g4simulation/g4detectors/PHG4HcalSteppingAction.cc

0001 #include "PHG4HcalSteppingAction.h"
0002 #include "PHG4HcalDetector.h"
0003
0004 #include <g4main/PHG4Hit.h>
0005 #include <g4main/PHG4HitContainer.h>
0006 #include <g4main/PHG4Hitv1.h>
0007 #include <g4main/PHG4Shower.h>
0008 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0009
0010 #include <g4main/PHG4TrackUserInfoV1.h>
0011
0012 #include <phool/getClass.h>
0013
0014 #include <Geant4/G4ParticleDefinition.hh>  // for G4ParticleDefinition
0015 #include <Geant4/G4Step.hh>
0016 #include <Geant4/G4StepPoint.hh>              // for G4StepPoint
0017 #include <Geant4/G4StepStatus.hh>             // for fGeomBoundary, fAtRestD...
0018 #include <Geant4/G4String.hh>                 // for G4String
0019 #include <Geant4/G4SystemOfUnits.hh>          // for cm, GeV, nanosecond
0020 #include <Geant4/G4ThreeVector.hh>            // for G4ThreeVector
0021 #include <Geant4/G4TouchableHandle.hh>        // for G4TouchableHandle
0022 #include <Geant4/G4Track.hh>                  // for G4Track
0023 #include <Geant4/G4TrackStatus.hh>            // for fStopAndKill
0024 #include <Geant4/G4Types.hh>                  // for G4double
0025 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0026 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0027
0028 #include <iostream>
0029 #include <string>  // for char_traits, operator<<
0030
0031 class G4VPhysicalVolume;
0032 class PHCompositeNode;
0033
0034 //____________________________________________________________________________..
0035 PHG4HcalSteppingAction::PHG4HcalSteppingAction(PHG4HcalDetector* detector)
0036   : PHG4SteppingAction(detector->GetName())
0037   , detector_(detector)
0038 {
0039 }
0040
0041 //____________________________________________________________________________..
0042 PHG4HcalSteppingAction::~PHG4HcalSteppingAction()
0043 {
0044   // if the last hit was a zero energie deposit hit, it is just reset
0045   // and the memory is still allocated, so we need to delete it here
0046   // if the last hit was saved, hit is a nullptr pointer which are
0047   // legal to delete (it results in a no operation)
0048   delete m_Hit;
0049 }
0050
0051 //____________________________________________________________________________..
0052 bool PHG4HcalSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0053 {
0054   // get volume of the current step
0055   G4VPhysicalVolume* volume = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
0056
0057   // collect energy and track length step by step
0058   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0059   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0060   G4double light_yield = 0;
0061
0062   const G4Track* aTrack = aStep->GetTrack();
0063
0064   int layer_id = detector_->get_Layer();
0065   // make sure we are in a volume
0066   // IsInCylinderActive returns the number of the scintillator
0067   // slat which has fired
0068   int isactive = detector_->IsInCylinderActive(volume);
0069   if (isactive > PHG4HcalDetector::INACTIVE)
0070   {
0071     bool geantino = false;
0072     // the check for the pdg code speeds things up, I do not want to make
0073     // an expensive string compare for every track when we know
0074     // geantino or chargedgeantino has pid=0
0075     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0076         aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
0077     {
0078       geantino = true;
0079     }
0080     G4StepPoint* prePoint = aStep->GetPreStepPoint();
0081     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0082     //       std::cout << "track id " << aTrack->GetTrackID() << std::endl;
0083     //        std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
0084     //        std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
0085     switch (prePoint->GetStepStatus())
0086     {
0087     case fGeomBoundary:
0088     case fUndefined:
0089
0090       if (!m_Hit)
0091       {
0092         m_Hit = new PHG4Hitv1();
0093       }
0094       m_Hit->set_layer((unsigned int) layer_id);
0095       m_Hit->set_scint_id(isactive);  // isactive contains the scintillator slat id
0096       // here we set the entrance values in cm
0097       m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0098       m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0099       m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0100
0101       // time in ns
0102       m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0103       // set the track ID
0104       m_Hit->set_trkid(aTrack->GetTrackID());
0105       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0106       {
0107         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0108         {
0109           m_Hit->set_trkid(pp->GetUserTrackId());
0110           m_Hit->set_shower_id(pp->GetShower()->get_id());
0111           m_SaveShower = pp->GetShower();
0112         }
0113       }
0114
0115       // set the initial energy deposit
0116       m_Hit->set_edep(0);
0117       m_Hit->set_eion(0);  // only implemented for v5 otherwise empty
0118       m_Hit->set_light_yield(0);
0119       //      m_Hit->print();
0120       // Now add the hit
0121       if (isactive >= 0)  // the slat ids start with zero
0122       {
0123         m_SaveHitContainer = m_HitContainer;
0124         //        unsigned int shift_layer_id = layer_id << (phg4hitdefs::keybits - 3);
0125       }
0126       else
0127       {
0128         m_SaveHitContainer = m_AbsorberHits;
0129       }
0130       //      m_SaveHitContainer->AddHit(layer_id, m_Hit);
0131
0132       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0133       {
0134         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0135         {
0136           pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
0137         }
0138       }
0139       if (m_Hit->get_z(0) > zmax || m_Hit->get_z(0) < zmin)
0140       {
0141         std::cout << "PHG4HcalSteppingAction: hit outside acceptance, layer: " << layer_id << std::endl;
0142         m_Hit->identify();
0143       }
0144       break;
0145     default:
0146       break;
0147     }
0148     // here we just update the exit values, it will be overwritten
0149     // for every step until we leave the volume or the particle
0150     // ceases to exist
0151     m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0152     m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0153     m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0154
0155     m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0156
0157     if (isactive >= 0)  // the slat ids start with zero
0158     {
0159       if (light_scint_model_)
0160       {
0161         light_yield = GetVisibleEnergyDeposition(aStep);
0162       }
0163       else
0164       {
0165         light_yield = eion;
0166       }
0167
0168       if (ValidCorrection())
0169       {
0170         light_yield = light_yield * GetLightCorrection(postPoint->GetPosition().x() / cm, (postPoint->GetPosition().y() / cm));
0171       }
0172     }
0173
0174     // sum up the energy to get total deposited
0175     m_Hit->set_edep(m_Hit->get_edep() + edep);
0176     m_Hit->set_eion(m_Hit->get_eion() + eion);
0177     m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
0178     m_Hit->set_path_length(aTrack->GetTrackLength() / cm);
0179
0180     if (m_Hit->get_z(1) > zmax || m_Hit->get_z(1) < zmin)
0181     {
0182       std::cout << "PHG4HcalSteppingAction: hit outside acceptance zmin " << zmin << ", zmax " << zmax << " at exit" << std::endl;
0183       m_Hit->identify();
0184     }
0185     if (geantino)
0186     {
0187       m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0188       m_Hit->set_eion(-1);
0189     }
0190     if (edep > 0)
0191     {
0192       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0193       {
0194         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0195         {
0196           pp->SetKeep(1);  // we want to keep the track
0197         }
0198       }
0199     }
0200     // if any of these conditions is true this is the last step in
0201     // this volume and we need to save the hit
0202     // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0203     // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0204     // (happens when your detector goes outside world volume)
0205     // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0206     // aTrack->GetTrackStatus() == fStopAndKill is also set)
0207     // aTrack->GetTrackStatus() == fStopAndKill: track ends
0208     if (postPoint->GetStepStatus() == fGeomBoundary ||
0209         postPoint->GetStepStatus() == fWorldBoundary ||
0210         postPoint->GetStepStatus() == fAtRestDoItProc ||
0211         aTrack->GetTrackStatus() == fStopAndKill)
0212     {
0213       // save only hits with energy deposit (or -1 for geantino)
0214       if (m_Hit->get_edep())
0215       {
0216         m_SaveHitContainer->AddHit(layer_id, m_Hit);
0217         if (m_SaveShower)
0218         {
0219           m_SaveShower->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
0220         }
0221         // ownership has been transferred to container, set to null
0222         // so we will create a new hit for the next track
0223         m_Hit = nullptr;
0224       }
0225       else
0226       {
0227         // if this hit has no energy deposit, just reset it for reuse
0228         // this means we have to delete it in the dtor. If this was
0229         // the last hit we processed the memory is still allocated
0230         m_Hit->Reset();
0231       }
0232     }
0233     // return true to indicate the hit was used
0234     //       hit->identify();
0235     // return true to indicate the hit was used
0236     return true;
0237   }
0238   else
0239   {
0240     return false;
0241   }
0242 }
0243
0244 //____________________________________________________________________________..
0245 void PHG4HcalSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0246 {
0247   std::string hitnodename;
0248   std::string absorbernodename;
0249   if (detector_->SuperDetector() != "NONE")
0250   {
0251     hitnodename = "G4HIT_" + detector_->SuperDetector();
0252     absorbernodename = "G4HIT_ABSORBER_" + detector_->SuperDetector();
0253   }
0254   else
0255   {
0256     hitnodename = "G4HIT_" + detector_->GetName();
0257     absorbernodename = "G4HIT_ABSORBER_" + detector_->GetName();
0258   }
0259
0260   // now look for the map and grab a pointer to it.
0261   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0262   m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
0263   // if we do not find the node we need to make it.
0264   if (!m_HitContainer)
0265   {
0266     std::cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
0267   }
0268   if (!m_AbsorberHits)
0269   {
0270     if (Verbosity() > 0)
0271     {
0272       std::cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << std::endl;
0273     }
0274   }
0275 }