g4simulation/g4caloprototype/PHG4HcalPrototypeSteppingAction.cc

0001 // This is the steppingaction code for the hcal prototype
0002 // created on 1/27/2014, Liang, HeXC
0003 //
0004 #include "PHG4HcalPrototypeSteppingAction.h"
0005 #include "PHG4HcalPrototypeDetector.h"
0006
0007 #include <g4main/PHG4Hit.h>
0008 #include <g4main/PHG4HitContainer.h>
0009 #include <g4main/PHG4Hitv1.h>
0010 #include <g4main/PHG4Shower.h>
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, fUndefined
0020 #include <Geant4/G4String.hh>                 // for G4String
0021 #include <Geant4/G4SystemOfUnits.hh>          // for cm, GeV, nanosecond
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 G4int, 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 operator+, string, ope...
0033
0034 class PHCompositeNode;
0035
0036 using namespace std;
0037 //____________________________________________________________________________..
0038 PHG4HcalPrototypeSteppingAction::PHG4HcalPrototypeSteppingAction(PHG4HcalPrototypeDetector* detector)
0039   : PHG4SteppingAction(detector->GetName())
0040   , detector_(detector)
0041   , hits_(nullptr)
0042   , absorberhits_(nullptr)
0043   , hit(nullptr)
0044 {
0045 }
0046
0047 //____________________________________________________________________________..
0048 bool PHG4HcalPrototypeSteppingAction::UserSteppingAction(const G4Step* aStep, bool)
0049 {
0050   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0051   // get volume of the current step
0052   G4VPhysicalVolume* volume = touch->GetVolume();
0053
0054   // We simply use ourown scintID to test on this condition
0055   //  int whichactive = detector_->IsInHcalPrototype(volume);
0056
0057   /*
0058     if (!whichactive)
0059     {
0060     return false;
0061     }
0062   */
0063
0064   // collect energy and track length step by step
0065   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0066   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0067
0068   const G4Track* aTrack = aStep->GetTrack();
0069
0070   // if this block stops everything, just put all kinetic energy into edep
0071   if (detector_->IsBlackHole())
0072   {
0073     edep = aTrack->GetKineticEnergy() / GeV;
0074     G4Track* killtrack = const_cast<G4Track*>(aTrack);
0075     killtrack->SetTrackStatus(fStopAndKill);
0076   }
0077
0078   // Add detector element ID into to the hit
0079
0080   G4int scintSheetCopyNumber, scintSheetLayerNumber;
0081
0082   G4int sectionID = 0;
0083   G4int scintID = 0x0;
0084   //  G4int absID = 0;
0085
0086   if (volume->GetName() == "hcal1AbsLayer")
0087   {
0088     sectionID = 1;
0089     //      absID = volume->GetCopyNo()+100;   // Inner section absorber layer ID's are 100, 101, ..., 115
0090     scintID = -1;  // scintID is not valid, we set it to -1
0091   }
0092   else if (volume->GetName() == "hcal2AbsLayer")
0093   {
0094     sectionID = 2;
0095     //      absID = volume->GetCopyNo()+200;   // Outer section absorber layer ID's are 200, 201, ..., 215
0096     scintID = -1;  // scintID is not valid, we set it to -1
0097   }
0098   else if (volume->GetName() == "inner1USheet")
0099   {
0100     sectionID = 1;
0101     scintSheetCopyNumber = volume->GetCopyNo();  // either 1 or 0
0102     scintSheetLayerNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(1)->GetCopyNo();
0103     scintID = (scintSheetLayerNumber << 2) + (1 << 1) + scintSheetCopyNumber;
0104     //      absID = -1;                        // absID is not valid, we set it to -1
0105     //G4cout << " **************** scintID: " << scintID << G4endl;
0106   }
0107   else if (volume->GetName() == "inner2USheet")
0108   {
0109     sectionID = 1;
0110     scintSheetCopyNumber = volume->GetCopyNo();  // either 1 or 0
0111     scintSheetLayerNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(1)->GetCopyNo();
0112     scintID = (scintSheetLayerNumber << 2) + (0 << 1) + scintSheetCopyNumber;
0113     //      absID = -1;                        // absID is not valid, we set it to -1
0114     //G4cout << " **************** scintID: " << scintID << G4endl;
0115   }
0116   else if (volume->GetName() == "outer1USheet")
0117   {
0118     sectionID = 2;
0119     scintSheetCopyNumber = volume->GetCopyNo();  // either 1 or 0
0120     scintSheetLayerNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(1)->GetCopyNo();
0121     scintID = (scintSheetLayerNumber << 2) + (1 << 1) + scintSheetCopyNumber;
0122     //      absID = -1;                        // absID is not valid, we set it to -1
0123     //G4cout << " **************** scintID: " << scintID << G4endl;
0124   }
0125   else if (volume->GetName() == "outer2USheet")
0126   {
0127     sectionID = 2;
0128     scintSheetCopyNumber = volume->GetCopyNo();  // either 1 or 0
0129     scintSheetLayerNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(1)->GetCopyNo();
0130     scintID = (scintSheetLayerNumber << 2) + (0 << 1) + scintSheetCopyNumber;
0131     //      absID = -1;                        // absID is not valid, we set it to -1
0132     //G4cout << " **************** scintID: " << scintID << G4endl;
0133   }
0134   else
0135   {
0136     return false;
0137   }
0138
0139   // make sure we are in a volume
0140   if (detector_->IsActive())
0141   {
0142     bool geantino = false;
0143     // the check for the pdg code speeds things up, I do not want to make
0144     // an expensive string compare for every track when we know
0145     // geantino or chargedgeantino has pid=0
0146     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0147         aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
0148     {
0149       geantino = true;
0150     }
0151
0152     G4StepPoint* prePoint = aStep->GetPreStepPoint();
0153     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0154     //       cout << "track id " << aTrack->GetTrackID() << endl;
0155     //       cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
0156     //       cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
0157     switch (prePoint->GetStepStatus())
0158     {
0159     case fGeomBoundary:
0160     case fUndefined:
0161       hit = new PHG4Hitv1();
0162       hit->set_layer((unsigned int) sectionID);
0163       hit->set_scint_id(scintID);
0164       //here we set the entrance values in cm
0165       hit->set_x(0, prePoint->GetPosition().x() / cm);
0166       hit->set_y(0, prePoint->GetPosition().y() / cm);
0167       hit->set_z(0, prePoint->GetPosition().z() / cm);
0168       // time in ns
0169       hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0170       //set the track ID
0171       {
0172         hit->set_trkid(aTrack->GetTrackID());
0173         if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0174         {
0175           if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0176           {
0177             hit->set_trkid(pp->GetUserTrackId());
0178             hit->set_shower_id(pp->GetShower()->get_id());
0179           }
0180         }
0181       }
0182
0183       //set the initial energy deposit
0184       hit->set_edep(0);
0185       hit->set_eion(0);
0186       if (scintID > 0)  // return of isinHcalTestDetector, > 0 hit in scintillator, < 0 hit in absorber
0187       {
0188         // Now add the hit
0189         hits_->AddHit(sectionID, hit);
0190
0191         {
0192           if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0193           {
0194             if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0195             {
0196               pp->GetShower()->add_g4hit_id(hits_->GetID(), hit->get_hit_id());
0197             }
0198           }
0199         }
0200       }
0201       else
0202       {
0203         absorberhits_->AddHit(sectionID, hit);
0204
0205         {
0206           if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0207           {
0208             if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0209             {
0210               pp->GetShower()->add_g4hit_id(absorberhits_->GetID(), hit->get_hit_id());
0211             }
0212           }
0213         }
0214       }
0215       break;
0216     default:
0217       break;
0218     }
0219
0220     // here we just update the exit values, it will be overwritten
0221     // for every step until we leave the volume or the particle
0222     // ceases to exist
0223     hit->set_x(1, postPoint->GetPosition().x() / cm);
0224     hit->set_y(1, postPoint->GetPosition().y() / cm);
0225     hit->set_z(1, postPoint->GetPosition().z() / cm);
0226
0227     hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0228     //sum up the energy to get total deposited
0229     hit->set_edep(hit->get_edep() + edep);
0230     hit->set_eion(hit->get_eion() + eion);
0231     if (geantino)
0232     {
0233       hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0234       hit->set_eion(-1);
0235     }
0236     if (edep > 0)
0237     {
0238       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0239       {
0240         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0241         {
0242           pp->SetKeep(1);  // we want to keep the track
0243         }
0244       }
0245     }
0246
0247     //       hit->identify();
0248     // return true to indicate the hit was used
0249     return true;
0250   }
0251   else
0252   {
0253     return false;
0254   }
0255 }
0256
0257 //____________________________________________________________________________..
0258 void PHG4HcalPrototypeSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0259 {
0260   string hitnodename;
0261   string absorbernodename;
0262   if (detector_->SuperDetector() != "NONE")
0263   {
0264     hitnodename = "G4HIT_" + detector_->SuperDetector();
0265     absorbernodename = "G4HIT_ABSORBER_" + detector_->SuperDetector();
0266   }
0267   else
0268   {
0269     hitnodename = "G4HIT_" + detector_->GetName();
0270     absorbernodename = "G4HIT_ABSORBER_" + detector_->GetName();
0271   }
0272
0273   //now look for the map and grab a pointer to it.
0274   hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
0275   absorberhits_ = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
0276
0277   // if we do not find the node it's messed up.
0278   if (!hits_)
0279   {
0280     std::cout << "PHG4HcalPrototypeSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
0281   }
0282   if (!absorberhits_)
0283   {
0284     if (Verbosity() > 0)
0285     {
0286       cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
0287     }
0288   }
0289 }