Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4OHCalSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4OHCalSteppingAction.cc
2 
3 #include "PHG4OHCalDetector.h"
4 
5 // our own headers in alphabetical order
6 
9 
10 #include <phparameter/PHParameters.h>
11 
12 #include <g4main/PHG4Hit.h>
14 #include <g4main/PHG4Hitv1.h>
15 #include <g4main/PHG4Shower.h>
16 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
18 
20 #include <fun4all/Fun4AllServer.h>
21 
22 #include <phool/PHCompositeNode.h>
23 #include <phool/PHIODataNode.h> // for PHIODataNode
24 #include <phool/PHNode.h> // for PHNode
25 #include <phool/PHNodeIterator.h> // for PHNodeIterator
26 #include <phool/PHObject.h> // for PHObject
27 #include <phool/getClass.h>
28 
29 #include <calobase/TowerInfo.h>
30 #include <calobase/TowerInfoContainer.h>
31 #include <calobase/TowerInfoContainerv1.h>
32 #include <calobase/TowerInfoDefs.h>
33 
34 // Root headers
35 #include <TAxis.h> // for TAxis
36 #include <TFile.h>
37 #include <TH2.h>
38 #include <TNamed.h> // for TNamed
39 #include <TSystem.h>
40 
41 // Geant4 headers
42 
43 #include <Geant4/G4AffineTransform.hh> // for G4AffineTransform
44 #include <Geant4/G4Field.hh>
45 #include <Geant4/G4FieldManager.hh>
46 #include <Geant4/G4LogicalVolume.hh> // for G4LogicalVolume
47 #include <Geant4/G4NavigationHistory.hh> // for G4NavigationHistory
48 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
49 #include <Geant4/G4PropagatorInField.hh>
50 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
51 #include <Geant4/G4Step.hh>
52 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
53 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
54 #include <Geant4/G4String.hh> // for G4String
55 #include <Geant4/G4SystemOfUnits.hh>
56 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
57 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
58 #include <Geant4/G4Track.hh> // for G4Track
59 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
60 #include <Geant4/G4TransportationManager.hh>
61 #include <Geant4/G4Types.hh> // for G4double
62 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
63 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
64 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
65 
66 // finally system headers
67 #include <cassert>
68 #include <cmath> // for isfinite, sqrt
69 #include <filesystem>
70 #include <iostream>
71 #include <string> // for operator<<, string
72 #include <tuple> // for get, tuple
73 
74 class PHCompositeNode;
75 
76 //____________________________________________________________________________..
78  : PHG4SteppingAction(detector->GetName())
79  , m_Detector(detector)
80  , m_Params(parameters)
81  , m_EnableFieldCheckerFlag(m_Params->get_int_param("field_check"))
82  , m_IsActiveFlag(m_Params->get_int_param("active"))
83  , m_IsBlackHoleFlag(m_Params->get_int_param("blackhole"))
84  , m_NScintiPlates(m_Params->get_int_param(PHG4HcalDefs::scipertwr) * m_Params->get_int_param("n_towers"))
85  , m_LightScintModelFlag(m_Params->get_int_param("light_scint_model"))
86  , m_doG4Hit(m_Params->get_int_param("saveg4hit"))
87  , m_tmin(m_Params->get_double_param("tmin"))
88  , m_tmax(m_Params->get_double_param("tmax"))
89  , m_dt(m_Params->get_double_param("dt"))
90 {
91  SetLightCorrection(m_Params->get_double_param("light_balance_inner_radius") * cm,
92  m_Params->get_double_param("light_balance_inner_corr"),
93  m_Params->get_double_param("light_balance_outer_radius") * cm,
94  m_Params->get_double_param("light_balance_outer_corr"));
95 }
96 
98 {
99  // if the last hit was a zero energie deposit hit, it is just reset
100  // and the memory is still allocated, so we need to delete it here
101  // if the last hit was saved, hit is a nullptr pointer which are
102  // legal to delete (it results in a no operation)
103  delete m_Hit;
104 
105  // since we have a copy in memory of this one - we need to delete it
106  for (auto it : m_MapCorrHist)
107  {
108  delete it;
109  }
110  for (int j = 0; j < 4; j++)
111  {
112  delete m_MapCorrHistChim[j];
113  }
114 }
115 
117 {
119 
121  {
122  std::string mappingfilename(m_Params->get_string_param("MapFileName"));
123  if (mappingfilename.empty())
124  {
125  return 0;
126  }
127  if (!std::filesystem::exists(m_Params->get_string_param("MapFileName")))
128  {
129  std::cout << PHWHERE << " Could not locate " << m_Params->get_string_param("MapFileName") << std::endl;
130  std::cout << "use empty filename to ignore mapfile" << std::endl;
131  gSystem->Exit(1);
132  }
133 
134  TFile* file = TFile::Open(mappingfilename.c_str());
135  std::string Tilehist = m_Params->get_string_param("MapHistoName");
136 
137  for (int i = 0; i < 24; i++)
138  {
139  std::string str2 = std::to_string(i);
140  Tilehist += str2;
141  file->GetObject(Tilehist.c_str(), m_MapCorrHist[i]);
142  if (i < 4)
143  {
144  Tilehist += "_chimney";
145  }
146  file->GetObject(Tilehist.c_str(), m_MapCorrHistChim[i]);
147  Tilehist = m_Params->get_string_param("MapHistoName");
148 
149  if ((!m_MapCorrHist[i]) || (!m_MapCorrHistChim[i]))
150  {
151  std::cout << "ERROR: could not find Histogram " << Tilehist << i << " in " << m_Params->get_string_param("MapFileName") << std::endl;
152  gSystem->Exit(1);
153  }
154 
155  m_MapCorrHist[i]->SetDirectory(nullptr); // rootism: this needs to be set otherwise histo vanished when closing the file
156  m_MapCorrHistChim[i]->SetDirectory(nullptr);
157  }
158 
159  file->Close();
160  delete file;
161  }
162  if (!m_doG4Hit)
163  {
164  try
165  {
166  CreateNodeTree(topNode);
167  }
168  catch (std::exception& e)
169  {
170  std::cout << e.what() << std::endl;
172  }
173  if (Verbosity() > 1) topNode->print();
174  }
175  return 0;
176 }
177 
178 //____________________________________________________________________________..
180 {
181  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
182  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
183  // get volume of the current step
184  G4VPhysicalVolume* volume = touch->GetVolume();
185 
186  // m_Detector->IsInIHCal(volume)
187  // returns
188  // 0 is outside of IHCal
189  // 1 is inside scintillator
190  // -1 is steel absorber
191 
192  int whichactive = m_Detector->IsInOHCal(volume);
193 
194  if (!whichactive)
195  {
196  return false;
197  }
199  {
200  FieldChecker(aStep);
201  }
202  int layer_id = -1;
203  int tower_id = -1;
204  int sector_id = -1;
205  if (whichactive > 0) // scintillator
206  {
207  std::tuple<int, int, int> layer_tower = m_Detector->GetRowColumnId(volume);
208  sector_id = std::get<0>(layer_tower);
209  layer_id = std::get<1>(layer_tower);
210  tower_id = std::get<2>(layer_tower);
211 
212  // std::cout<<"******** Outer HCal\t"<<volume->GetName()<<"\t"<<layer_id<<"\t"<<tower_id<<std::endl;
213  }
214  else
215  {
216  // absorber hit, do nothing
217  return false;
218  }
219 
220  if (!m_IsActiveFlag) return false;
221 
222  G4StepPoint* prePoint = aStep->GetPreStepPoint();
223  G4StepPoint* postPoint = aStep->GetPostStepPoint();
224  // time window cut
225  double pretime = prePoint->GetGlobalTime() / nanosecond;
226  double posttime = postPoint->GetGlobalTime() / nanosecond;
227  if (posttime < m_tmin || pretime > m_tmax) return false;
228  if ((posttime - pretime) > m_dt) return false;
229  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
230  const G4Track* aTrack = aStep->GetTrack();
231  // we only need visible energy here
232  double light_yield = eion;
233 
234  // correct evis using light map
236  {
237  light_yield = GetVisibleEnergyDeposition(aStep);
238  if ((m_MapCorrHistChim[tower_id]) || (m_MapCorrHist[tower_id]))
239  {
240  const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
241  const G4ThreeVector& worldPosition = postPoint->GetPosition();
242  G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
243  float lx = (localPosition.x() / cm);
244  float ly = (localPosition.y() / cm);
245 
246  // convert to the map bin coordinates:
247  int lcx = (int) (2.0 * lx) + 1;
248  int lcy = (int) (2.0 * (ly + 0.5)) + 1;
249 
250  if ((sector_id == 29) || (sector_id == 30) || (sector_id == 31))
251  {
252  if ((lcy >= 1) && (lcy <= m_MapCorrHistChim[tower_id]->GetNbinsY()) &&
253  (lcx >= 1) && (lcx <= m_MapCorrHistChim[tower_id]->GetNbinsX()))
254  {
255  light_yield *= (double) (m_MapCorrHistChim[tower_id]->GetBinContent(lcx, lcy));
256  }
257  else
258  {
259  light_yield = 0.0;
260  }
261  }
262  else
263  {
264  if ((lcy >= 1) && (lcy <= m_MapCorrHist[tower_id]->GetNbinsY()) &&
265  (lcx >= 1) && (lcx <= m_MapCorrHist[tower_id]->GetNbinsX()))
266  {
267  light_yield *= (double) (m_MapCorrHist[tower_id]->GetBinContent(lcx, lcy));
268  }
269  else
270  {
271  light_yield = 0.0;
272  }
273  }
274  }
275  else
276  {
277  // old correction (linear ligh yield dependence along r), never tested
278  light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
279  }
280  }
281  // find the tower index for this step, tower_id is ieta, layer_id/5 is iphi
282  unsigned int ieta = tower_id;
283  unsigned int iphi = (unsigned int) layer_id / 5;
284  unsigned int tower_key = TowerInfoDefs::encode_hcal(ieta, iphi);
286  // set keep for the track
287  if (light_yield > 0)
288  {
289  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
290  {
291  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
292  {
293  pp->SetKeep(1); // we want to keep the track
294  }
295  }
296  }
297 
298  return true;
299 }
300 
301 //____________________________________________________________________________..
302 bool PHG4OHCalSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
303 {
304  if ((!m_doG4Hit) && (!m_IsBlackHoleFlag))
305  {
306  return NoHitSteppingAction(aStep);
307  }
308 
309  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
310  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
311  // get volume of the current step
312  G4VPhysicalVolume* volume = touch->GetVolume();
313 
314  // m_Detector->IsInOHCal(volume)
315  // returns
316  // 0 is outside of OHCal
317  // 1 is inside scintillator
318  // -1 is steel absorber (if absorber set to active)
319 
320  int whichactive = m_Detector->IsInOHCal(volume);
321 
322  if (!whichactive)
323  {
324  return false;
325  }
326 
328  {
329  FieldChecker(aStep);
330  }
331 
332  int layer_id = -1;
333  int tower_id = -1;
334  int sector_id = -1;
335  if (whichactive > 0) // scintillator
336  {
337  std::tuple<int, int, int> layer_tower = m_Detector->GetRowColumnId(volume);
338  sector_id = std::get<0>(layer_tower); // calorimeter sector (0-31)
339  layer_id = std::get<1>(layer_tower); // calorimeter scintillator row (0-319)
340  tower_id = std::get<2>(layer_tower); // calorimeter scintillator tower in row (0-23)
341  }
342  else
343  {
344  layer_id = touch->GetCopyNumber(); // steel plate id
345  }
346 
347  // collect energy and track length step by step
348  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
349  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
350  G4double light_yield = 0;
351 
352  const G4Track* aTrack = aStep->GetTrack();
353 
354  // if this block stops everything, just put all kinetic energy into edep
355  if (m_IsBlackHoleFlag)
356  {
357  edep = aTrack->GetKineticEnergy() / GeV;
358  G4Track* killtrack = const_cast<G4Track*>(aTrack);
359  killtrack->SetTrackStatus(fStopAndKill);
360  }
361 
362  // make sure we are in a volume
363  if (m_IsActiveFlag)
364  {
365  bool geantino = false;
366 
367  // the check for the pdg code speeds things up, I do not want to make
368  // an expensive string compare for every track when we know
369  // geantino or chargedgeantino has pid=0
370  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
371  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
372  {
373  geantino = true;
374  }
375  G4StepPoint* prePoint = aStep->GetPreStepPoint();
376  G4StepPoint* postPoint = aStep->GetPostStepPoint();
377  // std::cout << "track id " << aTrack->GetTrackID() << std::endl;
378  // std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
379  // std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
380 
381  switch (prePoint->GetStepStatus())
382  {
383  case fPostStepDoItProc:
384  if (m_SavePostStepStatus != fGeomBoundary)
385  {
386  if (m_SavePostStepStatus != fAtRestDoItProc)
387  {
388  break;
389  }
390  else
391  {
392  if (aTrack->GetTrackID() == m_SaveTrackId)
393  {
394  std::cout << GetName() << ": Bad step status combination for the same track " << std::endl;
395  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
396  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
397  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
398  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
399  std::cout << "last track: " << m_SaveTrackId
400  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
401  std::cout << "phys pre vol: " << volume->GetName()
402  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
403  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
404  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
405  gSystem->Exit(1);
406  }
407  }
408  }
409  else
410  {
411  std::cout << GetName() << ": New Hit for " << std::endl;
412  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
413  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
414  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
415  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
416  std::cout << "last track: " << m_SaveTrackId
417  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
418  std::cout << "phys pre vol: " << volume->GetName()
419  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
420  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
421  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
422  }
423  [[fallthrough]];
424  case fGeomBoundary:
425  case fUndefined:
426  if (!m_Hit)
427  {
428  m_Hit = new PHG4Hitv1();
429  }
430  // here we set the entrance values in cm
431  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
432  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
433  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
434 
435  // time in ns
436  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
437  // set the track ID
438  m_Hit->set_trkid(aTrack->GetTrackID());
439  m_SaveTrackId = aTrack->GetTrackID();
440  // set the initial energy deposit
441  m_Hit->set_edep(0);
442  if (whichactive > 0) // return of IsInOHCalDetector, > 0 hit in scintillator, < 0 hit in absorber
443  {
444  m_Hit->set_sector(sector_id); // the sector id
445  m_Hit->set_scint_id(tower_id); // the slat id
446  m_Hit->set_eion(0);
447  m_Hit->set_raw_light_yield(0); // for scintillator only, initialize light yields
448  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
449  // Now save the container we want to add this hit to
451  }
452  else
453  {
455  }
456  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
457  {
458  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
459  {
460  m_Hit->set_trkid(pp->GetUserTrackId());
461  m_Hit->set_shower_id(pp->GetShower()->get_id());
462  m_SaveShower = pp->GetShower();
463  }
464  }
465  break;
466  default:
467  break;
468  }
469  // some sanity checks for inconsistencies
470  // check if this hit was created, if not print out last post step status
471  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
472  {
473  std::cout << GetName() << ": hit was not created" << std::endl;
474  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
475  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
476  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
477  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
478  std::cout << "last track: " << m_SaveTrackId
479  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
480  std::cout << "phys pre vol: " << volume->GetName()
481  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
482  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
483  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
484  gSystem->Exit(1);
485  }
486  m_SavePostStepStatus = postPoint->GetStepStatus();
487  // check if track id matches the initial one when the hit was created
488  if (aTrack->GetTrackID() != m_SaveTrackId)
489  {
490  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
491  std::cout << "saved track: " << m_SaveTrackId
492  << ", current trackid: " << aTrack->GetTrackID()
493  << std::endl;
494  gSystem->Exit(1);
495  }
496  m_SavePreStepStatus = prePoint->GetStepStatus();
497  m_SavePostStepStatus = postPoint->GetStepStatus();
498  m_SaveVolPre = volume;
499  m_SaveVolPost = touchpost->GetVolume();
500  // here we just update the exit values, it will be overwritten
501  // for every step until we leave the volume or the particle
502  // ceases to exist
503  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
504  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
505  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
506 
507  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
508 
509  // sum up the energy to get total deposited
510  m_Hit->set_edep(m_Hit->get_edep() + edep);
511 
512  if (whichactive > 0)
513  {
514  m_Hit->set_eion(m_Hit->get_eion() + eion);
515  light_yield = eion;
517  {
518  light_yield = GetVisibleEnergyDeposition(aStep);
519  m_Hit->set_raw_light_yield(m_Hit->get_raw_light_yield() + light_yield); // save raw Birks light yield
520  if ((m_MapCorrHistChim[tower_id]) || (m_MapCorrHist[tower_id]))
521  {
522  const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
523  const G4ThreeVector& worldPosition = postPoint->GetPosition();
524  G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
525  float lx = (localPosition.x() / cm);
526  float ly = (localPosition.y() / cm);
527 
528  // convert to the map bin coordinates:
529  int lcx = (int) (2.0 * lx) + 1;
530  int lcy = (int) (2.0 * (ly + 0.5)) + 1;
531 
532  if ((sector_id == 29) || (sector_id == 30) || (sector_id == 31))
533  {
534  if ((lcy >= 1) && (lcy <= m_MapCorrHistChim[tower_id]->GetNbinsY()) &&
535  (lcx >= 1) && (lcx <= m_MapCorrHistChim[tower_id]->GetNbinsX()))
536  {
537  light_yield *= (double) (m_MapCorrHistChim[tower_id]->GetBinContent(lcx, lcy));
538  }
539  else
540  {
541  light_yield = 0.0;
542  }
543  }
544  else
545  {
546  if ((lcy >= 1) && (lcy <= m_MapCorrHist[tower_id]->GetNbinsY()) &&
547  (lcx >= 1) && (lcx <= m_MapCorrHist[tower_id]->GetNbinsX()))
548  {
549  light_yield *= (double) (m_MapCorrHist[tower_id]->GetBinContent(lcx, lcy));
550  }
551  else
552  {
553  light_yield = 0.0;
554  }
555  }
556  }
557  else
558  {
559  // old correction (linear ligh yield dependence along r), never tested
560  light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
561  }
562  }
563  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
564  }
565 
566  if (geantino)
567  {
568  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
569  m_Hit->set_eion(-1);
570  }
571  if (edep > 0)
572  {
573  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
574  {
575  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
576  {
577  pp->SetKeep(1); // we want to keep the track
578  }
579  }
580  }
581 
582  // if any of these conditions is true this is the last step in
583  // this volume and we need to save the hit
584  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
585  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
586  // (happens when your detector goes outside world volume)
587  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
588  // aTrack->GetTrackStatus() == fStopAndKill is also set)
589  // aTrack->GetTrackStatus() == fStopAndKill: track ends
590  if (postPoint->GetStepStatus() == fGeomBoundary ||
591  postPoint->GetStepStatus() == fWorldBoundary ||
592  postPoint->GetStepStatus() == fAtRestDoItProc ||
593  aTrack->GetTrackStatus() == fStopAndKill)
594  {
595  // save only hits with energy deposit (or -1 for geantino)
596  if (m_Hit->get_edep() != 0)
597  {
598  m_SaveHitContainer->AddHit(layer_id, m_Hit);
599  if (m_SaveShower)
600  {
602  }
603  // ownership has been transferred to container, set to null
604  // so we will create a new hit for the next track
605  m_Hit = nullptr;
606  }
607  else
608  {
609  // if this hit has no energy deposit, just reset it for reuse
610  // this means we have to delete it in the dtor. If this was
611  // the last hit we processed the memory is still allocated
612  m_Hit->Reset();
613  }
614  }
615  // return true to indicate the hit was used
616  return true;
617  }
618  else
619  {
620  return false;
621  }
622 }
623 
624 //____________________________________________________________________________..
626 {
627  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
628  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
629 
630  // if we do not find the node it's messed up.
631  if (!m_HitContainer)
632  {
633  std::cout << "PHG4OHCalSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
634  }
635  if (!m_AbsorberHitContainer)
636  {
637  if (Verbosity() > 1)
638  {
639  std::cout << "PHG4OHcalSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
640  }
641  }
642 }
643 
644 void PHG4OHCalSteppingAction::FieldChecker(const G4Step* aStep)
645 {
647  assert(se);
648 
649  static const std::string h_field_name = "hOHCalField";
650 
651  if (!se->isHistoRegistered(h_field_name))
652  {
653  TH2F* h = new TH2F(h_field_name.c_str(), "Magnetic field (Tesla) in HCal;X (cm);Y (cm)", 2400,
654  -300, 300, 2400, -300, 300);
655 
656  se->registerHisto(h, 1);
657 
658  std::cout << "PHG4OHCalSteppingAction::FieldChecker - make a histograme to check outer Hcal field map."
659  << " Saved to Fun4AllServer Histo with name " << h_field_name << std::endl;
660  }
661 
662  TH2F* h = dynamic_cast<TH2F*>(se->getHisto(h_field_name));
663  assert(h);
664 
665  assert(aStep);
666  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
667  assert(touch);
668  // get volume of the current step
669  G4VPhysicalVolume* volume = touch->GetVolume();
670  assert(volume);
671 
672  G4ThreeVector globPosition = aStep->GetPreStepPoint()->GetPosition();
673 
674  G4double globPosVec[4] = {0};
675  G4double FieldValueVec[6] = {0};
676 
677  globPosVec[0] = globPosition.x();
678  globPosVec[1] = globPosition.y();
679  globPosVec[2] = globPosition.z();
680  globPosVec[3] = aStep->GetPreStepPoint()->GetGlobalTime();
681 
682  const Int_t binx = h->GetXaxis()->FindBin(globPosVec[0] / cm);
683  const Int_t biny = h->GetYaxis()->FindBin(globPosVec[1] / cm);
684 
685  if (h->GetBinContent(binx, binx) == 0)
686  { // only fille unfilled bins
687 
688  G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager();
689  assert(transportMgr);
690 
691  G4PropagatorInField* fFieldPropagator = transportMgr->GetPropagatorInField();
692  assert(fFieldPropagator);
693 
694  G4FieldManager* fieldMgr = fFieldPropagator->FindAndSetFieldManager(volume);
695  assert(fieldMgr);
696 
697  const G4Field* pField = fieldMgr->GetDetectorField();
698  assert(pField);
699 
700  pField->GetFieldValue(globPosVec, FieldValueVec);
701 
702  G4ThreeVector FieldValue = G4ThreeVector(FieldValueVec[0],
703  FieldValueVec[1], FieldValueVec[2]);
704 
705  const double Bz = FieldValue.z() / tesla;
706 
707  h->SetBinContent(binx, biny, Bz);
708 
709  std::cout << "PHG4OHCalSteppingAction::FieldChecker - "
710  << "volume " << volume->GetName() << " / " << volume->GetLogicalVolume()->GetName()
711  << "\t bin " << binx
712  << ", " << biny << " : Bz= " << Bz << " B = " << FieldValue.mag() / tesla
713  << " Tesla @ x,y = " << globPosVec[0] / cm
714  << "," << globPosVec[1] / cm << " cm" << std::endl;
715  }
716 }
717 
719 {
720  if (type == "G4HIT")
721  {
723  return;
724  }
725  else if (type == "G4HIT_ABSORBER")
726  {
728  return;
729  }
730  std::cout << "Invalid output hit node type " << type << std::endl;
731  gSystem->Exit(1);
732  return;
733 }
734 
736 {
737  PHNodeIterator nodeItr(topNode);
738  PHCompositeNode* dst_node = dynamic_cast<PHCompositeNode*>(
739  nodeItr.findFirst("PHCompositeNode", "DST"));
740  if (!dst_node)
741  {
742  std::cout << "PHComposite node created: DST" << std::endl;
743  dst_node = new PHCompositeNode("DST");
744  topNode->addNode(dst_node);
745  }
746  PHNodeIterator dstiter(dst_node);
747  PHCompositeNode* DetNode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", m_Detector->SuperDetector()));
748  if (!DetNode)
749  {
750  DetNode = new PHCompositeNode(m_Detector->SuperDetector());
751  dst_node->addNode(DetNode);
752  }
753  m_CaloInfoContainer = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::HCAL);
754  PHIODataNode<PHObject>* towerNode = new PHIODataNode<PHObject>(m_CaloInfoContainer, "TOWERINFO_SIM_" + m_Detector->SuperDetector(), "PHObject");
755  DetNode->addNode(towerNode);
756 }