Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4BaierKatkov.cc
Go to the documentation of this file.
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
27#include "G4BaierKatkov.hh"
28
29#include "Randomize.hh"
30#include "G4Gamma.hh"
31#include "G4SystemOfUnits.hh"
32#include "G4PhysicalConstants.hh"
33#include "G4Log.hh"
34
35G4BaierKatkov::G4BaierKatkov()
36{
37 //sets the default spectrum energy range of integration and
38 //calls ResetRadIntegral()
39 SetSpectrumEnergyRange(0.1*MeV,1.*GeV,110);
40
41 //Do not worry if the maximal energy > particle energy
42 //this elements of spectrum with non-physical energies
43 //will not be processed (they will be 0)
44
45 G4cout << " "<< G4endl;
46 G4cout << "G4BaierKatkov model is activated."<< G4endl;
47 G4cout << " "<< G4endl;
48}
49
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51
52void G4BaierKatkov::ResetRadIntegral()
53{
54 fAccumSpectrum.clear();
55
56 //reinitialize intermediate integrals with zeros
57 fFa.resize(fNMCPhotons);
58 fSs.resize(fNMCPhotons);
59 fSc.resize(fNMCPhotons);
60 fSsx.resize(fNMCPhotons);
61 fSsy.resize(fNMCPhotons);
62 fScx.resize(fNMCPhotons);
63 fScy.resize(fNMCPhotons);
64 std::fill(fFa.begin(), fFa.end(), 0.);
65 std::fill(fSs.begin(), fSs.end(), 0.);
66 std::fill(fSc.begin(), fSc.end(), 0.);
67 std::fill(fSsx.begin(), fSsx.end(), 0.);
68 std::fill(fSsy.begin(), fSsy.end(), 0.);
69 std::fill(fScx.begin(), fScx.end(), 0.);
70 std::fill(fScy.begin(), fScy.end(), 0.);
71
72 //Reset radiation integral internal variables to defaults
73 fMeanPhotonAngleX =0.; //average angle of
74 //radiated photon direction in sampling, x-plane
75 fParamPhotonAngleX=1.e-3*rad; //a parameter of
76 //radiated photon sampling distribution, x-plane
77 fMeanPhotonAngleY =0.; //average angle of
78 //radiated photon direction in sampling, y-plane
79 fParamPhotonAngleY=1.e-3*rad; //a parameter of
80 //radiated photon sampling distribution, y-plane
81
82 fImin0 = 0;//set the first vector element to 0
83
84 //reset the trajectory
85 fParticleAnglesX.clear();
86 fParticleAnglesY.clear();
87 fScatteringAnglesX.clear();
88 fScatteringAnglesY.clear();
89 fSteps.clear();
90 fGlobalTimes.clear();
91 fParticleCoordinatesXYZ.clear();
92
93 //resets the vector of element numbers at the trajectory start
94 fImax0.clear();
95 //sets 0 element of the vector of element numbers
96 fImax0.push_back(0.);
97
98 //resets the radiation probability
99 fTotalRadiationProbabilityAlongTrajectory.clear();
100 //sets the radiation probability at the trajectory start
101 fTotalRadiationProbabilityAlongTrajectory.push_back(0.);
102}
103
104//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105
106void G4BaierKatkov::SetSpectrumEnergyRange(G4double emin,
107 G4double emax,
108 G4int numberOfBins)
109{
110 fMinPhotonEnergy = emin;
111 fMaxPhotonEnergy = emax;
112 fNBinsSpectrum = numberOfBins;
113
114 fLogEmaxdEmin = G4Log(fMaxPhotonEnergy/fMinPhotonEnergy);
115
116 //in initializing fNPhotonsPerBin
117 fNPhotonsPerBin.resize(fNBinsSpectrum);
118 std::fill(fNPhotonsPerBin.begin(), fNPhotonsPerBin.end(), 0);
119
120 //initializing the Spectrum
121 fSpectrum.resize(fNBinsSpectrum);
122 std::fill(fSpectrum.begin(), fSpectrum.end(), 0);
123
124 //initializing the fAccumTotalSpectrum
125 fAccumTotalSpectrum.resize(fNBinsSpectrum);
126 std::fill(fAccumTotalSpectrum.begin(), fAccumTotalSpectrum.end(), 0);
127
128 //initializing the fTotalSpectrum
129 fTotalSpectrum.resize(fNBinsSpectrum);
130 std::fill(fTotalSpectrum.begin(), fTotalSpectrum.end(), 0);
131
132 fPhotonEnergyInSpectrum.clear();
133 for (G4int j=0;j<fNBinsSpectrum;j++)
134 {
135 //bin position (mean between 2 bin limits)
136 fPhotonEnergyInSpectrum.push_back(fMinPhotonEnergy*
137 (std::exp(fLogEmaxdEmin*j/fNBinsSpectrum)+
138 std::exp(fLogEmaxdEmin*(j+1)/fNBinsSpectrum))/2.);
139 }
140
141 fItrajectories = 0;
142
143 ResetRadIntegral();
144}
145
146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
148void G4BaierKatkov::AddStatisticsInPhotonEnergyRegion(G4double emin,
149 G4double emax,
150 G4int timesPhotonStatistics)
151{
152
153 if(timesPhotonStatistics<=1)
154 {
155 G4cout << "G4BaierKatkov model, "
156 "function AddStatisticsInPhotonEnergyRegion("
157 << emin/CLHEP::MeV << " MeV, "
158 << emax/CLHEP::MeV << " MeV, "
159 << timesPhotonStatistics << ")" << G4endl;
160 G4cout << "Warning: the statistics factor cannot be <=1." << G4endl;
161 G4cout << "The statistics was not added." << G4endl;
162 G4cout << " "<< G4endl;
163 }
164 else if(fMinPhotonEnergy>emin)
165 {
166 G4cout << "G4BaierKatkov model, "
167 "function AddStatisticsInPhotonEnergyRegion("
168 << emin/CLHEP::MeV << " MeV, "
169 << emax/CLHEP::MeV << " MeV, "
170 << timesPhotonStatistics << ")" << G4endl;
171 G4cout << "Warning: the minimal energy inserted is less then "
172 "the minimal energy cut of the spectrum: "
173 << fMinPhotonEnergy/CLHEP::MeV << " MeV." << G4endl;
174 G4cout << "The statistics was not added." << G4endl;
175 G4cout << " "<< G4endl;
176 }
177 else if(emax-emin<DBL_EPSILON)
178 {
179 G4cout << "G4BaierKatkov model, "
180 "function AddStatisticsInPhotonEnergyRegion("
181 << emin/CLHEP::MeV << " MeV, "
182 << emax/CLHEP::MeV << " MeV, "
183 << timesPhotonStatistics << ")" << G4endl;
184 G4cout << "Warning: the maximal energy <= the minimal energy." << G4endl;
185 G4cout << "The statistics was not added." << G4endl;
186 G4cout << " "<< G4endl;
187 }
188 else
189 {
190 G4bool setrange = true;
191 G4double logAddRangeEmindEmin = G4Log(emin/fMinPhotonEnergy);
192 G4double logAddRangeEmaxdEmin = G4Log(emax/fMinPhotonEnergy);
193
194 G4int nAddRange = (G4int)fTimesPhotonStatistics.size();
195 for (G4int j=0;j<nAddRange;j++)
196 {
197 if((logAddRangeEmindEmin>=fLogAddRangeEmindEmin[j]&&
198 logAddRangeEmindEmin< fLogAddRangeEmaxdEmin[j])||
199 (logAddRangeEmaxdEmin> fLogAddRangeEmindEmin[j]&&
200 logAddRangeEmaxdEmin<=fLogAddRangeEmaxdEmin[j])||
201 (logAddRangeEmindEmin<=fLogAddRangeEmindEmin[j]&&
202 logAddRangeEmaxdEmin>=fLogAddRangeEmaxdEmin[j]))
203 {
204 G4cout << "G4BaierKatkov model, "
205 "function AddStatisticsInPhotonEnergyRegion("
206 << emin/CLHEP::MeV << " MeV, "
207 << emax/CLHEP::MeV << " MeV, "
208 << timesPhotonStatistics << ")" << G4endl;
209 G4cout << "Warning: the energy range intersects another "
210 "added energy range." << G4endl;
211 G4cout << "The statistics was not added." << G4endl;
212 G4cout << " "<< G4endl;
213 setrange = false;
214 break;
215 }
216 }
217 if (setrange)
218 {
219 fLogAddRangeEmindEmin.push_back(logAddRangeEmindEmin);
220 fLogAddRangeEmaxdEmin.push_back(logAddRangeEmaxdEmin);
221 fTimesPhotonStatistics.push_back(timesPhotonStatistics);
222
223 G4cout << "G4BaierKatkov model: increasing the statistics of photon sampling "
224 "in Baier-Katkov with a factor of "
225 << timesPhotonStatistics << G4endl;
226 G4cout << "in the energy spectrum range: ("
227 << emin/CLHEP::MeV << " MeV, "
228 << emax/CLHEP::MeV << " MeV)" << G4endl;
229 }
230 }
231}
232
233//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
234
235void G4BaierKatkov::SetPhotonSamplingParameters(G4double ekin,
236 G4double minPhotonAngleX,
237 G4double maxPhotonAngleX,
238 G4double minPhotonAngleY,
239 G4double maxPhotonAngleY)
240{
241 fLogEdEmin = G4Log(ekin/fMinPhotonEnergy);
242 fMeanPhotonAngleX = (maxPhotonAngleX+minPhotonAngleX)/2.;
243 fParamPhotonAngleX = (maxPhotonAngleX-minPhotonAngleX)/2.;
244 fMeanPhotonAngleY = (maxPhotonAngleY+minPhotonAngleY)/2.;
245 fParamPhotonAngleY = (maxPhotonAngleY-minPhotonAngleY)/2.;
246}
247
248//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249
250void G4BaierKatkov::GeneratePhotonSampling()
251{
252 fPhotonEnergyInIntegral.clear();
253 fPhotonAngleInIntegralX.clear();
254 fPhotonAngleInIntegralY.clear();
255 fPhotonAngleNormCoef.clear();
256 fInsideVirtualCollimator.clear();
257 fIBinsSpectrum.clear();
258
259 G4double ksi=0.;
260 std::vector<G4int> moreStatistics;
261 moreStatistics.resize(fTimesPhotonStatistics.size());
262 std::fill(moreStatistics.begin(), moreStatistics.end(), 0);
263 G4int nAddRange = (G4int)fTimesPhotonStatistics.size();
264
265 //sampling of the energy of a photon emission
266 //(integration variables, Monte Carlo integration)
267 for (G4int j=0;j<fNMCPhotons;j++)
268 {
269 ksi = G4UniformRand()*fLogEdEmin;
270 fIBinsSpectrum.push_back((G4int)std::trunc(
271 ksi*fNBinsSpectrum/fLogEmaxdEmin));
272 //we consider also the energy outside the spectrum output range
273 //(E>Emax => fLogEdEmin>fLogEmaxdEmin)
274 //in this case we don't count the photon in the spectrum output
275 if(fIBinsSpectrum[j]<fNBinsSpectrum) {fNPhotonsPerBin[fIBinsSpectrum[j]]+=1;}
276
277 fPhotonEnergyInIntegral.push_back(fMinPhotonEnergy*std::exp(ksi));
278
279 fPhotonAngleNormCoef.push_back(1.);
280
281 for (G4int j2=0;j2<nAddRange;j2++)
282 {
283 if(ksi>fLogAddRangeEmindEmin[j2]&&
284 ksi<fLogAddRangeEmaxdEmin[j2])
285 {
286 //calculating the current statistics in this region
287 //to increase it proportionally
288 moreStatistics[j2]+=1;
289 fPhotonAngleNormCoef[j]/=fTimesPhotonStatistics[j2];
290 break;
291 }
292 }
293 }
294
295 for (G4int j2=0;j2<nAddRange;j2++)
296 {
297 G4int totalAddRangeStatistics = moreStatistics[j2]*fTimesPhotonStatistics[j2];
298 for (G4int j=moreStatistics[j2];j<totalAddRangeStatistics;j++)
299 {
300 ksi = fLogAddRangeEmindEmin[j2]+
301 G4UniformRand()*(std::min(fLogAddRangeEmaxdEmin[j2],fLogEdEmin)-
302 fLogAddRangeEmindEmin[j2]);
303 fIBinsSpectrum.push_back((G4int)std::trunc(
304 ksi*fNBinsSpectrum/fLogEmaxdEmin));
305 /* //we consider also the energy outside the spectrum output range
306 //(E>Emax => fLogEdEmin>fLogEmaxdEmin)
307 //in this case we don't count the photon in the spectrum output
308 if(fIBinsSpectrum.back()<fNBinsSpectrum)
309 {fNPhotonsPerBin[fIBinsSpectrum.back()]+=1;}*/
310
311 fPhotonEnergyInIntegral.push_back(fMinPhotonEnergy*std::exp(ksi));
312
313 fPhotonAngleNormCoef.push_back(1./fTimesPhotonStatistics[j2]);
314 }
315 }
316
317 G4double rho=1.;
318 const G4double rhocut=15.;//radial angular cut of the distribution
319 G4double norm=std::atan(rhocut*rhocut)*
320 CLHEP::pi*fParamPhotonAngleX*fParamPhotonAngleY;
321
322 //sampling of the angles of a photon emission
323 //(integration variables, Monte Carlo integration)
324 G4int nmctotal = (G4int)fPhotonEnergyInIntegral.size();
325 for (G4int j=0;j<nmctotal;j++)
326 {
327 //photon distribution with long tails (useful to not exclude particle angles
328 //after a strong single scattering)
329 //at ellipsescale < 1 => half of statistics of photons
330 do
331 {
332 rho = std::sqrt(std::tan(CLHEP::halfpi*G4UniformRand()));
333 }
334 while (rho>rhocut);
335
336 ksi = G4UniformRand();
337 fPhotonAngleInIntegralX.push_back(fMeanPhotonAngleX+
338 fParamPhotonAngleX*
339 rho*std::cos(CLHEP::twopi*ksi));
340 fPhotonAngleInIntegralY.push_back(fMeanPhotonAngleY+
341 fParamPhotonAngleY*
342 rho*std::sin(CLHEP::twopi*ksi));
343 fPhotonAngleNormCoef[j]*=(1.+rho*rho*rho*rho)*norm;
344
345 //test if the photon with these angles enter the virtual collimator
346 //(doesn't influence the Geant4 simulations,
347 //but only the accumulation of fTotalSpectrum
348 fInsideVirtualCollimator.push_back(fVirtualCollimatorAngularDiameter >
349 std::sqrt(fPhotonAngleInIntegralX[j]*
350 fPhotonAngleInIntegralX[j]+
351 fPhotonAngleInIntegralY[j]*
352 fPhotonAngleInIntegralY[j]));
353 }
354 //reinitialize the vector of radiation CDF for each photon
355 fPhotonProductionCDF.resize(nmctotal+1);//0 element equal to 0
356 std::fill(fPhotonProductionCDF.begin(), fPhotonProductionCDF.end(), 0.);
357
358 //if we have additional photons
359 if (nmctotal>fNMCPhotons)
360 {
361 //reinitialize intermediate integrals with zeros again
362 fFa.resize(nmctotal);
363 fSs.resize(nmctotal);
364 fSc.resize(nmctotal);
365 fSsx.resize(nmctotal);
366 fSsy.resize(nmctotal);
367 fScx.resize(nmctotal);
368 fScy.resize(nmctotal);
369 std::fill(fFa.begin(), fFa.end(), 0.);
370 std::fill(fSs.begin(), fSs.end(), 0.);
371 std::fill(fSc.begin(), fSc.end(), 0.);
372 std::fill(fSsx.begin(), fSsx.end(), 0.);
373 std::fill(fSsy.begin(), fSsy.end(), 0.);
374 std::fill(fScx.begin(), fScx.end(), 0.);
375 std::fill(fScy.begin(), fScy.end(), 0.);
376 }
377}
378
379//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
380
381G4double G4BaierKatkov::RadIntegral(G4double etotal, G4double mass,
382 std::vector<G4double> &vectorParticleAnglesX,
383 std::vector<G4double> &vectorParticleAnglesY,
384 std::vector<G4double> &vectorScatteringAnglesX,
385 std::vector<G4double> &vectorScatteringAnglesY,
386 std::vector<G4double> &vectorSteps,
387 G4int imin)
388{
389 //preliminary values are defined:
390
391 G4double om=0.;// photon energy
392 G4double eprime=0., eprime2=0.; //E'=E-omega eprime2=eprime*eprime
393 G4double omprime=0.,omprimed2=0.;//om'=(E*om/E'), omprimed2=omprime/2
394 G4double vxin =0.,vyin=0.,vxno=0.,vyno=0.;
395
396 G4double dzmod=0.;
397 G4double fa1=0.,faseBefore=0.,faseBeforedz=0.,faseBeforedzd2=0.;
398 G4double faseAfter=0.,fa2dfaseBefore2=0.;
399
400 G4double skJ=0, skIx=0., skIy=0.;
401 G4double sinfa1=0.,cosfa1=0.;
402 G4double i2=0.,j2=0.;// Ix^2+Iy^2 and of BK Jvector^2
403
404 std::size_t nparts=vectorParticleAnglesX.size();
405 G4int kmin = imin;
406 if(imin==0) {kmin=1;}//skipping 0 trajectory element
407
408 //total radiation probability for each photon
409 G4double totalRadiationProbabilityPhj = 0.;
410
411 //total radiation probability along this trajectory (fill with 0 only new elements)
412 fTotalRadiationProbabilityAlongTrajectory.resize(nparts);
413
414 //reset Spectrum
415 std::fill(fSpectrum.begin(), fSpectrum.end(), 0.);
416
417 //intermediate vectors to reduce calculations
418 std::vector<G4double> axt;//acceleration of a charged particle in a horizontal plane
419 axt.resize(nparts);
420 std::vector<G4double> ayt;//acceleration of a charged particle in a vertical plane
421 ayt.resize(nparts);
422 std::vector<G4double> dz;//step in in MeV^-1
423 dz.resize(nparts);
424 //setting values interesting for us
425 for(std::size_t k=kmin;k<nparts;k++)
426 {
427 dz[k]=vectorSteps[k]/CLHEP::hbarc;// dz in MeV^-1
428
429 // accelerations
430 axt[k]=(vectorParticleAnglesX[k]-vectorScatteringAnglesX[k]-
431 vectorParticleAnglesX[k-1])/dz[k];
432 ayt[k]=(vectorParticleAnglesY[k]-vectorScatteringAnglesY[k]-
433 vectorParticleAnglesY[k-1])/dz[k];
434 }
435
436 //intermediate variables to reduce calculations:
437 //the calculations inside for (G4int j=0;j<fNMCPhotons;j++)
438 //{for(G4int k=kmin;k<nparts;k++){...
439 //are the main cpu time consumption
440 G4double e2 = etotal*etotal;
441 G4double gammaInverse2 = mass*mass/(etotal*etotal);// 1/gamma^2 of
442 //the radiating charge particle
443 G4double coefNormLogdNMC = fLogEdEmin/fNMCPhotons;//here fNMCPhotons is correct,
444 //additional photons have been already
445 //taken into account in weights
446 G4double coefNorm = CLHEP::fine_structure_const/(8*(CLHEP::pi2))*coefNormLogdNMC;
447 G4double e2pluseprime2 = 0.;//e2pluseprime2 =e2+eprime2
448 G4double coefNormom2deprime2 = 0.; //coefNormom2deprime2 = coefNorm*om2/eprime2;
449 G4double gammaInverse2om = 0.; //gammaInverse2*om
450
451 std::size_t nmctotal = fPhotonEnergyInIntegral.size();
452 for (std::size_t j=0;j<nmctotal;j++)
453 //the final number of photons may be different from fNMCPhotons
454 {
455 om = fPhotonEnergyInIntegral[j];
456 eprime=etotal-om; //E'=E-omega
457 eprime2 = eprime*eprime;
458 e2pluseprime2 =e2+eprime2;
459 omprime=etotal*om/eprime;//om'=(E*om/E')
460 omprimed2=omprime/2;
461 coefNormom2deprime2 = coefNorm*om*om/eprime2;
462 gammaInverse2om = gammaInverse2*om;
463
464 for(std::size_t k=kmin;k<nparts;k++)
465 {
466 //the angles of a photon produced (with incoherent scattering)
467 vxin = vectorParticleAnglesX[k]-fPhotonAngleInIntegralX[j];
468 vyin = vectorParticleAnglesY[k]-fPhotonAngleInIntegralY[j];
469 //the angles of a photon produced (without incoherent scattering)
470 vxno = vxin-vectorScatteringAnglesX[k];
471 vyno = vyin-vectorScatteringAnglesY[k];
472
473 //phase difference before scattering
474 faseBefore=omprimed2*(gammaInverse2+vxno*vxno+vyno*vyno);//phi' t<ti//MeV
475
476 faseBeforedz = faseBefore*dz[k];
477 faseBeforedzd2 = faseBeforedz/2.;
478 fFa[j]+=faseBeforedz; //
479 fa1=fFa[j]-faseBeforedzd2;//
480 dzmod=2*std::sin(faseBeforedzd2)/faseBefore;//MeV^-1
481
482 //phi''/faseBefore^2
483 fa2dfaseBefore2 = omprime*(axt[k]*vxno+ayt[k]*vyno)/(faseBefore*faseBefore);
484
485 //phase difference after scattering
486 faseAfter=omprimed2*(gammaInverse2+vxin*vxin+vyin*vyin);//phi' ti+O//MeV
487
488 skJ=1/faseAfter-1/faseBefore-fa2dfaseBefore2*dzmod;//MeV^-1
489 skIx=vxin/faseAfter-vxno/faseBefore+dzmod*(axt[k]/faseBefore-
490 vxno*fa2dfaseBefore2);
491 skIy=vyin/faseAfter-vyno/faseBefore+dzmod*(ayt[k]/faseBefore-
492 vyno*fa2dfaseBefore2);
493
494 sinfa1 = std::sin(fa1);
495 cosfa1 = std::cos(fa1);
496
497 fSs[j]+=sinfa1*skJ;//sum sin integral J of BK
498 fSc[j]+=cosfa1*skJ;//sum cos integral J of BK
499 fSsx[j]+=sinfa1*skIx;// sum sin integral Ix of BK
500 fSsy[j]+=sinfa1*skIy;// sum sin integral Iy of BK
501 fScx[j]+=cosfa1*skIx;// sum cos integral Ix of BK
502 fScy[j]+=cosfa1*skIy;// sum cos integral Iy of BK
503
504 i2=fSsx[j]*fSsx[j]+fScx[j]*fScx[j]+fSsy[j]*fSsy[j]+fScy[j]*fScy[j];//MeV^-2
505 j2=fSs[j]*fSs[j]+fSc[j]*fSc[j];//MeV^-2
506
507 //updating the total radiation probability along the trajectory
508 totalRadiationProbabilityPhj = coefNormom2deprime2*fPhotonAngleNormCoef[j]*
509 (i2*e2pluseprime2+j2*gammaInverse2om);
510 fTotalRadiationProbabilityAlongTrajectory[k] += totalRadiationProbabilityPhj;
511 }
512
513 fPhotonProductionCDF[j+1] = fTotalRadiationProbabilityAlongTrajectory.back();
514
515 //filling spectrum (adding photon probabilities to a histogram)
516 //we consider also the energy outside the spectrum output range
517 //(E>Emax => fLogEdEmin>fLogEmaxdEmin)
518 //in this case we don't count the photon in the spectrum output
519 //we fill the spectrum only in case of the angles inside the virtual collimator
520 if((fIBinsSpectrum[j]<fNBinsSpectrum)&&fInsideVirtualCollimator[j])
521 {fSpectrum[fIBinsSpectrum[j]] += totalRadiationProbabilityPhj/
522 (om*coefNormLogdNMC);}
523
524 } // end cycle
525
526 fAccumSpectrum.push_back(fSpectrum);
527
528 return fTotalRadiationProbabilityAlongTrajectory.back();
529}
530
531//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
532
533G4int G4BaierKatkov::FindVectorIndex(std::vector<G4double> &myvector, G4double value)
534{
535 auto iteratorbegin = myvector.begin();
536 auto iteratorend = myvector.end();
537
538 //vector index (for non precise values lower_bound gives upper value)
539 auto loweriterator = std::lower_bound(iteratorbegin, iteratorend, value);
540 //return the index of the vector element
541 return (G4int)std::distance(iteratorbegin, loweriterator);
542}
543
544//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
545
546G4bool G4BaierKatkov::SetPhotonProductionParameters(G4double etotal, G4double mass)
547{
548 //flag if a photon was produced
549 G4bool flagPhotonProduced = false;
550
551 //how many small pieces of a trajectory are
552 //before radiation (all if not radiation)
553 std::size_t nodeNumber = fAccumSpectrum.size()-1;
554
555 //ksi is a random number
556 //Generally ksi = G4UniformRand() is ok, but
557 //we use as a correction for the case
558 //when the radiation probability becomes too high (> 0.1)
559 G4double ksi = -G4Log(G4UniformRand());
560
561 if (ksi< fTotalRadiationProbabilityAlongTrajectory.back()) // photon produced
562 {
563 G4double ksi1 = G4UniformRand()*fTotalRadiationProbabilityAlongTrajectory.back();
564
565 //randomly choosing the photon to be produced from the sampling list
566 //according to the probabilities calculated in the Baier-Katkov integral
567 G4int iphoton = FindVectorIndex(fPhotonProductionCDF,ksi1)-1;//index of
568 //a photon produced
569
570 //energy of the photon produced
571 fEph0 = fPhotonEnergyInIntegral[iphoton];
572 //anagles of the photon produced
573 G4double photonAngleX = fPhotonAngleInIntegralX[iphoton];
574 G4double photonAngleY = fPhotonAngleInIntegralY[iphoton];
575
576 G4double momentumDirectionZ = 1./
577 std::sqrt(1.+std::pow(std::tan(photonAngleX),2)+
578 std::pow(std::tan(photonAngleY),2));
579
580 //momentum direction vector of the photon produced
581 PhMomentumDirection.set(momentumDirectionZ*std::tan(photonAngleX),
582 momentumDirectionZ*std::tan(photonAngleY),
583 momentumDirectionZ);
584
585 //random calculation of the radiation point index (iNode)
586 //ksi = G4UniformRand()*fTotalRadiationProbabilityAlongTrajectory.back();
587
588 //sort fTotalRadiationProbabilityAlongTrajectory
589 //(increasing but oscillating function => non-monotonic)
590 std::vector<G4double> temporaryVector;
591 temporaryVector.assign(fTotalRadiationProbabilityAlongTrajectory.begin(),
592 fTotalRadiationProbabilityAlongTrajectory.end());
593 std::sort(temporaryVector.begin(), temporaryVector.end());
594
595 //index of the point of radiation ("poststep") in sorted vector
596 G4int iNode = FindVectorIndex(temporaryVector,ksi);
597
598 //index of the point of radiation ("poststep") in unsorted vector
599 auto it = std::find_if(fTotalRadiationProbabilityAlongTrajectory.begin(),
600 fTotalRadiationProbabilityAlongTrajectory.end(),
601 [&](G4double value)
602 {return std::abs(value - temporaryVector[iNode]) < DBL_EPSILON;});
603 iNode = (G4int)std::distance(fTotalRadiationProbabilityAlongTrajectory.begin(), it);
604
605 //the piece of trajectory number (necessary only for the spectrum output)
606 nodeNumber = FindVectorIndex(fImax0,iNode*1.)-1;
607
608 //set new parameters of the charged particle
609 //(returning the particle back to the radiation point, i.e.
610 //remembering the new parameters for corresponding get functions)
611 fNewParticleEnergy = etotal-fEph0;
612 fNewParticleAngleX = fParticleAnglesX[iNode];
613 fNewParticleAngleY = fParticleAnglesY[iNode];
614 fNewGlobalTime = fGlobalTimes[iNode];
615 fNewParticleCoordinateXYZ = fParticleCoordinatesXYZ[iNode];
616
617 //particle angle correction from momentum conservation
618 //(important for fEph0 comparable to E,
619 // may kick off a particle from channeling)
620 G4double pratio = fEph0/std::sqrt(etotal*etotal-mass*mass);
621 fNewParticleAngleX -= std::asin(pratio*std::sin(photonAngleX));
622 fNewParticleAngleY -= std::asin(pratio*std::sin(photonAngleY));
623
624 flagPhotonProduced = true;
625 }
626
627 //accumulation during entire code run
628 for (G4int j=0;j<fNBinsSpectrum;j++)
629 {
630 fAccumTotalSpectrum[j] += fAccumSpectrum[nodeNumber][j];
631 //in the case of missing photon in spectrum, probability is 0
632 //(may happen only at the beginning of simulations)
633 if(fNPhotonsPerBin[j]==0)
634 {
635 fTotalSpectrum[j] = 0;
636 }
637 else
638 {
639 fTotalSpectrum[j] = fAccumTotalSpectrum[j]/fNPhotonsPerBin[j]*fItrajectories;
640 }
641 }
642
643 return flagPhotonProduced;
644}
645
646//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
647
648G4bool G4BaierKatkov::DoRadiation(G4double etotal, G4double mass,
649 G4double angleX, G4double angleY,
650 G4double angleScatteringX, G4double angleScatteringY,
651 G4double step, G4double globalTime,
652 G4ThreeVector coordinateXYZ,
653 G4bool flagEndTrajectory)
654{
655 /**
656 DoRadiation description
657
658 The real trajectory is divided onto parts.
659 Every part is accumulated until the radiation cannot be considered as single photon,
660 i.e. the radiation probability exceeds some threshold
661 fSinglePhotonRadiationProbabilityLimit.
662 Typically fSinglePhotonRadiationProbabilityLimit should be between 0.01 and 0.1.
663 Then the photon generation as well as its parameters are simulated
664 in SetPhotonProductionParameters()
665
666 Finally if radiation took place, this function sets the new particle parameters
667 (the parameters at the point of radiation emission)
668 fNewParticleEnergy; fNewParticleAngleX;
669 fNewParticleAngleY; NewStep; fNewGlobalTime; fNewParticleCoordinateXYZ;
670
671 In order to control if the trajectory part reaches this limit,
672 one needs to divide it into smaller pieces (consisting of
673 fNSmallTrajectorySteps elements, typically several hundred) and to calculate the total
674 radiation probability along the trajectory part after each small piece accomplished.
675 If it exceeds the fSinglePhotonRadiationProbabilityLimit,
676 the trajectory part is over and the radiation can be generated.
677
678 In order to calculate the radiation probability the Baier-Katkov integral is called
679 after each small piece of the trajectory. It is summarized with the integral along
680 the previous small pieces for this part of the trajectory.
681
682 To speed-up the simulations, the photon angles are generated only once
683 at the start of the trajectory part.
684 */
685
686 //flag if a photon was produced
687 G4bool flagPhotonProduced = false;
688
689 //adding the next trajectory element to the vectors
690 fParticleAnglesX.push_back(angleX);
691 fParticleAnglesY.push_back(angleY);
692 fScatteringAnglesX.push_back(angleScatteringX);
693 fScatteringAnglesY.push_back(angleScatteringY);
694 fSteps.push_back(step);
695 fGlobalTimes.push_back(globalTime);
696 fParticleCoordinatesXYZ.push_back(coordinateXYZ);
697
698 G4double imax = fSteps.size();
699 if((imax==fImin0+fNSmallTrajectorySteps)||flagEndTrajectory)
700 {
701 //set the angular limits at the start of the trajectory part
702 if(fImin0==0)
703 {
704 //radiation within the angle = +-fRadiationAngleFactor/gamma
705 G4double radiationAngleLimit=fRadiationAngleFactor*mass/etotal;
706
707 SetPhotonSamplingParameters(etotal-mass,
708 *std::min_element(fParticleAnglesX.begin(),
709 fParticleAnglesX.end())-radiationAngleLimit,
710 *std::max_element(fParticleAnglesX.begin(),
711 fParticleAnglesX.end())+radiationAngleLimit,
712 *std::min_element(fParticleAnglesY.begin(),
713 fParticleAnglesY.end())-radiationAngleLimit,
714 *std::max_element(fParticleAnglesY.begin(),
715 fParticleAnglesY.end())+radiationAngleLimit);
716
717 //calculation of angles of photon emission
718 //(these angles are integration variables, Monte Carlo integration)
719 GeneratePhotonSampling();
720 }
721
722 //calculate Baier-Katkov integral after this
723 //small piece of trajectory (fNSmallTrajectorySteps elements)
724 fTotalRadiationProbability = RadIntegral(etotal,mass,
725 fParticleAnglesX,fParticleAnglesY,
726 fScatteringAnglesX,fScatteringAnglesY,
727 fSteps, fImin0);
728
729 //setting the last element of this small trajectory piece
730 fImin0 = imax;
731 fImax0.push_back(imax*1.);
732
733 //if the radiation probability is high enough or if we are finishing
734 //our trajectory => to simulate the photon emission
735 if(fTotalRadiationProbability>fSinglePhotonRadiationProbabilityLimit||
736 flagEndTrajectory)
737 {
738 fItrajectories += 1; //count this trajectory !!!correction 19.07.2023
739
740 flagPhotonProduced = SetPhotonProductionParameters(etotal,mass);
741
742 // correction 19.07.2023 fItrajectories += 1; //count this trajectory
743
744 //reinitialize intermediate integrals fFa, fSs, fSc, fSsx, fSsy, fScx, fScy;
745 //reset radiation integral internal variables to defaults;
746 //reset the trajectory and radiation probability along the trajectory
747 ResetRadIntegral();
748 }
749 }
750
751 return flagPhotonProduced;
752}
753
754//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
755
756void G4BaierKatkov::GeneratePhoton(G4FastStep &fastStep)
757{
758 const G4DynamicParticle theGamma = G4DynamicParticle(G4Gamma::Gamma(),
759 PhMomentumDirection,fEph0);
760 //generation of a secondary photon
761 fastStep.CreateSecondaryTrack(theGamma,fNewParticleCoordinateXYZ,fNewGlobalTime,true);
762}
G4BaierKatkov.hh
Definition of the G4BaierKatkov class This class is designed for the calculation of radiation probabi...
G4Log
G4double G4Log(G4double x)
Definition G4Log.hh:227
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4BaierKatkov::ResetRadIntegral
void ResetRadIntegral()
Definition G4BaierKatkov.cc:52
G4BaierKatkov::GeneratePhoton
void GeneratePhoton(G4FastStep &fastStep)
Definition G4BaierKatkov.cc:756
G4BaierKatkov::SetSpectrumEnergyRange
void SetSpectrumEnergyRange(G4double emin, G4double emax, G4int numberOfBins)
Definition G4BaierKatkov.cc:106
G4BaierKatkov::AddStatisticsInPhotonEnergyRegion
void AddStatisticsInPhotonEnergyRegion(G4double emin, G4double emax, G4int timesPhotonStatistics)
Definition G4BaierKatkov.cc:148
G4BaierKatkov::DoRadiation
G4bool DoRadiation(G4double etotal, G4double mass, G4double angleX, G4double angleY, G4double angleScatteringX, G4double angleScatteringY, G4double step, G4double globalTime, G4ThreeVector coordinateXYZ, G4bool flagEndTrajectory=false)
Definition G4BaierKatkov.cc:648
G4FastStep::CreateSecondaryTrack
G4Track * CreateSecondaryTrack(const G4DynamicParticle &, G4ThreeVector, G4ThreeVector, G4double, G4bool localCoordinates=true)
Definition G4FastStep.cc:182
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition G4Gamma.cc:81
DBL_EPSILON
#define DBL_EPSILON
Definition templates.hh:66