ATLAS Offline Software
LArBarrelPresamplerGeometry.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2017 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 
27 #include "LArBarrelPresamplerGeometry.h"
28 #include "LArG4Code/LArVG4DetectorParameters.h"
29 #include "LArG4Code/LArG4Identifier.h"
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  declareProperty("DetectorName",m_detectorName);
53  declareProperty("TestBeam", m_testbeam);
54  }
55 
56  //=====================================================================================
57  Geometry::~Geometry() {;}
58 
59  // ==================================================================
62  StatusCode Geometry::initialize()
63  {
66  const double eps=0.007*Athena::Units::mm;
67  m_end_module[0]=(m_mod[0][0]*m_cmm+2*eps)+m_zminPS+eps;
68  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;
69 #ifdef DEBUGHITS
70  for (int i=0;i<8;i++) {
71  ATH_MSG_VERBOSE("Module length " << m_mod[0][0]*m_cmm+2*eps);
72  ATH_MSG_VERBOSE("End of Module " << m_end_module[i]);
73  }
74 #endif
75 
76  m_cat_th=m_cathode_th*m_cmm;
77  m_first_cathod[0]=m_zminPS+m_mod[0][5]*m_cmm+m_cat_th/2.+2*eps;
78  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;
79 
80 #ifdef DEBUGHITS
81  for (int i=0;i<8;i++) ATH_MSG_VERBOSE("position of first cathode " << m_first_cathod[i]);
82 #endif
83 
84  // number of cells in eta per module
85  for (int i=0;i<7;i++) m_ncell_module[i]=8;
86  m_ncell_module[7]=5;
87 
88  // electrode tild in rad
89  for (int i=0;i<8;i++) m_tilt[i]=m_mod[i][3]*CLHEP::deg;
90 
91  // number of gaps per cell module 7 is somewhat pathological (last cell is shorter)
92  for (int i=0;i<7;i++) m_ngap_cell[i]=(int)((m_mod[i][1]+0.1)/m_ncell_module[i]);
93  m_ngap_cell[7]=18;
94 #ifdef DEBUGHITS
95  for (int i=0;i<8;i++) ATH_MSG_VERBOSE("ngap per cell " << m_ngap_cell[i]);
96 #endif
97 
98  // pitch in z of gaps
99  for (int i=0;i<8;i++) m_pitch[i]=m_mod[i][4]*m_cmm;
100 
101  return StatusCode::SUCCESS;
102  }
103 
104  //======================================================================================
114  //======================================================================================
115 
116  LArG4Identifier Geometry::CalculateIdentifier(const G4Step* a_step) const
117  {
118  const static G4String fullPSName = (m_detectorName.empty()) ? "LAr::Barrel::Presampler" : G4String(m_detectorName+"::LAr::Barrel::Presampler");
119  const static G4String fullCryoName = (m_detectorName.empty()) ? "LAr::TBBarrel::Cryostat::LAr" : G4String(m_detectorName+"::LAr::TBBarrel::Cryostat::LAr");
120  const static G4String fullModuleName = (m_detectorName.empty()) ? "LAr::Barrel::Presampler::Module" : G4String(m_detectorName+"::LAr::Barrel::Presampler::Module");
121 
123  const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
124  const G4int ndep = g4navigation->GetDepth();
125  bool iactive(false);
126  G4int idep(-999);
127 
129  for (G4int ii=0;ii<=ndep;ii++) {
130  // FIXME Need to find a way to avoid these string-comparisons
131  const G4String& vname = g4navigation->GetVolume(ii)->GetName();
132  if (idep<0 && vname==fullPSName) idep=ii; // half barrel
133  else if (!iactive && vname==fullModuleName) iactive=true;
134  }
135 
136  if (idep < 0) std::abort();
137 
138  if ( iactive ) {
139  return this->CalculatePSActiveIdentifier( a_step , idep );
140  }
141  return this->CalculatePS_DMIdentifier( a_step , idep );
142  }
143 
144  // ==========================================================================================
150  LArG4Identifier Geometry::CalculatePSActiveIdentifier(const G4Step* a_step, const G4int ind) const
151  {
152  LArG4Identifier PSActiveID = LArG4Identifier();
153 
154  const G4ThreeVector p = (a_step->GetPostStepPoint()->GetPosition() + a_step->GetPreStepPoint()->GetPosition()) * 0.5;
155 
156 #ifdef DEBUGHITS
157  ATH_MSG_VERBOSE("Position of the step in the ATLAS frame (x,y,z) --> " << p.x() << " " << p.y() << " " << p.z());
158  ATH_MSG_VERBOSE("Eta and Phi in the atlas frame --> " << p.eta() << " " << p.phi());
159 #endif
160 
161  /* to get coordinates in the Presampler frame */
162 
163  const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
164  const G4ThreeVector ploc = g4navigation->GetTransform(ind).TransformPoint(p);
165 
166  const G4int zSide(this->determineZSide(p.z()));
167 
168  CalcData currentCellData;
169  (void)this->findCell(currentCellData,ploc.x(),ploc.y(),ploc.z());
170 
171  if( zSide == -1 )
173  {
174  currentCellData.phiBin = 31 - currentCellData.phiBin;
175  if(currentCellData.phiBin < 0 ) currentCellData.phiBin += 64;
176  }
177 
179  PSActiveID << 4 // LArCalorimeter
180  << 1 // LArEM
181  << zSide
182  << currentCellData.sampling
183  << currentCellData.region
184  << currentCellData.etaBin
185  << currentCellData.phiBin;
186 
187 #ifdef DEBUGHITS
188  ATH_MSG_VERBOSE("Here the identifier for the presampler ACTIVE ----> ");
189  ATH_MSG_VERBOSE("m_zSide ----> " << zSide);
190  ATH_MSG_VERBOSE("m_sampling ----> " << currentCellData.sampling);
191  ATH_MSG_VERBOSE("m_region ----> " << currentCellData.region);
192  ATH_MSG_VERBOSE("currentCellData.etaBin ----> " << currentCellData.etaBin);
193  ATH_MSG_VERBOSE("currentCellData.phiBin ----> " << currentCellData.phiBin);
194 #endif
195 
196  return PSActiveID ;
197  }
198 
199  //==========================================================================================
200 
201  LArG4Identifier Geometry::CalculatePS_DMIdentifier(const G4Step* a_step, const G4int ind) const
202  {
203 
204  /******************************************************************************
205  CaloDM_ID identifier:
206 
207  detector_system/subdet/type/sampling/region/eta/phi
208  detector system = 10 -> Calorimeters
209  subdet = +/-4 -> LAr dead materials
210  type = 1 -> dead materials outside accordion and active presampler layers
211  sampling = 1 -> dead materials in front and in active LAr calorimeters
212  (starting from front warm cryostat walls)
213  regions: = 0 barrel warm wall and solenoid in front of the barrel presampler, 0 < |eta| < 1.5
214  = 1 barrel cryostat cold wall in front of the barrel presampler, 0 < |eta| < 1.5
215  = 2 all materials in front of the barrel presampler at radius larger than cold wall
216  outer radius, 0 < |eta| < 1.5
217  = 3 all materials from the active layer of the barrel presampler to the active layer
218  of accordion, 0 < |eta| < 1.5
219 
220  ---> Granularity : deta 0.1 granularity within region
221  dphi pi/32 ~ 0.1 granularity within region
222 
223  ***********************************************************************************/
224 
225  const G4ThreeVector p = (a_step->GetPostStepPoint()->GetPosition() + a_step->GetPreStepPoint()->GetPosition()) * 0.5;
226 
227 #ifdef DEBUGHITS
228  ATH_MSG_VERBOSE("Position of the step in the ATLAS frame (x,y,z) --> " << p.x() << " " << p.y() << " " << p.z());
229  ATH_MSG_VERBOSE("Eta and Phi in the atlas frame --> " << p.eta() << " " << p.phi());
230 #endif
231 
235  const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
236  const G4ThreeVector ploc = g4navigation->GetTransform(ind).TransformPoint(p);
237  const G4double radius=sqrt(ploc.x()*ploc.x() + ploc.y()*ploc.y());
238 
240  const G4ThreeVector ploc2(ploc.x(),ploc.y(),ploc.z()+m_zpres+m_zminPS);
241 
242 #ifdef DEBUGHITS
243  ATH_MSG_VERBOSE("Position of the step after traslation (x,y,z) --> " << ploc2.x() << " " << ploc2.y() << " " << ploc2.z());
244  ATH_MSG_VERBOSE("Eta and Phi after translation --> " << ploc2.eta() << " " << ploc2.phi());
245 #endif
246 
247  // 01-Feb-2001 WGS: Add zSide calculation.
248  const G4int zSide(this->determineZSide(p.z()));
249 
251  const G4double phi = (ploc2.phi() < 0.) ? ploc2.phi()+2.*M_PI : ploc2.phi();
252  const G4double eta = ploc2.eta();
253  //G4double z2=fabs(ploc2.z());
254 
263  const G4int numberPhiMod = 32;
264  const G4double dphi = ( 2.*M_PI ) / numberPhiMod;
265  const G4double inv_dphi = 1. / dphi;
266  const G4double PSModuleRmean = 1420 ;
267  const G4double phicheck = phi - int(phi * inv_dphi) * dphi - (dphi /2.);
268  const G4double Rcheck = PSModuleRmean / cos(phicheck);
269  CalcData currentCellData;
270  if (radius > Rcheck) {
271  currentCellData.region = 3;
272  } else {
273  currentCellData.region = 2;
274  }
275 
276  const G4double detaDM = 0.1 ;
277  const G4double dphiDM = ( 2 * M_PI ) / 64. ;
278 
279  currentCellData.phiBin = G4int( phi * (1. / dphiDM) );
280  currentCellData.etaBin = G4int( eta * (1. / detaDM) );
281 
282  if( zSide == -1 )
283  {
284  currentCellData.phiBin = 31 - currentCellData.phiBin;
285  if(currentCellData.phiBin < 0 ) currentCellData.phiBin += 64;
286  }
287 
288  // 07-Jul-2005 WGS: Handle an extremely rare rounding problem.
289  if ( currentCellData.phiBin == 64 ) currentCellData.phiBin = 0;
290 
292  LArG4Identifier PS_DM_ID = LArG4Identifier();
293  PS_DM_ID << 10 // ATLAS
294  << zSide*4 // LArEM
295  << 1
296  << 1
297  << currentCellData.region
298  << currentCellData.etaBin
299  << currentCellData.phiBin;
300 
301 #ifdef DEBUGHITS
302  ATH_MSG_VERBOSE("Here the identifier for the presampler DEAD materials ---->");
303  ATH_MSG_VERBOSE("m_zSide ----> " << zSide*4);
304  ATH_MSG_VERBOSE("region ----> " << currentCellData.region);
305  ATH_MSG_VERBOSE("etaBin ----> " << currentCellData.etaBin);
306  ATH_MSG_VERBOSE("phiBin ----> " << currentCellData.phiBin);
307 #endif
308 
309  return PS_DM_ID ;
310  }
311 
338  bool Geometry::findCell(CalcData & currentCellData, G4double xloc,G4double yloc,G4double zloc) const
339  {
340 
341  currentCellData.sampling = 0;
342  currentCellData.region = 0;
343 
345  G4double phi = atan2(yloc,xloc);
346  if ( phi < 0. ) phi += 2.*M_PI;
347  const G4double z2=fabs(zloc+m_zpres+m_zminPS);
348 
350  const G4int numberPhiBins = 64;
351  const G4double dphi = ( 2.*M_PI ) / numberPhiBins;
352  const G4double inv_dphi = 1. / dphi;
354  currentCellData.phiBin = G4int( phi * inv_dphi );
355  if (currentCellData.phiBin >63) currentCellData.phiBin=63;
356  if (currentCellData.phiBin <0) currentCellData.phiBin=0;
357 
361  if (z2 < m_zminPS ) {
362  currentCellData.etaBin=0;
363  return false;
364  }
365  if (z2 > m_end_module[7]) {
366  currentCellData.etaBin=60;
367  return false;
368  }
369 
371  currentCellData.module=0;
372  for (int i=1;i<8;i++) {
373  if (m_first_cathod[i]>=z2) break;
374  currentCellData.module++;
375  }
376  if (currentCellData.module <0 || currentCellData.module > 7)
377  {
378  G4cerr << "LArBarrelPresampler: invalid module/hit " << currentCellData.module << " " << z2 << G4endl;
379  if (currentCellData.module<0) currentCellData.etaBin=0;
380  if (currentCellData.module>7) currentCellData.etaBin=60;
381  return false;
382  }
383 
385  const G4int isect=G4int(phi*m_nsectors/(2.*M_PI));
386  const G4double phi0= ((double)isect+0.5)*2.*M_PI/((double) m_nsectors);
387  static const G4double r0=1420.4*CLHEP::mm; // FIXME should be recomputed from database
388  currentCellData.dist=(xloc*cos(phi0)+yloc*sin(phi0)-r0);
389 
390 #ifdef DEBUGHITS
391  ATH_MSG_VERBOSE("sector number, dist along height " << isect << " " << currentCellData.dist);
392  ATH_MSG_VERBOSE("z2,module number,m_first_cathod " << z2 << " " << currentCellData.module << " "
393  << m_first_cathod[currentCellData.module]);
394 #endif
395 
398  G4double deltaz=z2-(m_first_cathod[currentCellData.module]+currentCellData.dist*tan(m_tilt[currentCellData.module]));
399  if (deltaz<0 ) {
400  if (currentCellData.module>0) {
401  currentCellData.module=currentCellData.module-1;
402  deltaz=z2-(m_first_cathod[currentCellData.module]+currentCellData.dist*tan(m_tilt[currentCellData.module]));
403  }
404  else deltaz=0;
405  }
406 
408  currentCellData.gap = ((int)(deltaz/m_pitch[currentCellData.module]));
409 
410 #ifdef DEBUGHITS
411  ATH_MSG_VERBOSE("deltaz from first cathode,gap number " << deltaz << " " << currentCellData.gap);
412 #endif
413 
415  currentCellData.etaBin= currentCellData.gap/m_ngap_cell[currentCellData.module];
416 #ifdef DEBUGHITS
417  ATH_MSG_VERBOSE("etaBin inside module " << currentCellData.etaBin);
418 #endif
419  if (currentCellData.etaBin >= m_ncell_module[currentCellData.module]) currentCellData.etaBin=m_ncell_module[currentCellData.module]-1;
420 
421  for (int i=0;i<currentCellData.module;i++) currentCellData.etaBin=currentCellData.etaBin+m_ncell_module[i];
422 #ifdef DEBUGHITS
423  ATH_MSG_VERBOSE(" final etaBin " << currentCellData.etaBin);
424 #endif
425 
426  if (currentCellData.etaBin < 0 || currentCellData.etaBin > 60) {
427  ATH_MSG_WARNING("LArBarrelPresamplerGeometry::findCell etaBin outside range " << currentCellData.etaBin);
428  }
429 
431  const G4double zmiddle=m_first_cathod[currentCellData.module]+((double)(currentCellData.gap+0.5))*m_pitch[currentCellData.module];
432 
436  currentCellData.xElec=currentCellData.dist*cos(m_tilt[currentCellData.module])+(z2-zmiddle)*sin(m_tilt[currentCellData.module]);
437  currentCellData.distElec=(z2-zmiddle)*cos(m_tilt[currentCellData.module]) - currentCellData.dist*sin(m_tilt[currentCellData.module]);
438 #ifdef DEBUGHITS
439  ATH_MSG_VERBOSE("zmiddle,xloc,yloc " << zmiddle << " " << currentCellData.distElec << " " << currentCellData.xElec);
440 #endif
441 
442  bool status=true;
443  if (fabs(currentCellData.dist)>m_halfThickLAr) {
444 #ifdef DEBUGHITS
445  ATH_MSG_VERBOSE("Outside normal LAr 13mm gap "),
446 #endif
447  status=false;
448  }
449 
450  return status;
451  }
452 
453  }
455 }
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.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
LArG4::BarrelPresampler::CalcData
Definition: ILArBarrelPresamplerGeometry.h:26
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
LArG4::BarrelPresampler::CalcData::phiBin
G4int phiBin
Definition: ILArBarrelPresamplerGeometry.h:30
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
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
lumiFormat.i
int i
Definition: lumiFormat.py:92
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:20
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:195
LArBarrelPresamplerCalculator::m_detectorName
std::string m_detectorName
Definition: LArBarrelPresamplerCalculator.h:61
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
python.SystemOfUnits.mm
int mm
Definition: SystemOfUnits.py:83
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
LArBarrelPresamplerCalculator::m_testbeam
bool m_testbeam
Definition: LArBarrelPresamplerCalculator.h:63
declareProperty
#define declareProperty(n, p, h)
Definition: BaseFakeBkgTool.cxx:15
merge.status
status
Definition: merge.py:17
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
m_cathode_th
G4double m_cathode_th
Electrodes.
Definition: PresParameterDef.h:13
checkFileSG.ind
list ind
Definition: checkFileSG.py:118
Generated on Fri Jul 5 2024 21:20:03 for ATLAS Offline Software by   doxygen 1.8.18