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