Geant4  10.03
RMC01AnalysisManager.cc
Go to the documentation of this file.
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 // $Id: RMC01AnalysisManager.cc 101420 2016-11-17 10:37:14Z gcosmo $
30 //
32 // Class Name: RMC01AnalysisManager
33 // Author: L. Desorgher
34 // Organisation: SpaceIT GmbH
35 // Contract: ESA contract 21435/08/NL/AT
36 // Customer: ESA/ESTEC
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41 
42 #include "RMC01AnalysisManager.hh"
43 #include "G4AdjointSimManager.hh"
44 #include "G4SDManager.hh"
45 #include "RMC01SD.hh"
46 #include "G4THitsCollection.hh"
47 #include "G4Electron.hh"
48 #include "G4Proton.hh"
49 #include "G4Gamma.hh"
50 #include "G4Timer.hh"
51 #include "G4RunManager.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "RMC01AnalysisManagerMessenger.hh"
55 
56 RMC01AnalysisManager* RMC01AnalysisManager::fInstance = 0;
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 
60 RMC01AnalysisManager::RMC01AnalysisManager()
61  :fAccumulated_edep(0.), fAccumulated_edep2(0.), fMean_edep(0.),
62  fError_mean_edep(0.), fRelative_error(0.), fElapsed_time(0.),
63  fPrecision_to_reach(0.),fStop_run_if_precision_reached(true),
64  fNb_evt_modulo_for_convergence_test(5000),
65  fEdep_rmatrix_vs_electron_prim_energy(0),
66  fElectron_current_rmatrix_vs_electron_prim_energy(0),
67  fGamma_current_rmatrix_vs_electron_prim_energy(0),
68  fEdep_rmatrix_vs_gamma_prim_energy(0),
69  fElectron_current_rmatrix_vs_gamma_prim_energy(0),
70  fGamma_current_rmatrix_vs_gamma_prim_energy(0),
71  fEdep_rmatrix_vs_proton_prim_energy(0),
72  fElectron_current_rmatrix_vs_proton_prim_energy(0),
73  fProton_current_rmatrix_vs_proton_prim_energy(0),
74  fGamma_current_rmatrix_vs_proton_prim_energy(0),
75  fFactoryOn(false),
76  fPrimSpectrumType(EXPO),
77  fAlpha_or_E0(.5*MeV),fAmplitude_prim_spectrum (1.),
78  fEmin_prim_spectrum(1.*keV),fEmax_prim_spectrum (20.*MeV),
79  fAdjoint_sim_mode(true),fNb_evt_per_adj_evt(2)
80 {
81 
82  fMsg = new RMC01AnalysisManagerMessenger(this);
83 
84  //-------------
85  //Timer for convergence vector
86  //-------------
87 
88  fTimer = new G4Timer();
89 
90  //---------------------------------
91  //Primary particle ID for normalisation of adjoint results
92  //---------------------------------
93 
94  fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
95 
96  fFileName[0] = "sim";
97 
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 
102 RMC01AnalysisManager::~RMC01AnalysisManager()
103 {;
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 
108 RMC01AnalysisManager* RMC01AnalysisManager::GetInstance()
109 {
110  if (fInstance == 0) fInstance = new RMC01AnalysisManager;
111  return fInstance;
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115 
116 void RMC01AnalysisManager::BeginOfRun(const G4Run* aRun)
117 {
118 
119  fAccumulated_edep =0.;
120  fAccumulated_edep2 =0.;
121  fNentry = 0.0;
122  fRelative_error=1.;
123  fMean_edep=0.;
124  fError_mean_edep=0.;
125  fAdjoint_sim_mode =G4AdjointSimManager::GetInstance()->GetAdjointSimMode();
126 
127  if (fAdjoint_sim_mode){
128  fNb_evt_per_adj_evt=aRun->GetNumberOfEventToBeProcessed()/
129  G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
130  fConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt",
131  std::ios::out);
132  fConvergenceFileOutput<<
133  "Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
134  }
135  else {
136  fConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt",
137  std::ios::out);
138  fConvergenceFileOutput<<
139  "Edep per event [MeV]\terror[MeV]\tcomputing_time[s]"
140  <<std::endl;
141  }
142  fConvergenceFileOutput.setf(std::ios::scientific);
143  fConvergenceFileOutput.precision(6);
144 
145  fTimer->Start();
146  fElapsed_time=0.;
147 
148  Book();
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152 
153 void RMC01AnalysisManager::EndOfRun(const G4Run* aRun)
154 { fTimer->Stop();
155  G4int nb_evt=aRun->GetNumberOfEvent();
156  G4double factor =1./ nb_evt;
157  if (!fAdjoint_sim_mode){
158  G4cout<<"Results of forward simulation!"<<std::endl;
159  G4cout<<"edep per event [MeV] = "<<fMean_edep<<std::endl;
160  G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
161  }
162 
163  else {
164  G4cout<<"Results of reverse/adjoint simulation!"<<std::endl;
165  G4cout<<"normalised edep [MeV] = "<<fMean_edep<<std::endl;
166  G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
167  factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun()
168  *fNb_evt_per_adj_evt/aRun->GetNumberOfEvent();
169  }
170  Save(factor);
171  fConvergenceFileOutput.close();
172 }
173 
174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
175 
176 void RMC01AnalysisManager::BeginOfEvent(const G4Event* )
177 { ;
178 }
179 
180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
181 
182 void RMC01AnalysisManager::EndOfEvent(const G4Event* anEvent)
183 {
184  if (fAdjoint_sim_mode) EndOfEventForAdjointSimulation(anEvent);
185  else EndOfEventForForwardSimulation(anEvent);
186 
187  //Test convergence. The error is already computed
188  //--------------------------------------
189  G4int nb_event=anEvent->GetEventID()+1;
190  //G4double factor=1.;
191  if (fAdjoint_sim_mode) {
192  G4double n_adj_evt= nb_event/fNb_evt_per_adj_evt;
193  // nb_event/fNb_evt_per_adj_evt;
194  if (n_adj_evt*fNb_evt_per_adj_evt == nb_event) {
195  nb_event =static_cast<G4int>(n_adj_evt);
196  }
197  else nb_event=0;
198  }
199 
200  if (nb_event>100 && fStop_run_if_precision_reached &&
201  fPrecision_to_reach >fRelative_error) {
202  G4cout<<fPrecision_to_reach*100.<<"% Precision reached!"<<std::endl;
203  fTimer->Stop();
204  fElapsed_time+=fTimer->GetRealElapsed();
205  fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep
206  <<'\t'<<fElapsed_time<<std::endl;
207  G4RunManager::GetRunManager()->AbortRun(true);
208  }
209 
210  if (nb_event>0 && nb_event % fNb_evt_modulo_for_convergence_test == 0) {
211  fTimer->Stop();
212  fElapsed_time+=fTimer->GetRealElapsed();
213  fTimer->Start();
214  fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep<<'\t'
215  <<fElapsed_time<<std::endl;
216  }
217 }
218 
219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
220 
221 void RMC01AnalysisManager::EndOfEventForForwardSimulation(
222  const G4Event* anEvent)
223 {
224 
225  G4SDManager* SDman = G4SDManager::GetSDMpointer();
226  G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
227  RMC01DoubleWithWeightHitsCollection* edepCollection =
228  (RMC01DoubleWithWeightHitsCollection*)
229  (HCE->GetHC(SDman->GetCollectionID("edep")));
230 
231  RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
232  (RMC01DoubleWithWeightHitsCollection*)
233  (HCE->GetHC(SDman->GetCollectionID("current_electron")));
234 
235  RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
236  (RMC01DoubleWithWeightHitsCollection*)
237  (HCE->GetHC(SDman->GetCollectionID("current_proton")));
238 
239  RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
240  (RMC01DoubleWithWeightHitsCollection*)
241  (HCE->GetHC(SDman->GetCollectionID("current_gamma")));
242 
243  //Total energy deposited in Event
244  //-------------------------------
245  G4double totEdep=0;
246  G4int i;
247  for (i=0;i<edepCollection->entries();i++)
248  totEdep+=(*edepCollection)[i]->GetValue()
249  *(*edepCollection)[i]->GetWeight();
250 
251  if (totEdep>0.){
252  fAccumulated_edep +=totEdep ;
253  fAccumulated_edep2 +=totEdep*totEdep;
254  fNentry += 1.0;
255  G4PrimaryParticle* thePrimary=
256  anEvent->GetPrimaryVertex()->GetPrimary();
257  G4double E0= thePrimary->GetG4code()->GetPDGMass();
258  G4double P=thePrimary->GetMomentum().mag();
259  G4double prim_ekin =std::sqrt(E0*E0+P*P)-E0;
260  fEdep_vs_prim_ekin->fill(prim_ekin,totEdep);
261  }
262  ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
263  if (fError_mean_edep>0) fRelative_error= fError_mean_edep/fMean_edep;
264 
265  //Particle current on sensitive cylinder
266  //-------------------------------------
267 
268  for (i=0;i<electronCurrentCollection->entries();i++) {
269  G4double ekin =(*electronCurrentCollection)[i]->GetValue();
270  G4double weight=(*electronCurrentCollection)[i]->GetWeight();
271  fElectron_current->fill(ekin,weight);
272  }
273 
274  for (i=0;i<protonCurrentCollection->entries();i++) {
275  G4double ekin =(*protonCurrentCollection)[i]->GetValue();
276  G4double weight=(*protonCurrentCollection)[i]->GetWeight();
277  fProton_current->fill(ekin,weight);
278  }
279 
280  for (i=0;i<gammaCurrentCollection->entries();i++) {
281  G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
282  G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
283  fGamma_current->fill(ekin,weight);
284  }
285 
286 }
287 
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 
290 void RMC01AnalysisManager::EndOfEventForAdjointSimulation(
291  const G4Event* anEvent)
292 {
293  //Output from Sensitive volume computed during the forward tracking phase
294  //-----------------------------------------------------------------------
295  G4SDManager* SDman = G4SDManager::GetSDMpointer();
296  G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
297  RMC01DoubleWithWeightHitsCollection* edepCollection =
298  (RMC01DoubleWithWeightHitsCollection*)(
299  HCE->GetHC(SDman->GetCollectionID("edep")));
300 
301  RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
302  (RMC01DoubleWithWeightHitsCollection*)(
303  HCE->GetHC(SDman->GetCollectionID("current_electron")));
304 
305  RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
306  (RMC01DoubleWithWeightHitsCollection*)(
307  HCE->GetHC(SDman->GetCollectionID("current_proton")));
308 
309  RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
310  (RMC01DoubleWithWeightHitsCollection*)(
311  HCE->GetHC(SDman->GetCollectionID("current_gamma")));
312 
313  //Computation of total energy deposited in fwd tracking phase
314  //-------------------------------
315  G4double totEdep=0;
316  G4int i;
317  for (i=0;i<edepCollection->entries();i++)
318  totEdep+=(*edepCollection)[i]->GetValue()*
319  (*edepCollection)[i]->GetWeight();
320 
321  //Output from adjoint tracking phase
322  //----------------------------------------------------------------------------
323 
324  G4AdjointSimManager* theAdjointSimManager =
325  G4AdjointSimManager::GetInstance();
326 
327  size_t nb_adj_track =
328  theAdjointSimManager->GetNbOfAdointTracksReachingTheExternalSurface();
329  G4double total_normalised_weight = 0.;
330 
331  //We need to loop over the adjoint tracks that have reached the external
332  //surface.
333  for (size_t j=0;j<nb_adj_track;j++) {
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  G4AnaH1* edep_rmatrix =0;
354  G4AnaH2* electron_current_rmatrix =0;
355  G4AnaH2* gamma_current_rmatrix =0;
356  G4AnaH2* 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 , new_error;
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  if (!new_mean_computed){
448  ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
449  if (fError_mean_edep>0) fRelative_error= fError_mean_edep/fMean_edep;
450  }
451  }
452 }
453 
454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
455 
456 G4double RMC01AnalysisManager::PrimDiffAndDirFluxForAdjointSim(
457  G4double prim_energy)
458 {
459  G4double flux=fAmplitude_prim_spectrum;
460  if ( fPrimSpectrumType ==EXPO) flux*=std::exp(-prim_energy/fAlpha_or_E0);
461  else flux*=std::pow(prim_energy, -fAlpha_or_E0);
462  return flux;
463 }
464 
465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
466 /*
467 void RMC01AnalysisManager::WriteHisto(G4AnaH1* anHisto,
468  G4double scaling_factor, G4String fileName, G4String header_lines)
469 { std::fstream FileOutput(fileName, std::ios::out);
470  FileOutput<<header_lines;
471  FileOutput.setf(std::ios::scientific);
472  FileOutput.precision(6);
473 
474  for (G4int i =0;i<G4int(anHisto->axis().bins());i++) {
475  FileOutput<<anHisto->axis().bin_lower_edge(i)
476  <<'\t'<<anHisto->axis().bin_upper_edge(i)
477  <<'\t'<<anHisto->bin_height(i)*scaling_factor
478  <<'\t'<<anHisto->bin_error(i)*scaling_factor<<std::endl;
479  }
480 }
481 
482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
483 
484 void RMC01AnalysisManager::WriteHisto(G4AnaH2* anHisto,
485  G4double scaling_factor, G4String fileName, G4String header_lines)
486 { std::fstream FileOutput(fileName, std::ios::out);
487  FileOutput<<header_lines;
488 
489  FileOutput.setf(std::ios::scientific);
490  FileOutput.precision(6);
491 
492  for (G4int i =0;i<G4int(anHisto->axis_x().bins());i++) {
493  for (G4int j =0;j<G4int(anHisto->axis_y().bins());j++) {
494  FileOutput<<anHisto->axis_x().bin_lower_edge(i)
495  <<'\t'<<anHisto->axis_x().bin_upper_edge(i)
496  <<'\t'<<anHisto->axis_y().bin_lower_edge(i)
497  <<'\t'<<anHisto->axis_y().bin_upper_edge(i)
498  <<'\t'<<anHisto->bin_height(i,j)*scaling_factor
499  <<'\t'<<anHisto->bin_error(i,j)*scaling_factor
500  <<std::endl;
501  }
502  }
503 }
504 */
505 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
506 
507 void RMC01AnalysisManager::ComputeMeanEdepAndError(
508  const G4Event* anEvent,G4double& mean,G4double& error)
509 {
510  G4int nb_event=anEvent->GetEventID()+1;
511  G4double factor=1.;
512  if (fAdjoint_sim_mode) {
513  nb_event /=fNb_evt_per_adj_evt;
514  factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
515  }
516 
517  // VI: error computation now is based on number of entries and not
518  // number of events
519  // LD: This is wrong! With the use of fNentry the results were no longer
520  // correctly normalised. The mean and the error should be computed
521  // with nb_event. The old computation has been reset.
522  G4float nb_event_float = G4float(nb_event);
523  if (nb_event_float >1.) {
524  mean = fAccumulated_edep/nb_event_float;
525  G4double mean_x2 = fAccumulated_edep2/nb_event_float;
526  /*
527  G4cout << "Nevt= " << nb_event << " mean= " << mean
528  << " mean_x2= " << mean_x2 << " x2 - x*x= "
529  << mean_x2-mean*mean << G4endl;
530  */
531  error = factor*std::sqrt(mean_x2-mean*mean)/std::sqrt(nb_event_float);
532  mean *=factor;
533  }
534  else {
535  mean=0;
536  error=0;
537  }
538 }
539 
540 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
541 
542 void RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim(
543  const G4String& particle_name, G4double omni_fluence,
544  G4double E0, G4double Emin,
545  G4double Emax)
546 { fPrimSpectrumType = EXPO;
547  if (particle_name == "e-" ) fPrimPDG_ID =
548  G4Electron::Electron()->GetPDGEncoding();
549  else if (particle_name == "gamma") fPrimPDG_ID =
550  G4Gamma::Gamma()->GetPDGEncoding();
551  else if (particle_name == "proton") fPrimPDG_ID =
552  G4Proton::Proton()->GetPDGEncoding();
553  else {
554  G4cout<<"The particle that you did select is not in the candidate "<<
555  "list for primary [e-, gamma, proton]!"<<G4endl;
556  return;
557  }
558  fAlpha_or_E0 = E0 ;
559  fAmplitude_prim_spectrum = omni_fluence/E0/
560  (std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
561  fEmin_prim_spectrum = Emin ;
562  fEmax_prim_spectrum = Emax;
563 }
564 
565 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
566 
567 void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim(
568  const G4String& particle_name, G4double omni_fluence,
569  G4double alpha, G4double Emin,G4double Emax)
570 { fPrimSpectrumType =POWER;
571  if (particle_name == "e-" ) fPrimPDG_ID =
572  G4Electron::Electron()->GetPDGEncoding();
573  else if (particle_name == "gamma") fPrimPDG_ID =
574  G4Gamma::Gamma()->GetPDGEncoding();
575  else if (particle_name == "proton") fPrimPDG_ID =
576  G4Proton::Proton()->GetPDGEncoding();
577  else {
578  G4cout<<"The particle that you did select is not in the candidate"<<
579  " list for primary [e-, gamma, proton]!"<<G4endl;
580  return;
581  }
582 
583 
584  if (alpha ==1.) {
585  fAmplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
586  }
587  else {
588  G4double p=1.-alpha;
589  fAmplitude_prim_spectrum = omni_fluence/p/(std::pow(Emax,p)
590  -std::pow(Emin,p))/4./pi;
591  }
592 
593  fAlpha_or_E0 = alpha ;
594  fEmin_prim_spectrum = Emin ;
595  fEmax_prim_spectrum = Emax;
596 
597 }
598 
599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
600 
601 void RMC01AnalysisManager::Book()
602 {
603  //----------------------
604  //Creation of histograms
605  //----------------------
606 
607  //Energy binning of the histograms : 60 log bins over [1keV-1GeV]
608 
609  G4double emin=1.*keV;
610  G4double emax=1.*GeV;
611 
612  //file_name
613  fFileName[0]="forward_sim";
614  if (fAdjoint_sim_mode) fFileName[0]="adjoint_sim";
615 
616  //Histo manager
617  G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
618  G4String extension = theHistoManager->GetFileType();
619  fFileName[1] = fFileName[0] + "." + extension;
620  theHistoManager->SetFirstHistoId(1);
621 
622  G4bool fileOpen = theHistoManager->OpenFile(fFileName[0]);
623  if (!fileOpen) {
624  G4cout << "\n---> RMC01AnalysisManager::Book(): cannot open "
625  << fFileName[1]
626  << G4endl;
627  return;
628  }
629 
630  // Create directories
631  theHistoManager->SetHistoDirectoryName("histo");
632 
633  //Histograms for :
634  // 1)the forward simulation results
635  // 2)the Reverse MC simulation results normalised to a user spectrum
636  //------------------------------------------------------------------------
637 
638  G4int idHisto =
639  theHistoManager->CreateH1(G4String("Edep_vs_prim_ekin"),
640  G4String("edep vs e- primary energy"),60,emin,emax,
641  "none","none",G4String("log"));
642  fEdep_vs_prim_ekin = theHistoManager->GetH1(idHisto);
643 
644  idHisto = theHistoManager->CreateH1(G4String("elecron_current"),
645  G4String("electron"),60,emin,emax,
646  "none","none",G4String("log"));
647 
648  fElectron_current = theHistoManager->GetH1(idHisto);
649 
650  idHisto= theHistoManager->CreateH1(G4String("proton_current"),
651  G4String("proton"),60,emin,emax,
652  "none","none",G4String("log"));
653  fProton_current=theHistoManager->GetH1(idHisto);
654 
655  idHisto= theHistoManager->CreateH1(G4String("gamma_current"),
656  G4String("gamma"),60,emin,emax,
657  "none","none",G4String("log"));
658  fGamma_current=theHistoManager->GetH1(idHisto);
659 
660  //Response matrices for the adjoint simulation only
661  //-----------------------------------------------
662  if (fAdjoint_sim_mode){
663  //Response matrices for external isotropic e- source
664  //--------------------------------------------------
665 
666  idHisto =
667  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_electron_prim_energy"),
668  G4String("electron RM vs e- primary energy"),60,emin,emax,
669  "none","none",G4String("log"));
670  fEdep_rmatrix_vs_electron_prim_energy = theHistoManager->GetH1(idHisto);
671 
672  idHisto =
673  theHistoManager->
674  CreateH2(G4String("Electron_current_rmatrix_vs_electron_prim_energy"),
675  G4String("electron current RM vs e- primary energy"),
676  60,emin,emax,60,emin,emax,
677  "none","none","none","none",G4String("log"),G4String("log"));
678 
679  fElectron_current_rmatrix_vs_electron_prim_energy =
680  theHistoManager->GetH2(idHisto);
681 
682  idHisto =
683  theHistoManager->
684  CreateH2(G4String("Gamma_current_rmatrix_vs_electron_prim_energy"),
685  G4String("gamma current RM vs e- primary energy"),
686  60,emin,emax,60,emin,emax,
687  "none","none","none","none",G4String("log"),G4String("log"));
688 
689  fGamma_current_rmatrix_vs_electron_prim_energy =
690  theHistoManager->GetH2(idHisto);
691 
692  //Response matrices for external isotropic gamma source
693 
694  idHisto =
695  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_gamma_prim_energy"),
696  G4String("electron RM vs gamma primary energy"),60,emin,emax,
697  "none","none",G4String("log"));
698  fEdep_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH1(idHisto);
699 
700  idHisto =
701  theHistoManager->
702  CreateH2(G4String("Electron_current_rmatrix_vs_gamma_prim_energy"),
703  G4String("electron current RM vs gamma primary energy"),
704  60,emin,emax,60,emin,emax,
705  "none","none","none","none",G4String("log"),G4String("log"));
706 
707  fElectron_current_rmatrix_vs_gamma_prim_energy =
708  theHistoManager->GetH2(idHisto);
709 
710  idHisto =
711  theHistoManager->
712  CreateH2(G4String("Gamma_current_rmatrix_vs_gamma_prim_energy"),
713  G4String("gamma current RM vs gamma primary energy"),
714  60,emin,emax,60,emin,emax,
715  "none","none","none","none",G4String("log"),G4String("log"));
716 
717  fGamma_current_rmatrix_vs_gamma_prim_energy =
718  theHistoManager->GetH2(idHisto);
719 
720  //Response matrices for external isotropic proton source
721  idHisto =
722  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_proton_prim_energy"),
723  G4String("electron RM vs proton primary energy"),60,emin,emax,
724  "none","none",G4String("log"));
725  fEdep_rmatrix_vs_proton_prim_energy = theHistoManager->GetH1(idHisto);
726 
727  idHisto =
728  theHistoManager->
729  CreateH2(G4String("Electron_current_rmatrix_vs_proton_prim_energy"),
730  G4String("electron current RM vs proton primary energy"),
731  60,emin,emax,60,emin,emax,
732  "none","none","none","none",G4String("log"),G4String("log"));
733 
734  fElectron_current_rmatrix_vs_proton_prim_energy =
735  theHistoManager->GetH2(idHisto);
736 
737  idHisto =
738  theHistoManager->
739  CreateH2(G4String("Gamma_current_rmatrix_vs_proton_prim_energy"),
740  G4String("gamma current RM vs proton primary energy"),
741  60,emin,emax,60,emin,emax,
742  "none","none","none","none",G4String("log"),G4String("log"));
743 
744  fGamma_current_rmatrix_vs_proton_prim_energy =
745  theHistoManager->GetH2(idHisto);
746 
747  idHisto =
748  theHistoManager->
749  CreateH2(G4String("Proton_current_rmatrix_vs_proton_prim_energy"),
750  G4String("proton current RM vs proton primary energy"),
751  60,emin,emax,60,emin,emax,
752  "none","none","none","none",G4String("log"),G4String("log"));
753 
754  fProton_current_rmatrix_vs_proton_prim_energy =
755  theHistoManager->GetH2(idHisto);
756  }
757  fFactoryOn = true;
758  G4cout << "\n----> Histogram Tree is opened in " << fFileName[1] << G4endl;
759 }
760 
761 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
762 
763 void RMC01AnalysisManager::Save(G4double scaling_factor)
764 { if (fFactoryOn) {
765  G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
766  //scaling of results
767  //-----------------
768 
769  for (int ind=1; ind<=theHistoManager->GetNofH1s();ind++){
770  theHistoManager->SetH1Ascii(ind,true);
771  theHistoManager->ScaleH1(ind,scaling_factor);
772  }
773  for (int ind=1; ind<=theHistoManager->GetNofH2s();ind++)
774  theHistoManager->ScaleH2(ind,scaling_factor);
775 
776  theHistoManager->Write();
777  theHistoManager->CloseFile();
778  G4cout << "\n----> Histogram Tree is saved in " << fFileName[1] << G4endl;
779 
780  delete G4AnalysisManager::Instance();
781  fFactoryOn = false;
782  }
783 }
RMC01AnalysisManager::fEdep_rmatrix_vs_proton_prim_energy
G4AnaH1 * fEdep_rmatrix_vs_proton_prim_energy
Definition: RMC01AnalysisManager.hh:160
G4AdjointSimManager::GetInstance
static G4AdjointSimManager * GetInstance()
Definition: G4AdjointSimManager.cc:134
G4VAnalysisManager::SetHistoDirectoryName
G4bool SetHistoDirectoryName(const G4String &dirName)
Definition: G4VAnalysisManager.cc:250
G4RunManager::AbortRun
virtual void AbortRun(G4bool softAbort=false)
Definition: G4RunManager.cc:617
POWER
Definition: RMC01AnalysisManager.hh:71
RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim
void SetPrimaryPowerLawSpectrumForAdjointSim(const G4String &particle_name, G4double fluence, G4double alpha, G4double Emin, G4double Emax)
Definition: RMC01AnalysisManager.cc:567
G4Timer.hh
EXPO
Definition: RMC01AnalysisManager.hh:71
G4PrimaryParticle::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4PrimaryParticle.hh:219
RMC01AnalysisManager::fEdep_rmatrix_vs_electron_prim_energy
G4AnaH1 * fEdep_rmatrix_vs_electron_prim_energy
Definition: RMC01AnalysisManager.hh:152
G4VAnalysisManager::SetFirstHistoId
G4bool SetFirstHistoId(G4int firstId)
Definition: G4VAnalysisManager.cc:877
G4HCofThisEvent::GetHC
G4VHitsCollection * GetHC(G4int i)
Definition: G4HCofThisEvent.hh:83
G4SDManager::GetCollectionID
G4int GetCollectionID(G4String colName)
Definition: G4SDManager.cc:135
G4VAnalysisManager::CreateH1
G4int CreateH1(const G4String &name, const G4String &title, G4int nbins, G4double xmin, G4double xmax, const G4String &unitName="none", const G4String &fcnName="none", const G4String &binSchemeName="linear")
Definition: G4VAnalysisManager.cc:292
cm2
static constexpr double cm2
Definition: G4SIunits.hh:120
G4RunManager.hh
G4VAnalysisManager::ScaleH1
G4bool ScaleH1(G4int id, G4double factor)
Definition: G4VAnalysisManager.cc:520
G4VAnalysisManager::ScaleH2
G4bool ScaleH2(G4int id, G4double factor)
Definition: G4VAnalysisManager.cc:526
G4float
float G4float
Definition: G4Types.hh:77
RMC01AnalysisManager::fGamma_current_rmatrix_vs_proton_prim_energy
G4AnaH2 * fGamma_current_rmatrix_vs_proton_prim_energy
Definition: RMC01AnalysisManager.hh:163
RMC01AnalysisManager::fElectron_current_rmatrix_vs_electron_prim_energy
G4AnaH2 * fElectron_current_rmatrix_vs_electron_prim_energy
Definition: RMC01AnalysisManager.hh:153
G4THitsCollection::entries
G4int entries() const
Definition: G4THitsCollection.hh:114
RMC01AnalysisManager::fElectron_current_rmatrix_vs_gamma_prim_energy
G4AnaH2 * fElectron_current_rmatrix_vs_gamma_prim_energy
Definition: RMC01AnalysisManager.hh:157
RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim
void SetPrimaryExpSpectrumForAdjointSim(const G4String &particle_name, G4double fluence, G4double E0, G4double Emin, G4double Emax)
Definition: RMC01AnalysisManager.cc:542
RMC01AnalysisManager::BeginOfEvent
void BeginOfEvent(const G4Event *)
Definition: RMC01AnalysisManager.cc:176
RMC01AnalysisManager::fPrimSpectrumType
PRIM_SPECTRUM_TYPE fPrimSpectrumType
Definition: RMC01AnalysisManager.hh:173
G4int
int G4int
Definition: G4Types.hh:78
RMC01AnalysisManager::fProton_current_rmatrix_vs_proton_prim_energy
G4AnaH2 * fProton_current_rmatrix_vs_proton_prim_energy
Definition: RMC01AnalysisManager.hh:162
G4VAnalysisManager::OpenFile
G4bool OpenFile(const G4String &fileName="")
Definition: G4VAnalysisManager.cc:185
G4PrimaryParticle::GetG4code
G4ParticleDefinition * GetG4code() const
Definition: G4PrimaryParticle.hh:207
P
static double P[]
Definition: nf_gammaFunctions.cc:89
G4Event::GetEventID
G4int GetEventID() const
Definition: G4Event.hh:151
RMC01SD.hh
Definition of the RMC01SD class.
G4PrimaryVertex::GetPrimary
G4PrimaryParticle * GetPrimary(G4int i=0) const
Definition: G4PrimaryVertex.cc:140
G4VAnalysisManager::SetH1Ascii
void SetH1Ascii(G4int id, G4bool ascii)
Definition: G4VAnalysisManager.cc:1108
G4cout
G4GLOB_DLL std::ostream G4cout
G4SDManager.hh
G4Run::GetNumberOfEvent
G4int GetNumberOfEvent() const
Definition: G4Run.hh:79
G4bool
bool G4bool
Definition: G4Types.hh:79
G4VAnalysisManager::GetFileType
G4String GetFileType() const
G4Proton.hh
G4Run
Definition: G4Run.hh:46
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
RMC01AnalysisManager::EndOfRun
void EndOfRun(const G4Run *)
Definition: RMC01AnalysisManager.cc:153
G4Event::GetPrimaryVertex
G4PrimaryVertex * GetPrimaryVertex(G4int i=0) const
Definition: G4Event.hh:167
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
RMC01AnalysisManager.hh
Definition of the RMC01AnalysisManager class.
G4AdjointSimManager::GetWeightAtEndOfLastAdjointTrack
G4double GetWeightAtEndOfLastAdjointTrack(size_t i=0)
Definition: G4AdjointSimManager.cc:392
RMC01AnalysisManager::PrimDiffAndDirFluxForAdjointSim
G4double PrimDiffAndDirFluxForAdjointSim(G4double prim_energy)
Definition: RMC01AnalysisManager.cc:456
G4PhysicalConstants.hh
emax
static const G4double emax
Definition: G4ChatterjeeCrossSection.hh:39
RMC01AnalysisManager::fMsg
RMC01AnalysisManagerMessenger * fMsg
Definition: RMC01AnalysisManager.hh:122
RMC01AnalysisManager::Save
void Save(G4double scaling_factor)
Definition: RMC01AnalysisManager.cc:763
G4RunManager::GetRunManager
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
RMC01AnalysisManager::ComputeMeanEdepAndError
void ComputeMeanEdepAndError(const G4Event *anEvent, G4double &mean, G4double &error)
Definition: RMC01AnalysisManager.cc:507
G4Electron.hh
RMC01AnalysisManagerMessenger.hh
Definition of the RMC01AnalysisManagerMessenger class.
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
G4ToolsAnalysisManager::GetH2
tools::histo::h2d * GetH2(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
G4Timer::GetRealElapsed
G4double GetRealElapsed() const
Definition: G4Timer.cc:107
RMC01AnalysisManager::fGamma_current_rmatrix_vs_electron_prim_energy
G4AnaH2 * fGamma_current_rmatrix_vs_electron_prim_energy
Definition: RMC01AnalysisManager.hh:154
G4Timer::Stop
void Stop()
G4AdjointSimManager.hh
RMC01AnalysisManager::fElectron_current_rmatrix_vs_proton_prim_energy
G4AnaH2 * fElectron_current_rmatrix_vs_proton_prim_energy
Definition: RMC01AnalysisManager.hh:161
RMC01AnalysisManager::GetInstance
static RMC01AnalysisManager * GetInstance()
Definition: RMC01AnalysisManager.cc:108
G4SDManager::GetSDMpointer
static G4SDManager * GetSDMpointer()
Definition: G4SDManager.cc:40
RMC01AnalysisManager::EndOfEvent
void EndOfEvent(const G4Event *)
Definition: RMC01AnalysisManager.cc:182
RMC01AnalysisManager::BeginOfRun
void BeginOfRun(const G4Run *)
Definition: RMC01AnalysisManager.cc:116
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
Emin
static const G4double Emin
Definition: G4PhotoNuclearCrossSection.cc:50
G4AdjointSimManager::GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack
G4int GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack(size_t i=0)
Definition: G4AdjointSimManager.cc:407
G4Run::GetNumberOfEventToBeProcessed
G4int GetNumberOfEventToBeProcessed() const
Definition: G4Run.hh:83
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
Emax
static const G4double Emax
Definition: G4PhotoNuclearCrossSection.cc:51
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4Root::G4AnaH1
tools::histo::h1d G4AnaH1
Definition: g4root_defs.hh:46
pi
static constexpr double pi
Definition: G4SIunits.hh:75
error
static PROLOG_HANDLER error
Definition: xmlrole.cc:112
G4Root::G4AnaH2
tools::histo::h2d G4AnaH2
Definition: g4root_defs.hh:52
G4Timer::Start
void Start()
G4Event::GetHCofThisEvent
G4HCofThisEvent * GetHCofThisEvent() const
Definition: G4Event.hh:185
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
alpha
static const G4double alpha
Definition: G4LivermoreBremsstrahlungModel.cc:79
RMC01AnalysisManager::fGamma_current_rmatrix_vs_gamma_prim_energy
G4AnaH2 * fGamma_current_rmatrix_vs_gamma_prim_energy
Definition: RMC01AnalysisManager.hh:158
G4ToolsAnalysisManager::GetH1
tools::histo::h1d * GetH1(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
keV
static constexpr double keV
Definition: G4SIunits.hh:216
G4AdjointSimManager::GetEkinAtEndOfLastAdjointTrack
G4double GetEkinAtEndOfLastAdjointTrack(size_t i=0)
Definition: G4AdjointSimManager.cc:382
RMC01AnalysisManager::EndOfEventForAdjointSimulation
void EndOfEventForAdjointSimulation(const G4Event *anEvent)
Definition: RMC01AnalysisManager.cc:290
G4AdjointSimManager::GetNbOfAdointTracksReachingTheExternalSurface
size_t GetNbOfAdointTracksReachingTheExternalSurface()
Definition: G4AdjointSimManager.cc:418
RMC01AnalysisManager::EndOfEventForForwardSimulation
void EndOfEventForForwardSimulation(const G4Event *anEvent)
Definition: RMC01AnalysisManager.cc:221
G4Gamma.hh
RMC01AnalysisManager::fInstance
static RMC01AnalysisManager * fInstance
Definition: RMC01AnalysisManager.hh:106
RMC01AnalysisManager::fEdep_rmatrix_vs_gamma_prim_energy
G4AnaH1 * fEdep_rmatrix_vs_gamma_prim_energy
Definition: RMC01AnalysisManager.hh:156