Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
RMC01AnalysisManager Class Reference
biasing » ReverseMC01

#include <Doxymodules_biasing.h>

Public Member Functions

 ~RMC01AnalysisManager ()
 
void BeginOfRun (const G4Run *)
 
void EndOfRun (const G4Run *)
 
void BeginOfEvent (const G4Event *)
 
void EndOfEvent (const G4Event *)
 
void SetPrimaryExpSpectrumForAdjointSim (const G4String &particle_name, G4double fluence, G4double E0, G4double Emin, G4double Emax)
 
void SetPrimaryPowerLawSpectrumForAdjointSim (const G4String &particle_name, G4double fluence, G4double alpha, G4double Emin, G4double Emax)
 
void SetPrecision (G4double precision)
 
void Book ()
 
void Save (G4double scaling_factor)
 

Static Public Member Functions

static RMC01AnalysisManager * GetInstance ()
 

Private Member Functions

 RMC01AnalysisManager ()
 
void EndOfEventForForwardSimulation (const G4Event *anEvent)
 
void EndOfEventForAdjointSimulation (const G4Event *anEvent)
 
G4double PrimDiffAndDirFluxForAdjointSim (G4double prim_energy)
 
void ComputeMeanEdepAndError (const G4Event *anEvent, G4double &mean, G4double &error)
 

Private Attributes

RMC01AnalysisManagerMessenger * fMsg
 
G4H1 * fEdep_vs_prim_ekin
 
G4H1 * fElectron_current
 
G4H1 * fProton_current
 
G4H1 * fGamma_current
 
G4double fAccumulated_edep
 
G4double fAccumulated_edep2
 
G4double fNentry
 
G4double fMean_edep
 
G4double fError_mean_edep
 
G4double fRelative_error
 
G4double fElapsed_time
 
G4double fPrecision_to_reach
 
G4bool fStop_run_if_precision_reached
 
G4int fNb_evt_modulo_for_convergence_test
 
G4H1 * fEdep_rmatrix_vs_electron_prim_energy
 
G4H2 * fElectron_current_rmatrix_vs_electron_prim_energy
 
G4H2 * fGamma_current_rmatrix_vs_electron_prim_energy
 
G4H1 * fEdep_rmatrix_vs_gamma_prim_energy
 
G4H2 * fElectron_current_rmatrix_vs_gamma_prim_energy
 
G4H2 * fGamma_current_rmatrix_vs_gamma_prim_energy
 
G4H1 * fEdep_rmatrix_vs_proton_prim_energy
 
G4H2 * fElectron_current_rmatrix_vs_proton_prim_energy
 
G4H2 * fProton_current_rmatrix_vs_proton_prim_energy
 
G4H2 * fGamma_current_rmatrix_vs_proton_prim_energy
 
G4String fFileName [2]
 
G4bool fFactoryOn
 
PRIM_SPECTRUM_TYPE fPrimSpectrumType
 
G4int fPrimPDG_ID
 
G4double fAlpha_or_E0
 
G4double fAmplitude_prim_spectrum
 
G4double fEmin_prim_spectrum
 
G4double fEmax_prim_spectrum
 
G4bool fAdjoint_sim_mode
 
G4int fNb_evt_per_adj_evt
 
G4Timer * fTimer
 
std::fstream fConvergenceFileOutput
 

Static Private Attributes

static RMC01AnalysisManager * fInstance = 0
 

Detailed Description

Definition at line 170 of file Doxymodules_biasing.h.

Constructor & Destructor Documentation

◆ ~RMC01AnalysisManager()

RMC01AnalysisManager::~RMC01AnalysisManager ( )

Definition at line 101 of file RMC01AnalysisManager.cc.

102{;
103}

◆ RMC01AnalysisManager()

RMC01AnalysisManager::RMC01AnalysisManager ( )
private

Definition at line 59 of file RMC01AnalysisManager.cc.

