#include <G4NeutronHPInelasticCompFS.hh>

Inheritance diagram for G4NeutronHPInelasticCompFS:

Public Member Functions
	G4NeutronHPInelasticCompFS ()

virtual	~G4NeutronHPInelasticCompFS ()

void	Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aSFType)

void	InitGammas (G4double AR, G4double ZR)

virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &theTrack)=0

virtual G4NeutronHPFinalState *	New ()=0

virtual G4double	GetXsec (G4double anEnergy)

virtual G4NeutronHPVector *	GetXsec ()

G4int	SelectExitChannel (G4double eKinetic)

void	CompositeApply (const G4HadProjectile &theTrack, G4ParticleDefinition *aHadron)

void	InitDistributionInitialState (G4ReactionProduct &aNeutron, G4ReactionProduct &aTarget, G4int it)

Public Member Functions inherited from G4NeutronHPFinalState
	G4NeutronHPFinalState ()

virtual	~G4NeutronHPFinalState ()

void	Init (G4double A, G4double Z, G4String &dirName, G4String &aFSType)

G4bool	HasXsec ()

G4bool	HasFSData ()

G4bool	HasAnyData ()

void	SetA_Z (G4double anA, G4double aZ, G4int aM=0)

G4double	GetZ ()

G4double	GetN ()

G4int	GetM ()

Protected Attributes
G4NeutronHPVector *	theXsection [51]

G4NeutronHPEnergyDistribution *	theEnergyDistribution [51]

G4NeutronHPAngular *	theAngularDistribution [51]

G4NeutronHPEnAngCorrelation *	theEnergyAngData [51]

G4NeutronHPPhotonDist *	theFinalStatePhotons [51]

G4NeutronHPDeExGammas	theGammas

G4String	gammaPath

G4double	theCurrentA

G4double	theCurrentZ

std::vector< G4double >	QI

std::vector< G4int >	LR

Protected Attributes inherited from G4NeutronHPFinalState
G4bool	hasXsec

G4bool	hasFSData

G4bool	hasAnyData

G4NeutronHPNames	theNames

G4HadFinalState	theResult

G4double	theBaseA

G4double	theBaseZ

G4int	theBaseM

G4int	theNDLDataZ

G4int	theNDLDataA

G4int	theNDLDataM

Additional Inherited Members
Protected Member Functions inherited from G4NeutronHPFinalState
void	SetAZMs (G4double anA, G4double aZ, G4int aM, G4NeutronHPDataUsed used)

void	adjust_final_state (G4LorentzVector)

G4bool	DoNotAdjustFinalState ()

Detailed Description

Definition at line 41 of file G4NeutronHPInelasticCompFS.hh.

Constructor & Destructor Documentation

G4NeutronHPInelasticCompFS::G4NeutronHPInelasticCompFS ( )

inline

Definition at line 45 of file G4NeutronHPInelasticCompFS.hh.

References G4NeutronHPFinalState::hasXsec, LR, QI, theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, and theXsection.

   {
 
     QI.resize(51);
     LR.resize(51);
     for(G4int i=0; i<51; i++)
     {
       hasXsec = true; 
       theXsection[i] = 0;
       theEnergyDistribution[i] = 0;
       theAngularDistribution[i] = 0;
       theEnergyAngData[i] = 0;
       theFinalStatePhotons[i] = 0;
       QI[i]=0.0;
       LR[i]=0;
     }
 
   }

virtual G4NeutronHPInelasticCompFS::~G4NeutronHPInelasticCompFS ( )

inlinevirtual

Definition at line 63 of file G4NeutronHPInelasticCompFS.hh.

References theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, and theXsection.

   {
     for(G4int i=0; i<51; i++)
     {
       if(theXsection[i] != 0) delete theXsection[i];
       if(theEnergyDistribution[i] != 0) delete theEnergyDistribution[i];
       if(theAngularDistribution[i] != 0) delete theAngularDistribution[i];
       if(theEnergyAngData[i] != 0) delete theEnergyAngData[i];
       if(theFinalStatePhotons[i] != 0) delete theFinalStatePhotons[i];
     }
   }

