ATLAS Offline Software
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Private Member Functions | Static Private Member Functions | List of all members
PathLengthUtils Class Reference

#include <PathLengthUtils.h>

Collaboration diagram for PathLengthUtils:

Public Member Functions

 PathLengthUtils ()
 ~PathLengthUtils ()=default
double pathInsideCell (const CaloCell &cell, const CaloExtensionHelpers::EntryExitLayerMap &entryExitLayerMap) const
 Return the length(mm) of the path crossed inside the cell, given the parameters for the extrapolation at entrance and exit of the layer.

Static Public Member Functions

static double get3DPathLength (const CaloCell &cell, const Amg::Vector3D &entry, const Amg::Vector3D &exit, double drFix, double dzFix)
static double getPathLengthInTile (const CaloCell &cell, const Amg::Vector3D &entry, const Amg::Vector3D &exit)

Private Member Functions

double phiMean (double a, double b) const
bool crossedPhi (const CaloCell &cell, double phi_entrance, double phi_exit) const
 Return true if the cell crossed was crossed by the track in phi.
double getPathLengthInEta (const CaloCell &cell, double eta_entrance, double eta_exit) const
 Return the % of path length crossed by the track inside a cell in eta.
double getPathLengthInZ (double zMin, double zMax, double z_entrance, double z_exit) const
 Return the % of path length crossed by the track inside a cell in Z.
double getPathLengthInZ (const CaloCell &cell, double z_entrance, double z_exit) const
 Return the % of path length crossed by the track inside a cell in Z.

Static Private Member Functions

static CaloSampling::CaloSample tileEntrance (CaloSampling::CaloSample sample)
static CaloSampling::CaloSample tileExit (CaloSampling::CaloSample sample)
static bool crossingMatrix (const AmgMatrix(3, 3)&Matrix, const Amg::Vector3D &entry, Amg::Vector3D &path)

Detailed Description

Definition at line 28 of file PathLengthUtils.h.

Constructor & Destructor Documentation

◆ PathLengthUtils()

PathLengthUtils::PathLengthUtils ( )
default

◆ ~PathLengthUtils()

PathLengthUtils::~PathLengthUtils ( )
default

Member Function Documentation

◆ crossedPhi()

bool PathLengthUtils::crossedPhi ( const CaloCell & cell,
double phi_entrance,
double phi_exit ) const
inlineprivate

Return true if the cell crossed was crossed by the track in phi.

Definition at line 55 of file PathLengthUtils.h.

55 {
56 double mean_phi = phiMean(phi_entrance, phi_exit);
57 double dphi = fabs(CaloPhiRange::diff(phi_entrance, phi_exit));
58 double phi_min = mean_phi - dphi, phi_max = mean_phi + dphi;
59
60 return (CaloPhiRange::diff(cell.phi() + cell.caloDDE()->dphi() / 2., phi_min) > 0 &&
61 CaloPhiRange::diff(phi_max, cell.phi() - cell.caloDDE()->dphi() / 2.) > 0);
62}
CaloPhiRange::diff
static double diff(double phi1, double phi2)
simple phi1 - phi2 calculation, but result is fixed to respect range.
Definition CaloPhiRange.cxx:22
PathLengthUtils::phiMean
double phiMean(double a, double b) const
Definition PathLengthUtils.h:52

◆ crossingMatrix()

bool PathLengthUtils::crossingMatrix ( const AmgMatrix(3, 3)& Matrix,
const Amg::Vector3D & entry,
Amg::Vector3D & path )
staticprivate

Definition at line 337 of file PathLengthUtils.cxx.

337 {
338
339 if (Matrix.determinant() == 0) {
340 return false;
341 }
342 AmgMatrix(3,3) Minv = Matrix.inverse();
343
344 path = Minv * entry;
345
346 bool crossing = false;
347 //
348 // Inside middle of the plane ? range -0.5 to 0.5
349 //
350 if (fabs(path.x()) < 0.5 && fabs(path.y()) < 0.5) crossing = true;
351
352 return crossing;
353}
AmgMatrix
#define AmgMatrix(rows, cols)
Definition EventPrimitives.h:49
if
if(febId1==febId2)
Definition LArRodBlockPhysicsV0.cxx:567
y
#define y
x
#define x
RTTAlgmain.Matrix
list Matrix
Definition RTTAlgmain.py:19

◆ get3DPathLength()

double PathLengthUtils::get3DPathLength ( const CaloCell & cell,
const Amg::Vector3D & entry,
const Amg::Vector3D & exit,
double drFix,
double dzFix )
static

Definition at line 18 of file PathLengthUtils.cxx.

