Load in the reco tracks and the hits.
78{
79
82 continue;
83 }
84
85 const double Eclus = shower.getEnergy();
87 continue;
88 }
89
90
91 const double thetaClus = shower.getTheta();
92 const double phiClus = shower.getPhi();
93 const double rClus = shower.getR();
94
95 const double sinTheta = sin(thetaClus);
96 const double cosTheta = cos(thetaClus);
97 const double sinPhi = sin(phiClus);
98 const double cosPhi = cos(phiClus);
99
100 ROOT::Math::XYZVector pos(rClus * sinTheta * cosPhi, rClus * sinTheta * sinPhi, rClus * cosTheta);
101 const double tanLambda = pos.Z() / pos.Rho();
102
103
106 continue;
107 }
108 }
109
110
112
113 ROOT::Math::XYZVector mom = Eclus / pos.R() * pos;
114
115
116
117
119
120
121 if (2. * rad < pos.Rho()) {
122 rad = pos.Rho() / 2.0 + 1.0;
123 }
124
125
126 double q = pos.Rho();
127 double y3 = pos.Y() / 2.0;
128 double x3 = pos.X() / 2.0;
129
130 double basex =
sqrt(rad * rad - q * q / 4.0) * (-pos.Y()) / q;
131 double basey =
sqrt(rad * rad - q * q / 4.0) * pos.X() / q;
132
133 double centerx1 = x3 + basex;
134 double centery1 = y3 + basey;
135 double centerx2 = x3 - basex;
136 double centery2 = y3 - basey;
137
138
139 double momx1 = pos.Y() - centery1;
140 double momy1 = - (pos.X() - centerx1);
141 double momx2 = pos.Y() - centery2;
142 double momy2 = - (pos.X() - centerx2);
143
144
145
146 if (momx1 * pos.X() + momy1 * pos.Y() < 0) {
147 momx1 = -momx1;
148 momy1 = -momy1;
149 }
150 if (momx2 * pos.X() + momy2 * pos.Y() < 0) {
151 momx2 = -momx2;
152 momy2 = -momy2;
153 }
154
155
156 double mom1abs =
sqrt(momx1 * momx1 + momy1 * momy1);
157 double mom2abs =
sqrt(momx2 * momx2 + momy2 * momy2);
158 momx1 = momx1 / mom1abs;
159 momy1 = momy1 / mom1abs;
160 momx2 = momx2 / mom2abs;
161 momy2 = momy2 / mom2abs;
162
163 ROOT::Math::XYZVector mom1(momx1 * mom.Rho(), momy1 * mom.Rho(), mom.Z());
164 ROOT::Math::XYZVector mom2(momx2 * mom.Rho(), momy2 * mom.Rho(), mom.Z());
165
166
167
168 bool clockwise1 = true;
169 bool clockwise2 = true;
170 if (pos.Y() * mom1.X() - pos.X() * mom1.Y() > 0) {
171 clockwise1 = false;
172 }
173 if (pos.Y() * mom2.X() - pos.X() * mom2.Y() > 0) {
174 clockwise2 = false;
175 }
176
177 if (clockwise1 == clockwise2) {
178 B2WARNING("Something went wrong during helix extrapolation. Skipping track.");
179 continue;
180 }
181
182 ROOT::Math::XYZVector mompos;
183 ROOT::Math::XYZVector momneg;
184 if (clockwise1) {
185 mompos = mom2;
186 momneg = mom1;
187 } else {
188 mompos = mom1;
189 momneg = mom2;
190 }
191
192
197
198
200 genfit::MeasuredStateOnPlane msopNeg(repNeg);
202 genfit::MeasuredStateOnPlane msopPos(repPos);
203
204
205 TMatrixDSym cov(6);
206 double covArray[6];
207 shower.getCovarianceMatrixAsArray(covArray);
208
209
210 double dx2 = rClus * cosTheta * cosPhi * rClus * cosTheta * cosPhi * covArray[5]
211 + rClus * sinTheta * sinPhi * rClus * sinTheta * sinPhi * covArray[2];
212 double dy2 = rClus * cosTheta * sinPhi * rClus * cosTheta * sinPhi * covArray[5]
213 + rClus * sinTheta * cosPhi * rClus * sinTheta * cosPhi * covArray[2];
214 double dz2 = rClus * sinTheta * rClus * sinTheta * covArray[5];
215
216 double dpx2 = std::abs(mom1.X() - mom2.X()) / 4.0 * std::abs(mom1.X() - mom2.X()) / 4.0;
217 double dpy2 = std::abs(mom1.Y() - mom2.Y()) / 4.0 * std::abs(mom1.Y() - mom2.Y()) / 4.0;
218 double dpz2 = 0.25 * 0.25 * mom.Z() * mom.Z();
219
220 cov(0, 0) = dx2;
221 cov(1, 1) = dy2;
222 cov(2, 2) = dz2;
223 cov(3, 3) = dpx2;
224 cov(4, 4) = dpy2;
225 cov(5, 5) = dpz2;
226
227
231 msopNeg.setPlane(planeNeg);
233 msopPos.setPlane(planePos);
234
235
236 CDCCKFState seedStateNeg(eclSeedNeg, msopNeg);
237 seeds.push_back({seedStateNeg});
238 CDCCKFState seedStatePos(eclSeedPos, msopPos);
239 seeds.push_back({seedStatePos});
240 }
241}
StoreArray< ECLShower > m_inputECLshowers
Input ECL Showers Store Array.
bool m_param_restrictToForwardSeeds
Don't do Ecl seeding in central region to save computing time.
double m_param_minimalEnRequirement
Minimal pt requirement.
StoreArray< RecoTrack > m_eclSeedRecoTracks
Output Reco Tracks Store Array.
double m_param_tanLambdaForwardNeg
Up to which (neg) tanLambda value should the seeding be performed.
double m_param_tanLambdaForwardPos
Up to which (pos) tanLambda value should the seeding be performed.
double m_param_showerDepth
Correction if the shower is assumed to start in a certain depth.
static const ChargedStable pion
charged pion particle
static const double speedOfLight
[cm/ns]
@ c_nPhotons
CR is split into n photons (N1)
static genfit::AbsTrackRep * createOrReturnRKTrackRep(RecoTrack &recoTrack, int PDGcode)
Checks if a TrackRap for the PDG id of the RecoTrack (and its charge conjugate) does already exit and...
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
static constexpr auto XYZToTVector
Helper function to convert XYZVector to TVector3.
double sqrt(double a)
sqrt for double
1.13.2