ATLAS Offline Software
Loading...
Searching...
No Matches
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//=================================
12PathLengthUtils::PathLengthUtils() = default;
13
14//==================================
15
16
17//=========================================================================================================================
18double 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}
337bool 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
355double 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
586CaloSampling::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
598CaloSampling::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
612double 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}
mag
Scalar mag() const
mag method
Definition AmgMatrixBasePlugin.h:26
AmgMatrix
#define AmgMatrix(rows, cols)
Definition EventPrimitives.h:49
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
Phi
@ Phi
Definition RPCdef.h:8
S1
struct TBPatternUnitContext S1
CaloCell_Base_ID::TILE
@ TILE
Definition CaloCell_Base_ID.h:45
CaloCell
Data object for each calorimeter readout cell.
Definition CaloCell.h:57
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::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
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
PathLengthUtils::get3DPathLength
static double get3DPathLength(const CaloCell &cell, const Amg::Vector3D &entry, const Amg::Vector3D &exit, double drFix, double dzFix)
Definition PathLengthUtils.cxx:18
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
CaloExtensionHelpers::EntryExitLayerMap
std::map< CaloSampling::CaloSample, std::pair< Amg::Vector3D, Amg::Vector3D > > EntryExitLayerMap
Definition CaloExtensionHelpers.h:21
type
Generated on for ATLAS Offline Software by doxygen 1.14.0