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];
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,
226 dEdxMeanVector->PutValue(i,ionloss);
228 G4double dNdxCut = transferVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
229 G4double dNdxCutPhoton = photonVector->Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
230 G4double dNdxCutPlasmon = plasmonVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
232 G4double dEdxCut = dEdxVector->Value(cut)/cut;
235 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
236 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
237 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
239 dNdxCutVector->PutValue(i, dNdxCut);
240 dNdxCutPhotonVector->PutValue(i, dNdxCutPhoton);
241 dNdxCutPlasmonVector->PutValue(i, dNdxCutPlasmon);
243 dEdxCutVector->PutValue(i, dEdxCut);
245 PAItransferTable->insertAt(i,transferVector);
246 PAIphotonTable->insertAt(i,photonVector);
247 PAIplasmonTable->insertAt(i,plasmonVector);
248 PAIdEdxTable->insertAt(i,dEdxVector);
252 fPAIxscBank.push_back(PAItransferTable);
253 fPAIphotonBank.push_back(PAIphotonTable);
254 fPAIplasmonBank.push_back(PAIplasmonTable);
256 fPAIdEdxBank.push_back(PAIdEdxTable);
257 fdEdxTable.push_back(dEdxMeanVector);
259 fdNdxCutTable.push_back(dNdxCutVector);
260 fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
261 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
263 fdEdxCutTable.push_back(dEdxCutVector);
273 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
277 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
278 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
283 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
284 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
286 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
297 if( dEdx < 0.) { dEdx = 0.; }
308 G4double cross, xscEl, xscEl2, xscPh, xscPh2;
314 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
319 if( scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
320 else if( scaledTkin > fParticleEnergyVector->
Energy(0) ) one =
false;
323 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
324 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
329 cross = xscPh + xscEl;
333 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
345 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
347 E1 = fParticleEnergyVector->
Energy(iPlace);
348 E2 = fParticleEnergyVector->
Energy(iPlace+1);
351 W1 = (E2 - scaledTkin)*
W;
352 W2 = (scaledTkin - E1)*
W;
357 cross = xscEl + xscPh;
359 if( cross < 0.0) cross = 0.0;
369 G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
372 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
377 if( scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
378 else if( scaledTkin > fParticleEnergyVector->
Energy(0) ) one =
false;
381 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
382 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
387 cross = xscPh + xscEl;
391 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
403 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
405 E1 = fParticleEnergyVector->
Energy(iPlace);
406 E2 = fParticleEnergyVector->
Energy(iPlace+1);
409 W1 = (E2 - scaledTkin)*
W;
410 W2 = (scaledTkin - E1)*
W;
415 cross = xscEl + xscPh;
423 plRatio = xscEl/cross;
425 if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
441 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
446 if (scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
447 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) one =
false;
453 dNdxCut1 = (*vcut)[iPlace];
457 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
469 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
470 dNdxCut2 = (*vcut)[iPlace+1];
473 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
475 E1 = fParticleEnergyVector->
Energy(iPlace);
476 E2 = fParticleEnergyVector->
Energy(iPlace+1);
478 W1 = (E2 - scaledTkin)*
W;
479 W2 = (scaledTkin - E1)*
W;
480 meanNumber = W1*meanN1 + W2*meanN2;
485 if( meanNumber <= 0.0)
return 0.0;
491 if( 0 == numOfCollisions)
return 0.0;
493 for(
G4int i=0; i< numOfCollisions; ++i)
496 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
497 omega = GetEnergyTransfer(coupleIndex, iPlace,
position);
503 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
505 omega = omega*W1 + omega2*W2;
510 if( loss > kinEnergy) {
break; }
516 if ( loss > kinEnergy) loss = kinEnergy;
517 else if( loss < 0.) loss = 0.;
533 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
538 if (scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
539 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) one =
false;
545 dNdxCut1 = (*vcut)[iPlace];
549 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
561 v2 = (*(fPAIphotonBank[coupleIndex]))(iPlace+1);
562 dNdxCut2 = (*vcut)[iPlace+1];
565 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
567 E1 = fParticleEnergyVector->
Energy(iPlace);
568 E2 = fParticleEnergyVector->
Energy(iPlace+1);
570 W1 = (E2 - scaledTkin)*
W;
571 W2 = (scaledTkin - E1)*
W;
572 meanNumber = W1*meanN1 + W2*meanN2;
577 if( meanNumber <= 0.0)
return 0.0;
583 if( 0 == numOfCollisions)
return 0.0;
585 for(
G4int i=0; i< numOfCollisions; ++i)
588 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
589 omega = GetEnergyPhotonTransfer(coupleIndex, iPlace,
position);
595 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
597 omega = omega*W1 + omega2*W2;
602 if( loss > kinEnergy) {
break; }
608 if ( loss > kinEnergy) loss = kinEnergy;
609 else if( loss < 0.) loss = 0.;
625 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
630 if (scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) iPlace = nPlace;
631 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) one =
false;
637 dNdxCut1 = (*vcut)[iPlace];
641 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
653 v2 = (*(fPAIplasmonBank[coupleIndex]))(iPlace+1);
654 dNdxCut2 = (*vcut)[iPlace+1];
657 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
659 E1 = fParticleEnergyVector->
Energy(iPlace);
660 E2 = fParticleEnergyVector->
Energy(iPlace+1);
662 W1 = (E2 - scaledTkin)*
W;
663 W2 = (scaledTkin - E1)*
W;
664 meanNumber = W1*meanN1 + W2*meanN2;
669 if( meanNumber <= 0.0)
return 0.0;
675 if( 0 == numOfCollisions)
return 0.0;
677 for(
G4int i=0; i< numOfCollisions; ++i)
680 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
681 omega = GetEnergyPlasmonTransfer(coupleIndex, iPlace,
position);
687 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
688 G4double omega2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1,
position);
689 omega = omega*W1 + omega2*W2;
694 if( loss > kinEnergy) {
break; }
700 if ( loss > kinEnergy) loss = kinEnergy;
701 else if( loss < 0.) loss = 0.;
726 if( scaledTkin >= fParticleEnergyVector->
GetMaxEnergy())
729 transfer = GetEnergyTransfer(coupleIndex, nPlace,
position);
731 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
734 transfer = GetEnergyTransfer(coupleIndex, 0,
position);
738 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
740 dNdxCut1 = (*cutv)[iPlace];
741 dNdxCut2 = (*cutv)[iPlace+1];
743 E1 = fParticleEnergyVector->
Energy(iPlace);
744 E2 = fParticleEnergyVector->
Energy(iPlace+1);
746 W1 = (E2 - scaledTkin)*
W;
747 W2 = (scaledTkin - E1)*
W;
758 transfer = tr1*W1 + tr2*W2;
761 if(transfer < 0.0 ) { transfer = 0.0; }
783 if( scaledTkin >= fParticleEnergyVector->
GetMaxEnergy())
786 transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace,
position);
788 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
791 transfer = GetEnergyPhotonTransfer(coupleIndex, 0,
position);
795 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
797 dNdxCut1 = (*cutv)[iPlace];
798 dNdxCut2 = (*cutv)[iPlace+1];
800 E1 = fParticleEnergyVector->
Energy(iPlace);
801 E2 = fParticleEnergyVector->
Energy(iPlace+1);
803 W1 = (E2 - scaledTkin)*
W;
804 W2 = (scaledTkin - E1)*
W;
816 transfer = tr1*W1 + tr2*W2;
819 if(transfer < 0.0 ) { transfer = 0.0; }
840 if( scaledTkin >= fParticleEnergyVector->
GetMaxEnergy())
843 transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace,
position);
845 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
848 transfer = GetEnergyPlasmonTransfer(coupleIndex, 0,
position);
852 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
854 dNdxCut1 = (*cutv)[iPlace];
855 dNdxCut2 = (*cutv)[iPlace+1];
857 E1 = fParticleEnergyVector->
Energy(iPlace);
858 E2 = fParticleEnergyVector->
Energy(iPlace+1);
860 W1 = (E2 - scaledTkin)*
W;
861 W2 = (scaledTkin - E1)*
W;
872 transfer = tr1*W1 + tr2*W2;
876 if(transfer < 0.0 ) transfer = 0.0;
895 std::size_t iTransfer;
896 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
898 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
899 x2 = v->
Energy(iTransfer);
900 y2 = (*v)[iTransfer]/x2;
904 x1 = v->
Energy(iTransfer-1);
905 y1 = (*v)[iTransfer-1]/x1;
915 const G4int nbins = 5;
918 for(
G4int i=1; i<=nbins; ++i) {
920 y2 = v->
Value(x2)/x2;
926 energyTransfer = (y2 - y1)*x1*x2/(
position*(x1 - x2) - y1*x1 + y2*x2);
932 return energyTransfer;
937G4double G4PAIPhotData::GetEnergyPhotonTransfer(
G4int coupleIndex,
946 std::size_t iTransfer;
947 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
949 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer)
951 x2 = v->
Energy(iTransfer);
952 y2 = (*v)[iTransfer]/x2;
955 x1 = v->
Energy(iTransfer-1);
956 y1 = (*v)[iTransfer-1]/x1;
973 const G4int nbins = 5;
977 for(
G4int i=1; i<=nbins; ++i)
980 y2 = v->
Value(x2)/x2;
986 energyTransfer = (y2 - y1)*x1*x2/(
position*(x1 - x2) - y1*x1 + y2*x2);
993 return energyTransfer;
998G4double G4PAIPhotData::GetEnergyPlasmonTransfer(
G4int coupleIndex,
1008 std::size_t iTransfer;
1009 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
1011 for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer)
1013 x2 = v->
Energy(iTransfer);
1014 y2 = (*v)[iTransfer]/x2;
1017 x1 = v->
Energy(iTransfer-1);
1018 y1 = (*v)[iTransfer-1]/x1;
1023 energyTransfer = x1;
1035 const G4int nbins = 5;
1039 for(
G4int i = 1; i <= nbins; ++i )
1042 y2 = v->
Value(x2)/x2;
1050 energyTransfer = (y2 - y1)*x1*x2/(
position*(x1 - x2) - y1*x1 + y2*x2);
1057 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