Member Function Documentation

virtual G4HadFinalState* G4NeutronHPInelasticCompFS::ApplyYourself ( const G4HadProjectile & theTrack )

pure virtual

Reimplemented from G4NeutronHPFinalState.

Implemented in G4NeutronHPAInelasticFS, G4NeutronHPDInelasticFS, G4NeutronHPHe3InelasticFS, G4NeutronHPNInelasticFS, G4NeutronHPPInelasticFS, and G4NeutronHPTInelasticFS.

void G4NeutronHPInelasticCompFS::CompositeApply	(	const G4HadProjectile &	theTrack,
		G4ParticleDefinition *	aHadron
	)

Definition at line 219 of file G4NeutronHPInelasticCompFS.cc.

Referenced by G4NeutronHPAInelasticFS::ApplyYourself(), G4NeutronHPTInelasticFS::ApplyYourself(), G4NeutronHPPInelasticFS::ApplyYourself(), G4NeutronHPNInelasticFS::ApplyYourself(), G4NeutronHPHe3InelasticFS::ApplyYourself(), and G4NeutronHPDInelasticFS::ApplyYourself().

 {
 
 // prepare neutron
     theResult.Clear();
     G4double eKinetic = theTrack.GetKineticEnergy();
     const G4HadProjectile *incidentParticle = &theTrack;
     G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
     theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
     theNeutron.SetKineticEnergy( eKinetic );
 
 // prepare target
     G4int i;
     for(i=0; i<50; i++)
     { if(theXsection[i] != 0) { break; } } 
 
     G4double targetMass=0;
     G4double eps = 0.0001;
     targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
                    G4Neutron::Neutron()->GetPDGMass();
 //    if(theEnergyAngData[i]!=0)
 //        targetMass = theEnergyAngData[i]->GetTargetMass();
 //    else if(theAngularDistribution[i]!=0)
 //        targetMass = theAngularDistribution[i]->GetTargetMass();
 //    else if(theFinalStatePhotons[50]!=0)
 //        targetMass = theFinalStatePhotons[50]->GetTargetMass();
     G4Nucleus aNucleus;
     G4ReactionProduct theTarget; 
     G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
     theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
 
 // prepare the residual mass
     G4double residualMass=0;
     G4double residualZ = theBaseZ - aDefinition->GetPDGCharge();
     G4double residualA = theBaseA - aDefinition->GetBaryonNumber()+1;
     residualMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps)) ) /
                      G4Neutron::Neutron()->GetPDGMass();
 
 // prepare energy in target rest frame
     G4ReactionProduct boosted;
     boosted.Lorentz(theNeutron, theTarget);
     eKinetic = boosted.GetKineticEnergy();
 //    G4double momentumInCMS = boosted.GetTotalMomentum();
   
 // select exit channel for composite FS class.
     G4int it = SelectExitChannel( eKinetic );
    
 // set target and neutron in the relevant exit channel
     InitDistributionInitialState(theNeutron, theTarget, it);    
 
     G4ReactionProductVector * thePhotons = 0;
     G4ReactionProductVector * theParticles = 0;
     G4ReactionProduct aHadron;
     aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@    
     G4double availableEnergy = theNeutron.GetKineticEnergy() + theNeutron.GetMass() - aHadron.GetMass() +
                              (targetMass - residualMass)*G4Neutron::Neutron()->GetPDGMass();
 //080730c
     if ( availableEnergy < 0 )
     {
        //G4cout << "080730c Adjust availavleEnergy " << G4endl; 
        availableEnergy = 0; 
     }
     G4int nothingWasKnownOnHadron = 0;
     G4int dummy;
     G4double eGamm = 0;
     G4int iLevel=it-1;
 
 //  TK without photon has it = 0
     if( 50 == it ) 
     {
 
 //    TK Excitation level is not determined
       iLevel=-1;
       aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
                                (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
 
       //aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*
       //                  std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
       //                            aHadron.GetMass()*aHadron.GetMass()));
 
       //TK add safty 100909
       G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy() - aHadron.GetMass()*aHadron.GetMass() );
       G4double p = 0.0;
       if ( p2 > 0.0 ) p = std::sqrt( p ); 
 
       aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*p );
 
     }
     else
     {
       while( iLevel!=-1 && theGammas.GetLevel(iLevel) == 0 ) { iLevel--; }
     }
 
 
     if ( theAngularDistribution[it] != 0 ) // MF4
     {
       if(theEnergyDistribution[it]!=0) // MF5
       {
     //************************************************************
     /*
         aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
         G4double eSecN = aHadron.GetKineticEnergy();
     */
     //************************************************************
     //EMendoza --> maximum allowable energy should be taken into account.
         G4double dqi = 0.0;
         if ( QI[it] < 0 || 849 < QI[it] ) dqi = QI[it]; //For backword compatibility QI introduced since G4NDL3.15
     G4double MaxEne=eKinetic+dqi;
     G4double eSecN;
     do{
       eSecN=theEnergyDistribution[it]->Sample(eKinetic, dummy);
     }while(eSecN>MaxEne);
     aHadron.SetKineticEnergy(eSecN);
     //************************************************************
         eGamm = eKinetic-eSecN;
         for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
         {
           if(theGammas.GetLevelEnergy(iLevel)<eGamm) break;
         }
         G4double random = 2*G4UniformRand();
         iLevel+=G4int(random);
         if(iLevel>theGammas.GetNumberOfLevels()-1)iLevel = theGammas.GetNumberOfLevels()-1;
       }
       else
       {
         G4double eExcitation = 0;
         if(iLevel>=0) eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();    
         while (eKinetic-eExcitation < 0 && iLevel>0)
     {
       iLevel--;
       eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();    
     }
         //110610TK BEGIN
         //Use QI value for calculating excitation energy of residual.
         G4bool useQI=false;
         G4double dqi = QI[it]; 
         if ( dqi < 0 || 849 < dqi ) useQI = true; //Former libraies does not have values of this range
  
         if ( useQI ) 
         {
            // QI introudced since G4NDL3.15
            eExcitation = -QI[it];
            //Re-evluate iLevel based on this eExcitation 
            iLevel = 0;
            G4bool find = false;
            G4int imaxEx = 0;
            while( !theGammas.GetLevel(iLevel+1) == 0 ) 
            { 
               G4double maxEx = 0.0;
               if ( maxEx < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) 
               {
                  maxEx = theGammas.GetLevel(iLevel)->GetLevelEnergy();  
                  imaxEx = iLevel;
               }
               if ( eExcitation < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) 
               {
                  find = true; 
                  iLevel--; 
                  // very small eExcitation, iLevel becomes -1, this is protected below.
                  if ( iLevel == -1 ) iLevel = 0; // But cause energy trouble. 
                  break;
               }
               iLevel++; 
            }
            // In case, cannot find proper level, then use the maximum level. 
            if ( !find ) iLevel = imaxEx;
         }
         //110610TK END
     
     if(getenv("InelasticCompFSLogging") && eKinetic-eExcitation < 0) 
     {
       throw G4HadronicException(__FILE__, __LINE__, "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
     }
     if(eKinetic-eExcitation < 0) eExcitation = 0;
     if(iLevel!= -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
     
       }
       theAngularDistribution[it]->SampleAndUpdate(aHadron);
 
       if( theFinalStatePhotons[it] == 0 )
       {
         //G4cout << "110610 USE Gamma Level" << G4endl;
 // TK comment Most n,n* eneter to this  
     thePhotons = theGammas.GetDecayGammas(iLevel);
     eGamm -= theGammas.GetLevelEnergy(iLevel);
     if(eGamm>0) // @ ok for now, but really needs an efficient way of correllated sampling @
     {
           G4ReactionProduct * theRestEnergy = new G4ReactionProduct;
           theRestEnergy->SetDefinition(G4Gamma::Gamma());
           theRestEnergy->SetKineticEnergy(eGamm);
           G4double costh = 2.*G4UniformRand()-1.;
           G4double phi = twopi*G4UniformRand();
           theRestEnergy->SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi), 
                                      eGamm*std::sin(std::acos(costh))*std::sin(phi),
                                      eGamm*costh);
           if(thePhotons == 0) { thePhotons = new G4ReactionProductVector; }
           thePhotons->push_back(theRestEnergy);
     }
       }
     }
     else if(theEnergyAngData[it] != 0) // MF6  
     {
       theParticles = theEnergyAngData[it]->Sample(eKinetic);
     }
     else
     {
       // @@@ what to do, if we have photon data, but no info on the hadron itself
       nothingWasKnownOnHadron = 1;
     }
 
     //G4cout << "theFinalStatePhotons it " << it << G4endl;
     //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
     //G4cout << "theFinalStatePhotons it " << it << G4endl;
     //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
     //G4cout << "thePhotons " << thePhotons << G4endl;
 
     if ( theFinalStatePhotons[it] != 0 ) 
     {
        // the photon distributions are in the Nucleus rest frame.
        // TK residual rest frame
       G4ReactionProduct boosted_tmp;
       boosted_tmp.Lorentz(theNeutron, theTarget);
       G4double anEnergy = boosted_tmp.GetKineticEnergy();
       thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
       G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
       G4double testEnergy = 0;
       if(thePhotons!=0 && thePhotons->size()!=0)
       { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
       if(theFinalStatePhotons[it]->NeedsCascade())
       {
     while(aBaseEnergy>0.01*keV)
         {
           // cascade down the levels
       G4bool foundMatchingLevel = false;
           G4int closest = 2;
       G4double deltaEold = -1;
       for(G4int j=1; j<it; j++)
           {
             if(theFinalStatePhotons[j]!=0) 
             {
               testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
             }
             else
             {
               testEnergy = 0;
             }
         G4double deltaE = std::abs(testEnergy-aBaseEnergy);
             if(deltaE<0.1*keV)
             {
               G4ReactionProductVector * theNext = 
             theFinalStatePhotons[j]->GetPhotons(anEnergy);
               thePhotons->push_back(theNext->operator[](0));
               aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
               delete theNext;
           foundMatchingLevel = true;
               break; // ===>
             }
         if(theFinalStatePhotons[j]!=0 && ( deltaE<deltaEold||deltaEold<0.) )
         {
           closest = j;
           deltaEold = deltaE;     
         }
           } // <=== the break goes here.
       if(!foundMatchingLevel)
       {
             G4ReactionProductVector * theNext = 
                theFinalStatePhotons[closest]->GetPhotons(anEnergy);
             thePhotons->push_back(theNext->operator[](0));
             aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
             delete theNext;
       }
         } 
       }
     }
     unsigned int i0;
     if(thePhotons!=0)
     {
       for(i0=0; i0<thePhotons->size(); i0++)
       {
     // back to lab
     thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
       }
     }
     //G4cout << "nothingWasKnownOnHadron " << nothingWasKnownOnHadron << G4endl;
     if(nothingWasKnownOnHadron)
     {
 //    TKDB 100405
 //    In this case, hadron should be isotropic in CM
 //    mu and p should be correlated
 //
       G4double totalPhotonEnergy = 0.0;
       if ( thePhotons != 0 )
       {
          unsigned int nPhotons = thePhotons->size();
          unsigned int ii0;
          for ( ii0=0; ii0<nPhotons; ii0++)
          {
             //thePhotons has energies at LAB system 
             totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
          }
       }
 
       //isotropic distribution in CM 
       G4double mu = 1.0 - 2 * G4UniformRand();
 
       // need momentums in target rest frame;
       G4LorentzVector target_in_LAB ( theTarget.GetMomentum() , theTarget.GetTotalEnergy() );
       G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
       G4LorentzVector proj_in_LAB = incidentParticle->Get4Momentum();
 
       G4DynamicParticle* proj = new G4DynamicParticle( G4Neutron::Neutron() , proj_in_LAB.boost( boostToTargetRest ) ); 
       G4DynamicParticle* targ = new G4DynamicParticle( G4IonTable::GetIonTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , totalPhotonEnergy )  , G4ThreeVector(0) );
       G4DynamicParticle* hadron = new G4DynamicParticle( aHadron.GetDefinition() , G4ThreeVector(0) );  // will be fill momentum
 
       two_body_reaction ( proj , targ , hadron , mu );
 
       G4LorentzVector hadron_in_trag_rest = hadron->Get4Momentum();
       G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
       aHadron.SetMomentum( hadron_in_LAB.v() );
       aHadron.SetKineticEnergy ( hadron_in_LAB.e() - hadron_in_LAB.m() );
 
       delete proj;
       delete targ; 
       delete hadron;
 
 //TKDB 100405
 /*
       G4double totalPhotonEnergy = 0;
       if(thePhotons!=0)
       {
         unsigned int nPhotons = thePhotons->size();
     unsigned int i0;
     for(i0=0; i0<nPhotons; i0++)
         {
           totalPhotonEnergy += thePhotons->operator[](i0)->GetTotalEnergy();
         }
       }
       availableEnergy -= totalPhotonEnergy;
       residualMass += totalPhotonEnergy/G4Neutron::Neutron()->GetPDGMass();
       aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
                                (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
       G4double CosTheta = 1.0 - 2.0*G4UniformRand();
       G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
       G4double Phi = twopi*G4UniformRand();
       G4ThreeVector Vector(std::cos(Phi)*SinTheta, std::sin(Phi)*SinTheta, CosTheta);
       //aHadron.SetMomentum(Vector* std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
       //                                 aHadron.GetMass()*aHadron.GetMass()));
       G4double p2 = aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()- aHadron.GetMass()*aHadron.GetMass();
 
       G4double p = 0.0;
       if ( p2 > 0.0 )
          p = std::sqrt ( p2 ); 
 
       aHadron.SetMomentum( Vector*p ); 
 */
 
     }
 
 // fill the result
 // Beware - the recoil is not necessarily in the particles...
 // Can be calculated from momentum conservation?
 // The idea is that the particles ar emitted forst, and the gammas only once the
 // recoil is on the residual; assumption is that gammas do not contribute to 
 // the recoil.
 // This needs more design @@@
 
     G4int nSecondaries = 2; // the hadron and the recoil
     G4bool needsSeparateRecoil = false;
     G4int totalBaryonNumber = 0;
     G4int totalCharge = 0;
     G4ThreeVector totalMomentum(0);
     if(theParticles != 0) 
     {
       nSecondaries = theParticles->size();
       G4ParticleDefinition * aDef;
       unsigned int ii0;
       for(ii0=0; ii0<theParticles->size(); ii0++)
       {
         aDef = theParticles->operator[](ii0)->GetDefinition();
     totalBaryonNumber+=aDef->GetBaryonNumber();
     totalCharge+=G4int(aDef->GetPDGCharge()+eps);
         totalMomentum += theParticles->operator[](ii0)->GetMomentum();
       } 
       if(totalBaryonNumber!=G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber())) 
       {
         needsSeparateRecoil = true;
     nSecondaries++;
     residualA = G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber()
                       -totalBaryonNumber);
     residualZ = G4int(theBaseZ+eps+incidentParticle->GetDefinition()->GetPDGCharge()
                       -totalCharge);
       }
     }
     
     G4int nPhotons = 0;
     if(thePhotons!=0) { nPhotons = thePhotons->size(); }
     nSecondaries += nPhotons;
         
     G4DynamicParticle * theSec;
     
     if( theParticles==0 )
     {
       theSec = new G4DynamicParticle;   
       theSec->SetDefinition(aHadron.GetDefinition());
       theSec->SetMomentum(aHadron.GetMomentum());
       theResult.AddSecondary(theSec);    
  
     aHadron.Lorentz(aHadron, theTarget);
         G4ReactionProduct theResidual;   
         theResidual.SetDefinition(G4IonTable::GetIonTable()
                               ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));  
         theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass());
 
         //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
         //theResidual.SetMomentum(-1.*aHadron.GetMomentum());
     G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
         theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
 
     theResidual.Lorentz(theResidual, -1.*theTarget);
     G4ThreeVector totalPhotonMomentum(0,0,0);
     if(thePhotons!=0)
     {
           for(i=0; i<nPhotons; i++)
           {
             totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
           }
     }
         theSec = new G4DynamicParticle;   
         theSec->SetDefinition(theResidual.GetDefinition());
         theSec->SetMomentum(theResidual.GetMomentum()-totalPhotonMomentum);
         theResult.AddSecondary(theSec);    
     }
     else
     {
       for(i0=0; i0<theParticles->size(); i0++)
       {
         theSec = new G4DynamicParticle; 
         theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition());
         theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum());
         theResult.AddSecondary(theSec); 
         delete theParticles->operator[](i0); 
       } 
       delete theParticles;
       if(needsSeparateRecoil && residualZ!=0)
       {
         G4ReactionProduct theResidual;   
         theResidual.SetDefinition(G4IonTable::GetIonTable()
                               ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));  
         G4double resiualKineticEnergy  = theResidual.GetMass()*theResidual.GetMass();
                  resiualKineticEnergy += totalMomentum*totalMomentum;
              resiualKineticEnergy  = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
 //        cout << "Kinetic energy of the residual = "<<resiualKineticEnergy<<endl;
     theResidual.SetKineticEnergy(resiualKineticEnergy);
 
         //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
         //theResidual.SetMomentum(-1.*totalMomentum);
     //G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
         //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
 //080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
         theResidual.SetMomentum( theNeutron.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
 
         theSec = new G4DynamicParticle;   
         theSec->SetDefinition(theResidual.GetDefinition());
         theSec->SetMomentum(theResidual.GetMomentum());
         theResult.AddSecondary(theSec);  
       }  
     }
     if(thePhotons!=0)
     {
       for(i=0; i<nPhotons; i++)
       {
         theSec = new G4DynamicParticle;    
         //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 
         //theSec->SetDefinition(G4Gamma::Gamma());
         theSec->SetDefinition( thePhotons->operator[](i)->GetDefinition() );
         //But never cause real effect at least with G4NDL3.13 TK
         theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
         theResult.AddSecondary(theSec); 
         delete thePhotons->operator[](i);
       }
 // some garbage collection
       delete thePhotons;
     }
 
 //080721 
    G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
    G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
    G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
    G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
    adjust_final_state ( init_4p_lab ); 
 
 // clean up the primary neutron
     theResult.SetStatusChange(stopAndKill);
 }

