g4simulation/g4detectors/PHG4OuterHcalSteppingAction.cc

0001 // local headers in quotes (that is important when using include subdirs!)
0002 #include "PHG4OuterHcalSteppingAction.h"
0003
0004 #include "PHG4HcalDefs.h"
0005 #include "PHG4OuterHcalDetector.h"
0006 #include "PHG4StepStatusDecode.h"
0007
0008 // our own headers in alphabetical order
0009
0010 #include <phparameter/PHParameters.h>
0011
0012 #include <ffamodules/CDBInterface.h>
0013
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015 #include <fun4all/Fun4AllServer.h>
0016
0017 #include <g4main/PHG4Hit.h>
0018 #include <g4main/PHG4HitContainer.h>
0019 #include <g4main/PHG4Hitv1.h>
0020 #include <g4main/PHG4Shower.h>
0021 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0022 #include <g4main/PHG4TrackUserInfoV1.h>
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/PHIODataNode.h>    // for PHIODataNode
0026 #include <phool/PHNode.h>          // for PHNode
0027 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0028 #include <phool/PHObject.h>        // for PHObject
0029 #include <phool/getClass.h>
0030
0031 #include <calobase/TowerInfo.h>
0032 #include <calobase/TowerInfoContainer.h>
0033 #include <calobase/TowerInfoContainerv1.h>
0034 #include <calobase/TowerInfoDefs.h>
0035
0036 // Root headers
0037 #include <TAxis.h>  // for TAxis
0038 #include <TFile.h>
0039 #include <TH2.h>
0040 #include <TNamed.h>  // for TNamed
0041 #include <TSystem.h>
0042
0043 // Geant4 headers
0044
0045 #include <Geant4/G4Field.hh>
0046 #include <Geant4/G4FieldManager.hh>
0047 #include <Geant4/G4ParticleDefinition.hh>  // for G4ParticleDefinition
0048 #include <Geant4/G4PropagatorInField.hh>
0049 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHandle
0050 #include <Geant4/G4Step.hh>
0051 #include <Geant4/G4StepPoint.hh>   // for G4StepPoint
0052 #include <Geant4/G4StepStatus.hh>  // for fGeomBoundary, fAtRest...
0053 #include <Geant4/G4String.hh>      // for G4String
0054 #include <Geant4/G4SystemOfUnits.hh>
0055 #include <Geant4/G4ThreeVector.hh>      // for G4ThreeVector
0056 #include <Geant4/G4TouchableHandle.hh>  // for G4TouchableHandle
0057 #include <Geant4/G4Track.hh>            // for G4Track
0058 #include <Geant4/G4TrackStatus.hh>      // for fStopAndKill
0059 #include <Geant4/G4TransportationManager.hh>
0060 #include <Geant4/G4Types.hh>                  // for G4double
0061 #include <Geant4/G4VPhysicalVolume.hh>        // for G4VPhysicalVolume
0062 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0063 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0064
0065 // finally system headers
0066 #include <cassert>
0067 #include <cmath>  // for isfinite, sqrt
0068 #include <cstdlib>
0069 #include <filesystem>
0070 #include <iostream>
0071 #include <string>   // for operator<<, string
0072 #include <utility>  // for pair
0073
0074 class PHCompositeNode;
0075
0076 //____________________________________________________________________________..
0077 PHG4OuterHcalSteppingAction::PHG4OuterHcalSteppingAction(PHG4OuterHcalDetector* detector, const PHParameters* parameters)
0078   : PHG4SteppingAction(detector->GetName())
0079   , m_Detector(detector)
0080   , m_Params(parameters)
0081   , m_EnableFieldCheckerFlag(m_Params->get_int_param("field_check"))
0082   , m_IsActiveFlag(m_Params->get_int_param("active"))
0083   , m_IsBlackHoleFlag(m_Params->get_int_param("blackhole"))
0084   , m_NScintiPlates(m_Params->get_int_param(PHG4HcalDefs::scipertwr) * m_Params->get_int_param("n_towers"))
0085   , m_LightScintModelFlag(m_Params->get_int_param("light_scint_model"))
0086   , m_doG4Hit(m_Params->get_int_param("saveg4hit"))
0087   , m_tmin(m_Params->get_double_param("tmin"))
0088   , m_tmax(m_Params->get_double_param("tmax"))
0089   , m_dt(m_Params->get_double_param("dt"))
0090 {
0091   SetLightCorrection(m_Params->get_double_param("light_balance_inner_radius") * cm,
0092                      m_Params->get_double_param("light_balance_inner_corr"),
0093                      m_Params->get_double_param("light_balance_outer_radius") * cm,
0094                      m_Params->get_double_param("light_balance_outer_corr"));
0095 }
0096
0097 PHG4OuterHcalSteppingAction::~PHG4OuterHcalSteppingAction()
0098 {
0099   // if the last hit was a zero energie deposit hit, it is just reset
0100   // and the memory is still allocated, so we need to delete it here
0101   // if the last hit was saved, hit is a nullptr pointer which are
0102   // legal to delete (it results in a no operation)
0103   delete m_Hit;
0104   // since we have a copy in memory of this one - we need to delete it
0105   delete m_MapCorrHist;
0106 }
0107
0108 int PHG4OuterHcalSteppingAction::InitWithNode(PHCompositeNode* topNode)
0109 {
0110   m_EnableFieldCheckerFlag = m_Params->get_int_param("field_check");
0111   if (m_LightScintModelFlag)
0112   {
0113     std::string mappingfilename(m_Params->get_string_param("MapFileName"));
0114     if (mappingfilename.empty())
0115     {
0116       return 0;
0117     }
0118     std::string url = CDBInterface::instance()->getUrl("OLD_OUTER_HCAL_TILEMAP", mappingfilename);
0119     if (!std::filesystem::exists(url))
0120     {
0121       std::cout << PHWHERE << " Could not locate " << url << std::endl;
0122       std::cout << "use empty filename to ignore mapfile" << std::endl;
0123       gSystem->Exit(1);
0124     }
0125     TFile* file = TFile::Open(url.c_str());
0126     file->GetObject(m_Params->get_string_param("MapHistoName").c_str(), m_MapCorrHist);
0127     if (!m_MapCorrHist)
0128     {
0129       std::cout << "ERROR: could not find Histogram " << m_Params->get_string_param("MapHistoName") << " in " << m_Params->get_string_param("MapFileName") << std::endl;
0130       gSystem->Exit(1);
0131       exit(1);  // make code checkers which do not know gSystem->Exit() happy
0132     }
0133     m_MapCorrHist->SetDirectory(nullptr);  // rootism: this needs to be set otherwise histo vanished when closing the file
0134     file->Close();
0135     delete file;
0136   }
0137   if (!m_doG4Hit)
0138   {
0139     try
0140     {
0141       CreateNodeTree(topNode);
0142     }
0143     catch (std::exception& e)
0144     {
0145       std::cout << e.what() << std::endl;
0146       return Fun4AllReturnCodes::ABORTRUN;
0147     }
0148     if (Verbosity() > 1)
0149     {
0150       topNode->print();
0151     }
0152   }
0153   return 0;
0154 }
0155
0156 //____________________________________________________________________________..
0157 bool PHG4OuterHcalSteppingAction::NoHitSteppingAction(const G4Step* aStep)
0158 {
0159   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0160   G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
0161   // get volume of the current step
0162   G4VPhysicalVolume* volume = touch->GetVolume();
0163
0164   // m_Detector->IsInIHCal(volume)
0165   // returns
0166   //  0 is outside of IHCal
0167   //  1 is inside scintillator
0168   // -1 is steel absorber
0169
0170   int whichactive = m_Detector->IsInOuterHcal(volume);
0171
0172   if (!whichactive)
0173   {
0174     return false;
0175   }
0176   if (m_EnableFieldCheckerFlag)
0177   {
0178     FieldChecker(aStep);
0179   }
0180   int layer_id = -1;
0181   int tower_id = -1;
0182   if (whichactive > 0)  // scintillator
0183   {
0184     std::pair<int, int> layer_tower = m_Detector->GetLayerTowerId(volume);
0185     layer_id = layer_tower.first;
0186     tower_id = layer_tower.second;
0187
0188     // std::cout<<"******** Outer HCal\t"<<volume->GetName()<<"\t"<<layer_id<<"\t"<<tower_id<<std::endl;
0189   }
0190   else
0191   {
0192     // absorber hit, do nothing
0193     return false;
0194   }
0195
0196   if (!m_IsActiveFlag)
0197   {
0198     return false;
0199   }
0200
0201   G4StepPoint* prePoint = aStep->GetPreStepPoint();
0202   G4StepPoint* postPoint = aStep->GetPostStepPoint();
0203   // time window cut
0204   double pretime = prePoint->GetGlobalTime() / nanosecond;
0205   double posttime = postPoint->GetGlobalTime() / nanosecond;
0206   if (posttime < m_tmin || pretime > m_tmax)
0207   {
0208     return false;
0209   }
0210   if ((posttime - pretime) > m_dt)
0211   {
0212     return false;
0213   }
0214   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0215   const G4Track* aTrack = aStep->GetTrack();
0216   // we only need visible energy here
0217   double light_yield = eion;
0218
0219   // correct evis using light map
0220   if (m_LightScintModelFlag)
0221   {
0222     light_yield = GetVisibleEnergyDeposition(aStep);
0223     if (m_MapCorrHist)
0224     {
0225       const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
0226       const G4ThreeVector& worldPosition = postPoint->GetPosition();
0227       G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
0228       float lx = localPosition.x() / cm;
0229       float lz = fabs(localPosition.z() / cm);  // reverse the sense for towerid<12
0230
0231       // convert to the map bin coordinates:
0232       // map is in 0.5 cm bins
0233       int lcz = (int) (2.0 * lz) + 1;
0234       int lcx = (int) (2.0 * (lx + 42.75)) + 1;
0235
0236       if ((lcx >= 1) && (lcx <= m_MapCorrHist->GetNbinsY()) &&
0237           (lcz >= 1) && (lcz <= m_MapCorrHist->GetNbinsX()))
0238       {
0239         light_yield *= m_MapCorrHist->GetBinContent(lcz, lcx);
0240       }
0241       else
0242       {
0243         light_yield = 0.0;
0244       }
0245     }
0246     else
0247     {
0248       // old correction (linear ligh yield dependence along r), never tested
0249       light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
0250     }
0251   }
0252   // find the tower index for this step, tower_id is ieta, layer_id/5 is iphi
0253   unsigned int ieta = tower_id;
0254   unsigned int iphi = (unsigned int) layer_id / 5;
0255   unsigned int tower_key = TowerInfoDefs::encode_hcal(ieta, iphi);
0256   m_CaloInfoContainer->get_tower_at_key(tower_key)->set_energy(m_CaloInfoContainer->get_tower_at_key(tower_key)->get_energy() + light_yield);
0257   // set keep for the track
0258   if (light_yield > 0)
0259   {
0260     if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0261     {
0262       if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0263       {
0264         pp->SetKeep(1);  // we want to keep the track
0265       }
0266     }
0267   }
0268
0269   return true;
0270 }
0271
0272 //____________________________________________________________________________..
0273 bool PHG4OuterHcalSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0274 {
0275   if ((!m_doG4Hit) && (!m_IsBlackHoleFlag))
0276   {
0277     return NoHitSteppingAction(aStep);
0278   }
0279   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0280   G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
0281   // get volume of the current step
0282   G4VPhysicalVolume* volume = touch->GetVolume();
0283
0284   // m_Detector->IsInOuterHcal(volume)
0285   // returns
0286   //  0 is outside of OuterHcal
0287   //  1 is inside scintillator
0288   // -1 is steel absorber (if absorber set to active)
0289
0290   int whichactive = m_Detector->IsInOuterHcal(volume);
0291
0292   if (!whichactive)
0293   {
0294     return false;
0295   }
0296
0297   if (m_EnableFieldCheckerFlag)
0298   {
0299     FieldChecker(aStep);
0300   }
0301
0302   int layer_id = -1;
0303   int tower_id = -1;
0304   if (whichactive > 0)  // scintillator
0305   {
0306     std::pair<int, int> layer_tower = m_Detector->GetLayerTowerId(volume);
0307     layer_id = layer_tower.first;
0308     tower_id = layer_tower.second;
0309   }
0310   else
0311   {
0312     layer_id = touch->GetCopyNumber();  // steel plate id
0313   }
0314   // collect energy and track length step by step
0315   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0316   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0317   G4double light_yield = 0;
0318
0319   const G4Track* aTrack = aStep->GetTrack();
0320
0321   // if this block stops everything, just put all kinetic energy into edep
0322   if (m_IsBlackHoleFlag)
0323   {
0324     edep = aTrack->GetKineticEnergy() / GeV;
0325     G4Track* killtrack = const_cast<G4Track*>(aTrack);
0326     killtrack->SetTrackStatus(fStopAndKill);
0327   }
0328
0329   // make sure we are in a volume
0330   if (m_IsActiveFlag)
0331   {
0332     bool geantino = false;
0333
0334     // the check for the pdg code speeds things up, I do not want to make
0335     // an expensive string compare for every track when we know
0336     // geantino or chargedgeantino has pid=0
0337     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0338         aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
0339     {
0340       geantino = true;
0341     }
0342     G4StepPoint* prePoint = aStep->GetPreStepPoint();
0343     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0344     //       std::cout << "track id " << aTrack->GetTrackID() << std::endl;
0345     //       std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
0346     //       std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
0347
0348     switch (prePoint->GetStepStatus())
0349     {
0350     case fPostStepDoItProc:
0351       if (m_SavePostStepStatus != fGeomBoundary)
0352       {
0353         break;
0354       }
0355       else
0356       {
0357         std::cout << GetName() << ": New Hit for  " << std::endl;
0358         std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0359                   << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0360                   << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0361                   << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0362         std::cout << "last track: " << m_SaveTrackId
0363                   << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0364         std::cout << "phys pre vol: " << volume->GetName()
0365                   << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0366         std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0367                   << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0368       }
0369       [[fallthrough]];
0370     case fGeomBoundary:
0371     case fUndefined:
0372       if (!m_Hit)
0373       {
0374         m_Hit = new PHG4Hitv1();
0375       }
0376       // here we set the entrance values in cm
0377       m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0378       m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0379       m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0380
0381       // DEBUG
0382       // add the local coordinates
0383       // if(whichactive>0){
0384
0385       //    G4TouchableHandle theTouchable = prePoint->GetTouchableHandle();
0386       //    G4ThreeVector worldPosition = prePoint->GetPosition();
0387       //    G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
0388
0389       //    m_Hit->set_property(PHG4Hit::prop_local_x_0, (float)(localPosition.x()/cm));
0390       //    m_Hit->set_property(PHG4Hit::prop_local_y_0, (float)(localPosition.y()/cm));
0391       //    m_Hit->set_property(PHG4Hit::prop_local_z_0, (float)(localPosition.z()/cm));
0392
0393       //   m_Hit->set_property(PHG4Hit::prop_layer, (unsigned int) layer_id);
0394
0395       // }
0396
0397       // time in ns
0398       m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0399       // set the track ID
0400       m_Hit->set_trkid(aTrack->GetTrackID());
0401       m_SaveTrackId = aTrack->GetTrackID();
0402       // set the initial energy deposit
0403       m_Hit->set_edep(0);
0404       if (whichactive > 0)  // return of IsInOuterHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
0405       {
0406         m_Hit->set_scint_id(tower_id);  // the slat id
0407         m_Hit->set_eion(0);
0408         m_Hit->set_raw_light_yield(0);  //  for scintillator only, initialize light yields
0409         m_Hit->set_light_yield(0);      //  for scintillator only, initialize light yields
0410         // Now save the container we want to add this hit to
0411         m_SaveHitContainer = m_Hits;
0412       }
0413       else
0414       {
0415         m_SaveHitContainer = m_AbsorberHits;
0416       }
0417       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0418       {
0419         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0420         {
0421           m_Hit->set_trkid(pp->GetUserTrackId());
0422           m_Hit->set_shower_id(pp->GetShower()->get_id());
0423           m_SaveShower = pp->GetShower();
0424         }
0425       }
0426       break;
0427     default:
0428       break;
0429     }
0430     // some sanity checks for inconsistencies
0431     // check if this hit was created, if not print out last post step status
0432     if (!m_Hit || !isfinite(m_Hit->get_x(0)))
0433     {
0434       std::cout << GetName() << ": hit was not created" << std::endl;
0435       std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0436                 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0437                 << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0438                 << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0439       std::cout << "last track: " << m_SaveTrackId
0440                 << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0441       std::cout << "phys pre vol: " << volume->GetName()
0442                 << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0443       std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0444                 << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0445       gSystem->Exit(1);
0446     }
0447     m_SavePostStepStatus = postPoint->GetStepStatus();
0448     // check if track id matches the initial one when the hit was created
0449     if (aTrack->GetTrackID() != m_SaveTrackId)
0450     {
0451       std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
0452       std::cout << "saved track: " << m_SaveTrackId
0453                 << ", current trackid: " << aTrack->GetTrackID()
0454                 << std::endl;
0455       gSystem->Exit(1);
0456     }
0457     m_SavePreStepStatus = prePoint->GetStepStatus();
0458     m_SavePostStepStatus = postPoint->GetStepStatus();
0459     m_SaveVolPre = volume;
0460     m_SaveVolPost = touchpost->GetVolume();
0461     // here we just update the exit values, it will be overwritten
0462     // for every step until we leave the volume or the particle
0463     // ceases to exist
0464     m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0465     m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0466     m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0467
0468     m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0469
0470     // sum up the energy to get total deposited
0471     m_Hit->set_edep(m_Hit->get_edep() + edep);
0472
0473     if (whichactive > 0)
0474     {
0475       m_Hit->set_eion(m_Hit->get_eion() + eion);
0476       light_yield = eion;
0477       if (m_LightScintModelFlag)
0478       {
0479         light_yield = GetVisibleEnergyDeposition(aStep);
0480         m_Hit->set_raw_light_yield(m_Hit->get_raw_light_yield() + light_yield);
0481         if (m_MapCorrHist)
0482         {
0483           const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
0484           const G4ThreeVector& worldPosition = postPoint->GetPosition();
0485           G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
0486           float lx = localPosition.x() / cm;
0487           float lz = fabs(localPosition.z() / cm);  // reverse the sense for towerid<12
0488
0489           // convert to the map bin coordinates:
0490           // map is in 0.5 cm bins
0491           int lcz = (int) (2.0 * lz) + 1;
0492           int lcx = (int) (2.0 * (lx + 42.75)) + 1;
0493
0494           if ((lcx >= 1) && (lcx <= m_MapCorrHist->GetNbinsY()) &&
0495               (lcz >= 1) && (lcz <= m_MapCorrHist->GetNbinsX()))
0496           {
0497             light_yield *= m_MapCorrHist->GetBinContent(lcz, lcx);
0498           }
0499           else
0500           {
0501             light_yield = 0.0;
0502           }
0503         }
0504         else
0505         {
0506           // old correction (linear ligh yield dependence along r), never tested
0507           light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
0508         }
0509       }
0510       m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
0511     }
0512     if (geantino)
0513     {
0514       m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0515       m_Hit->set_eion(-1);
0516     }
0517     if (edep > 0)
0518     {
0519       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0520       {
0521         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0522         {
0523           pp->SetKeep(1);  // we want to keep the track
0524         }
0525       }
0526     }
0527
0528     // if any of these conditions is true this is the last step in
0529     // this volume and we need to save the hit
0530     // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0531     // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0532     // (happens when your detector goes outside world volume)
0533     // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0534     // aTrack->GetTrackStatus() == fStopAndKill is also set)
0535     // aTrack->GetTrackStatus() == fStopAndKill: track ends
0536     if (postPoint->GetStepStatus() == fGeomBoundary ||
0537         postPoint->GetStepStatus() == fWorldBoundary ||
0538         postPoint->GetStepStatus() == fAtRestDoItProc ||
0539         aTrack->GetTrackStatus() == fStopAndKill)
0540     {
0541       // save only hits with energy deposit (or -1 for geantino)
0542       if (m_Hit->get_edep() != 0)
0543       {
0544         m_SaveHitContainer->AddHit(layer_id, m_Hit);
0545         if (m_SaveShower)
0546         {
0547           m_SaveShower->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
0548         }
0549         // ownership has been transferred to container, set to null
0550         // so we will create a new hit for the next track
0551         m_Hit = nullptr;
0552       }
0553       else
0554       {
0555         // if this hit has no energy deposit, just reset it for reuse
0556         // this means we have to delete it in the dtor. If this was
0557         // the last hit we processed the memory is still allocated
0558         m_Hit->Reset();
0559       }
0560     }
0561     // return true to indicate the hit was used
0562     return true;
0563   }
0564   else
0565   {
0566     return false;
0567   }
0568 }
0569
0570 //____________________________________________________________________________..
0571 void PHG4OuterHcalSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0572 {
0573   std::string hitnodename;
0574   std::string absorbernodename;
0575   if (m_Detector->SuperDetector() != "NONE")
0576   {
0577     hitnodename = "G4HIT_" + m_Detector->SuperDetector();
0578     absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
0579   }
0580   else
0581   {
0582     hitnodename = "G4HIT_" + m_Detector->GetName();
0583     absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
0584   }
0585
0586   // now look for the map and grab a pointer to it.
0587   m_Hits = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0588   m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
0589
0590   // if we do not find the node it's messed up.
0591   if (!m_Hits)
0592   {
0593     std::cout << "PHG4OuterHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
0594   }
0595   if (!m_AbsorberHits)
0596   {
0597     if (Verbosity() > 1)
0598     {
0599       std::cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << std::endl;
0600     }
0601   }
0602 }
0603
0604 void PHG4OuterHcalSteppingAction::FieldChecker(const G4Step* aStep)
0605 {
0606   Fun4AllServer* se = Fun4AllServer::instance();
0607   assert(se);
0608
0609   static const std::string h_field_name = "hOuterHcalField";
0610
0611   if (not se->isHistoRegistered(h_field_name))
0612   {
0613     TH2* h = new TH2F(h_field_name.c_str(), "Magnetic field (Tesla) in HCal;X (cm);Y (cm)", 2400,
0614                       -300, 300, 2400, -300, 300);
0615
0616     se->registerHisto(h, 1);
0617
0618     std::cout << "PHG4OuterHcalSteppingAction::FieldChecker - make a histograme to check outer Hcal field map."
0619               << " Saved to Fun4AllServer Histo with name " << h_field_name << std::endl;
0620   }
0621
0622   TH2* h = dynamic_cast<TH2F*>(se->getHisto(h_field_name));
0623   assert(h);
0624
0625   assert(aStep);
0626   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0627   assert(touch);
0628   // get volume of the current step
0629   G4VPhysicalVolume* volume = touch->GetVolume();
0630   assert(volume);
0631
0632   G4ThreeVector globPosition = aStep->GetPreStepPoint()->GetPosition();
0633
0634   G4double globPosVec[4] = {0};
0635   G4double FieldValueVec[6] = {0};
0636
0637   globPosVec[0] = globPosition.x();
0638   globPosVec[1] = globPosition.y();
0639   globPosVec[2] = globPosition.z();
0640   globPosVec[3] = aStep->GetPreStepPoint()->GetGlobalTime();
0641
0642   const Int_t binx = h->GetXaxis()->FindBin(globPosVec[0] / cm);
0643   const Int_t biny = h->GetYaxis()->FindBin(globPosVec[1] / cm);
0644
0645   if (h->GetBinContent(binx, binx) == 0)
0646   {  // only fille unfilled bins
0647
0648     G4TransportationManager* transportMgr =
0649         G4TransportationManager::GetTransportationManager();
0650     assert(transportMgr);
0651
0652     G4PropagatorInField* fFieldPropagator =
0653         transportMgr->GetPropagatorInField();
0654     assert(fFieldPropagator);
0655
0656     G4FieldManager* fieldMgr = fFieldPropagator->FindAndSetFieldManager(volume);
0657     assert(fieldMgr);
0658
0659     const G4Field* pField = fieldMgr->GetDetectorField();
0660     assert(pField);
0661
0662     pField->GetFieldValue(globPosVec, FieldValueVec);
0663
0664     G4ThreeVector FieldValue = G4ThreeVector(FieldValueVec[0],
0665                                              FieldValueVec[1], FieldValueVec[2]);
0666
0667     const double B = FieldValue.mag() / tesla;
0668
0669     h->SetBinContent(binx, biny, B);
0670
0671     std::cout << "PHG4OuterHcalSteppingAction::FieldChecker - "
0672               << "bin " << binx
0673               << ", " << biny << " := " << B << " Tesla @ x,y = " << globPosVec[0] / cm
0674               << "," << globPosVec[1] / cm << " cm" << std::endl;
0675   }
0676 }
0677
0678 void PHG4OuterHcalSteppingAction::CreateNodeTree(PHCompositeNode* topNode)
0679 {
0680   PHNodeIterator nodeItr(topNode);
0681   PHCompositeNode* dst_node = dynamic_cast<PHCompositeNode*>(
0682       nodeItr.findFirst("PHCompositeNode", "DST"));
0683   if (!dst_node)
0684   {
0685     std::cout << "PHComposite node created: DST" << std::endl;
0686     dst_node = new PHCompositeNode("DST");
0687     topNode->addNode(dst_node);
0688   }
0689   PHNodeIterator dstiter(dst_node);
0690   PHCompositeNode* DetNode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", m_Detector->SuperDetector()));
0691   if (!DetNode)
0692   {
0693     DetNode = new PHCompositeNode(m_Detector->SuperDetector());
0694     dst_node->addNode(DetNode);
0695   }
0696   m_CaloInfoContainer = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::HCAL);
0697   PHIODataNode<PHObject>* towerNode = new PHIODataNode<PHObject>(m_CaloInfoContainer, "TOWERINFO_SIM_" + m_Detector->SuperDetector(), "PHObject");
0698   DetNode->addNode(towerNode);
0699 }