Load in the reco tracks and the hits.
76{
77
80 continue;
81 }
82
83 const double Eclus = shower.getEnergy();
85 continue;
86 }
87
88
89 const double thetaClus = shower.getTheta();
90 const double phiClus = shower.getPhi();
91 const double rClus = shower.getR();
92
93 const double sinTheta = sin(thetaClus);
94 const double cosTheta = cos(thetaClus);
95 const double sinPhi = sin(phiClus);
96 const double cosPhi = cos(phiClus);
97
98 ROOT::Math::XYZVector pos(rClus * sinTheta * cosPhi, rClus * sinTheta * sinPhi, rClus * cosTheta);
99 const double tanLambda = pos.Z() / pos.Rho();
100
101
104 continue;
105 }
106 }
107
108
110
111 ROOT::Math::XYZVector mom = Eclus / pos.R() * pos;
112
113
114
115
117
118
119 if (2. * rad < pos.Rho()) {
120 rad = pos.Rho() / 2.0 + 1.0;
121 }
122
123
124 double q = pos.Rho();
125 double y3 = pos.Y() / 2.0;
126 double x3 = pos.X() / 2.0;
127
128 double basex =
sqrt(rad * rad - q * q / 4.0) * (-pos.Y()) / q;
129 double basey =
sqrt(rad * rad - q * q / 4.0) * pos.X() / q;
130
131 double centerx1 = x3 + basex;
132 double centery1 = y3 + basey;
133 double centerx2 = x3 - basex;
134 double centery2 = y3 - basey;
135
136
137 double momx1 = pos.Y() - centery1;
138 double momy1 = - (pos.X() - centerx1);
139 double momx2 = pos.Y() - centery2;
140 double momy2 = - (pos.X() - centerx2);
141
142
143
144 if (momx1 * pos.X() + momy1 * pos.Y() < 0) {
145 momx1 = -momx1;
146 momy1 = -momy1;
147 }
148 if (momx2 * pos.X() + momy2 * pos.Y() < 0) {
149 momx2 = -momx2;
150 momy2 = -momy2;
151 }
152
153
154 double mom1abs =
sqrt(momx1 * momx1 + momy1 * momy1);
155 double mom2abs =
sqrt(momx2 * momx2 + momy2 * momy2);
156 momx1 = momx1 / mom1abs;
157 momy1 = momy1 / mom1abs;
158 momx2 = momx2 / mom2abs;
159 momy2 = momy2 / mom2abs;
160
161 ROOT::Math::XYZVector mom1(momx1 * mom.Rho(), momy1 * mom.Rho(), mom.Z());
162 ROOT::Math::XYZVector mom2(momx2 * mom.Rho(), momy2 * mom.Rho(), mom.Z());
163
164
165
166 bool clockwise1 = true;
167 bool clockwise2 = true;
168 if (pos.Y() * mom1.X() - pos.X() * mom1.Y() > 0) {
169 clockwise1 = false;
170 }
171 if (pos.Y() * mom2.X() - pos.X() * mom2.Y() > 0) {
172 clockwise2 = false;
173 }
174
175 if (clockwise1 == clockwise2) {
176 B2WARNING("Something went wrong during helix extrapolation. Skipping track.");
177 continue;
178 }
179
180 ROOT::Math::XYZVector mompos;
181 ROOT::Math::XYZVector momneg;
182 if (clockwise1) {
183 mompos = mom2;
184 momneg = mom1;
185 } else {
186 mompos = mom1;
187 momneg = mom2;
188 }
189
190
195
196
198 genfit::MeasuredStateOnPlane msopNeg(repNeg);
200 genfit::MeasuredStateOnPlane msopPos(repPos);
201
202
203 TMatrixDSym cov(6);
204 double covArray[6];
205 shower.getCovarianceMatrixAsArray(covArray);
206
207
208 double dx2 = rClus * cosTheta * cosPhi * rClus * cosTheta * cosPhi * covArray[5]
209 + rClus * sinTheta * sinPhi * rClus * sinTheta * sinPhi * covArray[2];
210 double dy2 = rClus * cosTheta * sinPhi * rClus * cosTheta * sinPhi * covArray[5]
211 + rClus * sinTheta * cosPhi * rClus * sinTheta * cosPhi * covArray[2];
212 double dz2 = rClus * sinTheta * rClus * sinTheta * covArray[5];
213
214 double dpx2 = std::abs(mom1.X() - mom2.X()) / 4.0 * std::abs(mom1.X() - mom2.X()) / 4.0;
215 double dpy2 = std::abs(mom1.Y() - mom2.Y()) / 4.0 * std::abs(mom1.Y() - mom2.Y()) / 4.0;
216 double dpz2 = 0.25 * 0.25 * mom.Z() * mom.Z();
217
218 cov(0, 0) = dx2;
219 cov(1, 1) = dy2;
220 cov(2, 2) = dz2;
221 cov(3, 3) = dpx2;
222 cov(4, 4) = dpy2;
223 cov(5, 5) = dpz2;
224
225
229 msopNeg.setPlane(planeNeg);
231 msopPos.setPlane(planePos);
232
233
234 CDCCKFState seedStateNeg(eclSeedNeg, msopNeg);
235 seeds.push_back({seedStateNeg});
236 CDCCKFState seedStatePos(eclSeedPos, msopPos);
237 seeds.push_back({seedStatePos});
238 }
239}
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