g4simulation/g4intt/PHG4InttSteppingAction.cc

0001 #include "PHG4InttSteppingAction.h"
0002 #include "PHG4InttDefs.h"
0003 #include "PHG4InttDetector.h"
0004
0005 #include <phparameter/PHParameters.h>
0006 #include <phparameter/PHParametersContainer.h>
0007
0008 #include <g4main/PHG4Hit.h>
0009 #include <g4main/PHG4HitContainer.h>
0010 #include <g4main/PHG4Hitv1.h>
0011 #include <g4main/PHG4Shower.h>
0012 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0013 #include <g4main/PHG4TrackUserInfoV1.h>
0014
0015 #include <phool/getClass.h>
0016
0017 #include <TSystem.h>
0018
0019 #include <Geant4/G4ParticleDefinition.hh>      // for G4ParticleDefinition
0020 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHandle
0021 #include <Geant4/G4Step.hh>
0022 #include <Geant4/G4StepPoint.hh>   // for G4StepPoint
0023 #include <Geant4/G4StepStatus.hh>  // for fGeomBoundary, fAtRes...
0024 #include <Geant4/G4String.hh>      // for G4String
0025 #include <Geant4/G4SystemOfUnits.hh>
0026 #include <Geant4/G4ThreeVector.hh>
0027 #include <Geant4/G4TouchableHandle.hh>
0028 #include <Geant4/G4Track.hh>                  // for G4Track
0029 #include <Geant4/G4TrackStatus.hh>            // for fStopAndKill
0030 #include <Geant4/G4Types.hh>                  // for G4double
0031 #include <Geant4/G4VPhysicalVolume.hh>        // for G4VPhysicalVolume
0032 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0033 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0034
0035 #include <cassert>  // for assert
0036 #include <iostream>
0037 #include <set>     // for set
0038 #include <string>  // for operator<<, string
0039 #include <tuple>   // for tie, tuple
0040
0041 class PHCompositeNode;
0042
0043 //____________________________________________________________________________..
0044 PHG4InttSteppingAction::PHG4InttSteppingAction(PHG4InttDetector* detector, const PHParametersContainer* parameters, const std::pair<std::vector<std::pair<int, int>>::const_iterator, std::vector<std::pair<int, int>>::const_iterator>& layer_begin_end)
0045   : PHG4SteppingAction(detector->GetName())
0046   , m_Detector(detector)
0047   , m_ParamsContainer(parameters)
0048 {
0049   // loop over layers to get laddertype nd active status for each layer
0050   for (auto layeriter = layer_begin_end.first; layeriter != layer_begin_end.second; ++layeriter)
0051   {
0052     int layer = layeriter->second;
0053     const PHParameters* par = m_ParamsContainer->GetParameters(layer);
0054     m_IsActiveMap[layer] = par->get_int_param("active");
0055     m_IsBlackHoleMap[layer] = par->get_int_param("blackhole");
0056     m_LadderTypeMap.insert(std::make_pair(layer, par->get_int_param("laddertype")));
0057     m_InttToTrackerLayerMap.insert(std::make_pair(layeriter->second, layeriter->first));
0058   }
0059
0060   // Get the parameters for each laddertype
0061   for (int iter : PHG4InttDefs::m_SensorSegmentationSet)
0062   {
0063     const PHParameters* par = m_ParamsContainer->GetParameters(iter);
0064     m_StripYMap.insert(std::make_pair(iter, par->get_double_param("strip_y") * cm));
0065     m_StripZMap.insert(std::make_pair(iter, std::make_pair(par->get_double_param("strip_z_0") * cm, par->get_double_param("strip_z_1") * cm)));
0066     m_nStripsPhiCell.insert(std::make_pair(iter, par->get_int_param("nstrips_phi_cell")));
0067     m_nStripsZSensor.insert(std::make_pair(iter, std::make_pair(par->get_int_param("nstrips_z_sensor_0"), par->get_int_param("nstrips_z_sensor_1"))));
0068   }
0069 }
0070
0071 PHG4InttSteppingAction::~PHG4InttSteppingAction()
0072 {
0073   // if the last hit was a zero energie deposit hit, it is just reset
0074   // and the memory is still allocated, so we need to delete it here
0075   // if the last hit was saved, hit is a nullptr pointer which are
0076   // legal to delete (it results in a no operation)
0077   delete m_Hit;
0078 }
0079
0080 //____________________________________________________________________________..
0081 bool PHG4InttSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0082 {
0083   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0084   // get volume of the current step
0085   G4VPhysicalVolume* volume = touch->GetVolume();
0086   const G4Track* aTrack = aStep->GetTrack();
0087   G4StepPoint* prePoint = aStep->GetPreStepPoint();
0088   G4StepPoint* postPoint = aStep->GetPostStepPoint();
0089
0090   const int whichactive = m_Detector->IsInIntt(volume);
0091
0092   if (!whichactive)
0093   {
0094     return false;
0095   }
0096
0097   if (Verbosity() > 0)
0098   {
0099     std::cout << std::endl
0100               << "PHG4SilicoTrackerSteppingAction::UserSteppingAction for volume name (pre) " << touch->GetVolume()->GetName()
0101               << " volume name (1) " << touch->GetVolume(1)->GetName()
0102               << " volume->GetTranslation " << touch->GetVolume()->GetTranslation()
0103               << std::endl;
0104   }
0105
0106   // set ladder index
0107   int sphxlayer = 0;
0108   int inttlayer = 0;
0109   int ladderz = 0;
0110   int ladderphi = 0;
0111   int zposneg = 0;
0112
0113   if (whichactive > 0)  // silicon active sensor
0114   {
0115     // Get the layer and ladder information which is step up in the volume hierarchy
0116     // the ladder also contains inactive volumes but we check in m_Detector->IsInIntt(volume)
0117     // if we are in an active logical volume whioch is located in this ladder
0118     auto iter = m_Detector->get_ActiveVolumeTuple(touch->GetVolume(1));
0119     std::tie(inttlayer, ladderz, ladderphi, zposneg) = iter->second;
0120     if (Verbosity() > 0)
0121     {
0122       std::cout << "     inttlayer " << inttlayer << " ladderz_base " << ladderz << " ladderphi " << ladderphi << " zposneg " << zposneg << std::endl;
0123     }
0124     if (inttlayer < 0 || inttlayer > 7)
0125     {
0126       assert(!"PHG4InttSteppingAction: check Intt ladder layer.");
0127     }
0128     sphxlayer = m_InttToTrackerLayerMap.find(inttlayer)->second;
0129     std::map<int, int>::const_iterator activeiter = m_IsActiveMap.find(inttlayer);
0130     if (activeiter == m_IsActiveMap.end())
0131     {
0132       std::cout << "PHG4InttSteppingAction: could not find active flag for layer " << inttlayer << std::endl;
0133       gSystem->Exit(1);
0134     }
0135     if (activeiter->second == 0)
0136     {
0137       return false;
0138     }
0139   }
0140   else  // whichactive < 0, silicon inactive area, FPHX, stabe etc. as absorbers
0141   {
0142     auto iter = m_Detector->get_PassiveVolumeTuple(touch->GetVolume(0)->GetLogicalVolume());
0143     std::tie(inttlayer, ladderz) = iter->second;
0144     sphxlayer = inttlayer;  // for absorber we use the Intt layer, not the tracking layer in sPHENIX
0145   }                         // end of si inactive area block
0146
0147   // collect energy and track length step by step
0148   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0149   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0150
0151   // if this block stops everything, just put all kinetic energy into edep
0152   if ((m_IsBlackHoleMap.find(inttlayer))->second == 1)
0153   {
0154     edep = aTrack->GetKineticEnergy() / GeV;
0155     G4Track* killtrack = const_cast<G4Track*>(aTrack);
0156     killtrack->SetTrackStatus(fStopAndKill);
0157   }
0158
0159   bool geantino = false;
0160
0161   // the check for the pdg code speeds things up, I do not want to make
0162   // an expensive string compare for every track when we know
0163   // geantino or chargedgeantino has pid=0
0164   if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
0165   {
0166     geantino = true;
0167   }
0168
0169   if (Verbosity() > 1)
0170   {
0171     std::cout << " aTrack->GetTrackID " << aTrack->GetTrackID() << " aTrack->GetParentID " << aTrack->GetParentID()
0172               << " Intt layer " << inttlayer
0173               << " prePoint step status = " << prePoint->GetStepStatus() << " postPoint step status = " << postPoint->GetStepStatus()
0174               << std::endl;
0175   }
0176   switch (prePoint->GetStepStatus())
0177   {
0178   case fGeomBoundary:
0179   case fUndefined:
0180
0181     if (Verbosity() > 1)
0182     {
0183       std::cout << "  start new g4hit for:  aTrack->GetParentID " << aTrack->GetParentID() << " aTrack->GetTrackID " << aTrack->GetTrackID() << " Intt layer " << inttlayer
0184                 << "   prePoint step status  = " << prePoint->GetStepStatus() << " postPoint->GetStepStatus = " << postPoint->GetStepStatus() << std::endl;
0185     }
0186     // if previous hit was saved, hit pointer was set to nullptr
0187     // and we have to make a new one
0188     if (!m_Hit)
0189     {
0190       m_Hit = new PHG4Hitv1();
0191     }
0192
0193     // set the index values needed to locate the sensor strip
0194     if (zposneg == 1)
0195     {
0196       ladderz += 2;  // ladderz = 0, 1 for negative z and = 2, 3 for positive z
0197     }
0198     if (Verbosity() > 0)
0199     {
0200       std::cout << "     ladderz = " << ladderz << std::endl;
0201     }
0202
0203     m_Hit->set_ladder_z_index(ladderz);
0204
0205     if (whichactive > 0)
0206     {
0207       m_Hit->set_layer(sphxlayer);
0208       m_Hit->set_ladder_phi_index(ladderphi);
0209       m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
0210       m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
0211       m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
0212       m_Hit->set_eion(0);
0213     }
0214
0215     // here we set the entrance values in cm
0216     m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0217     m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0218     m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0219
0220     StoreLocalCoordinate(m_Hit, aStep, true, false);
0221
0222     if (Verbosity() > 0)
0223     {
0224       std::cout << "     prePoint hit position x,y,z = " << prePoint->GetPosition().x() / cm
0225                 << "    " << prePoint->GetPosition().y() / cm
0226                 << "     " << prePoint->GetPosition().z() / cm
0227                 << std::endl;
0228     }
0229
0230     // time in ns
0231     m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0232
0233     // set the track ID
0234     m_Hit->set_trkid(aTrack->GetTrackID());
0235
0236     // set the initial energy deposit
0237     m_Hit->set_edep(0);
0238
0239     if (whichactive > 0)  // return of IsInIntt, > 0 hit in si-strip, < 0 hit in absorber
0240     {
0241       // Now save the container we want to add this hit to
0242       m_SaveHitContainer = m_HitContainer;
0243     }
0244     else
0245     {
0246       m_SaveHitContainer = m_AbsorberHitContainer;
0247     }
0248
0249     if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0250     {
0251       if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0252       {
0253         m_Hit->set_trkid(pp->GetUserTrackId());
0254         m_Hit->set_shower_id(pp->GetShower()->get_id());
0255         m_SaveShower = pp->GetShower();
0256       }
0257     }
0258     break;
0259
0260   default:
0261     break;
0262   }
0263
0264   // std::cout << "  Update exit values for prePoint->GetStepStatus = " << prePoint->GetStepStatus() << " and postPoint->GetStepStatus = " << postPoint->GetStepStatus() << std::endl;
0265
0266   // here we just update the exit values, it will be overwritten
0267   // for every step until we leave the volume or the particle
0268   // ceases to exist
0269   m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0270   m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0271   m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0272
0273   StoreLocalCoordinate(m_Hit, aStep, false, true);
0274
0275   if (whichactive > 0)
0276   {
0277     m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
0278     m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
0279     m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
0280     m_Hit->set_eion(m_Hit->get_eion() + eion);
0281   }
0282
0283   m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0284
0285   // sum up the energy to get total deposited
0286   m_Hit->set_edep(m_Hit->get_edep() + edep);
0287
0288   if (geantino)
0289   {
0290     m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0291     if (whichactive > 0)
0292     {
0293       m_Hit->set_eion(-1);
0294     }
0295   }
0296
0297   if (edep > 0)
0298   {
0299     if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0300     {
0301       if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0302       {
0303         pp->SetKeep(1);  // we want to keep the track
0304       }
0305     }
0306   }
0307
0308   // if any of these conditions is true this is the last step in
0309   // this volume and we need to save the hit
0310   // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0311   // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0312   // (happens when your detector goes outside world volume)
0313   // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0314   // aTrack->GetTrackStatus() == fStopAndKill is also set)
0315   // aTrack->GetTrackStatus() == fStopAndKill: track ends
0316   if (postPoint->GetStepStatus() == fGeomBoundary ||
0317       postPoint->GetStepStatus() == fWorldBoundary ||
0318       postPoint->GetStepStatus() == fAtRestDoItProc ||
0319       aTrack->GetTrackStatus() == fStopAndKill)
0320   {
0321     if (Verbosity() > 0)
0322     {
0323       std::cout << " postPoint step status changed to " << postPoint->GetStepStatus() << " or aTrack status changed to "
0324                 << aTrack->GetTrackStatus() << std::endl;
0325       std::cout << "  end g4hit for:  aTrack->GetParentID " << aTrack->GetParentID() << " aTrack->GetTrackID " << aTrack->GetTrackID() << " eion = " << eion << std::endl;
0326       std::cout << "     end hit position x,y,z = " << postPoint->GetPosition().x() / cm
0327                 << "    " << postPoint->GetPosition().y() / cm
0328                 << "     " << postPoint->GetPosition().z() / cm
0329                 << std::endl;
0330     }
0331
0332     // save only hits with energy deposit (or -1 for geantino)
0333     if (m_Hit->get_edep())
0334     {
0335       m_SaveHitContainer->AddHit(sphxlayer, m_Hit);
0336       if (m_SaveShower)
0337       {
0338         m_SaveShower->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
0339       }
0340       if (Verbosity() > 0)
0341       {
0342         m_Hit->print();
0343       }
0344       // ownership has been transferred to container, set to null
0345       // so we will create a new hit for the next track
0346       m_Hit = nullptr;
0347     }
0348     else
0349     {
0350       // if this hit has no energy deposit, just reset it for reuse
0351       // this means we have to delete it in the dtor. If this was
0352       // the last hit we processed the memory is still allocated
0353       m_Hit->Reset();
0354     }
0355   }
0356   return true;
0357 }
0358
0359 //____________________________________________________________________________..
0360 void PHG4InttSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0361 {
0362   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0363   m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
0364
0365   // if we do not find the node it's messed up.
0366   if (!m_HitContainer)
0367   {
0368     std::cout << "PHG4InttSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
0369     gSystem->Exit(1);
0370   }
0371
0372   if (!m_AbsorberHitContainer && Verbosity() > 1)
0373   {
0374     std::cout << "PHG4InttSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
0375   }
0376 }
0377
0378 void PHG4InttSteppingAction::SetHitNodeName(const std::string& type, const std::string& name)
0379 {
0380   if (type == "G4HIT")
0381   {
0382     m_HitNodeName = name;
0383     return;
0384   }
0385   else if (type == "G4HIT_ABSORBER")
0386   {
0387     m_AbsorberNodeName = name;
0388     return;
0389   }
0390   std::cout << "Invalid output hit node type " << type << std::endl;
0391   gSystem->Exit(1);
0392   return;
0393 }