Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-13 08:28:45 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file RMC01AnalysisManager.cc
0027 /// \brief Implementation of the RMC01AnalysisManager class
0028 
0029 //      Class Name:        RMC01AnalysisManager
0030 //        Author:               L. Desorgher
0031 //         Organisation:         SpaceIT GmbH
0032 //        Contract:        ESA contract 21435/08/NL/AT
0033 //         Customer:             ESA/ESTEC
0034 //////////////////////////////////////////////////////////////
0035 
0036 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0037 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0038 
0039 #include "RMC01AnalysisManager.hh"
0040 #include "G4AdjointSimManager.hh"
0041 #include "G4SDManager.hh"
0042 #include "RMC01SD.hh"
0043 #include "G4THitsCollection.hh"
0044 #include "G4Electron.hh"
0045 #include "G4Proton.hh"
0046 #include "G4Gamma.hh"
0047 #include "G4Timer.hh"
0048 #include "G4RunManager.hh"
0049 #include "G4MTRunManager.hh"
0050 #include "G4PhysicalConstants.hh"
0051 #include "G4SystemOfUnits.hh"
0052 #include "RMC01AnalysisManagerMessenger.hh"
0053 #include "G4AccumulableManager.hh"
0054 #include "G4Threading.hh"
0055 #include "G4PhysicsOrderedFreeVector.hh"
0056 #include "G4Event.hh"
0057 #include "globals.hh"
0058 #include "G4SPSEneDistribution.hh"
0059 #include "RMC01PrimaryGeneratorAction.hh"
0060 
0061 #include "G4AutoLock.hh"
0062 
0063 namespace
0064 {
0065   G4Mutex fMergeMutex = G4MUTEX_INITIALIZER;
0066   G4double fGlobalEdep = 0.;
0067   G4double fGlobalEdep2 = 0.;
0068   std::fstream fConvergenceFileOutput;
0069   G4Timer fTimer;
0070   G4double elapsed_time=0.;
0071   G4int nb_global_evt_processed=0;
0072 }
0073 
0074 RMC01AnalysisManager *RMC01AnalysisManager::GetInstance()
0075 {
0076   static G4ThreadLocal RMC01AnalysisManager theInstance;
0077   return &theInstance;
0078 }
0079 
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0081 
0082 RMC01AnalysisManager::RMC01AnalysisManager()
0083  :fAccumulated_edep("Accumulated_edep", 0.),
0084   fAccumulated_edep2("Accumulated_edep2", 0.), fMean_edep(0.),
0085   fError_mean_edep(0.), fRelative_error(0.),
0086   fPrecision_to_reach(0.),fStop_run_if_precision_reached(true),
0087   fNb_evt_modulo_for_convergence_test(5000),
0088   fEdep_rmatrix_vs_electron_prim_energy(0),
0089   fElectron_current_rmatrix_vs_electron_prim_energy(0),
0090   fGamma_current_rmatrix_vs_electron_prim_energy(0),
0091   fEdep_rmatrix_vs_gamma_prim_energy(0),
0092   fElectron_current_rmatrix_vs_gamma_prim_energy(0),
0093   fGamma_current_rmatrix_vs_gamma_prim_energy(0),
0094   fEdep_rmatrix_vs_proton_prim_energy(0),
0095   fElectron_current_rmatrix_vs_proton_prim_energy(0),
0096   fProton_current_rmatrix_vs_proton_prim_energy(0),
0097   fGamma_current_rmatrix_vs_proton_prim_energy(0),
0098   fFactoryOn(false),
0099   fPrimSpectrumType(EXPO),
0100   fAlpha_or_E0(.5*MeV),fAmplitude_prim_spectrum (1.),
0101   fEmin_prim_spectrum(1.*keV),fEmax_prim_spectrum (20.*MeV),
0102   fAdjoint_sim_mode(true),fNb_evt_per_adj_evt(2)
0103 { 
0104   
0105   fMsg = new RMC01AnalysisManagerMessenger(this);
0106 
0107   //-------------
0108   
0109 
0110   //---------------------------------
0111   //Primary particle ID for normalisation of adjoint results
0112   //---------------------------------
0113   
0114   fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
0115   
0116   fFileName[0] = "sim";
0117 
0118   G4AccumulableManager *accumulableManager = G4AccumulableManager::Instance();
0119   accumulableManager->Register(fAccumulated_edep);
0120   accumulableManager->Register(fAccumulated_edep2);
0121 }
0122 
0123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0124 
0125 RMC01AnalysisManager::~RMC01AnalysisManager() 
0126 {;
0127 }
0128 
0129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0130 
0131 void RMC01AnalysisManager::BeginOfRun(const G4Run* aRun)
0132 {
0133   fIsEndOfRun =false;
0134   G4AccumulableManager *accumulableManager = G4AccumulableManager::Instance();
0135   accumulableManager->Reset();
0136   fRelative_error=1.;
0137   fMean_edep=0.;
0138   fError_mean_edep=0.;
0139 
0140   fAdjoint_sim_mode =G4AdjointSimManager::GetInstance()->GetAdjointSimMode();
0141 
0142   if (fAdjoint_sim_mode)
0143   {
0144     fNb_evt_per_adj_evt=aRun->GetNumberOfEventToBeProcessed()/
0145                        G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
0146     if (G4Threading::IsMasterThread()
0147       || !G4Threading::IsMultithreadedApplication())
0148     {
0149       fConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt",
0150                                                                 std::ios::out);
0151       fConvergenceFileOutput << "Normalised Edep[MeV]\terror[MeV]"
0152                              << "\tcomputing_time[s]\tnb_adjoint_event"
0153                              << std::endl;
0154       fConvergenceFileOutput.setf(std::ios::scientific);
0155       fConvergenceFileOutput.precision(6);
0156     }
0157   }
0158   else
0159   {
0160     if (G4Threading::IsMasterThread()
0161       || !G4Threading::IsMultithreadedApplication())
0162     {
0163       fConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt",
0164                                                                std::ios::out);
0165       fConvergenceFileOutput
0166         << "Edep per event [MeV]\terror[MeV]\tcomputing_time[s]\tnb_event"
0167         <<std::endl;
0168       fConvergenceFileOutput.setf(std::ios::scientific);
0169       fConvergenceFileOutput.precision(6);
0170     }
0171   }
0172   Book();
0173   if (G4Threading::IsMasterThread()
0174      || !G4Threading::IsMultithreadedApplication())
0175   {
0176     fGlobalEdep = 0.;
0177     fGlobalEdep2 = 0.;
0178     nb_global_evt_processed=0;
0179     fTimer.Start();
0180     elapsed_time=0.;
0181   }
0182 }
0183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0184 
0185 void RMC01AnalysisManager::EndOfRun(const G4Run* aRun)
0186 {
0187   G4int nb_evt=aRun->GetNumberOfEvent();
0188   G4double factor =1./ nb_evt;
0189   G4AccumulableManager *accumulableManager = G4AccumulableManager::Instance();
0190   accumulableManager->Merge();
0191   if (G4Threading::IsMasterThread()
0192     || !G4Threading::IsMultithreadedApplication())
0193   {
0194     fTimer.Stop();
0195     elapsed_time+=fTimer.GetRealElapsed();
0196   }
0197 
0198   fIsEndOfRun =true;
0199   ComputeMeanEdepAndError(fMean_edep,fError_mean_edep,nb_evt);
0200 
0201   if (!fAdjoint_sim_mode)
0202   {
0203     if (G4Threading::IsMasterThread()
0204       || !G4Threading::IsMultithreadedApplication())
0205     {
0206       G4cout << "Results of forward simulation!" << std::endl;
0207       G4cout << "edep per event [MeV] = " << fMean_edep << std::endl;
0208       G4cout << "precision[MeV] = " << fError_mean_edep << std::endl;
0209     }
0210   }
0211   else
0212   {
0213     if (G4Threading::IsMasterThread())
0214     {
0215       G4cout << "Results of reverse/adjoint simulation!" << std::endl;
0216       G4cout << "normalised edep [MeV] = " << fMean_edep << std::endl;
0217       G4cout << "precision[MeV] = " << fError_mean_edep << std::endl;
0218     }
0219     factor = 1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun()
0220            * fNb_evt_per_adj_evt/aRun->GetNumberOfEvent();
0221     nb_evt=nb_evt/fNb_evt_per_adj_evt;
0222   }
0223   Save(factor);
0224   if (G4Threading::IsMasterThread()
0225      || !G4Threading::IsMultithreadedApplication())
0226   {
0227     fConvergenceFileOutput << fMean_edep << '\t' << fError_mean_edep
0228                            << '\t' << elapsed_time << '\t' << nb_evt
0229                            << std::endl;
0230     fConvergenceFileOutput.close();
0231   }
0232 }
0233 
0234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0235 
0236 void RMC01AnalysisManager::BeginOfEvent(const G4Event* )
0237 {
0238 }
0239 
0240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0241 
0242 void RMC01AnalysisManager::EndOfEvent(const G4Event* anEvent)
0243 {  
0244   if (fAdjoint_sim_mode) EndOfEventForAdjointSimulation(anEvent);
0245   else EndOfEventForForwardSimulation(anEvent);
0246 
0247   //Test convergence. The error is already computed
0248   //--------------------------------------
0249   //Nb event is the nb event already performed by the Thread
0250   G4int nb_event = G4RunManager::GetRunManager()
0251                  ->GetCurrentRun()->GetNumberOfEvent();
0252   nb_event = anEvent->GetEventID()+1;
0253 
0254   if (fAdjoint_sim_mode)
0255   {
0256     G4double  n_adj_evt= nb_event/fNb_evt_per_adj_evt;
0257     if (n_adj_evt*fNb_evt_per_adj_evt == nb_event)
0258     {
0259       nb_event =static_cast<G4int>(n_adj_evt);
0260     }        
0261     else nb_event=0;
0262   }
0263 
0264   if (nb_event>100 && fStop_run_if_precision_reached
0265      && fPrecision_to_reach >fRelative_error)
0266   {
0267     G4cout << fPrecision_to_reach*100. << "%  Precision reached!" << std::endl;
0268     G4RunManager::GetRunManager()->AbortRun(true);
0269   }
0270    
0271   if (nb_event>0 && nb_event % fNb_evt_modulo_for_convergence_test == 0)
0272   {
0273     G4AutoLock lock(&fMergeMutex);
0274     fTimer.Stop();
0275     elapsed_time+=fTimer.GetRealElapsed();
0276     fTimer.Start();
0277     fConvergenceFileOutput << fMean_edep << '\t' << fError_mean_edep
0278                            << '\t' << elapsed_time << '\t' << nb_event
0279                            << std::endl;
0280   }
0281 }   
0282 
0283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0284 
0285 void  RMC01AnalysisManager::EndOfEventForForwardSimulation(
0286                                                      const G4Event* anEvent)
0287 {  
0288    
0289    G4SDManager* SDman = G4SDManager::GetSDMpointer();
0290    G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
0291    RMC01DoubleWithWeightHitsCollection* edepCollection =
0292          (RMC01DoubleWithWeightHitsCollection*)
0293               (HCE->GetHC(0));
0294 
0295    RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
0296              (RMC01DoubleWithWeightHitsCollection*)
0297                 (HCE->GetHC(SDman->GetCollectionID("current_electron")));
0298 
0299    RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
0300              (RMC01DoubleWithWeightHitsCollection*)
0301                    (HCE->GetHC(SDman->GetCollectionID("current_proton")));
0302 
0303    RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
0304              (RMC01DoubleWithWeightHitsCollection*)
0305                      (HCE->GetHC(SDman->GetCollectionID("current_gamma")));
0306    
0307    // Total energy deposited in Event
0308    // -------------------------------
0309    G4double totEdep=0; 
0310    size_t i;
0311    if(edepCollection)
0312      for (i=0;i<edepCollection->entries();i++)
0313        totEdep+=(*edepCollection)[i]->GetValue()
0314                *(*edepCollection)[i]->GetWeight();
0315    {
0316      G4AutoLock lock(&fMergeMutex);
0317      nb_global_evt_processed +=1;
0318    }
0319    if (totEdep>0.)
0320    {
0321      fAccumulated_edep +=totEdep ;
0322      fAccumulated_edep2 +=totEdep*totEdep ;
0323 
0324      G4AutoLock lock(&fMergeMutex);
0325      fGlobalEdep += totEdep ;
0326      fGlobalEdep2 += totEdep*totEdep ;
0327      G4PrimaryParticle* thePrimary= anEvent->GetPrimaryVertex()->GetPrimary();
0328      G4double E0= thePrimary->GetG4code()->GetPDGMass();
0329      G4double P=thePrimary->GetMomentum().mag();
0330      G4double prim_ekin =std::sqrt(E0*E0+P*P)-E0;
0331      fEdep_vs_prim_ekin->fill(prim_ekin,totEdep);
0332    } 
0333    ComputeMeanEdepAndError(fMean_edep,fError_mean_edep,nb_global_evt_processed);
0334    if (fError_mean_edep>0) fRelative_error= fError_mean_edep/fMean_edep;
0335                    
0336    // Particle current on sensitive cylinder
0337    // -------------------------------------
0338 
0339    if(electronCurrentCollection)
0340      for (i=0;i<electronCurrentCollection->entries();i++)
0341      {
0342        G4double ekin =(*electronCurrentCollection)[i]->GetValue();
0343        G4double weight=(*electronCurrentCollection)[i]->GetWeight();
0344        fElectron_current->fill(ekin,weight);
0345      }
0346    
0347    if(protonCurrentCollection)
0348      for (i=0;i<protonCurrentCollection->entries();i++)
0349      {
0350        G4double ekin =(*protonCurrentCollection)[i]->GetValue();
0351        G4double weight=(*protonCurrentCollection)[i]->GetWeight();
0352        fProton_current->fill(ekin,weight);
0353      }        
0354    
0355    if(gammaCurrentCollection)
0356      for (i=0;i<gammaCurrentCollection->entries();i++)
0357      {
0358        G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
0359        G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
0360        fGamma_current->fill(ekin,weight);
0361      }
0362 }
0363 
0364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0365 
0366 void  RMC01AnalysisManager::EndOfEventForAdjointSimulation(
0367                                                       const G4Event* anEvent)
0368 {  
0369   //Output from Sensitive volume computed during the forward tracking phase
0370   //-----------------------------------------------------------------------
0371   G4SDManager* SDman = G4SDManager::GetSDMpointer();
0372   G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
0373   RMC01DoubleWithWeightHitsCollection* edepCollection =
0374                   (RMC01DoubleWithWeightHitsCollection*)(
0375                     HCE->GetHC(SDman->GetCollectionID("edep")));
0376 
0377   RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
0378                   (RMC01DoubleWithWeightHitsCollection*)(
0379                  HCE->GetHC(SDman->GetCollectionID("current_electron")));
0380 
0381   RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
0382                   (RMC01DoubleWithWeightHitsCollection*)(
0383                  HCE->GetHC(SDman->GetCollectionID("current_proton")));
0384 
0385   RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
0386                   (RMC01DoubleWithWeightHitsCollection*)(
0387                  HCE->GetHC(SDman->GetCollectionID("current_gamma")));
0388   
0389   // Computation of total energy deposited in fwd tracking phase
0390   // -------------------------------
0391   G4double totEdep=0;
0392   size_t i;
0393   if(edepCollection)
0394     for (i=0;i<edepCollection->entries();i++)
0395       totEdep+=(*edepCollection)[i]->GetValue()
0396               *(*edepCollection)[i]->GetWeight();
0397 
0398   // Output from adjoint tracking phase
0399   //----------------------------------------------------------------------------
0400   
0401   G4AdjointSimManager* theAdjointSimManager =
0402                    G4AdjointSimManager::GetInstance();
0403 
0404   size_t nb_adj_track =
0405        theAdjointSimManager->GetNbOfAdointTracksReachingTheExternalSurface();
0406   G4double total_normalised_weight = 0.;
0407 
0408   // We need to loop over the adjoint tracks that have reached the external
0409   // surface.
0410   for (size_t j=0;j<nb_adj_track;j++)
0411   {
0412     G4int pdg_nb =theAdjointSimManager
0413          ->GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack(j);
0414     G4double prim_ekin=theAdjointSimManager
0415                               ->GetEkinAtEndOfLastAdjointTrack(j);
0416     G4double adj_weight=theAdjointSimManager
0417                              ->GetWeightAtEndOfLastAdjointTrack(j);
0418  
0419     //Factor of normalisation to user defined prim spectrum (power law or exp)
0420     //------------------------------------------------------------------------
0421     G4double normalised_weight = 0.;
0422     if (pdg_nb== fPrimPDG_ID && prim_ekin>= fEmin_prim_spectrum
0423                                          && prim_ekin<= fEmax_prim_spectrum)
0424       normalised_weight =
0425                 adj_weight*PrimDiffAndDirFluxForAdjointSim(prim_ekin);
0426     total_normalised_weight += normalised_weight;
0427 
0428     // Answer matrices
0429     // ---------------
0430     G4AnaH1* edep_rmatrix =0;
0431     G4AnaH2* electron_current_rmatrix =0;
0432     G4AnaH2* gamma_current_rmatrix =0;
0433     G4AnaH2* proton_current_rmatrix =0;
0434 
0435     if (pdg_nb == G4Electron::Electron()->GetPDGEncoding())  //e- matrices
0436     {
0437       edep_rmatrix = fEdep_rmatrix_vs_electron_prim_energy;
0438       electron_current_rmatrix =
0439                           fElectron_current_rmatrix_vs_electron_prim_energy;
0440       gamma_current_rmatrix = fGamma_current_rmatrix_vs_electron_prim_energy;
0441     }
0442     else if (pdg_nb == G4Gamma::Gamma()->GetPDGEncoding())
0443     {
0444       //gammma answer matrices
0445       edep_rmatrix = fEdep_rmatrix_vs_gamma_prim_energy;
0446       electron_current_rmatrix = fElectron_current_rmatrix_vs_gamma_prim_energy;
0447       gamma_current_rmatrix = fGamma_current_rmatrix_vs_gamma_prim_energy;
0448     }
0449     else if (pdg_nb == G4Proton::Proton()->GetPDGEncoding()){
0450       //proton answer matrices
0451       edep_rmatrix = fEdep_rmatrix_vs_proton_prim_energy;
0452       electron_current_rmatrix =
0453                          fElectron_current_rmatrix_vs_proton_prim_energy;
0454       gamma_current_rmatrix = fGamma_current_rmatrix_vs_proton_prim_energy;
0455       proton_current_rmatrix = fProton_current_rmatrix_vs_proton_prim_energy;
0456     }
0457     // Register histo edep vs prim ekin
0458     //----------------------------------
0459     if (normalised_weight>0) fEdep_vs_prim_ekin
0460                         ->fill(prim_ekin,totEdep*normalised_weight);
0461     // Registering answer matrix
0462     //---------------------------
0463     edep_rmatrix->fill(prim_ekin,totEdep*adj_weight/cm2);
0464     
0465     //Registering of current of particles on the sensitive volume
0466     //------------------------------------------------------------
0467    
0468     if(electronCurrentCollection)
0469       for (i=0;i<electronCurrentCollection->entries();i++)
0470       {
0471         G4double ekin =(*electronCurrentCollection)[i]->GetValue();
0472         G4double weight=(*electronCurrentCollection)[i]->GetWeight();
0473         fElectron_current->fill(ekin,weight*normalised_weight);
0474         if(electron_current_rmatrix)
0475           electron_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
0476       }
0477     if(protonCurrentCollection)
0478       for (i=0;i<protonCurrentCollection->entries();i++)
0479       {
0480         G4double ekin =(*protonCurrentCollection)[i]->GetValue();
0481         G4double weight=(*protonCurrentCollection)[i]->GetWeight();
0482         fProton_current->fill(ekin,weight*normalised_weight);
0483         if(proton_current_rmatrix)
0484           proton_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
0485       }
0486     if(gammaCurrentCollection)
0487       for (i=0;i<gammaCurrentCollection->entries();i++)
0488       {
0489         G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
0490         G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
0491         fGamma_current->fill(ekin,weight*normalised_weight);
0492         if(gamma_current_rmatrix)
0493           gamma_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
0494       }
0495   }
0496 
0497   // Registering of total energy deposited in Event
0498   // -------------------------------
0499   {
0500     G4AutoLock lock(&fMergeMutex);
0501     nb_global_evt_processed +=1;
0502   }
0503   G4bool new_mean_computed=false;
0504   if (totEdep>0.)
0505   {
0506     if (total_normalised_weight>0.)
0507     {
0508       G4double edep=totEdep* total_normalised_weight;
0509 
0510       // Check if the edep is not wrongly too high
0511       // -----------------------------------------
0512       G4double new_mean , new_error;
0513       fAccumulated_edep +=edep;
0514       fAccumulated_edep2 +=edep*edep;
0515 
0516 
0517       G4AutoLock lock(&fMergeMutex);
0518       fGlobalEdep += edep;
0519       fGlobalEdep2 += edep*edep;
0520 
0521       ComputeMeanEdepAndError(new_mean,new_error,nb_global_evt_processed);
0522       G4double new_relative_error = 1.;
0523       if ( new_error >0)
0524         new_relative_error = new_error/ new_mean;
0525       fMean_edep = new_mean;
0526       fError_mean_edep = new_error;
0527       fRelative_error =new_relative_error;
0528       new_mean_computed=true;
0529     }
0530     if (!new_mean_computed)
0531     {
0532       G4AutoLock lock(&fMergeMutex);
0533       ComputeMeanEdepAndError(fMean_edep,fError_mean_edep,
0534                               nb_global_evt_processed);
0535       if (fError_mean_edep>0)
0536         fRelative_error= fError_mean_edep/fMean_edep;
0537     }
0538   }
0539 }
0540 
0541 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0542 
0543 G4double RMC01AnalysisManager::
0544 PrimDiffAndDirFluxForAdjointSim(G4double prim_energy)
0545 { 
0546   G4double flux=fAmplitude_prim_spectrum;
0547   if (fPrimSpectrumType ==EXPO)
0548     flux*=std::exp(-prim_energy/fAlpha_or_E0);
0549   else
0550     flux*=std::pow(prim_energy, -fAlpha_or_E0);
0551   return flux;
0552 }
0553 
0554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0555 /*
0556 void  RMC01AnalysisManager::WriteHisto(G4AnaH1* anHisto,
0557             G4double scaling_factor, G4String fileName, G4String header_lines)
0558 { std::fstream FileOutput(fileName, std::ios::out);
0559   FileOutput<<header_lines;
0560   FileOutput.setf(std::ios::scientific);
0561   FileOutput.precision(6);
0562 
0563   for (G4int i =0;i<G4int(anHisto->axis().bins());i++) {
0564         FileOutput<<anHisto->axis().bin_lower_edge(i)
0565               <<'\t'<<anHisto->axis().bin_upper_edge(i)
0566               <<'\t'<<anHisto->bin_height(i)*scaling_factor
0567               <<'\t'<<anHisto->bin_error(i)*scaling_factor<<std::endl;
0568   }
0569 }
0570 
0571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0572 
0573 void  RMC01AnalysisManager::WriteHisto(G4AnaH2* anHisto,
0574             G4double scaling_factor, G4String fileName, G4String header_lines)
0575 { std::fstream FileOutput(fileName, std::ios::out);
0576   FileOutput<<header_lines;
0577   
0578   FileOutput.setf(std::ios::scientific);
0579   FileOutput.precision(6);
0580 
0581   for (G4int i =0;i<G4int(anHisto->axis_x().bins());i++) {
0582     for (G4int j =0;j<G4int(anHisto->axis_y().bins());j++) {
0583        FileOutput<<anHisto->axis_x().bin_lower_edge(i)
0584                      <<'\t'<<anHisto->axis_x().bin_upper_edge(i)
0585                    <<'\t'<<anHisto->axis_y().bin_lower_edge(i)
0586                    <<'\t'<<anHisto->axis_y().bin_upper_edge(i)
0587                    <<'\t'<<anHisto->bin_height(i,j)*scaling_factor
0588                 <<'\t'<<anHisto->bin_error(i,j)*scaling_factor
0589                                                                   <<std::endl;
0590         }
0591   }
0592 }
0593 */
0594 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0595 
0596 void RMC01AnalysisManager::
0597 ComputeMeanEdepAndError(G4double& mean,G4double& error,
0598                         G4int nb_of_global_evt_processed)
0599 {  
0600   G4double nb_event=G4double(nb_of_global_evt_processed);
0601   G4double factor=1.;
0602   G4double factor1=1.;
0603 
0604   if (fAdjoint_sim_mode)
0605   {
0606      nb_event /= fNb_evt_per_adj_evt;
0607      G4double nb_evt_sim = G4double(G4AdjointSimManager::GetInstance()
0608                          ->GetNbEvtOfLastRun());
0609      factor = nb_evt_sim;
0610      factor1 = factor*nb_evt_sim;
0611   }
0612 
0613   // VI: error computation now is based on number of entries and not 
0614   //     number of events
0615   // LD: This is wrong! With the use of fNentry the results were no longer
0616   //     correctly normalised. The mean and the error should be computed
0617   //     with nb_event. The old computation has been reset.
0618 
0619   if (nb_event>1.)
0620   {
0621     mean = fGlobalEdep*factor/nb_event;
0622     G4double mean_x2 = fGlobalEdep2*factor1/nb_event;
0623     if (fIsEndOfRun)
0624     {
0625       mean = fAccumulated_edep.GetValue()*factor/nb_event;
0626       mean_x2 = fAccumulated_edep2.GetValue()*factor1/nb_event;
0627     }
0628     G4double n_eff= nb_event;
0629     G4double var_weighted = (mean_x2 - mean*mean)*n_eff/(n_eff-1.);
0630     error = std::sqrt(std::max(var_weighted, 0.)/n_eff);
0631   }
0632   else
0633   {
0634     mean=0;
0635     error=0;
0636   }
0637 }
0638 
0639 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0640 
0641 void RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim(
0642                      const G4String& particle_name, G4double omni_fluence,
0643                      G4double E0, G4double Emin, G4double Emax)
0644 {
0645   fPrimSpectrumType = EXPO;
0646   if (particle_name == "e-" )
0647     fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
0648   else if (particle_name == "gamma")
0649     fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEncoding();
0650   else if (particle_name == "proton")
0651     fPrimPDG_ID = G4Proton::Proton()->GetPDGEncoding();
0652   else
0653   {
0654     G4cout << "The particle that you did select is not in the candidate "
0655            << "list for primary [e-, gamma, proton]!" << G4endl;
0656     return;
0657   }        
0658   fAlpha_or_E0 = E0 ;
0659   fAmplitude_prim_spectrum = omni_fluence/E0
0660                            / (std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
0661   fEmin_prim_spectrum = Emin ;
0662   fEmax_prim_spectrum = Emax;
0663 }
0664 
0665 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0666 
0667 void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim(
0668                            const G4String& particle_name, G4double omni_fluence,
0669                                  G4double alpha, G4double Emin,G4double Emax)
0670 {
0671   fPrimSpectrumType=POWER;
0672   if (particle_name == "e-" )
0673     fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
0674   else if (particle_name == "gamma")
0675     fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEncoding();
0676   else if (particle_name == "proton")
0677     fPrimPDG_ID = G4Proton::Proton()->GetPDGEncoding();
0678   else
0679   {
0680     G4cout << "The particle that you did select is not in the candidate"
0681            << " list for primary [e-, gamma, proton]!"<<G4endl;
0682     return;
0683   }        
0684 
0685   if (alpha ==1.)
0686   {
0687     fAmplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
0688   }
0689   else
0690   {
0691     G4double p=1.-alpha;
0692     fAmplitude_prim_spectrum = omni_fluence/p
0693                              / (std::pow(Emax,p)-std::pow(Emin,p))/4./pi;
0694   }
0695 
0696   fAlpha_or_E0 = alpha;
0697   fEmin_prim_spectrum = Emin;
0698   fEmax_prim_spectrum = Emax;
0699 }
0700 
0701 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0702 
0703 void RMC01AnalysisManager::Book()
0704 {
0705   //----------------------
0706   //Creation of histograms
0707   //----------------------
0708 
0709   //Energy binning of the histograms : 60 log bins over [1keV-1GeV]
0710 
0711   G4double emin=1.*keV;
0712   G4double emax=1.*GeV;
0713 
0714   //file_name
0715   fFileName[0]="forward_sim";
0716   if (fAdjoint_sim_mode) fFileName[0]="adjoint_sim";
0717 
0718   //Histo manager
0719   G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
0720   theHistoManager->SetDefaultFileType("root");
0721   G4String extension = theHistoManager->GetFileType();
0722   fFileName[1] = fFileName[0] + "." + extension;
0723   theHistoManager->SetFirstHistoId(1);
0724 
0725   G4bool fileOpen = theHistoManager->OpenFile(fFileName[0]);
0726   if (!fileOpen)
0727   {
0728     G4cout << "\n---> RMC01AnalysisManager::Book(): cannot open "
0729            << fFileName[1]
0730            << G4endl;
0731     return;
0732   }
0733 
0734   // Create directories
0735   //   theHistoManager->SetHistoDirectoryName("histo");
0736 
0737   // Histograms for :
0738   //        1)the forward simulation results
0739   //        2)the Reverse MC  simulation results normalised to a user spectrum
0740   //------------------------------------------------------------------------
0741 
0742   G4int idHisto =
0743       theHistoManager->CreateH1(G4String("Edep_vs_prim_ekin"),
0744       G4String("edep vs e- primary energy"),60,emin,emax,
0745       "none","none",G4String("log"));
0746   fEdep_vs_prim_ekin = theHistoManager->GetH1(idHisto);
0747 
0748   idHisto = theHistoManager->CreateH1(G4String("elecron_current"),
0749         G4String("electron"),60,emin,emax,
0750         "none","none",G4String("log"));
0751 
0752   fElectron_current  =  theHistoManager->GetH1(idHisto);
0753 
0754   idHisto= theHistoManager->CreateH1(G4String("proton_current"),
0755         G4String("proton"),60,emin,emax,
0756         "none","none",G4String("log"));
0757   fProton_current=theHistoManager->GetH1(idHisto);
0758 
0759   idHisto= theHistoManager->CreateH1(G4String("gamma_current"),
0760           G4String("gamma"),60,emin,emax,
0761           "none","none",G4String("log"));
0762   fGamma_current=theHistoManager->GetH1(idHisto);
0763 
0764   // Response matrices for the adjoint simulation only
0765   //-----------------------------------------------
0766   if (fAdjoint_sim_mode)
0767   {
0768   // Response matrices for external isotropic e- source
0769   //--------------------------------------------------
0770 
0771     idHisto =
0772       theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_electron_prim_energy"),
0773         G4String("electron RM vs e- primary energy"),60,emin,emax,
0774         "none","none",G4String("log"));
0775      fEdep_rmatrix_vs_electron_prim_energy = theHistoManager->GetH1(idHisto);
0776 
0777     idHisto =
0778       theHistoManager->
0779          CreateH2(G4String("Electron_current_rmatrix_vs_electron_prim_energy"),
0780          G4String("electron current  RM vs e- primary energy"),
0781          60,emin,emax,60,emin,emax,
0782          "none","none","none","none",G4String("log"),G4String("log"));
0783 
0784     fElectron_current_rmatrix_vs_electron_prim_energy =
0785                                              theHistoManager->GetH2(idHisto);
0786 
0787     idHisto =
0788       theHistoManager->
0789            CreateH2(G4String("Gamma_current_rmatrix_vs_electron_prim_energy"),
0790            G4String("gamma current  RM vs e- primary energy"),
0791            60,emin,emax,60,emin,emax,
0792            "none","none","none","none",G4String("log"),G4String("log"));
0793 
0794     fGamma_current_rmatrix_vs_electron_prim_energy =
0795                                            theHistoManager->GetH2(idHisto);
0796 
0797     //Response matrices for external isotropic gamma source
0798 
0799     idHisto =
0800       theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_gamma_prim_energy"),
0801          G4String("electron RM vs gamma primary energy"),60,emin,emax,
0802         "none","none",G4String("log"));
0803     fEdep_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH1(idHisto);
0804 
0805     idHisto =
0806       theHistoManager->
0807           CreateH2(G4String("Electron_current_rmatrix_vs_gamma_prim_energy"),
0808           G4String("electron current  RM vs gamma primary energy"),
0809           60,emin,emax,60,emin,emax,
0810           "none","none","none","none",G4String("log"),G4String("log"));
0811 
0812     fElectron_current_rmatrix_vs_gamma_prim_energy =
0813                                                theHistoManager->GetH2(idHisto);
0814 
0815     idHisto =
0816       theHistoManager->
0817           CreateH2(G4String("Gamma_current_rmatrix_vs_gamma_prim_energy"),
0818           G4String("gamma current  RM vs gamma primary energy"),
0819           60,emin,emax,60,emin,emax,
0820           "none","none","none","none",G4String("log"),G4String("log"));
0821 
0822     fGamma_current_rmatrix_vs_gamma_prim_energy= theHistoManager->GetH2(idHisto);
0823 
0824     // Response matrices for external isotropic proton source
0825 
0826     idHisto =
0827       theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_proton_prim_energy"),
0828          G4String("electron RM vs proton primary energy"),60,emin,emax,
0829          "none","none",G4String("log"));
0830     fEdep_rmatrix_vs_proton_prim_energy = theHistoManager->GetH1(idHisto);
0831 
0832     idHisto =
0833       theHistoManager->
0834          CreateH2(G4String("Electron_current_rmatrix_vs_proton_prim_energy"),
0835          G4String("electron current  RM vs proton primary energy"),
0836          60,emin,emax,60,emin,emax,
0837          "none","none","none","none",G4String("log"),G4String("log"));
0838 
0839     fElectron_current_rmatrix_vs_proton_prim_energy =
0840                                                 theHistoManager->GetH2(idHisto);
0841 
0842     idHisto =
0843       theHistoManager->
0844          CreateH2(G4String("Gamma_current_rmatrix_vs_proton_prim_energy"),
0845          G4String("gamma current  RM vs proton primary energy"),
0846          60,emin,emax,60,emin,emax,
0847          "none","none","none","none",G4String("log"),G4String("log"));
0848 
0849     fGamma_current_rmatrix_vs_proton_prim_energy =
0850                                               theHistoManager->GetH2(idHisto);
0851 
0852     idHisto =
0853       theHistoManager->
0854          CreateH2(G4String("Proton_current_rmatrix_vs_proton_prim_energy"),
0855          G4String("proton current  RM vs proton primary energy"),
0856          60,emin,emax,60,emin,emax,
0857          "none","none","none","none",G4String("log"),G4String("log"));
0858 
0859     fProton_current_rmatrix_vs_proton_prim_energy =
0860                                                 theHistoManager->GetH2(idHisto);
0861   }
0862   fFactoryOn = true;
0863   G4cout << "\n----> Histogram Tree is opened in " << fFileName[1] << G4endl;
0864 }
0865 
0866 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0867 
0868 void RMC01AnalysisManager::Save(G4double scaling_factor)
0869 {
0870   if (fFactoryOn)
0871   {
0872     G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
0873 
0874     //scaling of results
0875     //-----------------
0876     for (int ind=1; ind<=theHistoManager->GetNofH1s();ind++)
0877     {
0878        theHistoManager->SetH1Ascii(ind,true);
0879        theHistoManager->ScaleH1(ind,scaling_factor);
0880     }
0881     for (int ind=1; ind<=theHistoManager->GetNofH2s();ind++)
0882        theHistoManager->ScaleH2(ind,scaling_factor);
0883 
0884     theHistoManager->Write();
0885     theHistoManager->CloseFile();
0886     G4cout << "\n----> Histogram Tree is saved in " << fFileName[1] << G4endl;
0887 
0888     fFactoryOn = false;
0889   }
0890 }
0891 
0892 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0893 
0894 void RMC01AnalysisManager::
0895 SetUserDefinedSpectrumPointForAdjointSim(const  G4String& particle_name,
0896                                          G4double omni_fluence,
0897                                          G4bool is_arbitrary_point_wise)
0898 {
0899   fPrimSpectrumType = USER;
0900   if (particle_name == "e-" )
0901     fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
0902   else if (particle_name == "gamma")
0903     fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEncoding();
0904   else if (particle_name == "proton")
0905     fPrimPDG_ID = G4Proton::Proton()->GetPDGEncoding();
0906   else
0907   {
0908     G4cout << "The particle that you did select is not in the candidate "
0909            << "list for primary [e-, gamma, proton]!" << G4endl;
0910     return;
0911   }
0912   G4cout << "Define user spectrum" << G4endl;
0913   G4cout << dynamic_cast<const RMC01PrimaryGeneratorAction*>
0914             (G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction())
0915          << G4endl;
0916 
0917   G4cout << G4MTRunManager::GetMasterRunManager()
0918             ->GetUserPrimaryGeneratorAction() << G4endl;
0919   G4SPSEneDistribution* theSPSEneDistribution =
0920     dynamic_cast<const RMC01PrimaryGeneratorAction*>
0921     (G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction())
0922      ->GetParticleGun()->GetCurrentSource()->GetEneDist();
0923   G4PhysicsOrderedFreeVector* aVec = nullptr;
0924 
0925   f_energy_vec.clear();
0926   f_flux_vec.clear();
0927 
0928   if (!is_arbitrary_point_wise)
0929   {
0930     aVec = new G4PhysicsOrderedFreeVector(theSPSEneDistribution
0931                                           ->GetUserDefinedEnergyHisto());
0932     fEmin_prim_spectrum = aVec->GetMinEnergy();
0933     fEmax_prim_spectrum = aVec->GetMaxEnergy();
0934     for (size_t i=0;i<aVec->GetVectorLength();i++)
0935     {
0936       f_energy_vec.push_back(aVec->GetLowEdgeEnergy(i));
0937       f_flux_vec.push_back((*aVec)[i]);
0938     }
0939     f_flux_vec[0]=0.;
0940   }
0941   else
0942   {
0943     G4PhysicsOrderedFreeVector aVec1= theSPSEneDistribution->GetArbEnergyHisto();
0944 
0945     // Interpolation spline method does not  works properly
0946     // Therfore we use the linear one. We will use G4DataInterpolation
0947 
0948     size_t n = aVec1.GetVectorLength();
0949     G4double* evec=new G4double[n-1];
0950     G4double* fvec=new G4double[n-1];
0951     for (size_t i=0;i<n-1;i++)
0952     {
0953       evec[i]=aVec1.GetLowEdgeEnergy(i);
0954       fvec[i]=aVec1[i];
0955     }
0956 
0957     f_flux_vec.push_back(0);
0958     G4double  flux1=aVec1[0];
0959     G4double E1 =aVec1.GetLowEdgeEnergy(0);
0960     f_energy_vec.push_back(E1);
0961     fEmin_prim_spectrum = E1;
0962 
0963     for (size_t i=1;i<n-1;i++)
0964     {
0965       G4double flux2=aVec1[i];
0966       G4double E2 = aVec1.GetLowEdgeEnergy(i);
0967       size_t lastidx=0;
0968 
0969       if ((E2/E1)>1.05)  //rebin this interval
0970       {
0971         G4double E =E1*1.05;
0972         G4double flux=0.;
0973         G4double logE1,logE2,logf1,logf2,d;
0974         logE1=std::log(E1);
0975         logE2=std::log(E2);
0976         logf1=std::log(flux1);
0977         logf2=std::log(flux2);
0978         d=(logf2-logf1)/(logE2-logE1);
0979 
0980         while(E<E2)
0981         {
0982           flux=aVec1.Value(E,lastidx);
0983           f_energy_vec.push_back(E);
0984           flux=std::exp(logf1+d*(std::log(E)-logE1));
0985           f_flux_vec.push_back((E-E1)*(flux1+flux)/2.);
0986           E1=E;
0987           flux1=flux;
0988           E=E1*1.05;
0989         }
0990       }
0991 
0992       f_energy_vec.push_back(E2);
0993       f_flux_vec.push_back((E2-E1)*(flux1+flux2)/2.);
0994       flux1=flux2;
0995       E1=E2;
0996       fEmax_prim_spectrum = E2;
0997     }
0998   }
0999 
1000   G4double integral_flux =0.;
1001   for (size_t i=0;i<f_flux_vec.size();i++)
1002     integral_flux += f_flux_vec[i];
1003   for (size_t i=1;i<f_flux_vec.size();i++)
1004   {
1005     G4double dE=f_energy_vec[i]-f_energy_vec[i-1];
1006     f_flux_vec[i] = omni_fluence*f_flux_vec[i]/integral_flux/4/pi/dE;
1007   }
1008 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!