G4ParticleHPAngular Class Reference

#include <G4ParticleHPAngular.hh>

Collaboration diagram for G4ParticleHPAngular:

Classes
struct	toBeCached

Public Member Functions
	G4ParticleHPAngular ()
	~G4ParticleHPAngular ()
void	Init (std::istream &aDataFile)
void	SampleAndUpdate (G4ReactionProduct &anIncidentParticle)
void	SetTarget (const G4ReactionProduct &aTarget)
void	SetProjectileRP (const G4ReactionProduct &anIncidentParticleRP)
G4double	GetTargetMass ()

Private Attributes
G4int	theAngularDistributionType
G4int	frameFlag
G4bool	theIsoFlag
G4ParticleHPLegendreStore *	theCoefficients
G4ParticleHPPartial *	theProbArray
G4double	targetMass
G4Cache< toBeCached >	fCache

Detailed Description

Definition at line 44 of file G4ParticleHPAngular.hh.

Constructor & Destructor Documentation

G4ParticleHPAngular()

G4ParticleHPAngular::G4ParticleHPAngular ( )

inline

Definition at line 55 of file G4ParticleHPAngular.hh.

  {
//TKDB110511
    //theAngularDistributionType = 0;
    //theIsoFlag = false;
    theIsoFlag = true;
// TKDB
      theCoefficients = 0;
      theProbArray = 0;
      toBeCached val;
      fCache.Put( val );
  } 

~G4ParticleHPAngular()

G4ParticleHPAngular::~G4ParticleHPAngular ( )

inline

Definition at line 68 of file G4ParticleHPAngular.hh.

   {
// TKDB
      delete theCoefficients;
      delete theProbArray;
   }

Here is the call graph for this function:

Member Function Documentation

GetTargetMass()

G4double G4ParticleHPAngular::GetTargetMass ( )

inline

Definition at line 83 of file G4ParticleHPAngular.hh.

{ return targetMass; }

Init()

void G4ParticleHPAngular::Init

(

std::istream &

aDataFile

)

Definition at line 41 of file G4ParticleHPAngular.cc.

{
//  G4cout << "here we are entering the Angular Init"<<G4endl;
  aDataFile >> theAngularDistributionType >> targetMass;
  aDataFile >> frameFlag;
  if(theAngularDistributionType == 0 )
  {
    theIsoFlag = true; 
  }
  else if(theAngularDistributionType==1)
  {
    theIsoFlag = false;
    G4int nEnergy;
    aDataFile >> nEnergy;  
    theCoefficients = new G4ParticleHPLegendreStore(nEnergy);
    theCoefficients->InitInterpolation(aDataFile);
    G4double temp, energy;
    G4int tempdep, nLegendre;
    G4int i, ii;
    for (i=0; i<nEnergy; i++)
    {
      aDataFile >> temp >> energy >> tempdep >> nLegendre;
      energy *=eV;
      theCoefficients->Init(i, energy, nLegendre);
      theCoefficients->SetTemperature(i, temp);
      G4double coeff=0;
      for(ii=0; ii<nLegendre; ii++)
      {
        aDataFile >> coeff;
        theCoefficients->SetCoeff(i, ii+1, coeff);
      }
    }
  }
  else if (theAngularDistributionType==2)
  {
    theIsoFlag = false;
    G4int nEnergy;
    aDataFile >> nEnergy;
    theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
    theProbArray->InitInterpolation(aDataFile);
    G4double temp, energy;
    G4int tempdep;
    for(G4int i=0; i<nEnergy; i++)
    {
      aDataFile >> temp >> energy >> tempdep;
      energy *= eV;
      theProbArray->SetT(i, temp);
      theProbArray->SetX(i, energy);
      theProbArray->InitData(i, aDataFile);
    }
  }
  else
  {
    theIsoFlag = false;
    G4cout << "unknown distribution found for Angular: "<< theAngularDistributionType << G4endl;
    throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
  }    
}

Here is the call graph for this function:

Here is the caller graph for this function:

SampleAndUpdate()

void G4ParticleHPAngular::SampleAndUpdate

(

G4ReactionProduct &

anIncidentParticle

)

Definition at line 100 of file G4ParticleHPAngular.cc.

