ATLAS Offline Software
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. More...
 

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. More...
 
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. More...
 
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. More...
 
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. More...
 

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 }

◆ 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 }

◆ 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 }

◆ 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

◆ 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 }

◆ 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 }

◆ 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; }

◆ 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:
Matrix
Definition: Trigger/TrigT1/TrigT1RPChardware/TrigT1RPChardware/Matrix.h:15
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
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:128
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
Monitored::Z
@ Z
Definition: HistogramFillerUtils.h:24
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
CaloCell_ID_FCS::TileExt0
@ TileExt0
Definition: FastCaloSim_CaloCell_ID.h:37
CaloCell_ID_FCS::TileBar1
@ TileBar1
Definition: FastCaloSim_CaloCell_ID.h:32
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
M_PI
#define M_PI
Definition: ActiveFraction.h:11
CaloExtensionHelpers::entryExitLayerMap
void entryExitLayerMap(const Trk::CaloExtension &extension, EntryExitLayerMap &result, const LayersToSelect *selection=nullptr)
Definition: CaloExtensionHelpers.h:192
xAOD::etaMax
etaMax
Definition: HIEventShape_v2.cxx:46
Phi
@ Phi
Definition: RPCdef.h:8
python.TurnDataReader.dr
dr
Definition: TurnDataReader.py:112
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
PathLengthUtils::phiMean
double phiMean(double a, double b) const
Definition: PathLengthUtils.h:52
x
#define x
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:49
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::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:113
CaloSampling::CaloSample
CaloSample
Definition: Calorimeter/CaloGeoHelpers/CaloGeoHelpers/CaloSampling.h:22
CaloCell_ID_FCS::TileBar0
@ TileBar0
Definition: FastCaloSim_CaloCell_ID.h:31
MuonR4::inverse
CalibratedSpacePoint::Covariance_t inverse(const CalibratedSpacePoint::Covariance_t &mat)
Inverts the parsed matrix.
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/UtilFunctions.cxx:65
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
AnalysisUtils::Delta::R
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
Definition: AnalysisMisc.h:49
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
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
query_example.col
col
Definition: query_example.py:7
a
TList * a
Definition: liststreamerinfos.cxx:10
LArCellBinning.etaMin
etaMin
Definition: LArCellBinning.py:84
y
#define y
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:567
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:26
PathLengthUtils::tileExit
static CaloSampling::CaloSample tileExit(CaloSampling::CaloSample sample)
Definition: PathLengthUtils.cxx:598
CaloPhiRange::diff
static double diff(double phi1, double phi2)
simple phi1 - phi2 calculation, but result is fixed to respect range.
Definition: CaloPhiRange.cxx:22
Generated on Mon Dec 23 2024 21:42:38 for ATLAS Offline Software by   doxygen 1.8.18