include/Geant4/G4UniversalFluctuation.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 // GEANT4 Class header file
0029 //
0030 //
0031 // File name:     G4UniversalFluctuation
0032 //
0033 // Author:        V.Ivanchenko make a class with the Laszlo Urban model
0034 //
0035 // Creation date: 03.01.2002
0036 //
0037 // Modifications:
0038 //
0039 //
0040 // Class Description:
0041 //
0042 // Implementation of energy loss fluctuations made by L.Urban in 2021
0043
0044 // -------------------------------------------------------------------
0045 //
0046
0047 #ifndef G4UniversalFluctuation_h
0048 #define G4UniversalFluctuation_h 1
0049
0050 #include "G4VEmFluctuationModel.hh"
0051 #include "G4ParticleDefinition.hh"
0052 #include "G4Poisson.hh"
0053 #include <CLHEP/Random/RandomEngine.h>
0054
0055 class G4UniversalFluctuation : public G4VEmFluctuationModel
0056 {
0057
0058 public:
0059
0060   explicit G4UniversalFluctuation(const G4String& nam = "UniFluc");
0061
0062   ~G4UniversalFluctuation() override;
0063
0064   G4double SampleFluctuations(const G4MaterialCutsCouple*,
0065                   const G4DynamicParticle*,
0066                               const G4double, const G4double,
0067                   const G4double, const G4double) override;
0068
0069   G4double Dispersion(const G4Material*,
0070               const G4DynamicParticle*,
0071                       const G4double, const G4double,
0072                       const G4double) override;
0073
0074   // Initialisation for a new particle type
0075   void InitialiseMe(const G4ParticleDefinition*) override;
0076
0077   // Initialisation prestep
0078   void SetParticleAndCharge(const G4ParticleDefinition*,
0079                 G4double q2) override;
0080
0081   // hide assignment operator
0082   G4UniversalFluctuation & operator=
0083   (const G4UniversalFluctuation &right) = delete;
0084   G4UniversalFluctuation(const G4UniversalFluctuation&) = delete;
0085
0086 protected:
0087
0088   virtual G4double SampleGlandz(CLHEP::HepRandomEngine* rndm,
0089                                 const G4Material*, const G4double tcut);
0090
0091   inline void AddExcitation(CLHEP::HepRandomEngine* rndm,
0092                             const G4double ax, const G4double ex,
0093                             G4double& eav,
0094                             G4double& eloss, G4double& esig2);
0095
0096   inline void SampleGauss(CLHEP::HepRandomEngine* rndm,
0097                           const G4double eav, const G4double esig2,
0098                           G4double& eloss);
0099
0100   // particle properties
0101   G4double particleMass = 0.0;
0102   G4double m_Inv_particleMass = DBL_MAX;
0103   G4double m_massrate = DBL_MAX;
0104   G4double chargeSquare = 1.0;
0105
0106   // material properties
0107   G4double ipotFluct = 0.0;
0108   G4double ipotLogFluct = 0.0;
0109   G4double e0 = 0.0;
0110
0111   // model parameters
0112   G4double minNumberInteractionsBohr = 10.0;
0113   G4double minLoss;
0114   G4double nmaxCont = 8.0;
0115   G4double rate = 0.56;
0116   G4double fw = 4.0;
0117   G4double a0 = 42.0;
0118   G4double w2 = 0.0;
0119   G4double meanLoss = 0.0;
0120
0121   const G4ParticleDefinition* particle = nullptr;
0122   G4double* rndmarray = nullptr;
0123   G4int sizearray = 30;
0124 };
0125
0126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0127
0128 inline void
0129 G4UniversalFluctuation::AddExcitation(CLHEP::HepRandomEngine* rndm,
0130                                       const G4double ax, const G4double ex,
0131                                       G4double& eav,
0132                                       G4double& eloss, G4double& esig2)
0133 {
0134   if(ax > nmaxCont) {
0135     eav  += ax*ex;
0136     esig2 += ax*ex*ex;
0137   } else {
0138     const G4int p = (G4int)G4Poisson(ax);
0139     if(p > 0) { eloss += ((p + 1) - 2.*rndm->flat())*ex; }
0140   }
0141 }
0142
0143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0144
0145 inline void
0146 G4UniversalFluctuation::SampleGauss(CLHEP::HepRandomEngine* rndm,
0147                                     const G4double eav, const G4double esig2,
0148                                     G4double& eloss)
0149 {
0150   G4double x = eav;
0151   const G4double sig = std::sqrt(esig2);
0152   if(eav < 0.25*sig) {
0153     x += (2.*rndm->flat() - 1.)*eav;
0154   } else {
0155     do {
0156       x = G4RandGauss::shoot(rndm, eav, sig);
0157     } while (x < 0.0 || x > 2*eav);
0158     // Loop checking, 23-Feb-2016, Vladimir Ivanchenko
0159   }
0160   eloss += x;
0161 }
0162
0163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0164
0165 #endif
0166