include/Geant4/G4BaierKatkov.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 #ifndef G4BaierKatkov_h
0028 #define G4BaierKatkov_h 1
0029
0030 #include "globals.hh"
0031 #include <CLHEP/Units/SystemOfUnits.h>
0032 #include <vector>
0033 #include "G4ThreeVector.hh"
0034
0035 #include "G4VFastSimulationModel.hh"
0036
0037 /** \file G4BaierKatkov.hh
0038 * \brief Definition of the G4BaierKatkov class
0039 * This class is designed for the calculation of radiation probability, radiation point
0040 * and the parameters of the photon produced as well as spectrum accumulation using
0041 * the Baier-Katkov integral:
0042 * V. N. Baier, V. M. Katkov, and V. M. Strakhovenko,
0043 * Electromagnetic Processes at High Energies in Oriented Single Crystals
0044 * (World Scientific, Singapore, 1998).
0045 */
0046
0047 class G4BaierKatkov
0048 {
0049 public:
0050     // default constructor
0051     G4BaierKatkov();
0052
0053     // destructor
0054     ~G4BaierKatkov() = default;
0055
0056     /**
0057        You may call DoRadiation at each step of your trajectory
0058        CAUTION: please ensure that your steps are physically small enough for calculation
0059        of the radiation type you are interested in
0060        CAUTION: do ResetRadIntegral() before the start of a new trajectory
0061
0062        1) change some model defaults if necessary
0063          (SetSinglePhotonRadiationProbabilityLimit,
0064           SetNSmallTrajectorySteps, SetSpectrumEnergyRange)
0065        2) call DoRadiation at each step of your trajectory
0066        3) if DoRadiation returns TRUE, this means that a photon is produced (not added
0067           as a secondary yet) and its parameters are calculated.
0068        4) You may generate a new photon using GeneratePhoton either with
0069           the parameters calculated in DoRadiation or your own parameters.
0070           CAUTION: By now GeneratePhoton works only for a FastSim model
0071        5) Use GetPhotonEnergyInSpectrum() and GetTotalSpectrum() to return calculated
0072           total spectrum (all the photons altogether)
0073           Caution: is not normalized on the event number
0074        6) Get the charged particle parameters in the radiation point:
0075           GetParticleNewTotalEnergy(),
0076           GetParticleNewAngleX(), GetParticleNewAngleY(),
0077           GetNewGlobalTime(),
0078           GetParticleNewCoordinateXYZ()
0079      */
0080
0081     ///get functions
0082
0083     /// get maximal radiation probability to preserve single photon radiation
0084     G4double GetSinglePhotonRadiationProbabilityLimit()
0085     {return fSinglePhotonRadiationProbabilityLimit;}
0086
0087     ///CAUTION! : use the get functions below ONLY AFTER the call of DoRadiation
0088     /// and ONLY IF IT RETURNS true
0089
0090     ///total probability of radiation: needs calculation of DoRadiation first
0091     G4double GetTotalRadiationProbability(){return fTotalRadiationProbability;}
0092     ///get new parameters of the particle
0093     ///(the parameters at the point of radiation emission)
0094     ///needs calculation of DoRadiation first
0095     G4double GetParticleNewTotalEnergy(){return fNewParticleEnergy;}
0096     G4double GetParticleNewAngleX(){return fNewParticleAngleX;}
0097     G4double GetParticleNewAngleY(){return fNewParticleAngleY;}
0098     G4double GetNewGlobalTime(){return fNewGlobalTime;}
0099     const G4ThreeVector& GetParticleNewCoordinateXYZ(){return fNewParticleCoordinateXYZ;}
0100
0101   ///get photon energies (x-value in spectrum)
0102   const std::vector<G4double>& GetPhotonEnergyInSpectrum()
0103   {return fPhotonEnergyInSpectrum;}
0104
0105   ///get fTotalSpectrum after finishing the trajectory part with DoRadiation
0106   const std::vector<G4double>& GetTotalSpectrum(){return fTotalSpectrum;}
0107
0108     ///set functions
0109
0110     ///set maximal radiation probability to preserve single photon radiation
0111     void SetSinglePhotonRadiationProbabilityLimit(G4double wmax)
0112     {fSinglePhotonRadiationProbabilityLimit = wmax;}
0113
0114     ///number of steps in a trajectory small piece before
0115     ///the next call of the radiation integral
0116     void SetNSmallTrajectorySteps(G4double nSmallTrajectorySteps)
0117     {fNSmallTrajectorySteps = nSmallTrajectorySteps;}
0118
0119     ///reinitialize intermediate integrals fFa, fSs, fSc, fSsx, fSsy, fScx, fScy;
0120     ///reset radiation integral internal variables to defaults;
0121     ///reset the trajectory and radiation probability along the trajectory
0122     void ResetRadIntegral();
0123
0124     ///setting the number of photons in sampling of Baier-Katkov Integral
0125     ///(MC integration by photon energy and angles <=> photon momentum)
0126     void SetSamplingPhotonsNumber(G4int nPhotons){fNMCPhotons = nPhotons;}
0127
0128     ///setting the number of radiation angles 1/gamma, defining the width of
0129     ///the angular distribution of photon sampling in the Baier-Katkov Integral
0130     void SetRadiationAngleFactor(G4double radiationAngleFactor)
0131     {fRadiationAngleFactor = radiationAngleFactor;}
0132
0133     ///CAUTION, the bins width is logarithmic
0134     ///Do not worry if the maximal energy > particle energy.
0135     ///This elements of spectrum with non-physical energies
0136     ///will not be processed (they will be 0).
0137     void SetSpectrumEnergyRange(G4double emin,
0138                                 G4double emax,
0139                                 G4int numberOfBins);
0140     /// SetSpectrumEnergyRange also calls ResetRadIntegral()
0141
0142     void SetMinPhotonEnergy(G4double emin){SetSpectrumEnergyRange(emin,
0143                                                                   fMaxPhotonEnergy,
0144                                                                   fNBinsSpectrum);}
0145     void SetMaxPhotonEnergy(G4double emax){SetSpectrumEnergyRange(fMinPhotonEnergy,
0146                                                                   emax,
0147                                                                   fNBinsSpectrum);}
0148     void SetNBinsSpectrum(G4int nbin){SetSpectrumEnergyRange(fMinPhotonEnergy,
0149                                                              fMaxPhotonEnergy,
0150                                                              nbin);}
0151
0152     /// Increase the statistic of virtual photons in a certain energy region
0153     /// CAUTION! : don't do it before SetSpectrumEnergyRange or SetMinPhotonEnergy
0154     void AddStatisticsInPhotonEnergyRegion(G4double emin, G4double emax,
0155                                            G4int timesPhotonStatistics);
0156
0157     /// Virtual collimator masks the selection of photon angles in fTotalSpectrum
0158     /// Virtual collimator doesn't influence on Geant4 simulations.
0159     void SetVirtualCollimator(G4double virtualCollimatorAngularDiameter)
0160          {fVirtualCollimatorAngularDiameter=virtualCollimatorAngularDiameter;}
0161
0162     /// add the new elements of the trajectory, calculate radiation in a crystal
0163     /// see complete description in G4BaierKatkov::DoRadiation
0164     /// calls RadIntegral and all the necessary functions
0165     /// sets the parameters of a photon produced (if any)
0166     /// using SetPhotonProductionParameters()
0167     /// returns true in the case of photon generation, false if not
0168     G4bool DoRadiation(G4double etotal, G4double mass,
0169                        G4double angleX, G4double angleY,
0170                        G4double angleScatteringX, G4double angleScatteringY,
0171                        G4double step, G4double globalTime,
0172                        G4ThreeVector coordinateXYZ,
0173                        G4bool flagEndTrajectory=false);
0174
0175     /// generates secondary photon belonging to fastStep with variables
0176     /// photon energy, momentum direction, coordinates and global time
0177     /// CALCULATED IN DoRadiation => USE IT ONLY AFTER DoRadiation returns true
0178     void GeneratePhoton(G4FastStep &fastStep);
0179
0180 private:
0181
0182     ///set functions
0183
0184     ///function setting the photon sampling parameters in the Baier-Katkov integral;
0185     ///only the maximal energy is set, while fMinPhotonEnergy is used as a minimal energy;
0186     ///the angles set the angular distribution (the tails are infinite)
0187     void SetPhotonSamplingParameters(G4double ekin,
0188                                    G4double minPhotonAngleX, G4double maxPhotonAngleX,
0189                                    G4double minPhotonAngleY, G4double maxPhotonAngleY);
0190
0191     ///main functions:
0192
0193     ///generation of the photons in sampling of Baier-Katkov Integral
0194     ///(MC integration by photon energy and angles <=> by photon momentum)
0195     void GeneratePhotonSampling();
0196
0197     ///Baier-Katkov method: calculation of integral, spectrum, full probability;
0198     ///returns the total radiation probability;
0199     ///calculates the radiation spectrum on this trajectory piece
0200     G4double RadIntegral(G4double etotal, G4double mass,
0201                          std::vector<G4double> &vectorParticleAnglesX,
0202                          std::vector<G4double> &vectorParticleAnglesY,
0203                          std::vector<G4double> &vectorScatteringAnglesX,
0204                          std::vector<G4double> &vectorScatteringAnglesY,
0205                          std::vector<G4double> &vectorSteps,
0206                          G4int imin);
0207
0208     ///set photon production parameters (returns false if no photon produced)
0209     ///accumulates fTotalSpectrum
0210     ///CAUTION: it is an accessory function of DoRadiation, do not use it separately
0211     G4bool SetPhotonProductionParameters(G4double etotal, G4double mass);
0212
0213     G4int FindVectorIndex(std::vector<G4double> &myvector, G4double value);
0214
0215     G4double fTotalRadiationProbability = 0.;
0216     G4double fSinglePhotonRadiationProbabilityLimit=0.25;//Maximal radiation
0217                                          //probability to preserve single photon radiation
0218
0219     //number of steps in a trajectory piece before the next call of the radiation integral
0220     G4int fNSmallTrajectorySteps=10000;
0221     ///trajectory element No (the first element of the array feeded in RadIntegral)
0222     G4int fImin0 = 0;
0223     ///Monte Carlo statistics of photon sampling in Baier-Katkov with 1 trajectory
0224     G4int fNMCPhotons =150;
0225     ///the number of bins in photon spectrum
0226     G4int fNBinsSpectrum = 110;
0227     G4double fMinPhotonEnergy = 0.1*CLHEP::MeV;//min energy in spectrum output
0228     G4double fMaxPhotonEnergy = 1*CLHEP::GeV;  //max energy in spectrum output
0229     G4double fLogEmaxdEmin = 1.;// = log(fMaxPhotonEnergy/fMinPhotonEnergy),
0230                                // 1/normalizing coefficient in
0231                                // 1/E distribution between
0232                                // fMinPhotonEnergy and fMaxPhotonEnergy
0233                                // is used only for spectrum output, not for simulations
0234                                //(we take bremsstrahlung for photon sampling)
0235     G4double fLogEdEmin = 1.;   // = log(E/fMinPhotonEnergy), the same as fLogEmaxdEmin
0236                                // but with the particle energy as the maximal limit
0237
0238     G4double fVirtualCollimatorAngularDiameter=1.;//default, infinite angle
0239     std::vector<G4bool> fInsideVirtualCollimator;
0240
0241     ///data of the phootn energy range with additional statistics
0242     std::vector<G4double> fLogAddRangeEmindEmin;//=G4Log(emin/fMinPhotonEnergy)
0243     std::vector<G4double> fLogAddRangeEmaxdEmin;//=G4Log(emax/fMinPhotonEnergy)
0244     std::vector<G4int> fTimesPhotonStatistics;
0245
0246     ///number of trajectories
0247     //(at each of the Baier-Katkov Integral is calculated for the same photons)
0248     G4int fItrajectories = 0;
0249
0250     G4double fEph0=0;   //energy of the photon produced
0251     G4ThreeVector PhMomentumDirection;   //momentum direction of the photon produced
0252
0253     ///Radiation integral variables
0254     G4double fMeanPhotonAngleX =0.;        //average angle of radiated photon direction
0255                                            //in sampling, x-plane
0256     G4double fParamPhotonAngleX=1.e-3*CLHEP::rad; //a parameter radiated photon
0257                                            //sampling distribution, x-plane
0258     G4double fMeanPhotonAngleY =0.;        //average angle of radiated photon direction
0259                                            //in sampling, y-plane
0260     G4double fParamPhotonAngleY=1.e-3*CLHEP::rad; //a parameter radiated photon
0261                                            //sampling distribution, y-plane
0262     G4double fRadiationAngleFactor = 1.; // number of radiation angles 1/gamma:
0263                                          // more fRadiationAngleFactor =>
0264                                          // higher fParamPhotonAngleX and Y
0265
0266     ///new particle parameters (the parameters at the point of radiation emission)
0267     G4double fNewParticleEnergy=0;
0268     G4double fNewParticleAngleX=0;
0269     G4double fNewParticleAngleY=0;
0270     G4double fNewGlobalTime=0;
0271     G4ThreeVector fNewParticleCoordinateXYZ;
0272
0273     ///sampling of the energy and the angles of a photon emission
0274     ///(integration variables, Monte Carlo integration)
0275     std::vector<G4double> fPhotonEnergyInIntegral;
0276     std::vector<G4double> fPhotonAngleInIntegralX;
0277     std::vector<G4double> fPhotonAngleInIntegralY;
0278     std::vector<G4double> fPhotonAngleNormCoef;
0279     ///spectrum bin index for each photon
0280     std::vector<G4double> fIBinsSpectrum;
0281     ///the vector of the discrete CDF of the radiation of sampling photons
0282     std::vector<G4double> fPhotonProductionCDF;
0283
0284     ///vectors of the trajectory
0285     std::vector<G4double> fParticleAnglesX;
0286     std::vector<G4double> fParticleAnglesY;
0287     std::vector<G4double> fScatteringAnglesX;
0288     std::vector<G4double> fScatteringAnglesY;
0289     std::vector<G4double> fSteps;
0290     std::vector<G4double> fGlobalTimes;
0291     std::vector<G4ThreeVector> fParticleCoordinatesXYZ;
0292
0293     ///intermediate integrals (different for each photon energy value)!!!
0294     std::vector<G4double> fFa;//phase
0295     std::vector<G4double> fSs;
0296     std::vector<G4double> fSc;
0297     std::vector<G4double> fSsx;
0298     std::vector<G4double> fSsy;
0299     std::vector<G4double> fScx;
0300     std::vector<G4double> fScy;
0301
0302     ///output
0303     std::vector<G4double> fPhotonEnergyInSpectrum; //energy values in spectrum
0304
0305     std::vector<G4int> fNPhotonsPerBin; //number of photons per spectrum bin
0306                                    //(accumulating during total run)
0307
0308     std::vector<G4double> fSpectrum; //spectrum normalized by the total
0309                         //radiation probability of one particle at one call of RadIntegral
0310
0311     std::vector<std::vector<G4double>> fAccumSpectrum; //accumulate Spectrum during
0312                                              //the part of a trajectory
0313
0314     std::vector<G4double> fAccumTotalSpectrum; //spectrum normalized by the total
0315                                           //radiation probability summed
0316                                           //for all the particles (is not divided
0317                                           //of one particle number fNPhotonsPerBin)
0318
0319     std::vector<G4double> fTotalSpectrum; //spectrum normalized by
0320                                      //the total radiation probability summed
0321                                      //for all the particles
0322                                      //(is divided by the photon number fNPhotonsPerBin)
0323                                      //multiplied by the number of trajectories
0324                                      //(fItrajectories)
0325
0326     std::vector<G4double> fImax0; //trajectory element numbers at the end of each
0327                             //small piece; G4double just for security of some operations
0328     ///total radiation probability along this trajectory
0329     std::vector<G4double> fTotalRadiationProbabilityAlongTrajectory;
0330 };
0331
0332 #endif