64 fCache.Get()->currentMeanEnergy = -2;
65 fCache.Get()->fresh =
true;
67 if ( getenv(
"G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) )
adjustResult =
false;
82 if ( getenv(
"G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) )
adjustResult =
false;
106 if( getenv(
"G4PHPTEST") )
G4cout <<
" G4ParticleHPContAngularPar::Sample " << anEnergy <<
" " << massCode <<
" " << angularRep <<
G4endl;
137 if(Z!=2)
throw G4HadronicException(__FILE__, __LINE__,
"G4ParticleHPContAngularPar: Unknown ion case 1");
148 if( angularRep == 1 )
159 if (
fCache.Get()->fresh == true )
164 fCache.Get()->fresh =
false;
181 if (
theAngular[j].GetValue(0) != 0.0 )
break;
195 running[j+1] = running[j] + delta;
213 if ( j == nDiscreteEnergies ) {
224 if ( j != nEnergies-1 ) {
229 throw G4HadronicException(__FILE__, __LINE__,
"G4ParticleHPContAngularPar: Unexpected non zero value of theAngular[nEnergies-1].GetValue(0)");
234 running[j+1] = running[j] + ( ( e_high - e_low ) * delta );
250 random *= (tot_prob_DIS + tot_prob_CON);
252 if ( random <= ( tot_prob_DIS / ( tot_prob_DIS + tot_prob_CON ) ) || nDiscreteEnergies ==
nEnergies )
258 if ( random < running[ j+1 ] )
282 if ( random < running[ j ] )
293 if ( it != nDiscreteEnergies )
307 if ( it == nDiscreteEnergies ) itt = it+1;
336 if (
fCache.Get()->fresh == true )
339 fCache.Get()->fresh =
false;
361 running[i]=running[i-1];
374 if ( nEnergies == 1 || running[nEnergies-1] == 0 )
375 fCache.Get()->currentMeanEnergy = 0.0;
378 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
382 if ( nEnergies == 1 ) it = 0;
385 if ( running[nEnergies-1] != 0 )
390 if ( random < running [ i ] / running [ nEnergies-1 ] )
break;
395 if ( running [ nEnergies-1 ] == 0 ) it = 0;
419 running[it-1]/running[nEnergies-1],
420 running[it]/running[nEnergies-1],
438 G4double x1 = running[it-1]/running[nEnergies-1];
439 G4double x2 = running[it]/running[nEnergies-1];
463 else if(angularRep==2)
470 if( getenv(
"G4PHPTEST") )
G4cout <<
" G4ParticleHPContAngularPar::Sample nEnergies " <<
nEnergies <<
G4endl;
473 if(j!=0) running[j]=running[j-1];
480 if( getenv(
"G4PHPTEST") )
G4cout <<
" G4ParticleHPContAngularPar::Sample " << j <<
" running " << running[j]
486 if ( nEnergies == 1 )
487 fCache.Get()->currentMeanEnergy = 0.0;
489 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
498 if(randkal<running[j]/running[nEnergies-1])
break;
504 x = randkal*running[nEnergies-1];
513 if( getenv(
"G4PHPTEST") )
G4cout << itt <<
" G4particleHPContAngularPar fsEnergy " << fsEnergy <<
" " <<
theManager.
GetInverseScheme(itt-1) <<
" x " << x <<
" " << x1 <<
" " << x2 <<
" y " << y1 <<
" " << y2 <<
G4endl;
518 fsEnergy, y1, y2, cLow,cHigh);
519 if( getenv(
"G4PHPTEST") )
G4cout << itt <<
" G4particleHPContAngularPar compoundFraction " << compoundFraction <<
" E " << fsEnergy <<
" " <<
theManager.
GetScheme(itt) <<
" ener " << fsEnergy <<
" y " << y1 <<
" " << y2 <<
" cLH " << cLow <<
" " << cHigh <<
G4endl;
531 targetA =
G4int (
fCache.Get()->theTarget->GetMass()/amu_c2 + 0.5 );
533 G4int residualA = targetA+1-
A;
534 G4int residualZ = targetZ-Z;
539 incidentEnergy, incidentMass,
540 productEnergy, productMass,
541 residualMass, residualA, residualZ,
542 targetMass, targetA, targetZ);
543 cosTh = theKallbach.
Sample(anEnergy);
544 if( getenv(
"G4PHPTEST") )
G4cout <<
" G4ParticleHPKallbachMannSyst::Sample resulttest " << cosTh <<
G4endl;
546 else if(angularRep>10&&angularRep<16)
554 if(i!=0) running[i]=running[i-1];
564 if ( nEnergies == 1 )
565 fCache.Get()->currentMeanEnergy = 0.0;
567 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
570 if ( nEnergies == 1 ) it = 0;
575 if(random<running[i]/running[nEnergies-1])
break;
591 aMan.
Init(angularRep-10, nAngularParameters-1);
593 cosTh = theStore.
Sample();
609 running[it-1]/running[nEnergies-1],
610 running[it]/running[nEnergies-1],
615 cosTh = theStore.
Sample();
620 G4double x1 = running[it-1]/running[nEnergies-1];
621 G4double x2 = running[it]/running[nEnergies-1];
659 x = theBuff1.
GetX(i);
660 y1 = theBuff1.
GetY(i);
661 y2 = theBuff2.
GetY(i);
668 cosTh = theStore.
Sample();
674 throw G4HadronicException(__FILE__, __LINE__,
"G4ParticleHPContAngularPar::Sample: Unknown angular representation");
681 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
698 if( !angParPrev )
return;
705 for(ie=0; ie<nDiscreteEnergiesPrev; ie++) {
710 for(ie=nDiscreteEnergies; ie<
nEnergies; ie++) {
719 for(ie=nDiscreteEnergiesPrev; ie<nEnergiesPrev; ie++) {
721 G4double enerT = (ener-minEnerPrev)/(maxEnerPrev-minEnerPrev);
734 G4int ie,ie1,ie2, ie1Prev, ie2Prev;
740 std::set<G4double>::const_iterator itede;
743 std::map<G4double,G4int>::const_iterator itedeo;
747 itedeo = discEnerOwn1.find(discEner);
748 if( itedeo == discEnerOwn1.end() ) {
753 itedeo = discEnerOwn2.find(discEner);
754 if( itedeo == discEnerOwn2.end() ) {
778 if( getenv(
"G4PHPTEST2") )
G4cout << ie <<
" " << ip <<
" G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 <<
" val2 " << val2 <<
" value " << value <<
G4endl;
796 std::set<G4double>::const_iterator iteet = energiesTransformed.begin();
809 G4double e1 = (maxEner1-minEner1) * eT + minEner1;
811 for(ie1=nDiscreteEnergies1; ie1<nEnergies1; ie1++) {
815 if( ie1 == 0 ) ie1Prev++;
816 if( ie1 == nEnergies1 ) {
821 G4double e2 = (maxEner2-minEner2) * eT + minEner2;
823 for(ie2=nDiscreteEnergies2; ie2<nEnergies2; ie2++) {
828 if( ie2 == 0 ) ie2Prev++;
829 if( ie2 == nEnergies2 ) {
836 if( getenv(
"G4PHPTEST2") )
G4cout << ie <<
" " << ie1 <<
" " << ie2 <<
" G4ParticleHPContAngularPar::loop eT " << eT <<
" -> eN " << eN <<
" e1 " << e1 <<
" e2 " << e2 <<
G4endl;
861 }
else if ( value != 0 ) {
862 throw G4HadronicException(__FILE__, __LINE__,
"G4ParticleHPContAngularPar::PrepareTableInterpolation theMaxEner == theMinEner and value != 0.");
864 if( getenv(
"G4PHPTEST2") )
G4cout << ie <<
" " << ip <<
" G4ParticleHPContAngularPar::Merge val1 " << val1 <<
" val2 " << val2 <<
" value " << value <<
G4endl;
873 if( getenv(
"G4PHPTEST2") ) {
874 G4cout <<
" G4ParticleHPContAngularPar::Merge ANGPAR1 " <<
G4endl;
876 G4cout <<
" G4ParticleHPContAngularPar::Merge ANGPAR2 " <<
G4endl;
878 G4cout <<
" G4ParticleHPContAngularPar::Merge ANGPARNEW " <<
G4endl;
std::map< G4double, G4int > GetDiscreteEnergiesOwn() const
void Init(G4int i, G4double e, G4int n)
G4int GetVectorLength() const
G4double GetTotalMomentum() const
CLHEP::Hep3Vector G4ThreeVector
void Init(std::istream &aDataFile, G4int nPar, G4double unit=1.)
G4ReactionProduct * Sample(G4double anEnergy, G4double massCode, G4double mass, G4int angularRep, G4int interpol)
void BuildByInterpolation(G4double anEnergy, G4InterpolationScheme aScheme, G4ParticleHPContAngularPar &store1, G4ParticleHPContAngularPar &store2)
G4double Interpolate2(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
void Init(G4int aScheme, G4int aRange)
void SetKineticEnergy(const G4double en)
G4double GetMinEner() const
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4double SampleMax(G4double energy)
G4Cache< toBeCached * > fCache
G4ParticleDefinition * theProjectile
G4InterpolationManager theManager
void SetInterpolationManager(const G4InterpolationManager &aManager)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4ParticleHPContAngularPar()
static constexpr double twopi
std::map< G4double, G4int > theDiscreteEnergiesOwn
std::set< G4double > GetEnergiesTransformed() const
void SetY(G4int i, G4double x)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
void Init(std::istream &aDataFile, G4ParticleDefinition *projectile)
G4int GetNEnergiesTransformed() const
void SetValue(G4int i, G4double y)
G4InterpolationScheme GetScheme(G4int index) const
static G4Triton * Triton()
static G4Proton * Proton()
static constexpr double eV
G4double GetX(G4int i) const
static G4Neutron * Neutron()
G4double GetY(G4double x)
G4double GetMaxEner() const
G4double GetWeightedBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
static G4Deuteron * Deuteron()
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
void SetLabel(G4double aLabel)
static G4IonTable * GetIonTable()
void SetManager(G4InterpolationManager &aManager)
static G4Positron * Positron()
G4int GetNDiscreteEnergies() const
G4ParticleHPInterpolator theInt
G4double GetPDGMass() const
static G4double GetBindingEnergy(const G4int A, const G4int Z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ParticleHPList * theAngular
void SetCoeff(G4int i, G4int l, G4double coeff)
std::set< G4double > theDiscreteEnergies
std::set< G4double > theEnergiesTransformed
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4InterpolationScheme GetInverseScheme(G4int index)
static G4Electron * Electron()
G4double GetValue(G4int i)
G4double GetBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
void PrepareTableInterpolation(const G4ParticleHPContAngularPar *angularPrev)
void SetX(G4int i, G4double e)
G4double Sample(G4double anEnergy)
G4int GetNEnergies() const