ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ForwardDualReadoutSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ForwardDualReadoutSteppingAction.cc
3 
4 #include <g4main/PHG4Hit.h>
6 #include <g4main/PHG4Hitv1.h>
7 #include <g4main/PHG4Shower.h>
8 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
9 
11 
12 #include <phool/getClass.h>
13 
14 #include <Geant4/G4IonisParamMat.hh> // for G4IonisParamMat
15 #include <Geant4/G4Material.hh> // for G4Material
16 #include <Geant4/G4MaterialCutsCouple.hh>
17 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
18 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
19 #include <Geant4/G4Step.hh>
20 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
21 #include <Geant4/G4String.hh> // for G4String
22 #include <Geant4/G4SystemOfUnits.hh>
23 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
24 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
25 #include <Geant4/G4Track.hh> // for G4Track
26 #include <Geant4/G4VSolid.hh> // for G4Track
27 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
28 #include <Geant4/G4Types.hh> // for G4double
29 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
30 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
31 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
32 #include <Geant4/G4VSensitiveDetector.hh> // for G4VUserTrackInformation
33 #include <Geant4/G4OpticalPhoton.hh>
34 #include <Geant4/G4Scintillation.hh>
35 #include <Geant4/G4Cerenkov.hh>
36 
37 #include <boost/tokenizer.hpp>
38 // this is an ugly hack, the gcc optimizer has a bug which
39 // triggers the uninitialized variable warning which
40 // stops compilation because of our -Werror
41 #include <boost/version.hpp> // to get BOOST_VERSION
42 #if (__GNUC__ == 4 && __GNUC_MINOR__ == 4 && BOOST_VERSION == 105700)
43 #pragma GCC diagnostic ignored "-Wuninitialized"
44 #pragma message "ignoring bogus gcc warning in boost header lexical_cast.hpp"
45 #include <boost/lexical_cast.hpp>
46 #pragma GCC diagnostic warning "-Wuninitialized"
47 #else
48 #include <boost/lexical_cast.hpp>
49 #endif
50 
51 #include <iostream>
52 #include <string> // for basic_string, operator+
53 
54 class PHCompositeNode;
55 
56 using namespace std;
57 
58 //____________________________________________________________________________..
60  : PHG4SteppingAction(detector->GetName())
61  , detector_(detector)
62  , hits_(nullptr)
63  , absorberhits_(nullptr)
64  , hitcontainer(nullptr)
65  , hit(nullptr)
66  , saveshower(nullptr)
67  , _towerdivision(0.0)
68  , _tower_size(1.0)
69  , _readout_size(1.0)
70  , _detector_size(100)
71  , absorbertruth(absorberactive)
72  , light_scint_model(1)
73 {
74 }
75 
77 {
78  // if the last hit was a zero energie deposit hit, it is just reset
79  // and the memory is still allocated, so we need to delete it here
80  // if the last hit was saved, hit is a nullptr pointer which are
81  // legal to delete (it results in a no operation)
82  delete hit;
83 }
84 
85 //____________________________________________________________________________..
87 {
89  // G4TouchableHistory* theTouchable = (G4TouchableHistory*)( aStep->GetPreStepPoint()->GetTouchable() );
91  // G4LogicalVolume* volumeMother = volume->GetLogicalVolume();
92  // G4VSensitiveDetector* sensdet = volumeMother->GetSensitiveDetector();
93  // G4VSolid* volumeGMother = nullptr;
94  // if(volumeMother) volumeGMother = volumeMother->GetSolid();
95 
96  // detector_->IsInForwardDualReadout(volume)
97  // returns
98  // 0 is outside of Forward HCAL
99  // 1 is inside scintillator
100  // -1 is inside absorber (dead material)
101 
102  int whichactive = detector_->IsInForwardDualReadout(volume);
103 
104 
105  // if(whichactive){
106  // // // if(sensdet)
107  // // // cout << theTouchable->GetReplicaNumber() << "\t"<< theTouchable->GetReplicaNumber(1) << "\t"<< theTouchable->GetReplicaNumber(2) << "\t" << volume->GetName() << endl;
108  // // // cout << whichactive << "\t" << volume->GetName() << "\tSD:" << sensdet->GetName() << endl;
109  // // // else if(volumeGMother)
110  // // // cout << whichactive << "\t" << volume->GetName() << "\tGM:" << volumeGMother->GetName() << endl;
111  // // // else if(volumeMother)
112  // // // cout << whichactive << "\t" << volume->GetName() << "\tMM:" << volumeMother->GetName() << endl;
113  // // // else
114  // cout << whichactive << "\t" << volume->GetName() << endl;
115  // // // cout << whichactive << "\t" << touch->GetCopyNumber(2)<< "\t" << touch->GetCopyNumber(3) << endl;
116  // // // cout << "\tx: " << (aStep->GetPreStepPoint()->GetPosition()).x()<< "\ty: " << (aStep->GetPreStepPoint()->GetPosition()).y()<< "\tz: " << (aStep->GetPreStepPoint()->GetPosition()).z()<< endl;
117  // } else {
118  // cout << volume->GetCopyNo() << "\t" << volume->GetName() << endl;
119 
120  // }
121 
122 
123  if (!whichactive)
124  {
125  return false;
126  }
127 
128  int layer_id = detector_->get_Layer();
129  // unsigned int tower_id = -1;
130  int idx_j = -1;
131  int idx_k = -1;
132 
133  // if (whichactive > 0) // in sctintillator
134  // {
135  // /* Find indizes of sctintillator / tower containing this step */
136  // // FindTowerIndex(touch, idx_j, idx_k);
137  // tower_id = touch->GetVolume(1)->GetCopyNo();
138  // }
139  // else
140  // {
141  // tower_id = touch->GetVolume(1)->GetCopyNo();
142  // }
143 
144  /* Get energy deposited by this step */
145  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
146  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
147  G4double light_yield = 0;
148 
149  /* Get pointer to associated Geant4 track */
150  const G4Track* aTrack = aStep->GetTrack();
151 
152  // if this block stops everything, just put all kinetic energy into edep
153  if (detector_->IsBlackHole())
154  {
155  edep = aTrack->GetKineticEnergy() / GeV;
156  G4Track* killtrack = const_cast<G4Track*>(aTrack);
157  killtrack->SetTrackStatus(fStopAndKill);
158  }
159 
160  /* Make sure we are in a volume */
161  if (detector_->IsActive())
162  {
163  // int idx_l = -1;
164  /* Check if particle is 'geantino' */
165  bool geantino = false;
166  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
167  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
168  {
169  geantino = true;
170  }
171 
172  /* Get Geant4 pre- and post-step points */
173  G4StepPoint* prePoint = aStep->GetPreStepPoint();
174  G4StepPoint* postPoint = aStep->GetPostStepPoint();
175  if(abs( prePoint->GetPosition().y() / cm)>(_detector_size*1.1) || abs( prePoint->GetPosition().x() / cm)>(_detector_size*1.1)) return false;
176  // cout << "x: " << prePoint->GetPosition().x() << "\ty: " << prePoint->GetPosition().y() << "\tz: " << prePoint->GetPosition().z() << endl;
177  switch (prePoint->GetStepStatus())
178  {
179  case fGeomBoundary:
180  case fUndefined:
181  if (!hit)
182  {
183  hit = new PHG4Hitv1();
184  }
185  // hit->set_scint_id(tower_id);
186  FindTowerIndexFromPosition(prePoint, idx_j, idx_k);
187  // cout << idx_j << ":" << idx_k << endl;
188  /* Set hit location (tower index) */
189  hit->set_index_j(idx_j);
190  hit->set_index_k(idx_k);
191  // hit->set_index_l(idx_l);
192 
193 // cout << "\tidj: " << idx_j << "\tidk: " << idx_k << "\tz " << prePoint->GetPosition().z()/cm << endl;
194  // TODO use these positions for an index conversion
195  /* Set hit location (space point) */
196  hit->set_x(0, prePoint->GetPosition().x() / cm);
197  hit->set_y(0, prePoint->GetPosition().y() / cm);
198  hit->set_z(0, prePoint->GetPosition().z() / cm);
199 
200  /* Set hit time */
201  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
202 
203  //set the track ID
204  hit->set_trkid(aTrack->GetTrackID());
205  /* set intial energy deposit */
206  hit->set_edep(0);
207  hit->set_eion(0);
208  hit->set_property(PHG4Hit::PROPERTY::scint_gammas, (float)0.);
209  hit->set_property(PHG4Hit::PROPERTY::cerenkov_gammas, (float)0.);
210 
211  /* Now add the hit to the hit collection */
212  if (whichactive > 0)
213  {
215  hit->set_light_yield(0);
216  }
217  else
218  {
220  }
222  {
223  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
224  {
225  hit->set_trkid(pp->GetUserTrackId());
226  hit->set_shower_id(pp->GetShower()->get_id());
227  saveshower = pp->GetShower();
228  }
229  }
230  break;
231  default:
232  break;
233  }
234 
235  if (whichactive > 0)
236  {
237  if (light_scint_model)
238  {
239  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
240  static bool once = true;
241  if (once && edep > 0)
242  {
243  once = false;
244 
245  if (Verbosity() > 0)
246  {
247  cout << "PHG4ForwardDualReadoutSteppingAction::UserSteppingAction::"
248  //
249  << detector_->GetName() << " - "
250  << " use scintillating light model at each Geant4 steps. "
251  << "First step: "
252  << "Material = "
253  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
254  << ", "
255  << "Birk Constant = "
257  << ","
258  << "edep = " << edep << ", "
259  << "eion = " << eion
260  << ", "
261  << "light_yield = " << light_yield << endl;
262  }
263  }
264  }
265  else
266  {
267  light_yield = eion;
268  }
269  }
270 
271 
272 
273  // Double_t fE; //deposited energy in the cell
274  // ULong64_t fNphot = 0; // number of optical photons
275  float fNscin = 0; // scintillation photons
276  float fNcerenkov = 0; // Cerenkov photons
277 
278  G4int fScinType; // scintillation process type
279  G4int fScinSubType; // scintillation process subtype
280  G4int fCerenkovType; // Cerenkov process type
281  G4int fCerenkovSubType; // Cerenkov process subtype
282 
283  G4Scintillation scin;
284  fScinType = scin.GetProcessType();
285  fScinSubType = scin.GetProcessSubType();
286  G4Cerenkov cer;
287  fCerenkovType = cer.GetProcessType();
288  fCerenkovSubType = cer.GetProcessSubType();
289 
290  //number of optical photons in the event from secondary tracks
291  // const std::vector<const G4Track*> *sec = aStep->GetSecondaryInCurrentStep();
292  // std::vector<const G4Track*>::const_iterator ittr;
293  // for(ittr = sec->begin(); ittr != sec->end(); ittr++) {
294  // if((*ittr)->GetParentID() <= 0) continue;
295 
296  //all optical photons
298  {
299  // fNphot++;
300 
301  //identify the process
302  G4int ptype = aTrack->GetCreatorProcess()->GetProcessType();
303  G4int pstype = aTrack->GetCreatorProcess()->GetProcessSubType();
304  // cout << aTrack->GetCreatorProcess()->GetProcessName() << endl;
305  //scintillation photons
306  G4Material* prevMaterial = aStep->GetPreStepPoint()->GetMaterial();
307  if((ptype == fScinType) && (pstype == fScinSubType) && (prevMaterial->GetName().find("G4_POLYSTYRENE") != std::string::npos)){ fNscin++;}
308  //Cerenkov photons
309  if( (ptype == fCerenkovType) && (pstype == fCerenkovSubType) && ((prevMaterial->GetName().find("PMMA") != std::string::npos) || (prevMaterial->GetName().find("Quartz") != std::string::npos))){ fNcerenkov++;}
310  // if(aTrack->GetParentID() > 0)
311  // {
312  // if(aTrack->GetCreatorProcess()->GetProcessName().find("enkov") != string::npos)cout << aTrack->GetCreatorProcess()->GetProcessName() << endl;
313  // }
314 // if(fNscin>0 || fNcerenkov>0) cout << aTrack->GetCreatorProcess()->GetProcessName() << "\tfNscin: " << fNscin << "\tfNcerenkov: " << fNcerenkov << endl;
315  }//secondary tracks loop
316 // cout << __LINE__ << endl;
317 // cout << hit->get_property_float(PHG4Hit::PROPERTY::scint_gammas) << "\tadd fNscin: " << fNscin << "\t" << hit->get_property_float(PHG4Hit::PROPERTY::cerenkov_gammas) << "\t add fNcerenkov: " << fNcerenkov<< endl;
318  //G4Material* nextMaterial = aStep->GetPostStepPoint()->GetMaterial();
319  //string materialstr = nextMaterial->GetName();
320 // string materialstr2 = prevMaterial->GetName();
321 // if(materialstr.find("G4_AIR") == std::string::npos)cout << materialstr << endl;;
322 // if(materialstr2.find("G4_AIR") == std::string::npos)cout << "\t" << materialstr2 << endl;;
323  hit->set_property(PHG4Hit::PROPERTY::scint_gammas,hit->get_property_float(PHG4Hit::PROPERTY::scint_gammas)+ fNscin);
324  // cout << __LINE__ << endl;
325  hit->set_property(PHG4Hit::PROPERTY::cerenkov_gammas,hit->get_property_float(PHG4Hit::PROPERTY::cerenkov_gammas)+ fNcerenkov);
326 // cout << hit->get_property_float(PHG4Hit::PROPERTY::scint_gammas)<< "\t" << hit->get_property_float(PHG4Hit::PROPERTY::cerenkov_gammas) << endl;
327 // cout << __LINE__ << endl;
328 
329  // //number of optical photons in the event from secondary tracks
330  // const std::vector<const G4Track*> *sec = aStep->GetSecondaryInCurrentStep();
331  // std::vector<const G4Track*>::const_iterator ittr;
332  // for(ittr = sec->begin(); ittr != sec->end(); ittr++) {
333  // if((*ittr)->GetParentID() <= 0) continue;
334 
335  // //all optical photons
336  // if((*ittr)->GetDynamicParticle()->GetParticleDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) continue;
337  // fNphot++;
338 
339  // //identify the process
340  // G4int ptype = (*ittr)->GetCreatorProcess()->GetProcessType();
341  // G4int pstype = (*ittr)->GetCreatorProcess()->GetProcessSubType();
342  // //scintillation photons
343  // if(ptype == fScinType && pstype == fScinSubType) fNscin++;
344  // //Cerenkov photons
345  // if(ptype == fCerenkovType && pstype == fCerenkovSubType) fNcerenkov++;
346 
347  // }//secondary tracks loop
348  // if(sec->size()>0) cout << "fNphot: " << fNphot << "\tfNscin: " << fNscin << "\tfNcerenkov: " << fNcerenkov << endl;
349  // }
350 
351 
352 
353  /* Update exit values- will be overwritten with every step until
354  * we leave the volume or the particle ceases to exist */
355  hit->set_x(1, postPoint->GetPosition().x() / cm);
356  hit->set_y(1, postPoint->GetPosition().y() / cm);
357  hit->set_z(1, postPoint->GetPosition().z() / cm);
358 
359  hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
360 
361  /* sum up the energy to get total deposited */
362  hit->set_edep(hit->get_edep() + edep);
363  if (whichactive > 0)
364  {
365  hit->set_eion(hit->get_eion() + eion);
366  hit->set_light_yield(hit->get_light_yield() + light_yield);
367  }
368 
369  if (geantino)
370  {
371  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
372  if (whichactive > 0)
373  {
374  hit->set_eion(-1);
375  hit->set_light_yield(-1);
376  }
377  }
378  if (edep > 0 && (whichactive > 0 || absorbertruth > 0))
379  {
381  {
382  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
383  {
384  pp->SetKeep(1); // we want to keep the track
385  }
386  }
387  }
388  // if any of these conditions is true this is the last step in
389  // this volume and we need to save the hit
390  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
391  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
392  // (not sure if this will ever be the case)
393  // aTrack->GetTrackStatus() == fStopAndKill: track ends
394  if (postPoint->GetStepStatus() == fGeomBoundary ||
395  postPoint->GetStepStatus() == fWorldBoundary ||
396  postPoint->GetStepStatus() == fAtRestDoItProc ||
397  aTrack->GetTrackStatus() == fStopAndKill)
398  {
399  // save only hits with energy deposit (or -1 for geantino)
400  if (hit->get_edep())
401  {
402  hitcontainer->AddHit(layer_id, hit);
403  if (saveshower)
404  {
406  }
407  // ownership has been transferred to container, set to null
408  // so we will create a new hit for the next track
409  hit = nullptr;
410  }
411  else
412  {
413  // if this hit has no energy deposit, just reset it for reuse
414  // this means we have to delete it in the dtor. If this was
415  // the last hit we processed the memory is still allocated
416  hit->Reset();
417  }
418  }
419 
420  return true;
421  }
422  else
423  {
424  return false;
425  }
426 }
427 
428 //____________________________________________________________________________..
430 {
431  std::string hitnodename;
432  std::string absorbernodename;
433 
434  // std::cout << detector_->SuperDetector() << "\t" << detector_->GetName() << endl;
435  if (detector_->SuperDetector() != "NONE")
436  {
437  hitnodename = "G4HIT_" + detector_->SuperDetector();
438  absorbernodename = "G4HIT_ABSORBER_" + detector_->SuperDetector();
439  }
440  else
441  {
442  hitnodename = "G4HIT_" + detector_->GetName();
443  absorbernodename = "G4HIT_ABSORBER_" + detector_->GetName();
444  }
445 
446  //now look for the map and grab a pointer to it.
447  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
448  absorberhits_ = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
449 
450  // if we do not find the node it's messed up.
451  if (!hits_)
452  {
453  std::cout << "PHG4ForwardDualReadoutSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
454  }
455  if (!absorberhits_)
456  {
457  if (Verbosity() > 0)
458  {
459  cout << "PHG4ForwardDualReadoutSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
460  }
461  }
462 }
463 
465 {
466  int j_0 = 0; //The j and k indices for the scintillator / tower
467  int k_0 = 0; //The j and k indices for the scintillator / tower
468 
469  // G4VPhysicalVolume* tower = touch->GetVolume(1); //Get the tower solid
470  // ParseG4VolumeName(tower, j_0, k_0);
471  if(_towerdivision==0.){
472  int maxsubtow = (int) ( (_tower_size) / (_readout_size));
473  _towerdivision = (_tower_size - (maxsubtow * _readout_size))/maxsubtow;
475  }
476  j_0 = (int) ( ( _detector_size + ( prePoint->GetPosition().x() ) ) / _towerdivision ); //TODO DRCALO TOWER SIZE
477  k_0 = (int) ( ( _detector_size + ( prePoint->GetPosition().y() ) ) / _towerdivision ); //TODO DRCALO TOWER SIZE
478  // if(prePoint->GetPosition().x()>290) cout << prePoint->GetPosition().x() << "\t" << k_0 << endl;
479  j = (j_0 * 1);
480  k = (k_0 * 1);
481 
482  return 0;
483 }
484 
486 {
487  int j_0, k_0; //The j and k indices for the scintillator / tower
488 
489  G4VPhysicalVolume* tower = touch->GetVolume(1); //Get the tower solid
490  ParseG4VolumeName(tower, j_0, k_0);
491 
492  j = (j_0 * 1);
493  k = (k_0 * 1);
494 
495  return 0;
496 }
497 
499 {
500  // cout << volume->GetName() << endl;
501  boost::char_separator<char> sep("_");
502  boost::tokenizer<boost::char_separator<char> > tok(volume->GetName(), sep);
503  boost::tokenizer<boost::char_separator<char> >::const_iterator tokeniter;
504  for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
505  {
506  if (*tokeniter == "j")
507  {
508  ++tokeniter;
509  if (tokeniter == tok.end()) break;
510  j = boost::lexical_cast<int>(*tokeniter);
511  }
512  else if (*tokeniter == "k")
513  {
514  ++tokeniter;
515  if (tokeniter == tok.end()) break;
516  k = boost::lexical_cast<int>(*tokeniter);
517  }
518  }
519  return 0;
520 }