Analysis Software
Documentation for sPHENIX simulation 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 
13 
15 #include <fun4all/Fun4AllServer.h>
16 
17 #include <g4main/PHG4Hit.h>
19 #include <g4main/PHG4Hitv1.h>
20 #include <g4main/PHG4Shower.h>
21 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
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 <calobase/TowerInfo.h>
32 #include <calobase/TowerInfoContainer.h>
33 #include <calobase/TowerInfoContainerv1.h>
34 #include <calobase/TowerInfoDefs.h>
35 
36 // Root headers
37 #include <TAxis.h> // for TAxis
38 #include <TFile.h>
39 #include <TH2.h>
40 #include <TNamed.h> // for TNamed
41 #include <TSystem.h>
42 
43 // Geant4 headers
44 
45 #include <Geant4/G4Field.hh>
46 #include <Geant4/G4FieldManager.hh>
47 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
48 #include <Geant4/G4PropagatorInField.hh>
49 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
50 #include <Geant4/G4Step.hh>
51 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
52 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
53 #include <Geant4/G4String.hh> // for G4String
54 #include <Geant4/G4SystemOfUnits.hh>
55 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
56 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
57 #include <Geant4/G4Track.hh> // for G4Track
58 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
59 #include <Geant4/G4TransportationManager.hh>
60 #include <Geant4/G4Types.hh> // for G4double
61 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
62 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
63 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
64 
65 // finally system headers
66 #include <cassert>
67 #include <cmath> // for isfinite, sqrt
68 #include <cstdlib>
69 #include <filesystem>
70 #include <iostream>
71 #include <string> // for operator<<, string
72 #include <utility> // for pair
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  // since we have a copy in memory of this one - we need to delete it
105  delete m_MapCorrHist;
106 }
107 
109 {
112  {
113  std::string mappingfilename(m_Params->get_string_param("MapFileName"));
114  if (mappingfilename.empty())
115  {
116  return 0;
117  }
118  std::string url = CDBInterface::instance()->getUrl("OLD_OUTER_HCAL_TILEMAP", mappingfilename);
119  if (!std::filesystem::exists(url))
120  {
121  std::cout << PHWHERE << " Could not locate " << url << std::endl;
122  std::cout << "use empty filename to ignore mapfile" << std::endl;
123  gSystem->Exit(1);
124  }
125  TFile* file = TFile::Open(url.c_str());
126  file->GetObject(m_Params->get_string_param("MapHistoName").c_str(), m_MapCorrHist);
127  if (!m_MapCorrHist)
128  {
129  std::cout << "ERROR: could not find Histogram " << m_Params->get_string_param("MapHistoName") << " in " << m_Params->get_string_param("MapFileName") << std::endl;
130  gSystem->Exit(1);
131  exit(1); // make code checkers which do not know gSystem->Exit() happy
132  }
133  m_MapCorrHist->SetDirectory(nullptr); // rootism: this needs to be set otherwise histo vanished when closing the file
134  file->Close();
135  delete file;
136  }
137  if (!m_doG4Hit)
138  {
139  try
140  {
141  CreateNodeTree(topNode);
142  }
143  catch (std::exception& e)
144  {
145  std::cout << e.what() << std::endl;
147  }
148  if (Verbosity() > 1) topNode->print();
149  }
150  return 0;
151 }
152 
153 //____________________________________________________________________________..
155 {
156  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
157  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
158  // get volume of the current step
159  G4VPhysicalVolume* volume = touch->GetVolume();
160 
161  // m_Detector->IsInIHCal(volume)
162  // returns
163  // 0 is outside of IHCal
164  // 1 is inside scintillator
165  // -1 is steel absorber
166 
167  int whichactive = m_Detector->IsInOuterHcal(volume);
168 
169  if (!whichactive)
170  {
171  return false;
172  }
174  {
175  FieldChecker(aStep);
176  }
177  int layer_id = -1;
178  int tower_id = -1;
179  if (whichactive > 0) // scintillator
180  {
181  std::pair<int, int> layer_tower = m_Detector->GetLayerTowerId(volume);
182  layer_id = layer_tower.first;
183  tower_id = layer_tower.second;
184 
185  // std::cout<<"******** Outer HCal\t"<<volume->GetName()<<"\t"<<layer_id<<"\t"<<tower_id<<std::endl;
186  }
187  else
188  {
189  // absorber hit, do nothing
190  return false;
191  }
192 
193  if (!m_IsActiveFlag) return false;
194 
195  G4StepPoint* prePoint = aStep->GetPreStepPoint();
196  G4StepPoint* postPoint = aStep->GetPostStepPoint();
197  // time window cut
198  double pretime = prePoint->GetGlobalTime() / nanosecond;
199  double posttime = postPoint->GetGlobalTime() / nanosecond;
200  if (posttime < m_tmin || pretime > m_tmax) return false;
201  if ((posttime - pretime) > m_dt) return false;
202  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
203  const G4Track* aTrack = aStep->GetTrack();
204  // we only need visible energy here
205  double light_yield = eion;
206 
207  // correct evis using light map
209  {
210  light_yield = GetVisibleEnergyDeposition(aStep);
211  if (m_MapCorrHist)
212  {
213  const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
214  const G4ThreeVector& worldPosition = postPoint->GetPosition();
215  G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
216  float lx = localPosition.x() / cm;
217  float lz = fabs(localPosition.z() / cm); // reverse the sense for towerid<12
218 
219  // convert to the map bin coordinates:
220  // map is in 0.5 cm bins
221  int lcz = (int) (2.0 * lz) + 1;
222  int lcx = (int) (2.0 * (lx + 42.75)) + 1;
223 
224  if ((lcx >= 1) && (lcx <= m_MapCorrHist->GetNbinsY()) &&
225  (lcz >= 1) && (lcz <= m_MapCorrHist->GetNbinsX()))
226  {
227  light_yield *= m_MapCorrHist->GetBinContent(lcz, lcx);
228  }
229  else
230  {
231  light_yield = 0.0;
232  }
233  }
234  else
235  {
236  // old correction (linear ligh yield dependence along r), never tested
237  light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
238  }
239  }
240  // find the tower index for this step, tower_id is ieta, layer_id/5 is iphi
241  unsigned int ieta = tower_id;
242  unsigned int iphi = (unsigned int) layer_id / 5;
243  unsigned int tower_key = TowerInfoDefs::encode_hcal(ieta, iphi);
245  // set keep for the track
246  if (light_yield > 0)
247  {
248  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
249  {
250  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
251  {
252  pp->SetKeep(1); // we want to keep the track
253  }
254  }
255  }
256 
257  return true;
258 }
259 
260 //____________________________________________________________________________..
261 bool PHG4OuterHcalSteppingAction::UserSteppingAction(const G4Step* aStep, bool)
262 {
263  if ((!m_doG4Hit) && (!m_IsBlackHoleFlag))
264  {
265  return NoHitSteppingAction(aStep);
266  }
267  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
268  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
269  // get volume of the current step
270  G4VPhysicalVolume* volume = touch->GetVolume();
271 
272  // m_Detector->IsInOuterHcal(volume)
273  // returns
274  // 0 is outside of OuterHcal
275  // 1 is inside scintillator
276  // -1 is steel absorber (if absorber set to active)
277 
278  int whichactive = m_Detector->IsInOuterHcal(volume);
279 
280  if (!whichactive)
281  {
282  return false;
283  }
284 
286  {
287  FieldChecker(aStep);
288  }
289 
290  int layer_id = -1;
291  int tower_id = -1;
292  if (whichactive > 0) // scintillator
293  {
294  std::pair<int, int> layer_tower = m_Detector->GetLayerTowerId(volume);
295  layer_id = layer_tower.first;
296  tower_id = layer_tower.second;
297  }
298  else
299  {
300  layer_id = touch->GetCopyNumber(); // steel plate id
301  }
302  // collect energy and track length step by step
303  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
304  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
305  G4double light_yield = 0;
306 
307  const G4Track* aTrack = aStep->GetTrack();
308 
309  // if this block stops everything, just put all kinetic energy into edep
310  if (m_IsBlackHoleFlag)
311  {
312  edep = aTrack->GetKineticEnergy() / GeV;
313  G4Track* killtrack = const_cast<G4Track*>(aTrack);
314  killtrack->SetTrackStatus(fStopAndKill);
315  }
316 
317  // make sure we are in a volume
318  if (m_IsActiveFlag)
319  {
320  bool geantino = false;
321 
322  // the check for the pdg code speeds things up, I do not want to make
323  // an expensive string compare for every track when we know
324  // geantino or chargedgeantino has pid=0
325  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
326  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
327  {
328  geantino = true;
329  }
330  G4StepPoint* prePoint = aStep->GetPreStepPoint();
331  G4StepPoint* postPoint = aStep->GetPostStepPoint();
332  // std::cout << "track id " << aTrack->GetTrackID() << std::endl;
333  // std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
334  // std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
335 
336  switch (prePoint->GetStepStatus())
337  {
338  case fPostStepDoItProc:
339  if (m_SavePostStepStatus != fGeomBoundary)
340  {
341  break;
342  }
343  else
344  {
345  std::cout << GetName() << ": New Hit for " << std::endl;
346  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
347  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
348  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
349  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
350  std::cout << "last track: " << m_SaveTrackId
351  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
352  std::cout << "phys pre vol: " << volume->GetName()
353  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
354  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
355  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
356  }
357  [[fallthrough]];
358  case fGeomBoundary:
359  case fUndefined:
360  if (!m_Hit)
361  {
362  m_Hit = new PHG4Hitv1();
363  }
364  // here we set the entrance values in cm
365  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
366  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
367  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
368 
369  // DEBUG
370  // add the local coordinates
371  // if(whichactive>0){
372 
373  // G4TouchableHandle theTouchable = prePoint->GetTouchableHandle();
374  // G4ThreeVector worldPosition = prePoint->GetPosition();
375  // G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
376 
377  // m_Hit->set_property(PHG4Hit::prop_local_x_0, (float)(localPosition.x()/cm));
378  // m_Hit->set_property(PHG4Hit::prop_local_y_0, (float)(localPosition.y()/cm));
379  // m_Hit->set_property(PHG4Hit::prop_local_z_0, (float)(localPosition.z()/cm));
380 
381  // m_Hit->set_property(PHG4Hit::prop_layer, (unsigned int) layer_id);
382 
383  // }
384 
385  // time in ns
386  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
387  // set the track ID
388  m_Hit->set_trkid(aTrack->GetTrackID());
389  m_SaveTrackId = aTrack->GetTrackID();
390  // set the initial energy deposit
391  m_Hit->set_edep(0);
392  if (whichactive > 0) // return of IsInOuterHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
393  {
394  m_Hit->set_scint_id(tower_id); // the slat id
395  m_Hit->set_eion(0);
396  m_Hit->set_raw_light_yield(0); // for scintillator only, initialize light yields
397  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
398  // Now save the container we want to add this hit to
400  }
401  else
402  {
404  }
405  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
406  {
407  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
408  {
409  m_Hit->set_trkid(pp->GetUserTrackId());
410  m_Hit->set_shower_id(pp->GetShower()->get_id());
411  m_SaveShower = pp->GetShower();
412  }
413  }
414  break;
415  default:
416  break;
417  }
418  // some sanity checks for inconsistencies
419  // check if this hit was created, if not print out last post step status
420  if (!m_Hit || !isfinite(m_Hit->get_x(0)))
421  {
422  std::cout << GetName() << ": hit was not created" << std::endl;
423  std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
424  << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
425  << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
426  << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
427  std::cout << "last track: " << m_SaveTrackId
428  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
429  std::cout << "phys pre vol: " << volume->GetName()
430  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
431  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
432  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
433  gSystem->Exit(1);
434  }
435  m_SavePostStepStatus = postPoint->GetStepStatus();
436  // check if track id matches the initial one when the hit was created
437  if (aTrack->GetTrackID() != m_SaveTrackId)
438  {
439  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
440  std::cout << "saved track: " << m_SaveTrackId
441  << ", current trackid: " << aTrack->GetTrackID()
442  << std::endl;
443  gSystem->Exit(1);
444  }
445  m_SavePreStepStatus = prePoint->GetStepStatus();
446  m_SavePostStepStatus = postPoint->GetStepStatus();
447  m_SaveVolPre = volume;
448  m_SaveVolPost = touchpost->GetVolume();
449  // here we just update the exit values, it will be overwritten
450  // for every step until we leave the volume or the particle
451  // ceases to exist
452  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
453  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
454  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
455 
456  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
457 
458  // sum up the energy to get total deposited
459  m_Hit->set_edep(m_Hit->get_edep() + edep);
460 
461  if (whichactive > 0)
462  {
463  m_Hit->set_eion(m_Hit->get_eion() + eion);
464  light_yield = eion;
466  {
467  light_yield = GetVisibleEnergyDeposition(aStep);
469  if (m_MapCorrHist)
470  {
471  const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
472  const G4ThreeVector& worldPosition = postPoint->GetPosition();
473  G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
474  float lx = localPosition.x() / cm;
475  float lz = fabs(localPosition.z() / cm); // reverse the sense for towerid<12
476 
477  // convert to the map bin coordinates:
478  // map is in 0.5 cm bins
479  int lcz = (int) (2.0 * lz) + 1;
480  int lcx = (int) (2.0 * (lx + 42.75)) + 1;
481 
482  if ((lcx >= 1) && (lcx <= m_MapCorrHist->GetNbinsY()) &&
483  (lcz >= 1) && (lcz <= m_MapCorrHist->GetNbinsX()))
484  {
485  light_yield *= m_MapCorrHist->GetBinContent(lcz, lcx);
486  }
487  else
488  {
489  light_yield = 0.0;
490  }
491  }
492  else
493  {
494  // old correction (linear ligh yield dependence along r), never tested
495  light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
496  }
497  }
498  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
499  }
500  if (geantino)
501  {
502  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
503  m_Hit->set_eion(-1);
504  }
505  if (edep > 0)
506  {
507  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
508  {
509  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
510  {
511  pp->SetKeep(1); // we want to keep the track
512  }
513  }
514  }
515 
516  // if any of these conditions is true this is the last step in
517  // this volume and we need to save the hit
518  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
519  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
520  // (happens when your detector goes outside world volume)
521  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
522  // aTrack->GetTrackStatus() == fStopAndKill is also set)
523  // aTrack->GetTrackStatus() == fStopAndKill: track ends
524  if (postPoint->GetStepStatus() == fGeomBoundary ||
525  postPoint->GetStepStatus() == fWorldBoundary ||
526  postPoint->GetStepStatus() == fAtRestDoItProc ||
527  aTrack->GetTrackStatus() == fStopAndKill)
528  {
529  // save only hits with energy deposit (or -1 for geantino)
530  if (m_Hit->get_edep() != 0)
531  {
532  m_SaveHitContainer->AddHit(layer_id, m_Hit);
533  if (m_SaveShower)
534  {
536  }
537  // ownership has been transferred to container, set to null
538  // so we will create a new hit for the next track
539  m_Hit = nullptr;
540  }
541  else
542  {
543  // if this hit has no energy deposit, just reset it for reuse
544  // this means we have to delete it in the dtor. If this was
545  // the last hit we processed the memory is still allocated
546  m_Hit->Reset();
547  }
548  }
549  // return true to indicate the hit was used
550  return true;
551  }
552  else
553  {
554  return false;
555  }
556 }
557 
558 //____________________________________________________________________________..
560 {
561  std::string hitnodename;
562  std::string absorbernodename;
563  if (m_Detector->SuperDetector() != "NONE")
564  {
565  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
566  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
567  }
568  else
569  {
570  hitnodename = "G4HIT_" + m_Detector->GetName();
571  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
572  }
573 
574  // now look for the map and grab a pointer to it.
575  m_Hits = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
576  m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
577 
578  // if we do not find the node it's messed up.
579  if (!m_Hits)
580  {
581  std::cout << "PHG4OuterHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
582  }
583  if (!m_AbsorberHits)
584  {
585  if (Verbosity() > 1)
586  {
587  std::cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << std::endl;
588  }
589  }
590 }
591 
593 {
595  assert(se);
596 
597  static const std::string h_field_name = "hOuterHcalField";
598 
599  if (not se->isHistoRegistered(h_field_name))
600  {
601  TH2* h = new TH2F(h_field_name.c_str(), "Magnetic field (Tesla) in HCal;X (cm);Y (cm)", 2400,
602  -300, 300, 2400, -300, 300);
603 
604  se->registerHisto(h, 1);
605 
606  std::cout << "PHG4OuterHcalSteppingAction::FieldChecker - make a histograme to check outer Hcal field map."
607  << " Saved to Fun4AllServer Histo with name " << h_field_name << std::endl;
608  }
609 
610  TH2* h = dynamic_cast<TH2F*>(se->getHisto(h_field_name));
611  assert(h);
612 
613  assert(aStep);
614  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
615  assert(touch);
616  // get volume of the current step
617  G4VPhysicalVolume* volume = touch->GetVolume();
618  assert(volume);
619 
620  G4ThreeVector globPosition = aStep->GetPreStepPoint()->GetPosition();
621 
622  G4double globPosVec[4] = {0};
623  G4double FieldValueVec[6] = {0};
624 
625  globPosVec[0] = globPosition.x();
626  globPosVec[1] = globPosition.y();
627  globPosVec[2] = globPosition.z();
628  globPosVec[3] = aStep->GetPreStepPoint()->GetGlobalTime();
629 
630  const Int_t binx = h->GetXaxis()->FindBin(globPosVec[0] / cm);
631  const Int_t biny = h->GetYaxis()->FindBin(globPosVec[1] / cm);
632 
633  if (h->GetBinContent(binx, binx) == 0)
634  { // only fille unfilled bins
635 
636  G4TransportationManager* transportMgr =
637  G4TransportationManager::GetTransportationManager();
638  assert(transportMgr);
639 
640  G4PropagatorInField* fFieldPropagator =
641  transportMgr->GetPropagatorInField();
642  assert(fFieldPropagator);
643 
644  G4FieldManager* fieldMgr = fFieldPropagator->FindAndSetFieldManager(volume);
645  assert(fieldMgr);
646 
647  const G4Field* pField = fieldMgr->GetDetectorField();
648  assert(pField);
649 
650  pField->GetFieldValue(globPosVec, FieldValueVec);
651 
652  G4ThreeVector FieldValue = G4ThreeVector(FieldValueVec[0],
653  FieldValueVec[1], FieldValueVec[2]);
654 
655  const double B = FieldValue.mag() / tesla;
656 
657  h->SetBinContent(binx, biny, B);
658 
659  std::cout << "PHG4OuterHcalSteppingAction::FieldChecker - "
660  << "bin " << binx
661  << ", " << biny << " := " << B << " Tesla @ x,y = " << globPosVec[0] / cm
662  << "," << globPosVec[1] / cm << " cm" << std::endl;
663  }
664 }
665 
667 {
668  PHNodeIterator nodeItr(topNode);
669  PHCompositeNode* dst_node = dynamic_cast<PHCompositeNode*>(
670  nodeItr.findFirst("PHCompositeNode", "DST"));
671  if (!dst_node)
672  {
673  std::cout << "PHComposite node created: DST" << std::endl;
674  dst_node = new PHCompositeNode("DST");
675  topNode->addNode(dst_node);
676  }
677  PHNodeIterator dstiter(dst_node);
678  PHCompositeNode* DetNode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", m_Detector->SuperDetector()));
679  if (!DetNode)
680  {
681  DetNode = new PHCompositeNode(m_Detector->SuperDetector());
682  dst_node->addNode(DetNode);
683  }
684  m_CaloInfoContainer = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::HCAL);
685  PHIODataNode<PHObject>* towerNode = new PHIODataNode<PHObject>(m_CaloInfoContainer, "TOWERINFO_SIM_" + m_Detector->SuperDetector(), "PHObject");
686  DetNode->addNode(towerNode);
687 }