46 #include "G4THitsCollection.hh"
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),
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)
135 "Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
141 "Edep per event [MeV]\terror[MeV]\tcomputing_time[s]"
160 G4cout<<
"Results of forward simulation!"<<std::endl;
166 G4cout<<
"Results of reverse/adjoint simulation!"<<std::endl;
197 nb_event =
static_cast<G4int>(n_adj_evt);
250 for (i=0;i<edepCollection->
entries();i++)
251 totEdep+=(*edepCollection)[i]->GetValue()
252 *(*edepCollection)[i]->GetWeight();
261 G4double prim_ekin =std::sqrt(E0*E0+P*P)-E0;
270 for (i=0;i<electronCurrentCollection->
entries();i++) {
271 G4double ekin =(*electronCurrentCollection)[i]->GetValue();
272 G4double weight=(*electronCurrentCollection)[i]->GetWeight();
276 for (i=0;i<protonCurrentCollection->
entries();i++) {
277 G4double ekin =(*protonCurrentCollection)[i]->GetValue();
278 G4double weight=(*protonCurrentCollection)[i]->GetWeight();
282 for (i=0;i<gammaCurrentCollection->
entries();i++) {
283 G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
284 G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
320 G4int pdg_nb =theAdjointSimManager
338 G4AnaH2* electron_current_rmatrix =0;
339 G4AnaH2* gamma_current_rmatrix =0;
340 G4AnaH2* proton_current_rmatrix =0;
345 electron_current_rmatrix =
368 for (i=0;i<edepCollection->
entries();i++)
369 totEdep+=(*edepCollection)[i]->GetValue()*
370 (*edepCollection)[i]->GetWeight();
372 G4bool new_mean_computed=
false;
374 if (normalised_weight>0.){
375 G4double edep=totEdep* normalised_weight;
386 if ( new_error >0) new_relative_error = new_error/ new_mean;
388 G4cout<<
"Potential wrong adjoint weight!"<<std::endl;
389 G4cout<<
"The results of this event will not be registered!"
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
404 fRelative_error =new_relative_error;
405 new_mean_computed=
true;
413 edep_rmatrix->fill(prim_ekin,totEdep*adj_weight/
cm2);
415 if (!new_mean_computed){
424 for (i=0;i<electronCurrentCollection->
entries();i++) {
425 G4double ekin =(*electronCurrentCollection)[i]->GetValue();
426 G4double weight=(*electronCurrentCollection)[i]->GetWeight();
428 electron_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/
cm2);
431 for (i=0;i<protonCurrentCollection->
entries();i++) {
432 G4double ekin =(*protonCurrentCollection)[i]->GetValue();
433 G4double weight=(*protonCurrentCollection)[i]->GetWeight();
435 proton_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/
cm2);
438 for (i=0;i<gammaCurrentCollection->
entries();i++) {
439 G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
440 G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
442 gamma_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/
cm2);
519 error = factor*std::sqrt(mean_x2-mean*mean)/std::sqrt(
fNentry);
542 G4cout<<
"The particle that you did select is not in the candidate "<<
543 "list for primary [e-, gamma, proton]!"<<
G4endl;
548 (std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./
pi;
566 G4cout<<
"The particle that you did select is not in the candidate"<<
567 " list for primary [e-, gamma, proton]!"<<
G4endl;
578 -std::pow(Emin,p))/4./
pi;
612 G4cout <<
"\n---> RMC01AnalysisManager::book(): cannot open " <<
fFileName[1]
628 G4String(
"edep vs e- primary energy"),60,emin,emax,
658 G4String(
"electron RM vs e- primary energy"),60,emin,emax,
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,
670 theHistoManager->
GetH2(idHisto);
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,
681 theHistoManager->
GetH2(idHisto);
688 G4String(
"electron RM vs gamma primary energy"),60,emin,emax,
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,
700 theHistoManager->
GetH2(idHisto);
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,
710 theHistoManager->
GetH2(idHisto);
717 G4String(
"electron RM vs proton primary energy"),60,emin,emax,
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,
729 theHistoManager->
GetH2(idHisto);
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,
739 theHistoManager->
GetH2(idHisto);
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,
749 theHistoManager->
GetH2(idHisto);
763 for (
int ind=1; ind<=theHistoManager->
GetNofH1s();ind++){
765 theHistoManager->
ScaleH1(ind,scaling_factor);
767 for (
int ind=1; ind<=theHistoManager->
GetNofH2s();ind++)
768 theHistoManager->
ScaleH2(ind,scaling_factor);
770 theHistoManager->
Write();
774 delete G4AnalysisManager::Instance();
G4AnaH1 * fEdep_rmatrix_vs_proton_prim_energy
static G4AdjointSimManager * GetInstance()
G4bool SetHistoDirectoryName(const G4String &dirName)
virtual void AbortRun(G4bool softAbort=false)
void SetPrimaryPowerLawSpectrumForAdjointSim(const G4String &particle_name, G4double fluence, G4double alpha, G4double Emin, G4double Emax)
G4ThreeVector GetMomentum() const
G4double fEmax_prim_spectrum
G4AnaH1 * fEdep_rmatrix_vs_electron_prim_energy
G4bool SetFirstHistoId(G4int firstId)
G4VHitsCollection * GetHC(G4int i)
G4int GetCollectionID(G4String colName)
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")
G4int fNb_evt_per_adj_evt
G4bool ScaleH1(G4int id, G4double factor)
G4bool ScaleH2(G4int id, G4double factor)
G4bool GetAdjointSimMode()
G4double fError_mean_edep
G4AnaH2 * fGamma_current_rmatrix_vs_proton_prim_energy
G4AnaH2 * fElectron_current_rmatrix_vs_electron_prim_energy
G4int GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack()
G4int GetPDGEncoding() const
G4AnaH2 * fElectron_current_rmatrix_vs_gamma_prim_energy
void SetPrimaryExpSpectrumForAdjointSim(const G4String &particle_name, G4double fluence, G4double E0, G4double Emin, G4double Emax)
void BeginOfEvent(const G4Event *)
PRIM_SPECTRUM_TYPE fPrimSpectrumType
G4AnaH2 * fProton_current_rmatrix_vs_proton_prim_energy
G4bool OpenFile(const G4String &fileName="")
G4ParticleDefinition * GetG4code() const
G4int GetNbEvtOfLastRun()
Definition of the RMC01SD class.
void save(G4double scaling_factor)
G4AnaH1 * fElectron_current
G4bool fStop_run_if_precision_reached
G4double GetEkinAtEndOfLastAdjointTrack()
G4PrimaryParticle * GetPrimary(G4int i=0) const
G4double fPrecision_to_reach
void SetH1Ascii(G4int id, G4bool ascii)
G4GLOB_DLL std::ostream G4cout
G4int GetNumberOfEvent() const
G4AnaH1 * fEdep_vs_prim_ekin
G4String GetFileType() const
G4double fAmplitude_prim_spectrum
static G4Proton * Proton()
void EndOfRun(const G4Run *)
G4PrimaryVertex * GetPrimaryVertex(G4int i=0) const
Definition of the RMC01AnalysisManager class.
G4AnaH1 * fProton_current
G4double PrimDiffAndDirFluxForAdjointSim(G4double prim_energy)
static const G4double emax
RMC01AnalysisManagerMessenger * fMsg
static G4RunManager * GetRunManager()
void ComputeMeanEdepAndError(const G4Event *anEvent, G4double &mean, G4double &error)
Definition of the RMC01AnalysisManagerMessenger class.
static const G4double factor
G4double GetPDGMass() const
G4double GetRealElapsed() const
G4double fAccumulated_edep
G4AnaH2 * fGamma_current_rmatrix_vs_electron_prim_energy
G4AnaH2 * fElectron_current_rmatrix_vs_proton_prim_energy
static RMC01AnalysisManager * GetInstance()
static G4SDManager * GetSDMpointer()
void EndOfEvent(const G4Event *)
G4double GetWeightAtEndOfLastAdjointTrack()
G4int fNb_evt_modulo_for_convergence_test
void BeginOfRun(const G4Run *)
static const G4double Emin
G4int GetNumberOfEventToBeProcessed() const
static G4Electron * Electron()
static const G4double Emax
tools::histo::h1d G4AnaH1
G4double fEmin_prim_spectrum
static PROLOG_HANDLER error
tools::histo::h2d G4AnaH2
G4HCofThisEvent * GetHCofThisEvent() const
static const G4double alpha
G4AnaH2 * fGamma_current_rmatrix_vs_gamma_prim_energy
G4double fAccumulated_edep2
std::fstream fConvergenceFileOutput
void EndOfEventForAdjointSimulation(const G4Event *anEvent)
void EndOfEventForForwardSimulation(const G4Event *anEvent)
static RMC01AnalysisManager * fInstance
G4AnaH1 * fEdep_rmatrix_vs_gamma_prim_energy