Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4SpacalPrototypeSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4SpacalPrototypeSteppingAction.cc
1 
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, fUndefined
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 <iostream>
40 #include <string> // for operator+, operator<<
41 
42 class G4VPhysicalVolume;
43 class PHCompositeNode;
44 
45 using namespace std;
46 //____________________________________________________________________________..
48  PHG4SteppingAction(detector->GetName()),
49  detector_(detector),
50  hits_(nullptr),
51  absorberhits_(nullptr),
52  hit(nullptr),
53  savehitcontainer(nullptr),
54  saveshower(nullptr)
55 {}
56 
57 //____________________________________________________________________________..
58 bool
60 {
61 
62  // get volume of the current step
63  G4VPhysicalVolume* volume =
64  aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
65 
66  // collect energy and track length step by step
67  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
68  G4double eion = (aStep->GetTotalEnergyDeposit()
69  - aStep->GetNonIonizingEnergyDeposit()) / GeV;
70 
71  const G4Track* aTrack = aStep->GetTrack();
72 
73  const int layer_id = 0;
74  // make sure we are in a volume
75  // IsInCylinderActive returns the number of the scintillator
76  // slat which has fired
77  int isactive = detector_->IsInCylinderActive(volume);
79  {
80  bool geantino = false;
81  // the check for the pdg code speeds things up, I do not want to make
82  // an expensive string compare for every track when we know
83  // geantino or chargedgeantino has pid=0
84  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0
85  && aTrack->GetParticleDefinition()->GetParticleName().find("geantino")
86  != string::npos)
87  {
88  geantino = true;
89  }
90  G4StepPoint * prePoint = aStep->GetPreStepPoint();
91  G4StepPoint * postPoint = aStep->GetPostStepPoint();
92  int scint_id = -1;
93 
94  if (//
96  or //
98  )
99  {
100  //SPACAL ID that is associated with towers
101  int sector_ID =0;
102  int tower_ID = 0;
103  int fiber_ID = 0;
104 
106  {
107 
108  fiber_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
109  tower_ID = prePoint->GetTouchable()->GetReplicaNumber(2);
110  sector_ID = prePoint->GetTouchable()->GetReplicaNumber(3);
111 
112  }
113 
114  else if (isactive == PHG4SpacalPrototypeDetector::FIBER_CLADING)
115  {
116  fiber_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
117  tower_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
118  sector_ID = prePoint->GetTouchable()->GetReplicaNumber(2);
119  }
120 
121  else if (isactive == PHG4SpacalPrototypeDetector::ABSORBER)
122  {
123  tower_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
124  sector_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
125  }
126 
127  // compact the tower/sector/fiber ID into 32 bit scint_id, so we could save some space for SPACAL hits
128  scint_id = PHG4CylinderGeom_Spacalv3::scint_id_coder(sector_ID, tower_ID, fiber_ID).scint_ID;
129 
130  }
131  else
132  {
133  // other configuraitons
135  {
136  scint_id = prePoint->GetTouchable()->GetReplicaNumber(2);
137  }
138  else if (isactive == PHG4SpacalPrototypeDetector::FIBER_CLADING)
139  {
140  scint_id = prePoint->GetTouchable()->GetReplicaNumber(1);
141  }
142  else
143  {
144  scint_id = prePoint->GetTouchable()->GetReplicaNumber(0);
145  }
146  }
147 
148  // cout << "track id " << aTrack->GetTrackID() << endl;
149  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
150  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
151  switch (prePoint->GetStepStatus())
152  {
153  case fGeomBoundary:
154  case fUndefined:
155  // case fPostStepDoItProc: // from point like interaction (compton,decay,hard rad)
156 
157  if (! hit)
158  {
159  hit = new PHG4Hitv1();
160  }
161  hit->set_layer((unsigned int) layer_id);
162  hit->set_scint_id(scint_id); // isactive contains the scintillator slat id
163  //here we set the entrance values in cm
164  hit->set_x(0, prePoint->GetPosition().x() / cm);
165  hit->set_y(0, prePoint->GetPosition().y() / cm);
166  hit->set_z(0, prePoint->GetPosition().z() / cm);
167 
168  // time in ns
169  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
170 
171  if (isactive == PHG4SpacalPrototypeDetector::FIBER_CORE) // only for active areas
172  {
173  // store all pre local coordinates
174  StoreLocalCoordinate(hit, aStep, true, false);
175  }
176 
177  //set the track ID
178  hit->set_trkid(aTrack->GetTrackID());
179  if ( G4VUserTrackInformation* p = aTrack->GetUserInformation() )
180  {
181  if ( PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p) )
182  {
183  hit->set_trkid(pp->GetUserTrackId());
184  hit->set_shower_id(pp->GetShower()->get_id());
185  }
186  }
187  //set the initial energy deposit
188  hit->set_edep(0);
189  if (isactive == PHG4SpacalPrototypeDetector::FIBER_CORE) // only for active areas
190  {
191  hit->set_eion(0); // only for fiber core hits
192  hit->set_light_yield(0);
194  }
195  else
196  {
198  }
199 
200  if ( G4VUserTrackInformation* p = aTrack->GetUserInformation() )
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 << "PHG4SpacalPrototypeSteppingAction: hit outside acceptance, layer: "
213 // << layer_id << endl;
214 // hit->identify();
215 // }
216  break;
217  default:
218  break;
219  }
220  // here we just update the exit values, it will be overwritten
221  // for every step until we leave the volume or the particle
222  // ceases to exist
223  hit->set_x(1, postPoint->GetPosition().x() / cm);
224  hit->set_y(1, postPoint->GetPosition().y() / cm);
225  hit->set_z(1, postPoint->GetPosition().z() / cm);
226 
227  hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
228 
229 
230  //sum up the energy to get total deposited
231  hit->set_edep(hit->get_edep() + edep);
232 
233  if (isactive == PHG4SpacalPrototypeDetector::FIBER_CORE) // only for active areas
234  {
235  // store all pre local coordinates
236  StoreLocalCoordinate(hit, aStep, false, true);
237  hit->set_eion(hit->get_eion() + eion);
238 
239  double light_yield = GetVisibleEnergyDeposition(aStep);
240 
241  static bool once = true;
242  if (once and edep > 0)
243  {
244  once = false;
245 
246  if (Verbosity() > 0) {
247  cout << "PHG4SpacalPrototypeSteppingAction::UserSteppingAction::"
248  //
249  << detector_->GetName() << " - "
250  << " use scintillating light model at each Geant4 steps. "
251  << "First step: " << "Material = "
252  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
253  << ", " << "Birk Constant = "
254  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetIonisation()->GetBirksConstant()
255  << "," << "edep = " << edep << ", " << "eion = " << eion
256  << ", " << "light_yield = " << light_yield << endl;
257  }
258  }
259 
260  hit->set_light_yield(hit->get_light_yield() + light_yield);
261  }
262 
263 // if (hit->get_z(1) > get_zmax() || hit->get_z(1) < get_zmin())
264 // {
265 // cout << "PHG4SpacalPrototypeSteppingAction: hit outside acceptance get_zmin() "
266 // << get_zmin() << ", get_zmax() " << get_zmax() << " at exit"
267 // << endl;
268 // hit->identify();
269 // }
270  if (geantino)
271  {
272  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
273 // hit->set_eion(-1);
274  }
275  if (edep > 0)
276  {
277  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
278  {
279  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
280  {
281  pp->SetKeep(1); // we want to keep the track
282  }
283  }
284  }
285 
286  // if any of these conditions is true this is the last step in
287  // this volume and we need to save the hit
288  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
289  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
290  // (not sure if this will ever be the case)
291  // aTrack->GetTrackStatus() == fStopAndKill: track ends
292  if (postPoint->GetStepStatus() == fGeomBoundary ||
293  postPoint->GetStepStatus() == fWorldBoundary||
294  aTrack->GetTrackStatus() == fStopAndKill)
295  {
296  // save only hits with energy deposit (or -1 for geantino)
297  if (hit->get_edep())
298  {
299  savehitcontainer->AddHit(layer_id, hit);
300  if (saveshower)
301  {
302  // save shower under container id of active volume (for running without absorber hit container to save space)
304  }
305  // ownership has been transferred to container, set to null
306  // so we will create a new hit for the next track
307  hit = nullptr;
308  }
309  else
310  {
311  // if this hit has no energy deposit, just reset it for reuse
312  // this means we have to delete it in the dtor. If this was
313  // the last hit we processed the memory is still allocated
314  hit->Reset();
315  }
316  }
317  // return true to indicate the hit was used
318  return true;
319 
320  }
321  else
322  {
323  return false;
324  }
325 }
326 
327 //____________________________________________________________________________..
328 void
330 {
331 
332  string hitnodename;
333  string absorbernodename;
334  if (detector_->SuperDetector() != "NONE")
335  {
336  hitnodename = "G4HIT_" + detector_->SuperDetector();
337  absorbernodename = "G4HIT_ABSORBER_" + detector_->SuperDetector();
338  }
339  else
340  {
341  hitnodename = "G4HIT_" + detector_->GetName();
342  absorbernodename = "G4HIT_ABSORBER_" + detector_->GetName();
343  }
344 
345  //now look for the map and grab a pointer to it.
346  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
347  absorberhits_ = findNode::getClass<PHG4HitContainer>(topNode,
348  absorbernodename.c_str());
349  // if we do not find the node we need to make it.
350  if (!hits_)
351  {
352  std::cout << "PHG4SpacalPrototypeSteppingAction::SetTopNode - unable to find "
353  << hitnodename << std::endl;
354  }
355  if (!absorberhits_)
356  {
357  if (Verbosity() > 1)
358  {
359  std::cout << "PHG4SpacalPrototypeSteppingAction::SetTopNode - unable to find "
360  << absorbernodename << std::endl;
361  }
362  }
363 }
364 
365 double
367 {
368  if (!detector_)
369  return 0;
370  else
371  return detector_->get_geom()->get_zmin() - .0001;
372 }
373 
374 double
376 {
377  if (!detector_)
378  return 0;
379  else
380  return detector_->get_geom()->get_zmax() + .0001;
381 }