77 const unsigned int nSuper = 9;
78 double nLayer[nSuper] = { 8, 6, 6, 6, 6, 6, 6, 6, 6 };
79 double rInner[nSuper] = { 16.80, 25.70, 36.52, 47.69, 58.41, 69.53, 80.25,
81 double rOuter[nSuper] = { 23.80, 34.80, 45.57, 56.69, 67.41, 78.53, 89.25,
83 double stereo[nSuper] = { 0., 0.068, 0., -0.060, 0., 0.064, 0., -0.072, 0. };
84 double zStart[nSuper] = { -35.9, -51.4, -57.5, -59.6, -61.8, -63.9, -66.0,
86 double zEnd[nSuper] = { 67.9, 98.6, 132.9, 144.7, 146.8, 148.9, 151.0,
89 unsigned int nTry = 1000;
90 std::cout <<
" GblComposedKin $Id$ " << nTry <<
", " << nSuper << std::endl;
92 clock_t startTime = clock();
94 const double bfac = 0.0045;
96 double beamPos[] = { 0., 0., 0. };
97 double beamSize[] = { 0.005, 0.005, 0.3 };
98 const bool useBeamSpot =
true;
106 unsigned int iEvent = 0;
109 for (
unsigned int iTry = 0; iTry < nTry; ++iTry) {
111 std::vector<std::array<double, 6>> allGenPar;
117 const double qbyp = 0.2;
118 const double genDca = beamSize[0] *
unrm();
119 const double genZ0 = beamSize[2] *
unrm();
120 const double genPhi0 = M_PI * (2. *
unif() - 1.);
121 const double genDzds = 2. *
unif() - 1.;
122 const double genCurv = bfac * qbyp * sqrt(1. + genDzds * genDzds);
124 std::cout <<
" gen(inline) " << iEvent <<
": " << genCurv <<
" "
125 << genPhi0 <<
" " << genDca <<
" " << genDzds <<
" " << genZ0
127 std::array<double, 6> par = { genCurv, genPhi0, genDca, genDzds, genZ0,
129 allGenPar.push_back(par);
131 std::array<double, 6> par2 = { -genCurv, genPhi0 + M_PI, -genDca,
132 -genDzds, genZ0, -1. };
133 allGenPar.push_back(par2);
135 std::vector<std::vector<GblPoint>> listsOfPoints;
136 std::vector<Matrix<double, 5, 8>> listOfTrafos;
137 for (
unsigned int iTrack = 0; iTrack < nTrack; ++iTrack) {
141 std::array<double, 6> genPar = allGenPar[iTrack];
142 double curv(genPar[0]), phi0(genPar[1]), dca(genPar[2]), dzds(
143 genPar[3]), z0(genPar[4]);
147 std::vector<GblDetectorLayer> layers;
154 beamPos[2], phi0 + sArc * curv, dzds, beamSize[0],
155 beamSize[1], beamSize[2]));
157 unsigned int cLayer = 0;
158 for (
unsigned int iSuper = 0; iSuper < nSuper; ++iSuper) {
159 double radius = rInner[iSuper];
160 double step = (rOuter[iSuper] - radius) / (nLayer[iSuper] - 1);
162 for (
unsigned int iLayer = 0; iLayer < nLayer[iSuper];
165 if (fabs(dca) < radius) {
171 double phiPos = hlx.
getPhi(radius);
173 double xPos(cos(phiPos) * radius), yPos(
174 sin(phiPos) * radius), zPos(z0 + sArc * dzds);
176 if (zPos > zStart[iSuper] and zPos < zEnd[iSuper])
179 zPos, phi0 + sArc * curv, dzds,
180 stereo[iSuper], 0.015));
192 std::array<double, 6> seedPar = allGenPar[iTrack];
194 Matrix<double, 5, 8> innerTrafo = Matrix<double, 5, 8>::Zero();
195 if (allGenPar[iTrack][5] < 0.) {
197 innerTrafo(0, 5) = -1.;
198 innerTrafo(1, 6) = 1.;
199 innerTrafo(2, 7) = -1.;
200 innerTrafo(3, 3) = -1.;
201 innerTrafo(4, 4) = 1.;
204 innerTrafo(0, 0) = 1.;
205 innerTrafo(1, 1) = 1.;
206 innerTrafo(2, 2) = 1.;
207 innerTrafo(3, 3) = 1.;
208 innerTrafo(4, 4) = 1.;
210 listOfTrafos.push_back(innerTrafo);
214 double seedCurv(seedPar[0]), seedPhi0(seedPar[1]), seedDca(
215 seedPar[2]), seedDzds(seedPar[3]), seedZ0(seedPar[4]);
220 const double cosLambdaSeed = 1. / sqrt(1. + (seedDzds * seedDzds));
222 const double qbyp = seedCurv / bfac * cosLambdaSeed;
224 std::vector<GblPoint> listOfPoints;
225 for (
unsigned int iLayer = 0; iLayer < layers.size(); ++iLayer) {
235 Vector2d res = measPrediction;
239 res[0] += sigma[0] *
unrm();
241 std::cout <<
" impact par " << res[0] <<
" " << res[1]
242 <<
" " << seedCurv <<
" " << seedPhi0 <<
" "
243 << seedDca <<
" " << seedDzds <<
" " << seedZ0
244 <<
" " << sArc <<
" " << layers.size() << std::endl;
250 Matrix2d proM2l = transG2l * transM2g.block<3, 2>(0, 0);
252 Matrix2d proL2m = proM2l.inverse();
255 (sArc - sOld) / cosLambdaSeed, cosLambdaSeed, bfac);
259 if (lid > 0 or (iTrack == 0 and useBeamSpot))
262 std::vector<int> labGlobal(6);
273 for (
int p = 0; p < 6; p++)
274 labGlobal[p] = layerID * 1000 + p + 1;
276 pos, dir).block<2, 6>(0, 0);
283 Vector2d scat(0., 0.);
284 Vector2d scatPrec(1. / (errMs * errMs),
285 1. / (errMs * errMs));
289 listOfPoints.push_back(point);
291 listsOfPoints.push_back(listOfPoints);
295 unsigned int nAccepted = 0;
296 for (
unsigned int iTrack = 0; iTrack < nTrack; ++iTrack)
298 if (listsOfPoints[iTrack].size() >= 10)
304 std::vector<std::pair<std::vector<GblPoint>, MatrixXd> > aPointsAndTransList;
305 for (
unsigned int iTrack = 0; iTrack < nTrack; ++iTrack) {
306 aPointsAndTransList.push_back(
307 make_pair(listsOfPoints[iTrack], listOfTrafos[iTrack]));
319 unsigned int ierr = traj.
fit(Chi2, Ndf, lostWeight);
320 std::cout <<
" Fit " << iTry <<
": " << Chi2 <<
", " << Ndf <<
", "
321 << lostWeight << std::endl;
330 LostSum += lostWeight;
336 clock_t endTime = clock();
337 double diff = endTime - startTime;
338 double cps = CLOCKS_PER_SEC;
339 std::cout <<
" Time elapsed " << diff / cps <<
" s" << std::endl;
340 std::cout <<
" Chi2/Ndf = " << Chi2Sum / NdfSum << std::endl;
341 std::cout <<
" Tracks fitted " << numFit << std::endl;
343 std::cout <<
" Weight lost " << LostSum << std::endl;
GblTrajectory definition.
GblHelixPrediction intersectWithHelix(GblSimpleHelix hlx) const
Intersect with helix.
double getRadiationLength() const
Get radiation length.
unsigned int getLayerID() const
Get layer ID.
Eigen::Matrix3d getMeasSystemDirs() const
Get directions of measurement system.
Eigen::Matrix< double, 3, 6 > getRigidBodyDerGlobal(Eigen::Vector3d &position, Eigen::Vector3d &direction) const
Get rigid body derivatives in global frame.
Eigen::Vector2d getResolution() const
Get resolution.
Eigen::Vector2d getPrecision() const
Get precision.
Eigen::Matrix< double, 2, 3 > getCurvilinearDirs() const
Get curvilinear directions (U,V)
double getArcLength() const
Get arc-length.
const Eigen::Vector2d & getMeasPred() const
Get (measurement) prediction.
double getCosIncidence() const
Get cosine of incidence.
const Eigen::Vector3d & getPosition() const
Get position.
const Eigen::Vector3d & getDirection() const
Get position.
void addMeasurement(const Eigen::MatrixBase< Projection > &aProjection, const Eigen::MatrixBase< Residuals > &aResiduals, const Eigen::MatrixBase< Precision > &aPrecision, double minPrecision=0.)
Add a measurement to a point.
void addGlobals(const std::vector< int > &aLabels, const Eigen::MatrixBase< Derivative > &aDerivatives)
Add global derivatives to a point.
void addScatterer(const Eigen::Vector2d &aResiduals, const Eigen::Matrix2d &aPrecision)
Add a thin scatterer to a point.
double getArcLengthR(double aRadius) const
Get (2D) arc length for given radius (to ref. point)
double getPhi(double aRadius) const
Get phi (of point on circle) for given radius (to ref. point)
double getArcLengthXY(double xPos, double yPos) const
Get (2D) arc length for given point.
void printTrajectory(unsigned int level=0) const
Print GblTrajectory.
void milleOut(MilleBinary &aMille)
Write valid trajectory to Millepede-II binary file.
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, const std::string &optionList="", unsigned int aLabel=0)
Perform fit of (valid) trajectory.
Millepede-II (binary) record.
void exampleComposedKin()
Drift chamber example.
Definitions for exampleUtil(ities) for a Cylindrical Drift Chamber.
Namespace for the general broken lines package.
GblDetectorLayer CreateWireCdc(const std::string aName, unsigned int layer, double xPos, double yPos, double zPos, double phi, double tanLambda, double stereoAngle, double uRes)
Create a drift chamber wire with 1D measurement.
double unrm()
unit normal distribution, Box-Muller method, polar form
double gblMultipleScatteringError(double qbyp, double xbyx0)
Multiple scattering error.
GblDetectorLayer CreateImpactPar(const std::string aName, unsigned int layer, double xPos, double yPos, double zPos, double phi, double tanLambda, double xRes, double yRes, double zRes)
Create detector plane for impact parameters as 2D measurement.
double unif()
uniform distribution [0..1]
Matrix5d gblSimpleJacobian(double ds, double cosl, double bfac)
Simple jacobian.
Eigen::Matrix< double, 5, 5 > Matrix5d