19 {
20
21 double dmin = 10.;
22 double dphimin = dmin / cell.caloDDE()->r();
23
24 // Get cell parameters
25
26 double dphi = cell.caloDDE()->dphi();
27 if (dphi < dphimin) dphi = dphimin;
28 double CellPhimin = cell.caloDDE()->phi() - dphi * 0.5;
29 double CellPhimax = cell.caloDDE()->phi() + dphi * 0.5;
30
31 // Order according to IP to follow track convention with entry and exit planes
32
33 int isign = 1;
34 if (cell.caloDDE()->z() < 0) isign = -1;
35 double dr = cell.caloDDE()->dr();
36 //
37 // always use fixed values that correspond to the Calorimeter Tracking Geometry
38 // these are different from the CaloCell values
39 //
40 if (dr == 0) dr = drFix;
41 dr = drFix;
42 if (dr < dmin) dr = dmin;
43 double dz = cell.caloDDE()->dz();
44 if (dz == 0) dz = dzFix;
45 dz = dzFix;
46 if (dz < dmin) dz = dmin;
47
48 double CellZmin = cell.caloDDE()->z() - isign * dz * 0.5;
49 double CellZmax = cell.caloDDE()->z() + isign * dz * 0.5;
50 double CellRmin = cell.caloDDE()->r() - dr * 0.5;
51 double CellRmax = cell.caloDDE()->r() + dr * 0.5;
52
53 // For CPU reason do a fast check on crossing
54 //
55 // track direction
56 //
57 Amg::Vector3D dir = exit - entrance;
58 if (dir.mag() == 0) return 0.;
59 dir = dir / dir.mag();
60 //
61 // point of closest approach in x,y
62 //
63 double r0 = (entrance.x() - cell.caloDDE()->x()) * sin(dir.phi()) - (entrance.y() - cell.caloDDE()->y()) * cos(dir.phi());
64 bool crossing = true;
65 if (fabs(r0) > sqrt(cell.caloDDE()->r() * cell.caloDDE()->r() * dphi * dphi + dr * dr)) crossing = false;
66 //
67 // point of closest approach in r,z
68 //
69 double rz0 = (entrance.perp() - cell.caloDDE()->r()) * cos(dir.theta()) - (entrance.z() - cell.caloDDE()->z()) * sin(dir.theta());
70 if (fabs(rz0) > sqrt(dr * dr + dz * dz)) crossing = false;
71
72
73 if (!crossing) return 0.;
74
75 // construct 8 corners of the volume
76
77 Amg::Vector3D Corner0(CellRmin * cos(CellPhimin), CellRmin * sin(CellPhimin), CellZmin);
78 Amg::Vector3D Corner1(CellRmin * cos(CellPhimax), CellRmin * sin(CellPhimax), CellZmin);
79 Amg::Vector3D Corner2(CellRmin * cos(CellPhimin), CellRmin * sin(CellPhimin), CellZmax);
80 Amg::Vector3D Corner3(CellRmin * cos(CellPhimax), CellRmin * sin(CellPhimax), CellZmax);
81
82 Amg::Vector3D Corner4(CellRmax * cos(CellPhimin), CellRmax * sin(CellPhimin), CellZmin);
83 Amg::Vector3D Corner5(CellRmax * cos(CellPhimax), CellRmax * sin(CellPhimax), CellZmin);
84 Amg::Vector3D Corner6(CellRmax * cos(CellPhimin), CellRmax * sin(CellPhimin), CellZmax);
85 Amg::Vector3D Corner7(CellRmax * cos(CellPhimax), CellRmax * sin(CellPhimax), CellZmax);
86
87 // Entry plane vectors through Corner0
88
89 Amg::Vector3D dir01 = Corner1 - Corner0; // direction along Phi
90 Amg::Vector3D dir02 = Corner2 - Corner0; // direction along Z
91 // Second (side) plane through Corner0 has dir20 and dir40
92 Amg::Vector3D dir04 = Corner4 - Corner0; // direction along R
93 // Third plane through Corner0 has dir10 and dir40
94
95 // Exit plane vectors through Corner6
96
97 Amg::Vector3D dir67 = Corner7 - Corner6; // direction along Phi
98 Amg::Vector3D dir64 = Corner4 - Corner6; // direction along Z
99 // Second (side) plane through Corner dir62 and dir64
100 Amg::Vector3D dir62 = Corner2 - Corner6; // direction along R
101 // Third plane through Corner dir67 and dir62
102
103 // minus direction of track
104 Amg::Vector3D dirT = exit - entrance;
105 dirT = -dirT;
106
107
108
109 // dot products of track with plane vectors
110 double dotp1 = fabs(dirT.dot(dir01) / dir01.mag() / dirT.mag());
111 double dotp2 = fabs(dirT.dot(dir02) / dir02.mag() / dirT.mag());
112 double dotp4 = fabs(dirT.dot(dir04) / dir04.mag() / dirT.mag());
113 //
114 // order planes according to dotproduct and calculate Matrices
115 //mEntry matrix
116 // column(0) vector for first plane e.g. dir01
117 // column(1) second vector for first plane e.g. dir02
118 // column(2) = minus track direction
119 // type = 124 here plane made by 1 and 2 = dir01 and dir02
120 // later also plane made by 14 = dir01 and dir04
121 // can be tried (second try) or 24 = dir02 and
122 // dir04 (third try)
123 //
124 AmgMatrix(3,3) mEntry;
125 AmgMatrix(3,3) mExit;
126 mEntry.setZero();
127 mExit.setZero();
128 mEntry.col(2) = dirT;
129 mExit.col(2) = dirT;
130
131 int type = 0;
132 if (dotp1 <= dotp2 && dotp2 <= dotp4) {
133 type = 124;
134 mEntry.col(0) = dir01;
135 mEntry.col(1) = dir02;
136 mExit.col(0) = dir67;
137 mExit.col(1) = dir64;
138 } else if (dotp1 <= dotp4 && dotp4 <= dotp2) {
139 type = 142;
140 mEntry.col(0) = dir01;
141 mEntry.col(1) = dir04;
142 mExit.col(0) = dir67;
143 mExit.col(1) = dir62;
144 } else if (dotp2 <= dotp1 && dotp1 <= dotp4) {
145 type = 214;
146 mEntry.col(0) = dir02;
147 mEntry.col(1) = dir01;
148 mExit.col(0) = dir64;
149 mExit.col(1) = dir67;
150 } else if (dotp2 <= dotp4 && dotp4 <= dotp1) {
151 type = 241;
152 mEntry.col(0) = dir02;
153 mEntry.col(1) = dir04;
154 mExit.col(0) = dir64;
155 mExit.col(1) = dir62;
156 } else if (dotp4 <= dotp2 && dotp2 <= dotp1) {
157 type = 421;
158 mEntry.col(0) = dir04;
159 mEntry.col(1) = dir02;
160 mExit.col(0) = dir62;
161 mExit.col(1) = dir64;
162 } else if (dotp4 <= dotp1 && dotp1 <= dotp2) {
163 type = 412;
164 mEntry.col(0) = dir04;
165 mEntry.col(1) = dir01;
166 mExit.col(0) = dir62;
167 mExit.col(1) = dir67;
168 }
169
170 // position of track w.r.t. all planes
171
172 Amg::Vector3D posEn = entrance - (Corner1 + Corner2 + Corner4) / 2.;
173 Amg::Vector3D posEx = entrance - (Corner2 + Corner4 + Corner7) / 2.;
174 //
175 // put posEn en posEx in the middle of the plane
176 //
177 // this makes it possible to check that the Matrix solution lies within the bounds
178 // |x|<0.5 and |y|<0.5
179 // x corresponds to the first column(0) e.g. dir01
180 // y corresponds to the second column(1) e.g. dir02
181 //
182 int type0 = type - 10 * int(type / 10);
183 if (type0 == 1) {
184 // last plane 1 not used
185 posEn += Corner1 / 2.;
186 posEx += Corner7 / 2.;
187 }
188 if (type0 == 2) {
189 // last plane 2 not used
190 posEn += Corner2 / 2.;
191 posEx += Corner4 / 2.;
192 }
193 if (type0 == 4) {
194 // last plane 4 not used
195 posEn += Corner4 / 2.;
196 posEx += Corner2 / 2.;
197 }
198
199
200 // cross with entry plane
201 int typeCheck = type0;
202
203 Amg::Vector3D pathEn;
204 bool crossingEn = crossingMatrix(mEntry, posEn, pathEn);
205 // if dotmax near to 1 Matrix is singular (track is parallel to the plane)
206 if (!crossingEn) {
207 int type1 = int((type - 100 * int(type / 100)) / 10);
208 int type2 = int(type / 100);
209 if (type1 == 1) posEn += Corner1 / 2.;
210 if (type1 == 2) posEn += Corner2 / 2.;
211 if (type1 == 4) posEn += Corner4 / 2.;
212 if (type0 == 1) {
213 mEntry.col(1) = dir01;
214 posEn -= Corner1 / 2.;
215 }
216 if (type0 == 2) {
217 mEntry.col(1) = dir02;
218 posEn -= Corner2 / 2.;
219 }
220 if (type0 == 4) {
221 mEntry.col(1) = dir04;
222 posEn -= Corner4 / 2.;
223 }
224 crossingEn = crossingMatrix(mEntry, posEn, pathEn);
225 typeCheck = type1;
226 if (!crossingEn) {
227 if (type2 == 1) posEn += Corner1 / 2.;
228 if (type2 == 2) posEn += Corner2 / 2.;
229 if (type2 == 4) posEn += Corner4 / 2.;
230 if (type1 == 1) {
231 mEntry.col(0) = dir01;
232 posEn -= Corner1 / 2.;
233 }
234 if (type1 == 2) {
235 mEntry.col(0) = dir02;
236 posEn -= Corner2 / 2.;
237 }
238 if (type1 == 4) {
239 mEntry.col(0) = dir04;
240 posEn -= Corner4 / 2.;
241 }
242 crossingEn = crossingMatrix(mEntry, posEn, pathEn);
243 typeCheck = type2;
244 }
245 }
246
247 Amg::Vector3D posXEn = entrance;
248 if (crossingEn) {
249 posXEn = entrance - dirT * pathEn.z();
250 if (typeCheck == 1) {
251 double dphiCheck =
252 acos(cos(posXEn.phi()) * cos((CellPhimax + CellPhimin) / 2.) + sin(posXEn.phi()) * sin((CellPhimax + CellPhimin) / 2.));
253 if (dphiCheck > fabs(CellPhimax - CellPhimin)) { crossingEn = false; }
254 }
255 if (typeCheck == 2) {
256 if (posXEn.perp() < CellRmin) crossingEn = false;
257 if (posXEn.perp() > CellRmax) crossingEn = false;
258 }
259 if (typeCheck == 4) {
260 if (isign * posXEn.z() < isign * CellZmin) crossingEn = false;
261 if (isign * posXEn.z() > isign * CellZmax) crossingEn = false;
262 }
263 }
264
265 if (!crossingEn) return 0.;
266
267 // cross with exit plane
268
269 Amg::Vector3D pathEx;
270 bool crossingEx = crossingMatrix(mExit, posEx, pathEx);
271 typeCheck = type0;
272 if (!crossingEx) {
273 int type1 = int((type - 100 * int(type / 100)) / 10);
274 int type2 = int(type / 100);
275 if (type1 == 1) posEx += Corner7 / 2.;
276 if (type1 == 2) posEx += Corner4 / 2.;
277 if (type1 == 4) posEx += Corner2 / 2.;
278 if (type0 == 1) {
279 mExit.col(1) = dir67;
280 posEx -= Corner7 / 2.;
281 }
282 if (type0 == 2) {
283 mExit.col(1) = dir64;
284 posEx -= Corner4 / 2.;
285 }
286 if (type0 == 4) {
287 mExit.col(1) = dir62;
288 posEx -= Corner2 / 2.;
289 }
290 crossingEx = crossingMatrix(mExit, posEx, pathEx);
291 typeCheck = type1;
292 if (!crossingEx) {
293 if (type2 == 1) posEx += Corner7 / 2.;
294 if (type2 == 2) posEx += Corner4 / 2.;
295 if (type2 == 4) posEx += Corner2 / 2.;
296 if (type1 == 1) {
297 mExit.col(0) = dir67;
298 posEx -= Corner7 / 2.;
299 }
300 if (type1 == 2) {
301 mExit.col(0) = dir64;
302 posEx -= Corner4 / 2.;
303 }
304 if (type1 == 4) {
305 mExit.col(0) = dir62;
306 posEx -= Corner2 / 2.;
307 }
308 crossingEx = crossingMatrix(mExit, posEx, pathEx);
309 typeCheck = type2;
310 }
311 }
312
313 Amg::Vector3D posXEx = posXEn;
314 if (crossingEx) {
315 posXEx = entrance - dirT * pathEx.z();
316 if (typeCheck == 1) {
317 double dphiCheck =
318 acos(cos(posXEx.phi()) * cos((CellPhimax + CellPhimin) / 2.) + sin(posXEx.phi()) * sin((CellPhimax + CellPhimin) / 2.));
319 if (dphiCheck > fabs(CellPhimax - CellPhimin)) { crossingEx = false; }
320 }
321 if (typeCheck == 2) {
322 if (posXEx.perp() < CellRmin) crossingEx = false;
323 if (posXEx.perp() > CellRmax) crossingEx = false;
324 }
325 if (typeCheck == 4) {
326 if (isign * posXEx.z() < isign * CellZmin) crossingEx = false;
327 if (isign * posXEx.z() > isign * CellZmax) crossingEx = false;
328 }
329 }
330 double path = 0.;
331 if (crossingEn && crossingEx) {
332 path = (posXEx - posXEn).mag();
333 }
334
335 return path;
336}
mag
Scalar mag() const
mag method
Definition AmgMatrixBasePlugin.h:26
PathLengthUtils::crossingMatrix
static bool crossingMatrix(const AmgMatrix(3, 3)&Matrix, const Amg::Vector3D &entry, Amg::Vector3D &path)
Definition PathLengthUtils.cxx:337
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
TRT_PAI_physicsConstants::r0
const double r0
electron radius{cm}
Definition TRT_PAI_physicsConstants.h:22
athena.path
path
python interpreter configuration --------------------------------------—
Definition athena.py:128

