60 const G4int nPerDecade = 10;
66 fLowestKineticEnergy = std::max(tmin, lowestTkin);
67 fHighestKineticEnergy = tmax;
69 if(tmax < 10*fLowestKineticEnergy)
71 fHighestKineticEnergy = 10*fLowestKineticEnergy;
73 else if(tmax > highestTkin)
75 fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
77 fTotBin = (
G4int)(nPerDecade*
78 std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
81 fHighestKineticEnergy,
84 G4cout <<
"### G4PAIPhotData: Nbins= " << fTotBin
85 <<
" Tmin(MeV)= " << fLowestKineticEnergy/MeV
86 <<
" Tmax(GeV)= " << fHighestKineticEnergy/GeV
87 <<
" tmin(keV)= " << tmin/keV <<
G4endl;
96 std::size_t n = fPAIxscBank.size();
99 for(std::size_t i=0; i<n; ++i)
103 fPAIxscBank[i]->clearAndDestroy();
104 delete fPAIxscBank[i];
105 fPAIxscBank[i] =
nullptr;
109 fPAIdEdxBank[i]->clearAndDestroy();
110 delete fPAIdEdxBank[i];
111 fPAIdEdxBank[i] =
nullptr;
113 delete fdEdxTable[i];
114 delete fdNdxCutTable[i];
115 fdEdxTable[i] =
nullptr;
116 fdNdxCutTable[i] =
nullptr;
119 delete fParticleEnergyVector;
120 fParticleEnergyVector =
nullptr;
134 for (jMatCC = 0; jMatCC < numOfCouples; ++jMatCC )
138 if( jMatCC == numOfCouples && jMatCC > 0 ) --jMatCC;
142 G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
143 G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
145 G4cout<<
"G4PAIPhotData::Initialise: "<<
"cut = "<<cut/keV<<
" keV; cutEl = "
146 <<deltaCutInKineticEnergyNow/keV<<
" keV; cutPh = "
147 <<photonCutInKineticEnergyNow/keV<<
" keV"<<
G4endl;
153 fHighestKineticEnergy,
158 fHighestKineticEnergy,
160 auto dNdxCutPhotonVector =
162 fHighestKineticEnergy,
164 auto dNdxCutPlasmonVector =
166 fHighestKineticEnergy,
177 auto dEdxMeanVector =
179 fHighestKineticEnergy,
188 for (
G4int i = 0; i <= fTotBin; ++i)
192 G4double tau = kinEnergy/proton_mass_c2;
195 if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow;
197 fPAIxSection.
Initialize( mat, Tmax, bg2, &fSandia);
210 for(
G4int k = 0; k < n; k++ )
214 transferVector->PutValue(k , t,
216 photonVector->PutValue(k , t,
218 plasmonVector->PutValue(k , t,
227 if(ionloss < 0.0) ionloss = 0.0;
229 dEdxMeanVector->PutValue(i,ionloss);
231 G4double dNdxCut = transferVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
232 G4double dNdxCutPhoton = photonVector->Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
233 G4double dNdxCutPlasmon = plasmonVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
235 G4double dEdxCut = dEdxVector->Value(cut)/cut;
238 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
239 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
240 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
242 dNdxCutVector->PutValue(i, dNdxCut);
243 dNdxCutPhotonVector->PutValue(i, dNdxCutPhoton);
244 dNdxCutPlasmonVector->PutValue(i, dNdxCutPlasmon);
246 dEdxCutVector->PutValue(i, dEdxCut);
248 PAItransferTable->insertAt(i,transferVector);
249 PAIphotonTable->insertAt(i,photonVector);
250 PAIplasmonTable->insertAt(i,plasmonVector);
251 PAIdEdxTable->insertAt(i,dEdxVector);
255 fPAIxscBank.push_back(PAItransferTable);
256 fPAIphotonBank.push_back(PAIphotonTable);
257 fPAIplasmonBank.push_back(PAIplasmonTable);
259 fPAIdEdxBank.push_back(PAIdEdxTable);
260 fdEdxTable.push_back(dEdxMeanVector);
262 fdNdxCutTable.push_back(dNdxCutVector);
263 fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
264 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
266 fdEdxCutTable.push_back(dEdxCutVector);
276 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
280 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
281 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
286 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
287 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
289 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
300 if( dEdx < 0.) { dEdx = 0.; }
311 G4double cross, xscEl, xscEl2, xscPh, xscPh2;
317 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
322 if( scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
323 else if( scaledTkin > fParticleEnergyVector->
Energy(0) ) one =
false;
326 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
327 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
332 cross = xscPh + xscEl;
336 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
348 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
350 E1 = fParticleEnergyVector->
Energy(iPlace);
351 E2 = fParticleEnergyVector->
Energy(iPlace+1);
354 W1 = (E2 - scaledTkin)*
W;
355 W2 = (scaledTkin - E1)*
W;
360 cross = xscEl + xscPh;
362 if( cross < 0.0) cross = 0.0;
372 G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
375 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
380 if( scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
381 else if( scaledTkin > fParticleEnergyVector->
Energy(0) ) one =
false;
384 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
385 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
390 cross = xscPh + xscEl;
394 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
406 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
408 E1 = fParticleEnergyVector->
Energy(iPlace);
409 E2 = fParticleEnergyVector->
Energy(iPlace+1);
412 W1 = (E2 - scaledTkin)*
W;
413 W2 = (scaledTkin - E1)*
W;
418 cross = xscEl + xscPh;
426 plRatio = xscEl/cross;
428 if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
444 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
449 if (scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
450 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) one =
false;
456 dNdxCut1 = (*vcut)[iPlace];
460 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
472 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
473 dNdxCut2 = (*vcut)[iPlace+1];
476 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
478 E1 = fParticleEnergyVector->
Energy(iPlace);
479 E2 = fParticleEnergyVector->
Energy(iPlace+1);
481 W1 = (E2 - scaledTkin)*
W;
482 W2 = (scaledTkin - E1)*
W;
483 meanNumber = W1*meanN1 + W2*meanN2;
488 if( meanNumber <= 0.0)
return 0.0;
494 if( 0 == numOfCollisions)
return 0.0;
496 for(
G4int i=0; i< numOfCollisions; ++i)
499 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
500 omega = GetEnergyTransfer(coupleIndex, iPlace,
position);
506 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
508 omega = omega*W1 + omega2*W2;
513 if( loss > kinEnergy) {
break; }
519 if ( loss > kinEnergy) loss = kinEnergy;
520 else if( loss < 0.) loss = 0.;
536 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
541 if (scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
542 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) one =
false;
548 dNdxCut1 = (*vcut)[iPlace];
552 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
564 v2 = (*(fPAIphotonBank[coupleIndex]))(iPlace+1);
565 dNdxCut2 = (*vcut)[iPlace+1];
568 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
570 E1 = fParticleEnergyVector->
Energy(iPlace);
571 E2 = fParticleEnergyVector->
Energy(iPlace+1);
573 W1 = (E2 - scaledTkin)*
W;
574 W2 = (scaledTkin - E1)*
W;
575 meanNumber = W1*meanN1 + W2*meanN2;
580 if( meanNumber <= 0.0)
return 0.0;
586 if( 0 == numOfCollisions)
return 0.0;
588 for(
G4int i=0; i< numOfCollisions; ++i)
591 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
592 omega = GetEnergyPhotonTransfer(coupleIndex, iPlace,
position);
598 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
600 omega = omega*W1 + omega2*W2;
605 if( loss > kinEnergy) {
break; }
611 if ( loss > kinEnergy) loss = kinEnergy;
612 else if( loss < 0.) loss = 0.;
628 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
633 if (scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
634 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) one =
false;
640 dNdxCut1 = (*vcut)[iPlace];
644 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
656 v2 = (*(fPAIplasmonBank[coupleIndex]))(iPlace+1);
657 dNdxCut2 = (*vcut)[iPlace+1];
660 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
662 E1 = fParticleEnergyVector->
Energy(iPlace);
663 E2 = fParticleEnergyVector->
Energy(iPlace+1);
665 W1 = (E2 - scaledTkin)*
W;
666 W2 = (scaledTkin - E1)*
W;
667 meanNumber = W1*meanN1 + W2*meanN2;
672 if( meanNumber <= 0.0)
return 0.0;
678 if( 0 == numOfCollisions)
return 0.0;
680 for(
G4int i=0; i< numOfCollisions; ++i)
683 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
684 omega = GetEnergyPlasmonTransfer(coupleIndex, iPlace,
position);
690 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
691 G4double omega2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1,
position);
692 omega = omega*W1 + omega2*W2;
697 if( loss > kinEnergy) {
break; }
703 if ( loss > kinEnergy) loss = kinEnergy;
704 else if( loss < 0.) loss = 0.;
729 if( scaledTkin >= fParticleEnergyVector->
GetMaxEnergy())
732 transfer = GetEnergyTransfer(coupleIndex, nPlace,
position);
734 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
737 transfer = GetEnergyTransfer(coupleIndex, 0,
position);
741 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
743 dNdxCut1 = (*cutv)[iPlace];
744 dNdxCut2 = (*cutv)[iPlace+1];
746 E1 = fParticleEnergyVector->
Energy(iPlace);
747 E2 = fParticleEnergyVector->
Energy(iPlace+1);
749 W1 = (E2 - scaledTkin)*
W;
750 W2 = (scaledTkin - E1)*
W;
761 transfer = tr1*W1 + tr2*W2;
764 if(transfer < 0.0 ) { transfer = 0.0; }
786 if( scaledTkin >= fParticleEnergyVector->
GetMaxEnergy())
789 transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace,
position);
791 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
794 transfer = GetEnergyPhotonTransfer(coupleIndex, 0,
position);
798 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
800 dNdxCut1 = (*cutv)[iPlace];
801 dNdxCut2 = (*cutv)[iPlace+1];
803 E1 = fParticleEnergyVector->
Energy(iPlace);
804 E2 = fParticleEnergyVector->
Energy(iPlace+1);
806 W1 = (E2 - scaledTkin)*
W;
807 W2 = (scaledTkin - E1)*
W;
819 transfer = tr1*W1 + tr2*W2;
822 if(transfer < 0.0 ) { transfer = 0.0; }
843 if( scaledTkin >= fParticleEnergyVector->
GetMaxEnergy())
846 transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace,
position);
848 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
851 transfer = GetEnergyPlasmonTransfer(coupleIndex, 0,
position);
855 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
857 dNdxCut1 = (*cutv)[iPlace];
858 dNdxCut2 = (*cutv)[iPlace+1];
860 E1 = fParticleEnergyVector->
Energy(iPlace);
861 E2 = fParticleEnergyVector->
Energy(iPlace+1);
863 W1 = (E2 - scaledTkin)*
W;
864 W2 = (scaledTkin - E1)*
W;
875 transfer = tr1*W1 + tr2*W2;
879 if(transfer < 0.0 ) transfer = 0.0;
898 std::size_t iTransfer;
899 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
901 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
902 x2 = v->
Energy(iTransfer);
903 y2 = (*v)[iTransfer]/x2;
907 x1 = v->
Energy(iTransfer-1);
908 y1 = (*v)[iTransfer-1]/x1;
918 const G4int nbins = 5;
921 for(
G4int i=1; i<=nbins; ++i) {
923 y2 = v->
Value(x2)/x2;
929 energyTransfer = (y2 - y1)*x1*x2/(
position*(x1 - x2) - y1*x1 + y2*x2);
935 return energyTransfer;
940G4double G4PAIPhotData::GetEnergyPhotonTransfer(
G4int coupleIndex,
949 std::size_t iTransfer;
950 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
952 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer)
954 x2 = v->
Energy(iTransfer);
955 y2 = (*v)[iTransfer]/x2;
958 x1 = v->
Energy(iTransfer-1);
959 y1 = (*v)[iTransfer-1]/x1;
976 const G4int nbins = 5;
980 for(
G4int i=1; i<=nbins; ++i)
983 y2 = v->
Value(x2)/x2;
989 energyTransfer = (y2 - y1)*x1*x2/(
position*(x1 - x2) - y1*x1 + y2*x2);
996 return energyTransfer;
1001G4double G4PAIPhotData::GetEnergyPlasmonTransfer(
G4int coupleIndex,
1011 std::size_t iTransfer;
1012 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
1014 for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer)
1016 x2 = v->
Energy(iTransfer);
1017 y2 = (*v)[iTransfer]/x2;
1020 x1 = v->
Energy(iTransfer-1);
1021 y1 = (*v)[iTransfer-1]/x1;
1026 energyTransfer = x1;
1038 const G4int nbins = 5;
1042 for(
G4int i = 1; i <= nbins; ++i )
1045 y2 = v->
Value(x2)/x2;
1053 energyTransfer = (y2 - y1)*x1*x2/(
position*(x1 - x2) - y1*x1 + y2*x2);
1060 return energyTransfer;
G4long G4Poisson(G4double mean)
G4GLOB_DLL std::ostream G4cout
const G4Material * GetMaterial() const
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
G4double SamplePostStepPlasmonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double SampleAlongStepPlasmonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
G4double GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double SamplePostStepPhotonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4PAIPhotData(G4double tmin, G4double tmax, G4int verbose)
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
G4double SampleAlongStepPhotonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double ComputeMaxEnergy(G4double scaledEnergy)
G4double GetIntegralCerenkov(G4int i) const
G4int GetSplineSize() const
G4double GetIntegralPAIdEdx(G4int i) const
G4double GetIntegralPlasmon(G4int i) const
G4double GetSplineEnergy(G4int i) const
G4double GetIntegralPAIxSection(G4int i) const
G4double GetMeanEnergyLoss() const
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
G4double GetMaxEnergy() const
G4double Energy(const std::size_t index) const
G4double Value(const G4double energy, std::size_t &lastidx) const
std::size_t GetVectorLength() const
std::size_t FindBin(const G4double energy, std::size_t idx) const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
const std::vector< G4double > * GetEnergyCutsVector(std::size_t pcIdx) const
static G4ProductionCutsTable * GetProductionCutsTable()
void Initialize(const G4Material *)
G4double GetSandiaMatTablePAI(G4int, G4int) const