Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-02-23 09:22:23 

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