◆ getPathLengthInEta()

double PathLengthUtils::getPathLengthInEta ( const CaloCell & cell,
double eta_entrance,
double eta_exit ) const
inlineprivate

Return the % of path length crossed by the track inside a cell in eta.

Definition at line 65 of file PathLengthUtils.h.

65 {
66 double etaMin = cell.eta() - 0.5 * cell.caloDDE()->deta();
67 double etaMax = cell.eta() + 0.5 * cell.caloDDE()->deta();
68 if (fabs(eta_entrance - eta_exit) < 1e-6) // to avoid FPE
69 return eta_entrance > etaMin && eta_entrance < etaMax;
70
71 double etaMinTrack = std::min(eta_entrance, eta_exit);
72 double etaMaxTrack = std::max(eta_entrance, eta_exit);
73 return (std::min(etaMax, etaMaxTrack) - std::max(etaMin, etaMinTrack)) / (etaMaxTrack - etaMinTrack);
74}

◆ getPathLengthInTile()

double PathLengthUtils::getPathLengthInTile ( const CaloCell & cell,
const Amg::Vector3D & entry,
const Amg::Vector3D & exit )
static

Definition at line 355 of file PathLengthUtils.cxx.

355 {
356 //=========================================================================================================================
357 // OBTAIN LAYER INDICES FOR LINEAR INTERPOLATION
358 unsigned int SampleID = cell.caloDDE()->getSampling();
359
360 constexpr double CellZB[9] = {141.495, 424.49, 707.485, 999.605, 1300.855, 1615.8, 1949., 2300.46, 2651.52};
361 constexpr double CellDZB[9] = {282.99, 283., 282.99, 301.25, 301.25, 328.64, 337.76, 365.16, 336.96};
362 constexpr double CellZC[9] = {159.755, 483.83, 812.465, 1150.23, 1497.125, 1857.71, 2241.12, 2628.695, 0};
363 constexpr double CellDZC[9] = {319.51, 328.64, 328.63, 346.9, 346.89, 374.28, 392.54, 382.61, 0};
364
365 // SPECIAL CASE: BC9
366 bool isBC9 = false;
367 if (SampleID == 13 && (fabs(cell.caloDDE()->eta() - 0.85) < 0.001 || fabs(cell.caloDDE()->eta() + 0.85) < 0.001)) isBC9 = true;
368
369 // OBTAIN TRACK AND CELL PARAMETERS
370 double pathl = 0.;
371
372 double Layer1X(exit.x()), Layer1Y(exit.y()), Layer1Z(exit.z());
373 double Layer2X(entrance.x()), Layer2Y(entrance.y()), Layer2Z(entrance.z());
374
375 double CellPhimin = cell.caloDDE()->phi() - cell.caloDDE()->dphi() * 0.5;
376 double CellPhimax = cell.caloDDE()->phi() + cell.caloDDE()->dphi() * 0.5;
377 double CellZmin = cell.caloDDE()->z() - cell.caloDDE()->dz() * 0.5;
378 double CellZmax = cell.caloDDE()->z() + cell.caloDDE()->dz() * 0.5;
379 double CellRmin = cell.caloDDE()->r() - cell.caloDDE()->dr() * 0.5;
380 double CellRmax = cell.caloDDE()->r() + cell.caloDDE()->dr() * 0.5;
381
382 // FIX FOR CELLS WITH ZERO WIDTH IN R
383 if (cell.caloDDE()->dr() == 0) {
384 // V = dr*dphi*r*dz
385 double CellVolume = cell.caloDDE()->volume();
386 double product = cell.caloDDE()->dphi() * cell.caloDDE()->r() * cell.caloDDE()->dz();
387 double drFormula = 0;
388 if (product != 0) { drFormula = CellVolume / product; } // IF
389 CellRmin = cell.caloDDE()->r() - drFormula * 0.5;
390 CellRmax = cell.caloDDE()->r() + drFormula * 0.5;
391 } // IF
392
393 // TileCellDim *cell_dim = m_tileDDM->get_cell_dim(cell->caloDDE()->identify());
394
395 double CellXimp[2], CellYimp[2], CellZimp[2];
396 double X(0), Y(0), Z(0), Phi(0);
397 double DeltaPhi;
398
399 // COMPUTE PATH
400 bool compute = true;
401 int lBC(0);
402 // LOOP IS USUALLY RUN ONCE, EXCEPT FOR LADDER SHAPED TILECAL CELLS
403 while (compute) {
404 if (lBC == 1 && isBC9) break;
405 int Np = 0;
406 if (sqrt((Layer1X - Layer2X) * (Layer1X - Layer2X) + (Layer1Y - Layer2Y) * (Layer1Y - Layer2Y)) < 3818.5) {
407 if (SampleID == 13 && lBC == 0) {
408 int cpos = fabs(cell.caloDDE()->eta()) / 0.1;
409
410 CellRmin = 2600.; // cell_dim->getRMin(0);
411 CellRmax = 2990.; // cell_dim->getRMax(2);
412 CellZmin = CellZB[cpos] - 0.5 * CellDZB[cpos]; // cell_dim->getZMin(0);
413 CellZmax = CellZB[cpos] + 0.5 * CellDZB[cpos]; // cell_dim->getZMax(0);
414 } // ELSE IF
415 else if (SampleID == 13 && lBC == 1) {
416 int cpos = fabs(cell.caloDDE()->eta()) / 0.1;
417
418 CellRmin = 2990.; // cell_dim->getRMin(3);
419 CellRmax = 3440.; // cell_dim->getRMax(5);
420 CellZmin = CellZC[cpos] - 0.5 * CellDZC[cpos]; // cell_dim->getZMin(3);
421 CellZmax = CellZC[cpos] + 0.5 * CellDZC[cpos]; // cell_dim->getZMax(3);
422 } // ELSE IF
423
424 // CALCULATE GRADIENTS
425 double Sxy = (Layer2X - Layer1X) / (Layer2Y - Layer1Y);
426 double Sxz = (Layer2X - Layer1X) / (Layer2Z - Layer1Z);
427 double Syz = (Layer2Y - Layer1Y) / (Layer2Z - Layer1Z);
428
429 // CALCULATE POINTS OF INTERSECTION
430 // INTERSECTIONS R PLANES
431 double Radius(CellRmin);
432
433 double x0int(exit.x()), x1int(entrance.x()), y0int(exit.y()), y1int(entrance.y()), z0int(exit.z()), z1int(entrance.z());
434 double S((y1int - y0int) / (x1int - x0int));
435 double a(1 + S * S), b(2 * S * y0int - 2 * S * S * x0int),
436 c(y0int * y0int - Radius * Radius + S * S * x0int * x0int - 2 * y0int * S * x0int);
437 double x1((-b + sqrt(b * b - 4 * a * c)) / (2 * a)), x2((-b - sqrt(b * b - 4 * a * c)) / (2 * a));
438 double y1(y0int + S * (x1 - x0int)), y2(y0int + S * (x2 - x0int));
439 double S1 = ((z1int - z0int) / (x1int - x0int));
440 double z1(z0int + S1 * (x1 - x0int)), z2(z0int + S1 * (x2 - x0int));
441
442 X = x1;
443 Y = y1;
444 Z = z1;
445
446 if (((x1 - x0int) * (x1 - x0int) + (y1 - y0int) * (y1 - y0int) + (z1 - z0int) * (z1 - z0int)) >
447 ((x2 - x0int) * (x2 - x0int) + (y2 - y0int) * (y2 - y0int) + (z2 - z0int) * (z2 - z0int))) {
448 X = x2;
449 Y = y2;
450 Z = z2;
451 } // IF
452
453 if (CellRmin != CellRmax) {
454 Phi = acos(X / sqrt(X * X + Y * Y));
455 if (Y <= 0) Phi = -Phi;
456 // R = CellRmin;
457
458 if (Z >= CellZmin && Z <= CellZmax && Phi >= CellPhimin && Phi <= CellPhimax) {
459 CellXimp[Np] = X;
460 CellYimp[Np] = Y;
461 CellZimp[Np] = Z;
462 Np = Np + 1;
463
464 } // IF
465
466 Radius = CellRmax;
467
468 c = y0int * y0int - Radius * Radius + S * S * x0int * x0int - 2 * y0int * S * x0int;
469 x1 = ((-b + sqrt(b * b - 4 * a * c)) / (2 * a));
470 x2 = ((-b - sqrt(b * b - 4 * a * c)) / (2 * a));
471 y1 = (y0int + S * (x1 - x0int));
472 y2 = (y0int + S * (x2 - x0int));
473 z1 = (z0int + S1 * (x1 - x0int));
474 z2 = (z0int + S1 * (x2 - x0int));
475 S1 = ((z1int - z0int) / (x1int - x0int));
476
477 X = x1;
478 Y = y1;
479 Z = z1;
480
481 if (((x1 - x0int) * (x1 - x0int) + (y1 - y0int) * (y1 - y0int) + (z1 - z0int) * (z1 - z0int)) >
482 ((x2 - x0int) * (x2 - x0int) + (y2 - y0int) * (y2 - y0int) + (z2 - z0int) * (z2 - z0int))) {
483 X = x2;
484 Y = y2;
485 Z = z2;
486 } // IF
487
488 Phi = std::acos(X / sqrt(X * X + Y * Y));
489 if (Y <= 0) Phi = -Phi;
490 // R=CellRmax;
491
492 if (Z >= CellZmin && Z <= CellZmax && Phi >= CellPhimin && Phi <= CellPhimax) {
493 CellXimp[Np] = X;
494 CellYimp[Np] = Y;
495 CellZimp[Np] = Z;
496 Np = Np + 1;
497 } // IF
498 } // IF
499
500 // INTERSECTIONS Z PLANES
501 if (CellZmin != CellZmax) {
502 if (Np < 2) {
503 Z = CellZmin;
504 X = Layer1X + Sxz * (Z - Layer1Z);
505 Y = Layer1Y + Syz * (Z - Layer1Z);
506 const double R = sqrt(X * X + Y * Y);
507 Phi = std::acos(X / R);
508 if (Y <= 0) Phi = -Phi;
509 if (R >= CellRmin && R <= CellRmax && Phi >= CellPhimin && Phi <= CellPhimax) {
510 CellXimp[Np] = X;
511 CellYimp[Np] = Y;
512 CellZimp[Np] = Z;
513 Np = Np + 1;
514 } // IF
515 } // IF
516
517 if (Np < 2) {
518 Z = CellZmax;
519 X = Layer1X + Sxz * (Z - Layer1Z);
520 Y = Layer1Y + Syz * (Z - Layer1Z);
521 const double R = sqrt(X * X + Y * Y);
522 Phi = std::acos(X / R);
523 if (Y <= 0) Phi = -Phi;
524 if (R >= CellRmin && R <= CellRmax && Phi >= CellPhimin && Phi <= CellPhimax) {
525 CellXimp[Np] = X;
526 CellYimp[Np] = Y;
527 CellZimp[Np] = Z;
528 Np = Np + 1;
529 } // IF
530 } // IF
531 } // IF
532
533 // INTERSECTIONS PHI PLANES
534 if (CellPhimin != CellPhimax) {
535 if (Np < 2) {
536 X = (Layer1X - Sxy * Layer1Y) / (1 - Sxy * tan(CellPhimin));
537 Y = X * std::tan(CellPhimin);
538 Z = Layer1Z + (1 / Sxz) * (X - Layer1X);
539 const double R = sqrt(X * X + Y * Y);
540 Phi = std::acos(X / R);
541 if (Y <= 0) Phi = -Phi;
542 DeltaPhi = fabs(Phi - CellPhimin);
543 if (DeltaPhi > 3.141593) DeltaPhi = fabs(Phi + CellPhimin);
544 if (R >= CellRmin && R <= CellRmax && Z >= CellZmin && Z <= CellZmax && DeltaPhi < 0.0001) {
545 CellXimp[Np] = X;
546 CellYimp[Np] = Y;
547 CellZimp[Np] = Z;
548 Np = Np + 1;
549 } // IF
550 } // IF
551
552 if (Np < 2) {
553 X = (Layer1X - Sxy * Layer1Y) / (1 - Sxy * tan(CellPhimax));
554 Y = X * std::tan(CellPhimax);
555 Z = Layer1Z + (1 / Sxz) * (X - Layer1X);
556 const double R = sqrt(X * X + Y * Y);
557 Phi = acos(X / R);
558 if (Y <= 0) Phi = -Phi;
559 DeltaPhi = fabs(Phi - CellPhimax);
560 if (DeltaPhi > 3.141593) DeltaPhi = fabs(Phi + CellPhimax);
561 if (R >= CellRmin && R <= CellRmax && Z >= CellZmin && Z <= CellZmax && DeltaPhi < 0.0001) {
562 CellXimp[Np] = X;
563 CellYimp[Np] = Y;
564 CellZimp[Np] = Z;
565 Np = Np + 1;
566 } // IF
567 } // IF
568 } // IF
569
570 // CALCULATE PATH IF TWO INTERSECTIONS WERE FOUND
571 if (Np == 2) {
572 pathl += sqrt((CellXimp[0] - CellXimp[1]) * (CellXimp[0] - CellXimp[1]) +
573 (CellYimp[0] - CellYimp[1]) * (CellYimp[0] - CellYimp[1]) +
574 (CellZimp[0] - CellZimp[1]) * (CellZimp[0] - CellZimp[1]));
575 } // IF
576 } // IF
577 if (SampleID == 13 && lBC == 0)
578 ++lBC;
579 else
580 compute = false;
581 } // WHILE (FOR LBBC LAYER)
582
583 return pathl;
584} // PathLengthUtils::getPathLengthInTile
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
Phi
@ Phi
Definition RPCdef.h:8
S1
struct TBPatternUnitContext S1
AnalysisUtils::Delta::R
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
Definition AnalysisMisc.h:49