virtual G4double G4NeutronHPInelasticCompFS::GetXsec ( G4double anEnergy )

inlinevirtual

Reimplemented from G4NeutronHPFinalState.

Definition at line 78 of file G4NeutronHPInelasticCompFS.hh.

References G4INCL::Math::max(), and theXsection.

   {
     return std::max(0., theXsection[50]->GetY(anEnergy));
   }

virtual G4NeutronHPVector* G4NeutronHPInelasticCompFS::GetXsec ( )

inlinevirtual

Reimplemented from G4NeutronHPFinalState.

Definition at line 82 of file G4NeutronHPInelasticCompFS.hh.

References theXsection.

Referenced by SelectExitChannel().

82 { return theXsection[50]; }

G4NeutronHPInelasticCompFS::theXsection

G4NeutronHPVector * theXsection[51]

Definition: G4NeutronHPInelasticCompFS.hh:103

void G4NeutronHPInelasticCompFS::Init	(	G4double	A,
		G4double	Z,
		G4int	M,
		G4String &	dirName,
		G4String &	aSFType
	)

virtual

Implements G4NeutronHPFinalState.

Reimplemented in G4NeutronHPNInelasticFS, G4NeutronHPPInelasticFS, and G4NeutronHPTInelasticFS.

Definition at line 79 of file G4NeutronHPInelasticCompFS.cc.

