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