◆ getPathLengthInZ() [1/2]

double PathLengthUtils::getPathLengthInZ ( const CaloCell & cell,
double z_entrance,
double z_exit ) const
inlineprivate

Return the % of path length crossed by the track inside a cell in Z.

Definition at line 87 of file PathLengthUtils.h.

87 {
88 return getPathLengthInZ(cell.z() - 0.5 * cell.caloDDE()->dz(), cell.z() + 0.5 * cell.caloDDE()->dz(), z_entrance, z_exit);
89}
PathLengthUtils::getPathLengthInZ
double getPathLengthInZ(double zMin, double zMax, double z_entrance, double z_exit) const
Return the % of path length crossed by the track inside a cell in Z.
Definition PathLengthUtils.h:77

◆ getPathLengthInZ() [2/2]

double PathLengthUtils::getPathLengthInZ ( double zMin,
double zMax,
double z_entrance,
double z_exit ) const
inlineprivate

Return the % of path length crossed by the track inside a cell in Z.

Definition at line 77 of file PathLengthUtils.h.

77 {
78 if (fabs(z_entrance - z_exit) < 1e-6) // to avoid FPE
79 return z_entrance > zMin && z_entrance < zMax;
80
81 double zMinTrack = std::min(z_entrance, z_exit);
82 double zMaxTrack = std::max(z_entrance, z_exit);
83 return (std::min(zMax, zMaxTrack) - std::max(zMin, zMinTrack)) / (zMaxTrack - zMinTrack);
84}

