ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
WLSSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file WLSSteppingAction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
29 //
30 //
31 #include "G4Run.hh"
32 #include "G4Step.hh"
33 #include "G4Track.hh"
34 #include "G4StepPoint.hh"
35 #include "G4TrackStatus.hh"
36 #include "G4VPhysicalVolume.hh"
37 #include "G4ParticleDefinition.hh"
38 
39 #include "WLSSteppingAction.hh"
40 #include "WLSDetectorConstruction.hh"
41 #include "WLSSteppingActionMessenger.hh"
42 #include "WLSPhotonDetSD.hh"
43 
44 #include "G4ParticleTypes.hh"
45 
46 #include "WLSUserTrackInformation.hh"
47 
48 #include "G4ProcessManager.hh"
49 #include "G4OpBoundaryProcess.hh"
50 
51 #include "G4RunManager.hh"
52 #include "G4SDManager.hh"
53 #include "G4UImanager.hh"
54 
55 #include "G4ThreeVector.hh"
56 #include "G4ios.hh"
57 #include "G4SystemOfUnits.hh"
58 #include <sstream>
59 
60 // Purpose: Save relevant information into User Track Information
61 
62 static const G4ThreeVector ZHat = G4ThreeVector(0.0,0.0,1.0);
63 
64 G4int WLSSteppingAction::fMaxRndmSave = 10000;
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 
68 WLSSteppingAction::WLSSteppingAction(WLSDetectorConstruction* detector)
69  : fDetector(detector)
70 {
71  fSteppingMessenger = new WLSSteppingActionMessenger(this);
72 
73  fCounterEnd = 0;
74  fCounterMid = 0;
75  fBounceLimit = 100000;
76 
77  fOpProcess = NULL;
78 
79  ResetCounters();
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 
84 WLSSteppingAction::~WLSSteppingAction()
85 {
86  delete fSteppingMessenger;
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 void WLSSteppingAction::SetBounceLimit(G4int i) {fBounceLimit = i;}
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94 
95 G4int WLSSteppingAction::GetNumberOfBounces() {return fCounterBounce;}
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
99 G4int WLSSteppingAction::GetNumberOfClad1Bounces() {return fCounterClad1Bounce;}
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 
103 G4int WLSSteppingAction::GetNumberOfClad2Bounces() {return fCounterClad2Bounce;}
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106 
107 G4int WLSSteppingAction::GetNumberOfWLSBounces() {return fCounterWLSBounce;}
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
111 G4int WLSSteppingAction::ResetSuccessCounter() {
112  G4int temp = fCounterEnd; fCounterEnd = 0; return temp;
113 }
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116 
117 inline void WLSSteppingAction::SaveRandomStatus(G4String subDir)
118 // save the random status into a sub-directory
119 // Pre: subDir must be empty or ended with "/"
120 {
121 
122  // don't save if the maximum amount has been reached
123  if (WLSSteppingAction::fMaxRndmSave == 0) return;
124 
125  G4RunManager* theRunManager = G4RunManager::GetRunManager();
126  G4String randomNumberStatusDir = theRunManager->GetRandomNumberStoreDir();
127 
128  G4String fileIn = randomNumberStatusDir + "currentEvent.rndm";
129 
130  std::ostringstream os;
131 
132  os << "run" << theRunManager->GetCurrentRun()->GetRunID() << "evt"
133  << theRunManager->GetCurrentEvent()->GetEventID() << ".rndm" << '\0';
134 
135  G4String fileOut = randomNumberStatusDir + subDir + os.str();
136 
137  G4String copCmd = "/control/shell cp "+fileIn+" "+fileOut;
138  G4UImanager::GetUIpointer()->ApplyCommand(copCmd);
139 
140  WLSSteppingAction::fMaxRndmSave--;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
145 void WLSSteppingAction::UserSteppingAction(const G4Step* theStep)
146 {
147  G4Track* theTrack = theStep->GetTrack();
148  WLSUserTrackInformation* trackInformation
149  = (WLSUserTrackInformation*)theTrack->GetUserInformation();
150 
151  G4StepPoint* thePrePoint = theStep->GetPreStepPoint();
152  G4StepPoint* thePostPoint = theStep->GetPostStepPoint();
153 
154  G4VPhysicalVolume* thePrePV = thePrePoint->GetPhysicalVolume();
155  G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
156 
157  G4String thePrePVname = " ";
158  G4String thePostPVname = " ";
159 
160  if (thePostPV) {
161  thePrePVname = thePrePV->GetName();
162  thePostPVname = thePostPV->GetName();
163  }
164 
165  //Recording data for start
166  if (theTrack->GetParentID()==0) {
167  //This is a primary track
168  if ( theTrack->GetCurrentStepNumber() == 1 ) {
169 // G4double x = theTrack->GetVertexPosition().x();
170 // G4double y = theTrack->GetVertexPosition().y();
171 // G4double z = theTrack->GetVertexPosition().z();
172 // G4double pz = theTrack->GetVertexMomentumDirection().z();
173 // G4double fInitTheta =
174 // theTrack->GetVertexMomentumDirection().angle(ZHat);
175  }
176  }
177 
178  // Retrieve the status of the photon
179  G4OpBoundaryProcessStatus theStatus = Undefined;
180 
181  G4ProcessManager* OpManager =
182  G4OpticalPhoton::OpticalPhoton()->GetProcessManager();
183 
184  if (OpManager) {
185  G4int MAXofPostStepLoops =
186  OpManager->GetPostStepProcessVector()->entries();
187  G4ProcessVector* fPostStepDoItVector =
188  OpManager->GetPostStepProcessVector(typeDoIt);
189 
190  for ( G4int i=0; i<MAXofPostStepLoops; i++) {
191  G4VProcess* fCurrentProcess = (*fPostStepDoItVector)[i];
192  fOpProcess = dynamic_cast<G4OpBoundaryProcess*>(fCurrentProcess);
193  if (fOpProcess) { theStatus = fOpProcess->GetStatus(); break;}
194  }
195  }
196 
197  // Find the skewness of the ray at first change of boundary
198  if ( fInitGamma == -1 &&
199  (theStatus == TotalInternalReflection
200  || theStatus == FresnelReflection
201  || theStatus == FresnelRefraction)
202  && trackInformation->IsStatus(InsideOfFiber) ) {
203 
204  G4double px = theTrack->GetVertexMomentumDirection().x();
205  G4double py = theTrack->GetVertexMomentumDirection().y();
206  G4double x = theTrack->GetPosition().x();
207  G4double y = theTrack->GetPosition().y();
208 
209  fInitGamma = x * px + y * py;
210 
211  fInitGamma =
212  fInitGamma / std::sqrt(px*px + py*py) / std::sqrt(x*x + y*y);
213 
214  fInitGamma = std::acos(fInitGamma*rad);
215 
216  if ( fInitGamma / deg > 90.0) { fInitGamma = 180 * deg - fInitGamma;}
217  }
218  // Record Photons that missed the photon detector but escaped from readout
219  if ( !thePostPV && trackInformation->IsStatus(EscapedFromReadOut) ) {
220 // UpdateHistogramSuccess(thePostPoint,theTrack);
221  ResetCounters();
222 
223  return;
224  }
225 
226  // Assumed photons are originated at the fiber OR
227  // the fiber is the first material the photon hits
228  switch (theStatus) {
229 
230  // Exiting the fiber
231  case FresnelRefraction:
232  case SameMaterial:
233 
234  G4bool isFiber;
235  isFiber = thePostPVname == "WLSFiber"
236  || thePostPVname == "Clad1"
237  || thePostPVname == "Clad2";
238 
239  if ( isFiber ) {
240 
241  if (trackInformation->IsStatus(OutsideOfFiber))
242  trackInformation->AddStatusFlag(InsideOfFiber);
243 
244  // Set the Exit flag when the photon refracted out of the fiber
245  } else if (trackInformation->IsStatus(InsideOfFiber)) {
246 
247  // EscapedFromReadOut if the z position is the same as fiber's end
248  if (theTrack->GetPosition().z() == fDetector->GetWLSFiberEnd())
249  {
250  trackInformation->AddStatusFlag(EscapedFromReadOut);
251  fCounterEnd++;
252  }
253  else // Escaped from side
254  {
255  trackInformation->AddStatusFlag(EscapedFromSide);
256  trackInformation->SetExitPosition(theTrack->GetPosition());
257 
258 // UpdateHistogramEscape(thePostPoint,theTrack);
259 
260  fCounterMid++;
261  ResetCounters();
262  }
263 
264  trackInformation->AddStatusFlag(OutsideOfFiber);
265  trackInformation->SetExitPosition(theTrack->GetPosition());
266 
267  }
268 
269  return;
270 
271  // Internal Reflections
272  case TotalInternalReflection:
273 
274  // Kill the track if it's number of bounces exceeded the limit
275  if (fBounceLimit > 0 && fCounterBounce >= fBounceLimit)
276  {
277  theTrack->SetTrackStatus(fStopAndKill);
278  trackInformation->AddStatusFlag(murderee);
279  ResetCounters();
280  G4cout << "\n Bounce Limit Exceeded" << G4endl;
281  return;
282  }
283  break;
284 
285  case FresnelReflection:
286 
287  fCounterBounce++;
288 
289  if ( thePrePVname == "WLSFiber") fCounterWLSBounce++;
290 
291  else if ( thePrePVname == "Clad1") fCounterClad1Bounce++;
292 
293  else if ( thePrePVname == "Clad2") fCounterClad2Bounce++;
294 
295  // Determine if the photon has reflected off the read-out end
296  if (theTrack->GetPosition().z() == fDetector->GetWLSFiberEnd())
297  {
298  if (!trackInformation->IsStatus(ReflectedAtReadOut) &&
299  trackInformation->IsStatus(InsideOfFiber))
300  {
301  trackInformation->AddStatusFlag(ReflectedAtReadOut);
302 
303  if (fDetector->IsPerfectFiber() &&
304  theStatus == TotalInternalReflection)
305  {
306  theTrack->SetTrackStatus(fStopAndKill);
307  trackInformation->AddStatusFlag(murderee);
308 // UpdateHistogramReflect(thePostPoint,theTrack);
309  ResetCounters();
310  return;
311  }
312  }
313  }
314  return;
315 
316  // Reflection of the mirror
317  case LambertianReflection:
318  case LobeReflection:
319  case SpikeReflection:
320 
321  // Check if it hits the mirror
322  if ( thePostPVname == "Mirror" )
323  trackInformation->AddStatusFlag(ReflectedAtMirror);
324 
325  return;
326 
327  // Detected by a detector
328  case Detection:
329 
330  // Check if the photon hits the detector and process the hit if it does
331  if ( thePostPVname == "PhotonDet" ) {
332 
333  G4SDManager* SDman = G4SDManager::GetSDMpointer();
334  G4String SDname="WLS/PhotonDet";
335  WLSPhotonDetSD* mppcSD =
336  (WLSPhotonDetSD*)SDman->FindSensitiveDetector(SDname);
337 
338  if (mppcSD) mppcSD->ProcessHits_constStep(theStep,NULL);
339 
340  // Record Photons that escaped at the end
341 // if (trackInformation->IsStatus(EscapedFromReadOut))
342 // UpdateHistogramSuccess(thePostPoint,theTrack);
343 
344  // Stop Tracking when it hits the detector's surface
345  ResetCounters();
346  theTrack->SetTrackStatus(fStopAndKill);
347 
348  return;
349  }
350  break;
351 
352  default: break;
353 
354  }
355 
356  // Check for absorbed photons
357  if (theTrack->GetTrackStatus() != fAlive &&
358  trackInformation->IsStatus(InsideOfFiber))
359  {
360 // UpdateHistogramAbsorb(thePostPoint,theTrack);
361  ResetCounters();
362  return;
363  }
364 
365 }