Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/parameterisations/channeling/include/G4BaierKatkov.hh

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // Author:      Alexei Sytov
 27 // Co-author:   Gianfranco Paterno (modifications & testing)
 28 // On the base of the CRYSTALRAD realization of the Baier-Katkov integral:
 29 // A. I. Sytov, V. V. Tikhomirov, and L. Bandiera PRAB 22, 064601 (2019)
 30 
 31 #ifndef G4BaierKatkov_h
 32 #define G4BaierKatkov_h 1
 33 
 34 #include "globals.hh"
 35 #include <CLHEP/Units/SystemOfUnits.h>
 36 #include <vector>
 37 #include "G4ThreeVector.hh"
 38 
 39 #include "G4VFastSimulationModel.hh"
 40 
 41 /** \file G4BaierKatkov.hh
 42 * \brief Definition of the G4BaierKatkov class
 43 * This class is designed for the calculation of radiation probability, radiation point
 44 * and the parameters of the photon produced as well as spectrum accumulation using
 45 * the Baier-Katkov integral:
 46 * V. N. Baier, V. M. Katkov, and V. M. Strakhovenko,
 47 * Electromagnetic Processes at High Energies in Oriented Single Crystals
 48 * (World Scientific, Singapore, 1998).
 49 */
 50 
 51 class G4BaierKatkov
 52 {
 53 public:
 54     // default constructor
 55     G4BaierKatkov();
 56 
 57     // destructor
 58     ~G4BaierKatkov() = default;
 59 
 60     /**
 61        You may call DoRadiation at each step of your trajectory
 62        CAUTION: please ensure that your steps are physically small enough for calculation
 63        of the radiation type you are interested in
 64        CAUTION: do ResetRadIntegral() before the start of a new trajectory
 65 
 66        1) change some model defaults if necessary
 67          (SetSinglePhotonRadiationProbabilityLimit,
 68           SetNSmallTrajectorySteps, SetSpectrumEnergyRange)
 69        2) call DoRadiation at each step of your trajectory
 70        3) if DoRadiation returns TRUE, this means that a photon is produced (not added
 71           as a secondary yet) and its parameters are calculated.
 72        4) You may generate a new photon using GeneratePhoton either with
 73           the parameters calculated in DoRadiation or your own parameters.
 74           CAUTION: By now GeneratePhoton works only for a FastSim model
 75        5) Use GetPhotonEnergyInSpectrum() and GetTotalSpectrum() to return calculated
 76           total spectrum (all the photons altogether)
 77           Caution: is not normalized on the event number
 78        6) Get the charged particle parameters in the radiation point:
 79           GetParticleNewTotalEnergy(),
 80           GetParticleNewAngleX(), GetParticleNewAngleY(),
 81           GetNewGlobalTime(),
 82           GetParticleNewCoordinateXYZ()
 83      */
 84 
 85     ///get functions
 86 
 87     /// get maximal radiation probability to preserve single photon radiation
 88     G4double GetSinglePhotonRadiationProbabilityLimit()
 89     {return fSinglePhotonRadiationProbabilityLimit;}
 90 
 91     ///CAUTION! : use the get functions below ONLY AFTER the call of DoRadiation
 92     /// and ONLY IF IT RETURNS true
 93 
 94     ///total probability of radiation: needs calculation of DoRadiation first
 95     G4double GetTotalRadiationProbability(){return fTotalRadiationProbability;}
 96     ///get new parameters of the particle
 97     ///(the parameters at the point of radiation emission)
 98     ///needs calculation of DoRadiation first
 99     G4double GetParticleNewTotalEnergy(){return fNewParticleEnergy;}
100     G4double GetParticleNewAngleX(){return fNewParticleAngleX;}
101     G4double GetParticleNewAngleY(){return fNewParticleAngleY;}
102     G4double GetNewGlobalTime(){return fNewGlobalTime;}
103     const G4ThreeVector& GetParticleNewCoordinateXYZ(){return fNewParticleCoordinateXYZ;}
104 
105   ///get photon energies (x-value in spectrum)
106   const std::vector<G4double>& GetPhotonEnergyInSpectrum()
107   {return fPhotonEnergyInSpectrum;}
108   
109   ///get fTotalSpectrum after finishing the trajectory part with DoRadiation
110   const std::vector<G4double>& GetTotalSpectrum(){return fTotalSpectrum;}
111 
112     ///set functions
113 
114     ///set maximal radiation probability to preserve single photon radiation
115     void SetSinglePhotonRadiationProbabilityLimit(G4double wmax)
116     {fSinglePhotonRadiationProbabilityLimit = wmax;}
117 
118     ///number of steps in a trajectory small piece before
119     ///the next call of the radiation integral
120     void SetNSmallTrajectorySteps(G4int nSmallTrajectorySteps)
121     {fNSmallTrajectorySteps = nSmallTrajectorySteps;}
122 
123     ///reinitialize intermediate integrals fFa, fSs, fSc, fSsx, fSsy, fScx, fScy;
124     ///reset radiation integral internal variables to defaults;
125     ///reset the trajectory and radiation probability along the trajectory
126     void ResetRadIntegral();
127 
128     ///setting the number of photons in sampling of Baier-Katkov Integral
129     ///(MC integration by photon energy and angles <=> photon momentum)
130     void SetSamplingPhotonsNumber(G4int nPhotons){fNMCPhotons = nPhotons;}
131 
132     ///setting the number of radiation angles 1/gamma, defining the width of
133     ///the angular distribution of photon sampling in the Baier-Katkov Integral
134     void SetRadiationAngleFactor(G4double radiationAngleFactor)
135     {fRadiationAngleFactor = radiationAngleFactor;}
136 
137     ///CAUTION, the bins width is logarithmic
138     ///Do not worry if the maximal energy > particle energy.
139     ///This elements of spectrum with non-physical energies
140     ///will not be processed (they will be 0).
141     void SetSpectrumEnergyRange(G4double emin,
142                                 G4double emax,
143                                 G4int numberOfBins);
144     /// SetSpectrumEnergyRange also calls ResetRadIntegral()
145 
146     void SetMinPhotonEnergy(G4double emin){SetSpectrumEnergyRange(emin,
147                                                                   fMaxPhotonEnergy,
148                                                                   fNBinsSpectrum);}
149     void SetMaxPhotonEnergy(G4double emax){SetSpectrumEnergyRange(fMinPhotonEnergy,
150                                                                   emax,
151                                                                   fNBinsSpectrum);}
152     void SetNBinsSpectrum(G4int nbin){SetSpectrumEnergyRange(fMinPhotonEnergy,
153                                                              fMaxPhotonEnergy,
154                                                              nbin);}
155 
156     /// Increase the statistic of virtual photons in a certain energy region
157     /// CAUTION! : don't do it before SetSpectrumEnergyRange or SetMinPhotonEnergy
158     void AddStatisticsInPhotonEnergyRegion(G4double emin, G4double emax,
159                                            G4int timesPhotonStatistics);
160 
161     /// Virtual collimator masks the selection of photon angles in fTotalSpectrum
162     /// Virtual collimator doesn't influence on Geant4 simulations.
163     void SetVirtualCollimator(G4double virtualCollimatorAngularDiameter)
164          {fVirtualCollimatorAngularDiameter=virtualCollimatorAngularDiameter;}
165 
166     /// add the new elements of the trajectory, calculate radiation in a crystal
167     /// see complete description in G4BaierKatkov::DoRadiation
168     /// calls RadIntegral and all the necessary functions
169     /// sets the parameters of a photon produced (if any)
170     /// using SetPhotonProductionParameters()
171     /// returns true in the case of photon generation, false if not
172     G4bool DoRadiation(G4double etotal, G4double mass,
173                        G4double angleX, G4double angleY,
174                        G4double angleScatteringX, G4double angleScatteringY,
175                        G4double step, G4double globalTime,
176                        G4ThreeVector coordinateXYZ,
177                        G4bool flagEndTrajectory=false);
178 
179     /// generates secondary photon belonging to fastStep with variables
180     /// photon energy, momentum direction, coordinates and global time
181     /// CALCULATED IN DoRadiation => USE IT ONLY AFTER DoRadiation returns true
182     void GeneratePhoton(G4FastStep &fastStep);
183 
184 private:
185 
186     ///set functions
187 
188     ///function setting the photon sampling parameters in the Baier-Katkov integral;
189     ///only the maximal energy is set, while fMinPhotonEnergy is used as a minimal energy;
190     ///the angles set the angular distribution (the tails are infinite)
191     void SetPhotonSamplingParameters(G4double ekin,
192                                    G4double minPhotonAngleX, G4double maxPhotonAngleX,
193                                    G4double minPhotonAngleY, G4double maxPhotonAngleY);
194 
195     ///main functions:
196 
197     ///generation of the photons in sampling of Baier-Katkov Integral
198     ///(MC integration by photon energy and angles <=> by photon momentum)
199     void GeneratePhotonSampling();
200 
201     ///Baier-Katkov method: calculation of integral, spectrum, full probability;
202     ///returns the total radiation probability;
203     ///calculates the radiation spectrum on this trajectory piece
204     G4double RadIntegral(G4double etotal, G4double mass,
205                          std::vector<G4double> &vectorParticleAnglesX,
206                          std::vector<G4double> &vectorParticleAnglesY,
207                          std::vector<G4double> &vectorScatteringAnglesX,
208                          std::vector<G4double> &vectorScatteringAnglesY,
209                          std::vector<G4double> &vectorSteps,
210                          G4int imin);
211 
212     ///set photon production parameters (returns false if no photon produced)
213     ///accumulates fTotalSpectrum
214     ///CAUTION: it is an accessory function of DoRadiation, do not use it separately
215     G4bool SetPhotonProductionParameters(G4double etotal, G4double mass);
216 
217     G4int FindVectorIndex(std::vector<G4double> &myvector, G4double value);
218 
219     G4double fTotalRadiationProbability = 0.;
220     G4double fSinglePhotonRadiationProbabilityLimit=0.25;//Maximal radiation
221                                          //probability to preserve single photon radiation
222 
223     //number of steps in a trajectory piece before the next call of the radiation integral
224     G4int fNSmallTrajectorySteps=10000;
225     ///trajectory element No (the first element of the array feeded in RadIntegral)
226     G4int fImin0 = 0;
227     ///Monte Carlo statistics of photon sampling in Baier-Katkov with 1 trajectory
228     G4int fNMCPhotons =150;
229     ///the number of bins in photon spectrum
230     G4int fNBinsSpectrum = 110;
231     G4double fMinPhotonEnergy = 0.1*CLHEP::MeV;//min energy in spectrum output
232     G4double fMaxPhotonEnergy = 1*CLHEP::GeV;  //max energy in spectrum output
233     G4double fLogEmaxdEmin = 1.;// = log(fMaxPhotonEnergy/fMinPhotonEnergy),
234                                // 1/normalizing coefficient in
235                                // 1/E distribution between
236                                // fMinPhotonEnergy and fMaxPhotonEnergy
237                                // is used only for spectrum output, not for simulations
238                                //(we take bremsstrahlung for photon sampling)
239     G4double fLogEdEmin = 1.;   // = log(E/fMinPhotonEnergy), the same as fLogEmaxdEmin
240                                // but with the particle energy as the maximal limit
241 
242     G4double fVirtualCollimatorAngularDiameter=1.;//default, infinite angle
243     std::vector<G4bool> fInsideVirtualCollimator;
244 
245     ///data of the phootn energy range with additional statistics
246     std::vector<G4double> fLogAddRangeEmindEmin;//=G4Log(emin/fMinPhotonEnergy)
247     std::vector<G4double> fLogAddRangeEmaxdEmin;//=G4Log(emax/fMinPhotonEnergy)
248     std::vector<G4int> fTimesPhotonStatistics;
249 
250     ///number of trajectories
251     //(at each of the Baier-Katkov Integral is calculated for the same photons)
252     G4int fItrajectories = 0;
253 
254     G4double fEph0=0;   //energy of the photon produced
255     G4ThreeVector PhMomentumDirection;   //momentum direction of the photon produced
256 
257     ///Radiation integral variables
258     G4double fMeanPhotonAngleX =0.;        //average angle of radiated photon direction
259                                            //in sampling, x-plane
260     G4double fParamPhotonAngleX=1.e-3*CLHEP::rad; //a parameter radiated photon
261                                            //sampling distribution, x-plane
262     G4double fMeanPhotonAngleY =0.;        //average angle of radiated photon direction
263                                            //in sampling, y-plane
264     G4double fParamPhotonAngleY=1.e-3*CLHEP::rad; //a parameter radiated photon
265                                            //sampling distribution, y-plane
266     G4double fRadiationAngleFactor = 4.; // number of radiation angles 1/gamma:
267                                          // more fRadiationAngleFactor =>
268                                          // higher fParamPhotonAngleX and Y
269 
270     ///new particle parameters (the parameters at the point of radiation emission)
271     G4double fNewParticleEnergy=0;
272     G4double fNewParticleAngleX=0;
273     G4double fNewParticleAngleY=0;
274     G4double fNewGlobalTime=0;
275     G4ThreeVector fNewParticleCoordinateXYZ;
276 
277     ///sampling of the energy and the angles of a photon emission
278     ///(integration variables, Monte Carlo integration)
279     std::vector<G4double> fPhotonEnergyInIntegral;
280     std::vector<G4double> fPhotonAngleInIntegralX;
281     std::vector<G4double> fPhotonAngleInIntegralY;
282     std::vector<G4double> fPhotonAngleNormCoef;
283     ///spectrum bin index for each photon
284     std::vector<G4double> fIBinsSpectrum;
285     ///the vector of the discrete CDF of the radiation of sampling photons
286     std::vector<G4double> fPhotonProductionCDF;
287 
288     ///vectors of the trajectory
289     std::vector<G4double> fParticleAnglesX;
290     std::vector<G4double> fParticleAnglesY;
291     std::vector<G4double> fScatteringAnglesX;
292     std::vector<G4double> fScatteringAnglesY;
293     std::vector<G4double> fSteps;
294     std::vector<G4double> fGlobalTimes;
295     std::vector<G4ThreeVector> fParticleCoordinatesXYZ;
296 
297     ///intermediate integrals (different for each photon energy value)!!!
298     std::vector<G4double> fFa;//phase
299     std::vector<G4double> fSs;
300     std::vector<G4double> fSc;
301     std::vector<G4double> fSsx;
302     std::vector<G4double> fSsy;
303     std::vector<G4double> fScx;
304     std::vector<G4double> fScy;
305 
306     ///output
307     std::vector<G4double> fPhotonEnergyInSpectrum; //energy values in spectrum
308 
309     std::vector<G4int> fNPhotonsPerBin; //number of photons per spectrum bin
310                                    //(accumulating during total run)
311 
312     std::vector<G4double> fSpectrum; //spectrum normalized by the total
313                         //radiation probability of one particle at one call of RadIntegral
314 
315     std::vector<std::vector<G4double>> fAccumSpectrum; //accumulate Spectrum during
316                                              //the part of a trajectory
317 
318     std::vector<G4double> fAccumTotalSpectrum; //spectrum normalized by the total
319                                           //radiation probability summed
320                                           //for all the particles (is not divided
321                                           //of one particle number fNPhotonsPerBin)
322 
323     std::vector<G4double> fTotalSpectrum; //spectrum normalized by
324                                      //the total radiation probability summed
325                                      //for all the particles
326                                      //(is divided by the photon number fNPhotonsPerBin)
327                                      //multiplied by the number of trajectories
328                                      //(fItrajectories)
329 
330     std::vector<G4double> fImax0; //trajectory element numbers at the end of each
331                             //small piece; G4double just for security of some operations
332     ///total radiation probability along this trajectory
333     std::vector<G4double> fTotalRadiationProbabilityAlongTrajectory;
334 };
335 
336 #endif
337