60 :fAccumulated_edep(0.), fAccumulated_edep2(0.), fMean_edep(0.),
61 fError_mean_edep(0.), fRelative_error(0.), fElapsed_time(0.),
62 fPrecision_to_reach(0.),fStop_run_if_precision_reached(true),
63 fNb_evt_modulo_for_convergence_test(5000),
64 fEdep_rmatrix_vs_electron_prim_energy(0),
65 fElectron_current_rmatrix_vs_electron_prim_energy(0),
66 fGamma_current_rmatrix_vs_electron_prim_energy(0),
67 fEdep_rmatrix_vs_gamma_prim_energy(0),
68 fElectron_current_rmatrix_vs_gamma_prim_energy(0),
69 fGamma_current_rmatrix_vs_gamma_prim_energy(0),
70 fEdep_rmatrix_vs_proton_prim_energy(0),
71 fElectron_current_rmatrix_vs_proton_prim_energy(0),
72 fProton_current_rmatrix_vs_proton_prim_energy(0),
73 fGamma_current_rmatrix_vs_proton_prim_energy(0),
74 fFactoryOn(false),
75 fPrimSpectrumType(EXPO),
76 fAlpha_or_E0(.5*MeV),fAmplitude_prim_spectrum (1.),
77 fEmin_prim_spectrum(1.*keV),fEmax_prim_spectrum (20.*MeV),
78 fAdjoint_sim_mode(true),fNb_evt_per_adj_evt(2)
79{
80
81 fMsg = new RMC01AnalysisManagerMessenger(this);
82
83 //-------------
84 //Timer for convergence vector
85 //-------------
86
87 fTimer = new G4Timer();
88
89 //---------------------------------
90 //Primary particle ID for normalisation of adjoint results
91 //---------------------------------
92
93 fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
94
95 fFileName[0] = "sim";
96
97}
RMC01AnalysisManager::fGamma_current_rmatrix_vs_electron_prim_energy
G4H2 * fGamma_current_rmatrix_vs_electron_prim_energy
Definition RMC01AnalysisManager.hh:153
RMC01AnalysisManager::fMsg
RMC01AnalysisManagerMessenger * fMsg
Definition RMC01AnalysisManager.hh:121
RMC01AnalysisManager::fGamma_current_rmatrix_vs_proton_prim_energy
G4H2 * fGamma_current_rmatrix_vs_proton_prim_energy
Definition RMC01AnalysisManager.hh:162
RMC01AnalysisManager::fElectron_current_rmatrix_vs_proton_prim_energy
G4H2 * fElectron_current_rmatrix_vs_proton_prim_energy
Definition RMC01AnalysisManager.hh:160
RMC01AnalysisManager::fPrimSpectrumType
PRIM_SPECTRUM_TYPE fPrimSpectrumType
Definition RMC01AnalysisManager.hh:172
RMC01AnalysisManager::fElectron_current_rmatrix_vs_gamma_prim_energy
G4H2 * fElectron_current_rmatrix_vs_gamma_prim_energy
Definition RMC01AnalysisManager.hh:156
RMC01AnalysisManager::fProton_current_rmatrix_vs_proton_prim_energy
G4H2 * fProton_current_rmatrix_vs_proton_prim_energy
Definition RMC01AnalysisManager.hh:161
RMC01AnalysisManager::fElectron_current_rmatrix_vs_electron_prim_energy
G4H2 * fElectron_current_rmatrix_vs_electron_prim_energy
Definition RMC01AnalysisManager.hh:152

Member Function Documentation

◆ GetInstance()

RMC01AnalysisManager * RMC01AnalysisManager::GetInstance ( )
static

Definition at line 107 of file RMC01AnalysisManager.cc.

108{
109 if (fInstance == 0) fInstance = new RMC01AnalysisManager;
110 return fInstance;
111}
RMC01AnalysisManager::fInstance
static RMC01AnalysisManager * fInstance
Definition RMC01AnalysisManager.hh:105

◆ BeginOfRun()

void RMC01AnalysisManager::BeginOfRun ( const G4Run *  aRun)

Definition at line 115 of file RMC01AnalysisManager.cc.

116{
117
118 fAccumulated_edep =0.;
119 fAccumulated_edep2 =0.;
120 fNentry = 0.0;
121 fRelative_error=1.;
122 fMean_edep=0.;
123 fError_mean_edep=0.;
124 fAdjoint_sim_mode =G4AdjointSimManager::GetInstance()->GetAdjointSimMode();
125
126 if (fAdjoint_sim_mode){
127 fNb_evt_per_adj_evt=aRun->GetNumberOfEventToBeProcessed()/
128 G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
129 fConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt",
130 std::ios::out);
131 fConvergenceFileOutput<<
132 "Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
133 }
134 else {
135 fConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt",
136 std::ios::out);
137 fConvergenceFileOutput<<
138 "Edep per event [MeV]\terror[MeV]\tcomputing_time[s]"
139 <<std::endl;
140 }
141 fConvergenceFileOutput.setf(std::ios::scientific);
142 fConvergenceFileOutput.precision(6);
143
144 fTimer->Start();
145 fElapsed_time=0.;
146
147 Book();
148}

◆ EndOfRun()

void RMC01AnalysisManager::EndOfRun ( const G4Run *  aRun)

Definition at line 152 of file RMC01AnalysisManager.cc.

153{ fTimer->Stop();
154 G4int nb_evt=aRun->GetNumberOfEvent();
155 G4double factor =1./ nb_evt;
156 if (!fAdjoint_sim_mode){
157 G4cout<<"Results of forward simulation!"<<std::endl;
158 G4cout<<"edep per event [MeV] = "<<fMean_edep<<std::endl;
159 G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
160 }
161
162 else {
163 G4cout<<"Results of reverse/adjoint simulation!"<<std::endl;
164 G4cout<<"normalised edep [MeV] = "<<fMean_edep<<std::endl;
165 G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
166 factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun()
167 *fNb_evt_per_adj_evt/aRun->GetNumberOfEvent();
168 }
169 Save(factor);
170 fConvergenceFileOutput.close();
171}
RMC01AnalysisManager::Save
void Save(G4double scaling_factor)
Definition RMC01AnalysisManager.cc:767

◆ BeginOfEvent()

void RMC01AnalysisManager::BeginOfEvent ( const G4Event *  )

Definition at line 175 of file RMC01AnalysisManager.cc.

176{ ;
177}

◆ EndOfEvent()

void RMC01AnalysisManager::EndOfEvent ( const G4Event *  anEvent)

Definition at line 181 of file RMC01AnalysisManager.cc.

