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;
275 boosted.
Lorentz(incidReactionProduct, theTarget);
292 if ( availableEnergy < 0 )
297 G4int nothingWasKnownOnHadron = 0;
318 if ( p2 > 0.0 ) p = std::sqrt( p );
341 if (
QI[it] < 0 || 849 <
QI[it] ) dqi =
QI[it];
346 G4int icounter_max=1024;
349 if ( icounter > icounter_max ) {
350 G4cout <<
"Loop-counter exceeded the threshold value at " << __LINE__ <<
"th line of " << __FILE__ <<
"." <<
G4endl;
354 }
while(eSecN>MaxEne);
357 eGamm = eKinetic-eSecN;
363 iLevel+=
G4int(random);
370 while (eKinetic-eExcitation < 0 && iLevel>0)
379 if ( dqi < 0 || 849 < dqi ) useQI =
true;
384 eExcitation = -
QI[it];
402 if ( iLevel == -1 ) iLevel = 0;
408 if ( !find ) iLevel = imaxEx;
412 if(getenv(
"G4ParticleHPDebug") && eKinetic-eExcitation < 0)
414 throw G4HadronicException(__FILE__, __LINE__,
"SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
416 if(eKinetic-eExcitation < 0) eExcitation = 0;
435 theRestEnergy->
SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
436 eGamm*std::sin(std::acos(costh))*std::sin(phi),
439 thePhotons->push_back(theRestEnergy);
450 if ( theParticles != NULL ) {
455 for (
G4int j = 0 ; j != (
G4int)theParticles->size() ; j++ ) {
456 if ( theParticles->at(j)->GetDefinition()->GetBaryonNumber() > maxA ) {
457 maxA = theParticles->at(j)->GetDefinition()->GetBaryonNumber();
460 sumA += theParticles->at(j)->GetDefinition()->GetBaryonNumber();
461 sumZ +=
G4int( theParticles->at(j)->GetDefinition()->GetPDGCharge() +
eps );
465 if ( dA < 0 || dZ < 0 ) {
466 G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA ;
467 G4int newZ =
G4int( theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge() +
eps ) + dZ;
469 theParticles->at( jAtMaxA )->SetDefinition( pd );
478 nothingWasKnownOnHadron = 1;
492 boosted_tmp.
Lorentz(incidReactionProduct, theTarget);
497 if(thePhotons!=0 && thePhotons->size()!=0)
498 { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
504 G4bool foundMatchingLevel =
false;
507 for(
G4int j=1; j<it; j++)
517 G4double deltaE = std::abs(testEnergy-aBaseEnergy);
522 thePhotons->push_back(theNext->operator[](0));
523 aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
525 foundMatchingLevel =
true;
534 if(!foundMatchingLevel)
538 thePhotons->push_back(theNext->operator[](0));
539 aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
548 for(i0=0; i0<thePhotons->size(); i0++)
551 thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
555 if(nothingWasKnownOnHadron)
562 if ( thePhotons != 0 )
564 unsigned int nPhotons = thePhotons->size();
566 for ( ii0=0; ii0<nPhotons; ii0++)
569 totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
578 G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
588 G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
637 G4int nSecondaries = 2;
638 G4bool needsSeparateRecoil =
false;
639 G4int totalBaryonNumber = 0;
640 G4int totalCharge = 0;
642 if(theParticles != 0)
644 nSecondaries = theParticles->size();
647 for(ii0=0; ii0<theParticles->size(); ii0++)
649 aDef = theParticles->operator[](ii0)->GetDefinition();
652 totalMomentum += theParticles->operator[](ii0)->GetMomentum();
656 needsSeparateRecoil =
true;
666 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
667 nSecondaries += nPhotons;
671 if( theParticles==0 )
681 aHadron.
Lorentz(aHadron, theTarget);
684 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
692 theResidual.
Lorentz(theResidual, -1.*theTarget);
696 for(i=0; i<nPhotons; i++)
698 totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
711 for(i0=0; i0<theParticles->size(); i0++)
714 theSec->
SetDefinition(theParticles->operator[](i0)->GetDefinition());
715 theSec->
SetMomentum(theParticles->operator[](i0)->GetMomentum());
720 delete theParticles->operator[](i0);
723 if(needsSeparateRecoil && residualZ!=0)
727 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
729 resiualKineticEnergy += totalMomentum*totalMomentum;
730 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.
GetMass();
753 for(i=0; i<nPhotons; i++)
758 theSec->
SetDefinition( thePhotons->operator[](i)->GetDefinition() );
760 theSec->
SetMomentum(thePhotons->operator[](i)->GetMomentum());
766 delete thePhotons->operator[](i);
801 rot2.setTheta( d.theta() );
823 if ( ( 1 + (1+A)/A*Q/E1 ) < 0 )
829 G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
832 G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu);
833 if ( omega3 > 1.0 ) omega3 = 1.0;
836 G4double omega4 = (1-gamma*mu)/std::sqrt(gamma*gamma+1-2*gamma*mu);
837 if ( omega4 > 1.0 ) omega4 = 1.0;
842 G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ;
843 G4ThreeVector p ( 0 , pmag*std::sqrt(1-omega3*omega3), pmag*omega3 );
846 G4double pmag4 = std::sqrt ((E4+M4)*(E4+M4)-M4*M4) ;
847 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 two_body_reaction(G4DynamicParticle *, G4DynamicParticle *, G4DynamicParticle *, G4double mu)
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepRotation G4RotationMatrix
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)
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
G4double GetKineticEnergy() const
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aSFType, G4ParticleDefinition *)
static const double twopi
G4ErrorTarget * theTarget
G4ParticleHPLevel * GetLevel(G4int i)
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()
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)
G4int GetNumberOfSecondaries() const
G4ThreeVector GetMomentum() const
void Init(std::istream &aDataFile)
G4int GetBaryonNumber() const
G4double GetLevelEnergy()
CLHEP::HepLorentzVector G4LorentzVector