ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RMC01AnalysisManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RMC01AnalysisManager.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 //
31 // Class Name: RMC01AnalysisManager
32 // Author: L. Desorgher
33 // Organisation: SpaceIT GmbH
34 // Contract: ESA contract 21435/08/NL/AT
35 // Customer: ESA/ESTEC
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
40 
41 #include "RMC01AnalysisManager.hh"
42 #include "G4AdjointSimManager.hh"
43 #include "G4SDManager.hh"
44 #include "RMC01SD.hh"
45 #include "G4THitsCollection.hh"
46 #include "G4Electron.hh"
47 #include "G4Proton.hh"
48 #include "G4Gamma.hh"
49 #include "G4Timer.hh"
50 #include "G4RunManager.hh"
51 #include "G4PhysicalConstants.hh"
52 #include "G4SystemOfUnits.hh"
53 #include "RMC01AnalysisManagerMessenger.hh"
54 
55 RMC01AnalysisManager* RMC01AnalysisManager::fInstance = 0;
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
59 RMC01AnalysisManager::RMC01AnalysisManager()
60  :fAccumulated_edep(0.), fAccumulated_edep2(0.), fMean_edep(0.),
61  fError_mean_edep(0.), fRelative_error(0.), fElapsed_time(0.),
62  fPrecision_to_reach(0.),fStop_run_if_precision_reached(true),
63  fNb_evt_modulo_for_convergence_test(5000),
64  fEdep_rmatrix_vs_electron_prim_energy(0),
65  fElectron_current_rmatrix_vs_electron_prim_energy(0),
66  fGamma_current_rmatrix_vs_electron_prim_energy(0),
67  fEdep_rmatrix_vs_gamma_prim_energy(0),
68  fElectron_current_rmatrix_vs_gamma_prim_energy(0),
69  fGamma_current_rmatrix_vs_gamma_prim_energy(0),
70  fEdep_rmatrix_vs_proton_prim_energy(0),
71  fElectron_current_rmatrix_vs_proton_prim_energy(0),
72  fProton_current_rmatrix_vs_proton_prim_energy(0),
73  fGamma_current_rmatrix_vs_proton_prim_energy(0),
74  fFactoryOn(false),
75  fPrimSpectrumType(EXPO),
76  fAlpha_or_E0(.5*MeV),fAmplitude_prim_spectrum (1.),
77  fEmin_prim_spectrum(1.*keV),fEmax_prim_spectrum (20.*MeV),
78  fAdjoint_sim_mode(true),fNb_evt_per_adj_evt(2)
79 {
80 
81  fMsg = new RMC01AnalysisManagerMessenger(this);
82 
83  //-------------
84  //Timer for convergence vector
85  //-------------
86 
87  fTimer = new G4Timer();
88 
89  //---------------------------------
90  //Primary particle ID for normalisation of adjoint results
91  //---------------------------------
92 
93  fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
94 
95  fFileName[0] = "sim";
96 
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
101 RMC01AnalysisManager::~RMC01AnalysisManager()
102 {;
103 }
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106 
107 RMC01AnalysisManager* RMC01AnalysisManager::GetInstance()
108 {
109  if (fInstance == 0) fInstance = new RMC01AnalysisManager;
110  return fInstance;
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 
115 void RMC01AnalysisManager::BeginOfRun(const G4Run* aRun)
116 {
117 
118  fAccumulated_edep =0.;
119  fAccumulated_edep2 =0.;
120  fNentry = 0.0;
121  fRelative_error=1.;
122  fMean_edep=0.;
123  fError_mean_edep=0.;
124  fAdjoint_sim_mode =G4AdjointSimManager::GetInstance()->GetAdjointSimMode();
125 
126  if (fAdjoint_sim_mode){
127  fNb_evt_per_adj_evt=aRun->GetNumberOfEventToBeProcessed()/
128  G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
129  fConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt",
130  std::ios::out);
131  fConvergenceFileOutput<<
132  "Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
133  }
134  else {
135  fConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt",
136  std::ios::out);
137  fConvergenceFileOutput<<
138  "Edep per event [MeV]\terror[MeV]\tcomputing_time[s]"
139  <<std::endl;
140  }
141  fConvergenceFileOutput.setf(std::ios::scientific);
142  fConvergenceFileOutput.precision(6);
143 
144  fTimer->Start();
145  fElapsed_time=0.;
146 
147  Book();
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
151 
152 void RMC01AnalysisManager::EndOfRun(const G4Run* aRun)
153 { fTimer->Stop();
154  G4int nb_evt=aRun->GetNumberOfEvent();
155  G4double factor =1./ nb_evt;
156  if (!fAdjoint_sim_mode){
157  G4cout<<"Results of forward simulation!"<<std::endl;
158  G4cout<<"edep per event [MeV] = "<<fMean_edep<<std::endl;
159  G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
160  }
161 
162  else {
163  G4cout<<"Results of reverse/adjoint simulation!"<<std::endl;
164  G4cout<<"normalised edep [MeV] = "<<fMean_edep<<std::endl;
165  G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
166  factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun()
167  *fNb_evt_per_adj_evt/aRun->GetNumberOfEvent();
168  }
169  Save(factor);
170  fConvergenceFileOutput.close();
171 }
172 
173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
174 
175 void RMC01AnalysisManager::BeginOfEvent(const G4Event* )
176 { ;
177 }
178 
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180 
181 void RMC01AnalysisManager::EndOfEvent(const G4Event* anEvent)
182 {
183  if (fAdjoint_sim_mode) EndOfEventForAdjointSimulation(anEvent);
184  else EndOfEventForForwardSimulation(anEvent);
185 
186  //Test convergence. The error is already computed
187  //--------------------------------------
188  G4int nb_event=anEvent->GetEventID()+1;
189  //G4double factor=1.;
190  if (fAdjoint_sim_mode) {
191  G4double n_adj_evt= nb_event/fNb_evt_per_adj_evt;
192  // nb_event/fNb_evt_per_adj_evt;
193  if (n_adj_evt*fNb_evt_per_adj_evt == nb_event) {
194  nb_event =static_cast<G4int>(n_adj_evt);
195  }
196  else nb_event=0;
197  }
198 
199  if (nb_event>100 && fStop_run_if_precision_reached &&
200  fPrecision_to_reach >fRelative_error) {
201  G4cout<<fPrecision_to_reach*100.<<"% Precision reached!"<<std::endl;
202  fTimer->Stop();
203  fElapsed_time+=fTimer->GetRealElapsed();
204  fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep
205  <<'\t'<<fElapsed_time<<std::endl;
206  G4RunManager::GetRunManager()->AbortRun(true);
207  }
208 
209  if (nb_event>0 && nb_event % fNb_evt_modulo_for_convergence_test == 0) {
210  fTimer->Stop();
211  fElapsed_time+=fTimer->GetRealElapsed();
212  fTimer->Start();
213  fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep<<'\t'
214  <<fElapsed_time<<std::endl;
215  }
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219 
220 void RMC01AnalysisManager::EndOfEventForForwardSimulation(
221  const G4Event* anEvent)
222 {
223 
224  G4SDManager* SDman = G4SDManager::GetSDMpointer();
225  G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
226  RMC01DoubleWithWeightHitsCollection* edepCollection =
227  (RMC01DoubleWithWeightHitsCollection*)
228  (HCE->GetHC(SDman->GetCollectionID("edep")));
229 
230  RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
231  (RMC01DoubleWithWeightHitsCollection*)
232  (HCE->GetHC(SDman->GetCollectionID("current_electron")));
233 
234  RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
235  (RMC01DoubleWithWeightHitsCollection*)
236  (HCE->GetHC(SDman->GetCollectionID("current_proton")));
237 
238  RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
239  (RMC01DoubleWithWeightHitsCollection*)
240  (HCE->GetHC(SDman->GetCollectionID("current_gamma")));
241 
242  //Total energy deposited in Event
243  //-------------------------------
244  G4double totEdep=0;
245  std::size_t i;
246  for (i=0;i<edepCollection->entries();++i)
247  totEdep+=(*edepCollection)[i]->GetValue()
248  *(*edepCollection)[i]->GetWeight();
249 
250  if (totEdep>0.){
251  fAccumulated_edep +=totEdep ;
252  fAccumulated_edep2 +=totEdep*totEdep;
253  fNentry += 1.0;
254  G4PrimaryParticle* thePrimary=
255  anEvent->GetPrimaryVertex()->GetPrimary();
256  G4double E0= thePrimary->GetG4code()->GetPDGMass();
257  G4double P=thePrimary->GetMomentum().mag();
258  G4double prim_ekin =std::sqrt(E0*E0+P*P)-E0;
259  fEdep_vs_prim_ekin->fill(prim_ekin,totEdep);
260  }
261  ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
262  if (fError_mean_edep>0) fRelative_error= fError_mean_edep/fMean_edep;
263 
264  //Particle current on sensitive cylinder
265  //-------------------------------------
266 
267  for (i=0;i<electronCurrentCollection->entries();++i) {
268  G4double ekin =(*electronCurrentCollection)[i]->GetValue();
269  G4double weight=(*electronCurrentCollection)[i]->GetWeight();
270  fElectron_current->fill(ekin,weight);
271  }
272 
273  for (i=0;i<protonCurrentCollection->entries();++i) {
274  G4double ekin =(*protonCurrentCollection)[i]->GetValue();
275  G4double weight=(*protonCurrentCollection)[i]->GetWeight();
276  fProton_current->fill(ekin,weight);
277  }
278 
279  for (i=0;i<gammaCurrentCollection->entries();++i) {
280  G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
281  G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
282  fGamma_current->fill(ekin,weight);
283  }
284 
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 
289 void RMC01AnalysisManager::EndOfEventForAdjointSimulation(
290  const G4Event* anEvent)
291 {
292  //Output from Sensitive volume computed during the forward tracking phase
293  //-----------------------------------------------------------------------
294  G4SDManager* SDman = G4SDManager::GetSDMpointer();
295  G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
296  RMC01DoubleWithWeightHitsCollection* edepCollection =
297  (RMC01DoubleWithWeightHitsCollection*)(
298  HCE->GetHC(SDman->GetCollectionID("edep")));
299 
300  RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
301  (RMC01DoubleWithWeightHitsCollection*)(
302  HCE->GetHC(SDman->GetCollectionID("current_electron")));
303 
304  RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
305  (RMC01DoubleWithWeightHitsCollection*)(
306  HCE->GetHC(SDman->GetCollectionID("current_proton")));
307 
308  RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
309  (RMC01DoubleWithWeightHitsCollection*)(
310  HCE->GetHC(SDman->GetCollectionID("current_gamma")));
311 
312  //Computation of total energy deposited in fwd tracking phase
313  //-------------------------------
314  G4double totEdep=0;
315  std::size_t i;
316  for (i=0;i<edepCollection->entries();++i)
317  totEdep+=(*edepCollection)[i]->GetValue()*
318  (*edepCollection)[i]->GetWeight();
319 
320  //Output from adjoint tracking phase
321  //----------------------------------------------------------------------------
322 
323  G4AdjointSimManager* theAdjointSimManager =
324  G4AdjointSimManager::GetInstance();
325 
326  size_t nb_adj_track =
327  theAdjointSimManager->GetNbOfAdointTracksReachingTheExternalSurface();
328  G4double total_normalised_weight = 0.;
329 
330  //We need to loop over the adjoint tracks that have reached the external
331  //surface.
332  for (std::size_t j=0;j<nb_adj_track;++j) {
333 
334  G4int pdg_nb =theAdjointSimManager
335  ->GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack(j);
336  G4double prim_ekin=theAdjointSimManager
337  ->GetEkinAtEndOfLastAdjointTrack(j);
338  G4double adj_weight=theAdjointSimManager
339  ->GetWeightAtEndOfLastAdjointTrack(j);
340 
341 
342  //Factor of normalisation to user defined prim spectrum (power law or exp)
343  //------------------------------------------------------------------------
344  G4double normalised_weight = 0.;
345  if (pdg_nb== fPrimPDG_ID && prim_ekin>= fEmin_prim_spectrum
346  && prim_ekin<= fEmax_prim_spectrum)
347  normalised_weight =
348  adj_weight*PrimDiffAndDirFluxForAdjointSim(prim_ekin);
349  total_normalised_weight += normalised_weight;
350 
351  //Answer matrices
352  //-------------
353  G4H1* edep_rmatrix =0;
354  G4H2* electron_current_rmatrix =0;
355  G4H2* gamma_current_rmatrix =0;
356  G4H2* proton_current_rmatrix =0;
357 
358  if (pdg_nb == G4Electron::Electron()->GetPDGEncoding()){ //e- matrices
359  edep_rmatrix = fEdep_rmatrix_vs_electron_prim_energy;
360  electron_current_rmatrix =
361  fElectron_current_rmatrix_vs_electron_prim_energy;
362  gamma_current_rmatrix = fGamma_current_rmatrix_vs_electron_prim_energy;
363  }
364  else if (pdg_nb == G4Gamma::Gamma()->GetPDGEncoding()){
365  //gammma answer matrices
366  edep_rmatrix = fEdep_rmatrix_vs_gamma_prim_energy;
367  electron_current_rmatrix = fElectron_current_rmatrix_vs_gamma_prim_energy;
368  gamma_current_rmatrix = fGamma_current_rmatrix_vs_gamma_prim_energy;
369  }
370  else if (pdg_nb == G4Proton::Proton()->GetPDGEncoding()){
371  //proton answer matrices
372  edep_rmatrix = fEdep_rmatrix_vs_proton_prim_energy;
373  electron_current_rmatrix =
374  fElectron_current_rmatrix_vs_proton_prim_energy;
375  gamma_current_rmatrix = fGamma_current_rmatrix_vs_proton_prim_energy;
376  proton_current_rmatrix = fProton_current_rmatrix_vs_proton_prim_energy;
377  }
378  //Register histo edep vs prim ekin
379  //----------------------------------
380  if (normalised_weight>0) fEdep_vs_prim_ekin
381  ->fill(prim_ekin,totEdep*normalised_weight);
382  // Registering answer matrix
383  //---------------------------
384  edep_rmatrix->fill(prim_ekin,totEdep*adj_weight/cm2);
385 
386  //Registering of current of particles on the sensitive volume
387  //------------------------------------------------------------
388 
389  for (i=0;i<electronCurrentCollection->entries();++i) {
390  G4double ekin =(*electronCurrentCollection)[i]->GetValue();
391  G4double weight=(*electronCurrentCollection)[i]->GetWeight();
392  fElectron_current->fill(ekin,weight*normalised_weight);
393  electron_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
394  }
395  for (i=0;i<protonCurrentCollection->entries();++i) {
396  G4double ekin =(*protonCurrentCollection)[i]->GetValue();
397  G4double weight=(*protonCurrentCollection)[i]->GetWeight();
398  fProton_current->fill(ekin,weight*normalised_weight);
399  proton_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
400  }
401  for (i=0;i<gammaCurrentCollection->entries();++i) {
402  G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
403  G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
404  fGamma_current->fill(ekin,weight*normalised_weight);
405  gamma_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
406  }
407  }
408 
409  //Registering of total energy deposited in Event
410  //-------------------------------
411  G4bool new_mean_computed=false;
412  if (totEdep>0.){
413  if (total_normalised_weight>0.){
414  G4double edep=totEdep* total_normalised_weight;
415 
416  //Check if the edep is not wrongly too high
417  //-----------------------------------------
418  G4double new_mean(0.0), new_error(0.0);
419  fAccumulated_edep +=edep;
420  fAccumulated_edep2 +=edep*edep;
421  fNentry += 1.0;
422  ComputeMeanEdepAndError(anEvent,new_mean,new_error);
423  G4double new_relative_error = 1.;
424  if ( new_error >0) new_relative_error = new_error/ new_mean;
425  if (fRelative_error <0.10 && new_relative_error>1.5*fRelative_error) {
426  G4cout<<"Potential wrong adjoint weight!"<<std::endl;
427  G4cout<<"The results of this event will not be registered!"
428  <<std::endl;
429  G4cout<<"previous mean edep [MeV] "<< fMean_edep<<std::endl;
430  G4cout<<"previous relative error "<< fRelative_error<<std::endl;
431  G4cout<<"new rejected mean edep [MeV] "<< new_mean<<std::endl;
432  G4cout<<"new rejected relative error "<< new_relative_error
433  <<std::endl;
434  fAccumulated_edep -=edep;
435  fAccumulated_edep2 -=edep*edep;
436  fNentry -= 1.0;
437  return;
438  }
439  else { //accepted
440  fMean_edep = new_mean;
441  fError_mean_edep = new_error;
442  fRelative_error =new_relative_error;
443  new_mean_computed=true;
444  }
445 
446  }
447 
448  if (!new_mean_computed){
449  ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
450  fRelative_error = (fMean_edep > 0.0) ? fError_mean_edep/fMean_edep : 0.0;
451  }
452  }
453 }
454 
455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456 
457 G4double RMC01AnalysisManager::PrimDiffAndDirFluxForAdjointSim(
458  G4double prim_energy)
459 {
460  G4double flux=fAmplitude_prim_spectrum;
461  if ( fPrimSpectrumType ==EXPO) flux*=std::exp(-prim_energy/fAlpha_or_E0);
462  else flux*=std::pow(prim_energy, -fAlpha_or_E0);
463  return flux;
464 }
465 
466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
467 /*
468 void RMC01AnalysisManager::WriteHisto(G4H1* anHisto,
469  G4double scaling_factor, G4String fileName, G4String header_lines)
470 { std::fstream FileOutput(fileName, std::ios::out);
471  FileOutput<<header_lines;
472  FileOutput.setf(std::ios::scientific);
473  FileOutput.precision(6);
474 
475  for (G4int i =0;i<G4int(anHisto->axis().bins());++i) {
476  FileOutput<<anHisto->axis().bin_lower_edge(i)
477  <<'\t'<<anHisto->axis().bin_upper_edge(i)
478  <<'\t'<<anHisto->bin_height(i)*scaling_factor
479  <<'\t'<<anHisto->bin_error(i)*scaling_factor<<std::endl;
480  }
481 }
482 
483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
484 
485 void RMC01AnalysisManager::WriteHisto(G4H2* anHisto,
486  G4double scaling_factor, G4String fileName, G4String header_lines)
487 { std::fstream FileOutput(fileName, std::ios::out);
488  FileOutput<<header_lines;
489 
490  FileOutput.setf(std::ios::scientific);
491  FileOutput.precision(6);
492 
493  for (G4int i =0;i<G4int(anHisto->axis_x().bins());++i) {
494  for (G4int j =0;j<G4int(anHisto->axis_y().bins());++j) {
495  FileOutput<<anHisto->axis_x().bin_lower_edge(i)
496  <<'\t'<<anHisto->axis_x().bin_upper_edge(i)
497  <<'\t'<<anHisto->axis_y().bin_lower_edge(i)
498  <<'\t'<<anHisto->axis_y().bin_upper_edge(i)
499  <<'\t'<<anHisto->bin_height(i,j)*scaling_factor
500  <<'\t'<<anHisto->bin_error(i,j)*scaling_factor
501  <<std::endl;
502  }
503  }
504 }
505 */
506 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
507 
508 void RMC01AnalysisManager::ComputeMeanEdepAndError(
509  const G4Event* anEvent,G4double& mean,G4double& error)
510 {
511  G4int nb_event = anEvent->GetEventID()+1;
512  G4double factor=1.;
513  if (fAdjoint_sim_mode) {
514  nb_event /= fNb_evt_per_adj_evt;
515  factor = 1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
516  }
517 
518  // VI: error computation now is based on number of entries and not
519  // number of events
520  // LD: This is wrong! With the use of fNentry the results were no longer
521  // correctly normalised. The mean and the error should be computed
522  // with nb_event. The old computation has been reset.
523  // VI: OK, but let computations be double
524  if (nb_event > 0) {
525  G4double norm = 1.0/(G4double)nb_event;
526  mean = fAccumulated_edep*norm;
527  G4double mean_x2 = fAccumulated_edep2*norm;
528  G4double zz = mean_x2 - mean*mean;
529  /*
530  G4cout << "Nevt= " << nb_event << " mean= " << mean
531  << " mean_x2= " << mean_x2 << " x2 - x*x= "
532  << zz << G4endl;
533  */
534  error = factor*std::sqrt(std::max(zz, 0.)*norm);
535  mean *= factor;
536  } else {
537  mean=0;
538  error=0;
539  }
540  //G4cout << "Aend: " << mean << " " << error << G4endl;
541 }
542 
543 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
544 
545 void RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim(
546  const G4String& particle_name, G4double omni_fluence,
547  G4double E0, G4double Emin,
548  G4double Emax)
549 { fPrimSpectrumType = EXPO;
550  if (particle_name == "e-" ) fPrimPDG_ID =
551  G4Electron::Electron()->GetPDGEncoding();
552  else if (particle_name == "gamma") fPrimPDG_ID =
553  G4Gamma::Gamma()->GetPDGEncoding();
554  else if (particle_name == "proton") fPrimPDG_ID =
555  G4Proton::Proton()->GetPDGEncoding();
556  else {
557  G4cout<<"The particle that you did select is not in the candidate "<<
558  "list for primary [e-, gamma, proton]!"<<G4endl;
559  return;
560  }
561  fAlpha_or_E0 = E0 ;
562  fAmplitude_prim_spectrum = omni_fluence/E0/
563  (std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
564  fEmin_prim_spectrum = Emin ;
565  fEmax_prim_spectrum = Emax;
566 }
567 
568 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
569 
570 void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim(
571  const G4String& particle_name, G4double omni_fluence,
572  G4double alpha, G4double Emin,G4double Emax)
573 { fPrimSpectrumType =POWER;
574  if (particle_name == "e-" ) fPrimPDG_ID =
575  G4Electron::Electron()->GetPDGEncoding();
576  else if (particle_name == "gamma") fPrimPDG_ID =
577  G4Gamma::Gamma()->GetPDGEncoding();
578  else if (particle_name == "proton") fPrimPDG_ID =
579  G4Proton::Proton()->GetPDGEncoding();
580  else {
581  G4cout<<"The particle that you did select is not in the candidate"<<
582  " list for primary [e-, gamma, proton]!"<<G4endl;
583  return;
584  }
585 
586 
587  if (alpha ==1.) {
588  fAmplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
589  }
590  else {
591  G4double p=1.-alpha;
592  fAmplitude_prim_spectrum = omni_fluence/p/(std::pow(Emax,p)
593  -std::pow(Emin,p))/4./pi;
594  }
595 
596  fAlpha_or_E0 = alpha ;
597  fEmin_prim_spectrum = Emin ;
598  fEmax_prim_spectrum = Emax;
599 
600 }
601 
602 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
603 
604 void RMC01AnalysisManager::Book()
605 {
606  //----------------------
607  //Creation of histograms
608  //----------------------
609 
610  //Energy binning of the histograms : 60 log bins over [1keV-1GeV]
611 
612  G4double emin=1.*keV;
613  G4double emax=1.*GeV;
614 
615  //file_name
616  fFileName[0]="forward_sim";
617  if (fAdjoint_sim_mode) fFileName[0]="adjoint_sim";
618 
619  //Histo manager
620  G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
621  G4String extension = theHistoManager->GetFileType();
622  fFileName[1] = fFileName[0] + "." + extension;
623  theHistoManager->SetFirstHistoId(1);
624 
625  G4bool fileOpen = theHistoManager->OpenFile(fFileName[0]);
626  if (!fileOpen) {
627  G4cout << "\n---> RMC01AnalysisManager::Book(): cannot open "
628  << fFileName[1]
629  << G4endl;
630  return;
631  }
632 
633  // Create directories
634  theHistoManager->SetHistoDirectoryName("histo");
635 
636  //Histograms for :
637  // 1)the forward simulation results
638  // 2)the Reverse MC simulation results normalised to a user spectrum
639  //------------------------------------------------------------------------
640 
641  G4int idHisto =
642  theHistoManager->CreateH1(G4String("Edep_vs_prim_ekin"),
643  G4String("edep vs e- primary energy"),60,emin,emax,
644  "none","none",G4String("log"));
645  fEdep_vs_prim_ekin = theHistoManager->GetH1(idHisto);
646 
647  idHisto = theHistoManager->CreateH1(G4String("elecron_current"),
648  G4String("electron"),60,emin,emax,
649  "none","none",G4String("log"));
650 
651  fElectron_current = theHistoManager->GetH1(idHisto);
652 
653  idHisto= theHistoManager->CreateH1(G4String("proton_current"),
654  G4String("proton"),60,emin,emax,
655  "none","none",G4String("log"));
656  fProton_current=theHistoManager->GetH1(idHisto);
657 
658  idHisto= theHistoManager->CreateH1(G4String("gamma_current"),
659  G4String("gamma"),60,emin,emax,
660  "none","none",G4String("log"));
661  fGamma_current=theHistoManager->GetH1(idHisto);
662 
663  //Response matrices for the adjoint simulation only
664  //-----------------------------------------------
665  if (fAdjoint_sim_mode){
666  //Response matrices for external isotropic e- source
667  //--------------------------------------------------
668 
669  idHisto =
670  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_electron_prim_energy"),
671  G4String("electron RM vs e- primary energy"),60,emin,emax,
672  "none","none",G4String("log"));
673  fEdep_rmatrix_vs_electron_prim_energy = theHistoManager->GetH1(idHisto);
674 
675  idHisto =
676  theHistoManager->
677  CreateH2(G4String("Electron_current_rmatrix_vs_electron_prim_energy"),
678  G4String("electron current RM vs e- primary energy"),
679  60,emin,emax,60,emin,emax,
680  "none","none","none","none",G4String("log"),G4String("log"));
681 
682  fElectron_current_rmatrix_vs_electron_prim_energy =
683  theHistoManager->GetH2(idHisto);
684 
685  idHisto =
686  theHistoManager->
687  CreateH2(G4String("Gamma_current_rmatrix_vs_electron_prim_energy"),
688  G4String("gamma current RM vs e- primary energy"),
689  60,emin,emax,60,emin,emax,
690  "none","none","none","none",G4String("log"),G4String("log"));
691 
692  fGamma_current_rmatrix_vs_electron_prim_energy =
693  theHistoManager->GetH2(idHisto);
694 
695  //Response matrices for external isotropic gamma source
696 
697  idHisto =
698  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_gamma_prim_energy"),
699  G4String("electron RM vs gamma primary energy"),60,emin,emax,
700  "none","none",G4String("log"));
701  fEdep_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH1(idHisto);
702 
703  idHisto =
704  theHistoManager->
705  CreateH2(G4String("Electron_current_rmatrix_vs_gamma_prim_energy"),
706  G4String("electron current RM vs gamma primary energy"),
707  60,emin,emax,60,emin,emax,
708  "none","none","none","none",G4String("log"),G4String("log"));
709 
710  fElectron_current_rmatrix_vs_gamma_prim_energy =
711  theHistoManager->GetH2(idHisto);
712 
713  idHisto =
714  theHistoManager->
715  CreateH2(G4String("Gamma_current_rmatrix_vs_gamma_prim_energy"),
716  G4String("gamma current RM vs gamma primary energy"),
717  60,emin,emax,60,emin,emax,
718  "none","none","none","none",G4String("log"),G4String("log"));
719 
720  fGamma_current_rmatrix_vs_gamma_prim_energy =
721  theHistoManager->GetH2(idHisto);
722 
723  //Response matrices for external isotropic proton source
724  idHisto =
725  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_proton_prim_energy"),
726  G4String("electron RM vs proton primary energy"),60,emin,emax,
727  "none","none",G4String("log"));
728  fEdep_rmatrix_vs_proton_prim_energy = theHistoManager->GetH1(idHisto);
729 
730  idHisto =
731  theHistoManager->
732  CreateH2(G4String("Electron_current_rmatrix_vs_proton_prim_energy"),
733  G4String("electron current RM vs proton primary energy"),
734  60,emin,emax,60,emin,emax,
735  "none","none","none","none",G4String("log"),G4String("log"));
736 
737  fElectron_current_rmatrix_vs_proton_prim_energy =
738  theHistoManager->GetH2(idHisto);
739 
740  idHisto =
741  theHistoManager->
742  CreateH2(G4String("Gamma_current_rmatrix_vs_proton_prim_energy"),
743  G4String("gamma current RM vs proton primary energy"),
744  60,emin,emax,60,emin,emax,
745  "none","none","none","none",G4String("log"),G4String("log"));
746 
747  fGamma_current_rmatrix_vs_proton_prim_energy =
748  theHistoManager->GetH2(idHisto);
749 
750  idHisto =
751  theHistoManager->
752  CreateH2(G4String("Proton_current_rmatrix_vs_proton_prim_energy"),
753  G4String("proton current RM vs proton primary energy"),
754  60,emin,emax,60,emin,emax,
755  "none","none","none","none",G4String("log"),G4String("log"));
756 
757  fProton_current_rmatrix_vs_proton_prim_energy =
758  theHistoManager->GetH2(idHisto);
759  }
760  fFactoryOn = true;
761  G4cout << "\n----> Histogram Tree is opened in " << fFileName[1] << G4endl;
762 }
763 
764 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
765 
766 void RMC01AnalysisManager::Save(G4double scaling_factor)
767 { if (fFactoryOn) {
768  G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
769  //scaling of results
770  //-----------------
771 
772  for (G4int ind=1; ind<=theHistoManager->GetNofH1s();++ind){
773  theHistoManager->SetH1Ascii(ind,true);
774  theHistoManager->ScaleH1(ind,scaling_factor);
775  }
776  for (G4int ind=1; ind<=theHistoManager->GetNofH2s();++ind)
777  theHistoManager->ScaleH2(ind,scaling_factor);
778 
779  theHistoManager->Write();
780  theHistoManager->CloseFile();
781  G4cout << "\n----> Histogram Tree is saved in " << fFileName[1] << G4endl;
782 
783  delete G4AnalysisManager::Instance();
784  fFactoryOn = false;
785  }
786 }