65 std::ostringstream ost;
81 if(!getenv(
"G4NEUTRONHPDATA"))
82 throw G4HadronicException(__FILE__, __LINE__,
"Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
83 G4String tBase = getenv(
"G4NEUTRONHPDATA");
97 if(getenv(
"NeutronHPNamesLogging"))
G4cout <<
"Skipped = "<< filename <<
" "<<A<<
" "<<Z<<
G4endl;
115 G4int infoType, dataType, dummy;
118 while (theData >> infoType)
122 theData >> sfType >> dummy;
124 if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
132 theData >> dqi >> ilr;
157 else if(dataType==12)
162 else if(dataType==13)
167 else if(dataType==14)
171 else if(dataType==15)
177 throw G4HadronicException(__FILE__, __LINE__,
"Data-type unknown to G4NeutronHPInelasticCompFS");
192 if(i!=0) running[i]=running[i-1];
204 for(i0=0; i0<50; i0++)
207 if(random < running[i0]/sum)
break;
256 boosted.
Lorentz(theNeutron, theTarget);
273 if ( availableEnergy < 0 )
278 G4int nothingWasKnownOnHadron = 0;
299 if ( p2 > 0.0 ) p = std::sqrt( p );
322 if (
QI[it] < 0 || 849 <
QI[it] ) dqi =
QI[it];
327 }
while(eSecN>MaxEne);
330 eGamm = eKinetic-eSecN;
336 iLevel+=
G4int(random);
343 while (eKinetic-eExcitation < 0 && iLevel>0)
352 if ( dqi < 0 || 849 < dqi ) useQI =
true;
357 eExcitation = -
QI[it];
375 if ( iLevel == -1 ) iLevel = 0;
381 if ( !find ) iLevel = imaxEx;
385 if(getenv(
"InelasticCompFSLogging") && eKinetic-eExcitation < 0)
387 throw G4HadronicException(__FILE__, __LINE__,
"SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
389 if(eKinetic-eExcitation < 0) eExcitation = 0;
408 theRestEnergy->
SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
409 eGamm*std::sin(std::acos(costh))*std::sin(phi),
412 thePhotons->push_back(theRestEnergy);
423 nothingWasKnownOnHadron = 1;
437 boosted_tmp.
Lorentz(theNeutron, theTarget);
442 if(thePhotons!=0 && thePhotons->size()!=0)
443 { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
446 while(aBaseEnergy>0.01*
keV)
449 G4bool foundMatchingLevel =
false;
452 for(
G4int j=1; j<it; j++)
462 G4double deltaE = std::abs(testEnergy-aBaseEnergy);
467 thePhotons->push_back(theNext->operator[](0));
468 aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
470 foundMatchingLevel =
true;
479 if(!foundMatchingLevel)
483 thePhotons->push_back(theNext->operator[](0));
484 aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
493 for(i0=0; i0<thePhotons->size(); i0++)
496 thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
500 if(nothingWasKnownOnHadron)
507 if ( thePhotons != 0 )
509 unsigned int nPhotons = thePhotons->size();
511 for ( ii0=0; ii0<nPhotons; ii0++)
514 totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
523 G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
530 two_body_reaction ( proj , targ , hadron , mu );
582 G4int nSecondaries = 2;
583 G4bool needsSeparateRecoil =
false;
584 G4int totalBaryonNumber = 0;
585 G4int totalCharge = 0;
587 if(theParticles != 0)
589 nSecondaries = theParticles->size();
592 for(ii0=0; ii0<theParticles->size(); ii0++)
594 aDef = theParticles->operator[](ii0)->GetDefinition();
597 totalMomentum += theParticles->operator[](ii0)->GetMomentum();
601 needsSeparateRecoil =
true;
611 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
612 nSecondaries += nPhotons;
616 if( theParticles==0 )
623 aHadron.
Lorentz(aHadron, theTarget);
626 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
634 theResidual.
Lorentz(theResidual, -1.*theTarget);
638 for(i=0; i<nPhotons; i++)
640 totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
650 for(i0=0; i0<theParticles->size(); i0++)
653 theSec->
SetDefinition(theParticles->operator[](i0)->GetDefinition());
654 theSec->
SetMomentum(theParticles->operator[](i0)->GetMomentum());
656 delete theParticles->operator[](i0);
659 if(needsSeparateRecoil && residualZ!=0)
663 ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
665 resiualKineticEnergy += totalMomentum*totalMomentum;
666 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.
GetMass();
685 for(i=0; i<nPhotons; i++)
690 theSec->
SetDefinition( thePhotons->operator[](i)->GetDefinition() );
692 theSec->
SetMomentum(thePhotons->operator[](i)->GetMomentum());
694 delete thePhotons->operator[](i);
751 if ( ( 1 + (1+A)/A*Q/E1 ) < 0 )
754 Q = -( A/(1+A)*E1 ) + 1.0e-6*
eV;
757 G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
759 G4double E3 = AA/std::pow((1+A),2)*(beta*beta+1+2*beta*mu)*E1;
760 G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu);
761 if ( omega3 > 1.0 ) omega3 = 1.0;
763 G4double E4 = (A+1-AA)/std::pow((1+A),2)*(gamma*gamma+1-2*gamma*mu)*E1;
764 G4double omega4 = (1-gamma*mu)/std::sqrt(gamma*gamma+1-2*gamma*mu);
765 if ( omega4 > 1.0 ) omega4 = 1.0;
770 G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ;
771 G4ThreeVector p ( 0 , pmag*std::sqrt(1-omega3*omega3), pmag*omega3 );
774 G4double pmag4 = std::sqrt ((E4+M4)*(E4+M4)-M4*M4) ;
775 G4ThreeVector p4 ( 0 , -pmag4*std::sqrt(1-omega4*omega4), pmag4*omega4 );