61 const G4int nPerDecade = 10;
67 fLowestKineticEnergy =
std::max(tmin, lowestTkin);
68 fHighestKineticEnergy = tmax;
70 if(tmax < 10*fLowestKineticEnergy)
72 fHighestKineticEnergy = 10*fLowestKineticEnergy;
74 else if(tmax > highestTkin)
76 fHighestKineticEnergy =
std::max(highestTkin, 10*fLowestKineticEnergy);
78 fTotBin = (
G4int)(nPerDecade*
79 std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
82 fHighestKineticEnergy,
85 G4cout <<
"### G4PAIPhotData: Nbins= " << fTotBin
86 <<
" Tmin(MeV)= " << fLowestKineticEnergy/
MeV
87 <<
" Tmax(GeV)= " << fHighestKineticEnergy/
GeV
97 size_t n = fPAIxscBank.size();
100 for(
size_t i=0; i<
n; ++i)
104 fPAIxscBank[i]->clearAndDestroy();
105 delete fPAIxscBank[i];
110 fPAIdEdxBank[i]->clearAndDestroy();
111 delete fPAIdEdxBank[i];
114 delete fdEdxTable[i];
115 delete fdNdxCutTable[i];
117 fdNdxCutTable[i] = 0;
120 delete fParticleEnergyVector;
121 fParticleEnergyVector = 0;
135 for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
139 if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
143 G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
144 G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
146 G4cout<<
"G4PAIPhotData::Initialise: "<<
"cut = "<<cut/
keV<<
" keV; cutEl = "
147 <<deltaCutInKineticEnergyNow/
keV<<
" keV; cutPh = "
148 <<photonCutInKineticEnergyNow/
keV<<
" keV"<<
G4endl;
154 fHighestKineticEnergy,
159 fHighestKineticEnergy,
163 fHighestKineticEnergy,
167 fHighestKineticEnergy,
171 fSandia.Initialize(const_cast<G4Material*>(mat));
180 fHighestKineticEnergy,
184 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
189 for (
G4int i = 0; i <= fTotBin; ++i)
191 G4double kinEnergy = fParticleEnergyVector->Energy(i);
196 if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow;
198 fPAIxSection.Initialize( mat, Tmax, bg2, &fSandia);
203 G4int n = fPAIxSection.GetSplineSize();
211 for(
G4int k = 0; k <
n; k++ )
213 G4double t = fPAIxSection.GetSplineEnergy(k+1);
216 t*fPAIxSection.GetIntegralPAIxSection(k+1));
218 t*fPAIxSection.GetIntegralCerenkov(k+1));
220 t*fPAIxSection.GetIntegralPlasmon(k+1));
222 dEdxVector->
PutValue(k, t, fPAIxSection.GetIntegralPAIdEdx(k+1));
226 G4double ionloss = fPAIxSection.GetMeanEnergyLoss();
228 if(ionloss < 0.0) ionloss = 0.0;
230 dEdxMeanVector->
PutValue(i,ionloss);
232 G4double dNdxCut = transferVector->
Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
233 G4double dNdxCutPhoton = photonVector->
Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
234 G4double dNdxCutPlasmon = plasmonVector->
Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
239 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
240 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
241 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
243 dNdxCutVector->
PutValue(i, dNdxCut);
244 dNdxCutPhotonVector->
PutValue(i, dNdxCutPhoton);
245 dNdxCutPlasmonVector->
PutValue(i, dNdxCutPlasmon);
247 dEdxCutVector->
PutValue(i, dEdxCut);
249 PAItransferTable->
insertAt(i,transferVector);
250 PAIphotonTable->
insertAt(i,photonVector);
251 PAIplasmonTable->
insertAt(i,plasmonVector);
252 PAIdEdxTable->
insertAt(i,dEdxVector);
256 fPAIxscBank.push_back(PAItransferTable);
257 fPAIphotonBank.push_back(PAIphotonTable);
258 fPAIplasmonBank.push_back(PAIplasmonTable);
260 fPAIdEdxBank.push_back(PAIdEdxTable);
261 fdEdxTable.push_back(dEdxMeanVector);
263 fdNdxCutTable.push_back(dNdxCutVector);
264 fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
265 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
267 fdEdxCutTable.push_back(dEdxCutVector);
277 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
278 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
281 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
282 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
287 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
288 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
290 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
291 G4double E1 = fParticleEnergyVector->Energy(iPlace);
292 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
301 if( dEdx < 0.) { dEdx = 0.; }
312 G4double cross, xscEl, xscEl2, xscPh, xscPh2;
318 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
319 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
323 if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
324 else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one =
false;
327 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
328 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
333 cross = xscPh + xscEl;
337 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
339 G4double E1 = fParticleEnergyVector->Energy(iPlace);
340 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
349 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
351 E1 = fParticleEnergyVector->Energy(iPlace);
352 E2 = fParticleEnergyVector->Energy(iPlace+1);
355 W1 = (E2 - scaledTkin)*W;
356 W2 = (scaledTkin - E1)*W;
361 cross = xscEl + xscPh;
363 if( cross < 0.0) cross = 0.0;
373 G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
376 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
377 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
381 if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
382 else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one =
false;
385 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
386 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
391 cross = xscPh + xscEl;
395 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
397 G4double E1 = fParticleEnergyVector->Energy(iPlace);
398 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
407 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
409 E1 = fParticleEnergyVector->Energy(iPlace);
410 E2 = fParticleEnergyVector->Energy(iPlace+1);
413 W1 = (E2 - scaledTkin)*W;
414 W2 = (scaledTkin - E1)*W;
419 cross = xscEl + xscPh;
427 plRatio = xscEl/cross;
429 if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
445 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
446 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
450 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
451 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one =
false;
457 dNdxCut1 = (*vcut)[iPlace];
461 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
473 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
474 dNdxCut2 = (*vcut)[iPlace+1];
477 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
479 E1 = fParticleEnergyVector->Energy(iPlace);
480 E2 = fParticleEnergyVector->Energy(iPlace+1);
482 W1 = (E2 - scaledTkin)*W;
483 W2 = (scaledTkin - E1)*W;
484 meanNumber = W1*meanN1 + W2*meanN2;
489 if( meanNumber <= 0.0)
return 0.0;
495 if( 0 == numOfCollisions)
return 0.0;
497 for(
G4int i=0; i< numOfCollisions; ++i)
500 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
501 omega = GetEnergyTransfer(coupleIndex, iPlace, position);
507 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
508 G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
509 omega = omega*W1 + omega2*W2;
514 if( loss > kinEnergy) {
break; }
520 if ( loss > kinEnergy) loss = kinEnergy;
521 else if( loss < 0.) loss = 0.;
537 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
538 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
542 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
543 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one =
false;
549 dNdxCut1 = (*vcut)[iPlace];
553 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
565 v2 = (*(fPAIphotonBank[coupleIndex]))(iPlace+1);
566 dNdxCut2 = (*vcut)[iPlace+1];
569 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
571 E1 = fParticleEnergyVector->Energy(iPlace);
572 E2 = fParticleEnergyVector->Energy(iPlace+1);
574 W1 = (E2 - scaledTkin)*W;
575 W2 = (scaledTkin - E1)*W;
576 meanNumber = W1*meanN1 + W2*meanN2;
581 if( meanNumber <= 0.0)
return 0.0;
587 if( 0 == numOfCollisions)
return 0.0;
589 for(
G4int i=0; i< numOfCollisions; ++i)
592 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
593 omega = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
599 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
600 G4double omega2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
601 omega = omega*W1 + omega2*W2;
606 if( loss > kinEnergy) {
break; }
612 if ( loss > kinEnergy) loss = kinEnergy;
613 else if( loss < 0.) loss = 0.;
629 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
630 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
634 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
635 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one =
false;
641 dNdxCut1 = (*vcut)[iPlace];
645 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
657 v2 = (*(fPAIplasmonBank[coupleIndex]))(iPlace+1);
658 dNdxCut2 = (*vcut)[iPlace+1];
661 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
663 E1 = fParticleEnergyVector->Energy(iPlace);
664 E2 = fParticleEnergyVector->Energy(iPlace+1);
666 W1 = (E2 - scaledTkin)*W;
667 W2 = (scaledTkin - E1)*W;
668 meanNumber = W1*meanN1 + W2*meanN2;
673 if( meanNumber <= 0.0)
return 0.0;
679 if( 0 == numOfCollisions)
return 0.0;
681 for(
G4int i=0; i< numOfCollisions; ++i)
684 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
685 omega = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
691 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
692 G4double omega2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
693 omega = omega*W1 + omega2*W2;
698 if( loss > kinEnergy) {
break; }
704 if ( loss > kinEnergy) loss = kinEnergy;
705 else if( loss < 0.) loss = 0.;
722 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
730 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
732 position = (*cutv)[nPlace]*rand;
733 transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
735 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
737 position = (*cutv)[0]*rand;
738 transfer = GetEnergyTransfer(coupleIndex, 0, position);
742 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
744 dNdxCut1 = (*cutv)[iPlace];
745 dNdxCut2 = (*cutv)[iPlace+1];
747 E1 = fParticleEnergyVector->Energy(iPlace);
748 E2 = fParticleEnergyVector->Energy(iPlace+1);
750 W1 = (E2 - scaledTkin)*W;
751 W2 = (scaledTkin - E1)*W;
756 position = dNdxCut1*rand;
757 G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
759 position = dNdxCut2*rand;
760 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
762 transfer = tr1*W1 + tr2*W2;
765 if(transfer < 0.0 ) { transfer = 0.0; }
778 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
787 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
789 position = (*cutv)[nPlace]*rand;
790 transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace, position);
792 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
794 position = (*cutv)[0]*rand;
795 transfer = GetEnergyPhotonTransfer(coupleIndex, 0, position);
799 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
801 dNdxCut1 = (*cutv)[iPlace];
802 dNdxCut2 = (*cutv)[iPlace+1];
804 E1 = fParticleEnergyVector->Energy(iPlace);
805 E2 = fParticleEnergyVector->Energy(iPlace+1);
807 W1 = (E2 - scaledTkin)*W;
808 W2 = (scaledTkin - E1)*W;
813 position = dNdxCut1*rand;
815 G4double tr1 = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
817 position = dNdxCut2*rand;
818 G4double tr2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
820 transfer = tr1*W1 + tr2*W2;
823 if(transfer < 0.0 ) { transfer = 0.0; }
836 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
844 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
846 position = (*cutv)[nPlace]*rand;
847 transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace, position);
849 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
851 position = (*cutv)[0]*rand;
852 transfer = GetEnergyPlasmonTransfer(coupleIndex, 0, position);
856 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
858 dNdxCut1 = (*cutv)[iPlace];
859 dNdxCut2 = (*cutv)[iPlace+1];
861 E1 = fParticleEnergyVector->Energy(iPlace);
862 E2 = fParticleEnergyVector->Energy(iPlace+1);
864 W1 = (E2 - scaledTkin)*W;
865 W2 = (scaledTkin - E1)*W;
870 position = dNdxCut1*rand;
871 G4double tr1 = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
873 position = dNdxCut2*rand;
874 G4double tr2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
876 transfer = tr1*W1 + tr2*W2;
880 if(transfer < 0.0 ) transfer = 0.0;
895 if(position*v->
Energy(0) >= (*v)[0]) {
return v->
Energy(0); }
900 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
902 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
903 x2 = v->
Energy(iTransfer);
904 y2 = (*v)[iTransfer]/x2;
905 if(position >= y2) {
break; }
908 x1 = v->
Energy(iTransfer-1);
909 y1 = (*v)[iTransfer-1]/x1;
919 const G4int nbins = 5;
922 for(
G4int i=1; i<=nbins; ++i) {
924 y2 = v->
Value(x2)/x2;
925 if(position >= y2) {
break; }
930 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
936 return energyTransfer;
941 G4double G4PAIPhotData::GetEnergyPhotonTransfer(
G4int coupleIndex,
946 if(position*v->
Energy(0) >= (*v)[0])
return v->
Energy(0);
951 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
953 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer)
955 x2 = v->
Energy(iTransfer);
956 y2 = (*v)[iTransfer]/x2;
957 if(position >= y2)
break;
959 x1 = v->
Energy(iTransfer-1);
960 y1 = (*v)[iTransfer-1]/x1;
977 const G4int nbins = 5;
981 for(
G4int i=1; i<=nbins; ++i)
984 y2 = v->
Value(x2)/x2;
985 if(position >= y2) {
break; }
990 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
997 return energyTransfer;
1002 G4double G4PAIPhotData::GetEnergyPlasmonTransfer(
G4int coupleIndex,
1008 if( position*v->
Energy(0) >= (*v)[0] )
return v->
Energy(0);
1013 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
1015 for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer)
1017 x2 = v->
Energy(iTransfer);
1018 y2 = (*v)[iTransfer]/x2;
1019 if(position >= y2)
break;
1021 x1 = v->
Energy(iTransfer-1);
1022 y1 = (*v)[iTransfer-1]/x1;
1027 energyTransfer = x1;
1039 const G4int nbins = 5;
1043 for(
G4int i = 1; i <= nbins; ++i )
1046 y2 = v->
Value(x2)/x2;
1048 if(position >= y2)
break;
1054 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
1061 return energyTransfer;
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double SamplePostStepPhotonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4long G4Poisson(G4double mean)
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
static constexpr double proton_mass_c2
G4double GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
void PutValue(size_t index, G4double energy, G4double dataValue)
G4double SampleAlongStepPhotonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4PAIPhotData(G4double tmin, G4double tmax, G4int verbose)
size_t GetVectorLength() const
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
static constexpr double TeV
G4GLOB_DLL std::ostream G4cout
G4double ComputeMaxEnergy(G4double scaledEnergy)
size_t GetTableSize() const
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
void PutValue(size_t index, G4double theValue)
static constexpr double eV
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double SamplePostStepPlasmonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
G4double SampleAlongStepPlasmonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
static constexpr double GeV
void insertAt(size_t, G4PhysicsVector *)
static constexpr double MeV
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
static constexpr double keV
const XML_Char XML_Content * model
const G4Material * GetMaterial() const