66 std::ostringstream ost;
69 std::ifstream from(aName, std::ios::in);
73 std::ifstream theGammaData(aName, std::ios::in);
82 if(!getenv(
"G4NEUTRONHPDATA"))
83 throw G4HadronicException(__FILE__, __LINE__,
"Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
84 G4String tBase = getenv(
"G4NEUTRONHPDATA");
98 if(getenv(
"NeutronHPNamesLogging"))
G4cout <<
"Skipped = "<< filename <<
" "<<A<<
" "<<Z<<
G4endl;
107 std::istringstream theData(std::ios::in);
118 G4int infoType, dataType, dummy;
121 while (theData >> infoType)
125 theData >> sfType >> dummy;
127 if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
135 theData >> dqi >> ilr;
160 else if(dataType==12)
165 else if(dataType==13)
170 else if(dataType==14)
174 else if(dataType==15)
180 throw G4HadronicException(__FILE__, __LINE__,
"Data-type unknown to G4NeutronHPInelasticCompFS");
195 if(i!=0) running[i]=running[i-1];
207 for(i0=0; i0<50; i0++)
210 if(random < running[i0]/sum)
break;
229 theNeutron.SetKineticEnergy( eKinetic );
260 boosted.
Lorentz(theNeutron, theTarget);
274 G4double availableEnergy = theNeutron.GetKineticEnergy() + theNeutron.GetMass() - aHadron.
GetMass() +
277 if ( availableEnergy < 0 )
282 G4int nothingWasKnownOnHadron = 0;
303 if ( p2 > 0.0 ) p = std::sqrt( p );
305 aHadron.
SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*p );
326 if (
QI[it] < 0 || 849 <
QI[it] ) dqi =
QI[it];
331 }
while(eSecN>MaxEne);
334 eGamm = eKinetic-eSecN;
340 iLevel+=
G4int(random);
347 while (eKinetic-eExcitation < 0 && iLevel>0)
356 if ( dqi < 0 || 849 < dqi ) useQI =
true;
361 eExcitation = -
QI[it];
379 if ( iLevel == -1 ) iLevel = 0;
385 if ( !find ) iLevel = imaxEx;
389 if(getenv(
"InelasticCompFSLogging") && eKinetic-eExcitation < 0)
391 throw G4HadronicException(__FILE__, __LINE__,
"SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
393 if(eKinetic-eExcitation < 0) eExcitation = 0;
412 theRestEnergy->
SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
413 eGamm*std::sin(std::acos(costh))*std::sin(phi),
416 thePhotons->push_back(theRestEnergy);
426 if ( theParticles != NULL ) {
431 for (
G4int j = 0 ; j != (
G4int)theParticles->size() ; j++ ) {
432 if ( theParticles->at(j)->GetDefinition()->GetBaryonNumber() > maxA ) {
433 maxA = theParticles->at(j)->GetDefinition()->GetBaryonNumber();
436 sumA += theParticles->at(j)->GetDefinition()->GetBaryonNumber();
437 sumZ +=
G4int( theParticles->at(j)->GetDefinition()->GetPDGCharge() +
eps );
441 if ( dA < 0 || dZ < 0 ) {
442 G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA ;
443 G4int newZ =
G4int( theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge() +
eps ) + dZ;
445 theParticles->at( jAtMaxA )->SetDefinition( pd );
454 nothingWasKnownOnHadron = 1;
468 boosted_tmp.
Lorentz(theNeutron, theTarget);
473 if(thePhotons!=0 && thePhotons->size()!=0)
474 { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
477 while(aBaseEnergy>0.01*
keV)
480 G4bool foundMatchingLevel =
false;
483 for(
G4int j=1; j<it; j++)
493 G4double deltaE = std::abs(testEnergy-aBaseEnergy);
498 thePhotons->push_back(theNext->operator[](0));
499 aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
501 foundMatchingLevel =
true;
510 if(!foundMatchingLevel)
514 thePhotons->push_back(theNext->operator[](0));
515 aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
524 for(i0=0; i0<thePhotons->size(); i0++)
527 thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
531 if(nothingWasKnownOnHadron)
538 if ( thePhotons != 0 )
540 unsigned int nPhotons = thePhotons->size();
542 for ( ii0=0; ii0<nPhotons; ii0++)
545 totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
554 G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
564 G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
613 G4int nSecondaries = 2;
614 G4bool needsSeparateRecoil =
false;
615 G4int totalBaryonNumber = 0;
616 G4int totalCharge = 0;
618 if(theParticles != 0)
620 nSecondaries = theParticles->size();
623 for(ii0=0; ii0<theParticles->size(); ii0++)
625 aDef = theParticles->operator[](ii0)->GetDefinition();
628 totalMomentum += theParticles->operator[](ii0)->GetMomentum();
632 needsSeparateRecoil =
true;
642 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
643 nSecondaries += nPhotons;
647 if( theParticles==0 )
654 aHadron.
Lorentz(aHadron, theTarget);
657 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
662 G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
665 theResidual.
Lorentz(theResidual, -1.*theTarget);
669 for(i=0; i<nPhotons; i++)
671 totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
681 for(i0=0; i0<theParticles->size(); i0++)
684 theSec->
SetDefinition(theParticles->operator[](i0)->GetDefinition());
685 theSec->
SetMomentum(theParticles->operator[](i0)->GetMomentum());
687 delete theParticles->operator[](i0);
690 if(needsSeparateRecoil && residualZ!=0)
694 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
696 resiualKineticEnergy += totalMomentum*totalMomentum;
697 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.
GetMass();
716 for(i=0; i<nPhotons; i++)
721 theSec->
SetDefinition( thePhotons->operator[](i)->GetDefinition() );
723 theSec->
SetMomentum(thePhotons->operator[](i)->GetMomentum());
725 delete thePhotons->operator[](i);
758 rot1.setPhi(
pi/2 + d.phi() );
760 rot2.setTheta( d.theta() );
782 if ( ( 1 + (1+A)/A*Q/E1 ) < 0 )
785 Q = -( A/(1+
A)*E1 ) + 1.0e-6*
eV;
788 G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
790 G4double E3 = AA/std::pow((1+A),2)*(beta*beta+1+2*beta*mu)*E1;
791 G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu);
792 if ( omega3 > 1.0 ) omega3 = 1.0;
801 G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ;
802 G4ThreeVector p ( 0 , pmag*std::sqrt(1-omega3*omega3), pmag*omega3 );
G4ReactionProductVector * GetPhotons(G4double anEnergy)
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetLevelEnergy()
G4NeutronHPNames theNames
void SetMomentum(const G4ThreeVector &momentum)
void two_body_reaction(G4DynamicParticle *, G4DynamicParticle *, G4DynamicParticle *, G4double mu)
G4NeutronHPLevel * GetLevel(G4int i)
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
G4double GetKineticEnergy() const
G4Cache< G4HadFinalState * > theResult
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepRotation G4RotationMatrix
void SetKineticEnergy(const G4double en)
static G4NeutronHPManager * GetInstance()
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4NeutronHPVector * theXsection[51]
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
void Init(std::istream &aDataFile)
G4double Sample(G4double anEnergy, G4int &it)
virtual G4NeutronHPVector * GetXsec()
void Init(std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
static const G4double eps
G4NeutronHPDeExGammas theGammas
G4ParticleDefinition * GetDefinition() const
void GetDataStream(G4String, std::istringstream &iss)
G4NeutronHPDataUsed GetName(G4int A, G4int Z, G4String base, G4String rest, G4bool &active)
std::vector< G4double > QI
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void InitAngular(std::istream &aDataFile)
const G4ParticleDefinition * GetDefinition() const
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
G4GLOB_DLL std::ostream G4cout
G4NeutronHPEnergyDistribution * theEnergyDistribution[51]
const G4ParticleDefinition * GetDefinition() const
void CompositeApply(const G4HadProjectile &theTrack, G4ParticleDefinition *aHadron)
G4double GetKineticEnergy() const
G4int SelectExitChannel(G4double eKinetic)
G4ErrorTarget * theTarget
void InitEnergies(std::istream &aDataFile)
void Init(std::istream &aDataFile)
static G4Neutron * Neutron()
static const G4double A[nN]
const G4LorentzVector & Get4Momentum() const
void InitGammas(G4double AR, G4double ZR)
void SetKineticEnergy(G4double aEnergy)
G4NeutronHPEnAngCorrelation * theEnergyAngData[51]
G4LorentzVector Get4Momentum() const
static G4IonTable * GetIonTable()
G4double GetKineticEnergy() const
void SampleAndUpdate(G4ReactionProduct &aNeutron)
G4NeutronHPAngular * theAngularDistribution[51]
G4int GetNumberOfLevels()
G4double GetTotalEnergy() const
G4bool InitMean(std::istream &aDataFile)
G4double total(Particle const *const p1, Particle const *const p2)
G4double GetPDGMass() const
G4double GetLevelEnergy(G4int aLevel)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4NeutronHPPhotonDist * theFinalStatePhotons[51]
void Init(std::istream &theData)
G4ThreeVector GetMomentum() const
G4double GetTemperature() const
const G4Material * GetMaterial() const
G4ReactionProductVector * Sample(G4double anEnergy)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aSFType)
G4ReactionProductVector * GetDecayGammas(G4int aLevel)
void SetAZMs(G4double anA, G4double aZ, G4int aM, G4NeutronHPDataUsed used)
void InitDistributionInitialState(G4ReactionProduct &aNeutron, G4ReactionProduct &aTarget, G4int it)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void adjust_final_state(G4LorentzVector)
G4double GetPDGCharge() const
void Put(const value_type &val) const
G4double GetLevelEnergy()
void InitPartials(std::istream &aDataFile)
void Init(std::istream &aDataFile)
G4ThreeVector GetMomentum() const
G4int GetBaryonNumber() const
CLHEP::HepLorentzVector G4LorentzVector