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