Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 //
26 /// \file biasing/ReverseMC01/src/RMC01AnalysisManager.cc
27 /// \brief Implementation of the RMC01AnalysisManager class
28 //
29 // $Id: RMC01AnalysisManager.cc 76461 2013-11-11 10:15:51Z gcosmo $
30 //
31 //////////////////////////////////////////////////////////////
32 // Class Name: RMC01AnalysisManager
33 // Author: L. Desorgher
34 // Organisation: SpaceIT GmbH
35 // Contract: ESA contract 21435/08/NL/AT
36 // Customer: ESA/ESTEC
37 //////////////////////////////////////////////////////////////
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  fRelative_error=1.;
124  fMean_edep=0.;
125  fError_mean_edep=0.;
126  fAdjoint_sim_mode =G4AdjointSimManager::GetInstance()->GetAdjointSimMode();
127 
128  if (fAdjoint_sim_mode){
129  fNb_evt_per_adj_evt=aRun->GetNumberOfEventToBeProcessed()/
130  G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
131  fConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt",
132  std::ios::out);
133  fConvergenceFileOutput<<
134  "Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
135  }
136  else {
137  fConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt",
138  std::ios::out);
139  fConvergenceFileOutput<<
140  "Edep per event [MeV]\terror[MeV]\tcomputing_time[s]"
141  <<std::endl;
142  }
143  fConvergenceFileOutput.setf(std::ios::scientific);
144  fConvergenceFileOutput.precision(6);
145 
146  fTimer->Start();
147  fElapsed_time=0.;
148 
149  book();
150 }
151 
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
153 
154 void RMC01AnalysisManager::EndOfRun(const G4Run* aRun)
155 { fTimer->Stop();
156  G4int nb_evt=aRun->GetNumberOfEvent();
157  G4double factor =1./ nb_evt;
158  if (!fAdjoint_sim_mode){
159  G4cout<<"Results of forward simulation!"<<std::endl;
160  G4cout<<"edep per event [MeV] = "<<fMean_edep<<std::endl;
161  G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
162  }
163 
164  else {
165  G4cout<<"Results of reverse/adjoint simulation!"<<std::endl;
166  G4cout<<"normalised edep [MeV] = "<<fMean_edep<<std::endl;
167  G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
168  factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun()
169  *fNb_evt_per_adj_evt/aRun->GetNumberOfEvent();
170  }
171  save(factor);
172  fConvergenceFileOutput.close();
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176 
177 void RMC01AnalysisManager::BeginOfEvent(const G4Event* )
178 { ;
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 
183 void RMC01AnalysisManager::EndOfEvent(const G4Event* anEvent)
184 {
185  if (fAdjoint_sim_mode) EndOfEventForAdjointSimulation(anEvent);
186  else EndOfEventForForwardSimulation(anEvent);
187 
188  //Test convergence. The error is already computed
189  //--------------------------------------
190  G4int nb_event=anEvent->GetEventID()+1;
191  //G4double factor=1.;
192  if (fAdjoint_sim_mode) {
193  G4double n_adj_evt= nb_event/fNb_evt_per_adj_evt;
194  // nb_event/fNb_evt_per_adj_evt;
195  if (n_adj_evt*fNb_evt_per_adj_evt == nb_event) {
196  nb_event =static_cast<G4int>(n_adj_evt);
197  }
198  else nb_event=0;
199 
200  }
201 
202  if (nb_event>100 && fStop_run_if_precision_reached &&
203  fPrecision_to_reach >fRelative_error) {
204  G4cout<<fPrecision_to_reach*100.<<"% Precision reached!"<<std::endl;
205  fTimer->Stop();
206  fElapsed_time+=fTimer->GetRealElapsed();
207  fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep
208  <<'\t'<<fElapsed_time<<std::endl;
209  G4RunManager::GetRunManager()->AbortRun(true);
210  }
211 
212  if (nb_event>0 && nb_event % fNb_evt_modulo_for_convergence_test == 0) {
213  fTimer->Stop();
214  fElapsed_time+=fTimer->GetRealElapsed();
215  fTimer->Start();
216  fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep<<'\t'
217  <<fElapsed_time<<std::endl;
218  }
219 }
220 
221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
222 
223 void RMC01AnalysisManager::EndOfEventForForwardSimulation(
224  const G4Event* anEvent)
225 {
226 
227  G4SDManager* SDman = G4SDManager::GetSDMpointer();
228  G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
229  RMC01DoubleWithWeightHitsCollection* edepCollection =
230  (RMC01DoubleWithWeightHitsCollection*)
231  (HCE->GetHC(SDman->GetCollectionID("edep")));
232 
233  RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
234  (RMC01DoubleWithWeightHitsCollection*)
235  (HCE->GetHC(SDman->GetCollectionID("current_electron")));
236 
237  RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
238  (RMC01DoubleWithWeightHitsCollection*)
239  (HCE->GetHC(SDman->GetCollectionID("current_proton")));
240 
241  RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
242  (RMC01DoubleWithWeightHitsCollection*)
243  (HCE->GetHC(SDman->GetCollectionID("current_gamma")));
244 
245  //Total energy deposited in Event
246  //-------------------------------
247  G4double totEdep=0;
248  G4int i;
249  for (i=0;i<edepCollection->entries();i++)
250  totEdep+=(*edepCollection)[i]->GetValue()
251  *(*edepCollection)[i]->GetWeight();
252 
253  if (totEdep>0.){
254  fAccumulated_edep +=totEdep ;
255  fAccumulated_edep2 +=totEdep*totEdep;
256  G4PrimaryParticle* thePrimary=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  //Output from adjoint tracking phase
314  //----------------------------------------------------------------------------
315 
316  G4AdjointSimManager* theAdjointSimManager =
317  G4AdjointSimManager::GetInstance();
318  G4int pdg_nb =theAdjointSimManager
319  ->GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack();
320  G4double prim_ekin=theAdjointSimManager->GetEkinAtEndOfLastAdjointTrack();
321  G4double adj_weight=theAdjointSimManager->GetWeightAtEndOfLastAdjointTrack();
322 
323 
324  //Factor of normalisation to user defined prim spectrum (power law or exp)
325  //---------------------------------------------------------------------------
326 
327  G4double normalised_weight = 0.;
328  if (pdg_nb== fPrimPDG_ID && prim_ekin>= fEmin_prim_spectrum
329  && prim_ekin<= fEmax_prim_spectrum)
330  normalised_weight =
331  adj_weight*PrimDiffAndDirFluxForAdjointSim(prim_ekin);
332 
333  //Answer matrices
334  //-------------
335  G4AnaH1* edep_rmatrix =0;
336  G4AnaH2* electron_current_rmatrix =0;
337  G4AnaH2* gamma_current_rmatrix =0;
338  G4AnaH2* proton_current_rmatrix =0;
339 
340  if (pdg_nb == G4Electron::Electron()->GetPDGEncoding()){ //e- answer matrices
341  edep_rmatrix = fEdep_rmatrix_vs_electron_prim_energy;
342 
343  electron_current_rmatrix =
344  fElectron_current_rmatrix_vs_electron_prim_energy;
345 
346  gamma_current_rmatrix = fGamma_current_rmatrix_vs_electron_prim_energy;
347  }
348  else if (pdg_nb == G4Gamma::Gamma()->GetPDGEncoding()){
349  //gammma answer matrices
350  edep_rmatrix = fEdep_rmatrix_vs_gamma_prim_energy;
351  electron_current_rmatrix = fElectron_current_rmatrix_vs_gamma_prim_energy;
352  gamma_current_rmatrix = fGamma_current_rmatrix_vs_gamma_prim_energy;
353  }
354  else if (pdg_nb == G4Proton::Proton()->GetPDGEncoding()){
355  //proton answer matrices
356  edep_rmatrix = fEdep_rmatrix_vs_proton_prim_energy;
357  electron_current_rmatrix = fElectron_current_rmatrix_vs_proton_prim_energy;
358  gamma_current_rmatrix = fGamma_current_rmatrix_vs_proton_prim_energy;
359  proton_current_rmatrix = fProton_current_rmatrix_vs_proton_prim_energy;
360  }
361 
362  //Registering of total energy deposited in Event
363  //-------------------------------
364  G4double totEdep=0;
365  G4int i;
366  for (i=0;i<edepCollection->entries();i++)
367  totEdep+=(*edepCollection)[i]->GetValue()*
368  (*edepCollection)[i]->GetWeight();
369 
370  G4bool new_mean_computed=false;
371  if (totEdep>0.){
372  if (normalised_weight>0.){
373  G4double edep=totEdep* normalised_weight;
374 
375  //Check if the edep is not wrongly too high
376  //-----------------------------------------
377  G4double new_mean , new_error;
378 
379  fAccumulated_edep +=edep;
380  fAccumulated_edep2 +=edep*edep;
381 
382  ComputeMeanEdepAndError(anEvent,new_mean,new_error);
383  G4double new_relative_error = 1.;
384  if ( new_error >0) new_relative_error = new_error/ new_mean;
385  if (fRelative_error <0.10 && new_relative_error>1.5*fRelative_error) {
386  G4cout<<"Potential wrong adjoint weight!"<<std::endl;
387  G4cout<<"The results of this event will not be registered!"
388  <<std::endl;
389  G4cout<<"previous mean edep [MeV] "<< fMean_edep<<std::endl;
390  G4cout<<"previous relative error "<< fRelative_error<<std::endl;
391  G4cout<<"new rejected mean edep [MeV] "<< new_mean<<std::endl;
392  G4cout<<"new rejected relative error "<< new_relative_error
393  <<std::endl;
394  fAccumulated_edep -=edep;
395  fAccumulated_edep2 -=edep*edep;
396  return;
397  }
398  else { //accepted
399  fMean_edep = new_mean;
400  fError_mean_edep = new_error;
401  fRelative_error =new_relative_error;
402  new_mean_computed=true;
403  }
404  fEdep_vs_prim_ekin->fill(prim_ekin,edep);
405  }
406 
407  // Registering answer matrix
408  //---------------------------
409 
410  edep_rmatrix->fill(prim_ekin,totEdep*adj_weight/cm2);
411  }
412  if (!new_mean_computed){
413  ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
414  if (fError_mean_edep>0) fRelative_error= fError_mean_edep/fMean_edep;
415  }
416 
417 
418  //Registering of current of particles on the sensitive volume
419  //------------------------------------------------------------
420 
421  for (i=0;i<electronCurrentCollection->entries();i++) {
422  G4double ekin =(*electronCurrentCollection)[i]->GetValue();
423  G4double weight=(*electronCurrentCollection)[i]->GetWeight();
424  fElectron_current->fill(ekin,weight*normalised_weight);
425  electron_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
426  }
427 
428  for (i=0;i<protonCurrentCollection->entries();i++) {
429  G4double ekin =(*protonCurrentCollection)[i]->GetValue();
430  G4double weight=(*protonCurrentCollection)[i]->GetWeight();
431  fProton_current->fill(ekin,weight*normalised_weight);
432  proton_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
433  }
434 
435  for (i=0;i<gammaCurrentCollection->entries();i++) {
436  G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
437  G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
438  fGamma_current->fill(ekin,weight*normalised_weight);
439  gamma_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
440  }
441 }
442 
443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
444 
445 G4double RMC01AnalysisManager::PrimDiffAndDirFluxForAdjointSim(
446  G4double prim_energy)
447 {
448  G4double flux=fAmplitude_prim_spectrum;
449  if ( fPrimSpectrumType ==EXPO) flux*=std::exp(-prim_energy/fAlpha_or_E0);
450  else flux*=std::pow(prim_energy, -fAlpha_or_E0);
451  return flux;
452 }
453 
454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
455 /*
456 void RMC01AnalysisManager::WriteHisto(G4AnaH1* anHisto,
457  G4double scaling_factor, G4String fileName, G4String header_lines)
458 { std::fstream FileOutput(fileName, std::ios::out);
459  FileOutput<<header_lines;
460  FileOutput.setf(std::ios::scientific);
461  FileOutput.precision(6);
462 
463  for (G4int i =0;i<G4int(anHisto->axis().bins());i++) {
464  FileOutput<<anHisto->axis().bin_lower_edge(i)
465  <<'\t'<<anHisto->axis().bin_upper_edge(i)
466  <<'\t'<<anHisto->bin_height(i)*scaling_factor
467  <<'\t'<<anHisto->bin_error(i)*scaling_factor<<std::endl;
468  }
469 }
470 
471 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
472 
473 void RMC01AnalysisManager::WriteHisto(G4AnaH2* anHisto,
474  G4double scaling_factor, G4String fileName, G4String header_lines)
475 { std::fstream FileOutput(fileName, std::ios::out);
476  FileOutput<<header_lines;
477 
478  FileOutput.setf(std::ios::scientific);
479  FileOutput.precision(6);
480 
481  for (G4int i =0;i<G4int(anHisto->axis_x().bins());i++) {
482  for (G4int j =0;j<G4int(anHisto->axis_y().bins());j++) {
483  FileOutput<<anHisto->axis_x().bin_lower_edge(i)
484  <<'\t'<<anHisto->axis_x().bin_upper_edge(i)
485  <<'\t'<<anHisto->axis_y().bin_lower_edge(i)
486  <<'\t'<<anHisto->axis_y().bin_upper_edge(i)
487  <<'\t'<<anHisto->bin_height(i,j)*scaling_factor
488  <<'\t'<<anHisto->bin_error(i,j)*scaling_factor
489  <<std::endl;
490  }
491  }
492 }
493 */
494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
495 
496 void RMC01AnalysisManager::ComputeMeanEdepAndError(
497  const G4Event* anEvent,G4double& mean,G4double& error)
498 {
499  G4int nb_event=anEvent->GetEventID()+1;
500  G4double factor=1.;
501  if (fAdjoint_sim_mode) {
502  nb_event /=fNb_evt_per_adj_evt;
503  factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
504  }
505 
506  //error computation
507  if (nb_event>1) {
508  mean = fAccumulated_edep/nb_event;
509  G4double mean_x2 =fAccumulated_edep2/nb_event;
510  error = factor*std::sqrt(mean_x2-mean*mean)/std::sqrt(G4double(nb_event));
511  mean *=factor;
512  }
513  else {
514  mean=0;
515  error=0;
516  }
517 }
518 
519 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
520 
521 void RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim(
522  const G4String& particle_name, G4double omni_fluence,
523  G4double E0, G4double Emin,
524  G4double Emax)
525 { fPrimSpectrumType = EXPO;
526  if (particle_name == "e-" ) fPrimPDG_ID =
527  G4Electron::Electron()->GetPDGEncoding();
528  else if (particle_name == "gamma") fPrimPDG_ID =
529  G4Gamma::Gamma()->GetPDGEncoding();
530  else if (particle_name == "proton") fPrimPDG_ID =
531  G4Proton::Proton()->GetPDGEncoding();
532  else {
533  G4cout<<"The particle that you did select is not in the candidate "<<
534  "list for primary [e-, gamma, proton]!"<<G4endl;
535  return;
536  }
537  fAlpha_or_E0 = E0 ;
538  fAmplitude_prim_spectrum = omni_fluence/E0/
539  (std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
540  fEmin_prim_spectrum = Emin ;
541  fEmax_prim_spectrum = Emax;
542 }
543 
544 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
545 
546 void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim(
547  const G4String& particle_name, G4double omni_fluence,
548  G4double alpha, G4double Emin,G4double Emax)
549 { fPrimSpectrumType =POWER;
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 
562 
563  if (alpha ==1.) {
564  fAmplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
565  }
566  else {
567  G4double p=1.-alpha;
568  fAmplitude_prim_spectrum = omni_fluence/p/(std::pow(Emax,p)
569  -std::pow(Emin,p))/4./pi;
570  }
571 
572  fAlpha_or_E0 = alpha ;
573  fEmin_prim_spectrum = Emin ;
574  fEmax_prim_spectrum = Emax;
575 
576 }
577 
578 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
579 
580 void RMC01AnalysisManager::book()
581 {
582  //----------------------
583  //Creation of histograms
584  //----------------------
585 
586  //Energy binning of the histograms : 60 log bins over [1keV-1GeV]
587 
588  G4double emin=1.*keV;
589  G4double emax=1.*GeV;
590 
591  //file_name
592  fFileName[0]="forward_sim";
593  if (fAdjoint_sim_mode) fFileName[0]="adjoint_sim";
594 
595  //Histo manager
596  G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
597  G4String extension = theHistoManager->GetFileType();
598  fFileName[1] = fFileName[0] + "." + extension;
599  theHistoManager->SetFirstHistoId(1);
600 
601  G4bool fileOpen = theHistoManager->OpenFile(fFileName[0]);
602  if (!fileOpen) {
603  G4cout << "\n---> RMC01AnalysisManager::book(): cannot open " << fFileName[1]
604  << G4endl;
605  return;
606  }
607 
608  // Create directories
609  theHistoManager->SetHistoDirectoryName("histo");
610 
611  //Histograms for :
612  // 1)the forward simulation results
613  // 2)the Reverse MC simulation results normalised to a user spectrum
614  //------------------------------------------------------------------------
615 
616 
617  G4int idHisto =
618  theHistoManager->CreateH1(G4String("Edep_vs_prim_ekin"),
619  G4String("edep vs e- primary energy"),60,emin,emax,
620  "none","none",G4String("log"));
621  fEdep_vs_prim_ekin = theHistoManager->GetH1(idHisto);
622 
623  idHisto = theHistoManager->CreateH1(G4String("elecron_current"),
624  G4String("electron"),60,emin,emax,
625  "none","none",G4String("log"));
626 
627  fElectron_current = theHistoManager->GetH1(idHisto);
628 
629  idHisto= theHistoManager->CreateH1(G4String("proton_current"),
630  G4String("proton"),60,emin,emax,
631  "none","none",G4String("log"));
632  fProton_current=theHistoManager->GetH1(idHisto);
633 
634 
635  idHisto= theHistoManager->CreateH1(G4String("gamma_current"),
636  G4String("gamma"),60,emin,emax,
637  "none","none",G4String("log"));
638  fGamma_current=theHistoManager->GetH1(idHisto);
639 
640 
641  //Response matrices for the adjoint simulation only
642  //-----------------------------------------------
643  if (fAdjoint_sim_mode){
644  //Response matrices for external isotropic e- source
645  //--------------------------------------------------
646 
647  idHisto =
648  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_electron_prim_energy"),
649  G4String("electron RM vs e- primary energy"),60,emin,emax,
650  "none","none",G4String("log"));
651  fEdep_rmatrix_vs_electron_prim_energy = theHistoManager->GetH1(idHisto);
652 
653  idHisto =
654  theHistoManager->
655  CreateH2(G4String("Electron_current_rmatrix_vs_electron_prim_energy"),
656  G4String("electron current RM vs e- primary energy"),
657  60,emin,emax,60,emin,emax,
658  "none","none","none","none",G4String("log"),G4String("log"));
659 
660  fElectron_current_rmatrix_vs_electron_prim_energy =
661  theHistoManager->GetH2(idHisto);
662 
663  idHisto =
664  theHistoManager->
665  CreateH2(G4String("Gamma_current_rmatrix_vs_electron_prim_energy"),
666  G4String("gamma current RM vs e- primary energy"),
667  60,emin,emax,60,emin,emax,
668  "none","none","none","none",G4String("log"),G4String("log"));
669 
670 
671  fGamma_current_rmatrix_vs_electron_prim_energy =
672  theHistoManager->GetH2(idHisto);
673 
674 
675  //Response matrices for external isotropic gamma source
676 
677  idHisto =
678  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_gamma_prim_energy"),
679  G4String("electron RM vs gamma primary energy"),60,emin,emax,
680  "none","none",G4String("log"));
681  fEdep_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH1(idHisto);
682 
683  idHisto =
684  theHistoManager->
685  CreateH2(G4String("Electron_current_rmatrix_vs_gamma_prim_energy"),
686  G4String("electron current RM vs gamma primary energy"),
687  60,emin,emax,60,emin,emax,
688  "none","none","none","none",G4String("log"),G4String("log"));
689 
690  fElectron_current_rmatrix_vs_gamma_prim_energy =
691  theHistoManager->GetH2(idHisto);
692 
693  idHisto =
694  theHistoManager->
695  CreateH2(G4String("Gamma_current_rmatrix_vs_gamma_prim_energy"),
696  G4String("gamma current RM vs gamma primary energy"),
697  60,emin,emax,60,emin,emax,
698  "none","none","none","none",G4String("log"),G4String("log"));
699 
700  fGamma_current_rmatrix_vs_gamma_prim_energy =
701  theHistoManager->GetH2(idHisto);
702 
703 
704 
705  //Response matrices for external isotropic proton source
706  idHisto =
707  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_proton_prim_energy"),
708  G4String("electron RM vs proton primary energy"),60,emin,emax,
709  "none","none",G4String("log"));
710  fEdep_rmatrix_vs_proton_prim_energy = theHistoManager->GetH1(idHisto);
711 
712  idHisto =
713  theHistoManager->
714  CreateH2(G4String("Electron_current_rmatrix_vs_proton_prim_energy"),
715  G4String("electron current RM vs proton primary energy"),
716  60,emin,emax,60,emin,emax,
717  "none","none","none","none",G4String("log"),G4String("log"));
718 
719  fElectron_current_rmatrix_vs_proton_prim_energy =
720  theHistoManager->GetH2(idHisto);
721 
722  idHisto =
723  theHistoManager->
724  CreateH2(G4String("Gamma_current_rmatrix_vs_proton_prim_energy"),
725  G4String("gamma current RM vs proton primary energy"),
726  60,emin,emax,60,emin,emax,
727  "none","none","none","none",G4String("log"),G4String("log"));
728 
729  fGamma_current_rmatrix_vs_proton_prim_energy =
730  theHistoManager->GetH2(idHisto);
731 
732  idHisto =
733  theHistoManager->
734  CreateH2(G4String("Proton_current_rmatrix_vs_proton_prim_energy"),
735  G4String("proton current RM vs proton primary energy"),
736  60,emin,emax,60,emin,emax,
737  "none","none","none","none",G4String("log"),G4String("log"));
738 
739  fProton_current_rmatrix_vs_proton_prim_energy =
740  theHistoManager->GetH2(idHisto);
741  }
742  fFactoryOn = true;
743  G4cout << "\n----> Histogram Tree is opened in " << fFileName[1] << G4endl;
744 }
745 
746 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
747 
748 void RMC01AnalysisManager::save(G4double scaling_factor)
749 { if (fFactoryOn) {
750  G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
751  //scaling of results
752  //-----------------
753 
754  for (int ind=1; ind<=theHistoManager->GetNofH1s();ind++){
755  theHistoManager->SetH1Ascii(ind,true);
756  theHistoManager->ScaleH1(ind,scaling_factor);
757  }
758  for (int ind=1; ind<=theHistoManager->GetNofH2s();ind++)
759  theHistoManager->ScaleH2(ind,scaling_factor);
760 
761  theHistoManager->Write();
762  theHistoManager->CloseFile();
763  G4cout << "\n----> Histogram Tree is saved in " << fFileName[1] << G4endl;
764 
765  delete G4AnalysisManager::Instance();
766  fFactoryOn = false;
767  }
768 }
G4AdjointSimManager::GetInstance
static G4AdjointSimManager * GetInstance()
Definition: G4AdjointSimManager.cc:129
G4VAnalysisManager::SetHistoDirectoryName
G4bool SetHistoDirectoryName(const G4String &dirName)
Definition: G4VAnalysisManager.cc:185
G4RunManager::AbortRun
virtual void AbortRun(G4bool softAbort=false)
Definition: G4RunManager.cc:551
POWER
Definition: RMC01AnalysisManager.hh:71
RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim
void SetPrimaryPowerLawSpectrumForAdjointSim(const G4String &particle_name, G4double fluence, G4double alpha, G4double Emin, G4double Emax)
Definition: RMC01AnalysisManager.cc:546
G4Timer.hh
EXPO
Definition: RMC01AnalysisManager.hh:71
G4PrimaryParticle::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4PrimaryParticle.hh:218
G4VAnalysisManager::SetFirstHistoId
G4bool SetFirstHistoId(G4int firstId)
Definition: G4VAnalysisManager.cc:401
G4SDManager::GetCollectionID
G4int GetCollectionID(G4String colName)
Definition: G4SDManager.cc:131
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:215
G4RunManager.hh
G4Root::G4AnaH2
tools::histo::h2d G4AnaH2
Definition: g4root_defs.hh:41
G4VAnalysisManager::ScaleH1
G4bool ScaleH1(G4int id, G4double factor)
Definition: G4VAnalysisManager.cc:332
G4VAnalysisManager::ScaleH2
G4bool ScaleH2(G4int id, G4double factor)
Definition: G4VAnalysisManager.cc:338
p
const char * p
Definition: xmltok.h:285
G4AdjointSimManager::GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack
G4int GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack()
Definition: G4AdjointSimManager.cc:393
RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim
void SetPrimaryExpSpectrumForAdjointSim(const G4String &particle_name, G4double fluence, G4double E0, G4double Emin, G4double Emax)
Definition: RMC01AnalysisManager.cc:521
RMC01AnalysisManager::BeginOfEvent
void BeginOfEvent(const G4Event *)
Definition: RMC01AnalysisManager.cc:177
G4int
int G4int
Definition: G4Types.hh:78
G4PrimaryParticle::GetG4code
G4ParticleDefinition * GetG4code() const
Definition: G4PrimaryParticle.hh:206
G4VAnalysisManager::GetNofH1s
G4int GetNofH1s() const
Definition: G4VAnalysisManager.cc:487
G4Event::GetEventID
G4int GetEventID() const
Definition: G4Event.hh:140
RMC01SD.hh
Definition of the RMC01SD class.
RMC01AnalysisManager::save
void save(G4double scaling_factor)
Definition: RMC01AnalysisManager.cc:748
G4AdjointSimManager::GetEkinAtEndOfLastAdjointTrack
G4double GetEkinAtEndOfLastAdjointTrack()
Definition: G4AdjointSimManager.cc:368
G4PrimaryVertex::GetPrimary
G4PrimaryParticle * GetPrimary(G4int i=0) const
Definition: G4PrimaryVertex.cc:133
G4VAnalysisManager::SetH1Ascii
void SetH1Ascii(G4int id, G4bool ascii)
Definition: G4VAnalysisManager.cc:523
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:154
G4Event::GetPrimaryVertex
G4PrimaryVertex * GetPrimaryVertex(G4int i=0) const
Definition: G4Event.hh:156
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
RMC01AnalysisManager.hh
Definition of the RMC01AnalysisManager class.
G4PhysicalConstants.hh
G4RunManager::GetRunManager
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:74
G4Root::G4AnaH1
tools::histo::h1d G4AnaH1
Definition: g4root_defs.hh:40
G4Electron.hh
RMC01AnalysisManagerMessenger.hh
Definition of the RMC01AnalysisManagerMessenger class.
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4Timer::GetRealElapsed
G4double GetRealElapsed() const
Definition: G4Timer.cc:107
G4Timer::Stop
void Stop()
G4AdjointSimManager.hh
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:183
G4AdjointSimManager::GetWeightAtEndOfLastAdjointTrack
G4double GetWeightAtEndOfLastAdjointTrack()
Definition: G4AdjointSimManager.cc:378
RMC01AnalysisManager::BeginOfRun
void BeginOfRun(const G4Run *)
Definition: RMC01AnalysisManager.cc:117
G4Run::GetNumberOfEventToBeProcessed
G4int GetNumberOfEventToBeProcessed() const
Definition: G4Run.hh:83
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Timer::Start
void Start()
G4Event::GetHCofThisEvent
G4HCofThisEvent * GetHCofThisEvent() const
Definition: G4Event.hh:174
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
CLHEP::Hep3Vector::mag
double mag() const
G4RootAnalysisManager::GetH1
tools::histo::h1d * GetH1(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
G4RootAnalysisManager::GetH2
tools::histo::h2d * GetH2(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
G4VAnalysisManager::GetNofH2s
G4int GetNofH2s() const
Definition: G4VAnalysisManager.cc:493
G4Gamma.hh
Generated on Wed Apr 30 2014 15:55:28 for Geant4.10 by   doxygen 1.8.7