wls/src/WLSSteppingAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file optical/wls/src/WLSSteppingAction.cc
0028 /// \brief Implementation of the WLSSteppingAction class
0029 //
0030 //
0031 #include "WLSSteppingAction.hh"
0032
0033 #include "WLSDetectorConstruction.hh"
0034 #include "WLSEventAction.hh"
0035 #include "WLSPhotonDetSD.hh"
0036 #include "WLSSteppingActionMessenger.hh"
0037 #include "WLSUserTrackInformation.hh"
0038
0039 #include "G4OpBoundaryProcess.hh"
0040 #include "G4OpticalPhoton.hh"
0041 #include "G4ProcessManager.hh"
0042 #include "G4Run.hh"
0043 #include "G4SDManager.hh"
0044 #include "G4Step.hh"
0045 #include "G4StepPoint.hh"
0046 #include "G4SystemOfUnits.hh"
0047 #include "G4ThreeVector.hh"
0048 #include "G4Track.hh"
0049 #include "G4TrackStatus.hh"
0050 #include "G4UImanager.hh"
0051 #include "G4VPhysicalVolume.hh"
0052 #include "G4ios.hh"
0053
0054 // Purpose: Save relevant information into User Track Information
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 WLSSteppingAction::WLSSteppingAction(WLSDetectorConstruction* detector, WLSEventAction* event)
0059   : fDetector(detector), fEventAction(event)
0060 {
0061   fSteppingMessenger = new WLSSteppingActionMessenger(this);
0062   ResetCounters();
0063 }
0064
0065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0066
0067 WLSSteppingAction::~WLSSteppingAction()
0068 {
0069   delete fSteppingMessenger;
0070 }
0071
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073
0074 void WLSSteppingAction::SetBounceLimit(G4int i)
0075 {
0076   fBounceLimit = i;
0077 }
0078
0079 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0080
0081 G4int WLSSteppingAction::GetNumberOfBounces()
0082 {
0083   return fCounterBounce;
0084 }
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087
0088 G4int WLSSteppingAction::GetNumberOfClad1Bounces()
0089 {
0090   return fCounterClad1Bounce;
0091 }
0092
0093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0094
0095 G4int WLSSteppingAction::GetNumberOfClad2Bounces()
0096 {
0097   return fCounterClad2Bounce;
0098 }
0099
0100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0101
0102 G4int WLSSteppingAction::GetNumberOfWLSBounces()
0103 {
0104   return fCounterWLSBounce;
0105 }
0106
0107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0108
0109 G4int WLSSteppingAction::ResetSuccessCounter()
0110 {
0111   G4int temp = fCounterEnd;
0112   fCounterEnd = 0;
0113   return temp;
0114 }
0115
0116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0117
0118 void WLSSteppingAction::UserSteppingAction(const G4Step* theStep)
0119 {
0120   G4Track* theTrack = theStep->GetTrack();
0121   auto trackInformation = (WLSUserTrackInformation*)theTrack->GetUserInformation();
0122
0123   G4StepPoint* thePrePoint = theStep->GetPreStepPoint();
0124   G4StepPoint* thePostPoint = theStep->GetPostStepPoint();
0125
0126   G4VPhysicalVolume* thePrePV = thePrePoint->GetPhysicalVolume();
0127   G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
0128
0129   G4String thePrePVname = " ";
0130   G4String thePostPVname = " ";
0131
0132   if (thePostPV) {
0133     thePrePVname = thePrePV->GetName();
0134     thePostPVname = thePostPV->GetName();
0135   }
0136
0137   // Recording data for start
0138   // static const G4ThreeVector ZHat = G4ThreeVector(0.0,0.0,1.0);
0139   if (theTrack->GetParentID() == 0) {
0140     // This is a primary track
0141     if (theTrack->GetCurrentStepNumber() == 1) {
0142       //        G4double x  = theTrack->GetVertexPosition().x();
0143       //        G4double y  = theTrack->GetVertexPosition().y();
0144       //        G4double z  = theTrack->GetVertexPosition().z();
0145       //        G4double pz = theTrack->GetVertexMomentumDirection().z();
0146       //        G4double fInitTheta =
0147       //                         theTrack->GetVertexMomentumDirection().angle(ZHat);
0148     }
0149   }
0150
0151   // Retrieve the status of the photon
0152   G4OpBoundaryProcessStatus theStatus = Undefined;
0153
0154   static G4ThreadLocal G4ProcessManager* OpManager =
0155     G4OpticalPhoton::OpticalPhoton()->GetProcessManager();
0156
0157   if (OpManager) {
0158     G4int nproc = OpManager->GetPostStepProcessVector()->entries();
0159     G4ProcessVector* fPostStepDoItVector = OpManager->GetPostStepProcessVector(typeDoIt);
0160
0161     for (G4int i = 0; i < nproc; ++i) {
0162       G4VProcess* fCurrentProcess = (*fPostStepDoItVector)[i];
0163       fOpProcess = dynamic_cast<G4OpBoundaryProcess*>(fCurrentProcess);
0164       if (fOpProcess) {
0165         theStatus = fOpProcess->GetStatus();
0166         break;
0167       }
0168     }
0169   }
0170
0171   // Find the skewness of the ray at first change of boundary
0172   if (fInitGamma == -1
0173       && (theStatus == TotalInternalReflection || theStatus == FresnelReflection
0174           || theStatus == FresnelRefraction)
0175       && trackInformation->IsStatus(InsideOfFiber))
0176   {
0177     G4double px = theTrack->GetVertexMomentumDirection().x();
0178     G4double py = theTrack->GetVertexMomentumDirection().y();
0179     G4double x = theTrack->GetPosition().x();
0180     G4double y = theTrack->GetPosition().y();
0181
0182     fInitGamma = x * px + y * py;
0183
0184     fInitGamma = fInitGamma / std::sqrt(px * px + py * py) / std::sqrt(x * x + y * y);
0185
0186     fInitGamma = std::acos(fInitGamma * rad);
0187
0188     if (fInitGamma / deg > 90.0) {
0189       fInitGamma = 180 * deg - fInitGamma;
0190     }
0191   }
0192   // Record Photons that missed the photon detector but escaped from readout
0193   if (!thePostPV && trackInformation->IsStatus(EscapedFromReadOut)) {
0194     // G4cout << "SteppingAction: status = EscapedFromReadOut" << G4endl;
0195     fEventAction->AddEscaped();
0196     // UpdateHistogramSuccess(thePostPoint,theTrack);
0197     ResetCounters();
0198
0199     return;
0200   }
0201
0202   // Assumed photons are originated at the fiber OR
0203   // the fiber is the first material the photon hits
0204   switch (theStatus) {
0205     // Exiting the fiber
0206     case FresnelRefraction:
0207     case SameMaterial:
0208       fEventAction->AddExiting();
0209
0210       if (thePostPVname == "WLSFiber" || thePostPVname == "Clad1" || thePostPVname == "Clad2") {
0211         if (trackInformation->IsStatus(OutsideOfFiber))
0212           trackInformation->AddStatusFlag(InsideOfFiber);
0213
0214         // Set the Exit flag when the photon refracted out of the fiber
0215       }
0216       else if (trackInformation->IsStatus(InsideOfFiber)) {
0217         // EscapedFromReadOut if the z position is the same as fiber's end
0218         if (theTrack->GetPosition().z() == fDetector->GetWLSFiberEnd()) {
0219           trackInformation->AddStatusFlag(EscapedFromReadOut);
0220           fCounterEnd++;
0221           fEventAction->AddEscapedEnd();
0222         }
0223         else  // Escaped from side
0224         {
0225           trackInformation->AddStatusFlag(EscapedFromSide);
0226           trackInformation->SetExitPosition(theTrack->GetPosition());
0227           //  UpdateHistogramEscape(thePostPoint,theTrack);
0228
0229           fCounterMid++;
0230           fEventAction->AddEscapedMid();
0231           ResetCounters();
0232         }
0233
0234         trackInformation->AddStatusFlag(OutsideOfFiber);
0235         trackInformation->SetExitPosition(theTrack->GetPosition());
0236       }
0237
0238       return;
0239
0240     // Internal Reflections
0241     case TotalInternalReflection:
0242
0243       fEventAction->AddTIR();
0244
0245       // Kill the track if it's number of bounces exceeded the limit
0246       if (fBounceLimit > 0 && fCounterBounce >= fBounceLimit) {
0247         theTrack->SetTrackStatus(fStopAndKill);
0248         trackInformation->AddStatusFlag(murderee);
0249         ResetCounters();
0250         G4cout << "\n Bounce Limit Exceeded" << G4endl;
0251         return;
0252       }
0253       break;
0254
0255     case FresnelReflection:
0256
0257       fCounterBounce++;
0258       fEventAction->AddBounce();
0259
0260       if (thePrePVname == "WLSFiber") {
0261         fCounterWLSBounce++;
0262         fEventAction->AddWLSBounce();
0263       }
0264       else if (thePrePVname == "Clad1") {
0265         fCounterClad1Bounce++;
0266         fEventAction->AddClad1Bounce();
0267       }
0268       else if (thePrePVname == "Clad2") {
0269         fCounterClad2Bounce++;
0270         fEventAction->AddClad1Bounce();
0271       }
0272
0273       // Determine if the photon has reflected off the read-out end
0274       if (theTrack->GetPosition().z() == fDetector->GetWLSFiberEnd()) {
0275         if (!trackInformation->IsStatus(ReflectedAtReadOut)
0276             && trackInformation->IsStatus(InsideOfFiber))
0277         {
0278           trackInformation->AddStatusFlag(ReflectedAtReadOut);
0279
0280           if (fDetector->IsPerfectFiber() && theStatus == TotalInternalReflection) {
0281             theTrack->SetTrackStatus(fStopAndKill);
0282             trackInformation->AddStatusFlag(murderee);
0283             // UpdateHistogramReflect(thePostPoint,theTrack);
0284             ResetCounters();
0285             return;
0286           }
0287         }
0288       }
0289       return;
0290
0291     // Reflection off the mirror
0292     case LambertianReflection:
0293     case LobeReflection:
0294     case SpikeReflection:
0295
0296       fEventAction->AddReflected();
0297       // Check if it hits the mirror
0298       if (thePostPVname == "Mirror") {
0299         trackInformation->AddStatusFlag(ReflectedAtMirror);
0300         fEventAction->AddMirror();
0301       }
0302       return;
0303
0304     // Detected by a detector
0305     case Detection:
0306       // Detected automatically with G4OpBoundaryProcess->InvokeSD set true
0307
0308       // Stop Tracking when it hits the detector's surface
0309       ResetCounters();
0310       theTrack->SetTrackStatus(fStopAndKill);
0311       return;
0312
0313     default:
0314       break;
0315   }
0316
0317   // Check for absorbed photons
0318   if (theTrack->GetTrackStatus() != fAlive && trackInformation->IsStatus(InsideOfFiber)) {
0319     // UpdateHistogramAbsorb(thePostPoint,theTrack);
0320     ResetCounters();
0321     return;
0322   }
0323 }