#include <G4ParticleHPInelasticCompFS.hh>

Inheritance diagram for G4ParticleHPInelasticCompFS:

Collaboration diagram for G4ParticleHPInelasticCompFS:

Public Member Functions
	G4ParticleHPInelasticCompFS ()

virtual	~G4ParticleHPInelasticCompFS ()

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

void	InitGammas (G4double AR, G4double ZR)

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

virtual G4ParticleHPFinalState *	New ()=0

virtual G4double	GetXsec (G4double anEnergy)

virtual G4ParticleHPVector *	GetXsec ()

G4int	SelectExitChannel (G4double eKinetic)

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

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

Public Member Functions inherited from G4ParticleHPFinalState
	G4ParticleHPFinalState ()

virtual	~G4ParticleHPFinalState ()

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

G4bool	HasXsec ()

G4bool	HasFSData ()

G4bool	HasAnyData ()

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

G4double	GetZ ()

G4double	GetN ()

G4double	GetA ()

G4int	GetM ()

void	SetAZMs (G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)

void	SetProjectile (G4ParticleDefinition *projectile)

Protected Attributes
G4ParticleHPVector *	theXsection [51]

G4ParticleHPEnergyDistribution *	theEnergyDistribution [51]

G4ParticleHPAngular *	theAngularDistribution [51]

G4ParticleHPEnAngCorrelation *	theEnergyAngData [51]

G4ParticleHPPhotonDist *	theFinalStatePhotons [51]

G4ParticleHPDeExGammas	theGammas

G4String	gammaPath

std::vector< G4double >	QI

std::vector< G4int >	LR

Protected Attributes inherited from G4ParticleHPFinalState
G4bool	hasXsec

G4bool	hasFSData

G4bool	hasAnyData

G4ParticleHPNames	theNames

G4Cache< G4HadFinalState * >	theResult

G4double	theBaseA

G4double	theBaseZ

G4int	theBaseM

G4int	theNDLDataZ

G4int	theNDLDataA

G4int	theNDLDataM

G4ParticleDefinition *	theProjectile

Additional Inherited Members
Protected Member Functions inherited from G4ParticleHPFinalState
void	adjust_final_state (G4LorentzVector)

G4bool	DoNotAdjustFinalState ()

Detailed Description

Definition at line 43 of file G4ParticleHPInelasticCompFS.hh.

Constructor & Destructor Documentation

G4ParticleHPInelasticCompFS::G4ParticleHPInelasticCompFS ( )

inline

