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
98 size_t n = fPAIxscBank.size();
101 for(
size_t i=0; i<
n; ++i)
106 delete fPAIxscBank[i];
111 delete fPAIdEdxBank[i];
129 for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
133 if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
137 G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
138 G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
140 G4cout<<
"G4PAIPhotData::Initialise: "<<
"cut = "<<cut/
keV<<
" keV; cutEl = "
141 <<deltaCutInKineticEnergyNow/
keV<<
" keV; cutPh = "
142 <<photonCutInKineticEnergyNow/
keV<<
" keV"<<
G4endl;
147 fSandia.Initialize(const_cast<G4Material*>(mat));
156 fHighestKineticEnergy,
160 fHighestKineticEnergy,
165 fHighestKineticEnergy,
169 fHighestKineticEnergy,
173 fHighestKineticEnergy,
177 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
182 for (
G4int i = 0; i <= fTotBin; ++i)
184 G4double kinEnergy = fParticleEnergyVector->Energy(i);
189 if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow;
191 fPAIxSection.Initialize( mat, Tmax, bg2, &fSandia);
196 G4int n = fPAIxSection.GetSplineSize();
204 for(
G4int k = 0; k <
n; k++ )
206 G4double t = fPAIxSection.GetSplineEnergy(k+1);
209 t*fPAIxSection.GetIntegralPAIxSection(k+1));
211 t*fPAIxSection.GetIntegralCerenkov(k+1));
213 t*fPAIxSection.GetIntegralPlasmon(k+1));
215 dEdxVector->
PutValue(k, t, fPAIxSection.GetIntegralPAIdEdx(k+1));
219 G4double ionloss = fPAIxSection.GetMeanEnergyLoss();
221 if(ionloss < 0.0) ionloss = 0.0;
223 dEdxMeanVector->
PutValue(i,ionloss);
225 G4double dNdxCut = transferVector->
Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
226 G4double dNdxCutPhoton = photonVector->
Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
227 G4double dNdxCutPlasmon = plasmonVector->
Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
232 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
233 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
234 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
236 dNdxCutVector->
PutValue(i, dNdxCut);
237 dNdxCutPhotonVector->
PutValue(i, dNdxCutPhoton);
238 dNdxCutPlasmonVector->
PutValue(i, dNdxCutPlasmon);
240 dEdxCutVector->
PutValue(i, dEdxCut);
242 PAItransferTable->
insertAt(i,transferVector);
243 PAIphotonTable->
insertAt(i,photonVector);
244 PAIplasmonTable->
insertAt(i,plasmonVector);
245 PAIdEdxTable->
insertAt(i,dEdxVector);
249 fPAIxscBank.push_back(PAItransferTable);
250 fPAIphotonBank.push_back(PAIphotonTable);
251 fPAIplasmonBank.push_back(PAIplasmonTable);
253 fPAIdEdxBank.push_back(PAIdEdxTable);
254 fdEdxTable.push_back(dEdxMeanVector);
256 fdNdxCutTable.push_back(dNdxCutVector);
257 fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
258 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
260 fdEdxCutTable.push_back(dEdxCutVector);
270 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
271 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
274 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
275 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
280 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
281 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
283 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
284 G4double E1 = fParticleEnergyVector->Energy(iPlace);
285 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
294 if( dEdx < 0.) { dEdx = 0.; }
305 G4double cross, xscEl, xscEl2, xscPh, xscPh2;
311 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
312 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
316 if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
317 else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one =
false;
320 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
321 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
326 cross = xscPh + xscEl;
330 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
332 G4double E1 = fParticleEnergyVector->Energy(iPlace);
333 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
342 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
344 E1 = fParticleEnergyVector->Energy(iPlace);
345 E2 = fParticleEnergyVector->Energy(iPlace+1);
348 W1 = (E2 - scaledTkin)*W;
349 W2 = (scaledTkin - E1)*W;
354 cross = xscEl + xscPh;
356 if( cross < 0.0) cross = 0.0;
366 G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
369 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
370 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
374 if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
375 else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one =
false;
378 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
379 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
384 cross = xscPh + xscEl;
388 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
390 G4double E1 = fParticleEnergyVector->Energy(iPlace);
391 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
400 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
402 E1 = fParticleEnergyVector->Energy(iPlace);
403 E2 = fParticleEnergyVector->Energy(iPlace+1);
406 W1 = (E2 - scaledTkin)*W;
407 W2 = (scaledTkin - E1)*W;
412 cross = xscEl + xscPh;
420 plRatio = xscEl/cross;
422 if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
438 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
439 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
443 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
444 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one =
false;
450 dNdxCut1 = (*vcut)[iPlace];
454 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
466 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
467 dNdxCut2 = (*vcut)[iPlace+1];
470 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
472 E1 = fParticleEnergyVector->Energy(iPlace);
473 E2 = fParticleEnergyVector->Energy(iPlace+1);
475 W1 = (E2 - scaledTkin)*W;
476 W2 = (scaledTkin - E1)*W;
477 meanNumber = W1*meanN1 + W2*meanN2;
482 if( meanNumber <= 0.0)
return 0.0;
488 if( 0 == numOfCollisions)
return 0.0;
490 for(
G4int i=0; i< numOfCollisions; ++i)
493 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
494 omega = GetEnergyTransfer(coupleIndex, iPlace, position);
500 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
501 G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
502 omega = omega*W1 + omega2*W2;
507 if( loss > kinEnergy) {
break; }
513 if ( loss > kinEnergy) loss = kinEnergy;
514 else if( loss < 0.) loss = 0.;
530 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
531 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
535 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
536 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one =
false;
542 dNdxCut1 = (*vcut)[iPlace];
546 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
558 v2 = (*(fPAIphotonBank[coupleIndex]))(iPlace+1);
559 dNdxCut2 = (*vcut)[iPlace+1];
562 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
564 E1 = fParticleEnergyVector->Energy(iPlace);
565 E2 = fParticleEnergyVector->Energy(iPlace+1);
567 W1 = (E2 - scaledTkin)*W;
568 W2 = (scaledTkin - E1)*W;
569 meanNumber = W1*meanN1 + W2*meanN2;
574 if( meanNumber <= 0.0)
return 0.0;
580 if( 0 == numOfCollisions)
return 0.0;
582 for(
G4int i=0; i< numOfCollisions; ++i)
585 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
586 omega = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
592 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
593 G4double omega2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
594 omega = omega*W1 + omega2*W2;
599 if( loss > kinEnergy) {
break; }
605 if ( loss > kinEnergy) loss = kinEnergy;
606 else if( loss < 0.) loss = 0.;
622 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
623 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
627 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
628 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one =
false;
634 dNdxCut1 = (*vcut)[iPlace];
638 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
650 v2 = (*(fPAIplasmonBank[coupleIndex]))(iPlace+1);
651 dNdxCut2 = (*vcut)[iPlace+1];
654 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
656 E1 = fParticleEnergyVector->Energy(iPlace);
657 E2 = fParticleEnergyVector->Energy(iPlace+1);
659 W1 = (E2 - scaledTkin)*W;
660 W2 = (scaledTkin - E1)*W;
661 meanNumber = W1*meanN1 + W2*meanN2;
666 if( meanNumber <= 0.0)
return 0.0;
672 if( 0 == numOfCollisions)
return 0.0;
674 for(
G4int i=0; i< numOfCollisions; ++i)
677 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
678 omega = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
684 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
685 G4double omega2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
686 omega = omega*W1 + omega2*W2;
691 if( loss > kinEnergy) {
break; }
697 if ( loss > kinEnergy) loss = kinEnergy;
698 else if( loss < 0.) loss = 0.;
715 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
723 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
725 position = (*cutv)[nPlace]*rand;
726 transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
728 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
730 position = (*cutv)[0]*rand;
731 transfer = GetEnergyTransfer(coupleIndex, 0, position);
735 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
737 dNdxCut1 = (*cutv)[iPlace];
738 dNdxCut2 = (*cutv)[iPlace+1];
740 E1 = fParticleEnergyVector->Energy(iPlace);
741 E2 = fParticleEnergyVector->Energy(iPlace+1);
743 W1 = (E2 - scaledTkin)*W;
744 W2 = (scaledTkin - E1)*W;
749 position = dNdxCut1*rand;
750 G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
752 position = dNdxCut2*rand;
753 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
755 transfer = tr1*W1 + tr2*W2;
758 if(transfer < 0.0 ) { transfer = 0.0; }
771 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
780 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
782 position = (*cutv)[nPlace]*rand;
783 transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace, position);
785 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
787 position = (*cutv)[0]*rand;
788 transfer = GetEnergyPhotonTransfer(coupleIndex, 0, position);
792 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
794 dNdxCut1 = (*cutv)[iPlace];
795 dNdxCut2 = (*cutv)[iPlace+1];
797 E1 = fParticleEnergyVector->Energy(iPlace);
798 E2 = fParticleEnergyVector->Energy(iPlace+1);
800 W1 = (E2 - scaledTkin)*W;
801 W2 = (scaledTkin - E1)*W;
806 position = dNdxCut1*rand;
808 G4double tr1 = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
810 position = dNdxCut2*rand;
811 G4double tr2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
813 transfer = tr1*W1 + tr2*W2;
816 if(transfer < 0.0 ) { transfer = 0.0; }
829 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
837 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
839 position = (*cutv)[nPlace]*rand;
840 transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace, position);
842 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
844 position = (*cutv)[0]*rand;
845 transfer = GetEnergyPlasmonTransfer(coupleIndex, 0, position);
849 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
851 dNdxCut1 = (*cutv)[iPlace];
852 dNdxCut2 = (*cutv)[iPlace+1];
854 E1 = fParticleEnergyVector->Energy(iPlace);
855 E2 = fParticleEnergyVector->Energy(iPlace+1);
857 W1 = (E2 - scaledTkin)*W;
858 W2 = (scaledTkin - E1)*W;
863 position = dNdxCut1*rand;
864 G4double tr1 = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
866 position = dNdxCut2*rand;
867 G4double tr2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
869 transfer = tr1*W1 + tr2*W2;
873 if(transfer < 0.0 ) transfer = 0.0;
888 if(position*v->
Energy(0) >= (*v)[0]) {
return v->
Energy(0); }
895 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
897 y2 = (*v)[iTransfer]/
x2;
898 if(position >=
y2) {
break; }
902 y1 = (*v)[iTransfer-1]/
x1;
912 const G4int nbins = 5;
915 for(
G4int i=1; i<=nbins; ++i) {
918 if(position >=
y2) {
break; }
929 return energyTransfer;
934 G4double G4PAIPhotData::GetEnergyPhotonTransfer(
G4int coupleIndex,
939 if(position*v->
Energy(0) >= (*v)[0])
return v->
Energy(0);
946 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer)
949 y2 = (*v)[iTransfer]/
x2;
950 if(position >=
y2)
break;
953 y1 = (*v)[iTransfer-1]/
x1;
970 const G4int nbins = 5;
974 for(
G4int i=1; i<=nbins; ++i)
978 if(position >=
y2) {
break; }
990 return energyTransfer;
995 G4double G4PAIPhotData::GetEnergyPlasmonTransfer(
G4int coupleIndex,
1001 if( position*v->
Energy(0) >= (*v)[0] )
return v->
Energy(0);
1008 for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer)
1011 y2 = (*v)[iTransfer]/
x2;
1012 if(position >=
y2)
break;
1015 y1 = (*v)[iTransfer-1]/
x1;
1020 energyTransfer =
x1;
1032 const G4int nbins = 5;
1036 for(
G4int i = 1; i <= nbins; ++i )
1041 if(position >=
y2)
break;
1054 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
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
G4double GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
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 *)
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)
const XML_Char XML_Content * model
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
void insertAt(size_t, G4PhysicsVector *)
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
const G4Material * GetMaterial() const