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