◆ pathInsideCell()

double PathLengthUtils::pathInsideCell ( const CaloCell & cell,
const CaloExtensionHelpers::EntryExitLayerMap & entryExitLayerMap ) const

Return the length(mm) of the path crossed inside the cell, given the parameters for the extrapolation at entrance and exit of the layer.

Definition at line 612 of file PathLengthUtils.cxx.

612 {
613 CaloSampling::CaloSample sample = cell.caloDDE()->getSampling();
614
615 //------------------------------------------------------------------------------------------
616 // special treatment for tile calo cells
617 CaloSampling::CaloSample tileEntranceID = tileEntrance(sample);
618 CaloSampling::CaloSample tileExitID = tileExit(sample);
619
620 auto tileEntrancePair = entryExitLayerMap.find(tileEntranceID);
621 auto tileExitPair = entryExitLayerMap.find(tileExitID);
622 if (tileEntrancePair != entryExitLayerMap.end() && tileExitPair != entryExitLayerMap.end()) {
623 return getPathLengthInTile(cell, tileEntrancePair->second.first, tileExitPair->second.second);
624 }
625 //------------------------------------------------------------------------------------------
626
627 auto entryExitPair = entryExitLayerMap.find(sample);
628 if (entryExitPair == entryExitLayerMap.end()) return -1.;
629
630 const Amg::Vector3D& entry = entryExitPair->second.first;
631 const Amg::Vector3D& exit = entryExitPair->second.second;
632
633 if (!crossedPhi(cell, entry.phi(), exit.phi())) return -1.;
634 double pathCrossed = -1;
635 if (cell.caloDDE()->getSubCalo() != CaloCell_ID::TILE) {
636 pathCrossed = getPathLengthInEta(cell, entry.eta(), exit.eta());
637 } else {
638 pathCrossed = getPathLengthInZ(cell, entry.z(), exit.z());
639 }
640 if (pathCrossed <= 0) return -1.;
641 return pathCrossed * (exit - entry).mag();
642}
CaloCell_Base_ID::TILE
@ TILE
Definition CaloCell_Base_ID.h:45
PathLengthUtils::getPathLengthInTile
static double getPathLengthInTile(const CaloCell &cell, const Amg::Vector3D &entry, const Amg::Vector3D &exit)
Definition PathLengthUtils.cxx:355
PathLengthUtils::crossedPhi
bool crossedPhi(const CaloCell &cell, double phi_entrance, double phi_exit) const
Return true if the cell crossed was crossed by the track in phi.
Definition PathLengthUtils.h:55
PathLengthUtils::tileEntrance
static CaloSampling::CaloSample tileEntrance(CaloSampling::CaloSample sample)
Definition PathLengthUtils.cxx:586
PathLengthUtils::getPathLengthInEta
double getPathLengthInEta(const CaloCell &cell, double eta_entrance, double eta_exit) const
Return the % of path length crossed by the track inside a cell in eta.
Definition PathLengthUtils.h:65
PathLengthUtils::tileExit
static CaloSampling::CaloSample tileExit(CaloSampling::CaloSample sample)
Definition PathLengthUtils.cxx:598
CaloExtensionHelpers::entryExitLayerMap
void entryExitLayerMap(const Trk::CaloExtension &extension, EntryExitLayerMap &result, const LayersToSelect *selection=nullptr)
Definition CaloExtensionHelpers.h:192
FullCPAlgorithmsTest_eljob.sample
sample
Definition FullCPAlgorithmsTest_eljob.py:116

