ATLAS Offline Software
PathLengthUtils.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "PathLengthUtils.h"
6 
7 #include <iostream>
8 
9 #include "CaloGeoHelpers/CaloSampling.h"
10 
11 //=================================
12 PathLengthUtils::PathLengthUtils() = default;
13 
14 //==================================
15 
16 
17 //=========================================================================================================================
18 double PathLengthUtils::get3DPathLength(const CaloCell& cell, const Amg::Vector3D& entrance, const Amg::Vector3D& exit, double drFix,
19  double dzFix) {
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 }
337 bool PathLengthUtils::crossingMatrix(const AmgMatrix(3,3)& Matrix, const Amg::Vector3D& entry, Amg::Vector3D& path) {
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 }
354 
355 double PathLengthUtils::getPathLengthInTile(const CaloCell& cell, const Amg::Vector3D& entrance, const Amg::Vector3D& exit) {
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
585 
586 CaloSampling::CaloSample PathLengthUtils::tileEntrance(CaloSampling::CaloSample sample) {
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
597 
598 CaloSampling::CaloSample PathLengthUtils::tileExit(CaloSampling::CaloSample sample) {
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
609 
612 double PathLengthUtils::pathInsideCell(const CaloCell& cell, const CaloExtensionHelpers::EntryExitLayerMap& entryExitLayerMap) const {
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 }
Matrix
Definition: Trigger/TrigT1/TrigT1RPChardware/TrigT1RPChardware/Matrix.h:15
GetLCDefs::Unknown
@ Unknown
Definition: GetLCDefs.h:21
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
S1
struct TBPatternUnitContext S1
CaloCell_ID_FCS::TileExt2
@ TileExt2
Definition: FastCaloSim_CaloCell_ID.h:39
ReadCellNoiseFromCool.cell
cell
Definition: ReadCellNoiseFromCool.py:53
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:126
IDTPM::R
float R(const U &p)
Definition: TrackParametersHelper.h:101
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
Monitored::Z
@ Z
Definition: HistogramFillerUtils.h:24
CaloCell_ID_FCS::TileExt0
@ TileExt0
Definition: FastCaloSim_CaloCell_ID.h:37
CaloCell_ID_FCS::TileBar1
@ TileBar1
Definition: FastCaloSim_CaloCell_ID.h:32
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
CaloExtensionHelpers::entryExitLayerMap
void entryExitLayerMap(const Trk::CaloExtension &extension, EntryExitLayerMap &result, const LayersToSelect *selection=nullptr)
Definition: CaloExtensionHelpers.h:192
Phi
@ Phi
Definition: RPCdef.h:8
python.TurnDataReader.dr
dr
Definition: TurnDataReader.py:112
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
JetTiledMap::S
@ S
Definition: TiledEtaPhiMap.h:44
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
CaloCell_ID_FCS::TileGap3
@ TileGap3
Definition: FastCaloSim_CaloCell_ID.h:36
AmgMatrix
#define AmgMatrix(rows, cols)
Definition: EventPrimitives.h:51
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
PathLengthUtils::crossingMatrix
static bool crossingMatrix(const AmgMatrix(3, 3)&Matrix, const Amg::Vector3D &entry, Amg::Vector3D &path)
Definition: PathLengthUtils.cxx:337
Monitored::X
@ X
Definition: HistogramFillerUtils.h:24
PathLengthUtils.h
PathLengthUtils::getPathLengthInTile
static double getPathLengthInTile(const CaloCell &cell, const Amg::Vector3D &entry, const Amg::Vector3D &exit)
Definition: PathLengthUtils.cxx:355
FullCPAlgorithmsTest_eljob.sample
sample
Definition: FullCPAlgorithmsTest_eljob.py:100
CaloSampling::CaloSample
CaloSample
Definition: Calorimeter/CaloGeoHelpers/CaloGeoHelpers/CaloSampling.h:22
CaloCell_ID_FCS::TileBar0
@ TileBar0
Definition: FastCaloSim_CaloCell_ID.h:31
CaloCell_ID_FCS::TileGap2
@ TileGap2
Definition: FastCaloSim_CaloCell_ID.h:35
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
TRT_PAI_physicsConstants::r0
const double r0
electron radius{cm}
Definition: TRT_PAI_physicsConstants.h:20
CaloCell_Base_ID::TILE
@ TILE
Definition: CaloCell_Base_ID.h:46
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
CaloExtensionHelpers::EntryExitLayerMap
std::map< CaloSampling::CaloSample, std::pair< Amg::Vector3D, Amg::Vector3D > > EntryExitLayerMap
Definition: CaloExtensionHelpers.h:21
CaloCell_ID_FCS::TileGap1
@ TileGap1
Definition: FastCaloSim_CaloCell_ID.h:34
calibdata.exit
exit
Definition: calibdata.py:236
beamspotman.dir
string dir
Definition: beamspotman.py:623
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
Monitored::Y
@ Y
Definition: HistogramFillerUtils.h:24
GetAllXsec.entry
list entry
Definition: GetAllXsec.py:132
PathLengthUtils::get3DPathLength
static double get3DPathLength(const CaloCell &cell, const Amg::Vector3D &entry, const Amg::Vector3D &exit, double drFix, double dzFix)
Definition: PathLengthUtils.cxx:18
CaloCell_ID_FCS::TileExt1
@ TileExt1
Definition: FastCaloSim_CaloCell_ID.h:38
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
PathLengthUtils::PathLengthUtils
PathLengthUtils()
a
TList * a
Definition: liststreamerinfos.cxx:10
CaloCell
Data object for each calorimeter readout cell.
Definition: CaloCell.h:57
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
PathLengthUtils::tileEntrance
static CaloSampling::CaloSample tileEntrance(CaloSampling::CaloSample sample)
Definition: PathLengthUtils.cxx:586
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
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
python.compressB64.c
def c
Definition: compressB64.py:93
CaloCell_ID_FCS::TileBar2
@ TileBar2
Definition: FastCaloSim_CaloCell_ID.h:33
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:25
PathLengthUtils::pathInsideCell
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...
Definition: PathLengthUtils.cxx:612
PathLengthUtils::tileExit
static CaloSampling::CaloSample tileExit(CaloSampling::CaloSample sample)
Definition: PathLengthUtils.cxx:598
Generated on Sun Jun 30 2024 21:24:00 for ATLAS Offline Software by   doxygen 1.8.18