G4NeutronHPDiscreteTwoBody Class Reference

#include <G4NeutronHPDiscreteTwoBody.hh>

Inheritance diagram for G4NeutronHPDiscreteTwoBody:

G4VNeutronHPEnergyAngular

Public Member Functions

 G4NeutronHPDiscreteTwoBody ()
 ~G4NeutronHPDiscreteTwoBody ()
void Init (std::ifstream &aDataFile)
G4ReactionProduct * Sample (G4double anEnergy, G4double massCode, G4double mass)
G4double MeanEnergyOfThisInteraction ()

Detailed Description

Definition at line 42 of file G4NeutronHPDiscreteTwoBody.hh.

Constructor & Destructor Documentation

G4NeutronHPDiscreteTwoBody::G4NeutronHPDiscreteTwoBody (  )  [inline]

Definition at line 46 of file G4NeutronHPDiscreteTwoBody.hh.

00047   {
00048     theCoeff = 0;
00049   }

G4NeutronHPDiscreteTwoBody::~G4NeutronHPDiscreteTwoBody (  )  [inline]

Definition at line 50 of file G4NeutronHPDiscreteTwoBody.hh.

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

Member Function Documentation

void G4NeutronHPDiscreteTwoBody::Init ( std::ifstream &  aDataFile  )  [inline, virtual]

Implements G4VNeutronHPEnergyAngular.

Definition at line 55 of file G4NeutronHPDiscreteTwoBody.hh.

References G4NeutronHPLegendreTable::Init(), G4InterpolationManager::Init(), G4NeutronHPLegendreTable::SetCoeff(), and G4NeutronHPLegendreTable::SetRepresentation().

00056   {
00057     aDataFile >> nEnergy;
00058     theManager.Init(aDataFile);
00059     theCoeff = new G4NeutronHPLegendreTable[nEnergy];
00060     for(G4int i=0; i<nEnergy; i++)
00061     {
00062       G4double energy;
00063       G4int aRep, nCoeff;
00064       aDataFile >> energy >> aRep >> nCoeff;
00065       energy*=CLHEP::eV;
00066       G4int nPoints=nCoeff;
00067       if(aRep>0) nPoints*=2;
00068       //theCoeff[i].Init(energy, nPoints);
00069 
00070       theCoeff[i].Init(energy, nPoints-1);
00071       theCoeff[i].SetRepresentation(aRep);
00072       for(G4int ii=0; ii<nPoints; ii++)
00073       {
00074         G4double y;
00075         aDataFile >> y;
00076         theCoeff[i].SetCoeff(ii, y);
00077       }
00078     }
00079   }

G4double G4NeutronHPDiscreteTwoBody::MeanEnergyOfThisInteraction (  )  [inline, virtual]

Implements G4VNeutronHPEnergyAngular.

Definition at line 82 of file G4NeutronHPDiscreteTwoBody.hh.

00082 { return -1; }

G4ReactionProduct * G4NeutronHPDiscreteTwoBody::Sample ( G4double  anEnergy,
G4double  massCode,
G4double  mass 
) [virtual]

Implements G4VNeutronHPEnergyAngular.

Definition at line 48 of file G4NeutronHPDiscreteTwoBody.cc.

References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4Electron::Electron(), G4UniformRand, G4Gamma::Gamma(), G4NeutronHPLegendreTable::GetCoeff(), G4NeutronHPLegendreTable::GetEnergy(), G4ReactionProduct::GetMass(), G4VNeutronHPEnergyAngular::GetNeutron(), G4NeutronHPLegendreTable::GetNumberOfPoly(), G4VNeutronHPEnergyAngular::GetQValue(), G4InterpolationManager::GetScheme(), G4VNeutronHPEnergyAngular::GetTarget(), G4ReactionProduct::GetTotalMomentum(), G4NeutronHPVector::GetVectorLength(), G4NeutronHPVector::GetX(), G4NeutronHPVector::GetY(), G4He3::He3(), G4InterpolationManager::Init(), G4NeutronHPInterpolator::Interpolate(), LINLIN, LOGLIN, G4NeutronHPVector::Merge(), G4Neutron::Neutron(), G4Positron::Positron(), G4Proton::Proton(), G4NeutronHPVector::Sample(), G4NeutronHPLegendreStore::SampleDiscreteTwoBody(), G4NeutronHPLegendreStore::SetCoeff(), G4NeutronHPVector::SetData(), G4ReactionProduct::SetDefinition(), G4NeutronHPVector::SetInterpolationManager(), G4ReactionProduct::SetKineticEnergy(), G4NeutronHPLegendreStore::SetManager(), G4ReactionProduct::SetMomentum(), G4NeutronHPVector::SetX(), G4NeutronHPVector::SetY(), and G4Triton::Triton().

