include/Geant4/G4LowPAIH2O.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 // G4LowPAIxs.hh -- header file
0027 //
0028 // GEANT 4 class header file --- Copyright CERN 1995
0029 // CERB Geneva Switzerland
0030 //
0031 // for information related to this code, please, contact
0032 // CERN, CN Division, ASD Group
0033 //
0034 // Preparation of ionizing collision cross section according to Photo Absorption
0035 // Ionization (PAI) model for simulation of ionization energy losses in very thin
0036 // layers of water for low energy protons and electrons.
0037 //
0038 // Author: Vladimir.Grichine@cern.ch
0039 //
0040 // History:
0041 //
0042
0043 // 2.12.24, V. Grichine: 1st version
0044
0045 #ifndef G4LOWPAIH2O_HH
0046 #define G4LOWPAIH2O_HH
0047
0048 #include "G4ios.hh"
0049 #include "globals.hh"
0050 #include "Randomize.hh"
0051 #include "G4PhysicsLogVector.hh"
0052 #include "G4DataVector.hh"
0053 #include "G4PhysicsTable.hh"
0054 #include "G4VEmModel.hh"
0055 #include "G4VEmFluctuationModel.hh"
0056
0057 class G4MaterialCutsCouple;
0058 class G4Material;
0059 class G4SandiaTable;
0060 class G4PhysicsLogVector;
0061 class G4PhysicsTable;
0062 class G4ParticleDefinition;
0063 class G4ParticleChangeForLoss;
0064
0065 class G4LowPAIH2O :
0066    public G4VEmModel  , public G4VEmFluctuationModel
0067 {
0068 public:
0069   // Constructors
0070
0071   explicit G4LowPAIH2O( const G4ParticleDefinition* p = nullptr,
0072            const G4String& nam = "lowpaih2o");
0073
0074  ~G4LowPAIH2O() override;
0075
0076   G4LowPAIH2O & operator=(const G4LowPAIH2O &right) = delete;
0077   G4LowPAIH2O(const G4LowPAIH2O&) = delete;
0078
0079   // methods
0080
0081   void Initialise(const G4ParticleDefinition*, const G4DataVector&) override;
0082
0083   void InitialiseLocal(const G4ParticleDefinition*,
0084                        G4VEmModel* masterModel) override;
0085
0086   G4double CrossSectionPerVolume(const G4Material*,
0087                                  const G4ParticleDefinition*,
0088                                  G4double kineticEnergy,
0089                                  G4double cutEnergy,
0090                                  G4double maxEnergy) override;
0091
0092   G4double CrossSectionPerAtom(
0093                                  const G4ParticleDefinition*,
0094                                  G4double kineticEnergy, G4double Z,
0095                                       G4double A,
0096                                  G4double cutEnergy,
0097                  G4double maxEnergy);
0098
0099   virtual G4double ComputeCrossSectionPerElectron(
0100                                  const G4ParticleDefinition*,
0101                                  G4double kineticEnergy,
0102                                  G4double cutEnergy,
0103                                  G4double maxEnergy);
0104
0105   void SampleSecondaries(std::vector<G4DynamicParticle*>*,
0106                                  const G4MaterialCutsCouple*,
0107                                  const G4DynamicParticle*,
0108                                  G4double tmin,
0109                      G4double maxEnergy) override;
0110
0111   G4double SampleFluctuations(const G4MaterialCutsCouple*,
0112                                       const G4DynamicParticle*,
0113                                       const G4double tcut,
0114                                       const G4double tmax,
0115                                       const G4double length,
0116                                       const G4double meanLoss) override;
0117   void CorrectionsAlongStep(const G4Material*,
0118                 const G4ParticleDefinition*,
0119                 const G4double kinEnergy,
0120                 const G4double cutEnergy,
0121                 const G4double& length,
0122                 G4double& eloss) override;
0123
0124   G4double Dispersion(const G4Material*, const G4DynamicParticle*,
0125            const G4double, const G4double, const G4double) override {return 0.;};
0126
0127   inline G4double ComputeMeanFreePath( const G4ParticleDefinition*,
0128                                      G4double kineticEnergy,
0129                                      const G4Material*,
0130                                      G4double cutEnergy = 0.0,
0131                                      G4double maxEnergy = DBL_MAX);
0132   void Initialize();
0133   void InitRuthELF();
0134
0135   void BuildPhysicsTable(const G4ParticleDefinition* pd);
0136   void BuildPrEnergyTable();
0137   void BuildElEnergyTable();
0138   G4double GetPrTransfer( G4double Tkin);
0139   G4double CorrectPrTransfer( G4double Tkin);
0140   G4double CorrectElTransfer( G4double Tkin);
0141
0142   G4double GetElTransfer( G4double Tkin);
0143   G4double GetPrMFP( G4double Tkin);
0144   G4double GetPrdNdx( G4double Tkin);
0145   G4double GetElMFP( G4double Tkin);
0146   G4double GetEldNdx( G4double Tkin);
0147   G4double PrPAId2Ndxdw( G4double omega );
0148   G4double ElPAId2Ndxdw( G4double omega );
0149
0150   void     SetBe2( G4double be2 ){ fBe2 = be2; };
0151   G4double GetBe2(){ return fBe2; };
0152   void     SetOmega( G4double ww ){ fOmega = ww; };
0153   G4double GetOmega(){ return fOmega;};
0154   void     SetBias( G4double bb ){ fBias = bb; };
0155   G4double GetBias(){ return fBias;};
0156
0157   G4double GetElectronTmax( G4double Tkin );
0158   G4double GetProtonTmax( G4double Tkin );
0159
0160   inline G4double GetSumELF( G4double energy );
0161   inline G4double GetSumRuth( G4double energy );
0162
0163 private:
0164
0165   // Local class members
0166   G4int fTotBin{0};
0167   G4int fBinTr{0};
0168   G4int fBias{0};
0169   G4double fCof{0.0};
0170   G4double fBeta{0.0};
0171   G4double fBe2{0.0};
0172   G4double fTkin{0.0};
0173
0174   G4double fBmin{0.0};
0175   G4double fBmax{0.0};
0176   G4double fWmin{0.0};
0177   G4double fWmax{0.0};
0178   G4double fOmega{0.0};
0179   G4double fElectronDensity{0.0};
0180   G4double fNat{0.0};
0181   G4double fNel{0.0};
0182   G4double fMass{0.0};
0183
0184   static const G4int theBin;
0185   static const G4double theEsum[1007], theELFsum[1007], theRuthSum[1007];
0186
0187   G4DataVector fEsum;
0188   G4DataVector fELFsum, fRuthSum;
0189   G4DataVector fPrWmaxVector;
0190
0191   G4Material* fMat{nullptr};
0192
0193   G4PhysicsLogVector*  fBetaVector{nullptr};
0194   G4PhysicsTable*      fPrEnergyTable{nullptr};
0195   G4PhysicsTable*      fElEnergyTable{nullptr};
0196   G4PhysicsLogVector*  fTransferVector{nullptr};
0197
0198   G4ParticleDefinition* theElectron{nullptr};
0199   G4ParticleDefinition* theProton{nullptr};
0200   G4ParticleChangeForLoss* fParticleChange{nullptr};
0201 };
0202
0203 ///////////////////////////////////////////////////////////////////
0204 ////////////////  Inline methods //////////////////////////////////
0205 ////////////////////////////////////////////////////////////////////
0206
0207
0208 /////////////////// fast (STL) get ELF /////////////////
0209
0210 G4double G4LowPAIH2O::GetSumELF( G4double energy )
0211 {
0212   G4double ee =  energy; // /CLHEP::eV;
0213   G4double elf(0.), y1(0.), y2(0.), x1(0.), x2(0.), aa(0.);
0214
0215   std::size_t nlow = std::lower_bound( fEsum.begin(), fEsum.end(), ee ) - fEsum.begin();
0216
0217   x1 = fEsum[nlow-1];
0218   x2 = fEsum[nlow];
0219   y1 = fELFsum[nlow-1];
0220   y2 = fELFsum[nlow];
0221   aa = (y2-y1)/(x2-x1);
0222   elf = y1 + aa*(ee-x1);
0223   return elf;
0224 }
0225
0226 /////////////////////// fast (STL) get Rutherford /////////////////////
0227
0228 G4double G4LowPAIH2O::GetSumRuth( G4double energy )
0229 {
0230   G4double ee =  energy; // /CLHEP::eV;
0231   G4double ruth(0.), y1(0.), y2(0.), x1(0.), x2(0.), aa(0.);
0232
0233   std::size_t nlow = std::lower_bound( fEsum.begin(), fEsum.end(), ee ) - fEsum.begin();
0234
0235   x1 = fEsum[nlow-1];
0236   x2 = fEsum[nlow];
0237   y1 = fRuthSum[nlow-1];
0238   y2 = fRuthSum[nlow];
0239   aa = (y2-y1)/(x2-x1);
0240   ruth = y1 + aa*(ee-x1);
0241
0242   return ruth;
0243 }
0244
0245 ///////////////////////////////////////////////
0246
0247 G4double G4LowPAIH2O::ComputeMeanFreePath( const G4ParticleDefinition* pd,
0248                                       G4double Tkin,
0249                                       const G4Material*,
0250                        G4double, // cutEnergy, // = 0.0,
0251                        G4double ) //maxEnergy ) //  = DBL_MAX)
0252 {
0253   G4double dndx(0.), mfp(DBL_MAX);
0254
0255   if     ( pd == theProton )    dndx = GetPrdNdx(Tkin);
0256   else if( pd == theElectron )  dndx = GetEldNdx(Tkin);
0257   else                          return DBL_MAX;
0258
0259   if( dndx > 0.) mfp = 1./dndx;
0260   else           mfp = DBL_MAX;
0261
0262   mfp *= fBias;
0263   return mfp;
0264 }
0265
0266 #endif
0267
0268 /////////////////   end of G4LowPAIH2O header file  //////////////////