{

  //********************************************************************
  //EMendoza -> sampling can be isotropic in LAB or in CMS
  /*
  if(theIsoFlag)
  {
//  G4cout << "Angular result "<<aHadron.GetTotalMomentum()<<" ";
// @@@ add code for isotropic emission in CMS.
      G4double costheta = 2.*G4UniformRand()-1;
      G4double theta = std::acos(costheta);
      G4double phi = twopi*G4UniformRand();
      G4double sinth = std::sin(theta);
      G4double en = aHadron.GetTotalMomentum();
      G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
      aHadron.SetMomentum( temp );
      aHadron.Lorentz(aHadron, -1.*theTarget);
  }
  else
  {
  */
  //********************************************************************
    if(frameFlag == 1) // LAB
    {
      G4double en = aHadron.GetTotalMomentum();
      G4ReactionProduct boosted;
      boosted.Lorentz(*fCache.Get().theProjectileRP, *fCache.Get().theTarget);
      G4double kineticEnergy = boosted.GetKineticEnergy();
      G4double cosTh = 0.0; 
      //********************************************************************
      //EMendoza --> sampling can be also isotropic
      /*
      if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy); 
      if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy); 
      */
      //********************************************************************
      if(theIsoFlag){cosTh =2.*G4UniformRand()-1;}
      else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
      else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
      else{ 
    G4cout << "unknown distribution found for Angular: "<< theAngularDistributionType << G4endl;
        throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
      }
      //********************************************************************
      G4double theta = std::acos(cosTh);
      G4double phi = twopi*G4UniformRand();
      G4double sinth = std::sin(theta);
      G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
      aHadron.SetMomentum( temp );
    }
    else if(frameFlag == 2) // costh in CMS
    {
      G4ReactionProduct boostedN;
      boostedN.Lorentz(*fCache.Get().theProjectileRP, *fCache.Get().theTarget);
      G4double kineticEnergy = boostedN.GetKineticEnergy();

      G4double cosTh = 0.0; 
      //********************************************************************
      //EMendoza --> sampling can be also isotropic
      /*
      if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy); 
      if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy); 
      */
      //********************************************************************
      if(theIsoFlag){cosTh =2.*G4UniformRand()-1;}
      else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
      else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
      else{
    G4cout << "unknown distribution found for Angular: "<< theAngularDistributionType << G4endl;
        throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
      }
      //********************************************************************
//080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) 
/*
    if(theAngularDistributionType == 1) // LAB
    {
      G4double en = aHadron.GetTotalMomentum();
      G4ReactionProduct boosted;
      boosted.Lorentz(theProjectile, theTarget);
      G4double kineticEnergy = boosted.GetKineticEnergy();
      G4double cosTh = theCoefficients->SampleMax(kineticEnergy);
      G4double theta = std::acos(cosTh);
      G4double phi = twopi*G4UniformRand();
      G4double sinth = std::sin(theta);
      G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
      aHadron.SetMomentum( temp );
    }
    else if(theAngularDistributionType == 2) // costh in CMS {
    }
*/

//      G4ReactionProduct boostedN;
//      boostedN.Lorentz(theProjectile, theTarget);
//      G4double kineticEnergy = boostedN.GetKineticEnergy();
//      G4double cosTh = theProbArray->Sample(kineticEnergy); 

      G4double theta = std::acos(cosTh);
      G4double phi = twopi*G4UniformRand();
      G4double sinth = std::sin(theta);      
      
      G4ThreeVector temp(sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) ); //CMS

//080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
/*
      G4double en = aHadron.GetTotalEnergy(); // Target rest
      
      // get trafo from Target rest frame to CMS
      G4ReactionProduct boostedT;
      boostedT.Lorentz(theTarget, theTarget);
      
      G4ThreeVector the3IncidentParticle = boostedN.GetMomentum();
      G4double nEnergy = boostedN.GetTotalEnergy();
      G4ThreeVector the3Target = boostedT.GetMomentum();
      G4double tEnergy = boostedT.GetTotalEnergy();
      G4double totE = nEnergy+tEnergy;
      G4ThreeVector the3trafo = -the3Target-the3IncidentParticle;
      G4ReactionProduct trafo; // for transformation from CMS to target rest frame
      trafo.SetMomentum(the3trafo);
      G4double cmsMom = std::sqrt(the3trafo*the3trafo);
      G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
      trafo.SetMass(sqrts);
      trafo.SetTotalEnergy(totE);
      
      G4double gamma = trafo.GetTotalEnergy()/trafo.GetMass();
      G4double cosalpha = temp*trafo.GetMomentum()/trafo.GetTotalMomentum()/temp.mag();
      G4double fac = cosalpha*trafo.GetTotalMomentum()/trafo.GetMass();
      fac*=gamma;
      
      G4double mom;
//    For G4FPE_DEBUG ON
      G4double mom2 = ( en*fac*en*fac - 
                   (fac*fac - gamma*gamma)*
                   (en*en - gamma*gamma*aHadron.GetMass()*aHadron.GetMass())
                );
      if ( mom2 > 0.0 ) 
        mom = std::sqrt( mom2 );
      else
        mom = 0.0; 

      mom = -en*fac - mom;
      mom /= (fac*fac-gamma*gamma);
      temp = mom*temp;
      
      aHadron.SetMomentum( temp ); // now all in CMS
      aHadron.SetTotalEnergy( std::sqrt( mom*mom + aHadron.GetMass()*aHadron.GetMass() ) );
      aHadron.Lorentz(aHadron, trafo); // now in target rest frame
*/
      // Determination of the hadron kinetic energy in CMS
      // aHadron.GetKineticEnergy() is actually the residual kinetic energy in CMS (or target frame)
      // kineticEnergy is incident neutron kinetic energy  in CMS (or target frame)  
      G4double QValue = aHadron.GetKineticEnergy() - kineticEnergy;
      G4double A1     =   fCache.Get().theTarget->GetMass()/boostedN.GetMass(); 
      G4double A1prim =   aHadron.GetMass()/ boostedN.GetMass();
      G4double kinE   = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*kineticEnergy+(1+A1)*QValue);
      G4double totalE = kinE + aHadron.GetMass();
      G4double mom2   = totalE*totalE - aHadron.GetMass()*aHadron.GetMass();
      G4double mom;
      if ( mom2 > 0.0 ) mom = std::sqrt( mom2 );
      else mom = 0.0;     

      aHadron.SetMomentum( mom*temp ); // Set momentum in CMS
      aHadron.SetKineticEnergy(kinE);  // Set kinetic energy in CMS

      // get trafo from Target rest frame to CMS
      G4ReactionProduct boostedT;
      boostedT.Lorentz(*fCache.Get().theTarget, *fCache.Get().theTarget);
      
      G4ThreeVector the3IncidentParticle = boostedN.GetMomentum();
      G4double nEnergy = boostedN.GetTotalEnergy();
      G4ThreeVector the3Target = boostedT.GetMomentum();
      G4double tEnergy = boostedT.GetTotalEnergy();
      G4double totE = nEnergy+tEnergy;
      G4ThreeVector the3trafo = -the3Target-the3IncidentParticle;
      G4ReactionProduct trafo; // for transformation from CMS to target rest frame
      trafo.SetMomentum(the3trafo);
      G4double cmsMom = std::sqrt(the3trafo*the3trafo);
      G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
      trafo.SetMass(sqrts);
      trafo.SetTotalEnergy(totE);

      aHadron.Lorentz(aHadron, trafo);

    }
    else
    {
      throw G4HadronicException(__FILE__, __LINE__, "Tried to sample non isotropic neutron angular");
    }
  aHadron.Lorentz(aHadron, -1.*(*fCache.Get().theTarget)); 
