ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4OuterHcalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4OuterHcalSteppingAction.cc
1 // local headers in quotes (that is important when using include subdirs!)
3 
4 #include "PHG4HcalDefs.h"
6 #include "PHG4StepStatusDecode.h"
7 
8 // our own headers in alphabetical order
9 
10 #include <phparameter/PHParameters.h>
11 
12 #include <fun4all/Fun4AllServer.h>
13 
14 #include <g4main/PHG4Hit.h>
16 #include <g4main/PHG4Hitv1.h>
17 #include <g4main/PHG4Shower.h>
18 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
20 
21 #include <phool/getClass.h>
22 
23 // Root headers
24 #include <TAxis.h> // for TAxis
25 #include <TH2.h>
26 #include <TH2F.h>
27 #include <TNamed.h> // for TNamed
28 #include <TSystem.h>
29 #include <TFile.h>
30 
31 // Geant4 headers
32 
33 #include <Geant4/G4Field.hh>
34 #include <Geant4/G4FieldManager.hh>
35 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
36 #include <Geant4/G4PropagatorInField.hh>
37 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
38 #include <Geant4/G4Step.hh>
39 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
40 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
41 #include <Geant4/G4String.hh> // for G4String
42 #include <Geant4/G4SystemOfUnits.hh>
43 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
44 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
45 #include <Geant4/G4Track.hh> // for G4Track
46 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
47 #include <Geant4/G4TransportationManager.hh>
48 #include <Geant4/G4Types.hh> // for G4double
49 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
50 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
51 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
52 #include <Geant4/G4Transform3D.hh>
53 
54 // finally system headers
55 #include <cassert>
56 #include <cmath> // for isfinite, sqrt
57 #include <iostream>
58 #include <string> // for operator<<, string
59 #include <utility> // for pair
60 
61 class PHCompositeNode;
62 
63 using namespace std;
64 
65 TH2F *MapCorr = NULL;
66 
67 //____________________________________________________________________________..
69  : PHG4SteppingAction(detector->GetName())
70  , m_Detector(detector)
71  , m_Hits(nullptr)
72  , m_AbsorberHits(nullptr)
73  , m_Hit(nullptr)
74  , m_Params(parameters)
75  , m_SaveHitContainer(nullptr)
76  , m_SaveShower(nullptr)
77  , m_SaveVolPre(nullptr)
78  , m_SaveVolPost(nullptr)
79  , m_SaveTrackId(-1)
80  , m_SavePreStepStatus(-1)
81  , m_SavePostStepStatus(-1)
82  , m_EnableFieldCheckerFlag(m_Params->get_int_param("field_check"))
83  , m_IsActiveFlag(m_Params->get_int_param("active"))
84  , m_IsBlackHoleFlag(m_Params->get_int_param("blackhole"))
85  , m_NScintiPlates(m_Params->get_int_param(PHG4HcalDefs::scipertwr) * m_Params->get_int_param("n_towers"))
86  , m_LightScintModelFlag(m_Params->get_int_param("light_scint_model"))
87 {
89 }
90 
92 {
93  // if the last hit was a zero energie deposit hit, it is just reset
94  // and the memory is still allocated, so we need to delete it here
95  // if the last hit was saved, hit is a nullptr pointer which are
96  // legal to delete (it results in a no operation)
97  delete m_Hit;
98 }
99 
101 {
103 // method in base class for light correction
104  SetLightCorrection(m_Params->get_double_param("light_balance_inner_radius") * cm,
105  m_Params->get_double_param("light_balance_inner_corr"),
106  m_Params->get_double_param("light_balance_outer_radius") * cm,
107  m_Params->get_double_param("light_balance_outer_corr")
108  );
109 
110 
111  std::ostringstream mappingfilename;
112  const char* calibroot = getenv("CALIBRATIONROOT");
113  if (calibroot)
114  {
115  mappingfilename << calibroot;
116  }
117  else
118  {
119  std::cout << "no CALIBRATIONROOT environment variable" << std::endl;
120  gSystem->Exit(1);
121  }
122 
123  mappingfilename << "/HCALOUT/tilemap/oHCALMaps092021.root";
124  TFile *f = new TFile(mappingfilename.str().data());
125  MapCorr = (TH2F *)f->Get("hCombinedMap");
126  if(!MapCorr){
127  std::cout << "ERROR: MapCorr is NULL" << std::endl;
128  gSystem->Exit(1);
129  }
130 
131  return 0;
132 }
133 
134 //____________________________________________________________________________..
136 {
138  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
139  // get volume of the current step
140  G4VPhysicalVolume* volume = touch->GetVolume();
141 
142  // m_Detector->IsInOuterHcal(volume)
143  // returns
144  // 0 is outside of OuterHcal
145  // 1 is inside scintillator
146  // -1 is steel absorber (if absorber set to active)
147 
148  int whichactive = m_Detector->IsInOuterHcal(volume);
149 
150  if (!whichactive)
151  {
152  return false;
153  }
154 
156  {
157  FieldChecker(aStep);
158  }
159 
160  int layer_id = -1;
161  int tower_id = -1;
162  if (whichactive > 0) // scintillator
163  {
164  pair<int, int> layer_tower = m_Detector->GetLayerTowerId(volume);
165  layer_id = layer_tower.first;
166  tower_id = layer_tower.second;
167  }
168  else
169  {
170  layer_id = touch->GetCopyNumber(); // steel plate id
171  }
172  // collect energy and track length step by step
173  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
174  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
175  G4double light_yield = 0;
176 
177  const G4Track* aTrack = aStep->GetTrack();
178 
179  // if this block stops everything, just put all kinetic energy into edep
180  if (m_IsBlackHoleFlag)
181  {
182  edep = aTrack->GetKineticEnergy() / GeV;
183  G4Track* killtrack = const_cast<G4Track*>(aTrack);
184  killtrack->SetTrackStatus(fStopAndKill);
185  }
186 
187  // make sure we are in a volume
188  if (m_IsActiveFlag)
189  {
190  bool geantino = false;
191 
192  // the check for the pdg code speeds things up, I do not want to make
193  // an expensive string compare for every track when we know
194  // geantino or chargedgeantino has pid=0
195  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
196  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
197  {
198  geantino = true;
199  }
200  G4StepPoint* prePoint = aStep->GetPreStepPoint();
201  G4StepPoint* postPoint = aStep->GetPostStepPoint();
202  // cout << "track id " << aTrack->GetTrackID() << endl;
203  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
204  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
205 
206  switch (prePoint->GetStepStatus())
207  {
208  case fPostStepDoItProc:
210  {
211  break;
212  }
213  else
214  {
215  cout << GetName() << ": New Hit for " << endl;
216  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
217  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
218  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
219  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
220  cout << "last track: " << m_SaveTrackId
221  << ", current trackid: " << aTrack->GetTrackID() << endl;
222  cout << "phys pre vol: " << volume->GetName()
223  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
224  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
225  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
226  }
227  [[fallthrough]];
228  case fGeomBoundary:
229  case fUndefined:
230  if (!m_Hit)
231  {
232  m_Hit = new PHG4Hitv1();
233  }
234  //here we set the entrance values in cm
235  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
236  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
237  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
238 
239  // DEBUG
240  // add the local coordinates
241  // if(whichactive>0){
242 
243  // G4TouchableHandle theTouchable = prePoint->GetTouchableHandle();
244  // G4ThreeVector worldPosition = prePoint->GetPosition();
245  // G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
246 
247  // m_Hit->set_property(PHG4Hit::prop_local_x_0, (float)(localPosition.x()/cm));
248  // m_Hit->set_property(PHG4Hit::prop_local_y_0, (float)(localPosition.y()/cm));
249  // m_Hit->set_property(PHG4Hit::prop_local_z_0, (float)(localPosition.z()/cm));
250 
251  // m_Hit->set_property(PHG4Hit::prop_layer, (unsigned int) layer_id);
252 
253  // }
254 
255  // time in ns
256  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
257  //set the track ID
258  m_Hit->set_trkid(aTrack->GetTrackID());
259  m_SaveTrackId = aTrack->GetTrackID();
260  //set the initial energy deposit
261  m_Hit->set_edep(0);
262  if (whichactive > 0) // return of IsInOuterHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
263  {
264  m_Hit->set_scint_id(tower_id); // the slat id
265  m_Hit->set_eion(0);
266  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
267  // Now save the container we want to add this hit to
269  }
270  else
271  {
273  }
275  {
276  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
277  {
278  m_Hit->set_trkid(pp->GetUserTrackId());
279  m_Hit->set_shower_id(pp->GetShower()->get_id());
280  m_SaveShower = pp->GetShower();
281  }
282  }
283  break;
284  default:
285  break;
286  }
287  // some sanity checks for inconsistencies
288  // check if this hit was created, if not print out last post step status
289  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
290  {
291  cout << GetName() << ": hit was not created" << endl;
292  cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
293  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
294  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
295  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
296  cout << "last track: " << m_SaveTrackId
297  << ", current trackid: " << aTrack->GetTrackID() << endl;
298  cout << "phys pre vol: " << volume->GetName()
299  << " post vol : " << touchpost->GetVolume()->GetName() << endl;
300  cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
301  << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
302  gSystem->Exit(1);
303  }
304  m_SavePostStepStatus = postPoint->GetStepStatus();
305  // check if track id matches the initial one when the hit was created
306  if (aTrack->GetTrackID() != m_SaveTrackId)
307  {
308  cout << GetName() << ": hits do not belong to the same track" << endl;
309  cout << "saved track: " << m_SaveTrackId
310  << ", current trackid: " << aTrack->GetTrackID()
311  << endl;
312  gSystem->Exit(1);
313  }
314  m_SavePreStepStatus = prePoint->GetStepStatus();
315  m_SavePostStepStatus = postPoint->GetStepStatus();
317  m_SaveVolPost = touchpost->GetVolume();
318  // here we just update the exit values, it will be overwritten
319  // for every step until we leave the volume or the particle
320  // ceases to exist
321  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
322  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
323  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
324 
325  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
326 
327  if (whichactive > 0)
328  {
329 
330  // Local Coordinates:
331 
332  //G4TouchableHandle theTouchable = postPoint->GetTouchableHandle();
333  // Use prePoint; sometimes the end point can be on the boundary/out of the scintillator
334  G4TouchableHandle theTouchable = prePoint->GetTouchableHandle();
335  G4ThreeVector worldPosition = postPoint->GetPosition();
336  G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
337 
338  // DEBUG
339  // m_Hit->set_property(PHG4Hit::prop_local_x_1, (float)(localPosition.x()/cm));
340  // m_Hit->set_property(PHG4Hit::prop_local_y_1, (float)(localPosition.y()/cm));
341  // m_Hit->set_property(PHG4Hit::prop_local_z_1, (float)(localPosition.z()/cm));
342 
343  // m_Hit->set_property(PHG4Hit::prop_layer, (unsigned int) layer_id);
344 
346  {
347  light_yield = GetVisibleEnergyDeposition(aStep);
348 
349  if(MapCorr){
350 
351  float lx = (localPosition.x()/cm);
352  float lz = fabs(localPosition.z()/cm); // reverse the sense for towerid<12
353 
354  // convert to the map bin coordinates:
355  // map is in 0.5 cm bins
356  int lcz = (int)(2.0*lz) + 1;
357  int lcx = (int)(2.0*(lx+42.75)) + 1;
358 
359  if( (lcx>=1) && (lcx<=MapCorr->GetNbinsY()) &&
360  (lcz>=1) && (lcz<=MapCorr->GetNbinsX()) ){
361 
362  light_yield *= (double) (MapCorr->GetBinContent(lcz, lcx));
363 
364  }
365  else{
366  light_yield = 0.0;
367  }
368 
369  }
370 
371  }
372  else
373  {
374  light_yield = eion;
375  }
376 
377  if (ValidCorrection())
378  {
379  double cor = GetLightCorrection(postPoint->GetPosition().x() , (postPoint->GetPosition().y() ));
380  cout << "applying cor: " << cor << endl;
381  light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x() , (postPoint->GetPosition().y() ));
382  }
383  }
384 
385  //sum up the energy to get total deposited
387  if (whichactive > 0)
388  {
389  m_Hit->set_eion(m_Hit->get_eion() + eion);
390  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
391  }
392  if (geantino)
393  {
394  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
395  m_Hit->set_eion(-1);
396  }
397  if (edep > 0)
398  {
400  {
401  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
402  {
403  pp->SetKeep(1); // we want to keep the track
404  }
405  }
406  }
407 
408  // if any of these conditions is true this is the last step in
409  // this volume and we need to save the hit
410  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
411  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
412  // (happens when your detector goes outside world volume)
413  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
414  // aTrack->GetTrackStatus() == fStopAndKill is also set)
415  // aTrack->GetTrackStatus() == fStopAndKill: track ends
416  if (postPoint->GetStepStatus() == fGeomBoundary ||
417  postPoint->GetStepStatus() == fWorldBoundary ||
418  postPoint->GetStepStatus() == fAtRestDoItProc ||
419  aTrack->GetTrackStatus() == fStopAndKill)
420  {
421  // save only hits with energy deposit (or -1 for geantino)
422  if (m_Hit->get_edep())
423  {
424  m_SaveHitContainer->AddHit(layer_id, m_Hit);
425  if (m_SaveShower)
426  {
428  }
429  // ownership has been transferred to container, set to null
430  // so we will create a new hit for the next track
431  m_Hit = nullptr;
432  }
433  else
434  {
435  // if this hit has no energy deposit, just reset it for reuse
436  // this means we have to delete it in the dtor. If this was
437  // the last hit we processed the memory is still allocated
438  m_Hit->Reset();
439  }
440  }
441  // return true to indicate the hit was used
442  return true;
443  }
444  else
445  {
446  return false;
447  }
448 }
449 
450 //____________________________________________________________________________..
452 {
453  string hitnodename;
454  string absorbernodename;
455  if (m_Detector->SuperDetector() != "NONE")
456  {
457  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
458  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
459  }
460  else
461  {
462  hitnodename = "G4HIT_" + m_Detector->GetName();
463  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
464  }
465 
466  //now look for the map and grab a pointer to it.
467  m_Hits = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
468  m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
469 
470  // if we do not find the node it's messed up.
471  if (!m_Hits)
472  {
473  std::cout << "PHG4OuterHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
474  }
475  if (!m_AbsorberHits)
476  {
477  if (Verbosity() > 1)
478  {
479  cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
480  }
481  }
482 }
483 
485 {
487  assert(se);
488 
489  static const string h_field_name = "hOuterHcalField";
490 
491  if (not se->isHistoRegistered(h_field_name))
492  {
493  TH2F* h = new TH2F(h_field_name.c_str(), "Magnetic field (Tesla) in HCal;X (cm);Y (cm)", 2400,
494  -300, 300, 2400, -300, 300);
495 
496  se->registerHisto(h, 1);
497 
498  cout << "PHG4OuterHcalSteppingAction::FieldChecker - make a histograme to check outer Hcal field map."
499  << " Saved to Fun4AllServer Histo with name " << h_field_name << endl;
500  }
501 
502  TH2F* h = dynamic_cast<TH2F*>(se->getHisto(h_field_name));
503  assert(h);
504 
505  assert(aStep);
507  assert(touch);
508  // get volume of the current step
509  G4VPhysicalVolume* volume = touch->GetVolume();
510  assert(volume);
511 
512  G4ThreeVector globPosition = aStep->GetPreStepPoint()->GetPosition();
513 
514  G4double globPosVec[4] = {0};
515  G4double FieldValueVec[6] = {0};
516 
517  globPosVec[0] = globPosition.x();
518  globPosVec[1] = globPosition.y();
519  globPosVec[2] = globPosition.z();
520  globPosVec[3] = aStep->GetPreStepPoint()->GetGlobalTime();
521 
522  const Int_t binx = h->GetXaxis()->FindBin(globPosVec[0] / cm);
523  const Int_t biny = h->GetYaxis()->FindBin(globPosVec[1] / cm);
524 
525  if (h->GetBinContent(binx, binx) == 0)
526  { // only fille unfilled bins
527 
528  G4TransportationManager* transportMgr =
530  assert(transportMgr);
531 
532  G4PropagatorInField* fFieldPropagator =
533  transportMgr->GetPropagatorInField();
534  assert(fFieldPropagator);
535 
536  G4FieldManager* fieldMgr = fFieldPropagator->FindAndSetFieldManager(volume);
537  assert(fieldMgr);
538 
539  const G4Field* pField = fieldMgr->GetDetectorField();
540  assert(pField);
541 
542  pField->GetFieldValue(globPosVec, FieldValueVec);
543 
544  G4ThreeVector FieldValue = G4ThreeVector(FieldValueVec[0],
545  FieldValueVec[1], FieldValueVec[2]);
546 
547  const double B = FieldValue.mag() / tesla;
548 
549  h->SetBinContent(binx, biny, B);
550 
551  cout << "PHG4OuterHcalSteppingAction::FieldChecker - "
552  << "bin " << binx
553  << ", " << biny << " := " << B << " Tesla @ x,y = " << globPosVec[0] / cm
554  << "," << globPosVec[1] / cm << " cm" << endl;
555  }
556 }