References python.hepunit::eV, G4cout, G4endl, gammaPath, G4NeutronHPManager::GetDataStream(), G4NeutronHPManager::GetInstance(), G4NeutronHPNames::GetName(), G4NeutronHPDataUsed::GetName(), G4NeutronHPFinalState::hasAnyData, G4NeutronHPFinalState::hasFSData, G4NeutronHPFinalState::hasXsec, G4NeutronHPEnAngCorrelation::Init(), G4NeutronHPAngular::Init(), G4NeutronHPEnergyDistribution::Init(), G4NeutronHPVector::Init(), G4NeutronHPPhotonDist::InitAngular(), G4NeutronHPPhotonDist::InitEnergies(), G4NeutronHPPhotonDist::InitMean(), G4NeutronHPPhotonDist::InitPartials(), LR, QI, G4NeutronHPFinalState::SetAZMs(), theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, G4NeutronHPFinalState::theNames, G4NeutronHPFinalState::theNDLDataA, G4NeutronHPFinalState::theNDLDataZ, theXsection, and G4INCL::CrossSections::total().

Referenced by G4NeutronHPAInelasticFS::Init(), G4NeutronHPTInelasticFS::Init(), G4NeutronHPPInelasticFS::Init(), G4NeutronHPNInelasticFS::Init(), G4NeutronHPHe3InelasticFS::Init(), and G4NeutronHPDInelasticFS::Init().

 {
   gammaPath = "/Inelastic/Gammas/";
     if(!getenv("G4NEUTRONHPDATA")) 
        throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
   G4String tBase = getenv("G4NEUTRONHPDATA");
   gammaPath = tBase+gammaPath;
   G4String tString = dirName;
   G4bool dbool;
   G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aFSType, dbool);
   G4String filename = aFile.GetName();
   SetAZMs( A, Z, M, aFile ); 
   //theBaseA = aFile.GetA();
   //theBaseZ = aFile.GetZ();
   //theNDLDataA = (int)aFile.GetA();
   //theNDLDataZ = aFile.GetZ();
   //if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
   if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) )
   {
     if(getenv("NeutronHPNamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
     hasAnyData = false;
     hasFSData = false; 
     hasXsec = false;
     return;
   }
   //theBaseA = A;
   //theBaseZ = G4int(Z+.5);
    //std::ifstream theData(filename, std::ios::in);
    std::istringstream theData(std::ios::in);
    G4NeutronHPManager::GetInstance()->GetDataStream(filename,theData);
   if(!theData) //"!" is a operator of ios
   {
     hasAnyData = false;
     hasFSData = false; 
     hasXsec = false;
     //theData.close();
     return;
   }
   // here we go
   G4int infoType, dataType, dummy;
   G4int sfType, it;
   hasFSData = false; 
   while (theData >> infoType)
   {
     hasFSData = true; 
     theData >> dataType;
     theData >> sfType >> dummy;
     it = 50;
     if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
     if(dataType==3) 
     {
       //theData >> dummy >> dummy;
       //TK110430
       // QI and LR introudced since G4NDL3.15
       G4double dqi;
       G4int ilr;
       theData >> dqi >> ilr;
 
       QI[ it ] = dqi*eV;
       LR[ it ] = ilr;
       theXsection[it] = new G4NeutronHPVector;
       G4int total;
       theData >> total;
       theXsection[it]->Init(theData, total, eV);
       //std::cout << theXsection[it]->GetXsec(1*MeV) << std::endl;
     }
     else if(dataType==4)
     {
       theAngularDistribution[it] = new G4NeutronHPAngular;
       theAngularDistribution[it]->Init(theData);
     }
     else if(dataType==5)
     {
       theEnergyDistribution[it] = new G4NeutronHPEnergyDistribution;
       theEnergyDistribution[it]->Init(theData); 
     }
     else if(dataType==6)
     {
       theEnergyAngData[it] = new G4NeutronHPEnAngCorrelation;
       theEnergyAngData[it]->Init(theData);
     }
     else if(dataType==12)
     {
       theFinalStatePhotons[it] = new G4NeutronHPPhotonDist;
       theFinalStatePhotons[it]->InitMean(theData);
     }
     else if(dataType==13)
     {
       theFinalStatePhotons[it] = new G4NeutronHPPhotonDist;
       theFinalStatePhotons[it]->InitPartials(theData);
     }
     else if(dataType==14)
     {
       theFinalStatePhotons[it]->InitAngular(theData);
     }
     else if(dataType==15)
     {
       theFinalStatePhotons[it]->InitEnergies(theData);
     }
     else
     {
       throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4NeutronHPInelasticCompFS");
     }
   }
   //theData.close();
 }

void G4NeutronHPInelasticCompFS::InitDistributionInitialState	(	G4ReactionProduct &	aNeutron,
		G4ReactionProduct &	aTarget,
		G4int	it
	)

inline

Definition at line 85 of file G4NeutronHPInelasticCompFS.hh.

References G4NeutronHPAngular::SetNeutron(), G4NeutronHPEnAngCorrelation::SetNeutron(), G4NeutronHPAngular::SetTarget(), G4NeutronHPEnAngCorrelation::SetTarget(), theAngularDistribution, and theEnergyAngData.

Referenced by CompositeApply().

   {
     if(theAngularDistribution[it]!=0) 
     {
       theAngularDistribution[it]->SetTarget(aTarget);
       theAngularDistribution[it]->SetNeutron(aNeutron);
     }
     if(theEnergyAngData[it]!=0)
     {
       theEnergyAngData[it]->SetTarget(aTarget);
       theEnergyAngData[it]->SetNeutron(aNeutron);
     }
   }

void G4NeutronHPInelasticCompFS::InitGammas	(	G4double	AR,
		G4double	ZR
	)

Definition at line 58 of file G4NeutronHPInelasticCompFS.cc.

References gammaPath, G4NeutronHPDeExGammas::Init(), and theGammas.

Referenced by G4NeutronHPAInelasticFS::Init(), G4NeutronHPTInelasticFS::Init(), G4NeutronHPPInelasticFS::Init(), G4NeutronHPNInelasticFS::Init(), G4NeutronHPHe3InelasticFS::Init(), and G4NeutronHPDInelasticFS::Init().

 {
   //   char the[100] = {""};
   //   std::ostrstream ost(the, 100, std::ios::out);
   //   ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
   //   G4String * aName = new G4String(the);
   //   std::ifstream from(*aName, std::ios::in);
 
    std::ostringstream ost;
    ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
    G4String aName = ost.str();
    std::ifstream from(aName, std::ios::in);
 
    if(!from) return; // no data found for this isotope
    //   std::ifstream theGammaData(*aName, std::ios::in);
    std::ifstream theGammaData(aName, std::ios::in);
     
    theGammas.Init(theGammaData);
    //   delete aName;
 }

virtual G4NeutronHPFinalState* G4NeutronHPInelasticCompFS::New ( )

pure virtual

Implements G4NeutronHPFinalState.

Implemented in G4NeutronHPAInelasticFS, G4NeutronHPDInelasticFS, G4NeutronHPHe3InelasticFS, G4NeutronHPNInelasticFS, G4NeutronHPPInelasticFS, and G4NeutronHPTInelasticFS.

G4int G4NeutronHPInelasticCompFS::SelectExitChannel ( G4double eKinetic )

Definition at line 186 of file G4NeutronHPInelasticCompFS.cc.

References G4UniformRand, GetXsec(), G4INCL::Math::max(), and theXsection.

Referenced by CompositeApply().

 {
 
 // it = 0 has without Photon
   G4double running[50];
   running[0] = 0;
   unsigned int i;
   for(i=0; i<50; i++)
   {
     if(i!=0) running[i]=running[i-1];
     if(theXsection[i] != 0) 
     {
       running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
     }
   }
   G4double random = G4UniformRand();
   G4double sum = running[49];
   G4int it = 50;
   if(0!=sum)
   {
     G4int i0;
     for(i0=0; i0<50; i0++)
     {
       it = i0;
       if(random < running[i0]/sum) break;
     }
   }
 //debug:  it = 1;
   return it;
 }