Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-02 07:37:24 

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 Run.cc
0027 /// \brief Implementation of the Run class
0028 
0029 #include "Run.hh"
0030 
0031 #include "DetectorConstruction.hh"
0032 #include "HistoManager.hh"
0033 
0034 #include "G4OpBoundaryProcess.hh"
0035 #include "G4SystemOfUnits.hh"
0036 #include "G4UnitsTable.hh"
0037 
0038 #include <numeric>
0039 
0040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0041 Run::Run() : G4Run()
0042 {
0043   fBoundaryProcs.assign(43, 0);
0044 }
0045 
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy, G4bool polarized,
0048                      G4double polarization)
0049 {
0050   fParticle = particle;
0051   fEkin = energy;
0052   fPolarized = polarized;
0053   fPolarization = polarization;
0054 }
0055 
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057 void Run::Merge(const G4Run* run)
0058 {
0059   const Run* localRun = static_cast<const Run*>(run);
0060 
0061   // pass information about primary particle
0062   fParticle = localRun->fParticle;
0063   fEkin = localRun->fEkin;
0064   fPolarized = localRun->fPolarized;
0065   fPolarization = localRun->fPolarization;
0066 
0067   fCerenkovEnergy += localRun->fCerenkovEnergy;
0068   fScintEnergy += localRun->fScintEnergy;
0069   fWLSAbsorptionEnergy += localRun->fWLSAbsorptionEnergy;
0070   fWLSEmissionEnergy += localRun->fWLSEmissionEnergy;
0071   fWLS2AbsorptionEnergy += localRun->fWLS2AbsorptionEnergy;
0072   fWLS2EmissionEnergy += localRun->fWLS2EmissionEnergy;
0073 
0074   fCerenkovCount += localRun->fCerenkovCount;
0075   fScintCount += localRun->fScintCount;
0076   fWLSAbsorptionCount += localRun->fWLSAbsorptionCount;
0077   fWLSEmissionCount += localRun->fWLSEmissionCount;
0078   fWLS2AbsorptionCount += localRun->fWLS2AbsorptionCount;
0079   fWLS2EmissionCount += localRun->fWLS2EmissionCount;
0080   fRayleighCount += localRun->fRayleighCount;
0081 
0082   fTotalSurface += localRun->fTotalSurface;
0083 
0084   fOpAbsorption += localRun->fOpAbsorption;
0085   fOpAbsorptionPrior += localRun->fOpAbsorptionPrior;
0086 
0087   for (size_t i = 0; i < fBoundaryProcs.size(); ++i) {
0088     fBoundaryProcs[i] += localRun->fBoundaryProcs[i];
0089   }
0090 
0091   G4Run::Merge(run);
0092 }
0093 
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0095 void Run::EndOfRun()
0096 {
0097   if (numberOfEvent == 0) return;
0098   auto TotNbofEvents = (G4double)numberOfEvent;
0099 
0100   G4AnalysisManager* analysisMan = G4AnalysisManager::Instance();
0101   G4int id = analysisMan->GetH1Id("Cerenkov spectrum");
0102   analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
0103   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0104 
0105   id = analysisMan->GetH1Id("Scintillation spectrum");
0106   analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
0107   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0108 
0109   id = analysisMan->GetH1Id("Scintillation time");
0110   analysisMan->SetH1XAxisTitle(id, "Creation time [ns]");
0111   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0112 
0113   id = analysisMan->GetH1Id("WLS abs");
0114   analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
0115   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0116 
0117   id = analysisMan->GetH1Id("WLS em");
0118   analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
0119   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0120 
0121   id = analysisMan->GetH1Id("WLS time");
0122   analysisMan->SetH1XAxisTitle(id, "Creation time [ns]");
0123   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0124 
0125   id = analysisMan->GetH1Id("WLS2 abs");
0126   analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
0127   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0128 
0129   id = analysisMan->GetH1Id("WLS2 em");
0130   analysisMan->SetH1XAxisTitle(id, "Energy [eV]");
0131   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0132 
0133   id = analysisMan->GetH1Id("WLS2 time");
0134   analysisMan->SetH1XAxisTitle(id, "Creation time [ns]");
0135   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0136 
0137   id = analysisMan->GetH1Id("bdry status");
0138   analysisMan->SetH1XAxisTitle(id, "Status code");
0139   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0140 
0141   id = analysisMan->GetH1Id("x_backward");
0142   analysisMan->SetH1XAxisTitle(id, "Direction cosine");
0143   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0144 
0145   id = analysisMan->GetH1Id("y_backward");
0146   analysisMan->SetH1XAxisTitle(id, "Direction cosine");
0147   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0148 
0149   id = analysisMan->GetH1Id("z_backward");
0150   analysisMan->SetH1XAxisTitle(id, "Direction cosine");
0151   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0152 
0153   id = analysisMan->GetH1Id("x_forward");
0154   analysisMan->SetH1XAxisTitle(id, "Direction cosine");
0155   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0156 
0157   id = analysisMan->GetH1Id("y_forward");
0158   analysisMan->SetH1XAxisTitle(id, "Direction cosine");
0159   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0160 
0161   id = analysisMan->GetH1Id("z_forward");
0162   analysisMan->SetH1XAxisTitle(id, "Direction cosine");
0163   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0164 
0165   id = analysisMan->GetH1Id("x_fresnel");
0166   analysisMan->SetH1XAxisTitle(id, "Direction cosine");
0167   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0168 
0169   id = analysisMan->GetH1Id("y_fresnel");
0170   analysisMan->SetH1XAxisTitle(id, "Direction cosine");
0171   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0172 
0173   id = analysisMan->GetH1Id("z_fresnel");
0174   analysisMan->SetH1XAxisTitle(id, "Direction cosine");
0175   analysisMan->SetH1YAxisTitle(id, "Number of photons");
0176 
0177   id = analysisMan->GetH1Id("Fresnel reflection");
0178   analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
0179   analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
0180 
0181   id = analysisMan->GetH1Id("Fresnel refraction");
0182   analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
0183   analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
0184 
0185   id = analysisMan->GetH1Id("Total internal reflection");
0186   analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
0187   analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
0188 
0189   id = analysisMan->GetH1Id("Fresnel reflection plus TIR");
0190   analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
0191   analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
0192 
0193   id = analysisMan->GetH1Id("Absorption");
0194   analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
0195   analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
0196 
0197   id = analysisMan->GetH1Id("Transmitted");
0198   analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
0199   analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
0200 
0201   id = analysisMan->GetH1Id("Spike reflection");
0202   analysisMan->SetH1XAxisTitle(id, "Angle [deg]");
0203   analysisMan->SetH1YAxisTitle(id, "Fraction of photons");
0204 
0205   const auto det =
0206     (const DetectorConstruction*)(G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0207 
0208   std::ios::fmtflags mode = G4cout.flags();
0209   G4int prec = G4cout.precision(2);
0210 
0211   G4cout << "\n    Run Summary\n";
0212   G4cout << "---------------------------------\n";
0213   G4cout << "Primary particle was: " << fParticle->GetParticleName() << " with energy "
0214          << G4BestUnit(fEkin, "Energy") << "." << G4endl;
0215   G4cout << "Number of events: " << numberOfEvent << G4endl;
0216 
0217   G4cout << "Material of world: " << det->GetWorldMaterial()->GetName() << G4endl;
0218   G4cout << "Material of tank:  " << det->GetTankMaterial()->GetName() << G4endl << G4endl;
0219 
0220   if (fParticle->GetParticleName() != "opticalphoton") {
0221     G4cout << "Average energy of Cerenkov photons created per event: "
0222            << (fCerenkovEnergy / eV) / TotNbofEvents << " eV." << G4endl;
0223     G4cout << "Average number of Cerenkov photons created per event: "
0224            << fCerenkovCount / TotNbofEvents << G4endl;
0225     if (fCerenkovCount > 0) {
0226       G4cout << " Average energy per photon: " << (fCerenkovEnergy / eV) / fCerenkovCount << " eV."
0227              << G4endl;
0228     }
0229     G4cout << "Average energy of scintillation photons created per event: "
0230            << (fScintEnergy / eV) / TotNbofEvents << " eV." << G4endl;
0231     G4cout << "Average number of scintillation photons created per event: "
0232            << fScintCount / TotNbofEvents << G4endl;
0233     if (fScintCount > 0) {
0234       G4cout << " Average energy per photon: " << (fScintEnergy / eV) / fScintCount << " eV."
0235              << G4endl;
0236     }
0237   }
0238 
0239   G4cout << "Average number of photons absorbed by WLS per event: "
0240          << fWLSAbsorptionCount / G4double(TotNbofEvents) << " " << G4endl;
0241   if (fWLSAbsorptionCount > 0) {
0242     G4cout << " Average energy per photon: " << (fWLSAbsorptionEnergy / eV) / fWLSAbsorptionCount
0243            << " eV." << G4endl;
0244   }
0245   G4cout << "Average number of photons created by WLS per event: "
0246          << fWLSEmissionCount / TotNbofEvents << G4endl;
0247   if (fWLSEmissionCount > 0) {
0248     G4cout << " Average energy per photon: " << (fWLSEmissionEnergy / eV) / fWLSEmissionCount
0249            << " eV." << G4endl;
0250   }
0251   G4cout << "Average energy of WLS photons created per event: "
0252          << (fWLSEmissionEnergy / eV) / TotNbofEvents << " eV." << G4endl;
0253 
0254   G4cout << "Average number of photons absorbed by WLS2 per event: "
0255          << fWLS2AbsorptionCount / G4double(TotNbofEvents) << " " << G4endl;
0256   if (fWLS2AbsorptionCount > 0) {
0257     G4cout << " Average energy per photon: " << (fWLS2AbsorptionEnergy / eV) / fWLS2AbsorptionCount
0258            << " eV." << G4endl;
0259   }
0260   G4cout << "Average number of photons created by WLS2 per event: "
0261          << fWLS2EmissionCount / TotNbofEvents << G4endl;
0262   if (fWLS2EmissionCount > 0) {
0263     G4cout << " Average energy per photon: " << (fWLS2EmissionEnergy / eV) / fWLS2EmissionCount
0264            << " eV." << G4endl;
0265   }
0266   G4cout << "Average energy of WLS2 photons created per event: "
0267          << (fWLS2EmissionEnergy / eV) / TotNbofEvents << " eV." << G4endl;
0268 
0269   G4cout << "Average number of OpRayleigh per event:   " << fRayleighCount / TotNbofEvents
0270          << G4endl;
0271   G4cout << "Average number of OpAbsorption per event: " << fOpAbsorption / TotNbofEvents << G4endl;
0272   G4cout << "\nSurface events (on +X surface, maximum one per photon) this run:" << G4endl;
0273   G4cout << "# of primary particles:      " << std::setw(8) << TotNbofEvents << G4endl;
0274   G4cout << "OpAbsorption before surface: " << std::setw(8) << fOpAbsorptionPrior << G4endl;
0275   G4cout << "Total # of surface events:   " << std::setw(8) << fTotalSurface << G4endl;
0276   if (fParticle->GetParticleName() == "opticalphoton") {
0277     G4cout << "Unaccounted for:             " << std::setw(8)
0278            << fTotalSurface + fOpAbsorptionPrior - TotNbofEvents << G4endl;
0279   }
0280   G4cout << "\nSurface events by process:" << G4endl;
0281   if (fBoundaryProcs[Transmission] > 0) {
0282     G4cout << "  Transmission:              " << std::setw(8) << fBoundaryProcs[Transmission]
0283            << G4endl;
0284   }
0285   if (fBoundaryProcs[FresnelRefraction] > 0) {
0286     G4cout << "  Fresnel refraction:        " << std::setw(8) << fBoundaryProcs[FresnelRefraction]
0287            << G4endl;
0288   }
0289   if (fBoundaryProcs[FresnelReflection] > 0) {
0290     G4cout << "  Fresnel reflection:        " << std::setw(8) << fBoundaryProcs[FresnelReflection]
0291            << G4endl;
0292   }
0293   if (fBoundaryProcs[TotalInternalReflection] > 0) {
0294     G4cout << "  Total internal reflection: " << std::setw(8)
0295            << fBoundaryProcs[TotalInternalReflection] << G4endl;
0296   }
0297   if (fBoundaryProcs[LambertianReflection] > 0) {
0298     G4cout << "  Lambertian reflection:     " << std::setw(8)
0299            << fBoundaryProcs[LambertianReflection] << G4endl;
0300   }
0301   if (fBoundaryProcs[LobeReflection] > 0) {
0302     G4cout << "  Lobe reflection:           " << std::setw(8) << fBoundaryProcs[LobeReflection]
0303            << G4endl;
0304   }
0305   if (fBoundaryProcs[SpikeReflection] > 0) {
0306     G4cout << "  Spike reflection:          " << std::setw(8) << fBoundaryProcs[SpikeReflection]
0307            << G4endl;
0308   }
0309   if (fBoundaryProcs[BackScattering] > 0) {
0310     G4cout << "  Backscattering:            " << std::setw(8) << fBoundaryProcs[BackScattering]
0311            << G4endl;
0312   }
0313   if (fBoundaryProcs[Absorption] > 0) {
0314     G4cout << "  Absorption:                " << std::setw(8) << fBoundaryProcs[Absorption]
0315            << G4endl;
0316   }
0317   if (fBoundaryProcs[Detection] > 0) {
0318     G4cout << "  Detection:                 " << std::setw(8) << fBoundaryProcs[Detection]
0319            << G4endl;
0320   }
0321   if (fBoundaryProcs[NotAtBoundary] > 0) {
0322     G4cout << "  Not at boundary:           " << std::setw(8) << fBoundaryProcs[NotAtBoundary]
0323            << G4endl;
0324   }
0325   if (fBoundaryProcs[SameMaterial] > 0) {
0326     G4cout << "  Same material:             " << std::setw(8) << fBoundaryProcs[SameMaterial]
0327            << G4endl;
0328   }
0329   if (fBoundaryProcs[StepTooSmall] > 0) {
0330     G4cout << "  Step too small:            " << std::setw(8) << fBoundaryProcs[StepTooSmall]
0331            << G4endl;
0332   }
0333   if (fBoundaryProcs[NoRINDEX] > 0) {
0334     G4cout << "  No RINDEX:                 " << std::setw(8) << fBoundaryProcs[NoRINDEX] << G4endl;
0335   }
0336   // LBNL polished
0337   if (fBoundaryProcs[PolishedLumirrorAirReflection] > 0) {
0338     G4cout << "  Polished Lumirror Air reflection: " << std::setw(8)
0339            << fBoundaryProcs[PolishedLumirrorAirReflection] << G4endl;
0340   }
0341   if (fBoundaryProcs[PolishedLumirrorGlueReflection] > 0) {
0342     G4cout << "  Polished Lumirror Glue reflection: " << std::setw(8)
0343            << fBoundaryProcs[PolishedLumirrorGlueReflection] << G4endl;
0344   }
0345   if (fBoundaryProcs[PolishedAirReflection] > 0) {
0346     G4cout << "  Polished Air reflection: " << std::setw(8) << fBoundaryProcs[PolishedAirReflection]
0347            << G4endl;
0348   }
0349   if (fBoundaryProcs[PolishedTeflonAirReflection] > 0) {
0350     G4cout << "  Polished Teflon Air reflection: " << std::setw(8)
0351            << fBoundaryProcs[PolishedTeflonAirReflection] << G4endl;
0352   }
0353   if (fBoundaryProcs[PolishedTiOAirReflection] > 0) {
0354     G4cout << "  Polished TiO Air reflection: " << std::setw(8)
0355            << fBoundaryProcs[PolishedTiOAirReflection] << G4endl;
0356   }
0357   if (fBoundaryProcs[PolishedTyvekAirReflection] > 0) {
0358     G4cout << "  Polished Tyvek Air reflection: " << std::setw(8)
0359            << fBoundaryProcs[PolishedTyvekAirReflection] << G4endl;
0360   }
0361   if (fBoundaryProcs[PolishedVM2000AirReflection] > 0) {
0362     G4cout << "  Polished VM2000 Air reflection: " << std::setw(8)
0363            << fBoundaryProcs[PolishedVM2000AirReflection] << G4endl;
0364   }
0365   if (fBoundaryProcs[PolishedVM2000GlueReflection] > 0) {
0366     G4cout << "  Polished VM2000 Glue reflection: " << std::setw(8)
0367            << fBoundaryProcs[PolishedVM2000GlueReflection] << G4endl;
0368   }
0369   // LBNL etched
0370   if (fBoundaryProcs[EtchedLumirrorAirReflection] > 0) {
0371     G4cout << "  Etched Lumirror Air reflection: " << std::setw(8)
0372            << fBoundaryProcs[EtchedLumirrorAirReflection] << G4endl;
0373   }
0374   if (fBoundaryProcs[EtchedLumirrorGlueReflection] > 0) {
0375     G4cout << "  Etched Lumirror Glue reflection: " << std::setw(8)
0376            << fBoundaryProcs[EtchedLumirrorGlueReflection] << G4endl;
0377   }
0378   if (fBoundaryProcs[EtchedAirReflection] > 0) {
0379     G4cout << "  Etched Air reflection: " << std::setw(8) << fBoundaryProcs[EtchedAirReflection]
0380            << G4endl;
0381   }
0382   if (fBoundaryProcs[EtchedTeflonAirReflection] > 0) {
0383     G4cout << "  Etched Teflon Air reflection: " << std::setw(8)
0384            << fBoundaryProcs[EtchedTeflonAirReflection] << G4endl;
0385   }
0386   if (fBoundaryProcs[EtchedTiOAirReflection] > 0) {
0387     G4cout << "  Etched TiO Air reflection: " << std::setw(8)
0388            << fBoundaryProcs[EtchedTiOAirReflection] << G4endl;
0389   }
0390   if (fBoundaryProcs[EtchedTyvekAirReflection] > 0) {
0391     G4cout << "  Etched Tyvek Air reflection: " << std::setw(8)
0392            << fBoundaryProcs[EtchedTyvekAirReflection] << G4endl;
0393   }
0394   if (fBoundaryProcs[EtchedVM2000AirReflection] > 0) {
0395     G4cout << "  Etched VM2000 Air reflection: " << std::setw(8)
0396            << fBoundaryProcs[EtchedVM2000AirReflection] << G4endl;
0397   }
0398   if (fBoundaryProcs[EtchedVM2000GlueReflection] > 0) {
0399     G4cout << "  Etched VM2000 Glue reflection: " << std::setw(8)
0400            << fBoundaryProcs[EtchedVM2000GlueReflection] << G4endl;
0401   }
0402   // LBNL ground
0403   if (fBoundaryProcs[GroundLumirrorAirReflection] > 0) {
0404     G4cout << "  Ground Lumirror Air reflection: " << std::setw(8)
0405            << fBoundaryProcs[GroundLumirrorAirReflection] << G4endl;
0406   }
0407   if (fBoundaryProcs[GroundLumirrorGlueReflection] > 0) {
0408     G4cout << "  Ground Lumirror Glue reflection: " << std::setw(8)
0409            << fBoundaryProcs[GroundLumirrorGlueReflection] << G4endl;
0410   }
0411   if (fBoundaryProcs[GroundAirReflection] > 0) {
0412     G4cout << "  Ground Air reflection: " << std::setw(8) << fBoundaryProcs[GroundAirReflection]
0413            << G4endl;
0414   }
0415   if (fBoundaryProcs[GroundTeflonAirReflection] > 0) {
0416     G4cout << "  Ground Teflon Air reflection: " << std::setw(8)
0417            << fBoundaryProcs[GroundTeflonAirReflection] << G4endl;
0418   }
0419   if (fBoundaryProcs[GroundTiOAirReflection] > 0) {
0420     G4cout << "  Ground TiO Air reflection: " << std::setw(8)
0421            << fBoundaryProcs[GroundTiOAirReflection] << G4endl;
0422   }
0423   if (fBoundaryProcs[GroundTyvekAirReflection] > 0) {
0424     G4cout << "  Ground Tyvek Air reflection: " << std::setw(8)
0425            << fBoundaryProcs[GroundTyvekAirReflection] << G4endl;
0426   }
0427   if (fBoundaryProcs[GroundVM2000AirReflection] > 0) {
0428     G4cout << "  Ground VM2000 Air reflection: " << std::setw(8)
0429            << fBoundaryProcs[GroundVM2000AirReflection] << G4endl;
0430   }
0431   if (fBoundaryProcs[GroundVM2000GlueReflection] > 0) {
0432     G4cout << "  Ground VM2000 Glue reflection: " << std::setw(8)
0433            << fBoundaryProcs[GroundVM2000GlueReflection] << G4endl;
0434   }
0435   if (fBoundaryProcs[CoatedDielectricRefraction] > 0) {
0436     G4cout << "  CoatedDielectricRefraction: " << std::setw(8)
0437            << fBoundaryProcs[CoatedDielectricRefraction] << G4endl;
0438   }
0439   if (fBoundaryProcs[CoatedDielectricReflection] > 0) {
0440     G4cout << "  CoatedDielectricReflection: " << std::setw(8)
0441            << fBoundaryProcs[CoatedDielectricReflection] << G4endl;
0442   }
0443   if (fBoundaryProcs[CoatedDielectricFrustratedTransmission] > 0) {
0444     G4cout << "  CoatedDielectricFrustratedTransmission: " << std::setw(8)
0445            << fBoundaryProcs[CoatedDielectricFrustratedTransmission] << G4endl;
0446   }
0447 
0448   G4int sum = std::accumulate(fBoundaryProcs.begin(), fBoundaryProcs.end(), 0);
0449   G4cout << " Sum:                        " << std::setw(8) << sum << G4endl;
0450   G4cout << " Unaccounted for:            " << std::setw(8) << fTotalSurface - sum << G4endl;
0451 
0452   G4cout << "---------------------------------\n";
0453   G4cout.setf(mode, std::ios::floatfield);
0454   G4cout.precision(prec);
0455 
0456   G4int histo_id_refract = analysisMan->GetH1Id("Fresnel refraction");
0457   G4int histo_id_reflect = analysisMan->GetH1Id("Fresnel reflection plus TIR");
0458   G4int histo_id_spike = analysisMan->GetH1Id("Spike reflection");
0459   G4int histo_id_absorption = analysisMan->GetH1Id("Absorption");
0460 
0461   if (analysisMan->GetH1Activation(histo_id_refract)
0462       && analysisMan->GetH1Activation(histo_id_reflect))
0463   {
0464     G4double rindex1 =
0465       det->GetTankMaterial()->GetMaterialPropertiesTable()->GetProperty(kRINDEX)->Value(fEkin);
0466     G4double rindex2 =
0467       det->GetWorldMaterial()->GetMaterialPropertiesTable()->GetProperty(kRINDEX)->Value(fEkin);
0468 
0469     auto histo_refract = analysisMan->GetH1(histo_id_refract);
0470     auto histo_reflect = analysisMan->GetH1(histo_id_reflect);
0471     // std::vector<G4double> refract;
0472     std::vector<G4double> reflect;
0473     // std::vector<G4double> tir;
0474     std::vector<G4double> tot;
0475     for (size_t i = 0; i < histo_refract->axis().bins(); ++i) {
0476       // refract.push_back(histo_refract->bin_height(i));
0477       reflect.push_back(histo_reflect->bin_height(i));
0478       // tir.push_back(histo_TIR->bin_height(i));
0479       tot.push_back(histo_refract->bin_height(i) + histo_reflect->bin_height(i));
0480     }
0481 
0482     // find Brewster angle: Rp = 0
0483     //  need enough statistics for this method to work
0484     G4double min_angle = -1.;
0485     G4double min_val = DBL_MAX;
0486     G4double bin_width = 0.;
0487     for (size_t i = 0; i < reflect.size(); ++i) {
0488       if (reflect[i] < min_val) {
0489         min_val = reflect[i];
0490         min_angle = histo_reflect->axis().bin_lower_edge(i);
0491         bin_width =
0492           histo_reflect->axis().bin_upper_edge(i) - histo_reflect->axis().bin_lower_edge(i);
0493         min_angle += bin_width / 2.;
0494       }
0495     }
0496     G4cout << "Polarization of primary optical photons: " << fPolarization / deg << " deg."
0497            << G4endl;
0498     if (fPolarized && fPolarization == 0.0) {
0499       G4cout << "Reflectance shows a minimum at: " << min_angle << " +/- " << bin_width / 2;
0500       G4cout << " deg. Expected Brewster angle: "
0501              << (360. / CLHEP::twopi) * std::atan(rindex2 / rindex1) << " deg. " << G4endl;
0502     }
0503 
0504     // find angle of total internal reflection:  T -> 0
0505     //   last bin for T > 0
0506     min_angle = -1.;
0507     min_val = DBL_MAX;
0508     for (size_t i = 0; i < histo_refract->axis().bins() - 1; ++i) {
0509       if (histo_refract->bin_height(i) > 0. && histo_refract->bin_height(i + 1) == 0.) {
0510         min_angle = histo_refract->axis().bin_lower_edge(i);
0511         bin_width =
0512           histo_reflect->axis().bin_upper_edge(i) - histo_reflect->axis().bin_lower_edge(i);
0513         min_angle += bin_width / 2.;
0514         break;
0515       }
0516     }
0517     if (fPolarized) {
0518       G4cout << "Fresnel transmission goes to 0 at: " << min_angle << " +/- " << bin_width / 2.
0519              << " deg."
0520              << " Expected: " << (360. / CLHEP::twopi) * std::asin(rindex2 / rindex1) << " deg."
0521              << G4endl;
0522     }
0523 
0524     // Normalize the transmission/reflection histos so that max is 1.
0525     // Only if x values are the same
0526     if ((analysisMan->GetH1Nbins(histo_id_refract) == analysisMan->GetH1Nbins(histo_id_reflect))
0527         && (analysisMan->GetH1Xmin(histo_id_refract) == analysisMan->GetH1Xmin(histo_id_reflect))
0528         && (analysisMan->GetH1Xmax(histo_id_refract) == analysisMan->GetH1Xmax(histo_id_reflect)))
0529     {
0530       unsigned int ent;
0531       G4double sw;
0532       G4double sw2;
0533       G4double sx2;
0534       G4double sx2w;
0535       for (size_t bin = 0; bin < histo_refract->axis().bins(); ++bin) {
0536         // "bin+1" below because bin 0 is underflow bin
0537         // NB. We are ignoring underflow/overflow bins
0538         histo_refract->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0539         if (tot[bin] > 0) {
0540           sw /= tot[bin];
0541           // bin error is sqrt(sw2)
0542           sw2 /= (tot[bin] * tot[bin]);
0543           sx2 /= (tot[bin] * tot[bin]);
0544           sx2w /= (tot[bin] * tot[bin]);
0545           histo_refract->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0546         }
0547 
0548         histo_reflect->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0549         if (tot[bin] > 0) {
0550           sw /= tot[bin];
0551           // bin error is sqrt(sw2)
0552           sw2 /= (tot[bin] * tot[bin]);
0553           sx2 /= (tot[bin] * tot[bin]);
0554           sx2w /= (tot[bin] * tot[bin]);
0555           histo_reflect->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0556         }
0557 
0558         G4int histo_id_fresnelrefl = analysisMan->GetH1Id("Fresnel reflection");
0559         auto histo_fresnelreflect = analysisMan->GetH1(histo_id_fresnelrefl);
0560         histo_fresnelreflect->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0561         if (tot[bin] > 0) {
0562           sw /= tot[bin];
0563           // bin error is sqrt(sw2)
0564           sw2 /= (tot[bin] * tot[bin]);
0565           sx2 /= (tot[bin] * tot[bin]);
0566           sx2w /= (tot[bin] * tot[bin]);
0567           histo_fresnelreflect->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0568         }
0569 
0570         G4int histo_id_TIR = analysisMan->GetH1Id("Total internal reflection");
0571         auto histo_TIR = analysisMan->GetH1(histo_id_TIR);
0572         if (analysisMan->GetH1Activation(histo_id_TIR)) {
0573           histo_TIR->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0574           if (tot[bin] > 0) {
0575             sw /= tot[bin];
0576             // bin error is sqrt(sw2)
0577             sw2 /= (tot[bin] * tot[bin]);
0578             sx2 /= (tot[bin] * tot[bin]);
0579             sx2w /= (tot[bin] * tot[bin]);
0580             histo_TIR->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0581           }
0582         }
0583       }
0584     }
0585     else {
0586       G4cout << "Not going to normalize refraction and reflection "
0587              << "histograms because bins are not the same." << G4endl;
0588     }
0589   }
0590 
0591   // complex index of refraction; have spike reflection and absorption
0592   // Only works for polished surfaces. Ground surfaces neglected.
0593   else if (analysisMan->GetH1Activation(histo_id_absorption)
0594            && analysisMan->GetH1Activation(histo_id_spike))
0595   {
0596     auto histo_spike = analysisMan->GetH1(histo_id_spike);
0597     auto histo_absorption = analysisMan->GetH1(histo_id_absorption);
0598 
0599     std::vector<G4double> tot;
0600     for (size_t i = 0; i < histo_absorption->axis().bins(); ++i) {
0601       tot.push_back(histo_absorption->bin_height(i) + histo_spike->bin_height(i));
0602     }
0603 
0604     if ((analysisMan->GetH1Nbins(histo_id_absorption) == analysisMan->GetH1Nbins(histo_id_spike))
0605         && (analysisMan->GetH1Xmin(histo_id_absorption) == analysisMan->GetH1Xmin(histo_id_spike))
0606         && (analysisMan->GetH1Xmax(histo_id_absorption) == analysisMan->GetH1Xmax(histo_id_spike)))
0607     {
0608       unsigned int ent;
0609       G4double sw;
0610       G4double sw2;
0611       G4double sx2;
0612       G4double sx2w;
0613       for (size_t bin = 0; bin < histo_absorption->axis().bins(); ++bin) {
0614         // "bin+1" below because bin 0 is underflow bin
0615         // NB. We are ignoring underflow/overflow bins
0616         histo_absorption->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0617         if (tot[bin] > 0) {
0618           sw /= tot[bin];
0619           // bin error is sqrt(sw2)
0620           sw2 /= (tot[bin] * tot[bin]);
0621           sx2 /= (tot[bin] * tot[bin]);
0622           sx2w /= (tot[bin] * tot[bin]);
0623           histo_absorption->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0624         }
0625 
0626         histo_spike->get_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0627         if (tot[bin] > 0) {
0628           sw /= tot[bin];
0629           // bin error is sqrt(sw2)
0630           sw2 /= (tot[bin] * tot[bin]);
0631           sx2 /= (tot[bin] * tot[bin]);
0632           sx2w /= (tot[bin] * tot[bin]);
0633           histo_spike->set_bin_content(bin + 1, ent, sw, sw2, sx2, sx2w);
0634         }
0635       }
0636     }
0637     else {
0638       G4cout << "Not going to normalize spike reflection and absorption "
0639              << "histograms because bins are not the same." << G4endl;
0640     }
0641   }
0642 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!