g4simulation/g4ohcal/PHG4OHCalSteppingAction.cc

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