106 if( getenv(
"G4PHPTEST") )
G4cout <<
" G4ParticleHPContAngularPar::Sample " << anEnergy <<
" " << massCode <<
" " << angularRep <<
G4endl;
107 if ( fCache.
Get() == NULL ) cacheInit();
137 if(Z!=2)
throw G4HadronicException(__FILE__, __LINE__,
"G4ParticleHPContAngularPar: Unknown ion case 1");
148 if( angularRep == 1 )
154 if ( nDiscreteEnergies != 0 )
159 if ( fCache.
Get()->fresh == true )
163 fCache.
Get()->remaining_energy =
std::max ( theAngular[0].GetLabel() , theAngular[nEnergies-1].GetLabel() );
164 fCache.
Get()->fresh =
false;
169 if ( nDiscreteEnergies == nEnergies )
171 fCache.
Get()->remaining_energy =
std::max ( fCache.
Get()->remaining_energy , theAngular[nDiscreteEnergies-1].
GetLabel() );
178 for (
G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
180 cont_min = theAngular[j].
GetLabel();
181 if ( theAngular[j].GetValue(0) != 0.0 )
break;
183 fCache.
Get()->remaining_energy =
std::max ( fCache.
Get()->remaining_energy ,
std::min ( theAngular[nDiscreteEnergies-1].GetLabel() , cont_min ) );
191 for (
G4int j = 0 ; j < nDiscreteEnergies ; j++ )
194 if ( theAngular[j].GetLabel() <= fCache.
Get()->remaining_energy ) delta = theAngular[i].GetValue(0);
195 running[j+1] = running[j] + delta;
197 G4double tot_prob_DIS = running[ nDiscreteEnergies ];
199 for (
G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
204 if ( theAngular[j].GetLabel() <= fCache.
Get()->remaining_energy ) delta = theAngular[j].GetValue(0);
211 if ( theAngular[j].GetLabel() != 0 )
213 if ( j == nDiscreteEnergies ) {
222 if ( theAngular[j].GetLabel() == 0.0 ) {
224 if ( j != nEnergies-1 ) {
228 if ( theAngular[j].GetValue(0) != 0.0 ) {
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 );
236 G4double tot_prob_CON = running[ nEnergies ] - running[ nDiscreteEnergies ];
250 random *= (tot_prob_DIS + tot_prob_CON);
252 if ( random <= ( tot_prob_DIS / ( tot_prob_DIS + tot_prob_CON ) ) || nDiscreteEnergies == nEnergies )
255 for (
G4int j = 0 ; j < nDiscreteEnergies ; j++ )
258 if ( random < running[ j+1 ] )
264 fsEnergy = theAngular[ it ].
GetLabel();
267 theStore.Init(0,fsEnergy,nAngularParameters);
268 for (
G4int j=0;j<nAngularParameters;j++)
270 theStore.SetCoeff(0,j,theAngular[it].GetValue(j));
273 cosTh = theStore.SampleMax(fsEnergy);
279 for (
G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
282 if ( random < running[ j ] )
293 if ( it != nDiscreteEnergies )
301 theStore.Init(0,y1,nAngularParameters);
302 theStore.Init(1,y2,nAngularParameters);
303 theStore.SetManager(theManager);
304 for (
G4int j=0;j<nAngularParameters;j++)
307 if ( it == nDiscreteEnergies ) itt = it+1;
314 theStore.SetCoeff(0,j,theAngular[itt-1].GetValue(j));
315 theStore.SetCoeff(1,j,theAngular[itt].GetValue(j));
318 cosTh = theStore.SampleMax(fsEnergy);
324 if( adjustResult ) fCache.
Get()->remaining_energy -= fsEnergy;
336 if ( fCache.
Get()->fresh == true )
338 fCache.
Get()->remaining_energy = theAngular[ nEnergies-1 ].
GetLabel();
339 fCache.
Get()->fresh =
false;
346 for(i=1; i<nEnergies; i++)
361 running[i]=running[i-1];
362 if ( fCache.
Get()->remaining_energy >= theAngular[i].
GetLabel() )
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 )
387 for ( i = 1 ; i < nEnergies ; i++ )
390 if ( random < running [ i ] / running [ nEnergies-1 ] )
break;
395 if ( running [ nEnergies-1 ] == 0 ) it = 0;
398 if (it<nDiscreteEnergies||it==0)
402 fsEnergy = theAngular[it].
GetLabel();
404 theStore.Init(0,fsEnergy,nAngularParameters);
405 for(i=0;i<nAngularParameters;i++)
407 theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
410 cosTh = theStore.SampleMax(fsEnergy);
419 running[it-1]/running[nEnergies-1],
420 running[it]/running[nEnergies-1],
424 theStore.Init(0,e1,nAngularParameters);
425 theStore.Init(1,e2,nAngularParameters);
426 for(i=0;i<nAngularParameters;i++)
428 theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
429 theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
432 theStore.SetManager(theManager);
433 cosTh = theStore.SampleMax(fsEnergy);
438 G4double x1 = running[it-1]/running[nEnergies-1];
439 G4double x2 = running[it]/running[nEnergies-1];
445 theStore.Init(0,y1,nAngularParameters);
446 theStore.Init(1,y2,nAngularParameters);
447 theStore.SetManager(theManager);
448 for(i=0;i<nAngularParameters;i++)
450 theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
451 theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
454 cosTh = theStore.SampleMax(fsEnergy);
459 if( adjustResult ) fCache.
Get()->remaining_energy -= fsEnergy;
463 else if(angularRep==2)
470 if( getenv(
"G4PHPTEST") )
G4cout <<
" G4ParticleHPContAngularPar::Sample nEnergies " << nEnergies <<
G4endl;
471 for(j=1; j<nEnergies; j++)
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];
495 for(j=1; j<nEnergies; j++)
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;
528 G4int targetA =
G4int(fCache.
Get()->theTargetCode-1000*targetZ);
531 targetA =
G4int ( fCache.
Get()->theTarget->GetMass()/
amu_c2 + 0.5 );
532 G4double targetMass = fCache.
Get()->theTarget->GetMass();
533 G4int residualA = targetA+1-
A;
534 G4int residualZ = targetZ-
Z;
536 residualMass +=(residualA-residualZ)*theProjectile->
GetPDGMass();
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)
552 for(i=1; i<nEnergies; i++)
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;
572 for(i=1; i<nEnergies; i++)
575 if(random<running[i]/running[nEnergies-1])
break;
577 if(it<nDiscreteEnergies||it==0)
581 fsEnergy = theAngular[0].
GetLabel();
584 for(
G4int j=1; j<nAngularParameters; j+=2)
586 theStore.
SetX(aCounter, theAngular[0].GetValue(j));
587 theStore.
SetY(aCounter, theAngular[0].GetValue(j+1));
591 aMan.
Init(angularRep-10, nAngularParameters-1);
593 cosTh = theStore.
Sample();
597 fsEnergy = theAngular[it].
GetLabel();
600 aMan.
Init(angularRep-10, nAngularParameters-1);
604 for(
G4int j=1; j<nAngularParameters; j+=2)
606 theStore.
SetX(aCounter, theAngular[it].GetValue(j));
609 running[it-1]/running[nEnergies-1],
610 running[it]/running[nEnergies-1],
611 theAngular[it-1].GetValue(j+1),
615 cosTh = theStore.
Sample();
620 G4double x1 = running[it-1]/running[nEnergies-1];
621 G4double x2 = running[it]/running[nEnergies-1];
629 aMan.
Init(angularRep-10, nAngularParameters-1);
643 for(i=0,j=1; i<nAngularParameters; i++,j+=2)
645 theBuff1.
SetX(i, theAngular[it-1].GetValue(j));
646 theBuff1.
SetY(i, theAngular[it-1].GetValue(j+1));
647 theBuff2.
SetX(i, theAngular[it].GetValue(j));
648 theBuff2.
SetY(i, theAngular[it].GetValue(j+1));
659 x = theBuff1.
GetX(i);
660 y1 = theBuff1.
GetY(i);
661 y2 = theBuff2.
GetY(i);
663 fsEnergy, theAngular[it-1].
GetLabel(),
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) );
G4double G4ParticleHPJENDLHEData::G4double result
G4int GetVectorLength() const
G4double GetTotalMomentum() const
void Init(G4int aScheme, G4int aRange)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetInterpolationManager(const G4InterpolationManager &aManager)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static constexpr double twopi
void SetY(G4int i, G4double x)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4InterpolationScheme GetScheme(G4int index) const
static G4Triton * Triton()
static G4Proton * Proton()
static constexpr double eV
static constexpr double amu_c2
G4double GetX(G4int i) const
static G4Neutron * Neutron()
G4double GetY(G4double x)
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
static G4IonTable * GetIonTable()
static G4Positron * Positron()
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
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 SetX(G4int i, G4double e)