◆ phiMean()

double PathLengthUtils::phiMean ( double a,
double b ) const
inlineprivate

Definition at line 52 of file PathLengthUtils.h.

52{ return 0.5 * (a + b) + (a * b < 0) * M_PI; }
M_PI
#define M_PI
Definition ActiveFraction.h:11

◆ tileEntrance()

CaloSampling::CaloSample PathLengthUtils::tileEntrance ( CaloSampling::CaloSample sample)
staticprivate

Definition at line 586 of file PathLengthUtils.cxx.

586 {
587 if (sample == CaloSampling::TileBar0 || sample == CaloSampling::TileBar1 || sample == CaloSampling::TileBar2 ||
588 sample == CaloSampling::TileGap2)
589 return CaloSampling::TileBar0;
590 if (sample == CaloSampling::TileGap1) return CaloSampling::TileBar1;
591 if (sample == CaloSampling::TileGap3 || sample == CaloSampling::TileExt0 || sample == CaloSampling::TileExt1 ||
592 sample == CaloSampling::TileExt2)
593 return CaloSampling::TileExt0;
594 return CaloSampling::Unknown;
595 // return sample;
596} // PathLengthUtils::entrance

◆ tileExit()

CaloSampling::CaloSample PathLengthUtils::tileExit ( CaloSampling::CaloSample sample)
staticprivate

Definition at line 598 of file PathLengthUtils.cxx.

598 {
599 if (sample == CaloSampling::TileBar0 || sample == CaloSampling::TileBar1 || sample == CaloSampling::TileBar2 ||
600 sample == CaloSampling::TileGap1)
601 return CaloSampling::TileBar2;
602 if (sample == CaloSampling::TileGap2) return CaloSampling::TileBar1;
603 if (sample == CaloSampling::TileGap3) return CaloSampling::TileExt1;
604 if (sample == CaloSampling::TileExt0 || sample == CaloSampling::TileExt1 || sample == CaloSampling::TileExt2)
605 return CaloSampling::TileExt2;
606 return CaloSampling::Unknown;
607 // return sample;
608} // PathLengthUtils::exit
The documentation for this class was generated from the following files:
Generated on for ATLAS Offline Software by doxygen 1.14.0