include/Geant4/G4AdjointCSManager.hh

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 ////////////////////////////////////////////////////////////////////////////////
0027 //  Class:    G4AdjointCSManager
0028 //  Author:         L. Desorgher
0029 //  Organisation:   SpaceIT GmbH
0030 //
0031 // Class is responsible for the management of all adjoint cross section
0032 // matrices, and for the computation of the total forward and adjoint cross
0033 // sections. Total adjoint and forward cross sections are needed to correct the
0034 // weight of a particle after a tracking step or after the occurrence of a
0035 // reverse reaction. It is also used to sample an adjoint secondary from a
0036 // given adjoint cross section matrix.
0037 //
0038 ////////////////////////////////////////////////////////////////////////////////
0039
0040 #ifndef G4AdjointCSManager_h
0041 #define G4AdjointCSManager_h 1
0042
0043 #include "globals.hh"
0044 #include "G4AdjointCSMatrix.hh"
0045 #include "G4ThreadLocalSingleton.hh"
0046
0047 #include <vector>
0048
0049 class G4Element;
0050 class G4Material;
0051 class G4MaterialCutsCouple;
0052 class G4ParticleDefinition;
0053 class G4PhysicsTable;
0054 class G4VEmProcess;
0055 class G4VEmAdjointModel;
0056 class G4VEnergyLossProcess;
0057
0058 class G4AdjointCSManager
0059 {
0060   friend class G4ThreadLocalSingleton<G4AdjointCSManager>;
0061
0062  public:
0063   ~G4AdjointCSManager();
0064   static G4AdjointCSManager* GetAdjointCSManager();
0065
0066   G4int GetNbProcesses();
0067
0068   // Registration of the different models and processes
0069
0070   std::size_t RegisterEmAdjointModel(G4VEmAdjointModel*);
0071
0072   void RegisterEmProcess(G4VEmProcess* aProcess,
0073                          G4ParticleDefinition* aPartDef);
0074
0075   void RegisterEnergyLossProcess(G4VEnergyLossProcess* aProcess,
0076                                  G4ParticleDefinition* aPartDef);
0077
0078   void RegisterAdjointParticle(G4ParticleDefinition* aPartDef);
0079
0080   // Building of the CS Matrices and Total Forward and Adjoint LambdaTables
0081   void BuildCrossSectionMatrices();
0082
0083   void BuildTotalSigmaTables();
0084
0085   // Get TotalCrossSections form Total Lambda Tables, Needed for Weight
0086   // correction and scaling of the
0087   G4double GetTotalAdjointCS(G4ParticleDefinition* aPartDef, G4double Ekin,
0088                              const G4MaterialCutsCouple* aCouple);
0089
0090   G4double GetTotalForwardCS(G4ParticleDefinition* aPartDef, G4double Ekin,
0091                              const G4MaterialCutsCouple* aCouple);
0092
0093   G4double GetAdjointSigma(G4double Ekin_nuc, std::size_t index_model,
0094                            G4bool is_scat_proj_to_proj,
0095                            const G4MaterialCutsCouple* aCouple);
0096
0097   void GetEminForTotalCS(G4ParticleDefinition* aPartDef,
0098                          const G4MaterialCutsCouple* aCouple,
0099                          G4double& emin_adj, G4double& emin_fwd);
0100
0101   void GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
0102                         const G4MaterialCutsCouple* aCouple,
0103                         G4double& e_sigma_max, G4double& sigma_max);
0104
0105   void GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
0106                         const G4MaterialCutsCouple* aCouple,
0107                         G4double& e_sigma_max, G4double& sigma_max);
0108
0109   // CrossSection Correction 1 or FwdCS/AdjCS following the G4boolean value of
0110   // forward_CS_is_used and forward_CS_mode
0111   G4double GetCrossSectionCorrection(G4ParticleDefinition* aPartDef,
0112                                      G4double PreStepEkin,
0113                                      const G4MaterialCutsCouple* aCouple,
0114                                      G4bool& fwd_is_used);
0115
0116   // Cross section mode
0117   inline void SetFwdCrossSectionMode(G4bool aBool) { fForwardCSMode = aBool; }
0118
0119   // Weight correction
0120   G4double GetContinuousWeightCorrection(G4ParticleDefinition* aPartDef,
0121                                          G4double PreStepEkin,
0122                                          G4double AfterStepEkin,
0123                                          const G4MaterialCutsCouple* aCouple,
0124                                          G4double step_length);
0125
0126   G4double GetPostStepWeightCorrection();
0127
0128   // called by the adjoint model to get the CS, if not otherwise specified
0129   G4double ComputeAdjointCS(G4Material* aMaterial, G4VEmAdjointModel* aModel,
0130                             G4double PrimEnergy, G4double Tcut,
0131                             G4bool isScatProjToProj,
0132                             std::vector<G4double>& AdjointCS_for_each_element);
0133
0134   // called by the adjoint model to sample secondary energy from the CS matrix
0135   G4Element* SampleElementFromCSMatrices(G4Material* aMaterial,
0136                                          G4VEmAdjointModel* aModel,
0137                                          G4double PrimEnergy, G4double Tcut,
0138                                          G4bool isScatProjToProj);
0139
0140   // Total Adjoint CS is computed at initialisation phase
0141   G4double ComputeTotalAdjointCS(const G4MaterialCutsCouple* aMatCutCouple,
0142                                  G4ParticleDefinition* aPart,
0143                                  G4double PrimEnergy);
0144
0145   G4ParticleDefinition* GetAdjointParticleEquivalent(
0146     G4ParticleDefinition* theFwdPartDef);
0147
0148   G4ParticleDefinition* GetForwardParticleEquivalent(
0149     G4ParticleDefinition* theAdjPartDef);
0150
0151   // inline
0152   inline void SetIon(G4ParticleDefinition* adjIon, G4ParticleDefinition* fwdIon)
0153   {
0154     fAdjIon = adjIon;
0155     fFwdIon = fwdIon;
0156   }
0157
0158  private:
0159   G4AdjointCSManager();
0160
0161   void DefineCurrentMaterial(const G4MaterialCutsCouple* couple);
0162
0163   void DefineCurrentParticle(const G4ParticleDefinition* aPartDef);
0164
0165   G4double ComputeAdjointCS(G4double aPrimEnergy,
0166                             G4AdjointCSMatrix* anAdjointCSMatrix,
0167                             G4double Tcut);
0168
0169   std::vector<G4AdjointCSMatrix*> BuildCrossSectionsModelAndElement(
0170     G4VEmAdjointModel* aModel, G4int Z, G4int A, G4int nbin_pro_decade);
0171
0172   std::vector<G4AdjointCSMatrix*> BuildCrossSectionsModelAndMaterial(
0173     G4VEmAdjointModel* aModel, G4Material* aMaterial, G4int nbin_pro_decade);
0174
0175   static constexpr G4double fTmin = 0.1 * CLHEP::keV;
0176   static constexpr G4double fTmax = 100. * CLHEP::TeV;
0177   // fNbins chosen to avoid error
0178   // in the CS value close to CS jump. (For example at Tcut)
0179   static constexpr G4int fNbins = 320;
0180
0181   static G4ThreadLocal G4AdjointCSManager* fInstance;
0182
0183   // only one ion can be considered by simulation
0184   G4ParticleDefinition* fAdjIon = nullptr;
0185   G4ParticleDefinition* fFwdIon = nullptr;
0186
0187   G4MaterialCutsCouple* fCurrentCouple = nullptr;
0188   G4Material* fCurrentMaterial         = nullptr;
0189
0190   // x dim is for G4VAdjointEM*, y dim is for elements
0191   std::vector<std::vector<G4AdjointCSMatrix*>>
0192     fAdjointCSMatricesForScatProjToProj;
0193
0194   std::vector<std::vector<G4AdjointCSMatrix*>> fAdjointCSMatricesForProdToProj;
0195
0196   std::vector<G4VEmAdjointModel*> fAdjointModels;
0197
0198   std::vector<std::size_t> fIndexOfAdjointEMModelInAction;
0199   std::vector<G4bool> fIsScatProjToProj;
0200   std::vector<std::vector<G4double>> fLastAdjointCSVsModelsAndElements;
0201
0202   // total adjoint and total forward cross section table in function of material
0203   // and in function of adjoint particle type
0204   std::vector<G4PhysicsTable*> fTotalFwdSigmaTable;
0205   std::vector<G4PhysicsTable*> fTotalAdjSigmaTable;
0206
0207   // Sigma table for each G4VAdjointEMModel
0208   std::vector<G4PhysicsTable*> fSigmaTableForAdjointModelScatProjToProj;
0209   std::vector<G4PhysicsTable*> fSigmaTableForAdjointModelProdToProj;
0210
0211   std::vector<std::vector<G4double>> fEminForFwdSigmaTables;
0212   std::vector<std::vector<G4double>> fEminForAdjSigmaTables;
0213   std::vector<std::vector<G4double>> fEkinofFwdSigmaMax;
0214   std::vector<std::vector<G4double>> fEkinofAdjSigmaMax;
0215
0216   // list of forward G4VEmProcess and of G4VEnergyLossProcess for the different
0217   // adjoint particle
0218   std::vector<std::vector<G4VEmProcess*>*> fForwardProcesses;
0219   std::vector<std::vector<G4VEnergyLossProcess*>*> fForwardLossProcesses;
0220
0221   // list of adjoint particles considered
0222   std::vector<G4ParticleDefinition*> fAdjointParticlesInAction;
0223
0224   G4double fMassRatio              = 1.;  // ion
0225   G4double fLastCSCorrectionFactor = 1.;
0226
0227   std::size_t fCurrentParticleIndex = 0;
0228   std::size_t fCurrentMatIndex      = 0;
0229
0230   G4bool fCSMatricesBuilt = false;
0231   G4bool fSigmaTableBuilt = false;
0232   G4bool fForwardCSUsed   = true;
0233   G4bool fForwardCSMode   = true;
0234   // Two CS mode are possible:
0235   // 1) fForwardCSMode = false, the Adjoint CS are used as it is implying
0236   //    an AlongStep Weight Correction.
0237   // 2) fForwardCSMode = true, the Adjoint CS are scaled to have the total
0238   //    adjoint CS equal to the fwd one implying a PostStep Weight Correction.
0239   // For energies where the total Fwd CS or the total adjoint CS are zero,
0240   // the scaling is not possible and fForwardCSUsed is set to false
0241 };
0242 #endif