182{
183 if (fAdjoint_sim_mode) EndOfEventForAdjointSimulation(anEvent);
184 else EndOfEventForForwardSimulation(anEvent);
185
186 //Test convergence. The error is already computed
187 //--------------------------------------
188 G4int nb_event=anEvent->GetEventID()+1;
189 //G4double factor=1.;
190 if (fAdjoint_sim_mode) {
191 G4double n_adj_evt= nb_event/fNb_evt_per_adj_evt;
192 // nb_event/fNb_evt_per_adj_evt;
193 if (n_adj_evt*fNb_evt_per_adj_evt == nb_event) {
194 nb_event =static_cast<G4int>(n_adj_evt);
195 }
196 else nb_event=0;
197 }
198
199 if (nb_event>100 && fStop_run_if_precision_reached &&
200 fPrecision_to_reach >fRelative_error) {
201 G4cout<<fPrecision_to_reach*100.<<"% Precision reached!"<<std::endl;
202 fTimer->Stop();
203 fElapsed_time+=fTimer->GetRealElapsed();
204 fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep
205 <<'\t'<<fElapsed_time<<std::endl;
206 G4RunManager::GetRunManager()->AbortRun(true);
207 }
208
209 if (nb_event>0 && nb_event % fNb_evt_modulo_for_convergence_test == 0) {
210 fTimer->Stop();
211 fElapsed_time+=fTimer->GetRealElapsed();
212 fTimer->Start();
213 fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep<<'\t'
214 <<fElapsed_time<<std::endl;
215 }
216}
RMC01AnalysisManager::EndOfEventForForwardSimulation
void EndOfEventForForwardSimulation(const G4Event *anEvent)
Definition RMC01AnalysisManager.cc:220
RMC01AnalysisManager::EndOfEventForAdjointSimulation
void EndOfEventForAdjointSimulation(const G4Event *anEvent)
Definition RMC01AnalysisManager.cc:289

◆ SetPrimaryExpSpectrumForAdjointSim()

void RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim ( const G4String &  particle_name,
G4double  fluence,
G4double  E0,
G4double  Emin,
G4double  Emax 
)

Definition at line 545 of file RMC01AnalysisManager.cc.