Definition at line 47 of file G4ParticleHPInelasticCompFS.hh.

   {
 
     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 G4ParticleHPInelasticCompFS::~G4ParticleHPInelasticCompFS ( )

inlinevirtual

Definition at line 65 of file G4ParticleHPInelasticCompFS.hh.

   {
     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* G4ParticleHPInelasticCompFS::ApplyYourself ( const G4HadProjectile & theTrack )

pure virtual

Reimplemented from G4ParticleHPFinalState.

Implemented in G4ParticleHPAInelasticFS, G4ParticleHPDInelasticFS, G4ParticleHPHe3InelasticFS, G4ParticleHPNInelasticFS, G4ParticleHPPInelasticFS, and G4ParticleHPTInelasticFS.

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

Definition at line 230 of file G4ParticleHPInelasticCompFS.cc.

 {
 
 // prepare neutron
     if ( theResult.Get() == NULL ) theResult.Put( new G4HadFinalState );
     theResult.Get()->Clear();
     G4double eKinetic = theTrack.GetKineticEnergy();
     const G4HadProjectile *hadProjectile = &theTrack;
     G4ReactionProduct incidReactionProduct( const_cast<G4ParticleDefinition *>(hadProjectile->GetDefinition()) ); // incidReactionProduct
     incidReactionProduct.SetMomentum( hadProjectile->Get4Momentum().vect() );
     incidReactionProduct.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));
 #ifdef G4PHPDEBUG
     if( getenv("G4ParticleHPDebug"))  G4cout <<this <<" G4ParticleHPInelasticCompFS::CompositeApply A " <<theBaseA <<" Z " <<theBaseZ <<" incident " <<hadProjectile->GetDefinition()->GetParticleName() <<G4endl;
 #endif
 //    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();
     G4ReactionProduct theTarget; 
     G4Nucleus aNucleus;
     //G4ThreeVector neuVelo = (1./hadProjectile->GetDefinition()->GetPDGMass())*incidReactionProduct.GetMomentum();
     //theTarget = aNucleus.GetBiasedThermalNucleus( targetMass/hadProjectile->GetDefinition()->GetPDGMass() , neuVelo, theTrack.GetMaterial()->GetTemperature());
     //G4Nucleus::GetBiasedThermalNucleus requests normalization of mass and velocity in neutron mass
     G4ThreeVector neuVelo = ( 1./G4Neutron::Neutron()->GetPDGMass() )*incidReactionProduct.GetMomentum();
     theTarget = aNucleus.GetBiasedThermalNucleus( targetMass/G4Neutron::Neutron()->GetPDGMass()
                                                 , neuVelo, theTrack.GetMaterial()->GetTemperature() );
     
     theTarget.SetDefinition( G4IonTable::GetIonTable()->GetIon( G4int(theBaseZ), G4int(theBaseA) , 0.0 ) );  //XX
 
 // prepare the residual mass
     G4double residualMass=0;
     G4double residualZ = theBaseZ + theProjectile->GetPDGCharge() - aDefinition->GetPDGCharge();
     G4double residualA = theBaseA + theProjectile->GetBaryonNumber() - aDefinition->GetBaryonNumber();
     residualMass = G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps));
 
 // prepare energy in target rest frame
     G4ReactionProduct boosted;
     boosted.Lorentz(incidReactionProduct, 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(incidReactionProduct, theTarget, it);    
 
     G4ReactionProductVector * thePhotons = 0;
     G4ReactionProductVector * theParticles = 0;
     G4ReactionProduct aHadron;
     aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@    
     G4double availableEnergy = incidReactionProduct.GetKineticEnergy() + incidReactionProduct.GetMass() - aHadron.GetMass() +
                              (targetMass - residualMass);
 //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/
                                (aHadron.GetMass()+residualMass));
 
       //aHadron.SetMomentum(incidReactionProduct.GetMomentum()*(1./incidReactionProduct.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(incidReactionProduct.GetMomentum()*(1./incidReactionProduct.GetTotalMomentum())*p );
 
     }
     else
     {
       while ( iLevel!=-1 && theGammas.GetLevel(iLevel) == 0 ) { iLevel--; } // Loop checking, 11.05.2015, T. Koi
     }
 
 
     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;
 
         G4int icounter=0;
         G4int icounter_max=1024;
         do {
            icounter++;
            if ( icounter > icounter_max ) {
           G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
               break;
            }
        eSecN=theEnergyDistribution[it]->Sample(eKinetic, dummy);
     }while(eSecN>MaxEne); // Loop checking, 11.05.2015, T. Koi
     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) // Loop checking, 11.05.2015, T. Koi
     {
       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 ) // Loop checking, 11.05.2015, T. Koi
            { 
               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("G4ParticleHPDebug") && 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 = CLHEP::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);
 
       //141017 Fix BEGIN
       //Adjust A and Z in the case of miss much between selected data and target nucleus 
       if ( theParticles != NULL ) {
          G4int sumA = 0;
          G4int sumZ = 0;
          G4int maxA = 0;
          G4int jAtMaxA = 0;
          for ( G4int j = 0 ; j != (G4int)theParticles->size() ; j++ ) {
             if ( theParticles->at(j)->GetDefinition()->GetBaryonNumber() > maxA ) {
                maxA = theParticles->at(j)->GetDefinition()->GetBaryonNumber(); 
                jAtMaxA = j; 
             }
             sumA += theParticles->at(j)->GetDefinition()->GetBaryonNumber();
             sumZ += G4int( theParticles->at(j)->GetDefinition()->GetPDGCharge() + eps );
          }
          G4int dA = (G4int)theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
          G4int dZ = (G4int)theBaseZ + G4int( hadProjectile->GetDefinition()->GetPDGCharge() + eps ) - sumZ;
          if ( dA < 0 || dZ < 0 ) {
             G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA ;
             G4int newZ = G4int( theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge() + eps ) + dZ;
             G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon ( newZ , newA );
             theParticles->at( jAtMaxA )->SetDefinition( pd );
          }
       }
       //141017 Fix END
 
     }
     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(incidReactionProduct, 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*CLHEP::keV) // Loop checking, 11.05.2015, T. Koi
         {
           // 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*CLHEP::keV)
             {
               G4ReactionProductVector * theNext = 
             theFinalStatePhotons[j]->GetPhotons(anEnergy);
               if ( thePhotons != NULL ) 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);
             if ( thePhotons != NULL ) 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 = hadProjectile->Get4Momentum();
 
       G4DynamicParticle* proj = new G4DynamicParticle( theProjectile , 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/theProjectile->GetPDGMass();
       aHadron.SetKineticEnergy(availableEnergy*residualMass*theProjectile->GetPDGMass()/
                                (aHadron.GetMass()+residualMass*theProjectile->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();
       const 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+hadProjectile->GetDefinition()->GetBaryonNumber())) 
       {
         needsSeparateRecoil = true;
     nSecondaries++;
     residualA = G4int(theBaseA+eps+hadProjectile->GetDefinition()->GetBaryonNumber()
                       -totalBaryonNumber);
     residualZ = G4int(theBaseZ+eps+hadProjectile->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.Get()->AddSecondary(theSec);    
 #ifdef G4PHPDEBUG
       if( getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply  add secondary1 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
 #endif
       
       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 = incidReactionProduct.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.Get()->AddSecondary(theSec);    
 #ifdef G4PHPDEBUG
       if( getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary2 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
 #endif
     }
     else
     {
       for(i0=0; i0<theParticles->size(); i0++)
       {
         theSec = new G4DynamicParticle; 
         theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition());
         theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum());
         theResult.Get()->AddSecondary(theSec); 
 #ifdef G4PHPDEBUG
       if( getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary3 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
 #endif
         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 = incidReactionProduct.GetMomentum();
         //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
 //080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
         theResidual.SetMomentum( incidReactionProduct.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
 
         theSec = new G4DynamicParticle;   
         theSec->SetDefinition(theResidual.GetDefinition());
         theSec->SetMomentum(theResidual.GetMomentum());
         theResult.Get()->AddSecondary(theSec);  
 #ifdef G4PHPDEBUG
       if( getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary4 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
 #endif
 
       }  
     }
     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.Get()->AddSecondary(theSec); 
 #ifdef G4PHPDEBUG
       if( getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary5 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
 #endif
 
         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.Get()->SetStatusChange( stopAndKill );
 }

Here is the caller graph for this function:

virtual G4double G4ParticleHPInelasticCompFS::GetXsec ( G4double anEnergy )

inlinevirtual

Reimplemented from G4ParticleHPFinalState.

Definition at line 80 of file G4ParticleHPInelasticCompFS.hh.

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

Here is the call graph for this function:

virtual G4ParticleHPVector* G4ParticleHPInelasticCompFS::GetXsec ( )

inlinevirtual

Reimplemented from G4ParticleHPFinalState.

Definition at line 84 of file G4ParticleHPInelasticCompFS.hh.

84 { return theXsection[50]; }

G4ParticleHPInelasticCompFS::theXsection

G4ParticleHPVector * theXsection[51]

Definition: G4ParticleHPInelasticCompFS.hh:105

Here is the caller graph for this function:

void G4ParticleHPInelasticCompFS::Init	(	G4double	A,
		G4double	Z,
		G4int	M,
		G4String &	dirName,
		G4String &	aSFType,
		G4ParticleDefinition *
	)

virtual

Implements G4ParticleHPFinalState.

Reimplemented in G4ParticleHPNInelasticFS, G4ParticleHPPInelasticFS, and G4ParticleHPTInelasticFS.

Definition at line 81 of file G4ParticleHPInelasticCompFS.cc.

 {
   gammaPath = "/Inelastic/Gammas/";  //only in neutron data base 
     if(!getenv("G4NEUTRONHPDATA")) 
        throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files where Inelastic/Gammas data is found.");
   G4String tBase = getenv("G4NEUTRONHPDATA");
   gammaPath = tBase+gammaPath;
   G4String tString = dirName;
   G4bool dbool;
   G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aFSType, dbool);
   G4String filename = aFile.GetName();
 #ifdef G4PHPDEBUG
   if( getenv("G4ParticleHPDebug") ) G4cout << " G4ParticleHPInelasticCompFS::Init FILE " << filename << G4endl;
 #endif
 
   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)) )
   {
 #ifdef G4PHPDEBUG
     if(getenv("G4ParticleHPDebug_NamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
 #endif
     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);
    G4ParticleHPManager::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) // Loop checking, 11.05.2015, T. Koi
   {
     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*CLHEP::eV;
       LR[ it ] = ilr;
       theXsection[it] = new G4ParticleHPVector;
       G4int total;
       theData >> total;
       theXsection[it]->Init(theData, total, CLHEP::eV);
       //std::cout << theXsection[it]->GetXsec(1*MeV) << std::endl;
     }
     else if(dataType==4)
     {
       theAngularDistribution[it] = new G4ParticleHPAngular;
       theAngularDistribution[it]->Init(theData);
     }
     else if(dataType==5)
     {
       theEnergyDistribution[it] = new G4ParticleHPEnergyDistribution;
       theEnergyDistribution[it]->Init(theData); 
     }
     else if(dataType==6)
     {
       theEnergyAngData[it] = new G4ParticleHPEnAngCorrelation(theProjectile);
       //      G4cout << this << " CompFS theEnergyAngData " << it << theEnergyAngData[it] << G4endl; //GDEB
       theEnergyAngData[it]->Init(theData);
     }
     else if(dataType==12)
     {
       theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
       theFinalStatePhotons[it]->InitMean(theData);
     }
     else if(dataType==13)
     {
       theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
       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 G4ParticleHPInelasticCompFS");
     }
   }
   //  theData.close();
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPInelasticCompFS::InitDistributionInitialState	(	G4ReactionProduct &	anIncidentPart,
		G4ReactionProduct &	aTarget,
		G4int	it
	)

inline

Definition at line 87 of file G4ParticleHPInelasticCompFS.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 59 of file G4ParticleHPInelasticCompFS.cc.

 {
   //   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;
 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

virtual G4ParticleHPFinalState* G4ParticleHPInelasticCompFS::New ( )

pure virtual

Implements G4ParticleHPFinalState.

Implemented in G4ParticleHPAInelasticFS, G4ParticleHPDInelasticFS, G4ParticleHPHe3InelasticFS, G4ParticleHPNInelasticFS, G4ParticleHPPInelasticFS, and G4ParticleHPTInelasticFS.

G4int G4ParticleHPInelasticCompFS::SelectExitChannel ( G4double eKinetic )

Definition at line 195 of file G4ParticleHPInelasticCompFS.cc.

 {
 
 // 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;
       //       G4cout << " SelectExitChannel " << it << " " << random << " " << running[i0]/sum << " " << running[i0] << G4endl; //GDEB
       if(random < running[i0]/sum) break;
     }
   }
 //debug:  it = 1;
 //  G4cout << " SelectExitChannel " << it << " " << sum << G4endl; //GDEB
   return it;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

G4String G4ParticleHPInelasticCompFS::gammaPath

protected

Definition at line 113 of file G4ParticleHPInelasticCompFS.hh.

std::vector<G4int > G4ParticleHPInelasticCompFS::LR

protected

Definition at line 120 of file G4ParticleHPInelasticCompFS.hh.

std::vector< G4double > G4ParticleHPInelasticCompFS::QI

protected

Definition at line 119 of file G4ParticleHPInelasticCompFS.hh.

G4ParticleHPAngular* G4ParticleHPInelasticCompFS::theAngularDistribution[51]

protected

Definition at line 107 of file G4ParticleHPInelasticCompFS.hh.

G4ParticleHPEnAngCorrelation* G4ParticleHPInelasticCompFS::theEnergyAngData[51]

protected

Definition at line 108 of file G4ParticleHPInelasticCompFS.hh.

G4ParticleHPEnergyDistribution* G4ParticleHPInelasticCompFS::theEnergyDistribution[51]

protected

Definition at line 106 of file G4ParticleHPInelasticCompFS.hh.

G4ParticleHPPhotonDist* G4ParticleHPInelasticCompFS::theFinalStatePhotons[51]

protected

Definition at line 110 of file G4ParticleHPInelasticCompFS.hh.

G4ParticleHPDeExGammas G4ParticleHPInelasticCompFS::theGammas

protected

Definition at line 112 of file G4ParticleHPInelasticCompFS.hh.

G4ParticleHPVector* G4ParticleHPInelasticCompFS::theXsection[51]

protected

Definition at line 105 of file G4ParticleHPInelasticCompFS.hh.

The documentation for this class was generated from the following files:

geant4.10.03.p01/source/processes/hadronic/models/particle_hp/include/G4ParticleHPInelasticCompFS.hh
geant4.10.03.p01/source/processes/hadronic/models/particle_hp/src/G4ParticleHPInelasticCompFS.cc

Public Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation