| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 /// \file optical/wls/src/WLSSteppingAction.cc 27 /// \file optical/wls/src/WLSSteppingAction.cc
28 /// \brief Implementation of the WLSSteppingAc 28 /// \brief Implementation of the WLSSteppingAction class
29 // 29 //
30 // 30 //
31 #include "WLSSteppingAction.hh" 31 #include "WLSSteppingAction.hh"
32 32
33 #include "WLSDetectorConstruction.hh" 33 #include "WLSDetectorConstruction.hh"
34 #include "WLSEventAction.hh" 34 #include "WLSEventAction.hh"
35 #include "WLSPhotonDetSD.hh" 35 #include "WLSPhotonDetSD.hh"
36 #include "WLSSteppingActionMessenger.hh" 36 #include "WLSSteppingActionMessenger.hh"
37 #include "WLSUserTrackInformation.hh" 37 #include "WLSUserTrackInformation.hh"
38 38
>> 39 #include "G4ios.hh"
39 #include "G4OpBoundaryProcess.hh" 40 #include "G4OpBoundaryProcess.hh"
40 #include "G4OpticalPhoton.hh" 41 #include "G4OpticalPhoton.hh"
41 #include "G4ProcessManager.hh" 42 #include "G4ProcessManager.hh"
42 #include "G4Run.hh" 43 #include "G4Run.hh"
43 #include "G4SDManager.hh" 44 #include "G4SDManager.hh"
44 #include "G4Step.hh" 45 #include "G4Step.hh"
45 #include "G4StepPoint.hh" 46 #include "G4StepPoint.hh"
46 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
47 #include "G4ThreeVector.hh" 48 #include "G4ThreeVector.hh"
48 #include "G4Track.hh" 49 #include "G4Track.hh"
49 #include "G4TrackStatus.hh" 50 #include "G4TrackStatus.hh"
50 #include "G4UImanager.hh" 51 #include "G4UImanager.hh"
51 #include "G4VPhysicalVolume.hh" 52 #include "G4VPhysicalVolume.hh"
52 #include "G4ios.hh" <<
53 53
54 // Purpose: Save relevant information into Use 54 // Purpose: Save relevant information into User Track Information
55 55
56 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 57
58 WLSSteppingAction::WLSSteppingAction(WLSDetect << 58 WLSSteppingAction::WLSSteppingAction(WLSDetectorConstruction* detector,
59 : fDetector(detector), fEventAction(event) << 59 WLSEventAction* event)
>> 60 : fDetector(detector)
>> 61 , fEventAction(event)
60 { 62 {
61 fSteppingMessenger = new WLSSteppingActionMe 63 fSteppingMessenger = new WLSSteppingActionMessenger(this);
>> 64
>> 65 fCounterEnd = 0;
>> 66 fCounterMid = 0;
>> 67 fBounceLimit = 100000;
>> 68
>> 69 fOpProcess = nullptr;
>> 70
62 ResetCounters(); 71 ResetCounters();
63 } 72 }
64 73
65 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66 75
67 WLSSteppingAction::~WLSSteppingAction() << 76 WLSSteppingAction::~WLSSteppingAction() { delete fSteppingMessenger; }
68 { <<
69 delete fSteppingMessenger; <<
70 } <<
71 77
72 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 79
74 void WLSSteppingAction::SetBounceLimit(G4int i << 80 void WLSSteppingAction::SetBounceLimit(G4int i) { fBounceLimit = i; }
75 { <<
76 fBounceLimit = i; <<
77 } <<
78 81
79 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 83
81 G4int WLSSteppingAction::GetNumberOfBounces() << 84 G4int WLSSteppingAction::GetNumberOfBounces() { return fCounterBounce; }
82 { <<
83 return fCounterBounce; <<
84 } <<
85 85
86 //....oooOO0OOooo........oooOO0OOooo........oo 86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 87
88 G4int WLSSteppingAction::GetNumberOfClad1Bounc 88 G4int WLSSteppingAction::GetNumberOfClad1Bounces()
89 { 89 {
90 return fCounterClad1Bounce; 90 return fCounterClad1Bounce;
91 } 91 }
92 92
93 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94 94
95 G4int WLSSteppingAction::GetNumberOfClad2Bounc 95 G4int WLSSteppingAction::GetNumberOfClad2Bounces()
96 { 96 {
97 return fCounterClad2Bounce; 97 return fCounterClad2Bounce;
98 } 98 }
99 99
100 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 101
102 G4int WLSSteppingAction::GetNumberOfWLSBounces << 102 G4int WLSSteppingAction::GetNumberOfWLSBounces() { return fCounterWLSBounce; }
103 { <<
104 return fCounterWLSBounce; <<
105 } <<
106 103
107 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 105
109 G4int WLSSteppingAction::ResetSuccessCounter() 106 G4int WLSSteppingAction::ResetSuccessCounter()
110 { 107 {
111 G4int temp = fCounterEnd; << 108 G4int temp = fCounterEnd;
112 fCounterEnd = 0; 109 fCounterEnd = 0;
113 return temp; 110 return temp;
114 } 111 }
115 112
116 //....oooOO0OOooo........oooOO0OOooo........oo 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 114
118 void WLSSteppingAction::UserSteppingAction(con 115 void WLSSteppingAction::UserSteppingAction(const G4Step* theStep)
119 { 116 {
120 G4Track* theTrack = theStep->GetTrack(); 117 G4Track* theTrack = theStep->GetTrack();
121 auto trackInformation = (WLSUserTrackInforma << 118 WLSUserTrackInformation* trackInformation =
>> 119 (WLSUserTrackInformation*) theTrack->GetUserInformation();
122 120
123 G4StepPoint* thePrePoint = theStep->GetPreSt << 121 G4StepPoint* thePrePoint = theStep->GetPreStepPoint();
124 G4StepPoint* thePostPoint = theStep->GetPost 122 G4StepPoint* thePostPoint = theStep->GetPostStepPoint();
125 123
126 G4VPhysicalVolume* thePrePV = thePrePoint->G << 124 G4VPhysicalVolume* thePrePV = thePrePoint->GetPhysicalVolume();
127 G4VPhysicalVolume* thePostPV = thePostPoint- 125 G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
128 126
129 G4String thePrePVname = " "; << 127 G4String thePrePVname = " ";
130 G4String thePostPVname = " "; 128 G4String thePostPVname = " ";
131 129
132 if (thePostPV) { << 130 if(thePostPV)
133 thePrePVname = thePrePV->GetName(); << 131 {
>> 132 thePrePVname = thePrePV->GetName();
134 thePostPVname = thePostPV->GetName(); 133 thePostPVname = thePostPV->GetName();
135 } 134 }
136 135
137 // Recording data for start 136 // Recording data for start
138 // static const G4ThreeVector ZHat = G4Three 137 // static const G4ThreeVector ZHat = G4ThreeVector(0.0,0.0,1.0);
139 if (theTrack->GetParentID() == 0) { << 138 if(theTrack->GetParentID() == 0)
>> 139 {
140 // This is a primary track 140 // This is a primary track
141 if (theTrack->GetCurrentStepNumber() == 1) << 141 if(theTrack->GetCurrentStepNumber() == 1)
>> 142 {
142 // G4double x = theTrack->GetVer 143 // G4double x = theTrack->GetVertexPosition().x();
143 // G4double y = theTrack->GetVer 144 // G4double y = theTrack->GetVertexPosition().y();
144 // G4double z = theTrack->GetVer 145 // G4double z = theTrack->GetVertexPosition().z();
145 // G4double pz = theTrack->GetVer 146 // G4double pz = theTrack->GetVertexMomentumDirection().z();
146 // G4double fInitTheta = 147 // G4double fInitTheta =
147 // theTrack->Get 148 // theTrack->GetVertexMomentumDirection().angle(ZHat);
148 } 149 }
149 } 150 }
150 151
151 // Retrieve the status of the photon 152 // Retrieve the status of the photon
152 G4OpBoundaryProcessStatus theStatus = Undefi 153 G4OpBoundaryProcessStatus theStatus = Undefined;
153 154
154 static G4ThreadLocal G4ProcessManager* OpMan 155 static G4ThreadLocal G4ProcessManager* OpManager =
155 G4OpticalPhoton::OpticalPhoton()->GetProce 156 G4OpticalPhoton::OpticalPhoton()->GetProcessManager();
156 157
157 if (OpManager) { << 158 if(OpManager)
>> 159 {
158 G4int nproc = OpManager->GetPostStepProces 160 G4int nproc = OpManager->GetPostStepProcessVector()->entries();
159 G4ProcessVector* fPostStepDoItVector = OpM << 161 G4ProcessVector* fPostStepDoItVector =
>> 162 OpManager->GetPostStepProcessVector(typeDoIt);
160 163
161 for (G4int i = 0; i < nproc; ++i) { << 164 for(G4int i = 0; i < nproc; ++i)
>> 165 {
162 G4VProcess* fCurrentProcess = (*fPostSte 166 G4VProcess* fCurrentProcess = (*fPostStepDoItVector)[i];
163 fOpProcess = dynamic_cast<G4OpBoundaryPr 167 fOpProcess = dynamic_cast<G4OpBoundaryProcess*>(fCurrentProcess);
164 if (fOpProcess) { << 168 if(fOpProcess)
>> 169 {
165 theStatus = fOpProcess->GetStatus(); 170 theStatus = fOpProcess->GetStatus();
166 break; 171 break;
167 } 172 }
168 } 173 }
169 } 174 }
170 175
171 // Find the skewness of the ray at first cha 176 // Find the skewness of the ray at first change of boundary
172 if (fInitGamma == -1 << 177 if(fInitGamma == -1 &&
173 && (theStatus == TotalInternalReflection << 178 (theStatus == TotalInternalReflection || theStatus == FresnelReflection ||
174 || theStatus == FresnelRefraction) << 179 theStatus == FresnelRefraction) &&
175 && trackInformation->IsStatus(InsideOfFi << 180 trackInformation->IsStatus(InsideOfFiber))
176 { 181 {
177 G4double px = theTrack->GetVertexMomentumD 182 G4double px = theTrack->GetVertexMomentumDirection().x();
178 G4double py = theTrack->GetVertexMomentumD 183 G4double py = theTrack->GetVertexMomentumDirection().y();
179 G4double x = theTrack->GetPosition().x(); << 184 G4double x = theTrack->GetPosition().x();
180 G4double y = theTrack->GetPosition().y(); << 185 G4double y = theTrack->GetPosition().y();
181 186
182 fInitGamma = x * px + y * py; 187 fInitGamma = x * px + y * py;
183 188
184 fInitGamma = fInitGamma / std::sqrt(px * p << 189 fInitGamma =
>> 190 fInitGamma / std::sqrt(px * px + py * py) / std::sqrt(x * x + y * y);
185 191
186 fInitGamma = std::acos(fInitGamma * rad); 192 fInitGamma = std::acos(fInitGamma * rad);
187 193
188 if (fInitGamma / deg > 90.0) { << 194 if(fInitGamma / deg > 90.0)
>> 195 {
189 fInitGamma = 180 * deg - fInitGamma; 196 fInitGamma = 180 * deg - fInitGamma;
190 } 197 }
191 } 198 }
192 // Record Photons that missed the photon det 199 // Record Photons that missed the photon detector but escaped from readout
193 if (!thePostPV && trackInformation->IsStatus << 200 if(!thePostPV && trackInformation->IsStatus(EscapedFromReadOut))
194 // G4cout << "SteppingAction: status = Esc << 201 {
>> 202 //G4cout << "SteppingAction: status = EscapedFromReadOut" << G4endl;
195 fEventAction->AddEscaped(); 203 fEventAction->AddEscaped();
196 // UpdateHistogramSuccess(thePostPoint,the 204 // UpdateHistogramSuccess(thePostPoint,theTrack);
197 ResetCounters(); 205 ResetCounters();
198 206
199 return; 207 return;
200 } 208 }
201 209
202 // Assumed photons are originated at the fib 210 // Assumed photons are originated at the fiber OR
203 // the fiber is the first material the photo 211 // the fiber is the first material the photon hits
204 switch (theStatus) { << 212 switch(theStatus)
>> 213 {
205 // Exiting the fiber 214 // Exiting the fiber
206 case FresnelRefraction: 215 case FresnelRefraction:
207 case SameMaterial: 216 case SameMaterial:
208 fEventAction->AddExiting(); 217 fEventAction->AddExiting();
209 218
210 if (thePostPVname == "WLSFiber" || thePo << 219 if(thePostPVname == "WLSFiber" || thePostPVname == "Clad1" ||
211 if (trackInformation->IsStatus(Outside << 220 thePostPVname == "Clad2")
>> 221 {
>> 222 if(trackInformation->IsStatus(OutsideOfFiber))
212 trackInformation->AddStatusFlag(Insi 223 trackInformation->AddStatusFlag(InsideOfFiber);
213 224
214 // Set the Exit flag when the photon r 225 // Set the Exit flag when the photon refracted out of the fiber
215 } 226 }
216 else if (trackInformation->IsStatus(Insi << 227 else if(trackInformation->IsStatus(InsideOfFiber))
>> 228 {
217 // EscapedFromReadOut if the z positio 229 // EscapedFromReadOut if the z position is the same as fiber's end
218 if (theTrack->GetPosition().z() == fDe << 230 if(theTrack->GetPosition().z() == fDetector->GetWLSFiberEnd())
>> 231 {
219 trackInformation->AddStatusFlag(Esca 232 trackInformation->AddStatusFlag(EscapedFromReadOut);
220 fCounterEnd++; 233 fCounterEnd++;
221 fEventAction->AddEscapedEnd(); 234 fEventAction->AddEscapedEnd();
222 } 235 }
223 else // Escaped from side 236 else // Escaped from side
224 { 237 {
225 trackInformation->AddStatusFlag(Esca 238 trackInformation->AddStatusFlag(EscapedFromSide);
226 trackInformation->SetExitPosition(th 239 trackInformation->SetExitPosition(theTrack->GetPosition());
227 // UpdateHistogramEscape(thePostPoi 240 // UpdateHistogramEscape(thePostPoint,theTrack);
228 241
229 fCounterMid++; 242 fCounterMid++;
230 fEventAction->AddEscapedMid(); 243 fEventAction->AddEscapedMid();
231 ResetCounters(); 244 ResetCounters();
232 } 245 }
233 246
234 trackInformation->AddStatusFlag(Outsid 247 trackInformation->AddStatusFlag(OutsideOfFiber);
235 trackInformation->SetExitPosition(theT 248 trackInformation->SetExitPosition(theTrack->GetPosition());
236 } 249 }
237 250
238 return; 251 return;
239 252
240 // Internal Reflections 253 // Internal Reflections
241 case TotalInternalReflection: 254 case TotalInternalReflection:
242 255
243 fEventAction->AddTIR(); 256 fEventAction->AddTIR();
244 257
245 // Kill the track if it's number of boun 258 // Kill the track if it's number of bounces exceeded the limit
246 if (fBounceLimit > 0 && fCounterBounce > << 259 if(fBounceLimit > 0 && fCounterBounce >= fBounceLimit)
>> 260 {
247 theTrack->SetTrackStatus(fStopAndKill) 261 theTrack->SetTrackStatus(fStopAndKill);
248 trackInformation->AddStatusFlag(murder 262 trackInformation->AddStatusFlag(murderee);
249 ResetCounters(); 263 ResetCounters();
250 G4cout << "\n Bounce Limit Exceeded" < 264 G4cout << "\n Bounce Limit Exceeded" << G4endl;
251 return; 265 return;
252 } 266 }
253 break; 267 break;
254 268
255 case FresnelReflection: 269 case FresnelReflection:
256 270
257 fCounterBounce++; 271 fCounterBounce++;
258 fEventAction->AddBounce(); 272 fEventAction->AddBounce();
259 273
260 if (thePrePVname == "WLSFiber") { << 274 if(thePrePVname == "WLSFiber")
>> 275 {
261 fCounterWLSBounce++; 276 fCounterWLSBounce++;
262 fEventAction->AddWLSBounce(); 277 fEventAction->AddWLSBounce();
263 } 278 }
264 else if (thePrePVname == "Clad1") { << 279 else if(thePrePVname == "Clad1")
>> 280 {
265 fCounterClad1Bounce++; 281 fCounterClad1Bounce++;
266 fEventAction->AddClad1Bounce(); 282 fEventAction->AddClad1Bounce();
267 } 283 }
268 else if (thePrePVname == "Clad2") { << 284 else if(thePrePVname == "Clad2")
>> 285 {
269 fCounterClad2Bounce++; 286 fCounterClad2Bounce++;
270 fEventAction->AddClad1Bounce(); 287 fEventAction->AddClad1Bounce();
271 } 288 }
272 289
273 // Determine if the photon has reflected 290 // Determine if the photon has reflected off the read-out end
274 if (theTrack->GetPosition().z() == fDete << 291 if(theTrack->GetPosition().z() == fDetector->GetWLSFiberEnd())
275 if (!trackInformation->IsStatus(Reflec << 292 {
276 && trackInformation->IsStatus(Insi << 293 if(!trackInformation->IsStatus(ReflectedAtReadOut) &&
>> 294 trackInformation->IsStatus(InsideOfFiber))
277 { 295 {
278 trackInformation->AddStatusFlag(Refl 296 trackInformation->AddStatusFlag(ReflectedAtReadOut);
279 297
280 if (fDetector->IsPerfectFiber() && t << 298 if(fDetector->IsPerfectFiber() &&
>> 299 theStatus == TotalInternalReflection)
>> 300 {
281 theTrack->SetTrackStatus(fStopAndK 301 theTrack->SetTrackStatus(fStopAndKill);
282 trackInformation->AddStatusFlag(mu 302 trackInformation->AddStatusFlag(murderee);
283 // UpdateHistogramReflect(thePostP 303 // UpdateHistogramReflect(thePostPoint,theTrack);
284 ResetCounters(); 304 ResetCounters();
285 return; 305 return;
286 } 306 }
287 } 307 }
288 } 308 }
289 return; 309 return;
290 310
291 // Reflection off the mirror 311 // Reflection off the mirror
292 case LambertianReflection: 312 case LambertianReflection:
293 case LobeReflection: 313 case LobeReflection:
294 case SpikeReflection: 314 case SpikeReflection:
295 315
296 fEventAction->AddReflected(); 316 fEventAction->AddReflected();
297 // Check if it hits the mirror 317 // Check if it hits the mirror
298 if (thePostPVname == "Mirror") { << 318 if(thePostPVname == "Mirror")
>> 319 {
299 trackInformation->AddStatusFlag(Reflec 320 trackInformation->AddStatusFlag(ReflectedAtMirror);
300 fEventAction->AddMirror(); 321 fEventAction->AddMirror();
301 } 322 }
302 return; 323 return;
303 324
304 // Detected by a detector 325 // Detected by a detector
305 case Detection: 326 case Detection:
306 // Detected automatically with G4OpBound 327 // Detected automatically with G4OpBoundaryProcess->InvokeSD set true
307 328
308 // Stop Tracking when it hits the detect 329 // Stop Tracking when it hits the detector's surface
309 ResetCounters(); 330 ResetCounters();
310 theTrack->SetTrackStatus(fStopAndKill); 331 theTrack->SetTrackStatus(fStopAndKill);
311 return; 332 return;
312 333
313 default: 334 default:
314 break; 335 break;
315 } 336 }
316 337
317 // Check for absorbed photons 338 // Check for absorbed photons
318 if (theTrack->GetTrackStatus() != fAlive && << 339 if(theTrack->GetTrackStatus() != fAlive &&
>> 340 trackInformation->IsStatus(InsideOfFiber))
>> 341 {
319 // UpdateHistogramAbsorb(thePostPoint,theT 342 // UpdateHistogramAbsorb(thePostPoint,theTrack);
320 ResetCounters(); 343 ResetCounters();
321 return; 344 return;
322 } 345 }
323 } 346 }
324 347