549{ fPrimSpectrumType = EXPO;
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 fAlpha_or_E0 = E0 ;
562 fAmplitude_prim_spectrum = omni_fluence/E0/
563 (std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
564 fEmin_prim_spectrum = Emin ;
565 fEmax_prim_spectrum = Emax;
566}

◆ SetPrimaryPowerLawSpectrumForAdjointSim()

void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim ( const G4String &  particle_name,
G4double  fluence,
G4double  alpha,
G4double  Emin,
G4double  Emax 
)

Definition at line 570 of file RMC01AnalysisManager.cc.

573{ fPrimSpectrumType =POWER;
574 if (particle_name == "e-" ) fPrimPDG_ID =
575 G4Electron::Electron()->GetPDGEncoding();
576 else if (particle_name == "gamma") fPrimPDG_ID =
577 G4Gamma::Gamma()->GetPDGEncoding();
578 else if (particle_name == "proton") fPrimPDG_ID =
579 G4Proton::Proton()->GetPDGEncoding();
580 else {
581 G4cout<<"The particle that you did select is not in the candidate"<<
582 " list for primary [e-, gamma, proton]!"<<G4endl;
583 return;
584 }
585
586
587 if (alpha ==1.) {
588 fAmplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
589 }
590 else {
591 G4double p=1.-alpha;
592 fAmplitude_prim_spectrum = omni_fluence/p/(std::pow(Emax,p)
593 -std::pow(Emin,p))/4./pi;
594 }
595
596 fAlpha_or_E0 = alpha ;
597 fEmin_prim_spectrum = Emin ;
598 fEmax_prim_spectrum = Emax;
599
600}

◆ SetPrecision()

void RMC01AnalysisManager::SetPrecision ( G4double  precision)
inline

Definition at line 96 of file RMC01AnalysisManager.hh.

96 {
97 fPrecision_to_reach =precision/100.;};

◆ Book()

void RMC01AnalysisManager::Book ( )

Definition at line 604 of file RMC01AnalysisManager.cc.

605{
606 //----------------------
607 //Creation of histograms
608 //----------------------
609
610 //Energy binning of the histograms : 60 log bins over [1keV-1GeV]
611
612 G4double emin=1.*keV;
613 G4double emax=1.*GeV;
614
615 //file_name
616 fFileName[0]="forward_sim";
617 if (fAdjoint_sim_mode) fFileName[0]="adjoint_sim";
618
619 //Histo manager
620 G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
621 theHistoManager->SetDefaultFileType("root");
622 G4String extension = theHistoManager->GetFileType();
623 fFileName[1] = fFileName[0] + "." + extension;
624 theHistoManager->SetFirstHistoId(1);
625
626 G4bool fileOpen = theHistoManager->OpenFile(fFileName[0]);
627 if (!fileOpen) {
628 G4cout << "\n---> RMC01AnalysisManager::Book(): cannot open "
629 << fFileName[1]
630 << G4endl;
631 return;
632 }
633
634 // Create directories
635 theHistoManager->SetHistoDirectoryName("histo");
636
637 //Histograms for :
638 // 1)the forward simulation results
639 // 2)the Reverse MC simulation results normalised to a user spectrum
640 //------------------------------------------------------------------------
641
642 G4int idHisto =
643 theHistoManager->CreateH1(G4String("Edep_vs_prim_ekin"),
644 G4String("edep vs e- primary energy"),60,emin,emax,
645 "none","none",G4String("log"));
646 fEdep_vs_prim_ekin = theHistoManager->GetH1(idHisto);
647
648 idHisto = theHistoManager->CreateH1(G4String("elecron_current"),
649 G4String("electron"),60,emin,emax,
650 "none","none",G4String("log"));
651
652 fElectron_current = theHistoManager->GetH1(idHisto);
653
654 idHisto= theHistoManager->CreateH1(G4String("proton_current"),
655 G4String("proton"),60,emin,emax,
656 "none","none",G4String("log"));
657 fProton_current=theHistoManager->GetH1(idHisto);
658
659 idHisto= theHistoManager->CreateH1(G4String("gamma_current"),
660 G4String("gamma"),60,emin,emax,
661 "none","none",G4String("log"));
662 fGamma_current=theHistoManager->GetH1(idHisto);
663
664 //Response matrices for the adjoint simulation only
665 //-----------------------------------------------
666 if (fAdjoint_sim_mode){
667 //Response matrices for external isotropic e- source
668 //--------------------------------------------------
669
670 idHisto =
671 theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_electron_prim_energy"),
672 G4String("electron RM vs e- primary energy"),60,emin,emax,
673 "none","none",G4String("log"));
674 fEdep_rmatrix_vs_electron_prim_energy = theHistoManager->GetH1(idHisto);
675
676 idHisto =
677 theHistoManager->
678 CreateH2(G4String("Electron_current_rmatrix_vs_electron_prim_energy"),
679 G4String("electron current RM vs e- primary energy"),
680 60,emin,emax,60,emin,emax,
681 "none","none","none","none",G4String("log"),G4String("log"));
682
683 fElectron_current_rmatrix_vs_electron_prim_energy =
684 theHistoManager->GetH2(idHisto);
685
686 idHisto =
687 theHistoManager->
688 CreateH2(G4String("Gamma_current_rmatrix_vs_electron_prim_energy"),
689 G4String("gamma current RM vs e- primary energy"),
690 60,emin,emax,60,emin,emax,
691 "none","none","none","none",G4String("log"),G4String("log"));
692
693 fGamma_current_rmatrix_vs_electron_prim_energy =
694 theHistoManager->GetH2(idHisto);
695
696 //Response matrices for external isotropic gamma source
697
698 idHisto =
699 theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_gamma_prim_energy"),
700 G4String("electron RM vs gamma primary energy"),60,emin,emax,
701 "none","none",G4String("log"));
702 fEdep_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH1(idHisto);
703
704 idHisto =
705 theHistoManager->
706 CreateH2(G4String("Electron_current_rmatrix_vs_gamma_prim_energy"),
707 G4String("electron current RM vs gamma primary energy"),
708 60,emin,emax,60,emin,emax,
709 "none","none","none","none",G4String("log"),G4String("log"));
710
711 fElectron_current_rmatrix_vs_gamma_prim_energy =
712 theHistoManager->GetH2(idHisto);
713
714 idHisto =
715 theHistoManager->
716 CreateH2(G4String("Gamma_current_rmatrix_vs_gamma_prim_energy"),
717 G4String("gamma current RM vs gamma primary energy"),
718 60,emin,emax,60,emin,emax,
719 "none","none","none","none",G4String("log"),G4String("log"));
720
721 fGamma_current_rmatrix_vs_gamma_prim_energy =
722 theHistoManager->GetH2(idHisto);
723
724 //Response matrices for external isotropic proton source
725 idHisto =
726 theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_proton_prim_energy"),
727 G4String("electron RM vs proton primary energy"),60,emin,emax,
728 "none","none",G4String("log"));
729 fEdep_rmatrix_vs_proton_prim_energy = theHistoManager->GetH1(idHisto);
730
731 idHisto =
732 theHistoManager->
733 CreateH2(G4String("Electron_current_rmatrix_vs_proton_prim_energy"),
734 G4String("electron current RM vs proton primary energy"),
735 60,emin,emax,60,emin,emax,
736 "none","none","none","none",G4String("log"),G4String("log"));
737
738 fElectron_current_rmatrix_vs_proton_prim_energy =
739 theHistoManager->GetH2(idHisto);
740
741 idHisto =
742 theHistoManager->
743 CreateH2(G4String("Gamma_current_rmatrix_vs_proton_prim_energy"),
744 G4String("gamma current RM vs proton primary energy"),
745 60,emin,emax,60,emin,emax,
746 "none","none","none","none",G4String("log"),G4String("log"));
747
748 fGamma_current_rmatrix_vs_proton_prim_energy =
749 theHistoManager->GetH2(idHisto);
750
751 idHisto =
752 theHistoManager->
753 CreateH2(G4String("Proton_current_rmatrix_vs_proton_prim_energy"),
754 G4String("proton current RM vs proton primary energy"),
755 60,emin,emax,60,emin,emax,
756 "none","none","none","none",G4String("log"),G4String("log"));
757
758 fProton_current_rmatrix_vs_proton_prim_energy =
759 theHistoManager->GetH2(idHisto);
760 }
761 fFactoryOn = true;
762 G4cout << "\n----> Histogram Tree is opened in " << fFileName[1] << G4endl;
763}

◆ Save()

void RMC01AnalysisManager::Save ( G4double  scaling_factor)

