include/Geant4/G4UrbanMscModel.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 // -------------------------------------------------------------------
0028 //
0029 //
0030 // GEANT4 Class header file
0031 //
0032 //
0033 // File name:     G4UrbanMscModel
0034 //
0035 // Author:        Laszlo Urban
0036 //
0037 // Creation date: 19.02.2013
0038 //
0039 // Created from G4UrbanMscModel96
0040 //
0041 // New parametrization for theta0
0042 // Correction for very small step length
0043 //
0044 // Class Description:
0045 //
0046 // Implementation of the model of multiple scattering based on
0047 // H.W.Lewis Phys Rev 78 (1950) 526 and L.Urban model
0048
0049 // -------------------------------------------------------------------
0050 //
0051
0052 #ifndef G4UrbanMscModel_h
0053 #define G4UrbanMscModel_h 1
0054
0055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0056
0057 #include <CLHEP/Units/SystemOfUnits.h>
0058
0059 #include "G4VMscModel.hh"
0060 #include "G4MscStepLimitType.hh"
0061 #include "G4Log.hh"
0062 #include "G4Exp.hh"
0063
0064 class G4ParticleChangeForMSC;
0065 class G4SafetyHelper;
0066
0067 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0068
0069 class G4UrbanMscModel : public G4VMscModel
0070 {
0071
0072 public:
0073
0074   explicit G4UrbanMscModel(const G4String& nam = "UrbanMsc");
0075
0076   ~G4UrbanMscModel() override;
0077
0078   void Initialise(const G4ParticleDefinition*,
0079           const G4DataVector&) override;
0080
0081   void StartTracking(G4Track*) override;
0082
0083   G4double
0084   ComputeCrossSectionPerAtom(const G4ParticleDefinition* particle,
0085                  G4double KineticEnergy,
0086                  G4double AtomicNumber,
0087                  G4double AtomicWeight=0.,
0088                  G4double cut =0.,
0089                  G4double emax=DBL_MAX) override;
0090
0091   G4ThreeVector& SampleScattering(const G4ThreeVector&,
0092                   G4double safety) override;
0093
0094   G4double ComputeTruePathLengthLimit(const G4Track& track,
0095                           G4double& currentMinimalStep) override;
0096
0097   G4double ComputeGeomPathLength(G4double truePathLength) override;
0098
0099   G4double ComputeTrueStepLength(G4double geomStepLength) override;
0100
0101   G4double ComputeTheta0(G4double truePathLength, G4double KineticEnergy);
0102
0103   inline void SetDisplacementAlgorithm96(const G4bool val);
0104
0105   inline void SetPositronCorrection(const G4bool val);
0106
0107   //  hide assignment operator
0108   G4UrbanMscModel & operator=(const  G4UrbanMscModel &right) = delete;
0109   G4UrbanMscModel(const  G4UrbanMscModel&) = delete;
0110
0111 private:
0112
0113   G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy);
0114
0115   void SampleDisplacement(G4double sinTheta, G4double phi);
0116
0117   void SampleDisplacementNew(G4double sinTheta, G4double phi);
0118
0119   void InitialiseModelCache();
0120
0121   inline void SetParticle(const G4ParticleDefinition*);
0122
0123   inline G4double Randomizetlimit();
0124
0125   inline G4double SimpleScattering();
0126
0127   inline G4double ComputeStepmin();
0128
0129   inline G4double ComputeTlimitmin();
0130
0131   CLHEP::HepRandomEngine* rndmEngineMod;
0132
0133   const G4ParticleDefinition* particle = nullptr;
0134   const G4ParticleDefinition* positron;
0135   G4ParticleChangeForMSC* fParticleChange = nullptr;
0136   const G4MaterialCutsCouple* couple = nullptr;
0137
0138   G4double mass;
0139   G4double charge,chargeSquare;
0140   G4double masslimite,fr;
0141
0142   G4double taubig;
0143   G4double tausmall;
0144   G4double taulim;
0145   G4double currentTau;
0146   G4double tlimit;
0147   G4double tlimitmin;
0148   G4double tlimitminfix,tlimitminfix2;
0149   G4double tgeom;
0150
0151   G4double geombig;
0152   G4double geommin;
0153   G4double geomlimit;
0154   G4double skindepth;
0155   G4double smallstep;
0156
0157   G4double presafety;
0158
0159   G4double lambda0;
0160   G4double lambdaeff;
0161   G4double tPathLength;
0162   G4double zPathLength;
0163   G4double par1,par2,par3;
0164
0165   G4double stepmin;
0166
0167   G4double currentKinEnergy;
0168   G4double currentLogKinEnergy;
0169   G4double currentRange;
0170   G4double rangeinit;
0171   G4double currentRadLength;
0172
0173   G4double drr,finalr;
0174
0175   G4double tlow;
0176   G4double invmev;
0177   G4double xmeanth = 0.0;
0178   G4double x2meanth = 1./3.;
0179   G4double rndmarray[2];
0180
0181   struct mscData {
0182     G4double Z23, sqrtZ, factmin;
0183     G4double coeffth1, coeffth2;
0184     G4double coeffc1, coeffc2, coeffc3, coeffc4;
0185     G4double stepmina, stepminb;
0186     G4double doverra, doverrb;
0187     G4double posa, posb, posc, posd, pose;
0188   };
0189   static std::vector<mscData*> msc;
0190
0191   // index of G4MaterialCutsCouple
0192   G4int idx = 0;
0193
0194   G4bool firstStep = true;
0195   G4bool insideskin = false;
0196
0197   G4bool latDisplasmentbackup = false;
0198   G4bool dispAlg96 = true;
0199   G4bool fPosiCorrection = true;
0200   G4bool isFirstInstance = false;
0201 };
0202
0203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0205
0206 inline
0207 void G4UrbanMscModel::SetParticle(const G4ParticleDefinition* p)
0208 {
0209   if (p != particle) {
0210     particle = p;
0211     mass = p->GetPDGMass();
0212     charge = p->GetPDGCharge()/CLHEP::eplus;
0213     chargeSquare = charge*charge;
0214   }
0215 }
0216
0217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0218
0219 inline G4double G4UrbanMscModel::Randomizetlimit()
0220 {
0221   G4double res = tlimitmin;
0222   if(tlimit > tlimitmin)
0223   {
0224     res = G4RandGauss::shoot(rndmEngineMod,tlimit,0.1*(tlimit-tlimitmin));
0225     res = std::max(res, tlimitmin);
0226   }
0227   return res;
0228 }
0229
0230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0231
0232 inline G4double G4UrbanMscModel::SimpleScattering()
0233 {
0234   // 'large angle scattering'
0235   // 2 model functions with correct xmean and x2mean
0236   const G4double a = (2.*xmeanth+9.*x2meanth-3.)/(2.*xmeanth-3.*x2meanth+1.);
0237   const G4double prob = (a+2.)*xmeanth/a;
0238
0239   // sampling
0240   rndmEngineMod->flatArray(2, rndmarray);
0241   return (rndmarray[1] < prob) ?
0242     -1.+2.*G4Exp(G4Log(rndmarray[0])/(a+1.)) : -1.+2.*rndmarray[0];
0243 }
0244
0245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0246
0247 inline G4double G4UrbanMscModel::ComputeStepmin()
0248 {
0249   // define stepmin using estimation of the ratio
0250   // of lambda_elastic/lambda_transport
0251   const G4double rat = currentKinEnergy*invmev;
0252   return lambda0*msc[idx]->factmin/
0253     (0.002 + rat*(msc[idx]->stepmina + msc[idx]->stepminb*rat));
0254 }
0255
0256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0257
0258 inline G4double G4UrbanMscModel::ComputeTlimitmin()
0259 {
0260   G4double x = (particle == positron) ?
0261     0.7*msc[idx]->sqrtZ*stepmin : 0.87*msc[idx]->Z23*stepmin;
0262   if(currentKinEnergy < tlow) { x *= 0.5*(1.+currentKinEnergy/tlow); }
0263   return std::max(x, tlimitminfix);
0264 }
0265
0266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0267
0268 inline void G4UrbanMscModel::SetDisplacementAlgorithm96(const G4bool val)
0269 {
0270   dispAlg96 = val;
0271 }
0272
0273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0274
0275 inline void G4UrbanMscModel::SetPositronCorrection(const G4bool val)
0276 {
0277   fPosiCorrection = val;
0278 }
0279
0280 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0281
0282 #endif
0283