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;
153 fSandia.Initialize(const_cast<G4Material*>(mat));
162 fHighestKineticEnergy,
166 fHighestKineticEnergy,
171 fHighestKineticEnergy,
175 fHighestKineticEnergy,
179 fHighestKineticEnergy,
183 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
188 for (
G4int i = 0; i <= fTotBin; ++i)
190 G4double kinEnergy = fParticleEnergyVector->Energy(i);
192 G4double tau = kinEnergy/proton_mass_c2;
195 if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow;
197 fPAIxSection.Initialize( mat, Tmax, bg2, &fSandia);
202 G4int n = fPAIxSection.GetSplineSize();
210 for(
G4int k = 0; k <
n; k++ )
212 G4double t = fPAIxSection.GetSplineEnergy(k+1);
215 t*fPAIxSection.GetIntegralPAIxSection(k+1));
217 t*fPAIxSection.GetIntegralCerenkov(k+1));
219 t*fPAIxSection.GetIntegralPlasmon(k+1));
221 dEdxVector->
PutValue(k, t, fPAIxSection.GetIntegralPAIdEdx(k+1));
225 G4double ionloss = fPAIxSection.GetMeanEnergyLoss();
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;
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 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
277 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
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);
290 G4double E1 = fParticleEnergyVector->Energy(iPlace);
291 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
300 if( dEdx < 0.) { dEdx = 0.; }
311 G4double cross, xscEl, xscEl2, xscPh, xscPh2;
317 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
318 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
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);
338 G4double E1 = fParticleEnergyVector->Energy(iPlace);
339 G4double E2 = fParticleEnergyVector->Energy(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 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
376 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
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);
396 G4double E1 = fParticleEnergyVector->Energy(iPlace);
397 G4double E2 = fParticleEnergyVector->Energy(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 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
445 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
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;
507 G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
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 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
537 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
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;
599 G4double omega2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
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 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
629 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
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.;
721 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
729 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
731 position = (*cutv)[nPlace]*rand;
732 transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
734 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
736 position = (*cutv)[0]*rand;
737 transfer = GetEnergyTransfer(coupleIndex, 0, position);
741 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;
755 position = dNdxCut1*rand;
756 G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
758 position = dNdxCut2*rand;
759 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
761 transfer = tr1*W1 + tr2*W2;
764 if(transfer < 0.0 ) { transfer = 0.0; }
777 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
786 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
788 position = (*cutv)[nPlace]*rand;
789 transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace, position);
791 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
793 position = (*cutv)[0]*rand;
794 transfer = GetEnergyPhotonTransfer(coupleIndex, 0, position);
798 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;
812 position = dNdxCut1*rand;
814 G4double tr1 = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
816 position = dNdxCut2*rand;
817 G4double tr2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
819 transfer = tr1*W1 + tr2*W2;
822 if(transfer < 0.0 ) { transfer = 0.0; }
835 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
843 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
845 position = (*cutv)[nPlace]*rand;
846 transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace, position);
848 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
850 position = (*cutv)[0]*rand;
851 transfer = GetEnergyPlasmonTransfer(coupleIndex, 0, position);
855 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;
869 position = dNdxCut1*rand;
870 G4double tr1 = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
872 position = dNdxCut2*rand;
873 G4double tr2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
875 transfer = tr1*W1 + tr2*W2;
879 if(transfer < 0.0 ) transfer = 0.0;
894 if(position*v->
Energy(0) >= (*v)[0]) {
return v->
Energy(0); }
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;
904 if(position >= y2) {
break; }
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;
924 if(position >= y2) {
break; }
929 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
935 return energyTransfer;
945 if(position*v->
Energy(0) >= (*v)[0])
return v->
Energy(0);
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;
956 if(position >= y2)
break;
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;
984 if(position >= y2) {
break; }
989 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
996 return energyTransfer;
1007 if( position*v->
Energy(0) >= (*v)[0] )
return v->
Energy(0);
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;
1018 if(position >= y2)
break;
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;
1047 if(position >= y2)
break;
1053 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
1060 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 GetEnergyTransfer(G4int coupleIndex, size_t iPlace, G4double position) const
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 *)
G4double GetEnergyPlasmonTransfer(G4int coupleIndex, size_t iPlace, G4double position) const
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)
G4double GetEnergyPhotonTransfer(G4int coupleIndex, size_t iPlace, G4double position) const
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