Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4CEmcTestBeamSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4CEmcTestBeamSteppingAction.cc
2 
4 
5 #include <g4main/PHG4Hit.h>
7 #include <g4main/PHG4Hitv1.h>
8 #include <g4main/PHG4Shower.h>
9 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
11 
12 #include <phool/getClass.h>
13 
14 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
15 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
16 #include <Geant4/G4Step.hh>
17 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
18 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fUndefined
19 #include <Geant4/G4String.hh> // for G4String
20 #include <Geant4/G4SystemOfUnits.hh>
21 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
22 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
23 #include <Geant4/G4Track.hh> // for G4Track
24 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
25 #include <Geant4/G4Types.hh> // for G4double
26 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
27 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
28 
29 #include <iostream>
30 #include <string> // for operator+, string, ope...
31 
32 class G4VPhysicalVolume;
33 class PHCompositeNode;
34 
35 using namespace std;
36 //____________________________________________________________________________..
38  : PHG4SteppingAction(detector->GetName())
39  , detector_(detector)
40  , hits_(nullptr)
41  , absorberhits_(nullptr)
42  , hit(nullptr)
43 {
44 }
45 
46 //____________________________________________________________________________..
48 {
49  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
50  // get volume of the current step
51  G4VPhysicalVolume* volume = touch->GetVolume();
52 
53  int whichactive = detector_->IsInCEmcTestBeam(volume);
54 
55  if (!whichactive)
56  {
57  return false;
58  }
59 
60  // collect energy and track length step by step
61  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
62  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
63 
64  const G4Track* aTrack = aStep->GetTrack();
65 
66  // if this block stops everything, just put all kinetic energy into edep
67  if (detector_->IsBlackHole())
68  {
69  edep = aTrack->GetKineticEnergy() / GeV;
70  G4Track* killtrack = const_cast<G4Track*>(aTrack);
71  killtrack->SetTrackStatus(fStopAndKill);
72  }
73 
74  int tower_id = touch->GetCopyNumber(2); // tower number
75  // make sure we are in a volume
76  if (detector_->IsActive())
77  {
78  bool geantino = false;
79  // the check for the pdg code speeds things up, I do not want to make
80  // an expensive string compare for every track when we know
81  // geantino or chargedgeantino has pid=0
82  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
83  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
84  {
85  geantino = true;
86  }
87  G4StepPoint* prePoint = aStep->GetPreStepPoint();
88  G4StepPoint* postPoint = aStep->GetPostStepPoint();
89  // cout << "track id " << aTrack->GetTrackID() << endl;
90  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
91  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
92  switch (prePoint->GetStepStatus())
93  {
94  case fGeomBoundary:
95  case fUndefined:
96  hit = new PHG4Hitv1();
97  hit->set_layer((unsigned int) tower_id);
98  hit->set_scint_id(touch->GetCopyNumber(1)); // the copy number of the sandwich
99  //here we set the entrance values in cm
100  hit->set_x(0, prePoint->GetPosition().x() / cm);
101  hit->set_y(0, prePoint->GetPosition().y() / cm);
102  hit->set_z(0, prePoint->GetPosition().z() / cm);
103  // time in ns
104  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
105  //set the track ID
106  {
107  hit->set_trkid(aTrack->GetTrackID());
108  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
109  {
110  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
111  {
112  hit->set_trkid(pp->GetUserTrackId());
113  hit->set_shower_id(pp->GetShower()->get_id());
114  }
115  }
116  }
117 
118  //set the initial energy deposit
119  hit->set_edep(0);
120  hit->set_eion(0); // only implemented for v5 otherwise empty
121  PHG4HitContainer* hitcontainer;
122  // here we do things which are different between scintillator and absorber hits
123  if (whichactive > 0) // return of isinCEmcTestDetector, > 0 hit in scintillator, < 0 hit in absorber
124  {
125  hitcontainer = hits_;
126  }
127  else
128  {
129  hitcontainer = absorberhits_;
130  }
131  // here we set what is common for scintillator and absorber hits
132  hitcontainer->AddHit(tower_id, hit);
133  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
134  {
135  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
136  {
137  pp->GetShower()->add_g4hit_id(hitcontainer->GetID(), hit->get_hit_id());
138  }
139  }
140  break;
141  default:
142  break;
143  }
144  // here we just update the exit values, it will be overwritten
145  // for every step until we leave the volume or the particle
146  // ceases to exist
147  hit->set_x(1, postPoint->GetPosition().x() / cm);
148  hit->set_y(1, postPoint->GetPosition().y() / cm);
149  hit->set_z(1, postPoint->GetPosition().z() / cm);
150 
151  hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
152  //sum up the energy to get total deposited
153  hit->set_edep(hit->get_edep() + edep);
154  hit->set_eion(hit->get_eion() + eion);
155  if (geantino)
156  {
157  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
158  hit->set_eion(-1);
159  }
160  if (edep > 0)
161  {
162  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
163  {
164  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
165  {
166  pp->SetKeep(1); // we want to keep the track
167  }
168  }
169  }
170 
171  // hit->identify();
172  // return true to indicate the hit was used
173  return true;
174  }
175  else
176  {
177  return false;
178  }
179 }
180 
181 //____________________________________________________________________________..
183 {
184  string hitnodename;
185  string absorbernodename;
186  if (detector_->SuperDetector() != "NONE")
187  {
188  hitnodename = "G4HIT_" + detector_->SuperDetector();
189  absorbernodename = "G4HIT_ABSORBER_" + detector_->SuperDetector();
190  }
191  else
192  {
193  hitnodename = "G4HIT_" + detector_->GetName();
194  absorbernodename = "G4HIT_ABSORBER_" + detector_->GetName();
195  }
196 
197  //now look for the map and grab a pointer to it.
198  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
199  absorberhits_ = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
200 
201  // if we do not find the node it's messed up.
202  if (!hits_)
203  {
204  std::cout << "PHG4CEmcTestBeamSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
205  }
206  if (!absorberhits_)
207  {
208  if (Verbosity() > 0)
209  {
210  cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
211  }
212  }
213 }