G4ParticleHPDiscreteTwoBody Class Reference

#include <G4ParticleHPDiscreteTwoBody.hh>

Inheritance diagram for G4ParticleHPDiscreteTwoBody:

Public Member Functions
	G4ParticleHPDiscreteTwoBody ()
	~G4ParticleHPDiscreteTwoBody ()
void	Init (std::istream &aDataFile)
G4ReactionProduct *	Sample (G4double anEnergy, G4double massCode, G4double mass)
G4double	MeanEnergyOfThisInteraction ()
Public Member Functions inherited from G4VParticleHPEnergyAngular
	G4VParticleHPEnergyAngular ()
virtual	~G4VParticleHPEnergyAngular ()
virtual void	Init (std::istream &aDataFile)=0
virtual G4ReactionProduct *	Sample (G4double anEnergy, G4double massCode, G4double mass)=0
virtual G4double	MeanEnergyOfThisInteraction ()=0
void	SetProjectileRP (G4ReactionProduct *aIncidentParticleRP)
void	SetTarget (G4ReactionProduct *aTarget)
G4ReactionProduct *	GetTarget ()
G4ReactionProduct *	GetProjectileRP ()
G4ReactionProduct *	GetCMS ()
void	SetQValue (G4double aValue)
virtual void	ClearHistories ()

Additional Inherited Members
Protected Member Functions inherited from G4VParticleHPEnergyAngular
G4double	GetQValue ()

Detailed Description

Definition at line 42 of file G4ParticleHPDiscreteTwoBody.hh.

Constructor & Destructor Documentation

G4ParticleHPDiscreteTwoBody()

G4ParticleHPDiscreteTwoBody::G4ParticleHPDiscreteTwoBody ( )

inline

Definition at line 46 of file G4ParticleHPDiscreteTwoBody.hh.

  {
    theCoeff = 0;
    bCheckDiffCoeffRepr = true;
    if ( std::getenv( "G4PHP_DO_NOT_CHECK_DIFF_COEFF_REPR" ) ) bCheckDiffCoeffRepr = false;
    nEnergy = 0;
  }

~G4ParticleHPDiscreteTwoBody()

G4ParticleHPDiscreteTwoBody::~G4ParticleHPDiscreteTwoBody ( )

inline

Definition at line 53 of file G4ParticleHPDiscreteTwoBody.hh.

  {
    if(theCoeff!=0) delete [] theCoeff;
  }

Member Function Documentation

Init()

void G4ParticleHPDiscreteTwoBody::Init ( std::istream &  aDataFile )

inlinevirtual

Implements G4VParticleHPEnergyAngular.

Definition at line 58 of file G4ParticleHPDiscreteTwoBody.hh.

  {
    aDataFile >> nEnergy;
    theManager.Init(aDataFile);
    theCoeff = new G4ParticleHPLegendreTable[nEnergy];
    for(G4int i=0; i<nEnergy; i++)
    {
      G4double energy;
      G4int aRep, nCoeff;
      aDataFile >> energy >> aRep >> nCoeff;
      //-      G4cout << this << " " << i << " G4ParticleHPDiscreteTwoBody READ DATA " << energy << " " << aRep << " " << nCoeff << G4endl;
      energy*=CLHEP::eV;
      G4int nPoints=nCoeff;
      if(aRep>0) nPoints*=2;
      //theCoeff[i].Init(energy, nPoints);
 
      theCoeff[i].Init(energy, nPoints-1);
      theCoeff[i].SetRepresentation(aRep);
      for(G4int ii=0; ii<nPoints; ii++)
      {
        G4double y;
        aDataFile >> y;
        theCoeff[i].SetCoeff(ii, y);
       }
    }
  }

MeanEnergyOfThisInteraction()

G4double G4ParticleHPDiscreteTwoBody::MeanEnergyOfThisInteraction ( )

inlinevirtual

Implements G4VParticleHPEnergyAngular.

Definition at line 86 of file G4ParticleHPDiscreteTwoBody.hh.

{ return -1; }

Sample()

G4ReactionProduct * G4ParticleHPDiscreteTwoBody::Sample ( G4double  anEnergy, G4double  massCode, G4double  mass  )

virtual

Implements G4VParticleHPEnergyAngular.

Definition at line 49 of file G4ParticleHPDiscreteTwoBody.cc.

