ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4SpacalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4SpacalSteppingAction.cc
1 
11 #include "PHG4SpacalDetector.h"
12 
14 #include <g4main/PHG4Hit.h> // for PHG4Hit
15 #include <g4main/PHG4Hitv1.h>
16 #include <g4main/PHG4Shower.h>
17 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
19 
20 #include <phool/getClass.h>
21 
22 #include <Geant4/G4IonisParamMat.hh> // for G4IonisParamMat
23 #include <Geant4/G4Material.hh> // for G4Material
24 #include <Geant4/G4MaterialCutsCouple.hh>
25 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
26 #include <Geant4/G4Step.hh>
27 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
28 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRestD...
29 #include <Geant4/G4String.hh> // for G4String
30 #include <Geant4/G4SystemOfUnits.hh>
31 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
32 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
33 #include <Geant4/G4Track.hh> // for G4Track
34 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
35 #include <Geant4/G4Types.hh> // for G4double
36 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
37 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
38 
39 #include <cmath> // for isfinite
40 #include <cstdlib> // for exit
41 #include <iostream>
42 #include <string> // for operator<<, char_traits
43 
44 class G4VPhysicalVolume;
45 class PHCompositeNode;
46 
47 using namespace std;
48 //____________________________________________________________________________..
50  detector_(detector),
51  hits_(nullptr),
52  absorberhits_(nullptr),
53  hit(nullptr),
54  savehitcontainer(nullptr),
55  saveshower(nullptr),
56  savetrackid(-1),
57  savepoststepstatus(-1)
58 {
59 }
60 
62 {
63  // if the last hit was a zero energie deposit hit, it is just reset
64  // and the memory is still allocated, so we need to delete it here
65  // if the last hit was saved, hit is a nullptr pointer which are
66  // legal to delete (it results in a no operation)
67  delete hit;
68 }
69 
70 //____________________________________________________________________________..
72 {
73  // get volume of the current step
76 
77  // collect energy and track length step by step
79  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
80 
81  const G4Track* aTrack = aStep->GetTrack();
82 
83  int layer_id = detector_->get_Layer();
84  // make sure we are in a volume
85  // IsInCylinderActive returns the number of the scintillator
86  // slat which has fired
87  int isactive = detector_->IsInCylinderActive(volume);
88  if (isactive > PHG4SpacalDetector::INACTIVE)
89  {
90  bool geantino = false;
91  // the check for the pdg code speeds things up, I do not want to make
92  // an expensive string compare for every track when we know
93  // geantino or chargedgeantino has pid=0
94  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
95  {
96  geantino = true;
97  }
98  G4StepPoint* prePoint = aStep->GetPreStepPoint();
99  G4StepPoint* postPoint = aStep->GetPostStepPoint();
100  int scint_id = -1;
101 
102  if ( //
104  or //
106  or //
108  or //
110  )
111  {
112  //SPACAL ID that is associated with towers
113  int sector_ID = 0;
114  int tower_ID = 0;
115  int fiber_ID = 0;
116 
117  if (isactive == PHG4SpacalDetector::FIBER_CORE)
118  {
119  fiber_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
120  tower_ID = prePoint->GetTouchable()->GetReplicaNumber(2);
121  sector_ID = prePoint->GetTouchable()->GetReplicaNumber(3);
122  }
123 
124  else if (isactive == PHG4SpacalDetector::FIBER_CLADING)
125  {
126  fiber_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
127  tower_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
128  sector_ID = prePoint->GetTouchable()->GetReplicaNumber(2);
129  }
130 
131  else if (isactive == PHG4SpacalDetector::ABSORBER)
132  {
133  tower_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
134  sector_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
135  }
136 
137  else if (isactive == PHG4SpacalDetector::SUPPORT)
138  {
139  tower_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
140  sector_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
141  fiber_ID = (1 << (PHG4CylinderGeom_Spacalv3::scint_id_coder::kfiber_bit)) - 1; // use max fiber ID to flag for support strucrtures.
142 
143 // cout <<"PHG4SpacalSteppingAction::UserSteppingAction - SUPPORT tower_ID = "<<tower_ID<<endl;
144  }
145 
146  // compact the tower/sector/fiber ID into 32 bit scint_id, so we could save some space for SPACAL hits
147  scint_id = PHG4CylinderGeom_Spacalv3::scint_id_coder(sector_ID, tower_ID, fiber_ID).scint_ID;
148  }
149  else
150  {
151  // other configuraitons
152  if (isactive == PHG4SpacalDetector::FIBER_CORE)
153  scint_id = prePoint->GetTouchable()->GetReplicaNumber(2);
154  else if (isactive == PHG4SpacalDetector::FIBER_CLADING)
155  scint_id = prePoint->GetTouchable()->GetReplicaNumber(1);
156  else
157  scint_id = prePoint->GetTouchable()->GetReplicaNumber(0);
158  }
159 
160  // cout << "track id " << aTrack->GetTrackID() << endl;
161  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
162  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
163  switch (prePoint->GetStepStatus())
164  {
165  case fGeomBoundary:
166  case fUndefined:
167  // if previous hit was saved, hit pointer was set to nullptr
168  // and we have to make a new one
169  if (!hit)
170  {
171  hit = new PHG4Hitv1();
172  }
173  hit->set_layer((unsigned int) layer_id);
174  hit->set_scint_id(scint_id); // isactive contains the scintillator slat id
175  //here we set the entrance values in cm
176  hit->set_x(0, prePoint->GetPosition().x() / cm);
177  hit->set_y(0, prePoint->GetPosition().y() / cm);
178  hit->set_z(0, prePoint->GetPosition().z() / cm);
179 
180  // time in ns
181  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
182  //set the track ID
183  hit->set_trkid(aTrack->GetTrackID());
184  savetrackid = aTrack->GetTrackID();
185  //set the initial energy deposit
186  hit->set_edep(0);
187  // Now add the hit
188  if (isactive == PHG4SpacalDetector::FIBER_CORE) // the slat ids start with zero
189  {
190  // store all pre local coordinates
191  StoreLocalCoordinate(hit, aStep, true, false);
192  hit->set_eion(0); // only implemented for v5 otherwise empty
193  hit->set_light_yield(0);
195  }
196  else
197  {
199  }
201  {
202  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
203  {
204  hit->set_trkid(pp->GetUserTrackId());
205  hit->set_shower_id(pp->GetShower()->get_id());
206  saveshower = pp->GetShower();
207  }
208  }
209 
210  if (hit->get_z(0) > get_zmax() || hit->get_z(0) < get_zmin())
211  {
212  cout << "PHG4SpacalSteppingAction: hit outside acceptance, layer: "
213  << layer_id << endl;
214  hit->identify();
215  }
216  break;
217  default:
218  break;
219  }
220  // some sanity checks for inconsistencies
221  // check if this hit was created, if not print out last post step status
222  if (!hit || !isfinite(hit->get_x(0)))
223  {
224  cout << GetName() << ": hit was not created" << endl;
225  cout << "prestep status: " << prePoint->GetStepStatus()
226  << ", last post step status: " << savepoststepstatus << endl;
227  exit(1);
228  }
229  savepoststepstatus = postPoint->GetStepStatus();
230  // check if track id matches the initial one when the hit was created
231  if (aTrack->GetTrackID() != savetrackid)
232  {
233  cout << GetName() << ": hits do not belong to the same track" << endl;
234  cout << "saved track: " << savetrackid
235  << ", current trackid: " << aTrack->GetTrackID()
236  << endl;
237  exit(1);
238  }
239  // here we just update the exit values, it will be overwritten
240  // for every step until we leave the volume or the particle
241  // ceases to exist
242  hit->set_x(1, postPoint->GetPosition().x() / cm);
243  hit->set_y(1, postPoint->GetPosition().y() / cm);
244  hit->set_z(1, postPoint->GetPosition().z() / cm);
245 
246  hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
247  //sum up the energy to get total deposited
248  hit->set_edep(hit->get_edep() + edep);
249 
250  if (isactive == PHG4SpacalDetector::FIBER_CORE) // only for active areas
251  {
252  // store all pre local coordinates
253  StoreLocalCoordinate(hit, aStep, false, true);
254 
255  hit->set_eion(hit->get_eion() + eion);
256 
257  double light_yield = GetVisibleEnergyDeposition(aStep);
258 
259  static bool once = true;
260  if (once and edep > 0)
261  {
262  once = false;
263 
264  if (Verbosity() > 0)
265  {
266  cout << "PHG4SpacalSteppingAction::UserSteppingAction::"
267  //
268  << detector_->GetName() << " - "
269  << " use scintillating light model at each Geant4 steps. "
270  << "First step: "
271  << "Material = "
272  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
273  << ", "
274  << "Birk Constant = "
276  << ","
277  << "edep = " << edep << ", "
278  << "eion = " << eion
279  << ", "
280  << "light_yield = " << light_yield << endl;
281  }
282  }
283 
284  hit->set_light_yield(hit->get_light_yield() + light_yield);
285  }
286 
287  if (hit->get_z(1) > get_zmax() || hit->get_z(1) < get_zmin())
288  {
289  cout << "PHG4SpacalSteppingAction: hit outside acceptance get_zmin() "
290  << get_zmin() << ", get_zmax() " << get_zmax() << " at exit"
291  << endl;
292  hit->identify();
293  }
294  if (geantino)
295  {
296  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
297  // hit->set_eion(-1);
298  }
299  if (edep > 0)
300  {
302  {
303  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
304  {
305  pp->SetKeep(1); // we want to keep the track
306  }
307  }
308  }
309  // if any of these conditions is true this is the last step in
310  // this volume and we need to save the hit
311  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
312  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
313  // (happens when your detector goes outside world volume)
314  // aTrack->GetTrackStatus() == fStopAndKill: track ends
315  if (postPoint->GetStepStatus() == fGeomBoundary ||
316  postPoint->GetStepStatus() == fWorldBoundary ||
317  postPoint->GetStepStatus() == fAtRestDoItProc ||
318  aTrack->GetTrackStatus() == fStopAndKill)
319  {
320  // save only hits with energy deposit (or -1 for geantino)
321  if (hit->get_edep())
322  {
323  savehitcontainer->AddHit(layer_id, hit);
324  if (saveshower)
325  {
327  }
328  // ownership has been transferred to container, set to null
329  // so we will create a new hit for the next track
330  hit = nullptr;
331  }
332  else
333  {
334  // if this hit has no energy deposit, just reset it for reuse
335  // this means we have to delete it in the dtor. If this was
336  // the last hit we processed the memory is still allocated
337  hit->Reset();
338  }
339  }
340  // return true to indicate the hit was used
341  return true;
342  }
343  else
344  {
345  return false;
346  }
347 }
348 
349 //____________________________________________________________________________..
351 {
352  string hitnodename;
353  string absorbernodename;
354  if (detector_->SuperDetector() != "NONE")
355  {
356  hitnodename = "G4HIT_" + detector_->SuperDetector();
357  absorbernodename = "G4HIT_ABSORBER_" + detector_->SuperDetector();
358  }
359  else
360  {
361  hitnodename = "G4HIT_" + detector_->GetName();
362  absorbernodename = "G4HIT_ABSORBER_" + detector_->GetName();
363  }
364 
365  //now look for the map and grab a pointer to it.
366  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
367  absorberhits_ = findNode::getClass<PHG4HitContainer>(topNode,
368  absorbernodename.c_str());
369  // if we do not find the node we need to make it.
370  if (!hits_)
371  {
372  std::cout << "PHG4SpacalSteppingAction::SetTopNode - unable to find "
373  << hitnodename << std::endl;
374  }
375  if (!absorberhits_)
376  {
377  if (Verbosity() > 1)
378  {
379  std::cout << "PHG4SpacalSteppingAction::SetTopNode - unable to find "
380  << absorbernodename << std::endl;
381  }
382  }
383 }
384 
385 double
387 {
388  if (!detector_)
389  return 0;
390  else
391  return detector_->get_geom()->get_zmin() - .0001;
392 }
393 
394 double
396 {
397  if (!detector_)
398  return 0;
399  else
400  return detector_->get_geom()->get_zmax() + .0001;
401 }