g4simulation/g4caloprototype/PHG4Prototype3InnerHcalSteppingAction.cc

0001 #include "PHG4Prototype3InnerHcalSteppingAction.h"
0002 #include "PHG4Prototype3InnerHcalDetector.h"
0003
0004 #include <phparameter/PHParameters.h>
0005
0006 #include <g4main/PHG4Hit.h>
0007 #include <g4main/PHG4HitContainer.h>
0008 #include <g4main/PHG4Hitv1.h>
0009 #include <g4main/PHG4Shower.h>
0010 #include <g4main/PHG4SteppingAction.h>         // for PHG4SteppingAction
0011 #include <g4main/PHG4TrackUserInfoV1.h>
0012
0013 #include <phool/getClass.h>
0014
0015 #include <Geant4/G4ParticleDefinition.hh>      // for G4ParticleDefinition
0016 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHandle
0017 #include <Geant4/G4Step.hh>
0018 #include <Geant4/G4StepPoint.hh>               // for G4StepPoint
0019 #include <Geant4/G4StepStatus.hh>              // for fGeomBoundary, fAtRest...
0020 #include <Geant4/G4String.hh>                  // for G4String
0021 #include <Geant4/G4SystemOfUnits.hh>
0022 #include <Geant4/G4ThreeVector.hh>             // for G4ThreeVector
0023 #include <Geant4/G4TouchableHandle.hh>         // for G4TouchableHandle
0024 #include <Geant4/G4Track.hh>                   // for G4Track
0025 #include <Geant4/G4TrackStatus.hh>             // for fStopAndKill
0026 #include <Geant4/G4Types.hh>                   // for G4double
0027 #include <Geant4/G4VPhysicalVolume.hh>         // for G4VPhysicalVolume
0028 #include <Geant4/G4VTouchable.hh>              // for G4VTouchable
0029 #include <Geant4/G4VUserTrackInformation.hh>   // for G4VUserTrackInformation
0030
0031 #include <iostream>
0032 #include <string>                              // for string, operator+, ope...
0033
0034 class PHCompositeNode;
0035
0036 //#define TESTSINGLESLAT
0037 #ifdef TESTSINGLESLAT
0038 static const int nrow = 3;
0039 static const int nslat = 10;
0040 #endif
0041
0042 using namespace std;
0043 //____________________________________________________________________________..
0044 PHG4Prototype3InnerHcalSteppingAction::PHG4Prototype3InnerHcalSteppingAction(PHG4Prototype3InnerHcalDetector* detector, const PHParameters* parameters)
0045   :   PHG4SteppingAction(detector->GetName())
0046   , m_Detector(detector)
0047   , m_HitContainer(nullptr)
0048   , m_AbsorberHitContainer(nullptr)
0049   , m_Hit(nullptr)
0050   , m_Params(parameters)
0051   , m_SaveHitContainer(nullptr)
0052   , m_SaveShower(nullptr)
0053   , m_AbsorberTruth(m_Params->get_int_param("absorbertruth"))
0054   , m_IsActive(m_Params->get_int_param("active"))
0055   , m_IsBlackHole(m_Params->get_int_param("blackhole"))
0056   , m_LightScintModel(m_Params->get_int_param("light_scint_model"))
0057 {
0058   SetLightCorrection(m_Params->get_double_param("light_balance_inner_radius") * cm,
0059                      m_Params->get_double_param("light_balance_inner_corr"),
0060              m_Params->get_double_param("light_balance_outer_radius") * cm,
0061              m_Params->get_double_param("light_balance_outer_corr"));
0062 }
0063
0064 PHG4Prototype3InnerHcalSteppingAction::~PHG4Prototype3InnerHcalSteppingAction()
0065 {
0066   // if the last hit was a zero energie deposit hit, it is just reset
0067   // and the memory is still allocated, so we need to delete it here
0068   // if the last hit was saved, hit is a nullptr pointer which are
0069   // legal to delete (it results in a no operation)
0070   delete m_Hit;
0071 }
0072
0073 //____________________________________________________________________________..
0074 bool PHG4Prototype3InnerHcalSteppingAction::UserSteppingAction(const G4Step* aStep, bool)
0075 {
0076   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0077   // get volume of the current step
0078   G4VPhysicalVolume* volume = touch->GetVolume();
0079
0080   // m_Detector->IsInPrototype3InnerHcal(volume)
0081   // returns
0082   //  0 is outside of Prototype3InnerHcal
0083   //  1 is inside scintillator
0084   // -1 is steel absorber
0085
0086   int whichactive = m_Detector->IsInPrototype3InnerHcal(volume);
0087
0088   if (!whichactive)
0089   {
0090     return false;
0091   }
0092   int row_id = -1;
0093   int slat_id = -1;
0094   if (whichactive > 0)  // scintillator
0095   {
0096     // The slat id is just the copy number
0097     slat_id = volume->GetCopyNo();
0098     // the row id comes from saved info in the detector construction
0099     G4VPhysicalVolume* mothervolume = touch->GetVolume(1);
0100     row_id = m_Detector->get_scinti_row_id(mothervolume->GetName());
0101     // cout << "mother volume: " <<  mothervolume->GetName()
0102     //   << ", volume name " << volume->GetName() << ", row: " << row_id
0103     //   << ", column: " << slat_id << endl;
0104 #ifdef TESTSINGLESLAT
0105     if (row_id != nrow)
0106     {
0107       return false;
0108     }
0109     if (slat_id != nslat)
0110     {
0111       return false;
0112     }
0113 #endif
0114   }
0115   else
0116   {
0117     // the row id comes from saved info in the detector construction
0118     row_id = m_Detector->get_steel_plate_id(volume->GetName());
0119 #ifdef TESTSINGLESLAT
0120     return false;
0121 #endif
0122   }
0123   // collect energy and track length step by step
0124   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0125   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0126   G4double light_yield = 0;
0127   const G4Track* aTrack = aStep->GetTrack();
0128
0129   // if this block stops everything, just put all kinetic energy into edep
0130   if (m_IsBlackHole)
0131   {
0132     edep = aTrack->GetKineticEnergy() / GeV;
0133     G4Track* killtrack = const_cast<G4Track*>(aTrack);
0134     killtrack->SetTrackStatus(fStopAndKill);
0135   }
0136   int layer_id = m_Detector->get_Layer();
0137   // make sure we are in a volume
0138   if (m_IsActive)
0139   {
0140     bool geantino = false;
0141
0142     // the check for the pdg code speeds things up, I do not want to make
0143     // an expensive string compare for every track when we know
0144     // geantino or chargedgeantino has pid=0
0145     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0146         aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
0147     {
0148       geantino = true;
0149     }
0150     G4StepPoint* prePoint = aStep->GetPreStepPoint();
0151     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0152     //       cout << "track id " << aTrack->GetTrackID() << endl;
0153     //       cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
0154     //       cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
0155     switch (prePoint->GetStepStatus())
0156     {
0157     case fGeomBoundary:
0158     case fUndefined:
0159       if (!m_Hit)
0160       {
0161         m_Hit = new PHG4Hitv1();
0162       }
0163       m_Hit->set_row(row_id);
0164       if (whichactive > 0)  // only for scintillators
0165       {
0166         m_Hit->set_scint_id(slat_id);  // the slat id in the mother volume (or steel plate id), the column
0167       }
0168       //here we set the entrance values in cm
0169       m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0170       m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0171       m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0172       // time in ns
0173       m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0174       //set the track ID
0175       m_Hit->set_trkid(aTrack->GetTrackID());
0176       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0177       {
0178         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0179         {
0180           m_Hit->set_trkid(pp->GetUserTrackId());
0181           m_Hit->set_shower_id(pp->GetShower()->get_id());
0182         }
0183       }
0184
0185       //set the initial energy deposit
0186       m_Hit->set_edep(0);
0187
0188       m_Hit->set_hit_type(0);
0189       if ((aTrack->GetParticleDefinition()->GetParticleName().find("e+") != string::npos) ||
0190           (aTrack->GetParticleDefinition()->GetParticleName().find("e-") != string::npos))
0191         m_Hit->set_hit_type(1);
0192
0193       if (whichactive > 0)  // return of IsInPrototype3InnerHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
0194       {
0195         m_SaveHitContainer = m_HitContainer;
0196         m_Hit->set_light_yield(0);  // for scintillator only, initialize light yields
0197         m_Hit->set_eion(0);
0198       }
0199       else
0200       {
0201         m_SaveHitContainer = m_AbsorberHitContainer;
0202       }
0203
0204       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0205       {
0206         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0207         {
0208           m_Hit->set_trkid(pp->GetUserTrackId());
0209           m_Hit->set_shower_id(pp->GetShower()->get_id());
0210           m_SaveShower = pp->GetShower();
0211         }
0212       }
0213       break;
0214     default:
0215       break;
0216     }
0217     // here we just update the exit values, it will be overwritten
0218     // for every step until we leave the volume or the particle
0219     // ceases to exist
0220     m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0221     m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0222     m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0223
0224     m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0225
0226     if (whichactive > 0)  // return of IsInPrototype3InnerHcalDetector, > 0 hit in scintillator, < 0 hit in absorber
0227     {
0228       m_Hit->set_eion(m_Hit->get_eion() + eion);
0229       light_yield = eion;
0230       if (m_LightScintModel)
0231       {
0232         light_yield = GetVisibleEnergyDeposition(aStep);  // for scintillator only, calculate light yields
0233       }
0234 // position correction for light yield, code is in base class
0235 // returns 1 if not implemented
0236       light_yield = light_yield *GetLightCorrection(postPoint->GetPosition().x()/cm,postPoint->GetPosition().y()/cm) ;
0237
0238       m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
0239     }
0240
0241     //sum up the energy to get total deposited
0242     m_Hit->set_edep(m_Hit->get_edep() + edep);
0243     if (geantino)
0244     {
0245       m_Hit->set_edep(-1);    // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0246       if (whichactive > 0)  // add light yield for scintillators
0247       {
0248         m_Hit->set_light_yield(-1);
0249         m_Hit->set_eion(-1);
0250       }
0251     }
0252     if (edep > 0 && (whichactive > 0 || m_AbsorberTruth > 0))
0253     {
0254       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0255       {
0256         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0257         {
0258           pp->SetKeep(1);  // we want to keep the track
0259         }
0260       }
0261     }
0262     // if any of these conditions is true this is the last step in
0263     // this volume and we need to save the hit
0264     // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0265     // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0266     // (happens when your detector goes outside world volume)
0267     // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0268     // aTrack->GetTrackStatus() == fStopAndKill is also set)
0269     // aTrack->GetTrackStatus() == fStopAndKill: track ends
0270     if (postPoint->GetStepStatus() == fGeomBoundary ||
0271         postPoint->GetStepStatus() == fWorldBoundary ||
0272         postPoint->GetStepStatus() == fAtRestDoItProc ||
0273         aTrack->GetTrackStatus() == fStopAndKill)
0274     {
0275       // save only hits with energy deposit (or -1 for geantino)
0276       if (m_Hit->get_edep())
0277       {
0278         m_SaveHitContainer->AddHit(layer_id, m_Hit);
0279         if (m_SaveShower)
0280         {
0281           m_SaveShower->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
0282         }
0283         // ownership has been transferred to container, set to null
0284         // so we will create a new hit for the next track
0285         m_Hit = nullptr;
0286       }
0287       else
0288       {
0289         // if this hit has no energy deposit, just reset it for reuse
0290         // this means we have to delete it in the dtor. If this was
0291         // the last hit we processed the memory is still allocated
0292         m_Hit->Reset();
0293       }
0294     }
0295
0296     //       m_Hit->identify();
0297     // return true to indicate the hit was used
0298     return true;
0299   }
0300   else
0301   {
0302     return false;
0303   }
0304 }
0305
0306 //____________________________________________________________________________..
0307 void PHG4Prototype3InnerHcalSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0308 {
0309   string hitnodename;
0310   string absorbernodename;
0311   if (m_Detector->SuperDetector() != "NONE")
0312   {
0313     hitnodename = "G4HIT_" + m_Detector->SuperDetector();
0314     absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
0315   }
0316   else
0317   {
0318     hitnodename = "G4HIT_" + m_Detector->GetName();
0319     absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
0320   }
0321
0322   //now look for the map and grab a pointer to it.
0323   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
0324   m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
0325
0326   // if it is active and we do not find the node it's messed up.
0327   if (!m_HitContainer && m_IsActive)
0328   {
0329     std::cout << "PHG4Prototype3InnerHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
0330   }
0331   if (!m_AbsorberHitContainer)
0332   {
0333     if (Verbosity() > 1)
0334     {
0335       cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
0336     }
0337   }
0338 }