Definition at line 767 of file RMC01AnalysisManager.cc.

768{ if (fFactoryOn) {
769 G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
770 //scaling of results
771 //-----------------
772
773 for (G4int ind=1; ind<=theHistoManager->GetNofH1s();++ind){
774 theHistoManager->SetH1Ascii(ind,true);
775 theHistoManager->ScaleH1(ind,scaling_factor);
776 }
777 for (G4int ind=1; ind<=theHistoManager->GetNofH2s();++ind)
778 theHistoManager->ScaleH2(ind,scaling_factor);
779
780 theHistoManager->Write();
781 theHistoManager->CloseFile();
782 G4cout << "\n----> Histogram Tree is saved in " << fFileName[1] << G4endl;
783
784 theHistoManager->Clear();
785 fFactoryOn = false;
786 }
787}

◆ EndOfEventForForwardSimulation()

void RMC01AnalysisManager::EndOfEventForForwardSimulation ( const G4Event *  anEvent)
private

Definition at line 220 of file RMC01AnalysisManager.cc.

222{
223
224 G4SDManager* SDman = G4SDManager::GetSDMpointer();
225 G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
226 RMC01DoubleWithWeightHitsCollection* edepCollection =
227 (RMC01DoubleWithWeightHitsCollection*)
228 (HCE->GetHC(SDman->GetCollectionID("edep")));
229
230 RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
231 (RMC01DoubleWithWeightHitsCollection*)
232 (HCE->GetHC(SDman->GetCollectionID("current_electron")));
233
234 RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
235 (RMC01DoubleWithWeightHitsCollection*)
236 (HCE->GetHC(SDman->GetCollectionID("current_proton")));
237
238 RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
239 (RMC01DoubleWithWeightHitsCollection*)
240 (HCE->GetHC(SDman->GetCollectionID("current_gamma")));
241
242 //Total energy deposited in Event
243 //-------------------------------
244 G4double totEdep=0;
245 std::size_t i;
246 for (i=0;i<edepCollection->entries();++i)
247 totEdep+=(*edepCollection)[i]->GetValue()
248 *(*edepCollection)[i]->GetWeight();
249
250 if (totEdep>0.){
251 fAccumulated_edep +=totEdep ;
252 fAccumulated_edep2 +=totEdep*totEdep;
253 fNentry += 1.0;
254 G4PrimaryParticle* thePrimary=
255 anEvent->GetPrimaryVertex()->GetPrimary();
256 G4double E0= thePrimary->GetG4code()->GetPDGMass();
257 G4double P=thePrimary->GetMomentum().mag();
258 G4double prim_ekin =std::sqrt(E0*E0+P*P)-E0;
259 fEdep_vs_prim_ekin->fill(prim_ekin,totEdep);
260 }
261 ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
262 if (fError_mean_edep>0) fRelative_error= fError_mean_edep/fMean_edep;
263
264 //Particle current on sensitive cylinder
265 //-------------------------------------
266
267 for (i=0;i<electronCurrentCollection->entries();++i) {
268 G4double ekin =(*electronCurrentCollection)[i]->GetValue();
269 G4double weight=(*electronCurrentCollection)[i]->GetWeight();
270 fElectron_current->fill(ekin,weight);
271 }
272
273 for (i=0;i<protonCurrentCollection->entries();++i) {
274 G4double ekin =(*protonCurrentCollection)[i]->GetValue();
275 G4double weight=(*protonCurrentCollection)[i]->GetWeight();
276 fProton_current->fill(ekin,weight);
277 }
278
279 for (i=0;i<gammaCurrentCollection->entries();++i) {
280 G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
281 G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
282 fGamma_current->fill(ekin,weight);
283 }
284
285}
RMC01DoubleWithWeightHitsCollection
G4THitsCollection< RMC01DoubleWithWeightHit > RMC01DoubleWithWeightHitsCollection
Definition RMC01DoubleWithWeightHit.hh:87
RMC01AnalysisManager::ComputeMeanEdepAndError
void ComputeMeanEdepAndError(const G4Event *anEvent, G4double &mean, G4double &error)
Definition RMC01AnalysisManager.cc:508

◆ EndOfEventForAdjointSimulation()

void RMC01AnalysisManager::EndOfEventForAdjointSimulation ( const G4Event *  anEvent)
private

Definition at line 289 of file RMC01AnalysisManager.cc.

