67 std::ostringstream ost;
70 std::ifstream from(aName, std::ios::in);
74 std::ifstream theGammaData(aName, std::ios::in);
84 if(!getenv(
"G4NEUTRONHPDATA"))
85 throw G4HadronicException(__FILE__, __LINE__,
"Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files where Inelastic/Gammas data is found.");
86 G4String tBase = getenv(
"G4NEUTRONHPDATA");
93 if( getenv(
"G4ParticleHPDebug") )
G4cout <<
" G4ParticleHPInelasticCompFS::Init FILE " << filename <<
G4endl;
105 if(getenv(
"G4ParticleHPDebug_NamesLogging"))
G4cout <<
"Skipped = "<< filename <<
" "<<A<<
" "<<Z<<
G4endl;
115 std::istringstream theData(std::ios::in);
126 G4int infoType, dataType, dummy;
129 while (theData >> infoType)
133 theData >> sfType >> dummy;
135 if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
143 theData >> dqi >> ilr;
169 else if(dataType==12)
174 else if(dataType==13)
179 else if(dataType==14)
183 else if(dataType==15)
189 throw G4HadronicException(__FILE__, __LINE__,
"Data-type unknown to G4ParticleHPInelasticCompFS");
204 if(i!=0) running[i]=running[i-1];
216 for(i0=0; i0<50; i0++)
220 if(random < running[i0]/sum)
break;
278 boosted.
Lorentz(incidReactionProduct, theTarget);
293 (targetMass - residualMass);
295 if ( availableEnergy < 0 )
300 G4int nothingWasKnownOnHadron = 0;
312 (aHadron.
GetMass()+residualMass));
321 if ( p2 > 0.0 ) p = std::sqrt( p );
344 if (
QI[it] < 0 || 849 <
QI[it] ) dqi =
QI[it];
349 G4int icounter_max=1024;
352 if ( icounter > icounter_max ) {
353 G4cout <<
"Loop-counter exceeded the threshold value at " << __LINE__ <<
"th line of " << __FILE__ <<
"." <<
G4endl;
357 }
while(eSecN>MaxEne);
360 eGamm = eKinetic-eSecN;
366 iLevel+=
G4int(random);
373 while (eKinetic-eExcitation < 0 && iLevel>0)
382 if ( dqi < 0 || 849 < dqi ) useQI =
true;
387 eExcitation = -
QI[it];
405 if ( iLevel == -1 ) iLevel = 0;
411 if ( !find ) iLevel = imaxEx;
415 if(getenv(
"G4ParticleHPDebug") && eKinetic-eExcitation < 0)
417 throw G4HadronicException(__FILE__, __LINE__,
"SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
419 if(eKinetic-eExcitation < 0) eExcitation = 0;
438 theRestEnergy->
SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
439 eGamm*std::sin(std::acos(costh))*std::sin(phi),
442 thePhotons->push_back(theRestEnergy);
453 if ( theParticles != NULL ) {
458 for (
G4int j = 0 ; j != (
G4int)theParticles->size() ; j++ ) {
459 if ( theParticles->at(j)->GetDefinition()->GetBaryonNumber() > maxA ) {
460 maxA = theParticles->at(j)->GetDefinition()->GetBaryonNumber();
463 sumA += theParticles->at(j)->GetDefinition()->GetBaryonNumber();
464 sumZ +=
G4int( theParticles->at(j)->GetDefinition()->GetPDGCharge() +
eps );
468 if ( dA < 0 || dZ < 0 ) {
469 G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA ;
470 G4int newZ =
G4int( theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge() +
eps ) + dZ;
472 theParticles->at( jAtMaxA )->SetDefinition( pd );
481 nothingWasKnownOnHadron = 1;
495 boosted_tmp.
Lorentz(incidReactionProduct, theTarget);
500 if(thePhotons!=0 && thePhotons->size()!=0)
501 { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
507 G4bool foundMatchingLevel =
false;
510 for(
G4int j=1; j<it; j++)
520 G4double deltaE = std::abs(testEnergy-aBaseEnergy);
525 if ( thePhotons != NULL ) thePhotons->push_back(theNext->operator[](0));
526 aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
528 foundMatchingLevel =
true;
537 if(!foundMatchingLevel)
541 if ( thePhotons != NULL ) thePhotons->push_back(theNext->operator[](0));
542 aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
551 for(i0=0; i0<thePhotons->size(); i0++)
554 thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
558 if(nothingWasKnownOnHadron)
565 if ( thePhotons != 0 )
567 unsigned int nPhotons = thePhotons->size();
569 for ( ii0=0; ii0<nPhotons; ii0++)
572 totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
581 G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
588 two_body_reaction ( proj , targ , hadron , mu );
640 G4int nSecondaries = 2;
641 G4bool needsSeparateRecoil =
false;
642 G4int totalBaryonNumber = 0;
643 G4int totalCharge = 0;
645 if(theParticles != 0)
647 nSecondaries = theParticles->size();
650 for(ii0=0; ii0<theParticles->size(); ii0++)
652 aDef = theParticles->operator[](ii0)->GetDefinition();
655 totalMomentum += theParticles->operator[](ii0)->GetMomentum();
659 needsSeparateRecoil =
true;
669 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
670 nSecondaries += nPhotons;
674 if( theParticles==0 )
684 aHadron.
Lorentz(aHadron, theTarget);
687 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
695 theResidual.
Lorentz(theResidual, -1.*theTarget);
699 for(i=0; i<nPhotons; i++)
701 totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
714 for(i0=0; i0<theParticles->size(); i0++)
717 theSec->
SetDefinition(theParticles->operator[](i0)->GetDefinition());
718 theSec->
SetMomentum(theParticles->operator[](i0)->GetMomentum());
723 delete theParticles->operator[](i0);
726 if(needsSeparateRecoil && residualZ!=0)
730 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
732 resiualKineticEnergy += totalMomentum*totalMomentum;
733 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.
GetMass();
756 for(i=0; i<nPhotons; i++)
761 theSec->
SetDefinition( thePhotons->operator[](i)->GetDefinition() );
763 theSec->
SetMomentum(thePhotons->operator[](i)->GetMomentum());
769 delete thePhotons->operator[](i);
826 if ( ( 1 + (1+A)/A*Q/E1 ) < 0 )
832 G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
835 G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu);
836 if ( omega3 > 1.0 ) omega3 = 1.0;
839 G4double omega4 = (1-gamma*mu)/std::sqrt(gamma*gamma+1-2*gamma*mu);
840 if ( omega4 > 1.0 ) omega4 = 1.0;
845 G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ;
846 G4ThreeVector p ( 0 , pmag*std::sqrt(1-omega3*omega3), pmag*omega3 );
849 G4double pmag4 = std::sqrt ((E4+M4)*(E4+M4)-M4*M4) ;
850 G4ThreeVector p4 ( 0 , -pmag4*std::sqrt(1-omega4*omega4), pmag4*omega4 );
static G4ParticleHPManager * GetInstance()
static G4Pow * GetInstance()
G4ReactionProductVector * GetDecayGammas(G4int aLevel)
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4ParticleHPEnergyDistribution * theEnergyDistribution[51]
G4Cache< G4HadFinalState * > theResult
void SetMomentum(const G4ThreeVector &momentum)
G4double GetTotalMomentum() const
G4double powN(G4double x, G4int n) const
G4double GetLevelEnergy()
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
static constexpr double keV
G4ParticleHPPhotonDist * theFinalStatePhotons[51]
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
void SetAZMs(G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)
void GetDataStream(G4String, std::istringstream &iss)
static const G4double eps
G4ParticleDefinition * GetDefinition() const
virtual G4ParticleHPVector * GetXsec()
void Init(std::istream &theData)
std::vector< G4double > QI
G4double Sample(G4double anEnergy, G4int &it)
G4ReactionProductVector * Sample(G4double anEnergy)
const G4String & GetParticleName() const
void InitPartials(std::istream &aDataFile)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
HepRotation inverse() const
G4ParticleDefinition * theProjectile
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void adjust_final_state(G4LorentzVector)
void Init(std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
G4ParticleHPDeExGammas theGammas
const G4ParticleDefinition * GetDefinition() const
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
G4GLOB_DLL std::ostream G4cout
G4ParticleHPEnAngCorrelation * theEnergyAngData[51]
double A(double temperature)
void SampleAndUpdate(G4ReactionProduct &anIncidentParticle)
const G4ParticleDefinition * GetDefinition() const
G4bool InitMean(std::istream &aDataFile)
G4int GetNumberOfLevels()
G4ParticleHPNames theNames
HepLorentzVector & boost(double, double, double)
G4double GetKineticEnergy() const
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aSFType, G4ParticleDefinition *)
void setTheta(double theta)
G4ParticleHPLevel * GetLevel(G4int i)
static G4Neutron * Neutron()
void Init(std::istream &aDataFile)
const G4ParticleDefinition * GetParticleDefinition() const
void InitDistributionInitialState(G4ReactionProduct &anIncidentPart, G4ReactionProduct &aTarget, G4int it)
const G4LorentzVector & Get4Momentum() const
void SetKineticEnergy(G4double aEnergy)
G4LorentzVector Get4Momentum() const
static G4IonTable * GetIonTable()
static constexpr double eV
G4double GetKineticEnergy() const
G4ParticleHPVector * theXsection[51]
G4double GetTotalEnergy() const
G4double total(Particle const *const p1, Particle const *const p2)
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void CompositeApply(const G4HadProjectile &theTrack, G4ParticleDefinition *aHadron)
G4ReactionProductVector * GetPhotons(G4double anEnergy)
G4ThreeVector GetMomentum() const
G4double GetTemperature() const
G4ParticleHPDataUsed GetName(G4int A, G4int Z, G4String base, G4String rest, G4bool &active)
const G4Material * GetMaterial() const
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
void InitGammas(G4double AR, G4double ZR)
void InitEnergies(std::istream &aDataFile)
G4ParticleHPAngular * theAngularDistribution[51]
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4int SelectExitChannel(G4double eKinetic)
void Init(std::istream &aDataFile)
void InitAngular(std::istream &aDataFile)
G4double GetPDGCharge() const
void Put(const value_type &val) const
G4double GetLevelEnergy(G4int aLevel)
static constexpr double twopi
G4int GetNumberOfSecondaries() const
static constexpr double pi
G4ThreeVector GetMomentum() const
void Init(std::istream &aDataFile)
G4int GetBaryonNumber() const
G4double GetLevelEnergy()