ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4mRICHSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4mRICHSteppingAction.cc
1 /*===============================================================*
2  * March 19th 2017 *
3  mRICH Stepping Action created by Cheuk-Ping Wong @GSU *
4  *===============================================================*/
6 
7 #include "PHG4mRICHDetector.h"
8 
9 #include <phparameter/PHParameters.h>
10 
11 #include <fun4all/Fun4AllBase.h>
12 
13 #include <g4main/PHG4Hit.h>
15 #include <g4main/PHG4Hitv1.h>
16 #include <g4main/PHG4Shower.h>
17 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
19 
20 #include <TSystem.h>
21 
22 #include <phool/getClass.h>
23 #include <phool/phool.h> // for PHWHERE
24 
25 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
26 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
27 #include <Geant4/G4Step.hh>
28 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
29 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
30 #include <Geant4/G4String.hh> // for G4String
31 #include <Geant4/G4SystemOfUnits.hh>
32 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
33 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
34 #include <Geant4/G4Track.hh> // for G4Track
35 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
36 #include <Geant4/G4Types.hh> // for G4double
37 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
38 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
39 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
40 
41 #include <boost/tokenizer.hpp>
42 // this is an ugly hack, the gcc optimizer has a bug which
43 // triggers the uninitialized variable warning which
44 // stops compilation because of our -Werror
45 #include <boost/version.hpp> // to get BOOST_VERSION
46 #if (__GNUC__ == 4 && __GNUC_MINOR__ == 4 && BOOST_VERSION == 105700)
47 #pragma GCC diagnostic ignored "-Wuninitialized"
48 #pragma message "ignoring bogus gcc warning in boost header lexical_cast.hpp"
49 #include <boost/lexical_cast.hpp>
50 #pragma GCC diagnostic warning "-Wuninitialized"
51 #else
52 #include <boost/lexical_cast.hpp>
53 #endif
54 
55 #include <cmath> // for isfinite
56 #include <iostream>
57 
58 class PHCompositeNode;
59 
60 using namespace CLHEP;
61 
62 //____________________________________________________________________________..
64  : PHG4SteppingAction(detector->GetName())
65  , detector_(detector)
66  ,
67  // active(params->get_int_param("active")),
68  IsBlackHole(params->get_int_param("blackhole"))
69  ,
70  // use_g4_steps(params->get_int_param("use_g4steps")),
71  detectorname(params->get_string_param("detectorname"))
72  , superdetector(params->get_string_param("superdetector"))
73  , hits_(nullptr)
74  , absorberhits_(nullptr)
75  , hit(nullptr)
76  , savehitcontainer(nullptr)
77  , saveshower(nullptr)
78  , savetrackid(-1)
79  , savepoststepstatus(-1)
80 {
81 }
82 //____________________________________________________________________________..
84 {
85  delete hit;
86 }
87 //____________________________________________________________________________..
89 {
92 
93  int isactive = detector_->IsInmRICH(volume);
94  if (isactive > PHG4mRICHDetector::INACTIVE)
95  {
96  // collect energy and track length step by step
98  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
99 
100  const G4Track* aTrack = aStep->GetTrack();
101 
102  /* Check if particle is 'geantino' */
103  bool geantino = false;
104  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
105  {
106  geantino = true;
107  }
108 
109  //cout << "Name of volume: " << volume->GetName() << ", isactive = " << isactive << endl;
110  // G4VPhysicalVolume* v1 = touch->GetVolume(1);
111  // cout << "Name of mother volume: " << v1->GetName() << endl;
112  // G4VPhysicalVolume* v2 = touch->GetVolume(2);
113  // cout << "Name of grand mother volume: " << v2->GetName() << endl;
114 
115  // cout << "copyNum = " << touch->GetReplicaNumber() << ", Or = " << touch->GetReplicaNumber(1) << ", Or = " << touch->GetReplicaNumber(2) << ", Or = " << touch->GetReplicaNumber(3) << endl;
116  // int module_id=GetModuleID(touch->GetVolume(2)); // use mother volume to determine module_id
117  int module_id = touch->GetReplicaNumber(2) - 1; // use copy number of mother volume to determine module_id
118  int PID = aTrack->GetDefinition()->GetPDGEncoding();
119  std::string PName = aTrack->GetDefinition()->GetParticleName();
120 
121  //-----------------------------------------------------------------------------------//
122  // if this block stops everything, just put all kinetic energy into edep
123  if (IsBlackHole)
124  {
125  edep = aTrack->GetKineticEnergy() / GeV;
126  G4Track* killtrack = const_cast<G4Track*>(aTrack);
127  killtrack->SetTrackStatus(fStopAndKill);
128  }
129 
130  /* Get Geant4 pre- and post-step points */
131  G4StepPoint* prePoint = aStep->GetPreStepPoint();
132  G4StepPoint* postPoint = aStep->GetPostStepPoint();
133 
134  switch (prePoint->GetStepStatus())
135  {
136  //-----------------
137  case fGeomBoundary:
138  //-----------------
139  case fUndefined:
140  if (!hit) hit = new PHG4Hitv1();
141  /* Set hit location (space point) */
142  hit->set_x(0, prePoint->GetPosition().x() / cm);
143  hit->set_y(0, prePoint->GetPosition().y() / cm);
144  hit->set_z(0, prePoint->GetPosition().z() / cm);
145 
146  /* Set hit time */
147  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
148 
149  //set the track ID
150  hit->set_trkid(aTrack->GetTrackID());
151  savetrackid = aTrack->GetTrackID();
152 
153  /* set intial energy deposit */
154  hit->set_edep(0);
155  hit->set_eion(0);
156  if (isactive == PHG4mRICHDetector::SENSOR)
157  {
159  if (PID == 0)
160  {
161  edep = aTrack->GetKineticEnergy() / GeV;
162  G4Track* killtrack = const_cast<G4Track*>(aTrack);
163  killtrack->SetTrackStatus(fStopAndKill);
164  hit->set_scint_id(module_id);
165  }
166  }
167  else
168  {
170  }
171 
172  // here we set what is common for scintillator and absorber hits
174  {
175  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
176  {
177  hit->set_trkid(pp->GetUserTrackId());
178  hit->set_shower_id(pp->GetShower()->get_id());
179  saveshower = pp->GetShower();
180  }
181  }
182  break;
183  //-----------------
184  default:
185  break;
186  }
187 
188  // some sanity checks for inconsistencies
189  // check if this hit was created, if not print out last post step status
190  if (!hit || !isfinite(hit->get_x(0)))
191  {
192  std::cout << __FILE__ << "::" << __func__ << "::" << GetName() << ": hit was not created" << std::endl;
193  std::cout << "prestep status: " << prePoint->GetStepStatus()
194  << ", last post step status: " << savepoststepstatus << std::endl;
195  gSystem->Exit(1);
196  }
197  savepoststepstatus = postPoint->GetStepStatus();
198 
199  // check if track id matches the initial one when the hit was created
200  if (aTrack->GetTrackID() != savetrackid)
201  {
202  std::cout << __FILE__ << "::" << __func__ << "::" << GetName() << ": hits do not belong to the same track" << std::endl;
203  std::cout << "saved track: " << savetrackid
204  << ", current trackid: " << aTrack->GetTrackID()
205  << std::endl;
206  gSystem->Exit(1);
207  }
208 
209  //-----------------------------------------------------------------------------------//
210  /* Update exit values- will be overwritten with every step until
211  * we leave the volume or the particle ceases to exist */
212  //-m/s- hit->set_x( 1, postPoint->GetPosition().x() / cm );
213  //-m/s- hit->set_y( 1, postPoint->GetPosition().y() / cm );
214  //-m/s- hit->set_z( 1, postPoint->GetPosition().z() / cm );
215 
216  //-m/s- hit->set_t( 1, postPoint->GetGlobalTime() / nanosecond );
217 
218  /* sum up the energy to get total deposited */
219  hit->set_edep(hit->get_edep() + edep);
220  hit->set_eion(hit->get_eion() + eion);
221 
222  if (geantino)
223  {
224  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
225  hit->set_eion(-1);
226  }
227  if (hit->get_edep() && PID == 0)
228  {
230  {
231  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
232  {
233  pp->SetKeep(1); // we want to keep the track
234  }
235  }
236  }
237 
238  //-----------------------------------------------------------------------------------//
239  // if any of these conditions is true this is the last step in
240  // this volume and we need to save the hit
241  if (postPoint->GetStepStatus() == fGeomBoundary ||
242  postPoint->GetStepStatus() == fWorldBoundary ||
243  postPoint->GetStepStatus() == fAtRestDoItProc ||
244  aTrack->GetTrackStatus() == fStopAndKill)
245  {
246  // save only hits with energy deposit (or -1 for geantino)
247  if (hit->get_edep() && PID == 0)
248  {
249  savehitcontainer->AddHit(module_id, hit);
250  if (saveshower)
251  {
253  }
254  // ownership has been transferred to container, set to null
255  // so we will create a new hit for the next track
256  hit = nullptr;
257  }
258  else
259  {
260  // if this hit has no energy deposit, just reset it for reuse
261  // this means we have to delete it in the dtor. If this was
262  // the last hit we processed the memory is still allocated
263  hit->Reset();
264  }
265  }
266  return true;
267  }
268  else
269  {
270  return false;
271  }
272 }
273 
274 //____________________________________________________________________________..
276 {
277  std::string hitnodename;
278  std::string absorbernodename;
279 
280  if (superdetector != "NONE")
281  {
282  hitnodename = "G4HIT_" + superdetector;
283  absorbernodename = "G4HIT_ABSORBER_" + superdetector;
284  }
285  else
286  {
287  hitnodename = "G4HIT_" + detectorname;
288  absorbernodename = "G4HIT_ABSORBER_" + detectorname;
289  }
290 
291  //now look for the map and grab a pointer to it.
292  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
293  absorberhits_ = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
294 
295  // if we do not find the node it's messed up.
296  if (!hits_)
297  {
298  if (Verbosity() >= Fun4AllBase::VERBOSITY_SOME) std::cout << __FILE__ << "::" << __func__ << " - unable to find " << hitnodename << std::endl;
299  }
300  if (!absorberhits_)
301  {
302  if (Verbosity() >= Fun4AllBase::VERBOSITY_SOME) std::cout << __FILE__ << "::" << __func__ << " - unable to find " << absorbernodename << std::endl;
303  }
304 }
305 //____________________________________________________________________________..
307 {
308  // G4AssemblyVolumes naming convention:
309  // av_WWW_impr_XXX_YYY_ZZZ
310  // where:
311  // WWW - assembly volume instance number
312  // XXX - assembly volume imprint number
313  // YYY - the name of the placed logical volume
314  // ZZZ - the logical volume index inside the assembly volume
315  // e.g. av_1_impr_4_mRICH_module_pv_11
316  // 4 is the number of the mRICH sector volume
317  // mRICH_module: name of mRICH module
318  // 11: number of mRICH module logical volume
319  // use boost tokenizer to separate the _, then take value
320  // after "impr" for mRICH sector volume and after "pv" for mRICH module volume
321  // use boost lexical cast for string -> int conversion
322  int module_id = -1;
323  int sector_id = -1;
324  int mRICH_id = -1;
325 
326  boost::char_separator<char> sep("_");
327  boost::tokenizer<boost::char_separator<char>> tok(volume->GetName(), sep);
328  boost::tokenizer<boost::char_separator<char>>::const_iterator tokeniter;
329 
330  for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
331  {
332  if (*tokeniter == "impr")
333  {
334  ++tokeniter;
335  if (tokeniter != tok.end())
336  {
337  sector_id = boost::lexical_cast<int>(*tokeniter);
338  }
339  else
340  {
341  std::cout << PHWHERE << " Error parsing " << volume->GetName()
342  << " for mRICH sector id " << std::endl;
343  gSystem->Exit(1);
344  }
345  }
346  else if (*tokeniter == "pv")
347  {
348  ++tokeniter;
349  if (tokeniter != tok.end())
350  {
351  mRICH_id = boost::lexical_cast<int>(*tokeniter);
352  }
353  else
354  {
355  std::cout << PHWHERE << " Error parsing " << volume->GetName()
356  << " for mRICH id " << std::endl;
357  gSystem->Exit(1);
358  }
359  }
360  }
361 
362  module_id = (sector_id - 1) * 100 + mRICH_id;
363 
364  if (Verbosity() >= Fun4AllBase::VERBOSITY_SOME) std::cout << "name of volume: " << volume->GetName() << ", sector_id = " << sector_id << ", mRICH_id = " << mRICH_id << ", module_id = " << module_id << std::endl;
365  //cout << endl;
366 
367  return module_id;
368 }