Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4HcalPrototypeSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4HcalPrototypeSteppingAction.cc
1 // This is the steppingaction code for the hcal prototype
2 // created on 1/27/2014, Liang, HeXC
3 //
6 
7 #include <g4main/PHG4Hit.h>
9 #include <g4main/PHG4Hitv1.h>
10 #include <g4main/PHG4Shower.h>
12 
13 #include <phool/getClass.h>
14 
15 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
16 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
17 #include <Geant4/G4Step.hh>
18 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
19 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fUndefined
20 #include <Geant4/G4String.hh> // for G4String
21 #include <Geant4/G4SystemOfUnits.hh> // for cm, GeV, nanosecond
22 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
23 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
24 #include <Geant4/G4Track.hh> // for G4Track
25 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
26 #include <Geant4/G4Types.hh> // for G4int, G4double
27 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
28 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
29 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
30 
31 #include <iostream>
32 #include <string> // for operator+, string, ope...
33 
34 class PHCompositeNode;
35 
36 using namespace std;
37 //____________________________________________________________________________..
39  : PHG4SteppingAction(detector->GetName())
40  , detector_(detector)
41  , hits_(nullptr)
42  , absorberhits_(nullptr)
43  , hit(nullptr)
44 {
45 }
46 
47 //____________________________________________________________________________..
49 {
50  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
51  // get volume of the current step
52  G4VPhysicalVolume* volume = touch->GetVolume();
53 
54  // We simply use ourown scintID to test on this condition
55  // int whichactive = detector_->IsInHcalPrototype(volume);
56 
57  /*
58  if (!whichactive)
59  {
60  return false;
61  }
62  */
63 
64  // collect energy and track length step by step
65  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
66  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
67 
68  const G4Track* aTrack = aStep->GetTrack();
69 
70  // if this block stops everything, just put all kinetic energy into edep
71  if (detector_->IsBlackHole())
72  {
73  edep = aTrack->GetKineticEnergy() / GeV;
74  G4Track* killtrack = const_cast<G4Track*>(aTrack);
75  killtrack->SetTrackStatus(fStopAndKill);
76  }
77 
78  // Add detector element ID into to the hit
79 
80  G4int scintSheetCopyNumber, scintSheetLayerNumber;
81 
82  G4int sectionID = 0;
83  G4int scintID = 0x0;
84  // G4int absID = 0;
85 
86  if (volume->GetName() == "hcal1AbsLayer")
87  {
88  sectionID = 1;
89  // absID = volume->GetCopyNo()+100; // Inner section absorber layer ID's are 100, 101, ..., 115
90  scintID = -1; // scintID is not valid, we set it to -1
91  }
92  else if (volume->GetName() == "hcal2AbsLayer")
93  {
94  sectionID = 2;
95  // absID = volume->GetCopyNo()+200; // Outer section absorber layer ID's are 200, 201, ..., 215
96  scintID = -1; // scintID is not valid, we set it to -1
97  }
98  else if (volume->GetName() == "inner1USheet")
99  {
100  sectionID = 1;
101  scintSheetCopyNumber = volume->GetCopyNo(); // either 1 or 0
102  scintSheetLayerNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(1)->GetCopyNo();
103  scintID = (scintSheetLayerNumber << 2) + (1 << 1) + scintSheetCopyNumber;
104  // absID = -1; // absID is not valid, we set it to -1
105  //G4cout << " **************** scintID: " << scintID << G4endl;
106  }
107  else if (volume->GetName() == "inner2USheet")
108  {
109  sectionID = 1;
110  scintSheetCopyNumber = volume->GetCopyNo(); // either 1 or 0
111  scintSheetLayerNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(1)->GetCopyNo();
112  scintID = (scintSheetLayerNumber << 2) + (0 << 1) + scintSheetCopyNumber;
113  // absID = -1; // absID is not valid, we set it to -1
114  //G4cout << " **************** scintID: " << scintID << G4endl;
115  }
116  else if (volume->GetName() == "outer1USheet")
117  {
118  sectionID = 2;
119  scintSheetCopyNumber = volume->GetCopyNo(); // either 1 or 0
120  scintSheetLayerNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(1)->GetCopyNo();
121  scintID = (scintSheetLayerNumber << 2) + (1 << 1) + scintSheetCopyNumber;
122  // absID = -1; // absID is not valid, we set it to -1
123  //G4cout << " **************** scintID: " << scintID << G4endl;
124  }
125  else if (volume->GetName() == "outer2USheet")
126  {
127  sectionID = 2;
128  scintSheetCopyNumber = volume->GetCopyNo(); // either 1 or 0
129  scintSheetLayerNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(1)->GetCopyNo();
130  scintID = (scintSheetLayerNumber << 2) + (0 << 1) + scintSheetCopyNumber;
131  // absID = -1; // absID is not valid, we set it to -1
132  //G4cout << " **************** scintID: " << scintID << G4endl;
133  }
134  else
135  {
136  return false;
137  }
138 
139  // make sure we are in a volume
140  if (detector_->IsActive())
141  {
142  bool geantino = false;
143  // the check for the pdg code speeds things up, I do not want to make
144  // an expensive string compare for every track when we know
145  // geantino or chargedgeantino has pid=0
146  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
147  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
148  {
149  geantino = true;
150  }
151 
152  G4StepPoint* prePoint = aStep->GetPreStepPoint();
153  G4StepPoint* postPoint = aStep->GetPostStepPoint();
154  // cout << "track id " << aTrack->GetTrackID() << endl;
155  // cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
156  // cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
157  switch (prePoint->GetStepStatus())
158  {
159  case fGeomBoundary:
160  case fUndefined:
161  hit = new PHG4Hitv1();
162  hit->set_layer((unsigned int) sectionID);
163  hit->set_scint_id(scintID);
164  //here we set the entrance values in cm
165  hit->set_x(0, prePoint->GetPosition().x() / cm);
166  hit->set_y(0, prePoint->GetPosition().y() / cm);
167  hit->set_z(0, prePoint->GetPosition().z() / cm);
168  // time in ns
169  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
170  //set the track ID
171  {
172  hit->set_trkid(aTrack->GetTrackID());
173  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
174  {
175  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
176  {
177  hit->set_trkid(pp->GetUserTrackId());
178  hit->set_shower_id(pp->GetShower()->get_id());
179  }
180  }
181  }
182 
183  //set the initial energy deposit
184  hit->set_edep(0);
185  hit->set_eion(0);
186  if (scintID > 0) // return of isinHcalTestDetector, > 0 hit in scintillator, < 0 hit in absorber
187  {
188  // Now add the hit
189  hits_->AddHit(sectionID, hit);
190 
191  {
192  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
193  {
194  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
195  {
196  pp->GetShower()->add_g4hit_id(hits_->GetID(), hit->get_hit_id());
197  }
198  }
199  }
200  }
201  else
202  {
203  absorberhits_->AddHit(sectionID, hit);
204 
205  {
206  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
207  {
208  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
209  {
210  pp->GetShower()->add_g4hit_id(absorberhits_->GetID(), hit->get_hit_id());
211  }
212  }
213  }
214  }
215  break;
216  default:
217  break;
218  }
219 
220  // here we just update the exit values, it will be overwritten
221  // for every step until we leave the volume or the particle
222  // ceases to exist
223  hit->set_x(1, postPoint->GetPosition().x() / cm);
224  hit->set_y(1, postPoint->GetPosition().y() / cm);
225  hit->set_z(1, postPoint->GetPosition().z() / cm);
226 
227  hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
228  //sum up the energy to get total deposited
229  hit->set_edep(hit->get_edep() + edep);
230  hit->set_eion(hit->get_eion() + eion);
231  if (geantino)
232  {
233  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
234  hit->set_eion(-1);
235  }
236  if (edep > 0)
237  {
238  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
239  {
240  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
241  {
242  pp->SetKeep(1); // we want to keep the track
243  }
244  }
245  }
246 
247  // hit->identify();
248  // return true to indicate the hit was used
249  return true;
250  }
251  else
252  {
253  return false;
254  }
255 }
256 
257 //____________________________________________________________________________..
259 {
260  string hitnodename;
261  string absorbernodename;
262  if (detector_->SuperDetector() != "NONE")
263  {
264  hitnodename = "G4HIT_" + detector_->SuperDetector();
265  absorbernodename = "G4HIT_ABSORBER_" + detector_->SuperDetector();
266  }
267  else
268  {
269  hitnodename = "G4HIT_" + detector_->GetName();
270  absorbernodename = "G4HIT_ABSORBER_" + detector_->GetName();
271  }
272 
273  //now look for the map and grab a pointer to it.
274  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
275  absorberhits_ = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
276 
277  // if we do not find the node it's messed up.
278  if (!hits_)
279  {
280  std::cout << "PHG4HcalPrototypeSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
281  }
282  if (!absorberhits_)
283  {
284  if (Verbosity() > 0)
285  {
286  cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
287  }
288  }
289 }