g4simulation/g4mvtx/PHG4EICMvtxSteppingAction.cc

0001 #include "PHG4EICMvtxSteppingAction.h"
0002
0003 #include "PHG4EICMvtxDetector.h"
0004
0005 #include <g4main/PHG4Hit.h>
0006 #include <g4main/PHG4HitContainer.h>
0007 #include <g4main/PHG4Hitv1.h>
0008 #include <g4main/PHG4Shower.h>
0009 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0010 #include <g4main/PHG4TrackUserInfoV1.h>
0011
0012 #include <phool/getClass.h>
0013 #include <phool/phool.h>  // for PHWHERE
0014
0015 #include <Geant4/G4NavigationHistory.hh>
0016 #include <Geant4/G4ParticleDefinition.hh>      // for G4ParticleDefinition
0017 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHandle
0018 #include <Geant4/G4Step.hh>
0019 #include <Geant4/G4StepPoint.hh>   // for G4StepPoint
0020 #include <Geant4/G4StepStatus.hh>  // for fGeomBoundary, fAtRest...
0021 #include <Geant4/G4String.hh>      // for G4String
0022 #include <Geant4/G4SystemOfUnits.hh>
0023 #include <Geant4/G4ThreeVector.hh>            // for G4ThreeVector
0024 #include <Geant4/G4TouchableHandle.hh>        // for G4TouchableHandle
0025 #include <Geant4/G4Track.hh>                  // for G4Track
0026 #include <Geant4/G4TrackStatus.hh>            // for fStopAndKill
0027 #include <Geant4/G4Types.hh>                  // for G4double
0028 #include <Geant4/G4VPhysicalVolume.hh>        // for G4VPhysicalVolume
0029 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0030 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0031
0032 #include <boost/tokenizer.hpp>
0033 // this is an ugly hack, the gcc optimizer has a bug which
0034 // triggers the uninitialized variable warning which
0035 // stops compilation because of our -Werror
0036 #include <boost/version.hpp>  // to get BOOST_VERSION
0037 #if (__GNUC__ == 4 && __GNUC_MINOR__ == 4 && BOOST_VERSION == 105700)
0038 #pragma GCC diagnostic ignored "-Wuninitialized"
0039 #pragma message "ignoring bogus gcc warning in boost header lexical_cast.hpp"
0040 #include <boost/lexical_cast.hpp>
0041 #pragma GCC diagnostic warning "-Wuninitialized"
0042 #else
0043 #include <boost/lexical_cast.hpp>
0044 #endif
0045
0046 #include <cstdlib>  // for exit
0047 #include <iostream>
0048 #include <string>  // for operator<<, basic_string
0049
0050 class PHCompositeNode;
0051
0052 using namespace std;
0053 //____________________________________________________________________________..
0054 PHG4EICMvtxSteppingAction::PHG4EICMvtxSteppingAction(PHG4EICMvtxDetector* detector)
0055   : PHG4SteppingAction(detector->GetName())
0056   , m_Detector(detector)
0057   , m_HitContainer(nullptr)
0058   , m_AbsorberhitContainer(nullptr)
0059   , m_Hit(nullptr)
0060   , m_SaveShower(nullptr)
0061 {
0062 }
0063
0064 //____________________________________________________________________________..
0065 PHG4EICMvtxSteppingAction::~PHG4EICMvtxSteppingAction()
0066 {
0067   // if the last hit was a zero energie deposit hit, it is just reset
0068   // and the memory is still allocated, so we need to delete it here
0069   // if the last hit was saved, hit is a nullptr pointer which are
0070   // legal to delete (it results in a no operation)
0071   delete m_Hit;
0072 }
0073
0074 //____________________________________________________________________________..
0075 bool PHG4EICMvtxSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0076 {
0077   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0078   // get volume of the current step
0079   G4VPhysicalVolume* sensor_volume = touch->GetVolume();
0080
0081   // PHG4MvtxDetector->IsInMvtx(volume)
0082   // returns
0083   //  0 if outside of Mvtx
0084   //  1 if inside sensor
0085
0086   // This checks if the volume is a sensor (doesn't tell us unique layer)
0087   // PHG4MvtxDetector->IsSensor(volume)
0088   // returns
0089   //  1 if volume is a sensor
0090   //  0 if not
0091   int whichactive = m_Detector->IsSensor(sensor_volume);
0092
0093   if (!whichactive)
0094   {
0095     return false;
0096   }
0097
0098   // This tells us if the volume belongs to the right stave for this layer
0099   // From the GDML file the 3rd volume up should be the half-stave
0100   // PHG4MvtxTelescopeDetector_->IsInMvtx(volume)
0101   // returns
0102   //  1 if in ladder belonging to this layer
0103   //  0 if not
0104   int layer_id = -9999;
0105   int stave_id = -9999;
0106   // cout << endl << "  In UserSteppingAction for layer " << layer_id << endl;
0107   G4VPhysicalVolume* vstave = touch->GetVolume(3);
0108   whichactive = m_Detector->IsInMvtx(vstave, layer_id, stave_id);
0109   if (layer_id < 0 || stave_id < 0)
0110   {
0111     cout << PHWHERE << "invalid Mvtx's layer (" << layer_id << ") or stave (" << stave_id << ") index " << endl;
0112     exit(1);
0113   }
0114
0115   if (!whichactive)
0116   {
0117     return false;
0118   }
0119
0120   if (Verbosity() > 5)
0121   {
0122     // make sure we know where we are!
0123     G4VPhysicalVolume* vtest = touch->GetVolume();
0124     cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume " << vtest->GetName() << endl;
0125     G4VPhysicalVolume* vtest1 = touch->GetVolume(1);
0126     cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume 1 up " << vtest1->GetName() << endl;
0127     G4VPhysicalVolume* vtest2 = touch->GetVolume(2);
0128     cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume 2 up " << vtest2->GetName() << endl;
0129     G4VPhysicalVolume* vtest3 = touch->GetVolume(3);
0130     cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume 3 up " << vtest3->GetName() << endl;
0131     G4VPhysicalVolume* vtest4 = touch->GetVolume(4);
0132     cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume 4 up " << vtest4->GetName() << endl;
0133   }
0134
0135   //=======================================================================
0136   // We want the location of the hit
0137   // Here we will record in the hit object:
0138   //   The stave, half-stave, module and chip numbers
0139   //   The energy deposited
0140   //   The entry point and exit point in world coordinates
0141   //   The entry point and exit point in local (sensor) coordinates
0142   // The pixel number will be derived later from the entry and exit points in the sensor local coordinates
0143   //=======================================================================
0144
0145   // int stave_number = -1; // stave_number from stave map values
0146   int half_stave_number = -1;
0147   int module_number = -1;
0148   int chip_number = -1;
0149
0150   if (Verbosity() > 0) {
0151     cout << endl
0152          << "  UserSteppingAction: layer " << layer_id;
0153 }
0154   boost::char_separator<char> sep("_");
0155   boost::tokenizer<boost::char_separator<char> >::const_iterator tokeniter;
0156
0157   // OLD ITS.gdml: chip number is from  "ITSUChip[layer number]_[chip number]
0158   // NEW: chip number is from  "MVTXChip_[chip number]
0159   G4VPhysicalVolume* v1 = touch->GetVolume(1);
0160   boost::tokenizer<boost::char_separator<char> > tok1(v1->GetName(), sep);
0161   tokeniter = tok1.begin();
0162   ++tokeniter;
0163   chip_number = boost::lexical_cast<int>(*tokeniter);
0164   if (Verbosity() > 0) { cout << " chip  " << chip_number;
0165 }
0166   G4VPhysicalVolume* v2 = touch->GetVolume(2);
0167   boost::tokenizer<boost::char_separator<char> > tok2(v2->GetName(), sep);
0168   tokeniter = tok2.begin();
0169   ++tokeniter;
0170   module_number = boost::lexical_cast<int>(*tokeniter);
0171   if (Verbosity() > 0) { cout << " module " << module_number;
0172 }
0173
0174   // The stave number  is the imprint number from the assembly volume imprint
0175   // The assembly volume history string format is (e.g.):
0176   // av_13_impr_4_ITSUHalfStave6_pv_1
0177   // where "impr_4" says stave number 4
0178   // and where "ITSUHalfStave6_pv_1" says hald stave number 1 in layer number 6
0179
0180   G4VPhysicalVolume* v3 = touch->GetVolume(3);
0181   boost::tokenizer<boost::char_separator<char> > tok3(v3->GetName(), sep);
0182   tokeniter = tok3.begin();
0183   ++tokeniter;
0184   ++tokeniter;
0185   ++tokeniter;
0186   // stave_number = boost::lexical_cast<int>(*tokeniter) - 1;  // starts counting imprints at 1, we count staves from 0!
0187   if (Verbosity() > 0) { cout << " stave " << stave_id;
0188 }
0189   ++tokeniter;
0190   ++tokeniter;
0191   ++tokeniter;
0192   half_stave_number = boost::lexical_cast<int>(*tokeniter);
0193   if (Verbosity() > 0) { cout << " half_stave " << half_stave_number;
0194 }
0195
0196   // FYI: doing string compares inside a stepping action sounds like a recipe
0197   // for failure inside a heavy ion event... we'll wait and see how badly
0198   // this profiles. -MPM
0199
0200   // Now we want to collect information about the hit
0201
0202   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0203   const G4Track* aTrack = aStep->GetTrack();
0204   if (Verbosity() > 0) { cout << " edep = " << edep << endl;
0205 }
0206
0207   // if this cylinder stops everything, just put all kinetic energy into edep
0208   if (m_Detector->IsBlackHole(layer_id))
0209   {
0210     edep = aTrack->GetKineticEnergy() / GeV;
0211     G4Track* killtrack = const_cast<G4Track*>(aTrack);
0212     killtrack->SetTrackStatus(fStopAndKill);
0213   }
0214
0215   // test if we are active
0216   if (m_Detector->IsActive(layer_id))
0217   {
0218     bool geantino = false;
0219     // the check for the pdg code speeds things up, I do not want to make
0220     // an expensive string compare for every track when we know
0221     // geantino or chargedgeantino has pid=0
0222     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0223         aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
0224     {
0225       geantino = true;
0226     }
0227
0228     G4ThreeVector worldPosition;  // localPosition;
0229     G4TouchableHandle theTouchable;
0230
0231     G4StepPoint* prePoint = aStep->GetPreStepPoint();
0232     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0233
0234     if (Verbosity() > 0)
0235     {
0236       G4ParticleDefinition* def = aTrack->GetDefinition();
0237       cout << " Particle: " << def->GetParticleName() << endl;
0238     }
0239
0240     switch (prePoint->GetStepStatus())
0241     {
0242     case fGeomBoundary:
0243     case fUndefined:
0244       if (!m_Hit)
0245       {
0246         m_Hit = new PHG4Hitv1();
0247       }
0248       m_Hit->set_layer((unsigned int) layer_id);
0249
0250       // set the index values needed to locate the sensor strip
0251       m_Hit->set_property(PHG4Hit::prop_stave_index, stave_id);
0252       m_Hit->set_property(PHG4Hit::prop_half_stave_index, half_stave_number);
0253       m_Hit->set_property(PHG4Hit::prop_module_index, module_number);
0254       m_Hit->set_property(PHG4Hit::prop_chip_index, chip_number);
0255
0256       worldPosition = prePoint->GetPosition();
0257
0258       if (Verbosity() > 0)
0259       {
0260         theTouchable = prePoint->GetTouchableHandle();
0261         cout << "entering: depth = " << theTouchable->GetHistory()->GetDepth() << endl;
0262         G4VPhysicalVolume* vol1 = theTouchable->GetVolume();
0263         cout << "entering volume name = " << vol1->GetName() << endl;
0264       }
0265
0266       /*
0267           localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
0268           */
0269
0270       // Store the local coordinates for the entry point
0271       StoreLocalCoordinate(m_Hit, aStep, true, false);
0272
0273       // Store the entrance values in cm in world coordinates
0274       m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0275       m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0276       m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0277
0278       m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
0279       m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
0280       m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
0281
0282       // time in ns
0283       m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0284       // set the track ID
0285       m_Hit->set_trkid(aTrack->GetTrackID());
0286       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0287       {
0288         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0289         {
0290           m_Hit->set_trkid(pp->GetUserTrackId());
0291           m_Hit->set_shower_id(pp->GetShower()->get_id());
0292           m_SaveShower = pp->GetShower();
0293         }
0294       }
0295       // set the initial energy deposit
0296       m_Hit->set_edep(0);
0297
0298       // Now add the hit
0299       //      m_Hit->print();
0300       // m_HitContainer->AddHit(layer_id, m_Hit);
0301       break;
0302     default:
0303       break;
0304     }
0305
0306     // here we just update the exit values, it will be overwritten
0307     // for every step until we leave the volume or the particle
0308     // ceases to exist
0309
0310     /*
0311       // Note that the after you reach the boundary the touchable for the postPoint points to the next volume, not the sensor!
0312       // This was given to me as the way to get back to the sensor volume, but it does not work
0313       theTouchable = postPoint->GetTouchableHandle();
0314       localPosition = history->GetTransform(history->GetDepth() - 1).TransformPoint(worldPosition);
0315       cout << "Exit local coords: x " <<  localPosition.x() / cm << " y " <<  localPosition.y() / cm << " z " <<  localPosition.z() / cm << endl;
0316       */
0317
0318     /*
0319       // Use the prePoint from the final step  for now, until I understand how to get the exit point in the sensor volume coordinates
0320       //======================================================================================
0321       theTouchable = prePoint->GetTouchableHandle();
0322       vol2 = theTouchable->GetVolume();
0323       if(Verbosity() > 0)
0324         cout << "exiting volume name = " << vol2->GetName() << endl;
0325       worldPosition = prePoint->GetPosition();
0326       if(Verbosity() > 0)
0327         cout << "Exit world coords prePoint: x " <<  worldPosition.x() / cm << " y " <<  worldPosition.y() / cm << " z " <<  worldPosition.z() / cm << endl;
0328
0329       // This is for consistency with the local coord position, the world coordinate exit position is correct
0330       m_Hit->set_x( 1, prePoint->GetPosition().x() / cm );
0331       m_Hit->set_y( 1, prePoint->GetPosition().y() / cm );
0332       m_Hit->set_z( 1, prePoint->GetPosition().z() / cm );
0333
0334       const G4NavigationHistory *history = theTouchable->GetHistory();
0335       //cout << "exiting: depth = " << history->GetDepth() <<  " volume name = " << history->GetVolume(history->GetDepth())->GetName() << endl;
0336       localPosition = history->GetTransform(history->GetDepth()).TransformPoint(worldPosition);
0337
0338       m_Hit->set_local_x(1, localPosition.x() / cm);
0339       m_Hit->set_local_y(1, localPosition.y() / cm);
0340       m_Hit->set_local_z(1, localPosition.z() / cm);
0341       */
0342
0343     // Store the local coordinates for the exit point
0344     StoreLocalCoordinate(m_Hit, aStep, false, true);
0345
0346     // Store world coordinates for the exit point
0347     m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0348     m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0349     m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0350
0351     m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
0352     m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
0353     m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
0354
0355     m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0356     // sum up the energy to get total deposited
0357     m_Hit->set_edep(m_Hit->get_edep() + edep);
0358     if (geantino)
0359     {
0360       m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0361     }
0362     if (edep > 0)
0363     {
0364       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0365       {
0366         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0367         {
0368           pp->SetKeep(1);  // we want to keep the track
0369         }
0370       }
0371     }
0372
0373     if (Verbosity() > 0)
0374     {
0375       G4StepPoint* prePointA = aStep->GetPreStepPoint();
0376       G4StepPoint* postPointA = aStep->GetPostStepPoint();
0377       cout << "----- PHg4MvtxSteppingAction::UserSteppingAction - active volume = " << sensor_volume->GetName() << endl;
0378       cout << "       layer = " << layer_id << endl;
0379       cout << "       stave number = " << stave_id << " half_stave_number = " << half_stave_number << endl;
0380       cout << "       module number  = " << module_number << endl;
0381       cout << "       chip number = " << chip_number << endl;
0382       cout << "       prepoint x position " << prePointA->GetPosition().x() / cm << endl;
0383       cout << "       prepoint y position " << prePointA->GetPosition().y() / cm << endl;
0384       cout << "       prepoint z position " << prePointA->GetPosition().z() / cm << endl;
0385       cout << "       postpoint x position " << postPointA->GetPosition().x() / cm << endl;
0386       cout << "       postpoint y position " << postPointA->GetPosition().y() / cm << endl;
0387       cout << "       postpoint z position " << postPointA->GetPosition().z() / cm << endl;
0388       cout << "       edep " << edep << endl;
0389     }
0390
0391     if (Verbosity() > 0)
0392     {
0393       cout << "  stepping action found hit:" << endl;
0394       m_Hit->print();
0395       cout << endl
0396            << endl;
0397     }
0398
0399     // if any of these conditions is true this is the last step in
0400     // this volume and we need to save the hit
0401     // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0402     // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0403     // (happens when your detector goes outside world volume)
0404     // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0405     // aTrack->GetTrackStatus() == fStopAndKill is also set)
0406     // aTrack->GetTrackStatus() == fStopAndKill: track ends
0407     if (postPoint->GetStepStatus() == fGeomBoundary ||
0408         postPoint->GetStepStatus() == fWorldBoundary ||
0409         postPoint->GetStepStatus() == fAtRestDoItProc ||
0410         aTrack->GetTrackStatus() == fStopAndKill)
0411     {
0412       // save only hits with energy deposit (or -1 for geantino)
0413       if (m_Hit->get_edep())
0414       {
0415         m_HitContainer->AddHit(layer_id, m_Hit);
0416         if (m_SaveShower)
0417         {
0418           m_SaveShower->add_g4hit_id(m_HitContainer->GetID(), m_Hit->get_hit_id());
0419         }
0420         // ownership has been transferred to container, set to null
0421         // so we will create a new hit for the next track
0422         m_Hit = nullptr;
0423       }
0424       else
0425       {
0426         // if this hit has no energy deposit, just reset it for reuse
0427         // this means we have to delete it in the dtor. If this was
0428         // the last hit we processed the memory is still allocated
0429         m_Hit->Reset();
0430       }
0431     }
0432
0433     // return true to indicate the hit was used
0434     return true;
0435   }
0436   else
0437   {
0438     return false;
0439   }
0440
0441   return false;
0442 }
0443
0444 //____________________________________________________________________________..
0445 void PHG4EICMvtxSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0446 {
0447   string hitnodename;
0448   string absorbernodename;
0449   if (m_Detector->SuperDetector() != "NONE")
0450   {
0451     hitnodename = "G4HIT_" + m_Detector->SuperDetector();
0452     absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
0453   }
0454   else
0455   {
0456     hitnodename = "G4HIT_" + m_Detector->GetName();
0457     absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
0458   }
0459
0460   // now look for the map and grab a pointer to it.
0461   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
0462   m_AbsorberhitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
0463
0464   // if we do not find the node it's messed up.
0465   if (!m_HitContainer)
0466   {
0467     std::cout << "PHG4MvtxSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
0468   }
0469   if (!m_AbsorberhitContainer)
0470   {
0471     if (Verbosity() > 0)
0472     {
0473       cout << "PHG4MvtxSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
0474     }
0475   }
0476 }