g4simulation/g4caloprototype/PHG4SpacalPrototype4SteppingAction.cc

0001 /*!
0002  * \file ${file_name}
0003  * \brief
0004  * \author Jin Huang <jhuang@bnl.gov>
0005  * \version $$Revision: 1.6 $$
0006  * \date $$Date: 2015/01/07 23:50:05 $$
0007  */
0008
0009 #include "PHG4SpacalPrototype4SteppingAction.h"
0010 #include <g4detectors/PHG4CylinderGeom_Spacalv3.h>
0011 #include "PHG4SpacalPrototype4Detector.h"
0012
0013 #include <g4main/PHG4Hit.h>  // for PHG4Hit
0014 #include <g4main/PHG4HitContainer.h>
0015 #include <g4main/PHG4Hitv1.h>
0016 #include <g4main/PHG4Shower.h>
0017 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0018 #include <g4main/PHG4TrackUserInfoV1.h>
0019
0020 #include <phool/getClass.h>
0021
0022 #include <Geant4/G4IonisParamMat.hh>  // for G4IonisParamMat
0023 #include <Geant4/G4Material.hh>       // for G4Material
0024 #include <Geant4/G4MaterialCutsCouple.hh>
0025 #include <Geant4/G4ParticleDefinition.hh>  // for G4ParticleDefinition
0026 #include <Geant4/G4Step.hh>
0027 #include <Geant4/G4StepPoint.hh>   // for G4StepPoint
0028 #include <Geant4/G4StepStatus.hh>  // for fGeomBoundary, fUndefined
0029 #include <Geant4/G4String.hh>      // for G4String
0030 #include <Geant4/G4SystemOfUnits.hh>
0031 #include <Geant4/G4ThreeVector.hh>            // for G4ThreeVector
0032 #include <Geant4/G4TouchableHandle.hh>        // for G4TouchableHandle
0033 #include <Geant4/G4Track.hh>                  // for G4Track
0034 #include <Geant4/G4TrackStatus.hh>            // for fStopAndKill
0035 #include <Geant4/G4Types.hh>                  // for G4double
0036 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0037 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0038
0039 #include <iostream>
0040 #include <string>  // for operator+, operator<<
0041
0042 class G4VPhysicalVolume;
0043 class PHCompositeNode;
0044
0045 using namespace std;
0046 //____________________________________________________________________________..
0047 PHG4SpacalPrototype4SteppingAction::PHG4SpacalPrototype4SteppingAction(PHG4SpacalPrototype4Detector* detector)
0048   : PHG4SteppingAction(detector->GetName())
0049   , detector_(detector)
0050   , hits_(nullptr)
0051   , absorberhits_(nullptr)
0052   , hit(nullptr)
0053   , savehitcontainer(nullptr)
0054   , saveshower(nullptr)
0055 {
0056 }
0057
0058 //____________________________________________________________________________..
0059 bool PHG4SpacalPrototype4SteppingAction::UserSteppingAction(const G4Step* aStep, bool)
0060 {
0061   // get volume of the current step
0062   G4VPhysicalVolume* volume =
0063       aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
0064
0065   // collect energy and track length step by step
0066   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0067   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0068
0069   const G4Track* aTrack = aStep->GetTrack();
0070
0071   const int layer_id = 0;
0072   // make sure we are in a volume
0073   // IsInCylinderActive returns the number of the scintillator
0074   // slat which has fired
0075   int isactive = detector_->IsInCylinderActive(volume);
0076   if (isactive > PHG4SpacalPrototype4Detector::INACTIVE)
0077   {
0078     bool geantino = false;
0079     // the check for the pdg code speeds things up, I do not want to make
0080     // an expensive string compare for every track when we know
0081     // geantino or chargedgeantino has pid=0
0082     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
0083     {
0084       geantino = true;
0085     }
0086     G4StepPoint* prePoint = aStep->GetPreStepPoint();
0087     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0088     int scint_id = -1;
0089
0090     //SPACAL ID that is associated with towers
0091     int sector_ID = 0;
0092     int tower_ID = 0;
0093     int fiber_ID = 0;
0094
0095     if (isactive == PHG4SpacalPrototype4Detector::FIBER_CORE)
0096     {
0097       fiber_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
0098       tower_ID = prePoint->GetTouchable()->GetReplicaNumber(2);
0099       sector_ID = prePoint->GetTouchable()->GetReplicaNumber(3);
0100
0101       if (Verbosity() >= 2)
0102         cout << __PRETTY_FUNCTION__ << " FIBER_CORE step with GetReplicaNumber[0-4] = "
0103              << prePoint->GetTouchable()->GetReplicaNumber(0) << ", "
0104              << prePoint->GetTouchable()->GetReplicaNumber(1) << ", "
0105              << prePoint->GetTouchable()->GetReplicaNumber(2) << ", "
0106              << prePoint->GetTouchable()->GetReplicaNumber(3) << ". edep = " << edep << ", eion = " << eion << endl;
0107     }
0108
0109     else if (isactive == PHG4SpacalPrototype4Detector::FIBER_CLADING)
0110     {
0111       fiber_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
0112       tower_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
0113       sector_ID = prePoint->GetTouchable()->GetReplicaNumber(2);
0114     }
0115
0116     else if (isactive == PHG4SpacalPrototype4Detector::ABSORBER)
0117     {
0118       tower_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
0119       sector_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
0120     }
0121
0122     // compact the tower/sector/fiber ID into 32 bit scint_id, so we could save some space for SPACAL hits
0123     scint_id = PHG4CylinderGeom_Spacalv3::scint_id_coder(sector_ID, tower_ID, fiber_ID).scint_ID;
0124
0125     //       cout << "track id " << aTrack->GetTrackID() << endl;
0126     //        cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
0127     //        cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
0128     switch (prePoint->GetStepStatus())
0129     {
0130     case fGeomBoundary:
0131     case fUndefined:
0132       //    case fPostStepDoItProc: // from point like interaction (compton,decay,hard rad)
0133
0134       if (!hit)
0135       {
0136         hit = new PHG4Hitv1();
0137       }
0138       hit->set_layer((unsigned int) layer_id);
0139       hit->set_scint_id(scint_id);  // isactive contains the scintillator slat id
0140       //here we set the entrance values in cm
0141       hit->set_x(0, prePoint->GetPosition().x() / cm);
0142       hit->set_y(0, prePoint->GetPosition().y() / cm);
0143       hit->set_z(0, prePoint->GetPosition().z() / cm);
0144
0145       // time in ns
0146       hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0147
0148       if (isactive == PHG4SpacalPrototype4Detector::FIBER_CORE)  // only for active areas
0149       {
0150         // store all pre local coordinates
0151         StoreLocalCoordinate(hit, aStep, true, false);
0152       }
0153
0154       //set the track ID
0155       hit->set_trkid(aTrack->GetTrackID());
0156       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0157       {
0158         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0159         {
0160           hit->set_trkid(pp->GetUserTrackId());
0161           hit->set_shower_id(pp->GetShower()->get_id());
0162         }
0163       }
0164       //set the initial energy deposit
0165       hit->set_edep(0);
0166       if (isactive == PHG4SpacalPrototype4Detector::FIBER_CORE)  // only for active areas
0167       {
0168         hit->set_eion(0);  // only for fiber core hits
0169         hit->set_light_yield(0);
0170         savehitcontainer = hits_;
0171       }
0172       else
0173       {
0174         savehitcontainer = absorberhits_;
0175       }
0176
0177       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0178       {
0179         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0180         {
0181           hit->set_trkid(pp->GetUserTrackId());
0182           hit->set_shower_id(pp->GetShower()->get_id());
0183           saveshower = pp->GetShower();
0184         }
0185       }
0186
0187       //        if (hit->get_z(0) > get_zmax() || hit->get_z(0) < get_zmin())
0188       //          {
0189       //            cout << "PHG4SpacalPrototype4SteppingAction: hit outside acceptance, layer: "
0190       //                << layer_id << endl;
0191       //            hit->identify();
0192       //          }
0193       break;
0194     default:
0195       break;
0196     }
0197     // here we just update the exit values, it will be overwritten
0198     // for every step until we leave the volume or the particle
0199     // ceases to exist
0200     hit->set_x(1, postPoint->GetPosition().x() / cm);
0201     hit->set_y(1, postPoint->GetPosition().y() / cm);
0202     hit->set_z(1, postPoint->GetPosition().z() / cm);
0203
0204     hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0205
0206     //sum up the energy to get total deposited
0207     hit->set_edep(hit->get_edep() + edep);
0208
0209     if (isactive == PHG4SpacalPrototype4Detector::FIBER_CORE)  // only for active areas
0210     {
0211       // store all pre local coordinates
0212       StoreLocalCoordinate(hit, aStep, false, true);
0213       hit->set_eion(hit->get_eion() + eion);
0214
0215       double light_yield = GetVisibleEnergyDeposition(aStep);
0216
0217       static bool once = true;
0218       if (once and edep > 0)
0219       {
0220         once = false;
0221
0222         if (Verbosity() > 0)
0223         {
0224           cout << "PHG4SpacalPrototype4SteppingAction::UserSteppingAction::"
0225                //
0226                << detector_->GetName() << " - "
0227                << " use scintillating light model at each Geant4 steps. "
0228                << "First step: "
0229                << "Material = "
0230                << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
0231                << ", "
0232                << "Birk Constant = "
0233                << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetIonisation()->GetBirksConstant()
0234                << ","
0235                << "edep = " << edep << ", "
0236                << "eion = " << eion
0237                << ", "
0238                << "light_yield = " << light_yield << endl;
0239         }
0240       }
0241
0242       hit->set_light_yield(hit->get_light_yield() + light_yield);
0243     }
0244
0245     //      if (hit->get_z(1) > get_zmax() || hit->get_z(1) < get_zmin())
0246     //        {
0247     //          cout << "PHG4SpacalPrototype4SteppingAction: hit outside acceptance get_zmin() "
0248     //              << get_zmin() << ", get_zmax() " << get_zmax() << " at exit"
0249     //              << endl;
0250     //          hit->identify();
0251     //        }
0252     if (geantino)
0253     {
0254       hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0255                           //          hit->set_eion(-1);
0256     }
0257     if (edep > 0)
0258     {
0259       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0260       {
0261         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0262         {
0263           pp->SetKeep(1);  // we want to keep the track
0264         }
0265       }
0266     }
0267
0268     // if any of these conditions is true this is the last step in
0269     // this volume and we need to save the hit
0270     // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0271     // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0272     // (not sure if this will ever be the case)
0273     // aTrack->GetTrackStatus() == fStopAndKill: track ends
0274     if (postPoint->GetStepStatus() == fGeomBoundary ||
0275         postPoint->GetStepStatus() == fWorldBoundary ||
0276         aTrack->GetTrackStatus() == fStopAndKill)
0277     {
0278       // save only hits with energy deposit (or -1 for geantino)
0279       if (hit->get_edep())
0280       {
0281         savehitcontainer->AddHit(layer_id, hit);
0282         if (saveshower)
0283         {
0284           // save shower under container id of active volume (for running without absorber hit container to save space)
0285           saveshower->add_g4hit_id(savehitcontainer->GetID(), hit->get_hit_id());
0286         }
0287         // ownership has been transferred to container, set to null
0288         // so we will create a new hit for the next track
0289         hit = nullptr;
0290       }
0291       else
0292       {
0293         // if this hit has no energy deposit, just reset it for reuse
0294         // this means we have to delete it in the dtor. If this was
0295         // the last hit we processed the memory is still allocated
0296         hit->Reset();
0297       }
0298     }
0299     // return true to indicate the hit was used
0300     return true;
0301   }
0302   else
0303   {
0304     return false;
0305   }
0306 }
0307
0308 //____________________________________________________________________________..
0309 void PHG4SpacalPrototype4SteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0310 {
0311   string hitnodename;
0312   string absorbernodename;
0313   if (detector_->SuperDetector() != "NONE")
0314   {
0315     hitnodename = "G4HIT_" + detector_->SuperDetector();
0316     absorbernodename = "G4HIT_ABSORBER_" + detector_->SuperDetector();
0317   }
0318   else
0319   {
0320     hitnodename = "G4HIT_" + detector_->GetName();
0321     absorbernodename = "G4HIT_ABSORBER_" + detector_->GetName();
0322   }
0323
0324   //now look for the map and grab a pointer to it.
0325   hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
0326   absorberhits_ = findNode::getClass<PHG4HitContainer>(topNode,
0327                                                        absorbernodename.c_str());
0328   // if we do not find the node we need to make it.
0329   if (!hits_)
0330   {
0331     std::cout << "PHG4SpacalPrototype4SteppingAction::SetTopNode - unable to find "
0332               << hitnodename << std::endl;
0333   }
0334   if (!absorberhits_)
0335   {
0336     if (Verbosity() > 1)
0337     {
0338       std::cout << "PHG4SpacalPrototype4SteppingAction::SetTopNode - unable to find "
0339                 << absorbernodename << std::endl;
0340     }
0341   }
0342 }
0343
0344 double
0345 PHG4SpacalPrototype4SteppingAction::get_zmin()
0346 {
0347   if (!detector_)
0348     return 0;
0349   else
0350     return detector_->get_geom()->get_zmin() - .0001;
0351 }
0352
0353 double
0354 PHG4SpacalPrototype4SteppingAction::get_zmax()
0355 {
0356   if (!detector_)
0357     return 0;
0358   else
0359     return detector_->get_geom()->get_zmax() + .0001;
0360 }