| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Author: Alexei Sytov
27 // Co-author: Gianfranco Paterno (modificati
28 // On the base of the CRYSTALRAD realization o
29 // A. I. Sytov, V. V. Tikhomirov, and L. Bandi
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 o
42 //calls ResetRadIntegral()
43 SetSpectrumEnergyRange(0.1*MeV,1.*GeV,110)
44
45 //Do not worry if the maximal energy > par
46 //this elements of spectrum with non-physi
47 //will not be processed (they will be 0)
48
49 G4cout << " "<< G4endl;
50 G4cout << "G4BaierKatkov model is activate
51 G4cout << " "<< G4endl;
52 }
53
54 //....oooOO0OOooo........oooOO0OOooo........oo
55
56 void G4BaierKatkov::ResetRadIntegral()
57 {
58 fAccumSpectrum.clear();
59
60 //reinitialize intermediate integrals with
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 variab
77 fMeanPhotonAngleX =0.; //average an
78 //radiated p
79 fParamPhotonAngleX=1.e-3*rad; //a paramete
80 //radiated p
81 fMeanPhotonAngleY =0.; //average an
82 //radiated p
83 fParamPhotonAngleY=1.e-3*rad; //a paramete
84 //radiated p
85
86 fImin0 = 0;//set the first vector element
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
98 fImax0.clear();
99 //sets 0 element of the vector of element
100 fImax0.push_back(0.);
101
102 //resets the radiation probability
103 fTotalRadiationProbabilityAlongTrajectory.
104 //sets the radiation probability at the tr
105 fTotalRadiationProbabilityAlongTrajectory.
106 }
107
108 //....oooOO0OOooo........oooOO0OOooo........oo
109
110 void G4BaierKatkov::SetSpectrumEnergyRange(G4d
111 G4d
112 G4i
113 {
114 fMinPhotonEnergy = emin;
115 fMaxPhotonEnergy = emax;
116 fNBinsSpectrum = numberOfBins;
117
118 fLogEmaxdEmin = G4Log(fMaxPhotonEnergy/fMi
119
120 //in initializing fNPhotonsPerBin
121 fNPhotonsPerBin.resize(fNBinsSpectrum);
122 std::fill(fNPhotonsPerBin.begin(), fNPhoto
123
124 //initializing the Spectrum
125 fSpectrum.resize(fNBinsSpectrum);
126 std::fill(fSpectrum.begin(), fSpectrum.end
127
128 //initializing the fAccumTotalSpectrum
129 fAccumTotalSpectrum.resize(fNBinsSpectrum)
130 std::fill(fAccumTotalSpectrum.begin(), fAc
131
132 //initializing the fTotalSpectrum
133 fTotalSpectrum.resize(fNBinsSpectrum);
134 std::fill(fTotalSpectrum.begin(), fTotalSp
135
136 fPhotonEnergyInSpectrum.clear();
137 for (G4int j=0;j<fNBinsSpectrum;j++)
138 {
139 //bin position (mean between 2 bin lim
140 fPhotonEnergyInSpectrum.push_back(fMin
141 (std::e
142 std::e
143 }
144
145 fItrajectories = 0;
146
147 ResetRadIntegral();
148 }
149
150 //....oooOO0OOooo........oooOO0OOooo........oo
151
152 void G4BaierKatkov::AddStatisticsInPhotonEnerg
153
154
155 {
156
157 if(timesPhotonStatistics<=1)
158 {
159 G4cout << "G4BaierKatkov model, "
160 "function AddStatisticsInPho
161 << emin/CLHEP::M
162 << emax/CLHEP::M
163 << timesPhotonSt
164 G4cout << "Warning: the statistics fac
165 G4cout << "The statistics was not adde
166 G4cout << " "<< G4endl;
167 }
168 else if(fMinPhotonEnergy>emin)
169 {
170 G4cout << "G4BaierKatkov model, "
171 "function AddStatisticsInPho
172 << emin/CLHEP::M
173 << emax/CLHEP::M
174 << timesPhotonSt
175 G4cout << "Warning: the minimal energy
176 "the minimal energy cut of t
177 << fMinPhotonEnergy/CLHEP::MeV
178 G4cout << "The statistics was not adde
179 G4cout << " "<< G4endl;
180 }
181 else if(emax-emin<DBL_EPSILON)
182 {
183 G4cout << "G4BaierKatkov model, "
184 "function AddStatisticsInPho
185 << emin/CLHEP::M
186 << emax/CLHEP::M
187 << timesPhotonSt
188 G4cout << "Warning: the maximal energy
189 G4cout << "The statistics was not adde
190 G4cout << " "<< G4endl;
191 }
192 else
193 {
194 G4bool setrange = true;
195 G4double logAddRangeEmindEmin = G4Log(
196 G4double logAddRangeEmaxdEmin = G4Log(
197
198 G4int nAddRange = (G4int)fTimesPhotonS
199 for (G4int j=0;j<nAddRange;j++)
200 {
201 if((logAddRangeEmindEmin>=fLogAddR
202 logAddRangeEmindEmin< fLogAddR
203 (logAddRangeEmaxdEmin> fLogAddR
204 logAddRangeEmaxdEmin<=fLogAddR
205 (logAddRangeEmindEmin<=fLogAddR
206 logAddRangeEmaxdEmin>=fLogAddR
207 {
208 G4cout << "G4BaierKatkov model
209 "function AddStatist
210 << emin/
211 << emax/
212 << times
213 G4cout << "Warning: the energy
214 "added energy range."
215 G4cout << "The statistics was n
216 G4cout << " "<< G4endl;
217 setrange = false;
218 break;
219 }
220 }
221 if (setrange)
222 {
223 fLogAddRangeEmindEmin.push_back(lo
224 fLogAddRangeEmaxdEmin.push_back(lo
225 fTimesPhotonStatistics.push_back(t
226
227 G4cout << "G4BaierKatkov model: in
228 "in Baier-Katkov with a
229 << timesPhotonStatistics <<
230 G4cout << "in the energy spectrum
231 << emin/CLHEP::MeV << " MeV
232 << emax/CLHEP::MeV << " MeV
233 }
234 }
235 }
236
237 //....oooOO0OOooo........oooOO0OOooo........oo
238
239 void G4BaierKatkov::SetPhotonSamplingParameter
240
241
242
243
244 {
245 fLogEdEmin = G4Log(ekin/fMinPhotonEnergy);
246 fMeanPhotonAngleX = (maxPhotonAngleX+minP
247 fParamPhotonAngleX = (maxPhotonAngleX-minP
248 fMeanPhotonAngleY = (maxPhotonAngleY+minP
249 fParamPhotonAngleY = (maxPhotonAngleY-minP
250 }
251
252 //....oooOO0OOooo........oooOO0OOooo........oo
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(fTimesPhotonStatisti
266 std::fill(moreStatistics.begin(), moreStat
267 G4int nAddRange = (G4int)fTimesPhotonStati
268
269 //sampling of the energy of a photon emiss
270 //(integration variables, Monte Carlo inte
271 for (G4int j=0;j<fNMCPhotons;j++)
272 {
273 ksi = G4UniformRand()*fLogEdEmin;
274 fIBinsSpectrum.push_back((G4int)std::t
275 ksi*fNBin
276 //we consider also the energy outside
277 //(E>Emax => fLogEdEmin>fLogEmaxdEmin)
278 //in this case we don't count the phot
279 if(fIBinsSpectrum[j]<fNBinsSpectrum) {
280
281 fPhotonEnergyInIntegral.push_back(fMin
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 stati
291 //to increase it proportionally
292 moreStatistics[j2]+=1;
293 fPhotonAngleNormCoef[j]/=fTimes
294 break;
295 }
296 }
297 }
298
299 for (G4int j2=0;j2<nAddRange;j2++)
300 {
301 G4int totalAddRangeStatistics = moreSt
302 for (G4int j=moreStatistics[j2];j<tota
303 {
304 ksi = fLogAddRangeEmindEmin[j2]+
305 G4UniformRand()*(std::min(fL
306 fLogAddRang
307 fIBinsSpectrum.push_back((G4int)st
308 ksi*f
309 /* //we consider also the energy ou
310 //(E>Emax => fLogEdEmin>fLogEmaxdE
311 //in this case we don't count the
312 if(fIBinsSpectrum.back()<fNBinsSpe
313 {fNPhotonsPerBin[fIBinsSpectru
314
315 fPhotonEnergyInIntegral.push_back(
316
317 fPhotonAngleNormCoef.push_back(1./
318 }
319 }
320
321 G4double rho=1.;
322 const G4double rhocut=15.;//radial angular
323 G4double norm=std::atan(rhocut*rhocut)*
324 CLHEP::pi*fParamPh
325
326 //sampling of the angles of a photon emiss
327 //(integration variables, Monte Carlo inte
328 G4int nmctotal = (G4int)fPhotonEnergyInInt
329 for (G4int j=0;j<nmctotal;j++)
330 {
331 //photon distribution with long tails
332 //after a strong single scattering)
333 //at ellipsescale < 1 => half of stati
334 do
335 {
336 rho = std::sqrt(std::tan(CLHEP::ha
337 }
338 while (rho>rhocut);
339
340 ksi = G4UniformRand();
341 fPhotonAngleInIntegralX.push_back(fMea
342 fPara
343 rho*
344 fPhotonAngleInIntegralY.push_back(fMea
345 fPara
346 rho*
347 fPhotonAngleNormCoef[j]*=(1.+rho*rho*r
348
349 //test if the photon with these angles
350 //(doesn't influence the Geant4 simula
351 //but only the accumulation of fTotalS
352 fInsideVirtualCollimator.push_back(fVi
353 std
354
355
356
357 }
358 //reinitialize the vector of radiation CDF
359 fPhotonProductionCDF.resize(nmctotal+1);//
360 std::fill(fPhotonProductionCDF.begin(), fP
361
362 //if we have additional photons
363 if (nmctotal>fNMCPhotons)
364 {
365 //reinitialize intermediate integrals
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........oo
384
385 G4double G4BaierKatkov::RadIntegral(G4double e
386 std::vecto
387 std::vecto
388 std::vecto
389 std::vecto
390 std::vecto
391 G4int imin
392 {
393 //preliminary values are defined:
394
395 G4double om=0.;// photon energy
396 G4double eprime=0., eprime2=0.; //E'=E-ome
397 G4double omprime=0.,omprimed2=0.;//om'=(E*
398 G4double vxin =0.,vyin=0.,vxno=0.,vyno=0.;
399
400 G4double dzmod=0.;
401 G4double fa1=0.,faseBefore=0.,faseBeforedz
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 B
407
408 std::size_t nparts=vectorParticleAnglesX.s
409 G4int kmin = imin;
410 if(imin==0) {kmin=1;}//skipping 0 trajecto
411
412 //total radiation probability for each pho
413 G4double totalRadiationProbabilityPhj = 0.
414
415 //total radiation probability along this t
416 fTotalRadiationProbabilityAlongTrajectory.
417
418 //reset Spectrum
419 std::fill(fSpectrum.begin(), fSpectrum.end
420
421 //intermediate vectors to reduce calculati
422 std::vector<G4double> axt;//acceleration o
423 axt.resize(nparts);
424 std::vector<G4double> ayt;//acceleration o
425 ayt.resize(nparts);
426 std::vector<G4double> dz;//step in in MeV^
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;// d
432
433 // accelerations
434 axt[k]=(vectorParticleAnglesX[k]-vecto
435 vectorParticleAnglesX[k-1])/dz
436 ayt[k]=(vectorParticleAnglesY[k]-vecto
437 vectorParticleAnglesY[k-1])/dz
438 }
439
440 //intermediate variables to reduce calcula
441 //the calculations inside for (G4int j=0;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
446
447 G4double coefNormLogdNMC = fLogEdEmin/fNMC
448 //ad
449 //ta
450 G4double coefNorm = CLHEP::fine_structure_
451 G4double e2pluseprime2 = 0.;//e2pluseprime
452 G4double coefNormom2deprime2 = 0.; //coefN
453 G4double gammaInverse2om2 = 0.; //gammaInv
454
455 std::size_t nmctotal = fPhotonEnergyInInte
456 for (std::size_t j=0;j<nmctotal;j++)
457 //the final number of photons may be d
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/e
466 gammaInverse2om2 = gammaInverse2*om*om
467
468 for(std::size_t k=kmin;k<nparts;k++)
469 {
470 //the angles of a photon produced
471 vxin = vectorParticleAnglesX[k]-fP
472 vyin = vectorParticleAnglesY[k]-fP
473 //the angles of a photon produced
474 vxno = vxin-vectorScatteringAngles
475 vyno = vyin-vectorScatteringAngles
476
477 //phase difference before scatteri
478 faseBefore=omprimed2*(gammaInverse
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)/f
485
486 //phi''/faseBefore^2
487 fa2dfaseBefore2 = omprime*(axt[k]*
488
489 //phase difference after scatterin
490 faseAfter=omprimed2*(gammaInverse2
491
492 skJ=1/faseAfter-1/faseBefore-fa2df
493 skIx=vxin/faseAfter-vxno/faseBefor
494
495 skIy=vyin/faseAfter-vyno/faseBefor
496
497
498 sinfa1 = std::sin(fa1);
499 cosfa1 = std::cos(fa1);
500
501 fSs[j]+=sinfa1*skJ;//sum sin integ
502 fSc[j]+=cosfa1*skJ;//sum cos integ
503 fSsx[j]+=sinfa1*skIx;// sum sin in
504 fSsy[j]+=sinfa1*skIy;// sum sin in
505 fScx[j]+=cosfa1*skIx;// sum cos in
506 fScy[j]+=cosfa1*skIy;// sum cos in
507
508 i2=fSsx[j]*fSsx[j]+fScx[j]*fScx[j]
509 j2=fSs[j]*fSs[j]+fSc[j]*fSc[j];//M
510
511 //updating the total radiation pro
512 totalRadiationProbabilityPhj = coe
513 (i2*e2pluseprime2+j2*gamma
514 fTotalRadiationProbabilityAlongTra
515 }
516
517 fPhotonProductionCDF[j+1] = fTotalRadi
518
519 //filling spectrum (adding photon prob
520 //we consider also the energy outside
521 //(E>Emax => fLogEdEmin>fLogEmaxdEmin)
522 //in this case we don't count the phot
523 //we fill the spectrum only in case of
524 if((fIBinsSpectrum[j]<fNBinsSpectrum)&
525 {fSpectrum[fIBinsSpectrum[j]] += t
526 (om*coefNormLogdNMC);}
527
528 } // end cycle
529
530 fAccumSpectrum.push_back(fSpectrum);
531
532 return fTotalRadiationProbabilityAlongTraj
533 }
534
535 //....oooOO0OOooo........oooOO0OOooo........oo
536
537 G4int G4BaierKatkov::FindVectorIndex(std::vect
538 {
539 auto iteratorbegin = myvector.begin();
540 auto iteratorend = myvector.end();
541
542 //vector index (for non precise values low
543 auto loweriterator = std::lower_bound(iter
544 //return the index of the vector element
545 return (G4int)std::distance(iteratorbegin,
546 }
547
548 //....oooOO0OOooo........oooOO0OOooo........oo
549
550 G4bool G4BaierKatkov::SetPhotonProductionParam
551 {
552 //flag if a photon was produced
553 G4bool flagPhotonProduced = false;
554
555 //how many small pieces of a trajectory ar
556 //before radiation (all if not radiation)
557 std::size_t nodeNumber = fAccumSpectrum.si
558
559 //ksi is a random number
560 //Generally ksi = G4UniformRand() is ok, b
561 //we use as a correction for the case
562 //when the radiation probability becomes t
563 G4double ksi = -G4Log(G4UniformRand());
564
565 if (ksi< fTotalRadiationProbabilityAlongTr
566 {
567 G4double ksi1 = G4UniformRand()*fTotalRa
568
569 //randomly choosing the photon to be pro
570 //according to the probabilities calcula
571 G4int iphoton = FindVectorIndex(fPhotonP
572
573
574 //energy of the photon produced
575 fEph0 = fPhotonEnergyInIntegral[iphoton]
576 //anagles of the photon produced
577 G4double photonAngleX = fPhotonAngleInIn
578 G4double photonAngleY = fPhotonAngleInIn
579
580 G4double momentumDirectionZ = 1./
581 std::sqrt(1.+std::pow(std::tan(photo
582 std::pow(std::tan(photo
583
584 //momentum direction vector of the photo
585 PhMomentumDirection.set(momentumDirectio
586 momentumDirectio
587 momentumDirectio
588
589 //random calculation of the radiation po
590 //ksi = G4UniformRand()*fTotalRadiationP
591
592 //sort fTotalRadiationProbabilityAlongTr
593 //(increasing but oscillating function =
594 std::vector<G4double> temporaryVector;
595 temporaryVector.assign(fTotalRadiationPr
596 fTotalRadiationPr
597 std::sort(temporaryVector.begin(), tempo
598
599 //index of the point of radiation ("post
600 G4int iNode = FindVectorIndex(temporaryV
601
602 //index of the point of radiation ("post
603 auto it = std::find_if(fTotalRadiationPr
604 fTotalRadiationPr
605 [&](G4double valu
606 {return std::abs(value
607 iNode = (G4int)std::distance(fTotalRadia
608
609 //the piece of trajectory number (necess
610 nodeNumber = FindVectorIndex(fImax0,iNod
611
612 //set new parameters of the charged part
613 //(returning the particle back to the ra
614 //remembering the new parameters for cor
615 fNewParticleEnergy = etotal-fEph0;
616 fNewParticleAngleX = fParticleAnglesX[iN
617 fNewParticleAngleY = fParticleAnglesY[iN
618 fNewGlobalTime = fGlobalTimes[iNode];
619 fNewParticleCoordinateXYZ = fParticleCoo
620
621 //particle angle correction from momentu
622 //(important for fEph0 comparable to E,
623 // may kick off a particle from channeli
624 G4double pratio = fEph0/std::sqrt(etotal
625 fNewParticleAngleX -= std::asin(pratio*s
626 fNewParticleAngleY -= std::asin(pratio*s
627
628 flagPhotonProduced = true;
629 }
630
631 //accumulation during entire code run
632 for (G4int j=0;j<fNBinsSpectrum;j++)
633 {
634 fAccumTotalSpectrum[j] += fAccumSpectr
635 //in the case of missing photon in spe
636 //(may happen only at the beginning of
637 if(fNPhotonsPerBin[j]==0)
638 {
639 fTotalSpectrum[j] = 0;
640 }
641 else
642 {
643 fTotalSpectrum[j] = fAccumTotalSpec
644 }
645 }
646
647 return flagPhotonProduced;
648 }
649
650 //....oooOO0OOooo........oooOO0OOooo........oo
651
652 G4bool G4BaierKatkov::DoRadiation(G4double eto
653 G4double ang
654 G4double ang
655 G4double ste
656 G4ThreeVecto
657 G4bool flagE
658 {
659 /**
660 DoRadiation description
661
662 The real trajectory is divided onto parts.
663 Every part is accumulated until the radiat
664 i.e. the radiation probability exceeds som
665 fSinglePhotonRadiationProbabilityLimit.
666 Typically fSinglePhotonRadiationProbabilit
667 Then the photon generation as well as its
668 in SetPhotonProductionParameters()
669
670 Finally if radiation took place, this func
671 (the parameters at the point of radiation
672 fNewParticleEnergy; fNewParticleAngleX;
673 fNewParticleAngleY; NewStep; fNewGlobalTim
674
675 In order to control if the trajectory part
676 one needs to divide it into smaller pieces
677 fNSmallTrajectorySteps elements, typically
678 radiation probability along the trajectory
679 If it exceeds the fSinglePhotonRadiationPr
680 the trajectory part is over and the radiat
681
682 In order to calculate the radiation probab
683 after each small piece of the trajectory.
684 the previous small pieces for this part of
685
686 To speed-up the simulations, the photon an
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 th
694 fParticleAnglesX.push_back(angleX);
695 fParticleAnglesY.push_back(angleY);
696 fScatteringAnglesX.push_back(angleScatteri
697 fScatteringAnglesY.push_back(angleScatteri
698 fSteps.push_back(step);
699 fGlobalTimes.push_back(globalTime);
700 fParticleCoordinatesXYZ.push_back(coordina
701
702 G4double imax = fSteps.size();
703 if((imax==fImin0+fNSmallTrajectorySteps)||
704 {
705 //set the angular limits at the start
706 if(fImin0==0)
707 {
708 //radiation within the angle = +-f
709 G4double radiationAngleLimit=fRadi
710
711 SetPhotonSamplingParameters(etotal
712 *std::min_element(fParticleAnglesX
713 fParticleAnglesX
714 *std::max_element(fParticleAnglesX
715 fParticleAnglesX
716 *std::min_element(fParticleAnglesY
717 fParticleAnglesY
718 *std::max_element(fParticleAnglesY
719 fParticleAnglesY
720
721 //calculation of angles of photon
722 //(these angles are integration va
723 GeneratePhotonSampling();
724 }
725
726 //calculate Baier-Katkov integral afte
727 //small piece of trajectory (fNSmallTr
728 fTotalRadiationProbability = RadIntegr
729
730
731
732
733 //setting the last element of this sma
734 fImin0 = imax;
735 fImax0.push_back(imax*1.);
736
737 //if the radiation probability is high
738 //our trajectory => to simulate the ph
739 if(fTotalRadiationProbability>fSingleP
740 flagEndTrajectory)
741 {
742 fItrajectories += 1; //count this
743
744 flagPhotonProduced = SetPhotonProd
745
746 // correction 19.07.2023 fItraject
747
748 //reinitialize intermediate integr
749 //reset radiation integral interna
750 //reset the trajectory and radiati
751 ResetRadIntegral();
752 }
753 }
754
755 return flagPhotonProduced;
756 }
757
758 //....oooOO0OOooo........oooOO0OOooo........oo
759
760 void G4BaierKatkov::GeneratePhoton(G4FastStep
761 {
762 const G4DynamicParticle theGamma = G4Dynam
763
764 //generation of a secondary photon
765 fastStep.CreateSecondaryTrack(theGamma,fNe
766 }
767