{ // Interpolation still only for the most used parts; rest to be Done @@@@@
   G4ReactionProduct * result = new G4ReactionProduct;
   G4int Z = static_cast<G4int>(massCode/1000);
   G4int A = static_cast<G4int>(massCode-1000*Z);
 
   if(massCode==0)
   {
     result->SetDefinition(G4Gamma::Gamma());
   }
   else if(A==0)
   {
     result->SetDefinition(G4Electron::Electron());     
     if(Z==1) result->SetDefinition(G4Positron::Positron());
   }
   else if(A==1)
   {
     result->SetDefinition(G4Neutron::Neutron());
     if(Z==1) result->SetDefinition(G4Proton::Proton());
   }
   else if(A==2)
   {
     result->SetDefinition(G4Deuteron::Deuteron());      
   }
   else if(A==3)
   {
     result->SetDefinition(G4Triton::Triton());  
     if(Z==2) result->SetDefinition(G4He3::He3());
   }
   else if(A==4)
   {
     result->SetDefinition(G4Alpha::Alpha());
     if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");    
   }
   else
   {
     throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPDiscreteTwoBody: Unknown ion case 2");
   }
   
// get cosine(theta)
   G4int i(0), it(0);
   G4double cosTh(0);
   for(i=0; i<nEnergy; i++)
   {
     it = i;
     if(theCoeff[i].GetEnergy()>anEnergy) break;
   }
   if(it==0||it==nEnergy-1)
   {
     if(theCoeff[it].GetRepresentation()==0)
     {
//TK Legendre expansion
       G4ParticleHPLegendreStore theStore(1);
       theStore.SetCoeff(0, theCoeff);
       theStore.SetManager(theManager);
       //cosTh = theStore.SampleMax(anEnergy);
       //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
       cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
     }
     else if(theCoeff[it].GetRepresentation()==12) // means LINLIN
     {
       G4ParticleHPVector theStore; 
       G4InterpolationManager aManager;
       aManager.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
       theStore.SetInterpolationManager(aManager);
       for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
       {
         //101110
         //theStore.SetX(i, theCoeff[it].GetCoeff(i));
         //theStore.SetY(i, theCoeff[it].GetCoeff(i));
         theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
         theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
       }
       cosTh = theStore.Sample();
     }
     else if(theCoeff[it].GetRepresentation()==14) //this is LOGLIN
     {
       G4ParticleHPVector theStore;
       G4InterpolationManager aManager;
       aManager.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
       theStore.SetInterpolationManager(aManager);
       for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
       {
         //101110
         //theStore.SetX(i, theCoeff[it].GetCoeff(i));
         //theStore.SetY(i, theCoeff[it].GetCoeff(i));
         theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
         theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
       }
       cosTh = theStore.Sample(); 
     }
     else
     {
       throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering");
     }
   }
   else
   {
     if(!bCheckDiffCoeffRepr || theCoeff[it].GetRepresentation() == theCoeff[it-1].GetRepresentation())
     {
       if(theCoeff[it].GetRepresentation()==0)
       {
//TK Legendre expansion
     G4ParticleHPLegendreStore theStore(2);
     theStore.SetCoeff(0, &(theCoeff[it-1]));
     theStore.SetCoeff(1, &(theCoeff[it]));
         G4InterpolationManager aManager;
         aManager.Init(theManager.GetScheme(it), 2);
         theStore.SetManager(aManager);
     //cosTh = theStore.SampleMax(anEnergy);
//080709 TKDB
         cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
       }
       else if(theCoeff[it].GetRepresentation()==12) // LINLIN
       {
     G4ParticleHPVector theBuff1;
         G4InterpolationManager aManager1;
         aManager1.Init(LINLIN, theCoeff[it-1].GetNumberOfPoly()/2);
         theBuff1.SetInterpolationManager(aManager1);
     for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i+=2)
     {
           //101110
           //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
           //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
           theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
           theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
     }
     G4ParticleHPVector theBuff2;
         G4InterpolationManager aManager2;
         aManager2.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
         theBuff2.SetInterpolationManager(aManager2);
     for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
     {
           //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
           //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
           theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
           theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
     }
 
     G4double x1 = theCoeff[it-1].GetEnergy();
     G4double x2 = theCoeff[it].GetEnergy();
     G4double x = anEnergy;
     G4double y1, y2, y, mu;
 
     G4ParticleHPVector theStore1;
         theStore1.SetInterpolationManager(aManager1);
     G4ParticleHPVector theStore2;
         theStore2.SetInterpolationManager(aManager2);
     G4ParticleHPVector theStore;
     
     // for fixed mu get p1, p2 and interpolate according to x
         for(i=0; i<theBuff1.GetVectorLength(); i++)
         {
           mu = theBuff1.GetX(i);
           y1 = theBuff1.GetY(i);
           y2 = theBuff2.GetY(mu);
           y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
           theStore1.SetData(i, mu, y);
         }
         for(i=0; i<theBuff2.GetVectorLength(); i++)
         {
           mu = theBuff2.GetX(i);
           y1 = theBuff2.GetY(i);
           y2 = theBuff1.GetY(mu);
           y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
           theStore2.SetData(i, mu, y);
         }
         theStore.Merge(&theStore1, &theStore2); // merge takes care of interpolationschemes
     cosTh = theStore.Sample();
       }
       else if(theCoeff[it].GetRepresentation()==14) //TK LOG_LIN
       {
     G4ParticleHPVector theBuff1;
         G4InterpolationManager aManager1;
         aManager1.Init(LOGLIN, theCoeff[it-1].GetNumberOfPoly()/2);
         theBuff1.SetInterpolationManager(aManager1);
     for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i+=2)
     {
           //101110
           //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
           //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
           theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
           theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
     }
     
     G4ParticleHPVector theBuff2;
         G4InterpolationManager aManager2;
         aManager2.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
         theBuff2.SetInterpolationManager(aManager2);
     for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
     {
           //101110
           //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
           //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
           theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
           theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
     }
 
     G4double x1 = theCoeff[it-1].GetEnergy();
     G4double x2 = theCoeff[it].GetEnergy();
     G4double x = anEnergy;
     G4double y1, y2, y, mu;
 
     G4ParticleHPVector theStore1;
         theStore1.SetInterpolationManager(aManager1);
     G4ParticleHPVector theStore2;
         theStore2.SetInterpolationManager(aManager2);
     G4ParticleHPVector theStore;
     
     // for fixed mu get p1, p2 and interpolate according to x
         for(i=0; i<theBuff1.GetVectorLength(); i++)
         {
           mu = theBuff1.GetX(i);
           y1 = theBuff1.GetY(i);
           y2 = theBuff2.GetY(mu);
           y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
           theStore1.SetData(i, mu, y);
         }
         for(i=0; i<theBuff2.GetVectorLength(); i++)
         {
           mu = theBuff2.GetX(i);
           y1 = theBuff2.GetY(i);
           y2 = theBuff1.GetY(mu);
           y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
           theStore2.SetData(i, mu, y);
         }
         theStore.Merge(&theStore1, &theStore2); 
     cosTh = theStore.Sample(); 
       }
       else
       {
         throw G4HadronicException(__FILE__, __LINE__, "Two neighbouring distributions with different interpolation");
       }
     }
     else
     {
       G4cout << " theCoeff[it].GetRepresent MEM " << it << " " << &theCoeff[it] <<  "  " << &theCoeff[it-1] << G4endl;
       G4cout << " theCoeff[it].GetRepresent " << it << " " << theCoeff[it].GetRepresentation() <<  " != " << theCoeff[it-1].GetRepresentation() << G4endl;
 
       throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering, case 2");
     }
   }
   
