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