Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4InnerHcalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4InnerHcalSteppingAction.cc
2 
4 #include "PHG4StepStatusDecode.h"
5 
6 #include <phparameter/PHParameters.h>
7 
8 #include <g4main/PHG4Hit.h>
10 #include <g4main/PHG4Hitv1.h>
11 #include <g4main/PHG4Shower.h>
12 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
14 
16 
18 
19 #include <calobase/TowerInfo.h>
20 #include <calobase/TowerInfoContainer.h>
21 #include <calobase/TowerInfoContainerv1.h>
22 #include <calobase/TowerInfoDefs.h>
23 
24 #include <phool/PHCompositeNode.h>
25 #include <phool/PHIODataNode.h> // for PHIODataNode
26 #include <phool/PHNode.h> // for PHNode
27 #include <phool/PHNodeIterator.h> // for PHNodeIterator
28 #include <phool/PHObject.h> // for PHObject
29 #include <phool/getClass.h>
30 
31 #include <TSystem.h>
32 
33 // Root headers
34 #include <TFile.h>
35 #include <TH2.h>
36 #include <TSystem.h>
37 
38 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
39 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
40 #include <Geant4/G4Step.hh>
41 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
42 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
43 #include <Geant4/G4String.hh> // for G4String
44 #include <Geant4/G4SystemOfUnits.hh>
45 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
46 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
47 #include <Geant4/G4Track.hh> // for G4Track
48 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
49 #include <Geant4/G4TransportationManager.hh>
50 #include <Geant4/G4Types.hh> // for G4double
51 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
52 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
53 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
54 
55 #include <cmath> // for isfinite
56 #include <cstdlib>
57 #include <filesystem>
58 #include <iostream>
59 #include <string> // for operator<<, operator+
60 #include <utility> // for pair
61 
62 class PHCompositeNode;
63 
64 //____________________________________________________________________________..
66  : PHG4SteppingAction(detector->GetName())
67  , m_Detector(detector)
68  , m_Params(parameters)
69  , m_IsActive(m_Params->get_int_param("active"))
70  , m_IsBlackHole(m_Params->get_int_param("blackhole"))
71  , m_LightScintModelFlag(m_Params->get_int_param("light_scint_model"))
72  , m_doG4Hit(m_Params->get_int_param("saveg4hit"))
73  , m_tmin(m_Params->get_double_param("tmin"))
74  , m_tmax(m_Params->get_double_param("tmax"))
75  , m_dt(m_Params->get_double_param("dt"))
76 {
77  SetLightCorrection(m_Params->get_double_param("light_balance_inner_radius") * cm,
78  m_Params->get_double_param("light_balance_inner_corr"),
79  m_Params->get_double_param("light_balance_outer_radius") * cm,
80  m_Params->get_double_param("light_balance_outer_corr"));
81 }
82 
84 {
85  // if the last hit was a zero energie deposit hit, it is just reset
86  // and the memory is still allocated, so we need to delete it here
87  // if the last hit was saved, hit is a nullptr pointer which are
88  // legal to delete (it results in a no operation)
89  delete m_Hit;
90  // since we have a copy in memory of this one - we need to delete it
91  delete m_MapCorrHist;
92 }
93 
94 //____________________________________________________________________________..
96 {
98  {
99  std::string ihcalmapname(m_Params->get_string_param("MapFileName"));
100  if (ihcalmapname.empty())
101  {
102  return 0;
103  }
104  std::string url = CDBInterface::instance()->getUrl("OLD_INNER_HCAL_TILEMAP", ihcalmapname);
105  if (!std::filesystem::exists(url))
106  {
107  std::cout << PHWHERE << " Could not locate " << url << std::endl;
108  std::cout << "use empty filename to ignore mapfile" << std::endl;
109  gSystem->Exit(1);
110  }
111  TFile* file = TFile::Open(url.c_str());
112  file->GetObject(m_Params->get_string_param("MapHistoName").c_str(), m_MapCorrHist);
113  if (!m_MapCorrHist)
114  {
115  std::cout << "ERROR: could not find Histogram " << m_Params->get_string_param("MapHistoName") << " in " << m_Params->get_string_param("MapFileName") << std::endl;
116  gSystem->Exit(1);
117  exit(1); // make code checkers which do not know gSystem->Exit() happy
118  }
119  m_MapCorrHist->SetDirectory(nullptr); // rootism: this needs to be set otherwise histo vanished when closing the file
120  file->Close();
121  delete file;
122  }
123  if (!m_doG4Hit)
124  {
125  try
126  {
127  CreateNodeTree(topNode);
128  }
129  catch (std::exception& e)
130  {
131  std::cout << e.what() << std::endl;
133  }
134  if (Verbosity() > 1) topNode->print();
135  }
136 
137  return 0;
138 }
139 
140 //____________________________________________________________________________..
142 {
143  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
144  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
145  // get volume of the current step
146  G4VPhysicalVolume* volume = touch->GetVolume();
147 
148  // m_Detector->IsInIHCal(volume)
149  // returns
150  // 0 is outside of IHCal
151  // 1 is inside scintillator
152  // -1 is steel absorber
153 
154  int whichactive = m_Detector->IsInInnerHcal(volume);
155 
156  if (!whichactive)
157  {
158  return false;
159  }
160  int layer_id = -1;
161  int tower_id = -1;
162  if (whichactive > 0) // scintillator
163  {
164  std::pair<int, int> layer_tower = m_Detector->GetLayerTowerId(volume);
165  layer_id = layer_tower.first;
166  tower_id = layer_tower.second;
167 
168  // std::cout<<"******** Inner HCal\t"<<volume->GetName()<<"\t"<<layer_id<<"\t"<<tower_id<<std::endl;
169  }
170  else
171  {
172  // absorber hit, do nothing
173  return false;
174  }
175 
176  if (!m_IsActive) return false;
177 
178  G4StepPoint* prePoint = aStep->GetPreStepPoint();
179  G4StepPoint* postPoint = aStep->GetPostStepPoint();
180  // time window cut
181  double pretime = prePoint->GetGlobalTime() / nanosecond;
182  double posttime = postPoint->GetGlobalTime() / nanosecond;
183  if (posttime < m_tmin || pretime > m_tmax) return false;
184  if ((posttime - pretime) > m_dt) return false;
185  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
186  const G4Track* aTrack = aStep->GetTrack();
187  // we only need visible energy here
188  double light_yield = eion;
189 
190  // correct evis using light map
192  {
193  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
194  if (m_MapCorrHist)
195  {
196  const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
197  const G4ThreeVector& worldPosition = postPoint->GetPosition();
198  G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
199  float lx = (localPosition.x() / cm);
200  float lz = fabs(localPosition.z() / cm);
201  // adjust to tilemap coordinates
202  int lcz = (int) (5.0 * lz) + 1;
203  int lcx = (int) (5.0 * (lx + 12.1)) + 1;
204 
205  if ((lcx >= 1) && (lcx <= m_MapCorrHist->GetNbinsY()) &&
206  (lcz >= 1) && (lcz <= m_MapCorrHist->GetNbinsX()))
207  {
208  light_yield *= (double) (m_MapCorrHist->GetBinContent(lcz, lcx));
209  }
210  else
211  {
212  light_yield = 0.0;
213  }
214  }
215  else
216  {
217  // old correction (linear ligh yield dependence along r), never tested
218  light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
219  }
220  }
221  // find the tower index for this step, tower_id is ieta, layer_id/4 is iphi
222  unsigned int ieta = tower_id;
223  unsigned int iphi = (unsigned int) layer_id / 4;
224  unsigned int tower_key = TowerInfoDefs::encode_hcal(ieta, iphi);
226  // set keep for the track
227  if (light_yield > 0)
228  {
229  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
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  return true;
239 }
240 
241 //____________________________________________________________________________..
242 bool PHG4InnerHcalSteppingAction::UserSteppingAction(const G4Step* aStep, bool)
243 {
244  if ((!m_doG4Hit) && (!m_IsBlackHole))
245  {
246  return NoHitSteppingAction(aStep);
247  }
248  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
249  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
250  // get volume of the current step
251  G4VPhysicalVolume* volume = touch->GetVolume();
252 
253  // m_Detector->IsInInnerHcal(volume)
254  // returns
255  // 0 is outside of InnerHcal
256  // 1 is inside scintillator
257  // -1 is steel absorber
258 
259  int whichactive = m_Detector->IsInInnerHcal(volume);
260 
261  if (!whichactive)
262  {
263  return false;
264  }
265  int layer_id = -1;
266  int tower_id = -1;
267  if (whichactive > 0) // scintillator
268  {
269  std::pair<int, int> layer_tower = m_Detector->GetLayerTowerId(volume);
270  layer_id = layer_tower.first;
271  tower_id = layer_tower.second;
272  // std::cout << "name " << volume->GetName() << ", mid: " << layer_id
273  // << ", twr: " << tower_id << std::endl;
274  }
275  else
276  {
277  layer_id = touch->GetCopyNumber(); // steel plate id
278  }
279  // collect energy and track length step by step
280  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
281  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
282  G4double light_yield = 0;
283  const G4Track* aTrack = aStep->GetTrack();
284 
285  // if this block stops everything, just put all kinetic energy into edep
286  if (m_IsBlackHole)
287  {
288  edep = aTrack->GetKineticEnergy() / GeV;
289  G4Track* killtrack = const_cast<G4Track*>(aTrack);
290  killtrack->SetTrackStatus(fStopAndKill);
291  }
292 
293  // make sure we are in a volume
294  if (m_IsActive)
295  {
296  bool geantino = false;
297 
298  // the check for the pdg code speeds things up, I do not want to make
299  // an expensive string compare for every track when we know
300  // geantino or chargedgeantino has pid=0
301  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
302  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
303  {
304  geantino = true;
305  }
306  G4StepPoint* prePoint = aStep->GetPreStepPoint();
307  G4StepPoint* postPoint = aStep->GetPostStepPoint();
308  // std::cout << "track id " << aTrack->GetTrackID() << std::endl;
309  // std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
310  // std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
311  switch (prePoint->GetStepStatus())
312  {
313  case fPostStepDoItProc:
314  if (m_SavePostStepStatus != fGeomBoundary)
315  {
316  break;
317  }
318  else
319  {
320  std::cout << GetName() << ": New Hit for " << std::endl;
321  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
322  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
323  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
324  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
325  std::cout << "last track: " << m_SaveTrackId
326  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
327  std::cout << "phys pre vol: " << volume->GetName()
328  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
329  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
330  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
331  }
332  [[fallthrough]];
333  case fGeomBoundary:
334  case fUndefined:
335  // if previous hit was saved, hit pointer was set to nullptr
336  // and we have to make a new one
337  if (!m_Hit)
338  {
339  m_Hit = new PHG4Hitv1();
340  }
341  // here we set the entrance values in cm
342  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
343  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
344  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
345  // time in ns
346  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
347  // set and save the track ID
348  m_Hit->set_trkid(aTrack->GetTrackID());
349  m_SaveTrackId = aTrack->GetTrackID();
350  // set the initial energy deposit
351  m_Hit->set_edep(0);
352  if (whichactive > 0) // return of IsInInnerHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
353  {
354  m_Hit->set_scint_id(tower_id); // the slat id
355  m_Hit->set_eion(0); // only implemented for v5 otherwise empty
356  m_Hit->set_raw_light_yield(0); // for scintillator only, initialize light yields
357  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
358  // Now save the container we want to add this hit to
360  }
361  else
362  {
364  }
365  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
366  {
367  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
368  {
369  m_Hit->set_trkid(pp->GetUserTrackId());
370  m_Hit->set_shower_id(pp->GetShower()->get_id());
371  m_SaveShower = pp->GetShower();
372  }
373  }
374  break;
375  default:
376  break;
377  }
378  // some sanity checks for inconsistencies
379  // check if this hit was created, if not print out last post step status
380  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
381  {
382  std::cout << GetName() << ": hit was not created" << std::endl;
383  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
384  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
385  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
386  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
387  std::cout << "last track: " << m_SaveTrackId
388  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
389  std::cout << "phys pre vol: " << volume->GetName()
390  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
391  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
392  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
393  gSystem->Exit(1);
394  }
395  // check if track id matches the initial one when the hit was created
396  if (aTrack->GetTrackID() != m_SaveTrackId)
397  {
398  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
399  std::cout << "saved track: " << m_SaveTrackId
400  << ", current trackid: " << aTrack->GetTrackID()
401  << ", prestep status: " << prePoint->GetStepStatus()
402  << ", previous post step status: " << m_SavePostStepStatus
403  << std::endl;
404 
405  gSystem->Exit(1);
406  }
407  m_SavePreStepStatus = prePoint->GetStepStatus();
408  m_SavePostStepStatus = postPoint->GetStepStatus();
409  m_SaveVolPre = volume;
410  m_SaveVolPost = touchpost->GetVolume();
411  // here we just update the exit values, it will be overwritten
412  // for every step until we leave the volume or the particle
413  // ceases to exist
414  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
415  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
416  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
417 
418  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
419 
420  // sum up the energy to get total deposited
421  m_Hit->set_edep(m_Hit->get_edep() + edep);
422  if (whichactive > 0) // return of IsInInnerHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
423  {
424  m_Hit->set_eion(m_Hit->get_eion() + eion);
425  light_yield = eion;
427  {
428  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
429  m_Hit->set_raw_light_yield(m_Hit->get_raw_light_yield() + light_yield); // save raw Birks light yield
430  if (m_MapCorrHist)
431  {
432  const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
433  const G4ThreeVector& worldPosition = postPoint->GetPosition();
434  G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
435  float lx = (localPosition.x() / cm);
436  float lz = fabs(localPosition.z() / cm);
437  // adjust to tilemap coordinates
438  int lcz = (int) (5.0 * lz) + 1;
439  int lcx = (int) (5.0 * (lx + 12.1)) + 1;
440 
441  if ((lcx >= 1) && (lcx <= m_MapCorrHist->GetNbinsY()) &&
442  (lcz >= 1) && (lcz <= m_MapCorrHist->GetNbinsX()))
443  {
444  light_yield *= (double) (m_MapCorrHist->GetBinContent(lcz, lcx));
445  }
446  else
447  {
448  light_yield = 0.0;
449  }
450  }
451  else
452  {
453  // old correction (linear ligh yield dependence along r), never tested
454  light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
455  }
456  }
457  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
458  }
459  if (geantino)
460  {
461  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
462  m_Hit->set_eion(-1);
463  }
464  if (edep > 0)
465  {
466  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
467  {
468  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
469  {
470  pp->SetKeep(1); // we want to keep the track
471  }
472  }
473  }
474 
475  // if any of these conditions is true this is the last step in
476  // this volume and we need to save the hit
477  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
478  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
479  // (happens when your detector goes outside world volume)
480  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
481  // aTrack->GetTrackStatus() == fStopAndKill is also set)
482  // aTrack->GetTrackStatus() == fStopAndKill: track ends
483  if (postPoint->GetStepStatus() == fGeomBoundary ||
484  postPoint->GetStepStatus() == fWorldBoundary ||
485  postPoint->GetStepStatus() == fAtRestDoItProc ||
486  aTrack->GetTrackStatus() == fStopAndKill)
487  {
488  // save only hits with energy deposit (or -1 for geantino)
489  if (m_Hit->get_edep() != 0)
490  {
491  m_SaveHitContainer->AddHit(layer_id, m_Hit);
492 
493  if (m_SaveShower)
494  {
496  }
497  // ownership has been transferred to container, set to null
498  // so we will create a new hit for the next track
499  m_Hit = nullptr;
500  }
501  else
502  {
503  // if this hit has no energy deposit, just reset it for reuse
504  // this means we have to delete it in the dtor. If this was
505  // the last hit we processed the memory is still allocated
506  m_Hit->Reset();
507  }
508  }
509  // return true to indicate the hit was used
510  return true;
511  }
512  else
513  {
514  return false;
515  }
516 }
517 
518 //____________________________________________________________________________..
520 {
521  std::string hitnodename;
522  std::string absorbernodename;
523  if (m_Detector->SuperDetector() != "NONE")
524  {
525  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
526  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
527  }
528  else
529  {
530  hitnodename = "G4HIT_" + m_Detector->GetName();
531  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
532  }
533 
534  // now look for the map and grab a pointer to it.
535  m_Hits = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
536  m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
537 
538  // if we do not find the node it's messed up.
539  if (!m_Hits)
540  {
541  std::cout << "PHG4InnerHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
542  }
543  if (!m_AbsorberHits)
544  {
545  if (Verbosity() > 1)
546  {
547  std::cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << std::endl;
548  }
549  }
550 }
551 
553 {
554  PHNodeIterator nodeItr(topNode);
555  PHCompositeNode* dst_node = dynamic_cast<PHCompositeNode*>(
556  nodeItr.findFirst("PHCompositeNode", "DST"));
557  if (!dst_node)
558  {
559  std::cout << "PHComposite node created: DST" << std::endl;
560  dst_node = new PHCompositeNode("DST");
561  topNode->addNode(dst_node);
562  }
563  PHNodeIterator dstiter(dst_node);
564  PHCompositeNode* DetNode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", m_Detector->SuperDetector()));
565  if (!DetNode)
566  {
567  DetNode = new PHCompositeNode(m_Detector->SuperDetector());
568  dst_node->addNode(DetNode);
569  }
570  m_CaloInfoContainer = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::HCAL);
571  PHIODataNode<PHObject>* towerNode = new PHIODataNode<PHObject>(m_CaloInfoContainer, "TOWERINFO_SIM_" + m_Detector->SuperDetector(), "PHObject");
572  DetNode->addNode(towerNode);
573 }