66 std::ostringstream ost;
69 std::ifstream from(aName, std::ios::in);
73 std::ifstream theGammaData(aName, std::ios::in);
83 if(!getenv(
"G4NEUTRONHPDATA"))
84 throw G4HadronicException(__FILE__, __LINE__,
"Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files where Inelastic/Gammas data is found.");
85 G4String tBase = getenv(
"G4NEUTRONHPDATA");
92 if( getenv(
"G4ParticleHPDebug") )
G4cout <<
" G4ParticleHPInelasticCompFS::Init FILE " << filename <<
G4endl;
104 if(getenv(
"G4ParticleHPDebug_NamesLogging"))
G4cout <<
"Skipped = "<< filename <<
" "<<A<<
" "<<Z<<
G4endl;
114 std::istringstream theData(std::ios::in);
125 G4int infoType, dataType, dummy;
128 while (theData >> infoType)
132 theData >> sfType >> dummy;
134 if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
142 theData >> dqi >> ilr;
168 else if(dataType==12)
173 else if(dataType==13)
178 else if(dataType==14)
182 else if(dataType==15)
188 throw G4HadronicException(__FILE__, __LINE__,
"Data-type unknown to G4ParticleHPInelasticCompFS");
203 if(i!=0) running[i]=running[i-1];
215 for(i0=0; i0<50; i0++)
219 if(random < running[i0]/sum)
break;
273 boosted.
Lorentz(incidReactionProduct, theTarget);
290 if ( availableEnergy < 0 )
295 G4int nothingWasKnownOnHadron = 0;
316 if ( p2 > 0.0 ) p = std::sqrt( p );
339 if (
QI[it] < 0 || 849 <
QI[it] ) dqi =
QI[it];
344 }
while(eSecN>MaxEne);
347 eGamm = eKinetic-eSecN;
353 iLevel+=
G4int(random);
360 while (eKinetic-eExcitation < 0 && iLevel>0)
369 if ( dqi < 0 || 849 < dqi ) useQI =
true;
374 eExcitation = -
QI[it];
392 if ( iLevel == -1 ) iLevel = 0;
398 if ( !find ) iLevel = imaxEx;
402 if(getenv(
"G4ParticleHPDebug") && eKinetic-eExcitation < 0)
404 throw G4HadronicException(__FILE__, __LINE__,
"SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
406 if(eKinetic-eExcitation < 0) eExcitation = 0;
425 theRestEnergy->
SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
426 eGamm*std::sin(std::acos(costh))*std::sin(phi),
429 thePhotons->push_back(theRestEnergy);
440 nothingWasKnownOnHadron = 1;
454 boosted_tmp.
Lorentz(incidReactionProduct, theTarget);
459 if(thePhotons!=0 && thePhotons->size()!=0)
460 { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
466 G4bool foundMatchingLevel =
false;
469 for(
G4int j=1; j<it; j++)
479 G4double deltaE = std::abs(testEnergy-aBaseEnergy);
484 thePhotons->push_back(theNext->operator[](0));
485 aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
487 foundMatchingLevel =
true;
496 if(!foundMatchingLevel)
500 thePhotons->push_back(theNext->operator[](0));
501 aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
510 for(i0=0; i0<thePhotons->size(); i0++)
513 thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
517 if(nothingWasKnownOnHadron)
524 if ( thePhotons != 0 )
526 unsigned int nPhotons = thePhotons->size();
528 for ( ii0=0; ii0<nPhotons; ii0++)
531 totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
540 G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
550 G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
599 G4int nSecondaries = 2;
600 G4bool needsSeparateRecoil =
false;
601 G4int totalBaryonNumber = 0;
602 G4int totalCharge = 0;
604 if(theParticles != 0)
606 nSecondaries = theParticles->size();
609 for(ii0=0; ii0<theParticles->size(); ii0++)
611 aDef = theParticles->operator[](ii0)->GetDefinition();
614 totalMomentum += theParticles->operator[](ii0)->GetMomentum();
618 needsSeparateRecoil =
true;
628 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
629 nSecondaries += nPhotons;
633 if( theParticles==0 )
643 aHadron.
Lorentz(aHadron, theTarget);
646 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
654 theResidual.
Lorentz(theResidual, -1.*theTarget);
658 for(i=0; i<nPhotons; i++)
660 totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
673 for(i0=0; i0<theParticles->size(); i0++)
676 theSec->
SetDefinition(theParticles->operator[](i0)->GetDefinition());
677 theSec->
SetMomentum(theParticles->operator[](i0)->GetMomentum());
682 delete theParticles->operator[](i0);
685 if(needsSeparateRecoil && residualZ!=0)
689 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
691 resiualKineticEnergy += totalMomentum*totalMomentum;
692 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.
GetMass();
715 for(i=0; i<nPhotons; i++)
720 theSec->
SetDefinition( thePhotons->operator[](i)->GetDefinition() );
722 theSec->
SetMomentum(thePhotons->operator[](i)->GetMomentum());
728 delete thePhotons->operator[](i);
763 rot2.setTheta( d.theta() );
785 if ( ( 1 + (1+A)/A*Q/E1 ) < 0 )
791 G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
793 G4double E3 = AA/std::pow((1+A),2)*(beta*beta+1+2*beta*mu)*E1;
794 G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu);
795 if ( omega3 > 1.0 ) omega3 = 1.0;
797 G4double E4 = (A+1-AA)/std::pow((1+A),2)*(gamma*gamma+1-2*gamma*mu)*E1;
798 G4double omega4 = (1-gamma*mu)/std::sqrt(gamma*gamma+1-2*gamma*mu);
799 if ( omega4 > 1.0 ) omega4 = 1.0;
804 G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ;
805 G4ThreeVector p ( 0 , pmag*std::sqrt(1-omega3*omega3), pmag*omega3 );
808 G4double pmag4 = std::sqrt ((E4+M4)*(E4+M4)-M4*M4) ;
809 G4ThreeVector p4 ( 0 , -pmag4*std::sqrt(1-omega4*omega4), pmag4*omega4 );
static G4ParticleHPManager * GetInstance()
G4ReactionProductVector * GetDecayGammas(G4int aLevel)
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4ParticleHPEnergyDistribution * theEnergyDistribution[51]
void SetMomentum(const G4ThreeVector &momentum)
G4double GetTotalMomentum() 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]
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 *)
G4ErrorTarget * theTarget
G4ParticleHPLevel * GetLevel(G4int i)
void Init(std::istream &aDataFile)
const G4ParticleDefinition * GetParticleDefinition() const
void InitDistributionInitialState(G4ReactionProduct &anIncidentPart, G4ReactionProduct &aTarget, G4int it)
static const G4double A[nN]
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)
G4HadFinalState theResult
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
G4double GetLevelEnergy(G4int aLevel)
G4int GetNumberOfSecondaries() const
G4ThreeVector GetMomentum() const
void Init(std::istream &aDataFile)
G4int GetBaryonNumber() const
G4double GetLevelEnergy()
CLHEP::HepLorentzVector G4LorentzVector