Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:19:29 

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         else
0412         {
0413           if (aTrack->GetTrackID() == m_SaveTrackId)
0414           {
0415             std::cout << GetName() << ": Bad step status combination for the same track " << std::endl;
0416             std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0417                       << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0418                       << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0419                       << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0420             std::cout << "last track: " << m_SaveTrackId
0421                       << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0422             std::cout << "phys pre vol: " << volume->GetName()
0423                       << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0424             std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0425                       << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0426             gSystem->Exit(1);
0427           }
0428         }
0429       }
0430       else
0431       {
0432         std::cout << GetName() << ": New Hit for  " << std::endl;
0433         std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0434                   << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0435                   << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0436                   << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0437         std::cout << "last track: " << m_SaveTrackId
0438                   << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0439         std::cout << "phys pre vol: " << volume->GetName()
0440                   << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0441         std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0442                   << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0443       }
0444       [[fallthrough]];
0445     case fGeomBoundary:
0446     case fUndefined:
0447       if (!m_Hit)
0448       {
0449         m_Hit = new PHG4Hitv1();
0450       }
0451       // here we set the entrance values in cm
0452       m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0453       m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0454       m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0455 
0456       // time in ns
0457       m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0458       // set the track ID
0459       m_Hit->set_trkid(aTrack->GetTrackID());
0460       m_SaveTrackId = aTrack->GetTrackID();
0461       // set the initial energy deposit
0462       m_Hit->set_edep(0);
0463       if (whichactive > 0)  // return of IsInOHCalDetector, > 0 hit in scintillator, < 0 hit in absorber
0464       {
0465         m_Hit->set_sector(sector_id);   // the sector id
0466         m_Hit->set_scint_id(tower_id);  // the slat id
0467         m_Hit->set_eion(0);
0468         m_Hit->set_raw_light_yield(0);  //  for scintillator only, initialize light yields
0469         m_Hit->set_light_yield(0);      //  for scintillator only, initialize light yields
0470         // Now save the container we want to add this hit to
0471         m_SaveHitContainer = m_HitContainer;
0472       }
0473       else
0474       {
0475         m_SaveHitContainer = m_AbsorberHitContainer;
0476       }
0477       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0478       {
0479         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0480         {
0481           m_Hit->set_trkid(pp->GetUserTrackId());
0482           m_Hit->set_shower_id(pp->GetShower()->get_id());
0483           m_SaveShower = pp->GetShower();
0484         }
0485       }
0486       break;
0487     default:
0488       break;
0489     }
0490     // some sanity checks for inconsistencies
0491     // check if this hit was created, if not print out last post step status
0492     if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
0493     {
0494       std::cout << GetName() << ": hit was not created" << std::endl;
0495       std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0496                 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0497                 << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0498                 << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0499       std::cout << "last track: " << m_SaveTrackId
0500                 << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0501       std::cout << "phys pre vol: " << volume->GetName()
0502                 << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0503       std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0504                 << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0505       gSystem->Exit(1);
0506     }
0507     m_SavePostStepStatus = postPoint->GetStepStatus();
0508     // check if track id matches the initial one when the hit was created
0509     if (aTrack->GetTrackID() != m_SaveTrackId)
0510     {
0511       std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
0512       std::cout << "saved track: " << m_SaveTrackId
0513                 << ", current trackid: " << aTrack->GetTrackID()
0514                 << std::endl;
0515       gSystem->Exit(1);
0516     }
0517     m_SavePreStepStatus = prePoint->GetStepStatus();
0518     m_SavePostStepStatus = postPoint->GetStepStatus();
0519     m_SaveVolPre = volume;
0520     m_SaveVolPost = touchpost->GetVolume();
0521     // here we just update the exit values, it will be overwritten
0522     // for every step until we leave the volume or the particle
0523     // ceases to exist
0524     m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0525     m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0526     m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0527 
0528     m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0529 
0530     // sum up the energy to get total deposited
0531     m_Hit->set_edep(m_Hit->get_edep() + edep);
0532 
0533     if (whichactive > 0)
0534     {
0535       m_Hit->set_eion(m_Hit->get_eion() + eion);
0536       light_yield = eion;
0537       if (m_LightScintModelFlag)
0538       {
0539         light_yield = GetVisibleEnergyDeposition(aStep);
0540         m_Hit->set_raw_light_yield(m_Hit->get_raw_light_yield() + light_yield);  // save raw Birks light yield
0541         if ((m_MapCorrHistChim[tower_id]) || (m_MapCorrHist[tower_id]))
0542         {
0543           const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
0544           const G4ThreeVector& worldPosition = postPoint->GetPosition();
0545           G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
0546           float lx = (localPosition.x() / cm);
0547           float ly = (localPosition.y() / cm);
0548 
0549           // convert to the map bin coordinates:
0550           int lcx = (int) (2.0 * lx) + 1;
0551           int lcy = (int) (2.0 * (ly + 0.5)) + 1;
0552 
0553           if ((sector_id == 29) || (sector_id == 30) || (sector_id == 31))
0554           {
0555             if ((lcy >= 1) && (lcy <= m_MapCorrHistChim[tower_id]->GetNbinsY()) &&
0556                 (lcx >= 1) && (lcx <= m_MapCorrHistChim[tower_id]->GetNbinsX()))
0557             {
0558               light_yield *= (double) (m_MapCorrHistChim[tower_id]->GetBinContent(lcx, lcy));
0559             }
0560             else
0561             {
0562               light_yield = 0.0;
0563             }
0564           }
0565           else
0566           {
0567             if ((lcy >= 1) && (lcy <= m_MapCorrHist[tower_id]->GetNbinsY()) &&
0568                 (lcx >= 1) && (lcx <= m_MapCorrHist[tower_id]->GetNbinsX()))
0569             {
0570               light_yield *= (double) (m_MapCorrHist[tower_id]->GetBinContent(lcx, lcy));
0571             }
0572             else
0573             {
0574               light_yield = 0.0;
0575             }
0576           }
0577         }
0578         else
0579         {
0580           // old correction (linear ligh yield dependence along r), never tested
0581           light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
0582         }
0583       }
0584       m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
0585     }
0586 
0587     if (geantino)
0588     {
0589       m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0590       m_Hit->set_eion(-1);
0591     }
0592     if (edep > 0)
0593     {
0594       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0595       {
0596         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0597         {
0598           pp->SetKeep(1);  // we want to keep the track
0599         }
0600       }
0601     }
0602 
0603     // if any of these conditions is true this is the last step in
0604     // this volume and we need to save the hit
0605     // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0606     // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0607     // (happens when your detector goes outside world volume)
0608     // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0609     // aTrack->GetTrackStatus() == fStopAndKill is also set)
0610     // aTrack->GetTrackStatus() == fStopAndKill: track ends
0611     if (postPoint->GetStepStatus() == fGeomBoundary ||
0612         postPoint->GetStepStatus() == fWorldBoundary ||
0613         postPoint->GetStepStatus() == fAtRestDoItProc ||
0614         aTrack->GetTrackStatus() == fStopAndKill)
0615     {
0616       // save only hits with energy deposit (or -1 for geantino)
0617       if (m_Hit->get_edep() != 0)
0618       {
0619         m_SaveHitContainer->AddHit(layer_id, m_Hit);
0620         if (m_SaveShower)
0621         {
0622           m_SaveShower->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
0623         }
0624         // ownership has been transferred to container, set to null
0625         // so we will create a new hit for the next track
0626         m_Hit = nullptr;
0627       }
0628       else
0629       {
0630         // if this hit has no energy deposit, just reset it for reuse
0631         // this means we have to delete it in the dtor. If this was
0632         // the last hit we processed the memory is still allocated
0633         m_Hit->Reset();
0634       }
0635     }
0636     // return true to indicate the hit was used
0637     return true;
0638   }
0639   else
0640   {
0641     return false;
0642   }
0643 }
0644 
0645 //____________________________________________________________________________..
0646 void PHG4OHCalSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0647 {
0648   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0649   m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
0650 
0651   // if we do not find the node it's messed up.
0652   if (!m_HitContainer)
0653   {
0654     std::cout << "PHG4OHCalSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
0655   }
0656   if (!m_AbsorberHitContainer)
0657   {
0658     if (Verbosity() > 1)
0659     {
0660       std::cout << "PHG4OHcalSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
0661     }
0662   }
0663 }
0664 
0665 void PHG4OHCalSteppingAction::FieldChecker(const G4Step* aStep)
0666 {
0667   Fun4AllServer* se = Fun4AllServer::instance();
0668   assert(se);
0669 
0670   static const std::string h_field_name = "hOHCalField";
0671 
0672   if (!se->isHistoRegistered(h_field_name))
0673   {
0674     TH2F* h = new TH2F(h_field_name.c_str(), "Magnetic field (Tesla) in HCal;X (cm);Y (cm)", 2400,
0675                        -300, 300, 2400, -300, 300);
0676 
0677     se->registerHisto(h, 1);
0678 
0679     std::cout << "PHG4OHCalSteppingAction::FieldChecker - make a histograme to check outer Hcal field map."
0680               << " Saved to Fun4AllServer Histo with name " << h_field_name << std::endl;
0681   }
0682 
0683   TH2F* h = dynamic_cast<TH2F*>(se->getHisto(h_field_name));
0684   assert(h);
0685 
0686   assert(aStep);
0687   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0688   assert(touch);
0689   // get volume of the current step
0690   G4VPhysicalVolume* volume = touch->GetVolume();
0691   assert(volume);
0692 
0693   G4ThreeVector globPosition = aStep->GetPreStepPoint()->GetPosition();
0694 
0695   G4double globPosVec[4] = {0};
0696   G4double FieldValueVec[6] = {0};
0697 
0698   globPosVec[0] = globPosition.x();
0699   globPosVec[1] = globPosition.y();
0700   globPosVec[2] = globPosition.z();
0701   globPosVec[3] = aStep->GetPreStepPoint()->GetGlobalTime();
0702 
0703   const int binx = h->GetXaxis()->FindBin(globPosVec[0] / cm);
0704   const int biny = h->GetYaxis()->FindBin(globPosVec[1] / cm);
0705 
0706   if (h->GetBinContent(binx, binx) == 0)
0707   {  // only fille unfilled bins
0708 
0709     G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager();
0710     assert(transportMgr);
0711 
0712     G4PropagatorInField* fFieldPropagator = transportMgr->GetPropagatorInField();
0713     assert(fFieldPropagator);
0714 
0715     G4FieldManager* fieldMgr = fFieldPropagator->FindAndSetFieldManager(volume);
0716     assert(fieldMgr);
0717 
0718     const G4Field* pField = fieldMgr->GetDetectorField();
0719     assert(pField);
0720 
0721     pField->GetFieldValue(globPosVec, FieldValueVec);
0722 
0723     G4ThreeVector FieldValue = G4ThreeVector(FieldValueVec[0],
0724                                              FieldValueVec[1], FieldValueVec[2]);
0725 
0726     const double Bz = FieldValue.z() / tesla;
0727 
0728     h->SetBinContent(binx, biny, Bz);
0729 
0730     std::cout << "PHG4OHCalSteppingAction::FieldChecker - "
0731               << "volume " << volume->GetName() << " / " << volume->GetLogicalVolume()->GetName()
0732               << "\t bin " << binx
0733               << ", " << biny << " : Bz= " << Bz << " B = " << FieldValue.mag() / tesla
0734               << " Tesla @ x,y = " << globPosVec[0] / cm
0735               << "," << globPosVec[1] / cm << " cm" << std::endl;
0736   }
0737 }
0738 
0739 void PHG4OHCalSteppingAction::SetHitNodeName(const std::string& type, const std::string& name)
0740 {
0741   if (type == "G4HIT")
0742   {
0743     m_HitNodeName = name;
0744     return;
0745   }
0746   else if (type == "G4HIT_ABSORBER")
0747   {
0748     m_AbsorberNodeName = name;
0749     return;
0750   }
0751   std::cout << "Invalid output hit node type " << type << std::endl;
0752   gSystem->Exit(1);
0753   return;
0754 }
0755 
0756 void PHG4OHCalSteppingAction::CreateNodeTree(PHCompositeNode* topNode)
0757 {
0758   PHNodeIterator nodeItr(topNode);
0759   PHCompositeNode* dst_node = dynamic_cast<PHCompositeNode*>(
0760       nodeItr.findFirst("PHCompositeNode", "DST"));
0761   if (!dst_node)
0762   {
0763     std::cout << "PHComposite node created: DST" << std::endl;
0764     dst_node = new PHCompositeNode("DST");
0765     topNode->addNode(dst_node);
0766   }
0767   PHNodeIterator dstiter(dst_node);
0768   PHCompositeNode* DetNode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", m_Detector->SuperDetector()));
0769   if (!DetNode)
0770   {
0771     DetNode = new PHCompositeNode(m_Detector->SuperDetector());
0772     dst_node->addNode(DetNode);
0773   }
0774   m_CaloInfoContainer = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::HCAL);
0775   PHIODataNode<PHObject>* towerNode = new PHIODataNode<PHObject>(m_CaloInfoContainer, "TOWERINFO_SIM_" + m_Detector->SuperDetector(), "PHObject");
0776   DetNode->addNode(towerNode);
0777 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!