00049 { // Interpolation still only for the most used parts; rest to be Done @@@@@
00050    G4ReactionProduct * result = new G4ReactionProduct;
00051    G4int Z = static_cast<G4int>(massCode/1000);
00052    G4int A = static_cast<G4int>(massCode-1000*Z);
00053 
00054    if(massCode==0)
00055    {
00056      result->SetDefinition(G4Gamma::Gamma());
00057    }
00058    else if(A==0)
00059    {
00060      result->SetDefinition(G4Electron::Electron());     
00061      if(Z==1) result->SetDefinition(G4Positron::Positron());
00062    }
00063    else if(A==1)
00064    {
00065      result->SetDefinition(G4Neutron::Neutron());
00066      if(Z==1) result->SetDefinition(G4Proton::Proton());
00067    }
00068    else if(A==2)
00069    {
00070      result->SetDefinition(G4Deuteron::Deuteron());      
00071    }
00072    else if(A==3)
00073    {
00074      result->SetDefinition(G4Triton::Triton());  
00075      if(Z==2) result->SetDefinition(G4He3::He3());
00076    }
00077    else if(A==4)
00078    {
00079      result->SetDefinition(G4Alpha::Alpha());
00080      if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");    
00081    }
00082    else
00083    {
00084      throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPDiscreteTwoBody: Unknown ion case 2");
00085    }
00086    
00087 // get cosine(theta)
00088    G4int i(0), it(0);
00089    G4double cosTh(0);
00090    for(i=0; i<nEnergy; i++)
00091    {
00092      it = i;
00093      if(theCoeff[i].GetEnergy()>anEnergy) break;
00094    }
00095    if(it==0||it==nEnergy-1)
00096    {
00097      if(theCoeff[it].GetRepresentation()==0)
00098      {
00099 //TK Legendre expansion
00100        G4NeutronHPLegendreStore theStore(1);
00101        theStore.SetCoeff(0, theCoeff);
00102        theStore.SetManager(theManager);
00103        //cosTh = theStore.SampleMax(anEnergy);
00104        //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
00105        cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
00106      }
00107      else if(theCoeff[it].GetRepresentation()==12) // means LINLIN
00108      {
00109        G4NeutronHPVector theStore; 
00110        G4InterpolationManager aManager;
00111        aManager.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
00112        theStore.SetInterpolationManager(aManager);
00113        for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
00114        {
00115          //101110
00116          //theStore.SetX(i, theCoeff[it].GetCoeff(i));
00117          //theStore.SetY(i, theCoeff[it].GetCoeff(i));
00118          theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
00119          theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
00120          i++;
00121        }
00122        cosTh = theStore.Sample();
00123      }
00124      else if(theCoeff[it].GetRepresentation()==14) //this is LOGLIN
00125      {
00126        G4NeutronHPVector theStore;
00127        G4InterpolationManager aManager;
00128        aManager.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
00129        theStore.SetInterpolationManager(aManager);
00130        for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
00131        {
00132          //101110
00133          //theStore.SetX(i, theCoeff[it].GetCoeff(i));
00134          //theStore.SetY(i, theCoeff[it].GetCoeff(i));
00135          theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
00136          theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
00137          i++;
00138        }
00139        cosTh = theStore.Sample(); 
00140      }
00141      else
00142      {
00143        throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering");
00144      }
00145    }
00146    else
00147    {
00148      if(theCoeff[it].GetRepresentation() == theCoeff[it-1].GetRepresentation())
00149      {
00150        if(theCoeff[it].GetRepresentation()==0)
00151        {
00152 //TK Legendre expansion
00153          G4NeutronHPLegendreStore theStore(2);
00154          theStore.SetCoeff(0, &(theCoeff[it-1]));
00155          theStore.SetCoeff(1, &(theCoeff[it]));
00156          G4InterpolationManager aManager;
00157          aManager.Init(theManager.GetScheme(it), 2);
00158          theStore.SetManager(aManager);
00159          //cosTh = theStore.SampleMax(anEnergy);
00160 //080709 TKDB
00161          cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
00162        }
00163        else if(theCoeff[it].GetRepresentation()==12) // LINLIN
00164        {
00165          G4NeutronHPVector theBuff1;
00166          G4InterpolationManager aManager1;
00167          aManager1.Init(LINLIN, theCoeff[it-1].GetNumberOfPoly()/2);
00168          theBuff1.SetInterpolationManager(aManager1);
00169          for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i++)
00170          {
00171            //101110
00172            //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
00173            //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
00174            theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
00175            theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
00176            i++;
00177          }
00178          G4NeutronHPVector theBuff2;
00179          G4InterpolationManager aManager2;
00180          aManager2.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
00181          theBuff2.SetInterpolationManager(aManager2);
00182          for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
00183          {
00184            //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
00185            //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
00186            theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
00187            theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
00188            i++;
00189          }
00190 
00191          G4double x1 = theCoeff[it-1].GetEnergy();
00192          G4double x2 = theCoeff[it].GetEnergy();
00193          G4double x = anEnergy;
00194          G4double y1, y2, y, mu;
00195 
00196          G4NeutronHPVector theStore1;
00197          theStore1.SetInterpolationManager(aManager1);
00198          G4NeutronHPVector theStore2;
00199          theStore2.SetInterpolationManager(aManager2);
00200          G4NeutronHPVector theStore;
00201          
00202          // for fixed mu get p1, p2 and interpolate according to x
00203          for(i=0; i<theBuff1.GetVectorLength(); i++)
00204          {
00205            mu = theBuff1.GetX(i);
00206            y1 = theBuff1.GetY(i);
00207            y2 = theBuff2.GetY(mu);
00208            y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00209            theStore1.SetData(i, mu, y);
00210          }
00211          for(i=0; i<theBuff2.GetVectorLength(); i++)
00212          {
00213            mu = theBuff2.GetX(i);
00214            y1 = theBuff2.GetY(i);
00215            y2 = theBuff1.GetY(mu);
00216            y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00217            theStore2.SetData(i, mu, y);
00218          }
00219          theStore.Merge(&theStore1, &theStore2); // merge takes care of interpolationschemes
00220          cosTh = theStore.Sample();
00221        }
00222        else if(theCoeff[it].GetRepresentation()==14) //TK LOG_LIN
00223        {
00224          G4NeutronHPVector theBuff1;
00225          G4InterpolationManager aManager1;
00226          aManager1.Init(LOGLIN, theCoeff[it-1].GetNumberOfPoly()/2);
00227          theBuff1.SetInterpolationManager(aManager1);
00228          for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i++)
00229          {
00230            //101110
00231            //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
00232            //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
00233            theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
00234            theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
00235            i++;
00236          }
00237          
00238          G4NeutronHPVector theBuff2;
00239          G4InterpolationManager aManager2;
00240          aManager2.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
00241          theBuff2.SetInterpolationManager(aManager2);
00242          for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
00243          {
00244            //101110
00245            //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
00246            //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
00247            theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
00248            theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
00249            i++;
00250          }
00251 
00252          G4double x1 = theCoeff[it-1].GetEnergy();
00253          G4double x2 = theCoeff[it].GetEnergy();
00254          G4double x = anEnergy;
00255          G4double y1, y2, y, mu;
00256 
00257          G4NeutronHPVector theStore1;
00258          theStore1.SetInterpolationManager(aManager1);
00259          G4NeutronHPVector theStore2;
00260          theStore2.SetInterpolationManager(aManager2);
00261          G4NeutronHPVector theStore;
00262          
00263          // for fixed mu get p1, p2 and interpolate according to x
00264          for(i=0; i<theBuff1.GetVectorLength(); i++)
00265          {
00266            mu = theBuff1.GetX(i);
00267            y1 = theBuff1.GetY(i);
00268            y2 = theBuff2.GetY(mu);
00269            y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00270            theStore1.SetData(i, mu, y);
00271          }
00272          for(i=0; i<theBuff2.GetVectorLength(); i++)
00273          {
00274            mu = theBuff2.GetX(i);
00275            y1 = theBuff2.GetY(i);
00276            y2 = theBuff1.GetY(mu);
00277            y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00278            theStore2.SetData(i, mu, y);
00279          }
00280          theStore.Merge(&theStore1, &theStore2); 
00281          cosTh = theStore.Sample(); 
00282        }
00283        else
00284        {
00285          throw G4HadronicException(__FILE__, __LINE__, "Two neighbouring distributions with different interpolation");
00286        }
00287      }
00288      else
00289      {
00290        throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering, case 2");
00291      }
00292    }
00293    
00294 // now get the energy from kinematics and Q-value.
00295 
00296    //G4double restEnergy = anEnergy+GetQValue();
00297    
00298 // assumed to be in CMS @@@@@@@@@@@@@@@@@
00299 
00300    //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2
00301    //G4double residualMass =   GetTarget()->GetMass() + GetNeutron()->GetMass()
00302    //                        - result->GetMass() - GetQValue();
00303    //G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
00304    G4double A1     =  GetTarget()->GetMass()/GetNeutron()->GetMass(); 
00305    G4double A1prim =  result->GetMass()/GetNeutron()->GetMass();
00306    G4double E1     =  (A1+1)*(A1+1)/A1/A1*anEnergy; 
00307    G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*E1+(1+A1)*GetQValue());
00308 
00309    result->SetKineticEnergy(kinE); // non relativistic @@
00310    G4double phi = twopi*G4UniformRand();
00311    G4double theta = std::acos(cosTh);
00312    G4double sinth = std::sin(theta);
00313    G4double mtot = result->GetTotalMomentum(); 
00314    G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
00315    result->SetMomentum(tempVector);
00316    
00317 // some garbage collection
00318    
00319 // return the result   
00320    return result;
00321 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:37 2013 for Geant4 by  doxygen 1.4.7