OpNovice2/src/SteppingAction.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/OpNovice2/src/SteppingAction.cc
0028 /// \brief Implementation of the SteppingAction class
0029 //
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0031
0032 #include "SteppingAction.hh"
0033
0034 #include "HistoManager.hh"
0035 #include "Run.hh"
0036 #include "SteppingMessenger.hh"
0037 #include "TrackInformation.hh"
0038
0039 #include "G4Cerenkov.hh"
0040 #include "G4Event.hh"
0041 #include "G4EventManager.hh"
0042 #include "G4OpBoundaryProcess.hh"
0043 #include "G4OpticalPhoton.hh"
0044 #include "G4ProcessManager.hh"
0045 #include "G4RunManager.hh"
0046 #include "G4Scintillation.hh"
0047 #include "G4Step.hh"
0048 #include "G4SteppingManager.hh"
0049 #include "G4SystemOfUnits.hh"
0050 #include "G4Track.hh"
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053 SteppingAction::SteppingAction() : G4UserSteppingAction()
0054 {
0055   fSteppingMessenger = new SteppingMessenger(this);
0056 }
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059 SteppingAction::~SteppingAction()
0060 {
0061   delete fSteppingMessenger;
0062 }
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065 void SteppingAction::UserSteppingAction(const G4Step* step)
0066 {
0067   static G4ParticleDefinition* opticalphoton = G4OpticalPhoton::OpticalPhotonDefinition();
0068
0069   G4AnalysisManager* analysisMan = G4AnalysisManager::Instance();
0070   Run* run = static_cast<Run*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
0071
0072   G4Track* track = step->GetTrack();
0073   G4StepPoint* endPoint = step->GetPostStepPoint();
0074   G4StepPoint* startPoint = step->GetPreStepPoint();
0075
0076   const G4DynamicParticle* theParticle = track->GetDynamicParticle();
0077   const G4ParticleDefinition* particleDef = theParticle->GetParticleDefinition();
0078
0079   auto trackInfo = (TrackInformation*)(track->GetUserInformation());
0080
0081   if (particleDef == opticalphoton) {
0082     const G4VProcess* pds = endPoint->GetProcessDefinedStep();
0083     G4String procname = pds->GetProcessName();
0084     if (procname == "OpAbsorption") {
0085       run->AddOpAbsorption();
0086       if (trackInfo->GetIsFirstTankX()) {
0087         run->AddOpAbsorptionPrior();
0088       }
0089     }
0090     else if (procname == "OpRayleigh") {
0091       run->AddRayleigh();
0092     }
0093     else if (procname == "OpWLS") {
0094       G4double en = track->GetKineticEnergy();
0095       run->AddWLSAbsorption();
0096       run->AddWLSAbsorptionEnergy(en);
0097       analysisMan->FillH1(4, en / eV);  // absorption energy
0098       // loop over secondaries, create statistics
0099       // const std::vector<const G4Track*>* secondaries =
0100       auto secondaries = step->GetSecondaryInCurrentStep();
0101       for (auto sec : *secondaries) {
0102         en = sec->GetKineticEnergy();
0103         run->AddWLSEmission();
0104         run->AddWLSEmissionEnergy(en);
0105         analysisMan->FillH1(5, en / eV);  // emission energy
0106         G4double time = sec->GetGlobalTime();
0107         analysisMan->FillH1(6, time / ns);
0108       }
0109     }
0110     else if (procname == "OpWLS2") {
0111       G4double en = track->GetKineticEnergy();
0112       run->AddWLS2Absorption();
0113       run->AddWLS2AbsorptionEnergy(en);
0114       analysisMan->FillH1(7, en / eV);  // absorption energy
0115       // loop over secondaries, create statistics
0116       // const std::vector<const G4Track*>* secondaries =
0117       auto secondaries = step->GetSecondaryInCurrentStep();
0118       for (auto sec : *secondaries) {
0119         en = sec->GetKineticEnergy();
0120         run->AddWLS2Emission();
0121         run->AddWLS2EmissionEnergy(en);
0122         analysisMan->FillH1(8, en / eV);  // emission energy
0123         G4double time = sec->GetGlobalTime();
0124         analysisMan->FillH1(9, time / ns);
0125       }
0126     }
0127
0128     // optical process has endpt on bdry,
0129     if (endPoint->GetStepStatus() == fGeomBoundary) {
0130       G4ThreeVector p0 = startPoint->GetMomentumDirection();
0131       G4ThreeVector p1 = endPoint->GetMomentumDirection();
0132
0133       G4OpBoundaryProcessStatus theStatus = Undefined;
0134
0135       G4ProcessManager* OpManager = opticalphoton->GetProcessManager();
0136       G4ProcessVector* postStepDoItVector = OpManager->GetPostStepProcessVector(typeDoIt);
0137       G4int n_proc = postStepDoItVector->entries();
0138
0139       if (trackInfo->GetIsFirstTankX()) {
0140         G4double px1 = p1.x();
0141         G4double py1 = p1.y();
0142         G4double pz1 = p1.z();
0143         // do not count Absorbed or Detected photons here
0144         if (track->GetTrackStatus() != fStopAndKill) {
0145           if (px1 < 0.) {
0146             analysisMan->FillH1(11, px1);
0147             analysisMan->FillH1(12, py1);
0148             analysisMan->FillH1(13, pz1);
0149           }
0150           else {
0151             analysisMan->FillH1(14, px1);
0152             analysisMan->FillH1(15, py1);
0153             analysisMan->FillH1(16, pz1);
0154           }
0155         }
0156
0157         trackInfo->SetIsFirstTankX(false);
0158         run->AddTotalSurface();
0159
0160         for (G4int i = 0; i < n_proc; ++i) {
0161           G4VProcess* currentProcess = (*postStepDoItVector)[i];
0162
0163           auto opProc = dynamic_cast<G4OpBoundaryProcess*>(currentProcess);
0164           if (opProc) {
0165             G4double angle = std::acos(p0.x());
0166             theStatus = opProc->GetStatus();
0167             analysisMan->FillH1(10, theStatus);
0168             switch (theStatus) {
0169               case Transmission:
0170                 run->AddTransmission();
0171                 analysisMan->FillH1(25, angle / deg);
0172                 break;
0173               case FresnelRefraction:
0174                 run->AddFresnelRefraction();
0175                 analysisMan->FillH1(17, px1);
0176                 analysisMan->FillH1(18, py1);
0177                 analysisMan->FillH1(19, pz1);
0178                 analysisMan->FillH1(20, angle / deg);
0179                 break;
0180               case FresnelReflection:
0181                 run->AddFresnelReflection();
0182                 analysisMan->FillH1(21, angle / deg);
0183                 analysisMan->FillH1(23, angle / deg);
0184                 break;
0185               case TotalInternalReflection:
0186                 run->AddTotalInternalReflection();
0187                 analysisMan->FillH1(22, angle / deg);
0188                 analysisMan->FillH1(23, angle / deg);
0189                 break;
0190               case LambertianReflection:
0191                 run->AddLambertianReflection();
0192                 break;
0193               case LobeReflection:
0194                 run->AddLobeReflection();
0195                 break;
0196               case SpikeReflection:
0197                 run->AddSpikeReflection();
0198                 analysisMan->FillH1(26, angle / deg);
0199                 break;
0200               case BackScattering:
0201                 run->AddBackScattering();
0202                 break;
0203               case Absorption:
0204                 analysisMan->FillH1(24, angle / deg);
0205                 run->AddAbsorption();
0206                 break;
0207               case Detection:
0208                 run->AddDetection();
0209                 break;
0210               case NotAtBoundary:
0211                 run->AddNotAtBoundary();
0212                 break;
0213               case SameMaterial:
0214                 run->AddSameMaterial();
0215                 break;
0216               case StepTooSmall:
0217                 run->AddStepTooSmall();
0218                 break;
0219               case NoRINDEX:
0220                 run->AddNoRINDEX();
0221                 break;
0222               case PolishedLumirrorAirReflection:
0223                 run->AddPolishedLumirrorAirReflection();
0224                 break;
0225               case PolishedLumirrorGlueReflection:
0226                 run->AddPolishedLumirrorGlueReflection();
0227                 break;
0228               case PolishedAirReflection:
0229                 run->AddPolishedAirReflection();
0230                 break;
0231               case PolishedTeflonAirReflection:
0232                 run->AddPolishedTeflonAirReflection();
0233                 break;
0234               case PolishedTiOAirReflection:
0235                 run->AddPolishedTiOAirReflection();
0236                 break;
0237               case PolishedTyvekAirReflection:
0238                 run->AddPolishedTyvekAirReflection();
0239                 break;
0240               case PolishedVM2000AirReflection:
0241                 run->AddPolishedVM2000AirReflection();
0242                 break;
0243               case PolishedVM2000GlueReflection:
0244                 run->AddPolishedVM2000AirReflection();
0245                 break;
0246               case EtchedLumirrorAirReflection:
0247                 run->AddEtchedLumirrorAirReflection();
0248                 break;
0249               case EtchedLumirrorGlueReflection:
0250                 run->AddEtchedLumirrorGlueReflection();
0251                 break;
0252               case EtchedAirReflection:
0253                 run->AddEtchedAirReflection();
0254                 break;
0255               case EtchedTeflonAirReflection:
0256                 run->AddEtchedTeflonAirReflection();
0257                 break;
0258               case EtchedTiOAirReflection:
0259                 run->AddEtchedTiOAirReflection();
0260                 break;
0261               case EtchedTyvekAirReflection:
0262                 run->AddEtchedTyvekAirReflection();
0263                 break;
0264               case EtchedVM2000AirReflection:
0265                 run->AddEtchedVM2000AirReflection();
0266                 break;
0267               case EtchedVM2000GlueReflection:
0268                 run->AddEtchedVM2000AirReflection();
0269                 break;
0270               case GroundLumirrorAirReflection:
0271                 run->AddGroundLumirrorAirReflection();
0272                 break;
0273               case GroundLumirrorGlueReflection:
0274                 run->AddGroundLumirrorGlueReflection();
0275                 break;
0276               case GroundAirReflection:
0277                 run->AddGroundAirReflection();
0278                 break;
0279               case GroundTeflonAirReflection:
0280                 run->AddGroundTeflonAirReflection();
0281                 break;
0282               case GroundTiOAirReflection:
0283                 run->AddGroundTiOAirReflection();
0284                 break;
0285               case GroundTyvekAirReflection:
0286                 run->AddGroundTyvekAirReflection();
0287                 break;
0288               case GroundVM2000AirReflection:
0289                 run->AddGroundVM2000AirReflection();
0290                 break;
0291               case GroundVM2000GlueReflection:
0292                 run->AddGroundVM2000AirReflection();
0293                 break;
0294               case Dichroic:
0295                 run->AddDichroic();
0296                 break;
0297               case CoatedDielectricReflection:
0298                 run->AddCoatedDielectricReflection();
0299                 break;
0300               case CoatedDielectricRefraction:
0301                 run->AddCoatedDielectricRefraction();
0302                 break;
0303               case CoatedDielectricFrustratedTransmission:
0304                 run->AddCoatedDielectricFrustratedTransmission();
0305                 break;
0306               default:
0307                 G4cout << "theStatus: " << theStatus << " was none of the above." << G4endl;
0308                 break;
0309             }
0310           }
0311         }
0312       }
0313       // when studying boundary scattering, it can be useful to only
0314       // visualize the photon before and after the first surface. If
0315       // selected, kill the photon when reaching the second surface
0316       // (note that there are 2 steps at the boundary, so the counter
0317       // equals 0 and 1 on the first surface)
0318       if (fKillOnSecondSurface) {
0319         if (trackInfo->GetReflectionNumber() >= 2) {
0320           track->SetTrackStatus(fStopAndKill);
0321         }
0322       }
0323       trackInfo->IncrementReflectionNumber();
0324     }
0325
0326     // This block serves to test that G4OpBoundaryProcess sets the group
0327     // velocity correctly. It is not necessary to include in user code.
0328     // Only steps where pre- and post- are the same material, to avoid
0329     // incorrect checks (so, in practice, set e.g. OpRayleigh low enough
0330     // for particles to step in the interior of each volume.
0331     if (endPoint->GetMaterial() == startPoint->GetMaterial()) {
0332       G4double trackVelocity = track->GetVelocity();
0333       G4double materialVelocity = CLHEP::c_light;
0334       G4MaterialPropertyVector* velVector =
0335         endPoint->GetMaterial()->GetMaterialPropertiesTable()->GetProperty(kGROUPVEL);
0336       if (velVector) {
0337         materialVelocity = velVector->Value(theParticle->GetTotalMomentum(), fIdxVelocity);
0338       }
0339
0340       if (std::abs(trackVelocity - materialVelocity) > 1e-9 * CLHEP::c_light) {
0341         G4ExceptionDescription ed;
0342         ed << "Optical photon group velocity: " << trackVelocity / (cm / ns)
0343            << " cm/ns is not what is expected from " << G4endl << "the material properties, "
0344            << materialVelocity / (cm / ns) << " cm/ns";
0345         G4Exception("OpNovice2 SteppingAction", "OpNovice2_1", FatalException, ed);
0346       }
0347     }
0348     // end of group velocity test
0349   }
0350
0351   else {  // particle != opticalphoton
0352     // print how many Cerenkov and scint photons produced this step
0353     // this demonstrates use of GetNumPhotons()
0354     auto proc_man = track->GetDynamicParticle()->GetParticleDefinition()->GetProcessManager();
0355     G4ProcessVector* proc_vec = proc_man->GetPostStepProcessVector(typeDoIt);
0356     G4int n_proc = proc_vec->entries();
0357
0358     G4int n_scint = 0;
0359     G4int n_cer = 0;
0360     for (G4int i = 0; i < n_proc; ++i) {
0361       G4String proc_name = (*proc_vec)[i]->GetProcessName();
0362       if (proc_name == "Cerenkov") {
0363         auto cer = (G4Cerenkov*)(*proc_vec)[i];
0364         n_cer = cer->GetNumPhotons();
0365       }
0366       else if (proc_name == "Scintillation") {
0367         auto scint = (G4Scintillation*)(*proc_vec)[i];
0368         n_scint = scint->GetNumPhotons();
0369       }
0370     }
0371     if (fVerbose > 0) {
0372       if (n_cer > 0 || n_scint > 0) {
0373         G4cout << "In this step, " << n_cer << " Cerenkov and " << n_scint
0374                << " scintillation photons were produced." << G4endl;
0375       }
0376     }
0377
0378     // loop over secondaries, create statistics
0379     const std::vector<const G4Track*>* secondaries = step->GetSecondaryInCurrentStep();
0380
0381     for (auto sec : *secondaries) {
0382       if (sec->GetDynamicParticle()->GetParticleDefinition() == opticalphoton) {
0383         G4String creator_process = sec->GetCreatorProcess()->GetProcessName();
0384         if (creator_process == "Cerenkov") {
0385           G4double en = sec->GetKineticEnergy();
0386           run->AddCerenkovEnergy(en);
0387           run->AddCerenkov();
0388           analysisMan->FillH1(1, en / eV);
0389         }
0390         else if (creator_process == "Scintillation") {
0391           G4double en = sec->GetKineticEnergy();
0392           run->AddScintillationEnergy(en);
0393           run->AddScintillation();
0394           analysisMan->FillH1(2, en / eV);
0395
0396           G4double time = sec->GetGlobalTime();
0397           analysisMan->FillH1(3, time / ns);
0398         }
0399       }
0400     }
0401   }
0402
0403   return;
0404 }
0405
0406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......