//  G4cout << aHadron.GetMomentum()<<" ";
//  G4cout << aHadron.GetTotalMomentum()<<G4endl;
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetProjectileRP()

void G4ParticleHPAngular::SetProjectileRP ( const G4ReactionProduct &  anIncidentParticleRP )

inline

Definition at line 81 of file G4ParticleHPAngular.hh.

{ fCache.Get().theProjectileRP = &anIncidentParticleRP; }

Here is the caller graph for this function:

SetTarget()

void G4ParticleHPAngular::SetTarget ( const G4ReactionProduct &  aTarget )

inline

Definition at line 79 of file G4ParticleHPAngular.hh.

{ fCache.Get().theTarget = &aTarget; }

Here is the caller graph for this function:

Member Data Documentation

fCache

G4Cache<toBeCached> G4ParticleHPAngular::fCache

private

Definition at line 108 of file G4ParticleHPAngular.hh.

frameFlag

G4int G4ParticleHPAngular::frameFlag

private

Definition at line 94 of file G4ParticleHPAngular.hh.

targetMass

G4double G4ParticleHPAngular::targetMass

private

Definition at line 104 of file G4ParticleHPAngular.hh.

theAngularDistributionType

G4int G4ParticleHPAngular::theAngularDistributionType

private

Definition at line 93 of file G4ParticleHPAngular.hh.

theCoefficients

G4ParticleHPLegendreStore* G4ParticleHPAngular::theCoefficients

private

Definition at line 98 of file G4ParticleHPAngular.hh.

theIsoFlag

G4bool G4ParticleHPAngular::theIsoFlag

private

Definition at line 96 of file G4ParticleHPAngular.hh.

theProbArray

G4ParticleHPPartial* G4ParticleHPAngular::theProbArray

private

Definition at line 100 of file G4ParticleHPAngular.hh.

---

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

• G4ParticleHPAngular.hh
• G4ParticleHPAngular.cc