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