291{
292 //Output from Sensitive volume computed during the forward tracking phase
293 //-----------------------------------------------------------------------
294 G4SDManager* SDman = G4SDManager::GetSDMpointer();
295 G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
296 RMC01DoubleWithWeightHitsCollection* edepCollection =
297 (RMC01DoubleWithWeightHitsCollection*)(
298 HCE->GetHC(SDman->GetCollectionID("edep")));
299
300 RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
301 (RMC01DoubleWithWeightHitsCollection*)(
302 HCE->GetHC(SDman->GetCollectionID("current_electron")));
303
304 RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
305 (RMC01DoubleWithWeightHitsCollection*)(
306 HCE->GetHC(SDman->GetCollectionID("current_proton")));
307
308 RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
309 (RMC01DoubleWithWeightHitsCollection*)(
310 HCE->GetHC(SDman->GetCollectionID("current_gamma")));
311
312 //Computation of total energy deposited in fwd tracking phase
313 //-------------------------------
314 G4double totEdep=0;
315 std::size_t i;
316 for (i=0;i<edepCollection->entries();++i)
317 totEdep+=(*edepCollection)[i]->GetValue()*
318 (*edepCollection)[i]->GetWeight();
319
320 //Output from adjoint tracking phase
321 //----------------------------------------------------------------------------
322
323 G4AdjointSimManager* theAdjointSimManager =
324 G4AdjointSimManager::GetInstance();
325
326 size_t nb_adj_track =
327 theAdjointSimManager->GetNbOfAdointTracksReachingTheExternalSurface();
328 G4double total_normalised_weight = 0.;
329
330 //We need to loop over the adjoint tracks that have reached the external
331 //surface.
332 for (std::size_t j=0;j<nb_adj_track;++j) {
333
334 G4int pdg_nb =theAdjointSimManager
335 ->GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack(j);
336 G4double prim_ekin=theAdjointSimManager
337 ->GetEkinAtEndOfLastAdjointTrack(j);
338 G4double adj_weight=theAdjointSimManager
339 ->GetWeightAtEndOfLastAdjointTrack(j);
340
341
342 //Factor of normalisation to user defined prim spectrum (power law or exp)
343 //------------------------------------------------------------------------
344 G4double normalised_weight = 0.;
345 if (pdg_nb== fPrimPDG_ID && prim_ekin>= fEmin_prim_spectrum
346 && prim_ekin<= fEmax_prim_spectrum)
347 normalised_weight =
348 adj_weight*PrimDiffAndDirFluxForAdjointSim(prim_ekin);
349 total_normalised_weight += normalised_weight;
350
351 //Answer matrices
352 //-------------
353 G4H1* edep_rmatrix =0;
354 G4H2* electron_current_rmatrix =0;
355 G4H2* gamma_current_rmatrix =0;
356 G4H2* proton_current_rmatrix =0;
357
358 if (pdg_nb == G4Electron::Electron()->GetPDGEncoding()){ //e- matrices
359 edep_rmatrix = fEdep_rmatrix_vs_electron_prim_energy;
360 electron_current_rmatrix =
361 fElectron_current_rmatrix_vs_electron_prim_energy;
362 gamma_current_rmatrix = fGamma_current_rmatrix_vs_electron_prim_energy;
363 }
364 else if (pdg_nb == G4Gamma::Gamma()->GetPDGEncoding()){
365 //gammma answer matrices
366 edep_rmatrix = fEdep_rmatrix_vs_gamma_prim_energy;
367 electron_current_rmatrix = fElectron_current_rmatrix_vs_gamma_prim_energy;
368 gamma_current_rmatrix = fGamma_current_rmatrix_vs_gamma_prim_energy;
369 }
370 else if (pdg_nb == G4Proton::Proton()->GetPDGEncoding()){
371 //proton answer matrices
372 edep_rmatrix = fEdep_rmatrix_vs_proton_prim_energy;
373 electron_current_rmatrix =
374 fElectron_current_rmatrix_vs_proton_prim_energy;
375 gamma_current_rmatrix = fGamma_current_rmatrix_vs_proton_prim_energy;
376 proton_current_rmatrix = fProton_current_rmatrix_vs_proton_prim_energy;
377 }
378 //Register histo edep vs prim ekin
379 //----------------------------------
380 if (normalised_weight>0) fEdep_vs_prim_ekin
381 ->fill(prim_ekin,totEdep*normalised_weight);
382 // Registering answer matrix
383 //---------------------------
384 edep_rmatrix->fill(prim_ekin,totEdep*adj_weight/cm2);
385
386 //Registering of current of particles on the sensitive volume
387 //------------------------------------------------------------
388
389 for (i=0;i<electronCurrentCollection->entries();++i) {
390 G4double ekin =(*electronCurrentCollection)[i]->GetValue();
391 G4double weight=(*electronCurrentCollection)[i]->GetWeight();
392 fElectron_current->fill(ekin,weight*normalised_weight);
393 electron_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
394 }
395 for (i=0;i<protonCurrentCollection->entries();++i) {
396 G4double ekin =(*protonCurrentCollection)[i]->GetValue();
397 G4double weight=(*protonCurrentCollection)[i]->GetWeight();
398 fProton_current->fill(ekin,weight*normalised_weight);
399 proton_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
400 }
401 for (i=0;i<gammaCurrentCollection->entries();++i) {
402 G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
403 G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
404 fGamma_current->fill(ekin,weight*normalised_weight);
405 gamma_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
406 }
407 }
408
409 //Registering of total energy deposited in Event
410 //-------------------------------
411 G4bool new_mean_computed=false;
412 if (totEdep>0.){
413 if (total_normalised_weight>0.){
414 G4double edep=totEdep* total_normalised_weight;
415
416 //Check if the edep is not wrongly too high
417 //-----------------------------------------
418 G4double new_mean(0.0), new_error(0.0);
419 fAccumulated_edep +=edep;
420 fAccumulated_edep2 +=edep*edep;
421 fNentry += 1.0;
422 ComputeMeanEdepAndError(anEvent,new_mean,new_error);
423 G4double new_relative_error = 1.;
424 if ( new_error >0) new_relative_error = new_error/ new_mean;
425 if (fRelative_error <0.10 && new_relative_error>1.5*fRelative_error) {
426 G4cout<<"Potential wrong adjoint weight!"<<std::endl;
427 G4cout<<"The results of this event will not be registered!"
428 <<std::endl;
429 G4cout<<"previous mean edep [MeV] "<< fMean_edep<<std::endl;
430 G4cout<<"previous relative error "<< fRelative_error<<std::endl;
431 G4cout<<"new rejected mean edep [MeV] "<< new_mean<<std::endl;
432 G4cout<<"new rejected relative error "<< new_relative_error
433 <<std::endl;
434 fAccumulated_edep -=edep;
435 fAccumulated_edep2 -=edep*edep;
436 fNentry -= 1.0;
437 return;
438 }
439 else { //accepted
440 fMean_edep = new_mean;
441 fError_mean_edep = new_error;
442 fRelative_error =new_relative_error;
443 new_mean_computed=true;
444 }
445
446 }
447
448 if (!new_mean_computed){
449 ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
450 fRelative_error = (fMean_edep > 0.0) ? fError_mean_edep/fMean_edep : 0.0;
451 }
452 }
453}
RMC01AnalysisManager::PrimDiffAndDirFluxForAdjointSim
G4double PrimDiffAndDirFluxForAdjointSim(G4double prim_energy)
Definition RMC01AnalysisManager.cc:457

