ATLAS Offline Software
PixelSensorSD.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 //
6 // Pixel Sensitive Detector.
7 // The Hits are processed here. For every hit I get the position and
8 // an information on the sensor in which the interaction happened
9 //
10 
11 // Class header
12 #include "PixelSensorSD.h"
13 
14 // Athena headers
16 #include "MCTruth/TrackHelper.h"
17 
18 // Geant4 headers
19 #include "G4ChargedGeantino.hh"
20 #include <G4EventManager.hh>
21 #include "G4Geantino.hh"
22 #include "G4SDManager.hh"
23 #include "G4Step.hh"
24 #include "G4ThreeVector.hh"
25 
26 // CLHEP headers
27 #include "CLHEP/Geometry/Transform3D.h"
28 #include "CLHEP/Units/SystemOfUnits.h"
29 
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 
32 PixelSensorSD::PixelSensorSD(const std::string& name, const std::string& hitCollectionName)
33  : G4VSensitiveDetector( name )
34  , m_HitCollName( hitCollectionName )
35 {
36 }
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
39 // Initialize from G4 - cache the hit collection for the current event
40 void PixelSensorSD::Initialize(G4HCofThisEvent *)
41 {
42  // ISF calls G4SDManager::PrepareNewEvent() before the Geant4 event loop starts...
43  if(auto* eventManger = G4EventManager::GetEventManager()){
44  if(auto* eventInfo = static_cast<AtlasG4EventUserInfo*>(eventManger->GetUserInformation())){
45  m_HitColl = eventInfo->GetHitCollectionMap()->Find<SiHitCollection>(m_HitCollName);
46  }
47  }
48 }
49 
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 
52 G4bool PixelSensorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* /*ROhist*/)
53 {
54  if (verboseLevel>5) G4cout << "Process Hit" << G4endl;
55 
56  G4double edep = aStep->GetTotalEnergyDeposit();
57  edep *= CLHEP::MeV;
58  if(edep==0.) {
59  if(aStep->GetTrack()->GetDefinition() != G4Geantino::GeantinoDefinition() &&
60  aStep->GetTrack()->GetDefinition() != G4ChargedGeantino::ChargedGeantinoDefinition())
61  return false;
62  }
63 
64  //use the global time. i.e. the time from the beginning of the event
65  //
66  // Get the Touchable History:
67  //
68  const G4TouchableHistory *myTouch = dynamic_cast<const G4TouchableHistory*>(aStep->GetPreStepPoint()->GetTouchable());
69  if (not myTouch) {
70  G4cout << "PixelSensorSD::ProcessHits bad dynamic_cast" << G4endl;
71  return false;
72  }
73  if(verboseLevel>5){
74  for (int i=0;i<myTouch->GetHistoryDepth();i++){
75  std::string detname=myTouch->GetVolume(i)->GetLogicalVolume()->GetName();
76  int copyno=myTouch->GetVolume(i)->GetCopyNo();
77  G4cout << "Volume " <<detname <<" Copy Nr. " << copyno << G4endl;
78  }
79  }
80  //
81  // Get the hit coordinates. Start and End Point
82  //
83  G4ThreeVector coord1 = aStep->GetPreStepPoint()->GetPosition();
84  G4ThreeVector coord2 = aStep->GetPostStepPoint()->GetPosition();
85 
86  // Calculate the local step begin and end position.
87  // From a G4 FAQ:
88  // http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/geometry/17/1.html
89  //
90  const G4AffineTransform transformation = myTouch->GetHistory()->GetTopTransform();
91  G4ThreeVector localPosition1 = transformation.TransformPoint(coord1);
92  G4ThreeVector localPosition2 = transformation.TransformPoint(coord2);
93 
94  HepGeom::Point3D<double> lP1,lP2;
95  lP1[SiHit::xEta] = localPosition1[2]*CLHEP::mm;
96  lP1[SiHit::xPhi] = localPosition1[1]*CLHEP::mm;
97  lP1[SiHit::xDep] = localPosition1[0]*CLHEP::mm;
98 
99  lP2[SiHit::xEta] = localPosition2[2]*CLHEP::mm;
100  lP2[SiHit::xPhi] = localPosition2[1]*CLHEP::mm;
101  lP2[SiHit::xDep] = localPosition2[0]*CLHEP::mm;
102 
103  //
104  // Get it into a vector in local coords and with the right units:
105 
106 
107  // Now Navigate the history to know in what detector the step is:
108  // and finally set the ID of det element in which the hit is.
109  //
110  //G4int History;
111  //
112  // Is it TB, barrel or endcap?
113  //
114  int BrlEcap = 0;
115  int LayerDisk = 0;
116  int etaMod = 0;
117  int phiMod = 0;
118  int side = 0;
119 
120  // G4VPhysicalVolume* TheVolume = myTouch->GetVolume();
121  //
122  // In the case of the TB the positioning integers won't be initialized
123  // and the identifying integer will be zero. There is no need to do
124  // anything else
125 
126  int BEcopyNo = myTouch->GetVolume()->GetCopyNo();
127 
128  // Barrel
129  if(BEcopyNo == 100) {
130  // Move up of one level in the tree to get the module number:
131  //
132  etaMod = myTouch->GetVolume(1)->GetCopyNo();
133  //
134  // Another level for the ladder
135  //
136  phiMod = myTouch->GetVolume(2)->GetCopyNo();
137  //
138  // Another one for the layer
139  //
140  LayerDisk = myTouch->GetVolume(3)->GetCopyNo();
141  //
142  // Some printout
143  //
144  if (verboseLevel>5){
145  G4cout << "In the Pixel Barrel" << G4endl;
146  G4cout << "----- Phi Module # " << phiMod << G4endl;
147  G4cout << "----- Eta Ladder # " << etaMod << G4endl;
148  G4cout << "----- Layer # " << LayerDisk << G4endl;
149  }
150 
151  // Standard ATLAS EndCap
152  } else if (BEcopyNo == 200) {
153 
154  // Use the copy no. to get the id:
155  phiMod = myTouch->GetVolume(1)->GetCopyNo();
156  if (phiMod == 48) phiMod = 0;
157 
158  // etaMod always 0.
159 
160  // Move up of one level in the tree to get the disk number:
161  LayerDisk = myTouch->GetVolume(2)->GetCopyNo();
162 
163  // Move up of one level in the tree to get the EC number:
164  BrlEcap = myTouch->GetVolume(3)->GetCopyNo();
165  //The following is no longer necessary
166  //if(BrlEcap == -2) BrlEcap = 1;
167 
168  // Workaround for bug in cosmic setup
169  if (BrlEcap == 0) {
170  BrlEcap = (coord1.z() > 0) ? 2 : -2;
171  }
172 
173  // Some printout
174  if (verboseLevel>5){
175  G4cout << "In the Pixel EndCap" << G4endl;
176  G4cout << "----- PhiModule # " << phiMod << G4endl;
177  G4cout << "----- Disk # " << LayerDisk << G4endl;
178  G4cout << "----- Endcap # " << BrlEcap << G4endl;
179  }
180 
181  // SLHC EndCap
182  } else if (BEcopyNo == 300) {
183 
184  // Use the copy no. to get the id:
185  phiMod = myTouch->GetVolume(1)->GetCopyNo();
186 
187  // Note: copyNo is 2*ieta for front face, 2*ieta+1 for back face
188  etaMod = myTouch->GetVolume(2)->GetCopyNo()/2;
189 
190  // Move up of one level in the tree to get the disk number:
191  LayerDisk = myTouch->GetVolume(3)->GetCopyNo();
192 
193  // Move up of one level in the tree to get the EC number:
194  BrlEcap = myTouch->GetVolume(4)->GetCopyNo();
195  //if(BrlEcap == -2) BrlEcap = 1;
196 
197  // Some printout
198  if (verboseLevel>5){
199  G4cout << "In the SLHC Pixel EndCap" << G4endl;
200  G4cout << "----- PhiModule # " << phiMod << G4endl;
201  G4cout << "----- Ring/Eta # " << etaMod << G4endl;
202  G4cout << "----- Disk # " << LayerDisk << G4endl;
203  G4cout << "----- Endcap # " << BrlEcap << G4endl;
204  }
205 
206  } else if(BEcopyNo == 400) { //DBM
207 
208  // Note: there is one volume (3-layers unit) contain the diamond+chip volume
209  // so one increment is needed after 'layerDisk'
210 
211  // only one eta module for the DBM
212  etaMod = 0;
213 
214  // Move up one level to get the layer number, 3 layers per telescope
215  LayerDisk = myTouch->GetVolume(1)->GetCopyNo();
216 
217  // Move up two level to get the phi module, azimuthal position of telescope
218  phiMod = myTouch->GetVolume(3)->GetCopyNo();
219 
220  // Move up one level to get the DBM side, A side or C side
221  BrlEcap = myTouch->GetVolume(4)->GetCopyNo();
222 
223  // Some printout
224  if (verboseLevel>5){
225  G4cout << "In the DBM" << G4endl;
226  G4cout << "----- PhiModule # " << phiMod << G4endl;
227  G4cout << "----- Ring/Eta # " << etaMod << G4endl;
228  G4cout << "----- Disk # " << LayerDisk << G4endl;
229  G4cout << "----- Endcap # " << BrlEcap << G4endl;
230  }
231 
232  } else if(BEcopyNo == 500) { //Inclined
233 
234  std::string volName = myTouch->GetVolume()->GetName();
235 
236  // split the volname vs '_'
237  std::replace(volName.begin(),volName.end(),'_',' ');
238 
239  std::vector<std::string> v;
240  std::istringstream s(volName);
241  std::string tmp;
242  while (s >> tmp) {
243  v.push_back(tmp);
244  }
245 
246  BrlEcap = 0; // no endcap defined in alpine
247  LayerDisk = atoi(v[1].c_str());
248  phiMod = atoi(v[2].c_str());
249  etaMod = atoi(v[3].c_str());
250 
251  if (verboseLevel>5){
252  G4cout << "Volume name " << volName <<G4endl;
253  double xpos = coord1.x();
254  double ypos = coord1.y();
255  double zpos = coord1.z();
256  double r = sqrt(xpos*xpos+ypos*ypos);
257  G4cout << "In the Alpine " << G4endl;
258  G4cout << "----- PhiModule # " << phiMod << G4endl;
259  G4cout << "----- Ring/Eta # " << etaMod << G4endl;
260  G4cout << "----- Disk # " << LayerDisk << G4endl;
261  G4cout << "----- Endcap # " << BrlEcap << G4endl;
262  G4cout << "----- Pos # " << r<<" "<<zpos << G4endl;
263  G4cout << "----- volume " << myTouch->GetVolume()->GetName()<< G4endl;
264  G4cout << " " << myTouch->GetVolume(1)->GetName()<<" "<<
265  myTouch->GetVolume(2)->GetName()<< " " << myTouch->GetVolume(3)->GetName()<<G4endl;
266  }
267 
268  } else if(BEcopyNo == 600) { //ITk ECring
269 
270  std::string volName = myTouch->GetVolume()->GetName();
271 
272  // split the volname vs '_'
273  std::replace(volName.begin(),volName.end(),'_',' ');
274 
275  std::vector<std::string> v;
276  std::istringstream s(volName);
277  std::string tmp;
278  while (s >> tmp) {
279  v.push_back(tmp);
280  }
281 
282  BrlEcap = atoi(v[1].c_str());
283  LayerDisk = atoi(v[2].c_str());
284  phiMod = atoi(v[3].c_str());
285  etaMod = atoi(v[4].c_str());
286 
287  if (verboseLevel>5){
288  double xpos = coord1.x();
289  double ypos = coord1.y();
290  double zpos = coord1.z();
291  double r = sqrt(xpos*xpos+ypos*ypos);
292  G4cout << "Volume name " << volName <<G4endl;
293  G4cout << "In the ITk EC ring " << G4endl;
294  G4cout << "----- PhiModule # " << phiMod << G4endl;
295  G4cout << "----- Ring/Eta # " << etaMod << G4endl;
296  G4cout << "----- Disk # " << LayerDisk << G4endl;
297  G4cout << "----- Endcap # " << BrlEcap << G4endl;
298  G4cout << "----- Pos # " << r<<" "<<zpos << G4endl;
299  G4cout << "----- volume " << myTouch->GetVolume()->GetName()<< G4endl;
300  G4cout << " " << myTouch->GetVolume(1)->GetName()<<" "<<
301  myTouch->GetVolume(2)->GetName()<< " " << myTouch->GetVolume(3)->GetName()<<G4endl;
302  }
303  }
304  else {
305  // Do not expect other numbers. Need to fix PixelGeoModel if this occurs.
306  G4ExceptionDescription description;
307  description << "ProcessHits: Unrecognized geometry in Pixel sensitive detector. Please contact the maintainer of the Pixel Detector Description.";
308  description << "If you are processing a GeoModelXML geometry, please ensure a PixelSensorGmxSD instance is used instead of this (PixelSensorSD).";
309  G4Exception("PixelSensorSD", "UnrecognizedPixelGeometry", FatalException, description);
310  abort();
311  }
312 
313  // get the HepMcParticleLink from the TrackHelper
314  TrackHelper trHelp(aStep->GetTrack());
315  m_HitColl->Emplace(lP1,
316  lP2,
317  edep,
318  aStep->GetPreStepPoint()->GetGlobalTime(),
319  trHelp.GenerateParticleLink(),
320  0,BrlEcap,LayerDisk,etaMod,phiMod,side);
321  return true;
322 }
replace
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition: hcg.cxx:307
beamspotman.r
def r
Definition: beamspotman.py:674
AtlasG4EventUserInfo
This class is attached to G4Event objects as UserInformation. It holds a pointer to the HepMC::GenEve...
Definition: AtlasG4EventUserInfo.h:23
SiHit::xDep
@ xDep
Definition: SiHit.h:162
python.SystemOfUnits.mm
float mm
Definition: SystemOfUnits.py:98
TrackHelper.h
SiHit::xPhi
@ xPhi
Definition: SiHit.h:162
AtlasHitsVector< SiHit >
PixelSensorSD::PixelSensorSD
PixelSensorSD(const std::string &name, const std::string &hitCollectionName)
Definition: PixelSensorSD.cxx:32
PixelSensorSD::Initialize
void Initialize(G4HCofThisEvent *) override final
Definition: PixelSensorSD.cxx:40
PixelSensorSD::m_HitCollName
std::string m_HitCollName
Definition: PixelSensorSD.h:47
python.SystemOfUnits.MeV
float MeV
Definition: SystemOfUnits.py:172
TRT::Hit::side
@ side
Definition: HitInfo.h:83
PixelSensorSD::m_HitColl
SiHitCollection * m_HitColl
Definition: PixelSensorSD.h:48
AtlasHitsVector::Emplace
void Emplace(Args &&... args)
Definition: AtlasHitsVector.h:80
TrackHelper
Definition: TrackHelper.h:14
LayerDisk
Definition: WaferTree.h:57
lumiFormat.i
int i
Definition: lumiFormat.py:85
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
PixelSensorSD.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
SiHit::xEta
@ xEta
Definition: SiHit.h:162
python.PyAthena.v
v
Definition: PyAthena.py:154
PixelSensorSD::ProcessHits
G4bool ProcessHits(G4Step *, G4TouchableHistory *) override final
Definition: PixelSensorSD.cxx:52
AtlasG4EventUserInfo.h
python.SystemOfUnits.s
float s
Definition: SystemOfUnits.py:147
CxxUtils::atoi
int atoi(std::string_view str)
Helper functions to unpack numbers decoded in string into integers and doubles The strings are requir...
Definition: Control/CxxUtils/Root/StringUtils.cxx:85
TrackHelper::GenerateParticleLink
HepMcParticleLink GenerateParticleLink()
Generates a creates new HepMcParticleLink object on the stack based on GetUniqueID(),...
Definition: TrackHelper.h:35
fillSCTHists.etaMod
etaMod
Definition: fillSCTHists.py:23
description
std::string description
glabal timer - how long have I taken so far?
Definition: hcg.cxx:88