68 fCache.Get()->currentMeanEnergy = -2;
69 fCache.Get()->fresh =
true;
118 if (std::getenv(
"G4PHPTEST") )
G4cout <<
" G4ParticleHPContAngularPar::Sample "
119 << anEnergy <<
" " << massCode <<
" "
141 "G4ParticleHPContAngularPar: Unknown ion case 1");
151 if (angularRep == 1) {
155 if (
fCache.Get()->fresh ==
true) {
160 fCache.Get()->fresh =
false;
161 if (std::getenv(
"G4PHPTEST") )
162 G4cout <<
" G4ParticleHPContAngularPar::Sample fresh remaining_energy "
171 if (std::getenv(
"G4PHPTEST") )
172 G4cout <<
" G4ParticleHPContAngularPar::Sample cheat case for small remaining_energy "
180 fCache.Get()->remaining_energy =
183 if (std::getenv(
"G4PHPTEST") )
184 G4cout <<
" G4ParticleHPContAngularPar::Sample max line or grid remaining_energy "
197 running[j+1] = running[j] + delta;
234 "G4ParticleHPContAngularPar: Unexpected non zero value of theAngular[nEnergies-1].GetValue(0)");
239 running[j+1] = running[j] + ( ( e_high - e_low ) * delta1);
244 random *= (tot_prob_DIS + tot_prob_CON);
252 if (random < running[ j+1 ] ) {
272 if (random < running[ j ] ) {
313 if (std::getenv(
"G4PHPTEST") )
314 G4cout <<
" G4ParticleHPContAngularPar::Sample remaining_energy - fsEnergy "
315 <<
fCache.Get()->remaining_energy <<
" - " << fsEnergy <<
G4endl;
316 fCache.Get()->remaining_energy -= fsEnergy;
323 if (
fCache.Get()->fresh ==
true) {
325 fCache.Get()->fresh =
false;
326 if (std::getenv(
"G4PHPTEST") )
327 G4cout <<
" G4ParticleHPContAngularPar::Sample 2nd fresh remaining_energy "
336 running[i]=running[i-1];
353 fCache.Get()->currentMeanEnergy = 0.0;
362 if (random < running [i]/running[
nEnergies-1] )
break;
424 if (std::getenv(
"G4PHPTEST") )
425 G4cout <<
" G4ParticleHPContAngularPar::Sample 2nd remaining_energy - fsEnergy "
426 <<
fCache.Get()->remaining_energy <<
" - " << fsEnergy <<
G4endl;
427 fCache.Get()->remaining_energy -= fsEnergy;
432 }
else if (angularRep == 2) {
438 if (std::getenv(
"G4PHPTEST") )
441 if (j != 0) running[j] = running[j-1];
448 if (std::getenv(
"G4PHPTEST") )
449 G4cout <<
" G4ParticleHPContAngularPar::Sample " << j <<
" running " << running[j]
457 fCache.Get()->currentMeanEnergy = 0.0;
465 if (randkal < running[j]/running[
nEnergies-1] )
break;
470 if (itt == 0) itt = 1;
480 if (std::getenv(
"G4PHPTEST") )
481 G4cout << itt <<
" G4particleHPContAngularPar fsEnergy " << fsEnergy <<
" "
483 << x2 <<
" y " << y1 <<
" " << y2 <<
G4endl;
488 fsEnergy, y1, y2, cLow, cHigh);
490 if (compoundFraction > 1.0) compoundFraction = 1.0;
492 if (std::getenv(
"G4PHPTEST") )
493 G4cout << itt <<
" G4particleHPContAngularPar compoundFraction " << compoundFraction
495 <<
" y " << y1 <<
" " << y2 <<
" cLH " << cLow <<
" " << cHigh <<
G4endl;
512 G4int residualA = targetA+incidentA-
A;
513 G4int residualZ = targetZ+incidentZ-
Z;
516 incidentEnergy, incidentMass,
517 productEnergy, productMass,
518 residualMass, residualA, residualZ,
519 targetMass, targetA, targetZ,
520 incidentA,incidentZ,
A,
Z);
521 cosTh = theKallbach.
Sample(anEnergy);
522 if (std::getenv(
"G4PHPTEST") )
523 G4cout <<
" G4ParticleHPKallbachMannSyst::Sample resulttest " << cosTh <<
G4endl;
526 }
else if (angularRep > 10 && angularRep < 16) {
532 if (i != 0) running[i] = running[i-1];
543 fCache.Get()->currentMeanEnergy = 0.0;
550 if(random<running[i]/running[
nEnergies-1])
break;
566 cosTh = theStore.
Sample();
577 theStore.
SetY(aCounter,
585 cosTh = theStore.
Sample();
614 x = theBuff1.
GetX(i);
615 y1 = theBuff1.
GetY(i);
616 y2 = theBuff2.
GetY(i);
623 cosTh = theStore.
Sample();
629 "G4ParticleHPContAngularPar::Sample: Unknown angular representation");
637 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
719 G4int ie, ie1, ie2, ie1Prev, ie2Prev;
723 if (
fCache.Get()->fresh !=
true)
return;
743 const std::map<G4double,G4int> discEnerOwn1 = copyAngpar1.GetDiscreteEnergiesOwn();
745 std::map<G4double,G4int>::const_iterator itedeo1;
746 std::map<G4double,G4int>::const_iterator itedeo2;
747 std::vector<G4ParticleHPList*> vAngular(discEnerOwn1.size() );
749 for (itedeo1 = discEnerOwn1.begin(); itedeo1 != discEnerOwn1.end(); ++itedeo1) {
750 discEner1 = itedeo1->first;
757 ie1 = itedeo1->second;
758 itedeo2 = discEnerOwn2.find(discEner1);
759 if( itedeo2 == discEnerOwn2.end() ) {
762 ie2 = itedeo2->second;
765 vAngular[ie1]->SetLabel(copyAngpar1.theAngular[ie1].GetLabel());
768 val1 = copyAngpar1.theAngular[ie1].GetValue(ip);
775 copyAngpar1.theEnergy, copyAngpar2.
theEnergy,
777 if (std::getenv(
"G4PHPTEST2") )
778 G4cout << ie1 <<
" " << ip
779 <<
" G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 "
780 << val1 <<
" val2 " << val2 <<
" value " << value <<
G4endl;
782 vAngular[ie1]->SetValue(ip, value);
787 std::vector<G4ParticleHPList*>::iterator itv;
789 for (itedeo2 = discEnerOwn2.begin(); itedeo2 != discEnerOwn2.end(); ++itedeo2) {
790 discEner2 = itedeo2->first;
791 ie2 = itedeo2->second;
793 itedeo1 = discEnerOwn1.find(discEner2);
794 if (itedeo1 != discEnerOwn1.end() ) {
806 G4bool isInserted =
false;
808 for (itv = vAngular.begin(); itv != vAngular.end(); ++itv,++ie) {
809 if (discEner2 > (*itv)->GetLabel() ) {
828 copyAngpar1.theEnergy,
831 if (std::getenv(
"G4PHPTEST2") )
833 <<
" G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 "
834 << val1 <<
" val2 " << val2 <<
" value " << value <<
G4endl;
835 vAngular[ie]->SetValue(ip, value);
866 G4double interMinEner = copyAngpar1.GetMinEner() + (
theEnergy-copyAngpar1.GetEnergy() )
867 * (copyAngpar2.
GetMinEner() - copyAngpar1.GetMinEner() )
868 / (copyAngpar2.
GetEnergy()-copyAngpar1.GetEnergy() );
869 G4double interMaxEner = copyAngpar1.GetMaxEner() + (
theEnergy-copyAngpar1.GetEnergy() )
870 * (copyAngpar2.
GetMaxEner()-copyAngpar1.GetMaxEner() )
871 / (copyAngpar2.
GetEnergy()-copyAngpar1.GetEnergy() );
873 if (std::getenv(
"G4PHPTEST2") )
874 G4cout <<
" G4ParticleHPContAngularPar::Merge E " << anEnergy <<
" minmax "
875 << interMinEner <<
" " << interMaxEner <<
G4endl;
880 G4int nEnergies1 = copyAngpar1.GetNEnergies();
881 G4int nDiscreteEnergies1 = copyAngpar1.GetNDiscreteEnergies();
882 G4double minEner1 = copyAngpar1.GetMinEner();
883 G4double maxEner1 = copyAngpar1.GetMaxEner();
896 for (ie1 = nDiscreteEnergies1; ie1 < nEnergies1; ++ie1) {
897 e1 = copyAngpar1.theAngular[ie1].GetLabel();
898 eTNorm1 = (
e1 - minEner1);
899 if (maxEner1 != minEner1) eTNorm1 /= (maxEner1-minEner1);
905 for (ie2 = nDiscreteEnergies2; ie2 < nEnergies2; ++ie2) {
907 eTNorm2 = (
e2 - minEner2);
908 if (maxEner2 != minEner2) eTNorm2 /= (maxEner2-minEner2);
940 e1Interp = (maxEner1 - minEner1) * eT + minEner1;
942 for (ie1 = nDiscreteEnergies1; ie1 < nEnergies1; ++ie1) {
943 if ((copyAngpar1.theAngular[ie1].GetLabel() - e1Interp) > 1.E-10*e1Interp)
break;
946 if (ie1 == 0) ie1Prev++;
947 if (ie1 == nEnergies1) {
953 e2Interp = (maxEner2-minEner2) * eT + minEner2;
955 for (ie2 = nDiscreteEnergies2; ie2 < nEnergies2; ++ie2) {
959 if (ie2 == 0) ie2Prev++;
960 if (ie2 == nEnergies2) {
966 G4double eN = (interMaxEner-interMinEner) * eT + interMinEner;
967 if (std::getenv(
"G4PHPTEST2") )
968 G4cout << ie <<
" " << ie1 <<
" " << ie2 <<
" G4ParticleHPContAngularPar::loop eT "
969 << eT <<
" -> eN " << eN <<
" e1Interp " << e1Interp <<
" e2Interp " << e2Interp <<
G4endl;
985 copyAngpar1.theAngular[ie1Prev].GetLabel(),
986 copyAngpar1.theAngular[ie1].GetLabel(),
987 copyAngpar1.theAngular[ie1Prev].GetValue(ip),
988 copyAngpar1.theAngular[ie1].GetValue(ip)) * (maxEner1-minEner1);
997 copyAngpar1.theEnergy, copyAngpar2.
theEnergy,
999 if (interMaxEner != interMinEner) {
1000 value /= (interMaxEner-interMinEner);
1001 }
else if (value != 0) {
1003 "G4ParticleHPContAngularPar::PrepareTableInterpolation interMaxEner == interMinEner and value != 0.");
1005 if (std::getenv(
"G4PHPTEST2") )
1006 G4cout << ie <<
" " << ip <<
" G4ParticleHPContAngularPar::Merge val1 " << val1
1007 <<
" val2 " << val2 <<
" value " << value <<
G4endl;
1012 for (itv = vAngular.begin(); itv != vAngular.end(); ++itv)
delete (*itv);
1014 if (std::getenv(
"G4PHPTEST2") ) {
1015 G4cout <<
" G4ParticleHPContAngularPar::Merge ANGPAR1 " <<
G4endl;
1017 G4cout <<
" G4ParticleHPContAngularPar::Merge ANGPAR2 " <<
G4endl;
1019 G4cout <<
" G4ParticleHPContAngularPar::Merge ANGPARNEW " <<
G4endl;
static const G4double e1[44]
static const G4double e2[44]
static constexpr double twopi
static constexpr double eV
G4GLOB_DLL std::ostream G4cout
static G4Deuteron * Deuteron()
static G4Electron * Electron()
void Init(G4int aScheme, G4int aRange)
G4InterpolationScheme GetScheme(G4int index) const
G4InterpolationScheme GetInverseScheme(G4int index)
static G4IonTable * GetIonTable()
static G4Neutron * Neutron()
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetPDGMass() const
G4double GetPDGCharge() const
G4ReactionProduct * Sample(G4double anEnergy, G4double massCode, G4double mass, G4int angularRep, G4int interpol)
void PrepareTableInterpolation()
std::set< G4double > theEnergiesTransformed
G4Cache< toBeCached * > fCache
G4InterpolationManager theManager
G4int GetNEnergies() const
void Init(std::istream &aDataFile, G4ParticleDefinition *projectile)
std::map< G4double, G4int > GetDiscreteEnergiesOwn() const
G4ParticleDefinition * theProjectile
std::set< G4double > theDiscreteEnergies
G4ParticleHPInterpolator theInt
G4int GetNDiscreteEnergies() const
G4double GetEnergy() const
G4double GetMinEner() const
std::map< G4double, G4int > theDiscreteEnergiesOwn
G4ParticleHPContAngularPar()
void BuildByInterpolation(G4double anEnergy, G4InterpolationScheme aScheme, G4ParticleHPContAngularPar &store1, G4ParticleHPContAngularPar &store2)
G4double GetMaxEner() const
G4ParticleHPList * theAngular
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
G4double GetWeightedBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
G4double Interpolate2(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
G4double GetBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
G4double Sample(G4double anEnergy)
void SetCoeff(G4int i, G4int l, G4double coeff)
void Init(G4int i, G4double e, G4int n)
G4double SampleMax(G4double energy)
void SetManager(G4InterpolationManager &aManager)
void SetLabel(G4double aLabel)
void SetValue(G4int i, G4double y)
void Init(std::istream &aDataFile, G4int nPar, G4double unit=1.)
G4double GetValue(G4int i)
static G4ParticleHPManager * GetInstance()
void SetY(G4int i, G4double x)
void SetX(G4int i, G4double e)
void SetInterpolationManager(const G4InterpolationManager &aManager)
G4double GetY(G4double x)
G4double GetX(G4int i) const
G4int GetVectorLength() const
static G4Positron * Positron()
static G4Proton * Proton()
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4double GetTotalMomentum() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetKineticEnergy(const G4double en)
static G4Triton * Triton()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments