Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/biasing/ReverseMC01/src/RMC01AnalysisManager.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/biasing/ReverseMC01/src/RMC01AnalysisManager.cc (Version 11.3.0) and /examples/extended/biasing/ReverseMC01/src/RMC01AnalysisManager.cc (Version 9.4.p4)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 /// \file biasing/ReverseMC01/src/RMC01Analysi <<  26 // $Id: RMC01AnalysisManager.cc,v 1.7 2010-11-11 14:39:42 ldesorgh Exp $
 27 /// \brief Implementation of the RMC01Analysis <<  27 // GEANT4 tag $Name: not supported by cvs2svn $
 28 //                                             << 
 29 //                                                 28 //
 30 //////////////////////////////////////////////     29 //////////////////////////////////////////////////////////////
 31 //      Class Name:        RMC01AnalysisManage <<  30 //      Class Name: RMC01AnalysisManager
 32 //        Author:               L. Desorgher   <<  31 //  Author:         L. Desorgher
 33 //         Organisation:         SpaceIT GmbH  <<  32 //  Organisation:   SpaceIT GmbH
 34 //        Contract:        ESA contract 21435/ <<  33 //  Contract: ESA contract 21435/08/NL/AT
 35 //         Customer:             ESA/ESTEC     <<  34 //  Customer:       ESA/ESTEC
 36 //////////////////////////////////////////////     35 //////////////////////////////////////////////////////////////
 37                                                << 
 38 //....oooOO0OOooo........oooOO0OOooo........oo << 
 39 //....oooOO0OOooo........oooOO0OOooo........oo << 
 40                                                << 
 41 #include "RMC01AnalysisManager.hh"                 36 #include "RMC01AnalysisManager.hh"
 42                                                << 
 43 #include "RMC01AnalysisManagerMessenger.hh"    << 
 44 #include "RMC01SD.hh"                          << 
 45                                                << 
 46 #include "G4AdjointSimManager.hh"                  37 #include "G4AdjointSimManager.hh"
 47 #include "G4Electron.hh"                       << 
 48 #include "G4Gamma.hh"                          << 
 49 #include "G4PhysicalConstants.hh"              << 
 50 #include "G4Proton.hh"                         << 
 51 #include "G4RunManager.hh"                     << 
 52 #include "G4SDManager.hh"                          38 #include "G4SDManager.hh"
 53 #include "G4SystemOfUnits.hh"                  <<  39 #include "RMC01SD.hh"
 54 #include "G4THitsCollection.hh"                    40 #include "G4THitsCollection.hh"
                                                   >>  41 #include "G4Electron.hh"
                                                   >>  42 #include "G4Proton.hh"
                                                   >>  43 #include "G4Gamma.hh"
                                                   >>  44 #include "Histo1DVar.hh"
                                                   >>  45 #include "Histo2DVar.hh"
 55 #include "G4Timer.hh"                              46 #include "G4Timer.hh"
                                                   >>  47 #include "G4RunManager.hh"
                                                   >>  48 #include "RMC01AnalysisManagerMessenger.hh"
                                                   >>  49 RMC01AnalysisManager* RMC01AnalysisManager::instance = 0;
 56                                                    50 
 57 RMC01AnalysisManager* RMC01AnalysisManager::fI <<  51 ////////////////////////////////////////////////////////////////////////////////
 58                                                <<  52 //
 59 //....oooOO0OOooo........oooOO0OOooo........oo <<  53 RMC01AnalysisManager::RMC01AnalysisManager()  
 60                                                <<  54 { 
 61 RMC01AnalysisManager::RMC01AnalysisManager()   <<  55   
 62   : fAccumulated_edep(0.),                     <<  56   theMsg = new RMC01AnalysisManagerMessenger(this);
 63     fAccumulated_edep2(0.),                    <<  57   
 64     fMean_edep(0.),                            <<  58   //----------------------
 65     fError_mean_edep(0.),                      <<  59   //Creation of histograms
 66     fRelative_error(0.),                       <<  60   //----------------------
 67     fElapsed_time(0.),                         <<  61   
 68     fPrecision_to_reach(0.),                   <<  62   //Energy binning of the histograms : 60 log bins over [1keV-1GeV]
 69     fStop_run_if_precision_reached(true),      <<  63   
 70     fNb_evt_modulo_for_convergence_test(5000), <<  64   G4double bins[61];
 71     fEdep_rmatrix_vs_electron_prim_energy(0),  <<  65   size_t nbin=60;
 72     fElectron_current_rmatrix_vs_electron_prim <<  66   G4double Emin=1.*keV;
 73     fGamma_current_rmatrix_vs_electron_prim_en <<  67   G4double Emax=1.*GeV;
 74     fEdep_rmatrix_vs_gamma_prim_energy(0),     <<  68   for ( size_t i=0; i <= nbin; i++) {
 75     fElectron_current_rmatrix_vs_gamma_prim_en <<  69     G4double val_bin;
 76     fGamma_current_rmatrix_vs_gamma_prim_energ <<  70   val_bin=Emin * std::pow(10., i * std::log10(Emax/Emin)/nbin);
 77     fEdep_rmatrix_vs_proton_prim_energy(0),    <<  71   G4double exp_10=4.-int(std::log10(val_bin));
 78     fElectron_current_rmatrix_vs_proton_prim_e <<  72   G4double factor =std::pow(10., exp_10);
 79     fProton_current_rmatrix_vs_proton_prim_ene <<  73   val_bin=int(factor*val_bin)/factor;
 80     fGamma_current_rmatrix_vs_proton_prim_ener <<  74   bins[i] = val_bin;
 81     fFactoryOn(false),                         <<  75   
 82     fPrimSpectrumType(EXPO),                   <<  76   }
 83     fAlpha_or_E0(.5 * MeV),                    <<  77   
 84     fAmplitude_prim_spectrum(1.),              <<  78   
 85     fEmin_prim_spectrum(1. * keV),             <<  79   //Histograms for :
 86     fEmax_prim_spectrum(20. * MeV),            <<  80   //  1)the forward simulation results 
 87     fAdjoint_sim_mode(true),                   <<  81   //  2)the Reverse MC  simulation results normalised to a user spectrum
 88     fNb_evt_per_adj_evt(2)                     <<  82   //-----------------------------------------------------------------------------------
 89 {                                              <<  83    
 90   fMsg = new RMC01AnalysisManagerMessenger(thi <<  84   edep_vs_prim_ekin = new Histo1DVar("edep_vs_prim_ekin",bins,  nbin+1, LEFT);  
 91                                                <<  85   electron_current  = new Histo1DVar("electron_current",bins,  nbin+1, LEFT); 
                                                   >>  86   proton_current= new Histo1DVar("proton_current",bins,  nbin+1, LEFT); 
                                                   >>  87   gamma_current= new Histo1DVar("gamma_current",bins,  nbin+1, LEFT); 
                                                   >>  88   
                                                   >>  89  
                                                   >>  90   
                                                   >>  91   //Answer matrices for the adjoint simulation only
                                                   >>  92   //-----------------------------------------------
                                                   >>  93   
                                                   >>  94   
                                                   >>  95   //answer matrices for external isotropic e- source
                                                   >>  96   
                                                   >>  97   edep_answer_matrix_vs_electron_prim_energy = new Histo1DVar("edep_answer_matrix_vs_electron_prim_energy", bins, nbin+1, LEFT); 
                                                   >>  98     
                                                   >>  99   electron_current_answer_matrix_vs_electron_prim_energy = new Histo2DVar("electron_current_answer_matrix_vs_electron_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT); 
                                                   >> 100     
                                                   >> 101   gamma_current_answer_matrix_vs_electron_prim_energy = new Histo2DVar("gamma_current_answer_matrix_vs_electron_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT); 
                                                   >> 102     
                                                   >> 103   
                                                   >> 104   //answer matrices for external isotropic gamma source
                                                   >> 105   
                                                   >> 106   edep_answer_matrix_vs_gamma_prim_energy = new Histo1DVar("edep_answer_matrix_vs_gamma_prim_energy", bins, nbin+1, LEFT); 
                                                   >> 107     
                                                   >> 108   electron_current_answer_matrix_vs_gamma_prim_energy = new Histo2DVar("electron_current_answer_matrix_vs_gamma_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT); 
                                                   >> 109     
                                                   >> 110   gamma_current_answer_matrix_vs_gamma_prim_energy = new Histo2DVar("gamma_current_answer_matrix_vs_gamma_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT); 
                                                   >> 111     
                                                   >> 112   
                                                   >> 113   
                                                   >> 114   //answer matrices for external isotropic proton source
                                                   >> 115   
                                                   >> 116   edep_answer_matrix_vs_proton_prim_energy = new Histo1DVar("edep_answer_matrix_vs_proton_prim_energy", bins, nbin+1, LEFT); 
                                                   >> 117     
                                                   >> 118   electron_current_answer_matrix_vs_proton_prim_energy = new Histo2DVar("electron_current_answer_matrix_vs_proton_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT); 
                                                   >> 119     
                                                   >> 120   gamma_current_answer_matrix_vs_proton_prim_energy = new Histo2DVar("gamma_current_answer_matrix_vs_proton_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT); 
                                                   >> 121     
                                                   >> 122   proton_current_answer_matrix_vs_proton_prim_energy = new Histo2DVar("proton_current_answer_matrix_vs_proton_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT); 
                                                   >> 123     
                                                   >> 124   
                                                   >> 125   
                                                   >> 126   //--------------------------
                                                   >> 127   // Convergence test variables
                                                   >> 128   //--------------------------
                                                   >> 129   
                                                   >> 130   precision_to_reach =0.;
                                                   >> 131   stop_run_if_precision_reached =true;
                                                   >> 132   nb_evt_modulo_for_convergence_test =5000;
                                                   >> 133   
                                                   >> 134   
                                                   >> 135   
 92   //-------------                                 136   //-------------
 93   // Timer for convergence vector              << 137   //Timer for convergence vector
 94   //-------------                                 138   //-------------
                                                   >> 139   
                                                   >> 140   theTimer = new G4Timer();
                                                   >> 141  
                                                   >> 142   
                                                   >> 143   
                                                   >> 144   
                                                   >> 145   
 95                                                   146 
 96   fTimer = new G4Timer();                      << 147   
 97                                                << 
 98   //---------------------------------             148   //---------------------------------
 99   // Primary particle ID for normalisation of  << 149   //Primary spectrum for normalisation of adjoint results
100   //---------------------------------             150   //---------------------------------
101                                                << 151   
102   fPrimPDG_ID = G4Electron::Electron()->GetPDG << 152   thePrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
103                                                << 153   alpha_or_E0=.5*MeV;
104   fFileName[0] = "sim";                        << 154   amplitude_prim_spectrum =1.;
                                                   >> 155   emin_prim_spectrum=1.*keV;
                                                   >> 156   emax_prim_spectrum = 20.*MeV;
                                                   >> 157               
                                                   >> 158   
                                                   >> 159 }
                                                   >> 160 ////////////////////////////////////////////////////////////////////////////////
                                                   >> 161 //  
                                                   >> 162 RMC01AnalysisManager::~RMC01AnalysisManager() 
                                                   >> 163 { 
                                                   >> 164   delete edep_vs_prim_ekin;
                                                   >> 165   delete electron_current;
                                                   >> 166   delete proton_current;
                                                   >> 167   delete gamma_current;
                                                   >> 168   
                                                   >> 169   delete edep_answer_matrix_vs_electron_prim_energy;
                                                   >> 170   delete electron_current_answer_matrix_vs_electron_prim_energy;
                                                   >> 171   delete gamma_current_answer_matrix_vs_electron_prim_energy;
                                                   >> 172   
                                                   >> 173   delete edep_answer_matrix_vs_gamma_prim_energy;
                                                   >> 174   delete electron_current_answer_matrix_vs_gamma_prim_energy;
                                                   >> 175   delete gamma_current_answer_matrix_vs_gamma_prim_energy;
                                                   >> 176   
                                                   >> 177   delete edep_answer_matrix_vs_proton_prim_energy;
                                                   >> 178   delete electron_current_answer_matrix_vs_proton_prim_energy;
                                                   >> 179   delete proton_current_answer_matrix_vs_proton_prim_energy;
                                                   >> 180   delete gamma_current_answer_matrix_vs_proton_prim_energy;
                                                   >> 181   
                                                   >> 182   
105 }                                                 183 }
106                                                << 184 ////////////////////////////////////////////////////////////////////////////////
107 //....oooOO0OOooo........oooOO0OOooo........oo << 185 //
108                                                << 
109 RMC01AnalysisManager::~RMC01AnalysisManager()  << 
110 {                                              << 
111   ;                                            << 
112 }                                              << 
113                                                << 
114 //....oooOO0OOooo........oooOO0OOooo........oo << 
115                                                << 
116 RMC01AnalysisManager* RMC01AnalysisManager::Ge    186 RMC01AnalysisManager* RMC01AnalysisManager::GetInstance()
117 {                                                 187 {
118   if (fInstance == 0) fInstance = new RMC01Ana << 188   if (instance == 0) instance = new RMC01AnalysisManager;
119   return fInstance;                            << 189   return instance;
120 }                                                 190 }
121                                                << 191 ////////////////////////////////////////////////////////////////////////////////
122 //....oooOO0OOooo........oooOO0OOooo........oo << 192 //  
123                                                << 
124 void RMC01AnalysisManager::BeginOfRun(const G4    193 void RMC01AnalysisManager::BeginOfRun(const G4Run* aRun)
125 {                                              << 194 {  accumulated_edep =0.;
126   fAccumulated_edep = 0.;                      << 195    accumulated_edep2 =0.;
127   fAccumulated_edep2 = 0.;                     << 196    relative_error=1.;
128   fNentry = 0.0;                               << 197    mean_edep=0.;
129   fRelative_error = 1.;                        << 198    error_mean_edep=0.;
130   fMean_edep = 0.;                             << 199    adjoint_sim_mode =G4AdjointSimManager::GetInstance()->GetAdjointSimMode();
131   fError_mean_edep = 0.;                       << 200    if (adjoint_sim_mode){
132   fAdjoint_sim_mode = G4AdjointSimManager::Get << 201     nb_evt_per_adj_evt=aRun->GetNumberOfEventToBeProcessed()/G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
133                                                << 202   ConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt", std::ios::out);
134   if (fAdjoint_sim_mode) {                     << 203     ConvergenceFileOutput<<"Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
135     fNb_evt_per_adj_evt = aRun->GetNumberOfEve << 204    }
136                           / G4AdjointSimManage << 205    else {
137     fConvergenceFileOutput.open("ConvergenceOf << 206     ConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt", std::ios::out);
138     fConvergenceFileOutput << "Normalised Edep << 207     ConvergenceFileOutput<<"Edep per event [MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
139   }                                            << 208    }
140   else {                                       << 209    ConvergenceFileOutput.setf(std::ios::scientific);
141     fConvergenceFileOutput.open("ConvergenceOf << 210    ConvergenceFileOutput.precision(6);  
142     fConvergenceFileOutput << "Edep per event  << 211    ResetHistograms();
143   }                                            << 212    theTimer->Start();
144   fConvergenceFileOutput.setf(std::ios::scient << 213    elapsed_time=0.;
145   fConvergenceFileOutput.precision(6);         << 
146                                                << 
147   fTimer->Start();                             << 
148   fElapsed_time = 0.;                          << 
149                                                << 
150   Book();                                      << 
151 }                                                 214 }
152                                                << 215 ////////////////////////////////////////////////////////////////////////////////
153 //....oooOO0OOooo........oooOO0OOooo........oo << 216 //  
154                                                << 
155 void RMC01AnalysisManager::EndOfRun(const G4Ru    217 void RMC01AnalysisManager::EndOfRun(const G4Run* aRun)
156 {                                              << 218 { theTimer->Stop();
157   fTimer->Stop();                              << 219  
158   G4int nb_evt = aRun->GetNumberOfEvent();     << 220   
159   G4double factor = 1. / nb_evt;               << 221   if (!adjoint_sim_mode){
160   if (!fAdjoint_sim_mode) {                    << 222         
161     G4cout << "Results of forward simulation!" << 223   G4cout<<"Results of forward simulation!"<<std::endl;
162     G4cout << "edep per event [MeV] = " << fMe << 224   G4cout<<"edep per event [MeV] = "<<mean_edep<<std::endl;
163     G4cout << "error[MeV] = " << fError_mean_e << 225   G4cout<<"error[MeV] = "<<error_mean_edep<<std::endl;
                                                   >> 226   G4int nb_evt=aRun->GetNumberOfEvent();
                                                   >> 227   WriteHisto(edep_vs_prim_ekin,1./nb_evt, G4String("Fwd_Edep_vs_EkinPrim.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tEdep[MeV]\terr_Edep[MeV]\n"));
                                                   >> 228   WriteHisto(electron_current,1./nb_evt, G4String("Fwd_ElectronCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering electron\terr\n"));
                                                   >> 229   WriteHisto(proton_current,1./nb_evt, G4String("Fwd_ProtonCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering proton\terr[]\n"));
                                                   >> 230     WriteHisto(gamma_current,1./nb_evt, G4String("Fwd_GammaCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering gamma\terr[]\n"));
164   }                                               231   }
165                                                << 232   
166   else {                                          233   else {
167     G4cout << "Results of reverse/adjoint simu << 234     G4cout<<"Results of reverse/adjoint simulation!"<<std::endl;
168     G4cout << "normalised edep [MeV] = " << fM << 235   G4cout<<"normalised edep [MeV] = "<<mean_edep<<std::endl;
169     G4cout << "error[MeV] = " << fError_mean_e << 236   G4cout<<"error[MeV] = "<<error_mean_edep<<std::endl;
170     factor = 1. * G4AdjointSimManager::GetInst << 237 
171              / aRun->GetNumberOfEvent();       << 238   
172   }                                            << 239   G4double factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun()*nb_evt_per_adj_evt/aRun->GetNumberOfEvent();
173   Save(factor);                                << 240     
174   fConvergenceFileOutput.close();              << 241   
                                                   >> 242   
                                                   >> 243   WriteHisto(edep_vs_prim_ekin,factor, G4String("Adj_Edep_vs_EkinPrim.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tEdep[MeV]\terr_Edep[MeV]\n"));
                                                   >> 244   WriteHisto(electron_current,factor, G4String("Adj_ElectronCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering electron\terr\n"));
                                                   >> 245   WriteHisto(proton_current,factor, G4String("Adj_ProtonCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering proton\terr[]\n"));
                                                   >> 246     WriteHisto(gamma_current,factor, G4String("Adj_GammaCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering gamma\terr[]\n"));
                                                   >> 247   
                                                   >> 248     
                                                   >> 249   WriteHisto(edep_answer_matrix_vs_electron_prim_energy,factor,
                                                   >> 250              G4String("Adj_Edep_vs_EkinPrimElectron_Answer.txt"),
                                                   >> 251              G4String("E1[MeV]\t\tE2[MeV]\t\tEdep Efficiency[MeV*cm2*MeV*str]\terr_Edep[MeV*cm2*MeV*str]\n"));
                                                   >> 252   
                                                   >> 253     WriteHisto(electron_current_answer_matrix_vs_electron_prim_energy,factor,
                                                   >> 254              G4String("Adj_ElectronCurrent_vs_EkinPrimElectron_Answer.txt"),
                                                   >> 255              G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
                                                   >> 256   
                                                   >> 257   WriteHisto(gamma_current_answer_matrix_vs_electron_prim_energy,factor,
                                                   >> 258              G4String("Adj_GammaCurrent_vs_EkinPrimElectron_Answer.txt"),
                                                   >> 259              G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
                                                   >> 260   
                                                   >> 261     
                                                   >> 262   
                                                   >> 263   WriteHisto(edep_answer_matrix_vs_gamma_prim_energy,factor,
                                                   >> 264              G4String("Adj_Edep_vs_EkinPrimGamma_Answer.txt"),
                                                   >> 265              G4String("E1[MeV]\t\tE2[MeV]\t\tEdep Efficiency[MeV*cm2*MeV*str]\terr_Edep[MeV*cm2*MeV*str]\n"));
                                                   >> 266   
                                                   >> 267     WriteHisto(electron_current_answer_matrix_vs_gamma_prim_energy,factor,
                                                   >> 268              G4String("Adj_ElectronCurrent_vs_EkinPrimGamma_Answer.txt"),
                                                   >> 269              G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
                                                   >> 270   
                                                   >> 271   WriteHisto(gamma_current_answer_matrix_vs_gamma_prim_energy,factor,
                                                   >> 272              G4String("Adj_GammaCurrent_vs_EkinPrimGamma_Answer.txt"),
                                                   >> 273              G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
                                                   >> 274   
                                                   >> 275   
                                                   >> 276   
                                                   >> 277   WriteHisto(edep_answer_matrix_vs_proton_prim_energy,factor,
                                                   >> 278              G4String("Adj_Edep_vs_EkinPrimProton_Answer.txt"),
                                                   >> 279              G4String("E1[MeV]\t\tE2[MeV]\t\tEdep Efficiency[MeV*cm2*MeV*str]\terr_Edep[MeV*cm2*MeV*str]\n"));
                                                   >> 280   
                                                   >> 281     WriteHisto(electron_current_answer_matrix_vs_proton_prim_energy,factor,
                                                   >> 282              G4String("Adj_ElectronCurrent_vs_EkinPrimProton_Answer.txt"),
                                                   >> 283              G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
                                                   >> 284   
                                                   >> 285   WriteHisto(gamma_current_answer_matrix_vs_proton_prim_energy,factor,
                                                   >> 286              G4String("Adj_GammaCurrent_vs_EkinPrimProton_Answer.txt"),
                                                   >> 287              G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
                                                   >> 288   
                                                   >> 289   WriteHisto(proton_current_answer_matrix_vs_proton_prim_energy,factor,
                                                   >> 290              G4String("Adj_ProtonCurrent_vs_EkinPrimProton_Answer.txt"),
                                                   >> 291              G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
                                                   >> 292   }
                                                   >> 293   ConvergenceFileOutput.close();
                                                   >> 294 }
                                                   >> 295 ////////////////////////////////////////////////////////////////////////////////
                                                   >> 296 //  
                                                   >> 297 void RMC01AnalysisManager::BeginOfEvent(const G4Event* )
                                                   >> 298 {;   
175 }                                                 299 }
176                                                << 300 ////////////////////////////////////////////////////////////////////////////////
177 //....oooOO0OOooo........oooOO0OOooo........oo << 301 //  
178                                                << 
179 void RMC01AnalysisManager::BeginOfEvent(const  << 
180 {                                              << 
181   ;                                            << 
182 }                                              << 
183                                                << 
184 //....oooOO0OOooo........oooOO0OOooo........oo << 
185                                                << 
186 void RMC01AnalysisManager::EndOfEvent(const G4    302 void RMC01AnalysisManager::EndOfEvent(const G4Event* anEvent)
187 {                                              << 303 {  
188   if (fAdjoint_sim_mode)                       << 
189     EndOfEventForAdjointSimulation(anEvent);   << 
190   else                                         << 
191     EndOfEventForForwardSimulation(anEvent);   << 
192                                                << 
193   // Test convergence. The error is already co << 
194   //--------------------------------------     << 
195   G4int nb_event = anEvent->GetEventID() + 1;  << 
196   // G4double factor=1.;                       << 
197   if (fAdjoint_sim_mode) {                     << 
198     G4double n_adj_evt = nb_event / fNb_evt_pe << 
199     // nb_event/fNb_evt_per_adj_evt;           << 
200     if (n_adj_evt * fNb_evt_per_adj_evt == nb_ << 
201       nb_event = static_cast<G4int>(n_adj_evt) << 
202     }                                          << 
203     else                                       << 
204       nb_event = 0;                            << 
205   }                                            << 
206                                                << 
207   if (nb_event > 100 && fStop_run_if_precision << 
208     G4cout << fPrecision_to_reach * 100. << "% << 
209     fTimer->Stop();                            << 
210     fElapsed_time += fTimer->GetRealElapsed(); << 
211     fConvergenceFileOutput << fMean_edep << '\ << 
212                            << std::endl;       << 
213     G4RunManager::GetRunManager()->AbortRun(tr << 
214   }                                            << 
215                                                   304 
216   if (nb_event > 0 && nb_event % fNb_evt_modul << 
217     fTimer->Stop();                            << 
218     fElapsed_time += fTimer->GetRealElapsed(); << 
219     fTimer->Start();                           << 
220     fConvergenceFileOutput << fMean_edep << '\ << 
221                            << std::endl;       << 
222   }                                            << 
223 }                                              << 
224                                                << 
225 //....oooOO0OOooo........oooOO0OOooo........oo << 
226                                                << 
227 void RMC01AnalysisManager::EndOfEventForForwar << 
228 {                                              << 
229   G4SDManager* SDman = G4SDManager::GetSDMpoin << 
230   G4HCofThisEvent* HCE = anEvent->GetHCofThisE << 
231   RMC01DoubleWithWeightHitsCollection* edepCol << 
232     (RMC01DoubleWithWeightHitsCollection*)(HCE << 
233                                                   305 
234   RMC01DoubleWithWeightHitsCollection* electro << 
235     (RMC01DoubleWithWeightHitsCollection*)(HCE << 
236                                                   306 
237   RMC01DoubleWithWeightHitsCollection* protonC << 307    if (adjoint_sim_mode) EndOfEventForAdjointSimulation(anEvent);
238     (RMC01DoubleWithWeightHitsCollection*)(HCE << 308    else EndOfEventForForwardSimulation(anEvent);
239                                                << 309    
240   RMC01DoubleWithWeightHitsCollection* gammaCu << 310    
241     (RMC01DoubleWithWeightHitsCollection*)(HCE << 311    //Test convergence. The error is already computed
242                                                << 312    //--------------------------------------
243   // Total energy deposited in Event           << 313    G4int nb_event=anEvent->GetEventID()+1;
244   //-------------------------------            << 314    G4double factor=1.;
245   G4double totEdep = 0;                        << 315    if (adjoint_sim_mode) {
246   std::size_t i;                               << 316     G4double  n_adj_evt= nb_event/nb_evt_per_adj_evt;
247   for (i = 0; i < edepCollection->entries(); + << 317   // nb_event/nb_evt_per_adj_evt;
248     totEdep += (*edepCollection)[i]->GetValue( << 318   if (n_adj_evt*nb_evt_per_adj_evt == nb_event) {
249                                                << 319     nb_event =static_cast<G4int>(n_adj_evt);
250   if (totEdep > 0.) {                          << 320     factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
251     fAccumulated_edep += totEdep;              << 321   } 
252     fAccumulated_edep2 += totEdep * totEdep;   << 322   else nb_event=0;
253     fNentry += 1.0;                            << 323   
254     G4PrimaryParticle* thePrimary = anEvent->G << 324    }
255     G4double E0 = thePrimary->GetG4code()->Get << 325    
256     G4double P = thePrimary->GetMomentum().mag << 326    
257     G4double prim_ekin = std::sqrt(E0 * E0 + P << 327    
258     fEdep_vs_prim_ekin->fill(prim_ekin, totEde << 328    if (nb_event>100 && stop_run_if_precision_reached && precision_to_reach >relative_error) {
259   }                                            << 329     G4cout<<precision_to_reach*100.<<"%  Precision reached!"<<std::endl;
260   ComputeMeanEdepAndError(anEvent, fMean_edep, << 330     theTimer->Stop();
261   if (fError_mean_edep > 0) fRelative_error =  << 331     elapsed_time+=theTimer->GetRealElapsed();
262                                                << 332     ConvergenceFileOutput<<mean_edep<<'\t'<<error_mean_edep<<'\t'<<elapsed_time<<std::endl;
263   // Particle current on sensitive cylinder    << 333     G4RunManager::GetRunManager()->AbortRun(true);
264   //-------------------------------------      << 334    }
265                                                << 335    
266   for (i = 0; i < electronCurrentCollection->e << 336    
267     G4double ekin = (*electronCurrentCollectio << 337    if (nb_event>0 && nb_event % nb_evt_modulo_for_convergence_test == 0) {
268     G4double weight = (*electronCurrentCollect << 338     theTimer->Stop();
269     fElectron_current->fill(ekin, weight);     << 339   elapsed_time+=theTimer->GetRealElapsed();
270   }                                            << 340   theTimer->Start();
271                                                << 341   ConvergenceFileOutput<<mean_edep<<'\t'<<error_mean_edep<<'\t'<<elapsed_time<<std::endl;
272   for (i = 0; i < protonCurrentCollection->ent << 342   
273     G4double ekin = (*protonCurrentCollection) << 343    }  
274     G4double weight = (*protonCurrentCollectio << 344     
275     fProton_current->fill(ekin, weight);       << 345   
276   }                                            << 346    
277                                                << 347   
278   for (i = 0; i < gammaCurrentCollection->entr << 348   
279     G4double ekin = (*gammaCurrentCollection)[ << 349    
280     G4double weight = (*gammaCurrentCollection << 350 }   
281     fGamma_current->fill(ekin, weight);        << 351 ////////////////////////////////////////////////////////////////////////////////
282   }                                            << 352 // 
283 }                                              << 353 void  RMC01AnalysisManager::EndOfEventForForwardSimulation(const G4Event* anEvent)
284                                                << 354 {  
285 //....oooOO0OOooo........oooOO0OOooo........oo << 355    
286                                                << 356    G4SDManager* SDman = G4SDManager::GetSDMpointer();
287 void RMC01AnalysisManager::EndOfEventForAdjoin << 357    G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
288 {                                              << 358    RMC01DoubleWithWeightHitsCollection* EdepCollection =  (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("edep")));
289   // Output from Sensitive volume computed dur << 359    RMC01DoubleWithWeightHitsCollection* ElectronCurrentCollection =  (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_electron")));
                                                   >> 360    RMC01DoubleWithWeightHitsCollection* ProtonCurrentCollection =  (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_proton")));
                                                   >> 361    RMC01DoubleWithWeightHitsCollection* GammaCurrentCollection =  (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_gamma")));
                                                   >> 362    
                                                   >> 363    //Total energy deposited in Event
                                                   >> 364    //-------------------------------
                                                   >> 365    G4double totEdep=0; 
                                                   >> 366    G4int i;
                                                   >> 367    for (i=0;i<EdepCollection->entries();i++) totEdep+=(*EdepCollection)[i]->GetValue()*(*EdepCollection)[i]->GetWeight();
                                                   >> 368    
                                                   >> 369    if (totEdep>0.){
                                                   >> 370     accumulated_edep +=totEdep ;
                                                   >> 371     accumulated_edep2 +=totEdep*totEdep;
                                                   >> 372         G4PrimaryParticle* thePrimary=anEvent->GetPrimaryVertex()->GetPrimary();
                                                   >> 373     G4double E0= thePrimary->GetG4code()->GetPDGMass();
                                                   >> 374     G4double P=thePrimary->GetMomentum().mag();
                                                   >> 375     G4double prim_ekin =std::sqrt(E0*E0+P*P)-E0;
                                                   >> 376     edep_vs_prim_ekin->fill(prim_ekin,totEdep);
                                                   >> 377    } 
                                                   >> 378    ComputeMeanEdepAndError(anEvent,mean_edep,error_mean_edep);
                                                   >> 379    if (error_mean_edep>0) relative_error= error_mean_edep/mean_edep;
                                                   >> 380       
                                                   >> 381    
                                                   >> 382    
                                                   >> 383    //Particle current on sensitive cylinder
                                                   >> 384    //-------------------------------------
                                                   >> 385    
                                                   >> 386    for (i=0;i<ElectronCurrentCollection->entries();i++) {
                                                   >> 387     G4double ekin =(*ElectronCurrentCollection)[i]->GetValue();
                                                   >> 388   G4double weight=(*ElectronCurrentCollection)[i]->GetWeight();
                                                   >> 389   electron_current->fill(ekin,weight);
                                                   >> 390    }
                                                   >> 391    
                                                   >> 392    for (i=0;i<ProtonCurrentCollection->entries();i++) {
                                                   >> 393     G4double ekin =(*ProtonCurrentCollection)[i]->GetValue();
                                                   >> 394   G4double weight=(*ProtonCurrentCollection)[i]->GetWeight();
                                                   >> 395   proton_current->fill(ekin,weight);
                                                   >> 396    }  
                                                   >> 397    
                                                   >> 398    for (i=0;i<GammaCurrentCollection->entries();i++) {
                                                   >> 399     G4double ekin =(*GammaCurrentCollection)[i]->GetValue();
                                                   >> 400   G4double weight=(*GammaCurrentCollection)[i]->GetWeight();
                                                   >> 401   gamma_current->fill(ekin,weight);
                                                   >> 402    }
                                                   >> 403    
                                                   >> 404    
                                                   >> 405 }
                                                   >> 406 ////////////////////////////////////////////////////////////////////////////////
                                                   >> 407 // 
                                                   >> 408 void  RMC01AnalysisManager::EndOfEventForAdjointSimulation(const G4Event* anEvent)
                                                   >> 409 {  
                                                   >> 410   //Output from Sensitive volume computed during the forward tracking phase
290   //------------------------------------------    411   //-----------------------------------------------------------------------
                                                   >> 412   
291   G4SDManager* SDman = G4SDManager::GetSDMpoin    413   G4SDManager* SDman = G4SDManager::GetSDMpointer();
292   G4HCofThisEvent* HCE = anEvent->GetHCofThisE    414   G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
293   RMC01DoubleWithWeightHitsCollection* edepCol << 415   RMC01DoubleWithWeightHitsCollection* EdepCollection =  (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("edep")));
294     (RMC01DoubleWithWeightHitsCollection*)(HCE << 416   RMC01DoubleWithWeightHitsCollection* ElectronCurrentCollection =  (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_electron")));
295                                                << 417   RMC01DoubleWithWeightHitsCollection* ProtonCurrentCollection =  (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_proton")));
296   RMC01DoubleWithWeightHitsCollection* electro << 418   RMC01DoubleWithWeightHitsCollection* GammaCurrentCollection =  (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_gamma")));
297     (RMC01DoubleWithWeightHitsCollection*)(HCE << 419   
298                                                << 420   //Output from adjoint tracking phase
299   RMC01DoubleWithWeightHitsCollection* protonC << 
300     (RMC01DoubleWithWeightHitsCollection*)(HCE << 
301                                                << 
302   RMC01DoubleWithWeightHitsCollection* gammaCu << 
303     (RMC01DoubleWithWeightHitsCollection*)(HCE << 
304                                                << 
305   // Computation of total energy deposited in  << 
306   //-------------------------------            << 
307   G4double totEdep = 0;                        << 
308   std::size_t i;                               << 
309   for (i = 0; i < edepCollection->entries(); + << 
310     totEdep += (*edepCollection)[i]->GetValue( << 
311                                                << 
312   // Output from adjoint tracking phase        << 
313   //------------------------------------------    421   //----------------------------------------------------------------------------
314                                                << 422   
315   G4AdjointSimManager* theAdjointSimManager =     423   G4AdjointSimManager* theAdjointSimManager = G4AdjointSimManager::GetInstance();
316                                                << 424   G4int pdg_nb =theAdjointSimManager->GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack();
317   size_t nb_adj_track = theAdjointSimManager-> << 425   G4double prim_ekin=theAdjointSimManager->GetEkinAtEndOfLastAdjointTrack();
318   G4double total_normalised_weight = 0.;       << 426   G4double adj_weight=theAdjointSimManager->GetWeightAtEndOfLastAdjointTrack();
319                                                << 427  
320   // We need to loop over the adjoint tracks t << 428   
321   // surface.                                  << 429   //Factor of normalisation to user selected primary spectrum (power law or exponential) 
322   for (std::size_t j = 0; j < nb_adj_track; ++ << 430   //------------------------------------------------------------------------------------
323     G4int pdg_nb = theAdjointSimManager->GetFw << 431   
324     G4double prim_ekin = theAdjointSimManager- << 432   G4double normalised_weight = 0.;
325     G4double adj_weight = theAdjointSimManager << 433   if (pdg_nb== thePrimPDG_ID && prim_ekin>= emin_prim_spectrum && prim_ekin<= emax_prim_spectrum) 
326                                                << 434         normalised_weight = adj_weight*PrimDiffAndDirectionalFluxForAdjointSim(prim_ekin);
327     // Factor of normalisation to user defined << 435   
328     //---------------------------------------- << 436   
329     G4double normalised_weight = 0.;           << 437   //Answer matrices
330     if (pdg_nb == fPrimPDG_ID && prim_ekin >=  << 438   //-------------
331         && prim_ekin <= fEmax_prim_spectrum)   << 439    Histo1DVar* edep_answer_matrix =0;
332       normalised_weight = adj_weight * PrimDif << 440    Histo2DVar* electron_current_answer_matrix =0;
333     total_normalised_weight += normalised_weig << 441    Histo2DVar* gamma_current_answer_matrix =0;
334                                                << 442    Histo2DVar* proton_current_answer_matrix =0;
335     // Answer matrices                         << 443   
336     //-------------                            << 444    
337     G4H1* edep_rmatrix = 0;                    << 445    if (pdg_nb == G4Electron::Electron()->GetPDGEncoding()){ //electron answer matrices
338     G4H2* electron_current_rmatrix = 0;        << 446   edep_answer_matrix = edep_answer_matrix_vs_electron_prim_energy;
339     G4H2* gamma_current_rmatrix = 0;           << 447   electron_current_answer_matrix = electron_current_answer_matrix_vs_electron_prim_energy;
340     G4H2* proton_current_rmatrix = 0;          << 448   gamma_current_answer_matrix = gamma_current_answer_matrix_vs_electron_prim_energy;
341                                                << 449    }
342     if (pdg_nb == G4Electron::Electron()->GetP << 450    else if (pdg_nb == G4Gamma::Gamma()->GetPDGEncoding()){ //electron answer matrices
343       edep_rmatrix = fEdep_rmatrix_vs_electron << 451   edep_answer_matrix = edep_answer_matrix_vs_gamma_prim_energy;
344       electron_current_rmatrix = fElectron_cur << 452   electron_current_answer_matrix = electron_current_answer_matrix_vs_gamma_prim_energy;
345       gamma_current_rmatrix = fGamma_current_r << 453   gamma_current_answer_matrix = gamma_current_answer_matrix_vs_gamma_prim_energy;
346     }                                          << 454    }
347     else if (pdg_nb == G4Gamma::Gamma()->GetPD << 455    else if (pdg_nb == G4Proton::Proton()->GetPDGEncoding()){ //electron answer matrices
348       // gammma answer matrices                << 456   edep_answer_matrix = edep_answer_matrix_vs_proton_prim_energy;
349       edep_rmatrix = fEdep_rmatrix_vs_gamma_pr << 457   electron_current_answer_matrix = electron_current_answer_matrix_vs_proton_prim_energy;
350       electron_current_rmatrix = fElectron_cur << 458   gamma_current_answer_matrix = gamma_current_answer_matrix_vs_proton_prim_energy;
351       gamma_current_rmatrix = fGamma_current_r << 459   proton_current_answer_matrix = proton_current_answer_matrix_vs_proton_prim_energy;
352     }                                          << 460    }
353     else if (pdg_nb == G4Proton::Proton()->Get << 461   
354       // proton answer matrices                << 462   
355       edep_rmatrix = fEdep_rmatrix_vs_proton_p << 463   
356       electron_current_rmatrix = fElectron_cur << 464   //Registering of total energy deposited in Event
357       gamma_current_rmatrix = fGamma_current_r << 
358       proton_current_rmatrix = fProton_current << 
359     }                                          << 
360     // Register histo edep vs prim ekin        << 
361     //----------------------------------       << 
362     if (normalised_weight > 0) fEdep_vs_prim_e << 
363     // Registering answer matrix               << 
364     //---------------------------              << 
365     edep_rmatrix->fill(prim_ekin, totEdep * ad << 
366                                                << 
367     // Registering of current of particles on  << 
368     //---------------------------------------- << 
369                                                << 
370     for (i = 0; i < electronCurrentCollection- << 
371       G4double ekin = (*electronCurrentCollect << 
372       G4double weight = (*electronCurrentColle << 
373       fElectron_current->fill(ekin, weight * n << 
374       electron_current_rmatrix->fill(prim_ekin << 
375     }                                          << 
376     for (i = 0; i < protonCurrentCollection->e << 
377       G4double ekin = (*protonCurrentCollectio << 
378       G4double weight = (*protonCurrentCollect << 
379       fProton_current->fill(ekin, weight * nor << 
380       proton_current_rmatrix->fill(prim_ekin,  << 
381     }                                          << 
382     for (i = 0; i < gammaCurrentCollection->en << 
383       G4double ekin = (*gammaCurrentCollection << 
384       G4double weight = (*gammaCurrentCollecti << 
385       fGamma_current->fill(ekin, weight * norm << 
386       gamma_current_rmatrix->fill(prim_ekin, e << 
387     }                                          << 
388   }                                            << 
389                                                << 
390   // Registering of total energy deposited in  << 
391   //-------------------------------               465   //-------------------------------
392   G4bool new_mean_computed = false;            << 466    G4double totEdep=0; 
393   if (totEdep > 0.) {                          << 467    G4int i;
394     if (total_normalised_weight > 0.) {        << 468    for (i=0;i<EdepCollection->entries();i++) totEdep+=(*EdepCollection)[i]->GetValue()*(*EdepCollection)[i]->GetWeight();
395       G4double edep = totEdep * total_normalis << 469    
396                                                << 470    G4bool new_mean_computed=false;
397       // Check if the edep is not wrongly too  << 471    if (totEdep>0.){
398       //-------------------------------------- << 472     if (normalised_weight>0.){
399       G4double new_mean(0.0), new_error(0.0);  << 473     G4double edep=totEdep* normalised_weight;
400       fAccumulated_edep += edep;               << 474     
401       fAccumulated_edep2 += edep * edep;       << 475     //Check if the edep is not wrongly too high
402       fNentry += 1.0;                          << 476     //----------------------------------------- 
403       ComputeMeanEdepAndError(anEvent, new_mea << 477     G4double new_mean , new_error;
404       G4double new_relative_error = 1.;        << 478     
405       if (new_error > 0) new_relative_error =  << 479     accumulated_edep +=edep;
406       if (fRelative_error < 0.10 && new_relati << 480       accumulated_edep2 +=edep*edep;
407         G4cout << "Potential wrong adjoint wei << 481     
408         G4cout << "The results of this event w << 482     ComputeMeanEdepAndError(anEvent,new_mean,new_error);
409         G4cout << "previous mean edep [MeV] "  << 483     G4double new_relative_error = 1.;
410         G4cout << "previous relative error " < << 484     if ( new_error >0) new_relative_error = new_error/ new_mean;
411         G4cout << "new rejected mean edep [MeV << 485     if (relative_error <0.10 && new_relative_error>1.5*relative_error) { // rejected
412         G4cout << "new rejected relative error << 486       G4cout<<"Potential wrong adjoint weight!"<<std::endl;
413         fAccumulated_edep -= edep;             << 487       G4cout<<"The results of this event will not be registered!"<<std::endl;
414         fAccumulated_edep2 -= edep * edep;     << 488       G4cout<<"previous mean edep [MeV] "<< mean_edep<<std::endl;
415         fNentry -= 1.0;                        << 489       G4cout<<"previous relative error "<< relative_error<<std::endl;
416         return;                                << 490       G4cout<<"new rejected mean edep [MeV] "<< new_mean<<std::endl;
417       }                                        << 491       G4cout<<"new rejected relative error "<< new_relative_error<<std::endl;
418       else {  // accepted                      << 492       accumulated_edep -=edep;
419         fMean_edep = new_mean;                 << 493       accumulated_edep2 -=edep*edep;
420         fError_mean_edep = new_error;          << 494       return; 
421         fRelative_error = new_relative_error;  << 495     }
422         new_mean_computed = true;              << 496     else { //accepted
423       }                                        << 497       mean_edep = new_mean;
424     }                                          << 498       error_mean_edep = new_error;
425                                                << 499       relative_error =new_relative_error;
426     if (!new_mean_computed) {                  << 500       new_mean_computed=true;
427       ComputeMeanEdepAndError(anEvent, fMean_e << 501     }  
428       fRelative_error = (fMean_edep > 0.0) ? f << 502     edep_vs_prim_ekin->fill(prim_ekin,edep); 
429     }                                          << 503   }
430   }                                            << 504   
                                                   >> 505   // Registering answer matrix
                                                   >> 506   //---------------------------
                                                   >> 507   
                                                   >> 508   edep_answer_matrix->fill(prim_ekin,totEdep*adj_weight/cm2);
                                                   >> 509   
                                                   >> 510 
                                                   >> 511    }
                                                   >> 512    if (!new_mean_computed){
                                                   >> 513      ComputeMeanEdepAndError(anEvent,mean_edep,error_mean_edep);
                                                   >> 514      if (error_mean_edep>0) relative_error= error_mean_edep/mean_edep;
                                                   >> 515    }   
                                                   >> 516     
                                                   >> 517    
                                                   >> 518   //Registering of current of particles on the sensitive volume
                                                   >> 519   //------------------------------------------------------------
                                                   >> 520    
                                                   >> 521    for (i=0;i<ElectronCurrentCollection->entries();i++) {
                                                   >> 522     G4double ekin =(*ElectronCurrentCollection)[i]->GetValue();
                                                   >> 523   G4double weight=(*ElectronCurrentCollection)[i]->GetWeight();
                                                   >> 524   electron_current->fill(ekin,weight*normalised_weight);
                                                   >> 525   electron_current_answer_matrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
                                                   >> 526    }
                                                   >> 527    
                                                   >> 528    for (i=0;i<ProtonCurrentCollection->entries();i++) {
                                                   >> 529     G4double ekin =(*ProtonCurrentCollection)[i]->GetValue();
                                                   >> 530   G4double weight=(*ProtonCurrentCollection)[i]->GetWeight();
                                                   >> 531   proton_current->fill(ekin,weight*normalised_weight);
                                                   >> 532   proton_current_answer_matrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
                                                   >> 533   } 
                                                   >> 534    
                                                   >> 535    for (i=0;i<GammaCurrentCollection->entries();i++) {
                                                   >> 536     G4double ekin =(*GammaCurrentCollection)[i]->GetValue();
                                                   >> 537   G4double weight=(*GammaCurrentCollection)[i]->GetWeight();
                                                   >> 538   gamma_current->fill(ekin,weight*normalised_weight);
                                                   >> 539   gamma_current_answer_matrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
                                                   >> 540 
                                                   >> 541    }
                                                   >> 542       
                                                   >> 543   
431 }                                                 544 }
432                                                << 545 ////////////////////////////////////////////////////////////////////////////////
433 //....oooOO0OOooo........oooOO0OOooo........oo << 546 //
434                                                << 547 G4double RMC01AnalysisManager::PrimDiffAndDirectionalFluxForAdjointSim(G4double prim_energy)
435 G4double RMC01AnalysisManager::PrimDiffAndDirF << 548 { 
436 {                                              << 549   G4double flux=amplitude_prim_spectrum;
437   G4double flux = fAmplitude_prim_spectrum;    << 550   if ( the_prim_spectrum_type ==EXPO) flux*=std::exp(-prim_energy/alpha_or_E0);
438   if (fPrimSpectrumType == EXPO)               << 551   else flux*=std::pow(prim_energy, -alpha_or_E0);
439     flux *= std::exp(-prim_energy / fAlpha_or_ << 
440   else                                         << 
441     flux *= std::pow(prim_energy, -fAlpha_or_E << 
442   return flux;                                    552   return flux;
443 }                                                 553 }
444                                                << 554 ////////////////////////////////////////////////////////////////////////////////
445 //....oooOO0OOooo........oooOO0OOooo........oo << 555 //
446 /*                                             << 556 void  RMC01AnalysisManager::WriteHisto(Histo1DVar* anHisto, G4double scaling_factor, G4String fileName, G4String header_lines)
447 void  RMC01AnalysisManager::WriteHisto(G4H1* a << 
448             G4double scaling_factor, G4String  << 
449 { std::fstream FileOutput(fileName, std::ios::    557 { std::fstream FileOutput(fileName, std::ios::out);
450   FileOutput<<header_lines;                       558   FileOutput<<header_lines;
                                                   >> 559   
451   FileOutput.setf(std::ios::scientific);          560   FileOutput.setf(std::ios::scientific);
452   FileOutput.precision(6);                        561   FileOutput.precision(6);
453                                                   562 
454   for (G4int i =0;i<G4int(anHisto->axis().bins << 563   size_t nxbins = anHisto->part.total_bins();
455         FileOutput<<anHisto->axis().bin_lower_ << 564   for (size_t i =0;i<nxbins;i++) {
456               <<'\t'<<anHisto->axis().bin_uppe << 565   FileOutput<<anHisto->part.get_bin_position(i)
457               <<'\t'<<anHisto->bin_height(i)*s << 566         <<'\t'<<anHisto->part.get_bin_position(i+1)
458               <<'\t'<<anHisto->bin_error(i)*sc << 567               <<'\t'<<anHisto->get_bin_value(int(i))*scaling_factor
                                                   >> 568         <<'\t'<<anHisto->get_bin_error(int(i))*scaling_factor<<std::endl;
459   }                                               569   }
460 }                                                 570 }
461                                                << 571 ////////////////////////////////////////////////////////////////////////////////
462 //....oooOO0OOooo........oooOO0OOooo........oo << 572 //
463                                                << 573 void  RMC01AnalysisManager::WriteHisto(Histo2DVar* anHisto, G4double scaling_factor, G4String fileName, G4String header_lines)
464 void  RMC01AnalysisManager::WriteHisto(G4H2* a << 
465             G4double scaling_factor, G4String  << 
466 { std::fstream FileOutput(fileName, std::ios::    574 { std::fstream FileOutput(fileName, std::ios::out);
467   FileOutput<<header_lines;                       575   FileOutput<<header_lines;
468                                                << 576   
469   FileOutput.setf(std::ios::scientific);          577   FileOutput.setf(std::ios::scientific);
470   FileOutput.precision(6);                        578   FileOutput.precision(6);
471                                                << 579   
472   for (G4int i =0;i<G4int(anHisto->axis_x().bi << 580   size_t nxbins = anHisto->x_part.total_bins();
473     for (G4int j =0;j<G4int(anHisto->axis_y(). << 581   size_t nybins = anHisto->y_part.total_bins();
474        FileOutput<<anHisto->axis_x().bin_lower << 582   for (size_t i =0;i<nxbins;i++) {
475                      <<'\t'<<anHisto->axis_x() << 583       for (size_t j =0;j<nybins;j++) {
476                    <<'\t'<<anHisto->axis_y().b << 584     FileOutput<<anHisto->x_part.get_bin_position(i)
477                    <<'\t'<<anHisto->axis_y().b << 585         <<'\t'<<anHisto->x_part.get_bin_position(i+1)
478                    <<'\t'<<anHisto->bin_height << 586         <<'\t'<<anHisto->y_part.get_bin_position(j)
479                 <<'\t'<<anHisto->bin_error(i,j << 587         <<'\t'<<anHisto->y_part.get_bin_position(j+1)
480                                                << 588         <<'\t'<<anHisto->get_bin_value(int(i),int(j))*scaling_factor
481         }                                      << 589               <<'\t'<<anHisto->get_bin_error(int(i),int(j))*scaling_factor<<std::endl;
                                                   >> 590     }
482   }                                               591   }
483 }                                                 592 }
484 */                                             << 593 ////////////////////////////////////////////////////////////////////////////////
485 //....oooOO0OOooo........oooOO0OOooo........oo << 594 //
486                                                << 595 void RMC01AnalysisManager::ResetHistograms()
487 void RMC01AnalysisManager::ComputeMeanEdepAndE << 596 { edep_vs_prim_ekin->reset();
488                                                << 597   electron_current->reset();
489 {                                              << 598   proton_current->reset();
490   G4int nb_event = anEvent->GetEventID() + 1;  << 599   gamma_current->reset();
491   G4double factor = 1.;                        << 600   
492   if (fAdjoint_sim_mode) {                     << 601   edep_answer_matrix_vs_electron_prim_energy->reset();
493     nb_event /= fNb_evt_per_adj_evt;           << 602   electron_current_answer_matrix_vs_electron_prim_energy->reset();
494     factor = 1. * G4AdjointSimManager::GetInst << 603   gamma_current_answer_matrix_vs_electron_prim_energy->reset();
495   }                                            << 604   
496                                                << 605   edep_answer_matrix_vs_gamma_prim_energy->reset();
497   // VI: error computation now is based on num << 606   electron_current_answer_matrix_vs_gamma_prim_energy->reset();
498   //     number of events                      << 607   gamma_current_answer_matrix_vs_gamma_prim_energy->reset();
499   // LD: This is wrong! With the use of fNentr << 608   
500   //     correctly normalised. The mean and th << 609   edep_answer_matrix_vs_proton_prim_energy->reset();
501   //     with nb_event. The old computation ha << 610   electron_current_answer_matrix_vs_proton_prim_energy->reset();
502   // VI: OK, but let computations be double    << 611   proton_current_answer_matrix_vs_proton_prim_energy->reset();
503   if (nb_event > 0) {                          << 612   gamma_current_answer_matrix_vs_proton_prim_energy->reset();
504     G4double norm = 1.0 / (G4double)nb_event;  << 
505     mean = fAccumulated_edep * norm;           << 
506     G4double mean_x2 = fAccumulated_edep2 * no << 
507     G4double zz = mean_x2 - mean * mean;       << 
508     /*                                         << 
509      G4cout << "Nevt= " << nb_event <<  " mean << 
510             << "  mean_x2= " <<  mean_x2 << "  << 
511             << zz << G4endl;                   << 
512     */                                         << 
513     error = factor * std::sqrt(std::max(zz, 0. << 
514     mean *= factor;                            << 
515   }                                            << 
516   else {                                       << 
517     mean = 0;                                  << 
518     error = 0;                                 << 
519   }                                            << 
520   // G4cout << "Aend: " << mean << " " << erro << 
521 }                                                 613 }
522                                                << 614 ////////////////////////////////////////////////////////////////////////////////
523 //....oooOO0OOooo........oooOO0OOooo........oo << 615 //
524                                                << 616 void RMC01AnalysisManager::ComputeMeanEdepAndError(const G4Event* anEvent,G4double& mean,G4double& error)
525 void RMC01AnalysisManager::SetPrimaryExpSpectr << 617 {  
526                                                << 618    G4int nb_event=anEvent->GetEventID()+1;
527                                                << 619    G4double factor=1.;
528 {                                              << 620    if (adjoint_sim_mode) {
529   fPrimSpectrumType = EXPO;                    << 621     nb_event /=nb_evt_per_adj_evt;
530   if (particle_name == "e-")                   << 622   factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
531     fPrimPDG_ID = G4Electron::Electron()->GetP << 623    }
532   else if (particle_name == "gamma")           << 624    
533     fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco << 625    //error computation
534   else if (particle_name == "proton")          << 626    if (nb_event>1) {
535     fPrimPDG_ID = G4Proton::Proton()->GetPDGEn << 627       mean = accumulated_edep/nb_event;
536   else {                                       << 628     G4double mean_x2 =accumulated_edep2/nb_event;
537     G4cout << "The particle that you did selec << 629       error = factor*std::sqrt(mean_x2-mean*mean)/std::sqrt(G4double(nb_event));
538            << "list for primary [e-, gamma, pr << 630     mean *=factor;
539     return;                                    << 631    } else {
                                                   >> 632           mean=0;
                                                   >> 633           error=0;
540   }                                               634   }
541   fAlpha_or_E0 = E0;                           << 
542   fAmplitude_prim_spectrum =                   << 
543     omni_fluence / E0 / (std::exp(-Emin / E0)  << 
544   fEmin_prim_spectrum = Emin;                  << 
545   fEmax_prim_spectrum = Emax;                  << 
546 }                                                 635 }
547                                                << 636 ////////////////////////////////////////////////////////////////////////////////
548 //....oooOO0OOooo........oooOO0OOooo........oo << 637 //
549                                                << 638 void RMC01AnalysisManager::SetPrimaryExponentialSpectrumForAdjointSim(const G4String& particle_name, G4double omni_fluence,G4double E0, G4double Emin,G4double Emax)
550 void RMC01AnalysisManager::SetPrimaryPowerLawS << 639 { the_prim_spectrum_type = EXPO;
551                                                << 640   if (particle_name == "e-" ) thePrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
552                                                << 641   else if (particle_name == "gamma") thePrimPDG_ID = G4Gamma::Gamma()->GetPDGEncoding();
553                                                << 642   else if (particle_name == "proton") thePrimPDG_ID = G4Proton::Proton()->GetPDGEncoding();
554 {                                              << 
555   fPrimSpectrumType = POWER;                   << 
556   if (particle_name == "e-")                   << 
557     fPrimPDG_ID = G4Electron::Electron()->GetP << 
558   else if (particle_name == "gamma")           << 
559     fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco << 
560   else if (particle_name == "proton")          << 
561     fPrimPDG_ID = G4Proton::Proton()->GetPDGEn << 
562   else {                                       << 
563     G4cout << "The particle that you did selec << 
564            << " list for primary [e-, gamma, p << 
565     return;                                    << 
566   }                                            << 
567                                                << 
568   if (alpha == 1.) {                           << 
569     fAmplitude_prim_spectrum = omni_fluence /  << 
570   }                                            << 
571   else {                                          643   else {
572     G4double p = 1. - alpha;                   << 644     G4cout<<"The particle that you did select is not in the candidate list for primary [e-, gamma, proton]!"<<G4endl;
573     fAmplitude_prim_spectrum = omni_fluence /  << 645     return;
574   }                                            << 646   } 
575                                                << 647   
576   fAlpha_or_E0 = alpha;                        << 648 
577   fEmin_prim_spectrum = Emin;                  << 649   alpha_or_E0 = E0 ;
578   fEmax_prim_spectrum = Emax;                  << 650   amplitude_prim_spectrum = omni_fluence/E0/(std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
                                                   >> 651   emin_prim_spectrum = Emin ;
                                                   >> 652   emax_prim_spectrum = Emax;
579 }                                                 653 }
580                                                << 654 ////////////////////////////////////////////////////////////////////////////////
581 //....oooOO0OOooo........oooOO0OOooo........oo << 655 //
582                                                << 656 void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim(const G4String& particle_name, G4double omni_fluence,G4double alpha, G4double Emin,G4double Emax)
583 void RMC01AnalysisManager::Book()              << 657 { the_prim_spectrum_type  =POWER;
584 {                                              << 658   if (particle_name == "e-" ) thePrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
585   //----------------------                     << 659   else if (particle_name == "gamma") thePrimPDG_ID = G4Gamma::Gamma()->GetPDGEncoding();
586   // Creation of histograms                    << 660   else if (particle_name == "proton") thePrimPDG_ID = G4Proton::Proton()->GetPDGEncoding();
587   //----------------------                     << 661   else {
588                                                << 662     G4cout<<"The particle that you did select is not in the candidate list for primary [e-, gamma, proton]!"<<G4endl;
589   // Energy binning of the histograms : 60 log << 663     return;
590                                                << 664   } 
591   G4double emin = 1. * keV;                    << 665   
592   G4double emax = 1. * GeV;                    << 666 
593                                                << 667  if (alpha ==1.) {
594   // file_name                                 << 668   amplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
595   fFileName[0] = "forward_sim";                << 669  }
596   if (fAdjoint_sim_mode) fFileName[0] = "adjoi << 670  else {
597                                                << 671   G4double p=1.-alpha;
598   // Histo manager                             << 672   amplitude_prim_spectrum = omni_fluence/p/(std::pow(Emax,p)-std::pow(Emin,p))/4./pi;
599   G4AnalysisManager* theHistoManager = G4Analy << 673   
600   theHistoManager->SetDefaultFileType("root"); << 674  }
601   G4String extension = theHistoManager->GetFil << 675 
602   fFileName[1] = fFileName[0] + "." + extensio << 676   alpha_or_E0 = alpha ;
603   theHistoManager->SetFirstHistoId(1);         << 677   emin_prim_spectrum = Emin ;
604                                                << 678   emax_prim_spectrum = Emax;
605   G4bool fileOpen = theHistoManager->OpenFile( << 
606   if (!fileOpen) {                             << 
607     G4cout << "\n---> RMC01AnalysisManager::Bo << 
608     return;                                    << 
609   }                                            << 
610                                                << 
611   // Create directories                        << 
612   theHistoManager->SetHistoDirectoryName("hist << 
613                                                << 
614   // Histograms for :                          << 
615   //         1)the forward simulation results  << 
616   //         2)the Reverse MC  simulation resu << 
617   //------------------------------------------ << 
618                                                << 
619   G4int idHisto =                              << 
620     theHistoManager->CreateH1(G4String("Edep_v << 
621                               60, emin, emax,  << 
622   fEdep_vs_prim_ekin = theHistoManager->GetH1( << 
623                                                << 
624   idHisto = theHistoManager->CreateH1(G4String << 
625                                       emax, "n << 
626                                                << 
627   fElectron_current = theHistoManager->GetH1(i << 
628                                                << 
629   idHisto = theHistoManager->CreateH1(G4String << 
630                                       emax, "n << 
631   fProton_current = theHistoManager->GetH1(idH << 
632                                                << 
633   idHisto = theHistoManager->CreateH1(G4String << 
634                                       "none",  << 
635   fGamma_current = theHistoManager->GetH1(idHi << 
636                                                << 
637   // Response matrices for the adjoint simulat << 
638   //------------------------------------------ << 
639   if (fAdjoint_sim_mode) {                     << 
640     // Response matrices for external isotropi << 
641     //---------------------------------------- << 
642                                                << 
643     idHisto = theHistoManager->CreateH1(G4Stri << 
644                                         G4Stri << 
645                                         emax,  << 
646     fEdep_rmatrix_vs_electron_prim_energy = th << 
647                                                << 
648     idHisto = theHistoManager->CreateH2(       << 
649       G4String("Electron_current_rmatrix_vs_el << 
650       G4String("electron current  RM vs e- pri << 
651       "none", "none", "none", G4String("log"), << 
652                                                << 
653     fElectron_current_rmatrix_vs_electron_prim << 
654                                                << 
655     idHisto = theHistoManager->CreateH2(G4Stri << 
656                                         G4Stri << 
657                                         emin,  << 
658                                         G4Stri << 
659                                                << 
660     fGamma_current_rmatrix_vs_electron_prim_en << 
661                                                << 
662     // Response matrices for external isotropi << 
663                                                << 
664     idHisto = theHistoManager->CreateH1(G4Stri << 
665                                         G4Stri << 
666                                         emax,  << 
667     fEdep_rmatrix_vs_gamma_prim_energy = theHi << 
668                                                << 
669     idHisto = theHistoManager->CreateH2(G4Stri << 
670                                         G4Stri << 
671                                         60, em << 
672                                         "none" << 
673                                                << 
674     fElectron_current_rmatrix_vs_gamma_prim_en << 
675                                                << 
676     idHisto = theHistoManager->CreateH2(G4Stri << 
677                                         G4Stri << 
678                                         emin,  << 
679                                         G4Stri << 
680                                                << 
681     fGamma_current_rmatrix_vs_gamma_prim_energ << 
682                                                << 
683     // Response matrices for external isotropi << 
684     idHisto = theHistoManager->CreateH1(G4Stri << 
685                                         G4Stri << 
686                                         emax,  << 
687     fEdep_rmatrix_vs_proton_prim_energy = theH << 
688                                                << 
689     idHisto = theHistoManager->CreateH2(G4Stri << 
690                                         G4Stri << 
691                                         60, em << 
692                                         "none" << 
693                                                << 
694     fElectron_current_rmatrix_vs_proton_prim_e << 
695                                                << 
696     idHisto = theHistoManager->CreateH2(G4Stri << 
697                                         G4Stri << 
698                                         emin,  << 
699                                         G4Stri << 
700                                                << 
701     fGamma_current_rmatrix_vs_proton_prim_ener << 
702                                                << 
703     idHisto = theHistoManager->CreateH2(G4Stri << 
704                                         G4Stri << 
705                                         emin,  << 
706                                         G4Stri << 
707                                                << 
708     fProton_current_rmatrix_vs_proton_prim_ene << 
709   }                                            << 
710   fFactoryOn = true;                           << 
711   G4cout << "\n----> Histogram Tree is opened  << 
712 }                                              << 
713                                                << 
714 //....oooOO0OOooo........oooOO0OOooo........oo << 
715                                                << 
716 void RMC01AnalysisManager::Save(G4double scali << 
717 {                                              << 
718   if (fFactoryOn) {                            << 
719     G4AnalysisManager* theHistoManager = G4Ana << 
720     // scaling of results                      << 
721     //-----------------                        << 
722                                                << 
723     for (G4int ind = 1; ind <= theHistoManager << 
724       theHistoManager->SetH1Ascii(ind, true);  << 
725       theHistoManager->ScaleH1(ind, scaling_fa << 
726     }                                          << 
727     for (G4int ind = 1; ind <= theHistoManager << 
728       theHistoManager->ScaleH2(ind, scaling_fa << 
729                                                << 
730     theHistoManager->Write();                  << 
731     theHistoManager->CloseFile();              << 
732     G4cout << "\n----> Histogram Tree is saved << 
733                                                << 
734     theHistoManager->Clear();                  << 
735     fFactoryOn = false;                        << 
736   }                                            << 
737 }                                                 679 }
738                                                   680