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); }
902 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
904 y2 = (*v)[iTransfer]/
x2;
905 if(position >=
y2) {
break; }
909 y1 = (*v)[iTransfer-1]/
x1;
919 const G4int nbins = 5;
922 for(
G4int i=1; i<=nbins; ++i) {
925 if(position >=
y2) {
break; }
936 return energyTransfer;
946 if(position*v->
Energy(0) >= (*v)[0])
return v->
Energy(0);
953 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer)
956 y2 = (*v)[iTransfer]/
x2;
957 if(position >=
y2)
break;
960 y1 = (*v)[iTransfer-1]/
x1;
977 const G4int nbins = 5;
981 for(
G4int i=1; i<=nbins; ++i)
985 if(position >=
y2) {
break; }
997 return energyTransfer;
1008 if( position*v->
Energy(0) >= (*v)[0] )
return v->
Energy(0);
1015 for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer)
1018 y2 = (*v)[iTransfer]/
x2;
1019 if(position >=
y2)
break;
1022 y1 = (*v)[iTransfer-1]/
x1;
1027 energyTransfer =
x1;
1039 const G4int nbins = 5;
1043 for(
G4int i = 1; i <= nbins; ++i )
1048 if(position >=
y2)
break;
1061 return energyTransfer;
G4double Energy(size_t index) const
G4double SamplePostStepPlasmonTransfer(G4int coupleIndex, G4double scaledTkin) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetEnergyPlasmonTransfer(G4int coupleIndex, size_t iPlace, G4double position) const
static constexpr double MeV
Float_t y1[n_points_granero]
static constexpr double keV
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin) const
Float_t x1[n_points_granero]
void insertAt(size_t, G4PhysicsVector *)
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
Float_t y2[n_points_geant4]
static constexpr double proton_mass_c2
static constexpr double TeV
G4PAIPhotData(G4double tmin, G4double tmax, G4int verbose)
static constexpr double eV
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
size_t GetTableSize() const
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
static G4ProductionCutsTable * GetProductionCutsTable()
G4double SamplePostStepPhotonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double SampleAlongStepPlasmonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double SampleAlongStepPhotonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
G4GLOB_DLL std::ostream G4cout
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4long G4Poisson(G4double mean)
G4double ComputeMaxEnergy(G4double scaledEnergy)
const G4Material * GetMaterial() const
void PutValue(size_t index, G4double energy, G4double dataValue)
Float_t x2[n_points_geant4]
const XML_Char XML_Content * model
static constexpr double GeV
void PutValue(size_t index, G4double theValue)
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
G4double GetEnergyPhotonTransfer(G4int coupleIndex, size_t iPlace, G4double position) const
size_t GetVectorLength() const
G4double GetEnergyTransfer(G4int coupleIndex, size_t iPlace, G4double position) const