◆ PrimDiffAndDirFluxForAdjointSim()

G4double RMC01AnalysisManager::PrimDiffAndDirFluxForAdjointSim ( G4double  prim_energy)
private

Definition at line 457 of file RMC01AnalysisManager.cc.

459{
460 G4double flux=fAmplitude_prim_spectrum;
461 if ( fPrimSpectrumType ==EXPO) flux*=std::exp(-prim_energy/fAlpha_or_E0);
462 else flux*=std::pow(prim_energy, -fAlpha_or_E0);
463 return flux;
464}

◆ ComputeMeanEdepAndError()

void RMC01AnalysisManager::ComputeMeanEdepAndError ( const G4Event *  anEvent,
G4double &  mean,
G4double &  error 
)
private

Definition at line 508 of file RMC01AnalysisManager.cc.

510{
511 G4int nb_event = anEvent->GetEventID()+1;
512 G4double factor=1.;
513 if (fAdjoint_sim_mode) {
514 nb_event /= fNb_evt_per_adj_evt;
515 factor = 1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
516 }
517
518 // VI: error computation now is based on number of entries and not
519 // number of events
520 // LD: This is wrong! With the use of fNentry the results were no longer
521 // correctly normalised. The mean and the error should be computed
522 // with nb_event. The old computation has been reset.
523 // VI: OK, but let computations be double
524 if (nb_event > 0) {
525 G4double norm = 1.0/(G4double)nb_event;
526 mean = fAccumulated_edep*norm;
527 G4double mean_x2 = fAccumulated_edep2*norm;
528 G4double zz = mean_x2 - mean*mean;
529 /*
530 G4cout << "Nevt= " << nb_event << " mean= " << mean
531 << " mean_x2= " << mean_x2 << " x2 - x*x= "
532 << zz << G4endl;
533 */
534 error = factor*std::sqrt(std::max(zz, 0.)*norm);
535 mean *= factor;
536 } else {
537 mean=0;
538 error=0;
539 }
540 //G4cout << "Aend: " << mean << " " << error << G4endl;
541}

Member Data Documentation

◆ fInstance

RMC01AnalysisManager * RMC01AnalysisManager::fInstance = 0
staticprivate

Definition at line 105 of file RMC01AnalysisManager.hh.

◆ fMsg

RMC01AnalysisManagerMessenger* RMC01AnalysisManager::fMsg
private

Definition at line 121 of file RMC01AnalysisManager.hh.

◆ fEdep_vs_prim_ekin

G4H1* RMC01AnalysisManager::fEdep_vs_prim_ekin
private

Definition at line 125 of file RMC01AnalysisManager.hh.

◆ fElectron_current

G4H1* RMC01AnalysisManager::fElectron_current
private

Definition at line 126 of file RMC01AnalysisManager.hh.

◆ fProton_current

G4H1* RMC01AnalysisManager::fProton_current
private

Definition at line 127 of file RMC01AnalysisManager.hh.

◆ fGamma_current

G4H1* RMC01AnalysisManager::fGamma_current
private

Definition at line 128 of file RMC01AnalysisManager.hh.

◆ fAccumulated_edep

G4double RMC01AnalysisManager::fAccumulated_edep
private

Definition at line 137 of file RMC01AnalysisManager.hh.

◆ fAccumulated_edep2

G4double RMC01AnalysisManager::fAccumulated_edep2
private

Definition at line 138 of file RMC01AnalysisManager.hh.

◆ fNentry

G4double RMC01AnalysisManager::fNentry
private

Definition at line 139 of file RMC01AnalysisManager.hh.

◆ fMean_edep

G4double RMC01AnalysisManager::fMean_edep
private

Definition at line 140 of file RMC01AnalysisManager.hh.

◆ fError_mean_edep

G4double RMC01AnalysisManager::fError_mean_edep
private

Definition at line 141 of file RMC01AnalysisManager.hh.

◆ fRelative_error

G4double RMC01AnalysisManager::fRelative_error
private

Definition at line 142 of file RMC01AnalysisManager.hh.

◆ fElapsed_time

G4double RMC01AnalysisManager::fElapsed_time
private