// now get the energy from kinematics and Q-value.
 
   //G4double restEnergy = anEnergy+GetQValue();
   
// assumed to be in CMS @@@@@@@@@@@@@@@@@
 
   //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2
   //G4double residualMass =   GetTarget()->GetMass() + GetNeutron()->GetMass()
   //                        - result->GetMass() - GetQValue();
   //G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
   G4double A1     =  GetTarget()->GetMass()/GetProjectileRP()->GetMass(); 
   G4double A1prim =  result->GetMass()/GetProjectileRP()->GetMass();
   //G4double E1     =  (A1+1)*(A1+1)/A1/A1*anEnergy; 
   //Bug fix Bugzilla #1815
   G4double E1     =  anEnergy; 
   G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*E1+(1+A1)*GetQValue());
 
   result->SetKineticEnergy(kinE); // non relativistic @@
   G4double phi = CLHEP::twopi*G4UniformRand();
   G4double theta = std::acos(cosTh);
   G4double sinth = std::sin(theta);
   G4double mtot = result->GetTotalMomentum(); 
   G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
   result->SetMomentum(tempVector);
   
// some garbage collection
   
// return the result   
   return result;
}

---

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

• G4ParticleHPDiscreteTwoBody.hh
• G4ParticleHPDiscreteTwoBody.cc