ATLAS Offline Software
LArBarrelPresamplerGeometry.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /********************************************************************
6 
7 NAME: LArBarrelPresamplerGeometry.cxx
8 PACKAGE: offline/LArCalorimeter/LArG4/LArG4Barrel
9 
10 AUTHORS: G. Unal, L. Carminati
11 CREATED: September, 2004
12 
13 PURPOSE: 'geometrical' methods used by the LArBarrelPresamplerCalculator.
14  These original methods (previously in LArBarrelPresampler Calculator) were
15  written by Bill Seligman
16 
17 UPDATES: - Calculate identifier method used by PresamplerCalibrationCalculator.
18  (sept 2004)
19  - Extended (added findCell method) to compute distance to electrode
20  (GU oct 2005). Also improved computation of etaBin taking into account
21  more correctly effect of tilted electrodes
22 
23 ********************************************************************/
24 
25 //#undef DEBUGHITS
26 
30 
31 #include "G4AffineTransform.hh"
32 #include "G4NavigationHistory.hh"
33 
34 #include "G4ThreeVector.hh"
35 #include "G4StepPoint.hh"
36 #include "G4Step.hh"
37 
38 #include "G4ios.hh"
39 #include "globals.hh"
40 
41 #include <cmath>
42 #include <limits.h>
43 #include <iostream>
44 
45 namespace LArG4 {
46 
47  namespace BarrelPresampler {
48 
49  Geometry::Geometry(const std::string& name, ISvcLocator *pSvcLocator)
50  : base_class(name, pSvcLocator)
51  {
52  }
53 
54  // ==================================================================
58  {
61  const double eps=0.007*Athena::Units::mm;
62  m_end_module[0]=(m_mod[0][0]*m_cmm+2*eps)+m_zminPS+eps;
63  for (int i=1;i<8;i++) m_end_module[i]=m_end_module[i-1]+(m_mod[i][0]*m_cmm+2*eps)+eps;
64 #ifdef DEBUGHITS
65  for (int i=0;i<8;i++) {
66  ATH_MSG_VERBOSE("Module length " << m_mod[0][0]*m_cmm+2*eps);
67  ATH_MSG_VERBOSE("End of Module " << m_end_module[i]);
68  }
69 #endif
70 
71  m_cat_th=m_cathode_th*m_cmm;
72  m_first_cathod[0]=m_zminPS+m_mod[0][5]*m_cmm+m_cat_th/2.+2*eps;
73  for (int i=1;i<8;i++) m_first_cathod[i]=m_end_module[i-1]+m_mod[i][5]*m_cmm+m_cat_th/2.+2*eps;
74 
75 #ifdef DEBUGHITS
76  for (int i=0;i<8;i++) ATH_MSG_VERBOSE("position of first cathode " << m_first_cathod[i]);
77 #endif
78 
79  // number of cells in eta per module
80  for (int i=0;i<7;i++) m_ncell_module[i]=8;
81  m_ncell_module[7]=5;
82 
83  // electrode tild in rad
84  for (int i=0;i<8;i++) m_tilt[i]=m_mod[i][3]*CLHEP::deg;
85 
86  // number of gaps per cell module 7 is somewhat pathological (last cell is shorter)
87  for (int i=0;i<7;i++) m_ngap_cell[i]=(int)((m_mod[i][1]+0.1)/m_ncell_module[i]);
88  m_ngap_cell[7]=18;
89 #ifdef DEBUGHITS
90  for (int i=0;i<8;i++) ATH_MSG_VERBOSE("ngap per cell " << m_ngap_cell[i]);
91 #endif
92 
93  // pitch in z of gaps
94  for (int i=0;i<8;i++) m_pitch[i]=m_mod[i][4]*m_cmm;
95 
96  return StatusCode::SUCCESS;
97  }
98 
99  //======================================================================================
109  //======================================================================================
110 
111  LArG4Identifier Geometry::CalculateIdentifier(const G4Step* a_step) const
112  {
113  const static G4String fullPSName = (m_detectorName.empty()) ? "LAr::Barrel::Presampler" : G4String(m_detectorName+"::LAr::Barrel::Presampler");
114  const static G4String fullCryoName = (m_detectorName.empty()) ? "LAr::TBBarrel::Cryostat::LAr" : G4String(m_detectorName+"::LAr::TBBarrel::Cryostat::LAr");
115  const static G4String fullModuleName = (m_detectorName.empty()) ? "LAr::Barrel::Presampler::Module" : G4String(m_detectorName+"::LAr::Barrel::Presampler::Module");
116 
118  const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
119  const G4int ndep = g4navigation->GetDepth();
120  bool iactive(false);
121  G4int idep(-999);
122 
124  for (G4int ii=0;ii<=ndep;ii++) {
125  // FIXME Need to find a way to avoid these string-comparisons
126  const G4String& vname = g4navigation->GetVolume(ii)->GetName();
127  if (idep<0 && vname==fullPSName) idep=ii; // half barrel
128  else if (!iactive && vname==fullModuleName) iactive=true;
129  }
130 
131  if (idep < 0) std::abort();
132 
133  if ( iactive ) {
134  return this->CalculatePSActiveIdentifier( a_step , idep );
135  }
136  return this->CalculatePS_DMIdentifier( a_step , idep );
137  }
138 
139  // ==========================================================================================
145  LArG4Identifier Geometry::CalculatePSActiveIdentifier(const G4Step* a_step, const G4int ind) const
146  {
148 
149  const G4ThreeVector p = (a_step->GetPostStepPoint()->GetPosition() + a_step->GetPreStepPoint()->GetPosition()) * 0.5;
150 
151 #ifdef DEBUGHITS
152  ATH_MSG_VERBOSE("Position of the step in the ATLAS frame (x,y,z) --> " << p.x() << " " << p.y() << " " << p.z());
153  ATH_MSG_VERBOSE("Eta and Phi in the atlas frame --> " << p.eta() << " " << p.phi());
154 #endif
155 
156  /* to get coordinates in the Presampler frame */
157 
158  const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
159  const G4ThreeVector ploc = g4navigation->GetTransform(ind).TransformPoint(p);
160 
161  const G4int zSide(this->determineZSide(p.z()));
162 
163  CalcData currentCellData;
164  (void)this->findCell(currentCellData,ploc.x(),ploc.y(),ploc.z());
165 
166  if( zSide == -1 )
168  {
169  currentCellData.phiBin = 31 - currentCellData.phiBin;
170  if(currentCellData.phiBin < 0 ) currentCellData.phiBin += 64;
171  }
172 
174  PSActiveID << 4 // LArCalorimeter
175  << 1 // LArEM
176  << zSide
177  << currentCellData.sampling
178  << currentCellData.region
179  << currentCellData.etaBin
180  << currentCellData.phiBin;
181 
182 #ifdef DEBUGHITS
183  ATH_MSG_VERBOSE("Here the identifier for the presampler ACTIVE ----> ");
184  ATH_MSG_VERBOSE("m_zSide ----> " << zSide);
185  ATH_MSG_VERBOSE("m_sampling ----> " << currentCellData.sampling);
186  ATH_MSG_VERBOSE("m_region ----> " << currentCellData.region);
187  ATH_MSG_VERBOSE("currentCellData.etaBin ----> " << currentCellData.etaBin);
188  ATH_MSG_VERBOSE("currentCellData.phiBin ----> " << currentCellData.phiBin);
189 #endif
190 
191  return PSActiveID ;
192  }
193 
194  //==========================================================================================
195 
196  LArG4Identifier Geometry::CalculatePS_DMIdentifier(const G4Step* a_step, const G4int ind) const
197  {
198 
199  /******************************************************************************
200  CaloDM_ID identifier:
201 
202  detector_system/subdet/type/sampling/region/eta/phi
203  detector system = 10 -> Calorimeters
204  subdet = +/-4 -> LAr dead materials
205  type = 1 -> dead materials outside accordion and active presampler layers
206  sampling = 1 -> dead materials in front and in active LAr calorimeters
207  (starting from front warm cryostat walls)
208  regions: = 0 barrel warm wall and solenoid in front of the barrel presampler, 0 < |eta| < 1.5
209  = 1 barrel cryostat cold wall in front of the barrel presampler, 0 < |eta| < 1.5
210  = 2 all materials in front of the barrel presampler at radius larger than cold wall
211  outer radius, 0 < |eta| < 1.5
212  = 3 all materials from the active layer of the barrel presampler to the active layer
213  of accordion, 0 < |eta| < 1.5
214 
215  ---> Granularity : deta 0.1 granularity within region
216  dphi pi/32 ~ 0.1 granularity within region
217 
218  ***********************************************************************************/
219 
220  const G4ThreeVector p = (a_step->GetPostStepPoint()->GetPosition() + a_step->GetPreStepPoint()->GetPosition()) * 0.5;
221 
222 #ifdef DEBUGHITS
223  ATH_MSG_VERBOSE("Position of the step in the ATLAS frame (x,y,z) --> " << p.x() << " " << p.y() << " " << p.z());
224  ATH_MSG_VERBOSE("Eta and Phi in the atlas frame --> " << p.eta() << " " << p.phi());
225 #endif
226 
230  const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
231  const G4ThreeVector ploc = g4navigation->GetTransform(ind).TransformPoint(p);
232  const G4double radius=sqrt(ploc.x()*ploc.x() + ploc.y()*ploc.y());
233 
235  const G4ThreeVector ploc2(ploc.x(),ploc.y(),ploc.z()+m_zpres+m_zminPS);
236 
237 #ifdef DEBUGHITS
238  ATH_MSG_VERBOSE("Position of the step after traslation (x,y,z) --> " << ploc2.x() << " " << ploc2.y() << " " << ploc2.z());
239  ATH_MSG_VERBOSE("Eta and Phi after translation --> " << ploc2.eta() << " " << ploc2.phi());
240 #endif
241 
242  // 01-Feb-2001 WGS: Add zSide calculation.
243  const G4int zSide(this->determineZSide(p.z()));
244 
246  const G4double phi = (ploc2.phi() < 0.) ? ploc2.phi()+2.*M_PI : ploc2.phi();
247  const G4double eta = ploc2.eta();
248  //G4double z2=fabs(ploc2.z());
249 
258  const G4int numberPhiMod = 32;
259  const G4double dphi = ( 2.*M_PI ) / numberPhiMod;
260  const G4double inv_dphi = 1. / dphi;
261  const G4double PSModuleRmean = 1420 ;
262  const G4double phicheck = phi - int(phi * inv_dphi) * dphi - (dphi /2.);
263  const G4double Rcheck = PSModuleRmean / cos(phicheck);
264  CalcData currentCellData;
265  if (radius > Rcheck) {
266  currentCellData.region = 3;
267  } else {
268  currentCellData.region = 2;
269  }
270 
271  const G4double detaDM = 0.1 ;
272  const G4double dphiDM = ( 2 * M_PI ) / 64. ;
273 
274  currentCellData.phiBin = G4int( phi * (1. / dphiDM) );
275  currentCellData.etaBin = G4int( eta * (1. / detaDM) );
276 
277  if( zSide == -1 )
278  {
279  currentCellData.phiBin = 31 - currentCellData.phiBin;
280  if(currentCellData.phiBin < 0 ) currentCellData.phiBin += 64;
281  }
282 
283  // 07-Jul-2005 WGS: Handle an extremely rare rounding problem.
284  if ( currentCellData.phiBin == 64 ) currentCellData.phiBin = 0;
285 
287  LArG4Identifier PS_DM_ID = LArG4Identifier();
288  PS_DM_ID << 10 // ATLAS
289  << zSide*4 // LArEM
290  << 1
291  << 1
292  << currentCellData.region
293  << currentCellData.etaBin
294  << currentCellData.phiBin;
295 
296 #ifdef DEBUGHITS
297  ATH_MSG_VERBOSE("Here the identifier for the presampler DEAD materials ---->");
298  ATH_MSG_VERBOSE("m_zSide ----> " << zSide*4);
299  ATH_MSG_VERBOSE("region ----> " << currentCellData.region);
300  ATH_MSG_VERBOSE("etaBin ----> " << currentCellData.etaBin);
301  ATH_MSG_VERBOSE("phiBin ----> " << currentCellData.phiBin);
302 #endif
303 
304  return PS_DM_ID ;
305  }
306 
333  bool Geometry::findCell(CalcData & currentCellData, G4double xloc,G4double yloc,G4double zloc) const
334  {
335 
336  currentCellData.sampling = 0;
337  currentCellData.region = 0;
338 
340  G4double phi = atan2(yloc,xloc);
341  if ( phi < 0. ) phi += 2.*M_PI;
342  const G4double z2=fabs(zloc+m_zpres+m_zminPS);
343 
345  const G4int numberPhiBins = 64;
346  const G4double dphi = ( 2.*M_PI ) / numberPhiBins;
347  const G4double inv_dphi = 1. / dphi;
349  currentCellData.phiBin = G4int( phi * inv_dphi );
350  if (currentCellData.phiBin >63) currentCellData.phiBin=63;
351  if (currentCellData.phiBin <0) currentCellData.phiBin=0;
352 
356  if (z2 < m_zminPS ) {
357  currentCellData.etaBin=0;
358  return false;
359  }
360  if (z2 > m_end_module[7]) {
361  currentCellData.etaBin=60;
362  return false;
363  }
364 
366  currentCellData.module=0;
367  for (int i=1;i<8;i++) {
368  if (m_first_cathod[i]>=z2) break;
369  currentCellData.module++;
370  }
371  if (currentCellData.module <0 || currentCellData.module > 7)
372  {
373  G4cerr << "LArBarrelPresampler: invalid module/hit " << currentCellData.module << " " << z2 << G4endl;
374  if (currentCellData.module<0) currentCellData.etaBin=0;
375  if (currentCellData.module>7) currentCellData.etaBin=60;
376  return false;
377  }
378 
380  const G4int isect=G4int(phi*m_nsectors/(2.*M_PI));
381  const G4double phi0= ((double)isect+0.5)*2.*M_PI/((double) m_nsectors);
382  static const G4double r0=1420.4*CLHEP::mm; // FIXME should be recomputed from database
383  currentCellData.dist=(xloc*cos(phi0)+yloc*sin(phi0)-r0);
384 
385 #ifdef DEBUGHITS
386  ATH_MSG_VERBOSE("sector number, dist along height " << isect << " " << currentCellData.dist);
387  ATH_MSG_VERBOSE("z2,module number,m_first_cathod " << z2 << " " << currentCellData.module << " "
388  << m_first_cathod[currentCellData.module]);
389 #endif
390 
393  G4double deltaz=z2-(m_first_cathod[currentCellData.module]+currentCellData.dist*tan(m_tilt[currentCellData.module]));
394  if (deltaz<0 ) {
395  if (currentCellData.module>0) {
396  currentCellData.module=currentCellData.module-1;
397  deltaz=z2-(m_first_cathod[currentCellData.module]+currentCellData.dist*tan(m_tilt[currentCellData.module]));
398  }
399  else deltaz=0;
400  }
401 
403  currentCellData.gap = ((int)(deltaz/m_pitch[currentCellData.module]));
404 
405 #ifdef DEBUGHITS
406  ATH_MSG_VERBOSE("deltaz from first cathode,gap number " << deltaz << " " << currentCellData.gap);
407 #endif
408 
410  currentCellData.etaBin= currentCellData.gap/m_ngap_cell[currentCellData.module];
411 #ifdef DEBUGHITS
412  ATH_MSG_VERBOSE("etaBin inside module " << currentCellData.etaBin);
413 #endif
414  if (currentCellData.etaBin >= m_ncell_module[currentCellData.module]) currentCellData.etaBin=m_ncell_module[currentCellData.module]-1;
415 
416  for (int i=0;i<currentCellData.module;i++) currentCellData.etaBin=currentCellData.etaBin+m_ncell_module[i];
417 #ifdef DEBUGHITS
418  ATH_MSG_VERBOSE(" final etaBin " << currentCellData.etaBin);
419 #endif
420 
421  if (currentCellData.etaBin < 0 || currentCellData.etaBin > 60) {
422  ATH_MSG_WARNING("LArBarrelPresamplerGeometry::findCell etaBin outside range " << currentCellData.etaBin);
423  }
424 
426  const G4double zmiddle=m_first_cathod[currentCellData.module]+((double)(currentCellData.gap+0.5))*m_pitch[currentCellData.module];
427 
431  currentCellData.xElec=currentCellData.dist*cos(m_tilt[currentCellData.module])+(z2-zmiddle)*sin(m_tilt[currentCellData.module]);
432  currentCellData.distElec=(z2-zmiddle)*cos(m_tilt[currentCellData.module]) - currentCellData.dist*sin(m_tilt[currentCellData.module]);
433 #ifdef DEBUGHITS
434  ATH_MSG_VERBOSE("zmiddle,xloc,yloc " << zmiddle << " " << currentCellData.distElec << " " << currentCellData.xElec);
435 #endif
436 
437  bool status=true;
438  if (fabs(currentCellData.dist)>m_halfThickLAr) {
439 #ifdef DEBUGHITS
440  ATH_MSG_VERBOSE("Outside normal LAr 13mm gap "),
441 #endif
442  status=false;
443  }
444 
445  return status;
446  }
447 
448  }
450 }
RoiUtil::phicheck
double phicheck(double phi)
basic range checkers
Definition: RoiUtil.cxx:123
LArG4Identifier
Definition: LArG4Identifier.h:121
LArG4::BarrelPresampler::CalcData::region
G4int region
Definition: ILArBarrelPresamplerGeometry.h:28
python.SystemOfUnits.mm
float mm
Definition: SystemOfUnits.py:98
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
LArG4::BarrelPresampler::CalcData
Definition: ILArBarrelPresamplerGeometry.h:26
LArG4::BarrelPresampler::CalcData::phiBin
G4int phiBin
Definition: ILArBarrelPresamplerGeometry.h:30
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
initialize
void initialize()
Definition: run_EoverP.cxx:894
InDetAccessor::phi0
@ phi0
Definition: InDetAccessor.h:33
LArBarrelPresamplerGeometry.h
M_PI
#define M_PI
Definition: ActiveFraction.h:11
LArVG4DetectorParameters.h
deg
#define deg
Definition: SbPolyhedron.cxx:17
LArG4::BarrelPresampler::Geometry::Geometry
Geometry(const std::string &name, ISvcLocator *pSvcLocator)
Definition: LArBarrelPresamplerGeometry.cxx:85
m_mod
G4double m_mod[8][6]
modules para [Length,NAnodes,NCathodes,elec.
Definition: PresParameterDef.h:19
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
m_nsectors
G4int m_nsectors
GU add phi min and phi span of overall volume.
Definition: PresParameterDef.h:60
LArG4
Definition: LArWheelCalculatorEnums.h:8
m_cmm
G4double m_cmm
default units : mm , deg.
Definition: PresParameterDef.h:7
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
lumiFormat.i
int i
Definition: lumiFormat.py:85
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
TRT_PAI_physicsConstants::r0
const double r0
electron radius{cm}
Definition: TRT_PAI_physicsConstants.h:22
LArG4Identifier.h
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
LArG4::BarrelPresampler::CalcData::sampling
G4int sampling
Definition: ILArBarrelPresamplerGeometry.h:27
LArG4::BarrelPresampler::CalcData::etaBin
G4int etaBin
Definition: ILArBarrelPresamplerGeometry.h:29
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
merge.status
status
Definition: merge.py:16
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
m_cathode_th
G4double m_cathode_th
Electrodes.
Definition: PresParameterDef.h:13