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