ReverseMC01/src/RMC01AnalysisManager.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 biasing/ReverseMC01/src/RMC01AnalysisManager.cc
0027 /// \brief Implementation of the RMC01AnalysisManager class
0028 //
0029 //
0030 //////////////////////////////////////////////////////////////
0031 //      Class Name:        RMC01AnalysisManager
0032 //        Author:               L. Desorgher
0033 //         Organisation:         SpaceIT GmbH
0034 //        Contract:        ESA contract 21435/08/NL/AT
0035 //         Customer:             ESA/ESTEC
0036 //////////////////////////////////////////////////////////////
0037
0038 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0039 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0040
0041 #include "RMC01AnalysisManager.hh"
0042
0043 #include "RMC01AnalysisManagerMessenger.hh"
0044 #include "RMC01SD.hh"
0045
0046 #include "G4AdjointSimManager.hh"
0047 #include "G4Electron.hh"
0048 #include "G4Gamma.hh"
0049 #include "G4PhysicalConstants.hh"
0050 #include "G4Proton.hh"
0051 #include "G4RunManager.hh"
0052 #include "G4SDManager.hh"
0053 #include "G4SystemOfUnits.hh"
0054 #include "G4THitsCollection.hh"
0055 #include "G4Timer.hh"
0056
0057 RMC01AnalysisManager* RMC01AnalysisManager::fInstance = 0;
0058
0059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0060
0061 RMC01AnalysisManager::RMC01AnalysisManager()
0062   : fAccumulated_edep(0.),
0063     fAccumulated_edep2(0.),
0064     fMean_edep(0.),
0065     fError_mean_edep(0.),
0066     fRelative_error(0.),
0067     fElapsed_time(0.),
0068     fPrecision_to_reach(0.),
0069     fStop_run_if_precision_reached(true),
0070     fNb_evt_modulo_for_convergence_test(5000),
0071     fEdep_rmatrix_vs_electron_prim_energy(0),
0072     fElectron_current_rmatrix_vs_electron_prim_energy(0),
0073     fGamma_current_rmatrix_vs_electron_prim_energy(0),
0074     fEdep_rmatrix_vs_gamma_prim_energy(0),
0075     fElectron_current_rmatrix_vs_gamma_prim_energy(0),
0076     fGamma_current_rmatrix_vs_gamma_prim_energy(0),
0077     fEdep_rmatrix_vs_proton_prim_energy(0),
0078     fElectron_current_rmatrix_vs_proton_prim_energy(0),
0079     fProton_current_rmatrix_vs_proton_prim_energy(0),
0080     fGamma_current_rmatrix_vs_proton_prim_energy(0),
0081     fFactoryOn(false),
0082     fPrimSpectrumType(EXPO),
0083     fAlpha_or_E0(.5 * MeV),
0084     fAmplitude_prim_spectrum(1.),
0085     fEmin_prim_spectrum(1. * keV),
0086     fEmax_prim_spectrum(20. * MeV),
0087     fAdjoint_sim_mode(true),
0088     fNb_evt_per_adj_evt(2)
0089 {
0090   fMsg = new RMC01AnalysisManagerMessenger(this);
0091
0092   //-------------
0093   // Timer for convergence vector
0094   //-------------
0095
0096   fTimer = new G4Timer();
0097
0098   //---------------------------------
0099   // Primary particle ID for normalisation of adjoint results
0100   //---------------------------------
0101
0102   fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
0103
0104   fFileName[0] = "sim";
0105 }
0106
0107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0108
0109 RMC01AnalysisManager::~RMC01AnalysisManager()
0110 {
0111   ;
0112 }
0113
0114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0115
0116 RMC01AnalysisManager* RMC01AnalysisManager::GetInstance()
0117 {
0118   if (fInstance == 0) fInstance = new RMC01AnalysisManager;
0119   return fInstance;
0120 }
0121
0122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0123
0124 void RMC01AnalysisManager::BeginOfRun(const G4Run* aRun)
0125 {
0126   fAccumulated_edep = 0.;
0127   fAccumulated_edep2 = 0.;
0128   fNentry = 0.0;
0129   fRelative_error = 1.;
0130   fMean_edep = 0.;
0131   fError_mean_edep = 0.;
0132   fAdjoint_sim_mode = G4AdjointSimManager::GetInstance()->GetAdjointSimMode();
0133
0134   if (fAdjoint_sim_mode) {
0135     fNb_evt_per_adj_evt = aRun->GetNumberOfEventToBeProcessed()
0136                           / G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
0137     fConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt", std::ios::out);
0138     fConvergenceFileOutput << "Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]" << std::endl;
0139   }
0140   else {
0141     fConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt", std::ios::out);
0142     fConvergenceFileOutput << "Edep per event [MeV]\terror[MeV]\tcomputing_time[s]" << std::endl;
0143   }
0144   fConvergenceFileOutput.setf(std::ios::scientific);
0145   fConvergenceFileOutput.precision(6);
0146
0147   fTimer->Start();
0148   fElapsed_time = 0.;
0149
0150   Book();
0151 }
0152
0153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0154
0155 void RMC01AnalysisManager::EndOfRun(const G4Run* aRun)
0156 {
0157   fTimer->Stop();
0158   G4int nb_evt = aRun->GetNumberOfEvent();
0159   G4double factor = 1. / nb_evt;
0160   if (!fAdjoint_sim_mode) {
0161     G4cout << "Results of forward simulation!" << std::endl;
0162     G4cout << "edep per event [MeV] = " << fMean_edep << std::endl;
0163     G4cout << "error[MeV] = " << fError_mean_edep << std::endl;
0164   }
0165
0166   else {
0167     G4cout << "Results of reverse/adjoint simulation!" << std::endl;
0168     G4cout << "normalised edep [MeV] = " << fMean_edep << std::endl;
0169     G4cout << "error[MeV] = " << fError_mean_edep << std::endl;
0170     factor = 1. * G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun() * fNb_evt_per_adj_evt
0171              / aRun->GetNumberOfEvent();
0172   }
0173   Save(factor);
0174   fConvergenceFileOutput.close();
0175 }
0176
0177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0178
0179 void RMC01AnalysisManager::BeginOfEvent(const G4Event*)
0180 {
0181   ;
0182 }
0183
0184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0185
0186 void RMC01AnalysisManager::EndOfEvent(const G4Event* anEvent)
0187 {
0188   if (fAdjoint_sim_mode)
0189     EndOfEventForAdjointSimulation(anEvent);
0190   else
0191     EndOfEventForForwardSimulation(anEvent);
0192
0193   // Test convergence. The error is already computed
0194   //--------------------------------------
0195   G4int nb_event = anEvent->GetEventID() + 1;
0196   // G4double factor=1.;
0197   if (fAdjoint_sim_mode) {
0198     G4double n_adj_evt = nb_event / fNb_evt_per_adj_evt;
0199     // nb_event/fNb_evt_per_adj_evt;
0200     if (n_adj_evt * fNb_evt_per_adj_evt == nb_event) {
0201       nb_event = static_cast<G4int>(n_adj_evt);
0202     }
0203     else
0204       nb_event = 0;
0205   }
0206
0207   if (nb_event > 100 && fStop_run_if_precision_reached && fPrecision_to_reach > fRelative_error) {
0208     G4cout << fPrecision_to_reach * 100. << "%  Precision reached!" << std::endl;
0209     fTimer->Stop();
0210     fElapsed_time += fTimer->GetRealElapsed();
0211     fConvergenceFileOutput << fMean_edep << '\t' << fError_mean_edep << '\t' << fElapsed_time
0212                            << std::endl;
0213     G4RunManager::GetRunManager()->AbortRun(true);
0214   }
0215
0216   if (nb_event > 0 && nb_event % fNb_evt_modulo_for_convergence_test == 0) {
0217     fTimer->Stop();
0218     fElapsed_time += fTimer->GetRealElapsed();
0219     fTimer->Start();
0220     fConvergenceFileOutput << fMean_edep << '\t' << fError_mean_edep << '\t' << fElapsed_time
0221                            << std::endl;
0222   }
0223 }
0224
0225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0226
0227 void RMC01AnalysisManager::EndOfEventForForwardSimulation(const G4Event* anEvent)
0228 {
0229   G4SDManager* SDman = G4SDManager::GetSDMpointer();
0230   G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
0231   RMC01DoubleWithWeightHitsCollection* edepCollection =
0232     (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("edep")));
0233
0234   RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
0235     (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_electron")));
0236
0237   RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
0238     (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_proton")));
0239
0240   RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
0241     (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_gamma")));
0242
0243   // Total energy deposited in Event
0244   //-------------------------------
0245   G4double totEdep = 0;
0246   std::size_t i;
0247   for (i = 0; i < edepCollection->entries(); ++i)
0248     totEdep += (*edepCollection)[i]->GetValue() * (*edepCollection)[i]->GetWeight();
0249
0250   if (totEdep > 0.) {
0251     fAccumulated_edep += totEdep;
0252     fAccumulated_edep2 += totEdep * totEdep;
0253     fNentry += 1.0;
0254     G4PrimaryParticle* thePrimary = anEvent->GetPrimaryVertex()->GetPrimary();
0255     G4double E0 = thePrimary->GetG4code()->GetPDGMass();
0256     G4double P = thePrimary->GetMomentum().mag();
0257     G4double prim_ekin = std::sqrt(E0 * E0 + P * P) - E0;
0258     fEdep_vs_prim_ekin->fill(prim_ekin, totEdep);
0259   }
0260   ComputeMeanEdepAndError(anEvent, fMean_edep, fError_mean_edep);
0261   if (fError_mean_edep > 0) fRelative_error = fError_mean_edep / fMean_edep;
0262
0263   // Particle current on sensitive cylinder
0264   //-------------------------------------
0265
0266   for (i = 0; i < electronCurrentCollection->entries(); ++i) {
0267     G4double ekin = (*electronCurrentCollection)[i]->GetValue();
0268     G4double weight = (*electronCurrentCollection)[i]->GetWeight();
0269     fElectron_current->fill(ekin, weight);
0270   }
0271
0272   for (i = 0; i < protonCurrentCollection->entries(); ++i) {
0273     G4double ekin = (*protonCurrentCollection)[i]->GetValue();
0274     G4double weight = (*protonCurrentCollection)[i]->GetWeight();
0275     fProton_current->fill(ekin, weight);
0276   }
0277
0278   for (i = 0; i < gammaCurrentCollection->entries(); ++i) {
0279     G4double ekin = (*gammaCurrentCollection)[i]->GetValue();
0280     G4double weight = (*gammaCurrentCollection)[i]->GetWeight();
0281     fGamma_current->fill(ekin, weight);
0282   }
0283 }
0284
0285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0286
0287 void RMC01AnalysisManager::EndOfEventForAdjointSimulation(const G4Event* anEvent)
0288 {
0289   // Output from Sensitive volume computed during the forward tracking phase
0290   //-----------------------------------------------------------------------
0291   G4SDManager* SDman = G4SDManager::GetSDMpointer();
0292   G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
0293   RMC01DoubleWithWeightHitsCollection* edepCollection =
0294     (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("edep")));
0295
0296   RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
0297     (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_electron")));
0298
0299   RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
0300     (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_proton")));
0301
0302   RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
0303     (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_gamma")));
0304
0305   // Computation of total energy deposited in fwd tracking phase
0306   //-------------------------------
0307   G4double totEdep = 0;
0308   std::size_t i;
0309   for (i = 0; i < edepCollection->entries(); ++i)
0310     totEdep += (*edepCollection)[i]->GetValue() * (*edepCollection)[i]->GetWeight();
0311
0312   // Output from adjoint tracking phase
0313   //----------------------------------------------------------------------------
0314
0315   G4AdjointSimManager* theAdjointSimManager = G4AdjointSimManager::GetInstance();
0316
0317   size_t nb_adj_track = theAdjointSimManager->GetNbOfAdointTracksReachingTheExternalSurface();
0318   G4double total_normalised_weight = 0.;
0319
0320   // We need to loop over the adjoint tracks that have reached the external
0321   // surface.
0322   for (std::size_t j = 0; j < nb_adj_track; ++j) {
0323     G4int pdg_nb = theAdjointSimManager->GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack(j);
0324     G4double prim_ekin = theAdjointSimManager->GetEkinAtEndOfLastAdjointTrack(j);
0325     G4double adj_weight = theAdjointSimManager->GetWeightAtEndOfLastAdjointTrack(j);
0326
0327     // Factor of normalisation to user defined prim spectrum (power law or exp)
0328     //------------------------------------------------------------------------
0329     G4double normalised_weight = 0.;
0330     if (pdg_nb == fPrimPDG_ID && prim_ekin >= fEmin_prim_spectrum
0331         && prim_ekin <= fEmax_prim_spectrum)
0332       normalised_weight = adj_weight * PrimDiffAndDirFluxForAdjointSim(prim_ekin);
0333     total_normalised_weight += normalised_weight;
0334
0335     // Answer matrices
0336     //-------------
0337     G4H1* edep_rmatrix = 0;
0338     G4H2* electron_current_rmatrix = 0;
0339     G4H2* gamma_current_rmatrix = 0;
0340     G4H2* proton_current_rmatrix = 0;
0341
0342     if (pdg_nb == G4Electron::Electron()->GetPDGEncoding()) {  // e- matrices
0343       edep_rmatrix = fEdep_rmatrix_vs_electron_prim_energy;
0344       electron_current_rmatrix = fElectron_current_rmatrix_vs_electron_prim_energy;
0345       gamma_current_rmatrix = fGamma_current_rmatrix_vs_electron_prim_energy;
0346     }
0347     else if (pdg_nb == G4Gamma::Gamma()->GetPDGEncoding()) {
0348       // gammma answer matrices
0349       edep_rmatrix = fEdep_rmatrix_vs_gamma_prim_energy;
0350       electron_current_rmatrix = fElectron_current_rmatrix_vs_gamma_prim_energy;
0351       gamma_current_rmatrix = fGamma_current_rmatrix_vs_gamma_prim_energy;
0352     }
0353     else if (pdg_nb == G4Proton::Proton()->GetPDGEncoding()) {
0354       // proton answer matrices
0355       edep_rmatrix = fEdep_rmatrix_vs_proton_prim_energy;
0356       electron_current_rmatrix = fElectron_current_rmatrix_vs_proton_prim_energy;
0357       gamma_current_rmatrix = fGamma_current_rmatrix_vs_proton_prim_energy;
0358       proton_current_rmatrix = fProton_current_rmatrix_vs_proton_prim_energy;
0359     }
0360     // Register histo edep vs prim ekin
0361     //----------------------------------
0362     if (normalised_weight > 0) fEdep_vs_prim_ekin->fill(prim_ekin, totEdep * normalised_weight);
0363     // Registering answer matrix
0364     //---------------------------
0365     edep_rmatrix->fill(prim_ekin, totEdep * adj_weight / cm2);
0366
0367     // Registering of current of particles on the sensitive volume
0368     //------------------------------------------------------------
0369
0370     for (i = 0; i < electronCurrentCollection->entries(); ++i) {
0371       G4double ekin = (*electronCurrentCollection)[i]->GetValue();
0372       G4double weight = (*electronCurrentCollection)[i]->GetWeight();
0373       fElectron_current->fill(ekin, weight * normalised_weight);
0374       electron_current_rmatrix->fill(prim_ekin, ekin, weight * adj_weight / cm2);
0375     }
0376     for (i = 0; i < protonCurrentCollection->entries(); ++i) {
0377       G4double ekin = (*protonCurrentCollection)[i]->GetValue();
0378       G4double weight = (*protonCurrentCollection)[i]->GetWeight();
0379       fProton_current->fill(ekin, weight * normalised_weight);
0380       proton_current_rmatrix->fill(prim_ekin, ekin, weight * adj_weight / cm2);
0381     }
0382     for (i = 0; i < gammaCurrentCollection->entries(); ++i) {
0383       G4double ekin = (*gammaCurrentCollection)[i]->GetValue();
0384       G4double weight = (*gammaCurrentCollection)[i]->GetWeight();
0385       fGamma_current->fill(ekin, weight * normalised_weight);
0386       gamma_current_rmatrix->fill(prim_ekin, ekin, weight * adj_weight / cm2);
0387     }
0388   }
0389
0390   // Registering of total energy deposited in Event
0391   //-------------------------------
0392   G4bool new_mean_computed = false;
0393   if (totEdep > 0.) {
0394     if (total_normalised_weight > 0.) {
0395       G4double edep = totEdep * total_normalised_weight;
0396
0397       // Check if the edep is not wrongly too high
0398       //-----------------------------------------
0399       G4double new_mean(0.0), new_error(0.0);
0400       fAccumulated_edep += edep;
0401       fAccumulated_edep2 += edep * edep;
0402       fNentry += 1.0;
0403       ComputeMeanEdepAndError(anEvent, new_mean, new_error);
0404       G4double new_relative_error = 1.;
0405       if (new_error > 0) new_relative_error = new_error / new_mean;
0406       if (fRelative_error < 0.10 && new_relative_error > 1.5 * fRelative_error) {
0407         G4cout << "Potential wrong adjoint weight!" << std::endl;
0408         G4cout << "The results of this event will not be registered!" << std::endl;
0409         G4cout << "previous mean edep [MeV] " << fMean_edep << std::endl;
0410         G4cout << "previous relative error " << fRelative_error << std::endl;
0411         G4cout << "new rejected mean edep [MeV] " << new_mean << std::endl;
0412         G4cout << "new rejected relative error " << new_relative_error << std::endl;
0413         fAccumulated_edep -= edep;
0414         fAccumulated_edep2 -= edep * edep;
0415         fNentry -= 1.0;
0416         return;
0417       }
0418       else {  // accepted
0419         fMean_edep = new_mean;
0420         fError_mean_edep = new_error;
0421         fRelative_error = new_relative_error;
0422         new_mean_computed = true;
0423       }
0424     }
0425
0426     if (!new_mean_computed) {
0427       ComputeMeanEdepAndError(anEvent, fMean_edep, fError_mean_edep);
0428       fRelative_error = (fMean_edep > 0.0) ? fError_mean_edep / fMean_edep : 0.0;
0429     }
0430   }
0431 }
0432
0433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0434
0435 G4double RMC01AnalysisManager::PrimDiffAndDirFluxForAdjointSim(G4double prim_energy)
0436 {
0437   G4double flux = fAmplitude_prim_spectrum;
0438   if (fPrimSpectrumType == EXPO)
0439     flux *= std::exp(-prim_energy / fAlpha_or_E0);
0440   else
0441     flux *= std::pow(prim_energy, -fAlpha_or_E0);
0442   return flux;
0443 }
0444
0445 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0446 /*
0447 void  RMC01AnalysisManager::WriteHisto(G4H1* anHisto,
0448             G4double scaling_factor, G4String fileName, G4String header_lines)
0449 { std::fstream FileOutput(fileName, std::ios::out);
0450   FileOutput<<header_lines;
0451   FileOutput.setf(std::ios::scientific);
0452   FileOutput.precision(6);
0453
0454   for (G4int i =0;i<G4int(anHisto->axis().bins());++i) {
0455         FileOutput<<anHisto->axis().bin_lower_edge(i)
0456               <<'\t'<<anHisto->axis().bin_upper_edge(i)
0457               <<'\t'<<anHisto->bin_height(i)*scaling_factor
0458               <<'\t'<<anHisto->bin_error(i)*scaling_factor<<std::endl;
0459   }
0460 }
0461
0462 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0463
0464 void  RMC01AnalysisManager::WriteHisto(G4H2* anHisto,
0465             G4double scaling_factor, G4String fileName, G4String header_lines)
0466 { std::fstream FileOutput(fileName, std::ios::out);
0467   FileOutput<<header_lines;
0468
0469   FileOutput.setf(std::ios::scientific);
0470   FileOutput.precision(6);
0471
0472   for (G4int i =0;i<G4int(anHisto->axis_x().bins());++i) {
0473     for (G4int j =0;j<G4int(anHisto->axis_y().bins());++j) {
0474        FileOutput<<anHisto->axis_x().bin_lower_edge(i)
0475                      <<'\t'<<anHisto->axis_x().bin_upper_edge(i)
0476                    <<'\t'<<anHisto->axis_y().bin_lower_edge(i)
0477                    <<'\t'<<anHisto->axis_y().bin_upper_edge(i)
0478                    <<'\t'<<anHisto->bin_height(i,j)*scaling_factor
0479                 <<'\t'<<anHisto->bin_error(i,j)*scaling_factor
0480                                                                   <<std::endl;
0481         }
0482   }
0483 }
0484 */
0485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0486
0487 void RMC01AnalysisManager::ComputeMeanEdepAndError(const G4Event* anEvent, G4double& mean,
0488                                                    G4double& error)
0489 {
0490   G4int nb_event = anEvent->GetEventID() + 1;
0491   G4double factor = 1.;
0492   if (fAdjoint_sim_mode) {
0493     nb_event /= fNb_evt_per_adj_evt;
0494     factor = 1. * G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
0495   }
0496
0497   // VI: error computation now is based on number of entries and not
0498   //     number of events
0499   // LD: This is wrong! With the use of fNentry the results were no longer
0500   //     correctly normalised. The mean and the error should be computed
0501   //     with nb_event. The old computation has been reset.
0502   // VI: OK, but let computations be double
0503   if (nb_event > 0) {
0504     G4double norm = 1.0 / (G4double)nb_event;
0505     mean = fAccumulated_edep * norm;
0506     G4double mean_x2 = fAccumulated_edep2 * norm;
0507     G4double zz = mean_x2 - mean * mean;
0508     /*
0509      G4cout << "Nevt= " << nb_event <<  " mean= " << mean
0510             << "  mean_x2= " <<  mean_x2 << " x2 - x*x= "
0511             << zz << G4endl;
0512     */
0513     error = factor * std::sqrt(std::max(zz, 0.) * norm);
0514     mean *= factor;
0515   }
0516   else {
0517     mean = 0;
0518     error = 0;
0519   }
0520   // G4cout << "Aend: " << mean << " " << error << G4endl;
0521 }
0522
0523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0524
0525 void RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim(const G4String& particle_name,
0526                                                               G4double omni_fluence, G4double E0,
0527                                                               G4double Emin, G4double Emax)
0528 {
0529   fPrimSpectrumType = EXPO;
0530   if (particle_name == "e-")
0531     fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
0532   else if (particle_name == "gamma")
0533     fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEncoding();
0534   else if (particle_name == "proton")
0535     fPrimPDG_ID = G4Proton::Proton()->GetPDGEncoding();
0536   else {
0537     G4cout << "The particle that you did select is not in the candidate "
0538            << "list for primary [e-, gamma, proton]!" << G4endl;
0539     return;
0540   }
0541   fAlpha_or_E0 = E0;
0542   fAmplitude_prim_spectrum =
0543     omni_fluence / E0 / (std::exp(-Emin / E0) - std::exp(-Emax / E0)) / 4. / pi;
0544   fEmin_prim_spectrum = Emin;
0545   fEmax_prim_spectrum = Emax;
0546 }
0547
0548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0549
0550 void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim(const G4String& particle_name,
0551                                                                    G4double omni_fluence,
0552                                                                    G4double alpha, G4double Emin,
0553                                                                    G4double Emax)
0554 {
0555   fPrimSpectrumType = POWER;
0556   if (particle_name == "e-")
0557     fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
0558   else if (particle_name == "gamma")
0559     fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEncoding();
0560   else if (particle_name == "proton")
0561     fPrimPDG_ID = G4Proton::Proton()->GetPDGEncoding();
0562   else {
0563     G4cout << "The particle that you did select is not in the candidate"
0564            << " list for primary [e-, gamma, proton]!" << G4endl;
0565     return;
0566   }
0567
0568   if (alpha == 1.) {
0569     fAmplitude_prim_spectrum = omni_fluence / std::log(Emax / Emin) / 4. / pi;
0570   }
0571   else {
0572     G4double p = 1. - alpha;
0573     fAmplitude_prim_spectrum = omni_fluence / p / (std::pow(Emax, p) - std::pow(Emin, p)) / 4. / pi;
0574   }
0575
0576   fAlpha_or_E0 = alpha;
0577   fEmin_prim_spectrum = Emin;
0578   fEmax_prim_spectrum = Emax;
0579 }
0580
0581 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0582
0583 void RMC01AnalysisManager::Book()
0584 {
0585   //----------------------
0586   // Creation of histograms
0587   //----------------------
0588
0589   // Energy binning of the histograms : 60 log bins over [1keV-1GeV]
0590
0591   G4double emin = 1. * keV;
0592   G4double emax = 1. * GeV;
0593
0594   // file_name
0595   fFileName[0] = "forward_sim";
0596   if (fAdjoint_sim_mode) fFileName[0] = "adjoint_sim";
0597
0598   // Histo manager
0599   G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
0600   theHistoManager->SetDefaultFileType("root");
0601   G4String extension = theHistoManager->GetFileType();
0602   fFileName[1] = fFileName[0] + "." + extension;
0603   theHistoManager->SetFirstHistoId(1);
0604
0605   G4bool fileOpen = theHistoManager->OpenFile(fFileName[0]);
0606   if (!fileOpen) {
0607     G4cout << "\n---> RMC01AnalysisManager::Book(): cannot open " << fFileName[1] << G4endl;
0608     return;
0609   }
0610
0611   // Create directories
0612   theHistoManager->SetHistoDirectoryName("histo");
0613
0614   // Histograms for :
0615   //         1)the forward simulation results
0616   //         2)the Reverse MC  simulation results normalised to a user spectrum
0617   //------------------------------------------------------------------------
0618
0619   G4int idHisto =
0620     theHistoManager->CreateH1(G4String("Edep_vs_prim_ekin"), G4String("edep vs e- primary energy"),
0621                               60, emin, emax, "none", "none", G4String("log"));
0622   fEdep_vs_prim_ekin = theHistoManager->GetH1(idHisto);
0623
0624   idHisto = theHistoManager->CreateH1(G4String("elecron_current"), G4String("electron"), 60, emin,
0625                                       emax, "none", "none", G4String("log"));
0626
0627   fElectron_current = theHistoManager->GetH1(idHisto);
0628
0629   idHisto = theHistoManager->CreateH1(G4String("proton_current"), G4String("proton"), 60, emin,
0630                                       emax, "none", "none", G4String("log"));
0631   fProton_current = theHistoManager->GetH1(idHisto);
0632
0633   idHisto = theHistoManager->CreateH1(G4String("gamma_current"), G4String("gamma"), 60, emin, emax,
0634                                       "none", "none", G4String("log"));
0635   fGamma_current = theHistoManager->GetH1(idHisto);
0636
0637   // Response matrices for the adjoint simulation only
0638   //-----------------------------------------------
0639   if (fAdjoint_sim_mode) {
0640     // Response matrices for external isotropic e- source
0641     //--------------------------------------------------
0642
0643     idHisto = theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_electron_prim_energy"),
0644                                         G4String("electron RM vs e- primary energy"), 60, emin,
0645                                         emax, "none", "none", G4String("log"));
0646     fEdep_rmatrix_vs_electron_prim_energy = theHistoManager->GetH1(idHisto);
0647
0648     idHisto = theHistoManager->CreateH2(
0649       G4String("Electron_current_rmatrix_vs_electron_prim_energy"),
0650       G4String("electron current  RM vs e- primary energy"), 60, emin, emax, 60, emin, emax, "none",
0651       "none", "none", "none", G4String("log"), G4String("log"));
0652
0653     fElectron_current_rmatrix_vs_electron_prim_energy = theHistoManager->GetH2(idHisto);
0654
0655     idHisto = theHistoManager->CreateH2(G4String("Gamma_current_rmatrix_vs_electron_prim_energy"),
0656                                         G4String("gamma current  RM vs e- primary energy"), 60,
0657                                         emin, emax, 60, emin, emax, "none", "none", "none", "none",
0658                                         G4String("log"), G4String("log"));
0659
0660     fGamma_current_rmatrix_vs_electron_prim_energy = theHistoManager->GetH2(idHisto);
0661
0662     // Response matrices for external isotropic gamma source
0663
0664     idHisto = theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_gamma_prim_energy"),
0665                                         G4String("electron RM vs gamma primary energy"), 60, emin,
0666                                         emax, "none", "none", G4String("log"));
0667     fEdep_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH1(idHisto);
0668
0669     idHisto = theHistoManager->CreateH2(G4String("Electron_current_rmatrix_vs_gamma_prim_energy"),
0670                                         G4String("electron current  RM vs gamma primary energy"),
0671                                         60, emin, emax, 60, emin, emax, "none", "none", "none",
0672                                         "none", G4String("log"), G4String("log"));
0673
0674     fElectron_current_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH2(idHisto);
0675
0676     idHisto = theHistoManager->CreateH2(G4String("Gamma_current_rmatrix_vs_gamma_prim_energy"),
0677                                         G4String("gamma current  RM vs gamma primary energy"), 60,
0678                                         emin, emax, 60, emin, emax, "none", "none", "none", "none",
0679                                         G4String("log"), G4String("log"));
0680
0681     fGamma_current_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH2(idHisto);
0682
0683     // Response matrices for external isotropic proton source
0684     idHisto = theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_proton_prim_energy"),
0685                                         G4String("electron RM vs proton primary energy"), 60, emin,
0686                                         emax, "none", "none", G4String("log"));
0687     fEdep_rmatrix_vs_proton_prim_energy = theHistoManager->GetH1(idHisto);
0688
0689     idHisto = theHistoManager->CreateH2(G4String("Electron_current_rmatrix_vs_proton_prim_energy"),
0690                                         G4String("electron current  RM vs proton primary energy"),
0691                                         60, emin, emax, 60, emin, emax, "none", "none", "none",
0692                                         "none", G4String("log"), G4String("log"));
0693
0694     fElectron_current_rmatrix_vs_proton_prim_energy = theHistoManager->GetH2(idHisto);
0695
0696     idHisto = theHistoManager->CreateH2(G4String("Gamma_current_rmatrix_vs_proton_prim_energy"),
0697                                         G4String("gamma current  RM vs proton primary energy"), 60,
0698                                         emin, emax, 60, emin, emax, "none", "none", "none", "none",
0699                                         G4String("log"), G4String("log"));
0700
0701     fGamma_current_rmatrix_vs_proton_prim_energy = theHistoManager->GetH2(idHisto);
0702
0703     idHisto = theHistoManager->CreateH2(G4String("Proton_current_rmatrix_vs_proton_prim_energy"),
0704                                         G4String("proton current  RM vs proton primary energy"), 60,
0705                                         emin, emax, 60, emin, emax, "none", "none", "none", "none",
0706                                         G4String("log"), G4String("log"));
0707
0708     fProton_current_rmatrix_vs_proton_prim_energy = theHistoManager->GetH2(idHisto);
0709   }
0710   fFactoryOn = true;
0711   G4cout << "\n----> Histogram Tree is opened in " << fFileName[1] << G4endl;
0712 }
0713
0714 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0715
0716 void RMC01AnalysisManager::Save(G4double scaling_factor)
0717 {
0718   if (fFactoryOn) {
0719     G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
0720     // scaling of results
0721     //-----------------
0722
0723     for (G4int ind = 1; ind <= theHistoManager->GetNofH1s(); ++ind) {
0724       theHistoManager->SetH1Ascii(ind, true);
0725       theHistoManager->ScaleH1(ind, scaling_factor);
0726     }
0727     for (G4int ind = 1; ind <= theHistoManager->GetNofH2s(); ++ind)
0728       theHistoManager->ScaleH2(ind, scaling_factor);
0729
0730     theHistoManager->Write();
0731     theHistoManager->CloseFile();
0732     G4cout << "\n----> Histogram Tree is saved in " << fFileName[1] << G4endl;
0733
0734     theHistoManager->Clear();
0735     fFactoryOn = false;
0736   }
0737 }