Definition at line 143 of file RMC01AnalysisManager.hh.

◆ fPrecision_to_reach

G4double RMC01AnalysisManager::fPrecision_to_reach
private

Definition at line 144 of file RMC01AnalysisManager.hh.

◆ fStop_run_if_precision_reached

G4bool RMC01AnalysisManager::fStop_run_if_precision_reached
private

Definition at line 145 of file RMC01AnalysisManager.hh.

◆ fNb_evt_modulo_for_convergence_test

G4int RMC01AnalysisManager::fNb_evt_modulo_for_convergence_test
private

Definition at line 146 of file RMC01AnalysisManager.hh.

◆ fEdep_rmatrix_vs_electron_prim_energy

G4H1* RMC01AnalysisManager::fEdep_rmatrix_vs_electron_prim_energy
private

Definition at line 151 of file RMC01AnalysisManager.hh.

◆ fElectron_current_rmatrix_vs_electron_prim_energy

G4H2* RMC01AnalysisManager::fElectron_current_rmatrix_vs_electron_prim_energy
private

Definition at line 152 of file RMC01AnalysisManager.hh.

◆ fGamma_current_rmatrix_vs_electron_prim_energy

G4H2* RMC01AnalysisManager::fGamma_current_rmatrix_vs_electron_prim_energy
private

Definition at line 153 of file RMC01AnalysisManager.hh.

◆ fEdep_rmatrix_vs_gamma_prim_energy

G4H1* RMC01AnalysisManager::fEdep_rmatrix_vs_gamma_prim_energy
private

Definition at line 155 of file RMC01AnalysisManager.hh.

◆ fElectron_current_rmatrix_vs_gamma_prim_energy

G4H2* RMC01AnalysisManager::fElectron_current_rmatrix_vs_gamma_prim_energy
private

Definition at line 156 of file RMC01AnalysisManager.hh.

◆ fGamma_current_rmatrix_vs_gamma_prim_energy

G4H2* RMC01AnalysisManager::fGamma_current_rmatrix_vs_gamma_prim_energy
private

Definition at line 157 of file RMC01AnalysisManager.hh.

◆ fEdep_rmatrix_vs_proton_prim_energy

G4H1* RMC01AnalysisManager::fEdep_rmatrix_vs_proton_prim_energy
private

Definition at line 159 of file RMC01AnalysisManager.hh.

◆ fElectron_current_rmatrix_vs_proton_prim_energy

G4H2* RMC01AnalysisManager::fElectron_current_rmatrix_vs_proton_prim_energy
private

Definition at line 160 of file RMC01AnalysisManager.hh.

◆ fProton_current_rmatrix_vs_proton_prim_energy

G4H2* RMC01AnalysisManager::fProton_current_rmatrix_vs_proton_prim_energy
private

Definition at line 161 of file RMC01AnalysisManager.hh.

◆ fGamma_current_rmatrix_vs_proton_prim_energy

G4H2* RMC01AnalysisManager::fGamma_current_rmatrix_vs_proton_prim_energy
private

Definition at line 162 of file RMC01AnalysisManager.hh.

◆ fFileName

G4String RMC01AnalysisManager::fFileName[2]
private

Definition at line 164 of file RMC01AnalysisManager.hh.

◆ fFactoryOn

G4bool RMC01AnalysisManager::fFactoryOn
private

Definition at line 165 of file RMC01AnalysisManager.hh.

◆ fPrimSpectrumType

PRIM_SPECTRUM_TYPE RMC01AnalysisManager::fPrimSpectrumType
private

Definition at line 172 of file RMC01AnalysisManager.hh.

◆ fPrimPDG_ID

G4int RMC01AnalysisManager::fPrimPDG_ID
private

Definition at line 173 of file RMC01AnalysisManager.hh.

◆ fAlpha_or_E0

G4double RMC01AnalysisManager::fAlpha_or_E0
private

Definition at line 174 of file RMC01AnalysisManager.hh.

◆ fAmplitude_prim_spectrum

G4double RMC01AnalysisManager::fAmplitude_prim_spectrum
private

Definition at line 175 of file RMC01AnalysisManager.hh.

◆ fEmin_prim_spectrum

G4double RMC01AnalysisManager::fEmin_prim_spectrum
private

Definition at line 176 of file RMC01AnalysisManager.hh.

◆ fEmax_prim_spectrum

G4double RMC01AnalysisManager::fEmax_prim_spectrum
private

Definition at line 177 of file RMC01AnalysisManager.hh.

◆ fAdjoint_sim_mode

G4bool RMC01AnalysisManager::fAdjoint_sim_mode
private

Definition at line 178 of file RMC01AnalysisManager.hh.

◆ fNb_evt_per_adj_evt

G4int RMC01AnalysisManager::fNb_evt_per_adj_evt
private

Definition at line 179 of file RMC01AnalysisManager.hh.

◆ fTimer

G4Timer* RMC01AnalysisManager::fTimer
private

Definition at line 183 of file RMC01AnalysisManager.hh.

◆ fConvergenceFileOutput

std::fstream RMC01AnalysisManager::fConvergenceFileOutput
private

Definition at line 184 of file RMC01AnalysisManager.hh.

The documentation for this class was generated from the following files:
Generated on Tue Jan 30 2024 12:18:47 for Geant4 examples by   doxygen 1.9.8
Applications | User Support | Publications | Collaboration