g4simulation/g4detectors/PHG4SpacalSteppingAction.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 "PHG4SpacalSteppingAction.h"
0010 #include "PHG4CylinderGeom_Spacalv3.h"
0011 #include "PHG4SpacalDetector.h"
0012
0013 #include "PHG4CellDefs.h"
0014 #include "PHG4CylinderCellGeom.h"
0015 #include "PHG4CylinderCellGeomContainer.h"
0016 #include "PHG4CylinderCellGeom_Spacalv1.h"
0017 #include "PHG4CylinderGeom.h"  // for PHG4CylinderGeom
0018 #include "PHG4CylinderGeomContainer.h"
0019 #include "PHG4CylinderGeom_Spacalv1.h"  // for PHG4CylinderGeom_Spaca...
0020
0021 #include <calobase/TowerInfo.h>
0022 #include <calobase/TowerInfoContainer.h>
0023 #include <calobase/TowerInfoContainerv1.h>
0024 #include <calobase/TowerInfoDefs.h>
0025
0026 #include <g4main/PHG4Hit.h>  // for PHG4Hit
0027 #include <g4main/PHG4HitContainer.h>
0028 #include <g4main/PHG4Hitv1.h>
0029 #include <g4main/PHG4Shower.h>
0030 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0031 #include <g4main/PHG4TrackUserInfoV1.h>
0032
0033 #include <phool/PHCompositeNode.h>
0034 #include <phool/PHIODataNode.h>
0035 #include <phool/PHNode.h>  // for PHNode
0036 #include <phool/PHNodeIterator.h>
0037 #include <phool/PHObject.h>  // for PHObject
0038 #include <phool/getClass.h>
0039 #include <phool/phool.h>  // for PHWHERE
0040 #include <phool/recoConsts.h>
0041
0042 #include <phparameter/PHParameters.h>
0043
0044 #include <Geant4/G4IonisParamMat.hh>  // for G4IonisParamMat
0045 #include <Geant4/G4Material.hh>       // for G4Material
0046 #include <Geant4/G4MaterialCutsCouple.hh>
0047 #include <Geant4/G4ParticleDefinition.hh>      // for G4ParticleDefinition
0048 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHandle
0049 #include <Geant4/G4Step.hh>
0050 #include <Geant4/G4StepPoint.hh>   // for G4StepPoint
0051 #include <Geant4/G4StepStatus.hh>  // for fGeomBoundary, fAtRestD...
0052 #include <Geant4/G4String.hh>      // for G4String
0053 #include <Geant4/G4SystemOfUnits.hh>
0054 #include <Geant4/G4ThreeVector.hh>      // for G4ThreeVector
0055 #include <Geant4/G4TouchableHandle.hh>  // for G4TouchableHandle
0056 #include <Geant4/G4Track.hh>            // for G4Track
0057 #include <Geant4/G4TrackStatus.hh>      // for fStopAndKill
0058 #include <Geant4/G4TransportationManager.hh>
0059 #include <Geant4/G4Types.hh>                  // for G4double
0060 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0061 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0062
0063 #include <TSystem.h>
0064
0065 #include <cmath>    // for isfinite
0066 #include <cstdlib>  // for exit
0067 #include <iostream>
0068 #include <string>  // for operator<<, char_traits
0069
0070 class G4VPhysicalVolume;
0071 class PHCompositeNode;
0072
0073 //____________________________________________________________________________..
0074 PHG4SpacalSteppingAction::PHG4SpacalSteppingAction(PHG4SpacalDetector *indetector, const PHParameters *parameters)
0075   : PHG4SteppingAction(indetector->GetName())
0076   , m_Detector(indetector)
0077   , m_Params(parameters)
0078   , m_doG4Hit(m_Params->get_int_param("saveg4hit"))
0079   , m_tmin(m_Params->get_double_param("tmin"))
0080   , m_tmax(m_Params->get_double_param("tmax"))
0081   , m_dt(m_Params->get_double_param("dt"))
0082 {
0083 }
0084
0085 PHG4SpacalSteppingAction::~PHG4SpacalSteppingAction()
0086 {
0087   // if the last hit was a zero energie deposit hit, it is just reset
0088   // and the memory is still allocated, so we need to delete it here
0089   // if the last hit was saved, hit is a nullptr pointer which are
0090   // legal to delete (it results in a no operation)
0091   delete m_Hit;
0092 }
0093
0094 int PHG4SpacalSteppingAction::InitWithNode(PHCompositeNode *topNode)
0095 {
0096   if (m_doG4Hit)
0097   {
0098     return 0;
0099   }
0100   PHNodeIterator iter(topNode);
0101   detector = m_Detector->SuperDetector();
0102   // Looking for the DST node
0103   PHCompositeNode *dstNode;
0104   dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0105   if (!dstNode)
0106   {
0107     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0108     exit(1);
0109   }
0110
0111   // add towerinfo here
0112   PHNodeIterator dstiter(dstNode);
0113   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", detector));
0114   if (!DetNode)
0115   {
0116     DetNode = new PHCompositeNode(detector);
0117     dstNode->addNode(DetNode);
0118   }
0119   m_CaloInfoContainer = new TowerInfoContainerv1(TowerInfoContainer::DETECTOR::EMCAL);
0120   PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(m_CaloInfoContainer, "TOWERINFO_SIM_" + detector, "PHObject");
0121   DetNode->addNode(towerNode);
0122
0123   return 0;
0124 }
0125
0126 int PHG4SpacalSteppingAction::SetUpGeomNode(PHCompositeNode *topNode)
0127 {
0128   if (m_geomsetup)
0129   {
0130     return 0;
0131   }
0132
0133   if (m_doG4Hit)
0134   {
0135     return 0;
0136   }
0137   PHNodeIterator iter(topNode);
0138   detector = m_Detector->SuperDetector();
0139
0140   geonodename = "CYLINDERGEOM_" + detector;
0141   seggeonodename = "CYLINDERCELLGEOM_" + detector;
0142   // get the private layergeo
0143   _layergeo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename);
0144   if (!_layergeo)
0145   {
0146     std::cout << "PHG4SpacalSteppingAction::InitWithNode - Fatal Error - Could not locate sim geometry node "
0147               << geonodename << std::endl;
0148     exit(1);
0149   }
0150
0151   _seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename);
0152   if (!_seggeo)
0153   {
0154     std::cout << "PHG4FullProjSpacalCellReco::process_event - Fatal Error - could not locate cell geometry node "
0155               << seggeonodename << std::endl;
0156     exit(1);
0157   }
0158   PHG4CylinderCellGeom *_geo_raw = _seggeo->GetFirstLayerCellGeom();
0159   _geo = dynamic_cast<PHG4CylinderCellGeom_Spacalv1 *>(_geo_raw);
0160   assert(_geo);
0161   const PHG4CylinderGeom *_layergeom_raw = _layergeo->GetFirstLayerGeom();
0162   assert(_layergeom_raw);
0163   // a special implimentation of PHG4CylinderGeom is required here.
0164   _layergeom = dynamic_cast<const PHG4CylinderGeom_Spacalv3 *>(_layergeom_raw);
0165   assert(_layergeom);
0166   m_geomsetup = true;
0167   return 0;
0168 }
0169
0170 bool PHG4SpacalSteppingAction::NoHitSteppingAction(const G4Step *aStep)
0171 {
0172   // get volume of the current step
0173   G4VPhysicalVolume *volume = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
0174   int isactive = m_Detector->IsInCylinderActive(volume);
0175   if (isactive > PHG4SpacalDetector::INACTIVE)
0176   {
0177     G4StepPoint *prePoint = aStep->GetPreStepPoint();
0178     G4StepPoint *postPoint = aStep->GetPostStepPoint();
0179     // time window cut
0180     double pretime = prePoint->GetGlobalTime() / nanosecond;
0181     double posttime = postPoint->GetGlobalTime() / nanosecond;
0182     if (posttime < m_tmin || pretime > m_tmax)
0183     {
0184       return false;
0185     }
0186     if ((posttime - pretime) > m_dt)
0187     {
0188       return false;
0189     }
0190
0191     int scint_id = -1;
0192
0193     if (  //
0194         m_Detector->get_geom()->get_config() == PHG4SpacalDetector::SpacalGeom_t::kFullProjective_2DTaper ||
0195         m_Detector->get_geom()->get_config() == PHG4SpacalDetector::SpacalGeom_t::kFullProjective_2DTaper_SameLengthFiberPerTower ||
0196         m_Detector->get_geom()->get_config() == PHG4SpacalDetector::SpacalGeom_t::kFullProjective_2DTaper_Tilted ||
0197         m_Detector->get_geom()->get_config() == PHG4SpacalDetector::SpacalGeom_t::kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower  //
0198     )
0199     {
0200       // SPACAL ID that is associated with towers
0201       int sector_ID = 0;
0202       int tower_ID = 0;
0203       int fiber_ID = 0;
0204
0205       if (isactive == PHG4SpacalDetector::FIBER_CORE)
0206       {
0207         fiber_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
0208         tower_ID = prePoint->GetTouchable()->GetReplicaNumber(2);
0209         sector_ID = prePoint->GetTouchable()->GetReplicaNumber(3);
0210       }
0211
0212       else
0213       {
0214         return false;
0215       }
0216
0217       // compact the tower/sector/fiber ID into 32 bit scint_id, so we could save some space for SPACAL hits
0218       scint_id = PHG4CylinderGeom_Spacalv3::scint_id_coder(sector_ID, tower_ID, fiber_ID).scint_ID;
0219     }
0220     else
0221     {
0222       // other configuraitons
0223       if (isactive == PHG4SpacalDetector::FIBER_CORE)
0224       {
0225         scint_id = prePoint->GetTouchable()->GetReplicaNumber(2);
0226       }
0227       else
0228       {
0229         return false;
0230       }
0231     }
0232     // decode scint_id
0233     PHG4CylinderGeom_Spacalv3::scint_id_coder decoder(scint_id);
0234
0235     // convert to z_ID, phi_ID
0236     std::pair<int, int> tower_z_phi_ID = _layergeom->get_tower_z_phi_ID(decoder.tower_ID, decoder.sector_ID);
0237     const int &tower_ID_z = tower_z_phi_ID.first;
0238     const int &tower_ID_phi = tower_z_phi_ID.second;
0239
0240     PHG4CylinderGeom_Spacalv3::tower_map_t::const_iterator it_tower =
0241         _layergeom->get_sector_tower_map().find(decoder.tower_ID);
0242     assert(it_tower != _layergeom->get_sector_tower_map().end());
0243
0244     // convert tower_ID_z to to eta bin number
0245     int etabin = -1;
0246     try
0247     {
0248       etabin = _geo->get_etabin_block(tower_ID_z);  // block eta bin
0249     }
0250     catch (std::exception &e)
0251     {
0252       std::cout << "Print cell geometry:" << std::endl;
0253       _geo->identify();
0254       std::cout << "Print scint_id_coder:" << std::endl;
0255       decoder.identify();
0256       std::cout << "PHG4SpacalSteppingAction::UserSteppingAction::"
0257                 << " - Fatal Error - " << e.what() << std::endl;
0258       exit(1);
0259     }
0260
0261     const int sub_tower_ID_x = it_tower->second.get_sub_tower_ID_x(decoder.fiber_ID);
0262     const int sub_tower_ID_y = it_tower->second.get_sub_tower_ID_y(decoder.fiber_ID);
0263     unsigned short etabinshort = etabin * _layergeom->get_n_subtower_eta() + sub_tower_ID_y;
0264     unsigned short phibin = tower_ID_phi * _layergeom->get_n_subtower_phi() + sub_tower_ID_x;
0265
0266     // get light yield
0267     double light_yield = GetVisibleEnergyDeposition(aStep);
0268
0269     if (light_collection_model.use_fiber_model())
0270     {
0271       const G4TouchableHandle &theTouchable0 = prePoint->GetTouchableHandle();
0272       const G4ThreeVector &worldPosition0 = prePoint->GetPosition();
0273       G4ThreeVector localPosition = theTouchable0->GetHistory()->GetTopTransform().TransformPoint(worldPosition0);
0274       const double localz0 = localPosition.z();
0275       // post point
0276       const G4TouchableHandle &theTouchable1 = postPoint->GetTouchableHandle();
0277       const G4ThreeVector &worldPosition1 = postPoint->GetPosition();
0278       localPosition = theTouchable1->GetHistory()->GetTopTransform().TransformPoint(worldPosition1);
0279       const double localz1 = localPosition.z();
0280
0281       const double z = 0.5 * (localz0 + localz1);
0282       assert(not std::isnan(z));
0283
0284       light_yield *= light_collection_model.get_fiber_transmission(z);
0285     }
0286
0287     // light yield correction from light guide collection efficiency:
0288     if (light_collection_model.use_fiber_model())
0289     {
0290       const double x = it_tower->second.get_position_fraction_x_in_sub_tower(decoder.fiber_ID);
0291       const double y = it_tower->second.get_position_fraction_y_in_sub_tower(decoder.fiber_ID);
0292
0293       light_yield *= light_collection_model.get_light_guide_efficiency(x, y);
0294     }
0295     unsigned int tower_key = TowerInfoDefs::encode_emcal(etabinshort, phibin);
0296     m_CaloInfoContainer->get_tower_at_key(tower_key)->set_energy(m_CaloInfoContainer->get_tower_at_key(tower_key)->get_energy() + light_yield);
0297
0298     // set keep for the track
0299     const G4Track *aTrack = aStep->GetTrack();
0300     if (light_yield > 0)
0301     {
0302       if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
0303       {
0304         if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
0305         {
0306           pp->SetKeep(1);  // we want to keep the track
0307         }
0308       }
0309     }
0310     return true;
0311   }
0312   else
0313   {
0314     return false;
0315   }
0316 }
0317
0318 //____________________________________________________________________________..
0319 bool PHG4SpacalSteppingAction::UserSteppingAction(const G4Step *aStep, bool /*was_used*/)
0320 {
0321   if (!m_doG4Hit)
0322   {
0323     return NoHitSteppingAction(aStep);
0324   }
0325
0326   // get volume of the current step
0327   G4VPhysicalVolume *volume = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
0328
0329   // collect energy and track length step by step
0330   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0331   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0332
0333   const G4Track *aTrack = aStep->GetTrack();
0334
0335   int layer_id = m_Detector->get_Layer();
0336   // make sure we are in a volume
0337   // IsInCylinderActive returns the number of the scintillator
0338   // slat which has fired
0339   int isactive = m_Detector->IsInCylinderActive(volume);
0340   if (isactive > PHG4SpacalDetector::INACTIVE)
0341   {
0342     bool geantino = false;
0343     // the check for the pdg code speeds things up, I do not want to make
0344     // an expensive string compare for every track when we know
0345     // geantino or chargedgeantino has pid=0
0346     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
0347     {
0348       geantino = true;
0349     }
0350     G4StepPoint *prePoint = aStep->GetPreStepPoint();
0351     G4StepPoint *postPoint = aStep->GetPostStepPoint();
0352     int scint_id = -1;
0353
0354     if (  //
0355         m_Detector->get_geom()->get_config() == PHG4SpacalDetector::SpacalGeom_t::kFullProjective_2DTaper ||
0356         m_Detector->get_geom()->get_config() == PHG4SpacalDetector::SpacalGeom_t::kFullProjective_2DTaper_SameLengthFiberPerTower ||
0357         m_Detector->get_geom()->get_config() == PHG4SpacalDetector::SpacalGeom_t::kFullProjective_2DTaper_Tilted ||
0358         m_Detector->get_geom()->get_config() == PHG4SpacalDetector::SpacalGeom_t::kFullProjective_2DTaper_Tilted_SameLengthFiberPerTower  //
0359     )
0360     {
0361       // SPACAL ID that is associated with towers
0362       int sector_ID = 0;
0363       int tower_ID = 0;
0364       int fiber_ID = 0;
0365
0366       if (isactive == PHG4SpacalDetector::FIBER_CORE)
0367       {
0368         fiber_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
0369         tower_ID = prePoint->GetTouchable()->GetReplicaNumber(2);
0370         sector_ID = prePoint->GetTouchable()->GetReplicaNumber(3);
0371       }
0372
0373       else if (isactive == PHG4SpacalDetector::FIBER_CLADING)
0374       {
0375         fiber_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
0376         tower_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
0377         sector_ID = prePoint->GetTouchable()->GetReplicaNumber(2);
0378       }
0379
0380       else if (isactive == PHG4SpacalDetector::ABSORBER)
0381       {
0382         tower_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
0383         sector_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
0384       }
0385
0386       else if (isactive == PHG4SpacalDetector::SUPPORT)
0387       {
0388         tower_ID = prePoint->GetTouchable()->GetReplicaNumber(0);
0389         sector_ID = prePoint->GetTouchable()->GetReplicaNumber(1);
0390 // NOLINTNEXTLINE(hicpp-signed-bitwise)
0391         fiber_ID = (1 << (PHG4CylinderGeom_Spacalv3::scint_id_coder::kfiber_bit)) - 1;  // use max fiber ID to flag for support strucrtures.
0392
0393         //        std::cout <<"PHG4SpacalSteppingAction::UserSteppingAction - SUPPORT tower_ID = "<<tower_ID<<std::endl;
0394       }
0395
0396       // compact the tower/sector/fiber ID into 32 bit scint_id, so we could save some space for SPACAL hits
0397       scint_id = PHG4CylinderGeom_Spacalv3::scint_id_coder(sector_ID, tower_ID, fiber_ID).scint_ID;
0398     }
0399     else
0400     {
0401       // other configuraitons
0402       if (isactive == PHG4SpacalDetector::FIBER_CORE)
0403       {
0404         scint_id = prePoint->GetTouchable()->GetReplicaNumber(2);
0405       }
0406       else if (isactive == PHG4SpacalDetector::FIBER_CLADING)
0407       {
0408         scint_id = prePoint->GetTouchable()->GetReplicaNumber(1);
0409       }
0410       else
0411       {
0412         scint_id = prePoint->GetTouchable()->GetReplicaNumber(0);
0413       }
0414     }
0415
0416     //       std::cout << "track id " << aTrack->GetTrackID() << std::endl;
0417     //        std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
0418     //        std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
0419     switch (prePoint->GetStepStatus())
0420     {
0421     case fGeomBoundary:
0422     case fUndefined:
0423       // if previous hit was saved, hit pointer was set to nullptr
0424       // and we have to make a new one
0425       if (!m_Hit)
0426       {
0427         m_Hit = new PHG4Hitv1();
0428       }
0429       m_Hit->set_layer((unsigned int) layer_id);
0430       m_Hit->set_scint_id(scint_id);  // isactive contains the scintillator slat id
0431       // here we set the entrance values in cm
0432       m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0433       m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0434       m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0435
0436       // time in ns
0437       m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0438       // set the track ID
0439       m_Hit->set_trkid(aTrack->GetTrackID());
0440       m_SaveTrackid = aTrack->GetTrackID();
0441       // set the initial energy deposit
0442       m_Hit->set_edep(0);
0443       // Now add the hit
0444       if (isactive == PHG4SpacalDetector::FIBER_CORE)  // the slat ids start with zero
0445       {
0446         // store all pre local coordinates
0447         StoreLocalCoordinate(m_Hit, aStep, true, false);
0448         m_Hit->set_eion(0);  // only implemented for v5 otherwise empty
0449         m_Hit->set_light_yield(0);
0450         m_CurrentHitContainer = m_HitContainer;
0451       }
0452       else
0453       {
0454         m_CurrentHitContainer = m_AbsorberHitContainer;
0455       }
0456       if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
0457       {
0458         if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
0459         {
0460           m_Hit->set_trkid(pp->GetUserTrackId());
0461           m_Hit->set_shower_id(pp->GetShower()->get_id());
0462           m_CurrentShower = pp->GetShower();
0463         }
0464       }
0465
0466       if (m_Hit->get_z(0) > get_zmax() || m_Hit->get_z(0) < get_zmin())
0467       {
0468         std::cout << "PHG4SpacalSteppingAction: hit outside acceptance, layer: "
0469                   << layer_id << std::endl;
0470         m_Hit->identify();
0471       }
0472       break;
0473     default:
0474       break;
0475     }
0476     // some sanity checks for inconsistencies
0477     // check if this hit was created, if not print out last post step status
0478     if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
0479     {
0480       std::cout << GetName() << ": hit was not created" << std::endl;
0481       std::cout << "prestep status: " << prePoint->GetStepStatus()
0482                 << ", last post step status: " << m_SavePostStepStatus << std::endl;
0483       exit(1);
0484     }
0485     m_SavePostStepStatus = postPoint->GetStepStatus();
0486     // check if track id matches the initial one when the hit was created
0487     if (aTrack->GetTrackID() != m_SaveTrackid)
0488     {
0489       std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
0490       std::cout << "saved track: " << m_SaveTrackid
0491                 << ", current trackid: " << aTrack->GetTrackID()
0492                 << std::endl;
0493       exit(1);
0494     }
0495     // here we just update the exit values, it will be overwritten
0496     // for every step until we leave the volume or the particle
0497     // ceases to exist
0498     m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0499     m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0500     m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0501
0502     m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0503     // sum up the energy to get total deposited
0504     m_Hit->set_edep(m_Hit->get_edep() + edep);
0505
0506     if (isactive == PHG4SpacalDetector::FIBER_CORE)  // only for active areas
0507     {
0508       // store all pre local coordinates
0509       StoreLocalCoordinate(m_Hit, aStep, false, true);
0510
0511       m_Hit->set_eion(m_Hit->get_eion() + eion);
0512
0513       double light_yield = GetVisibleEnergyDeposition(aStep);
0514
0515       static bool once = true;
0516       if (once and edep > 0)
0517       {
0518         once = false;
0519
0520         if (Verbosity() > 0)
0521         {
0522           std::cout << "PHG4SpacalSteppingAction::UserSteppingAction::"
0523                     //
0524                     << m_Detector->GetName() << " - "
0525                     << " use scintillating light model at each Geant4 steps. "
0526                     << "First step: "
0527                     << "Material = "
0528                     << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
0529                     << ", "
0530                     << "Birk Constant = "
0531                     << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetIonisation()->GetBirksConstant()
0532                     << ","
0533                     << "edep = " << edep << ", "
0534                     << "eion = " << eion
0535                     << ", "
0536                     << "light_yield = " << light_yield << std::endl;
0537         }
0538       }
0539
0540       m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
0541     }
0542
0543     if (m_Hit->get_z(1) > get_zmax() || m_Hit->get_z(1) < get_zmin())
0544     {
0545       std::cout << "PHG4SpacalSteppingAction: hit outside acceptance get_zmin() "
0546                 << get_zmin() << ", get_zmax() " << get_zmax() << " at exit"
0547                 << std::endl;
0548       m_Hit->identify();
0549     }
0550     if (geantino)
0551     {
0552       m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0553                             //          m_Hit->set_eion(-1);
0554     }
0555     if (edep > 0)
0556     {
0557       if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
0558       {
0559         if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
0560         {
0561           pp->SetKeep(1);  // we want to keep the track
0562         }
0563       }
0564     }
0565     // if any of these conditions is true this is the last step in
0566     // this volume and we need to save the hit
0567     // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0568     // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0569     // (happens when your detector goes outside world volume)
0570     // aTrack->GetTrackStatus() == fStopAndKill: track ends
0571     if (postPoint->GetStepStatus() == fGeomBoundary ||
0572         postPoint->GetStepStatus() == fWorldBoundary ||
0573         postPoint->GetStepStatus() == fAtRestDoItProc ||
0574         aTrack->GetTrackStatus() == fStopAndKill)
0575     {
0576       // save only hits with energy deposit (or -1 for geantino)
0577       if (m_Hit->get_edep())
0578       {
0579         m_CurrentHitContainer->AddHit(layer_id, m_Hit);
0580         if (m_CurrentShower)
0581         {
0582           m_CurrentShower->add_g4hit_id(m_CurrentHitContainer->GetID(), m_Hit->get_hit_id());
0583         }
0584         // ownership has been transferred to container, set to null
0585         // so we will create a new hit for the next track
0586         m_Hit = nullptr;
0587       }
0588       else
0589       {
0590         // if this hit has no energy deposit, just reset it for reuse
0591         // this means we have to delete it in the dtor. If this was
0592         // the last hit we processed the memory is still allocated
0593         m_Hit->Reset();
0594       }
0595     }
0596     // return true to indicate the hit was used
0597     return true;
0598   }
0599   else
0600   {
0601     return false;
0602   }
0603 }
0604
0605 //____________________________________________________________________________..
0606 void PHG4SpacalSteppingAction::SetInterfacePointers(PHCompositeNode *topNode)
0607 {
0608   // we can only play with the geometry node after ConstructMe is called
0609   if ((!m_geomsetup) && (!m_doG4Hit))
0610   {
0611     SetUpGeomNode(topNode);
0612   }
0613   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0614   m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
0615   // if we do not find the node it's messed up.
0616   if (!m_HitContainer)
0617   {
0618     std::cout << "PHG4SpacalSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
0619     gSystem->Exit(1);
0620   }
0621   // this is perfectly fine if absorber hits are disabled
0622   if (!m_AbsorberHitContainer)
0623   {
0624     if (Verbosity() > 0)
0625     {
0626       std::cout << "PHG4SpacalSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
0627     }
0628   }
0629 }
0630
0631 double
0632 PHG4SpacalSteppingAction::get_zmin() const
0633 {
0634   if (!m_Detector)
0635   {
0636     return 0;
0637   }
0638   else
0639   {
0640     return m_Detector->get_geom()->get_zmin() - .0001;
0641   }
0642 }
0643
0644 double
0645 PHG4SpacalSteppingAction::get_zmax() const
0646 {
0647   if (!m_Detector)
0648   {
0649     return 0;
0650   }
0651   else
0652   {
0653     return m_Detector->get_geom()->get_zmax() + .0001;
0654   }
0655 }
0656
0657 void PHG4SpacalSteppingAction::SetHitNodeName(const std::string &type, const std::string &name)
0658 {
0659   if (type == "G4HIT")
0660   {
0661     m_HitNodeName = name;
0662     return;
0663   }
0664   else if (type == "G4HIT_ABSORBER")
0665   {
0666     m_AbsorberNodeName = name;
0667     return;
0668   }
0669   std::cout << "Invalid output hit node type " << type << std::endl;
0670   gSystem->Exit(1);
0671   return;
0672 }