G4NeutronHPAngular Class Reference

#include <G4NeutronHPAngular.hh>


Public Member Functions

 G4NeutronHPAngular ()
 ~G4NeutronHPAngular ()
void Init (std::ifstream &aDataFile)
void SampleAndUpdate (G4ReactionProduct &aNeutron)
void SetTarget (const G4ReactionProduct &aTarget)
void SetNeutron (const G4ReactionProduct &aNeutron)
G4double GetTargetMass ()


Detailed Description

Definition at line 42 of file G4NeutronHPAngular.hh.


Constructor & Destructor Documentation

G4NeutronHPAngular::G4NeutronHPAngular (  )  [inline]

Definition at line 46 of file G4NeutronHPAngular.hh.

00047   {
00048 //TKDB110511
00049     //theAngularDistributionType = 0;
00050     //theIsoFlag = false;
00051     theIsoFlag = true;
00052 // TKDB
00053       theCoefficients = 0;
00054       theProbArray = 0;
00055   } 

G4NeutronHPAngular::~G4NeutronHPAngular (  )  [inline]

Definition at line 57 of file G4NeutronHPAngular.hh.

00058    {
00059 // TKDB
00060       delete theCoefficients;
00061       delete theProbArray;
00062    }


Member Function Documentation

G4double G4NeutronHPAngular::GetTargetMass (  )  [inline]

Definition at line 72 of file G4NeutronHPAngular.hh.

Referenced by G4NeutronHPInelasticBaseFS::BaseApply().

00072 { return targetMass; }

void G4NeutronHPAngular::Init ( std::ifstream &  aDataFile  ) 

Definition at line 39 of file G4NeutronHPAngular.cc.

References G4cout, G4endl, G4NeutronHPLegendreStore::Init(), G4NeutronHPPartial::InitData(), G4NeutronHPPartial::InitInterpolation(), G4NeutronHPLegendreStore::InitInterpolation(), G4NeutronHPLegendreStore::SetCoeff(), G4NeutronHPPartial::SetT(), G4NeutronHPLegendreStore::SetTemperature(), and G4NeutronHPPartial::SetX().

Referenced by G4NeutronHPInelasticCompFS::Init(), G4NeutronHPInelasticBaseFS::Init(), G4NeutronHPFSFissionFS::Init(), G4NeutronHPFissionBaseFS::Init(), and G4FissionLibrary::Init().

00040 {
00041 //  G4cout << "here we are entering the Angular Init"<<G4endl;
00042   aDataFile >> theAngularDistributionType >> targetMass;
00043   aDataFile >> frameFlag;
00044   if(theAngularDistributionType == 0 )
00045   {
00046     theIsoFlag = true; 
00047   }
00048   else if(theAngularDistributionType==1)
00049   {
00050     theIsoFlag = false;
00051     G4int nEnergy;
00052     aDataFile >> nEnergy;  
00053     theCoefficients = new G4NeutronHPLegendreStore(nEnergy);
00054     theCoefficients->InitInterpolation(aDataFile);
00055     G4double temp, energy;
00056     G4int tempdep, nLegendre;
00057     G4int i, ii;
00058     for (i=0; i<nEnergy; i++)
00059     {
00060       aDataFile >> temp >> energy >> tempdep >> nLegendre;
00061       energy *=eV;
00062       theCoefficients->Init(i, energy, nLegendre);
00063       theCoefficients->SetTemperature(i, temp);
00064       G4double coeff=0;
00065       for(ii=0; ii<nLegendre; ii++)
00066       {
00067         aDataFile >> coeff;
00068         theCoefficients->SetCoeff(i, ii+1, coeff);
00069       }
00070     }
00071   }
00072   else if (theAngularDistributionType==2)
00073   {
00074     theIsoFlag = false;
00075     G4int nEnergy;
00076     aDataFile >> nEnergy;
00077     theProbArray = new G4NeutronHPPartial(nEnergy, nEnergy);
00078     theProbArray->InitInterpolation(aDataFile);
00079     G4double temp, energy;
00080     G4int tempdep;
00081     for(G4int i=0; i<nEnergy; i++)
00082     {
00083       aDataFile >> temp >> energy >> tempdep;
00084       energy *= eV;
00085       theProbArray->SetT(i, temp);
00086       theProbArray->SetX(i, energy);
00087       theProbArray->InitData(i, aDataFile);
00088     }
00089   }
00090   else
00091   {
00092     theIsoFlag = false;
00093     G4cout << "unknown distribution found for Angular"<<G4endl;
00094     throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
00095   }    
00096 }

void G4NeutronHPAngular::SampleAndUpdate ( G4ReactionProduct aNeutron  ) 

Definition at line 98 of file G4NeutronHPAngular.cc.

References G4cout, G4endl, G4UniformRand, G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4ReactionProduct::Lorentz(), G4NeutronHPPartial::Sample(), G4NeutronHPLegendreStore::SampleMax(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), and G4ReactionProduct::SetTotalEnergy().

Referenced by G4NeutronHPFSFissionFS::ApplyYourself(), G4NeutronHPFissionBaseFS::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), and G4NeutronHPInelasticCompFS::CompositeApply().

00099 {
00100 
00101   //********************************************************************
00102   //EMendoza -> sampling can be isotropic in LAB or in CMS
00103   /*
00104   if(theIsoFlag)
00105   {
00106 //  G4cout << "Angular result "<<aHadron.GetTotalMomentum()<<" ";
00107 // @@@ add code for isotropic emission in CMS.
00108       G4double costheta = 2.*G4UniformRand()-1;
00109       G4double theta = std::acos(costheta);
00110       G4double phi = twopi*G4UniformRand();
00111       G4double sinth = std::sin(theta);
00112       G4double en = aHadron.GetTotalMomentum();
00113       G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00114       aHadron.SetMomentum( temp );
00115       aHadron.Lorentz(aHadron, -1.*theTarget);
00116   }
00117   else
00118   {
00119   */
00120   //********************************************************************
00121     if(frameFlag == 1) // LAB
00122     {
00123       G4double en = aHadron.GetTotalMomentum();
00124       G4ReactionProduct boosted;
00125       boosted.Lorentz(theNeutron, theTarget);
00126       G4double kineticEnergy = boosted.GetKineticEnergy();
00127       G4double cosTh = 0.0; 
00128       //********************************************************************
00129       //EMendoza --> sampling can be also isotropic
00130       /*
00131       if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy); 
00132       if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy); 
00133       */
00134       //********************************************************************
00135       if(theIsoFlag){cosTh =2.*G4UniformRand()-1;}
00136       else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
00137       else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
00138       else{
00139         G4cout << "unknown distribution found for Angular"<<G4endl;
00140         throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
00141       }
00142       //********************************************************************
00143       G4double theta = std::acos(cosTh);
00144       G4double phi = twopi*G4UniformRand();
00145       G4double sinth = std::sin(theta);
00146       G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00147       aHadron.SetMomentum( temp );
00148     }
00149     else if(frameFlag == 2) // costh in CMS
00150     {
00151       G4ReactionProduct boostedN;
00152       boostedN.Lorentz(theNeutron, theTarget);
00153       G4double kineticEnergy = boostedN.GetKineticEnergy();
00154 
00155       G4double cosTh = 0.0; 
00156       //********************************************************************
00157       //EMendoza --> sampling can be also isotropic
00158       /*
00159       if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy); 
00160       if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy); 
00161       */
00162       //********************************************************************
00163       if(theIsoFlag){cosTh =2.*G4UniformRand()-1;}
00164       else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
00165       else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
00166       else{
00167         G4cout << "unknown distribution found for Angular"<<G4endl;
00168         throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
00169       }
00170       //********************************************************************
00171 //080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) 
00172 /*
00173     if(theAngularDistributionType == 1) // LAB
00174     {
00175       G4double en = aHadron.GetTotalMomentum();
00176       G4ReactionProduct boosted;
00177       boosted.Lorentz(theNeutron, theTarget);
00178       G4double kineticEnergy = boosted.GetKineticEnergy();
00179       G4double cosTh = theCoefficients->SampleMax(kineticEnergy);
00180       G4double theta = std::acos(cosTh);
00181       G4double phi = twopi*G4UniformRand();
00182       G4double sinth = std::sin(theta);
00183       G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00184       aHadron.SetMomentum( temp );
00185     }
00186     else if(theAngularDistributionType == 2) // costh in CMS {
00187     }
00188 */
00189 
00190 //      G4ReactionProduct boostedN;
00191 //      boostedN.Lorentz(theNeutron, theTarget);
00192 //      G4double kineticEnergy = boostedN.GetKineticEnergy();
00193 //      G4double cosTh = theProbArray->Sample(kineticEnergy); 
00194 
00195       G4double theta = std::acos(cosTh);
00196       G4double phi = twopi*G4UniformRand();
00197       G4double sinth = std::sin(theta);      
00198       
00199       G4ThreeVector temp(sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) ); //CMS
00200 
00201 //080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
00202 /*
00203       G4double en = aHadron.GetTotalEnergy(); // Target rest
00204       
00205       // get trafo from Target rest frame to CMS
00206       G4ReactionProduct boostedT;
00207       boostedT.Lorentz(theTarget, theTarget);
00208       
00209       G4ThreeVector the3Neutron = boostedN.GetMomentum();
00210       G4double nEnergy = boostedN.GetTotalEnergy();
00211       G4ThreeVector the3Target = boostedT.GetMomentum();
00212       G4double tEnergy = boostedT.GetTotalEnergy();
00213       G4double totE = nEnergy+tEnergy;
00214       G4ThreeVector the3trafo = -the3Target-the3Neutron;
00215       G4ReactionProduct trafo; // for transformation from CMS to target rest frame
00216       trafo.SetMomentum(the3trafo);
00217       G4double cmsMom = std::sqrt(the3trafo*the3trafo);
00218       G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
00219       trafo.SetMass(sqrts);
00220       trafo.SetTotalEnergy(totE);
00221       
00222       G4double gamma = trafo.GetTotalEnergy()/trafo.GetMass();
00223       G4double cosalpha = temp*trafo.GetMomentum()/trafo.GetTotalMomentum()/temp.mag();
00224       G4double fac = cosalpha*trafo.GetTotalMomentum()/trafo.GetMass();
00225       fac*=gamma;
00226       
00227       G4double mom;
00228 //    For G4FPE_DEBUG ON
00229       G4double mom2 = ( en*fac*en*fac - 
00230                    (fac*fac - gamma*gamma)*
00231                    (en*en - gamma*gamma*aHadron.GetMass()*aHadron.GetMass())
00232                 );
00233       if ( mom2 > 0.0 ) 
00234         mom = std::sqrt( mom2 );
00235       else
00236         mom = 0.0; 
00237 
00238       mom = -en*fac - mom;
00239       mom /= (fac*fac-gamma*gamma);
00240       temp = mom*temp;
00241       
00242       aHadron.SetMomentum( temp ); // now all in CMS
00243       aHadron.SetTotalEnergy( std::sqrt( mom*mom + aHadron.GetMass()*aHadron.GetMass() ) );
00244       aHadron.Lorentz(aHadron, trafo); // now in target rest frame
00245 */
00246       // Determination of the hadron kinetic energy in CMS
00247       // aHadron.GetKineticEnergy() is actually the residual kinetic energy in CMS (or target frame)
00248       // kineticEnergy is incident neutron kinetic energy  in CMS (or target frame)  
00249       G4double QValue = aHadron.GetKineticEnergy() - kineticEnergy;
00250       G4double A1     =   theTarget.GetMass()/boostedN.GetMass(); 
00251       G4double A1prim =   aHadron.GetMass()/ boostedN.GetMass();
00252       G4double kinE   = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*kineticEnergy+(1+A1)*QValue);
00253       G4double totalE = kinE + aHadron.GetMass();
00254       G4double mom2   = totalE*totalE - aHadron.GetMass()*aHadron.GetMass();
00255       G4double mom;
00256       if ( mom2 > 0.0 ) mom = std::sqrt( mom2 );
00257       else mom = 0.0;     
00258 
00259       aHadron.SetMomentum( mom*temp ); // Set momentum in CMS
00260       aHadron.SetKineticEnergy(kinE);  // Set kinetic energy in CMS
00261 
00262       // get trafo from Target rest frame to CMS
00263       G4ReactionProduct boostedT;
00264       boostedT.Lorentz(theTarget, theTarget);
00265       
00266       G4ThreeVector the3Neutron = boostedN.GetMomentum();
00267       G4double nEnergy = boostedN.GetTotalEnergy();
00268       G4ThreeVector the3Target = boostedT.GetMomentum();
00269       G4double tEnergy = boostedT.GetTotalEnergy();
00270       G4double totE = nEnergy+tEnergy;
00271       G4ThreeVector the3trafo = -the3Target-the3Neutron;
00272       G4ReactionProduct trafo; // for transformation from CMS to target rest frame
00273       trafo.SetMomentum(the3trafo);
00274       G4double cmsMom = std::sqrt(the3trafo*the3trafo);
00275       G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
00276       trafo.SetMass(sqrts);
00277       trafo.SetTotalEnergy(totE);
00278 
00279       aHadron.Lorentz(aHadron, trafo);
00280 
00281     }
00282     else
00283     {
00284       throw G4HadronicException(__FILE__, __LINE__, "Tried to sample non isotropic neutron angular");
00285     }
00286   aHadron.Lorentz(aHadron, -1.*theTarget); 
00287 //  G4cout << aHadron.GetMomentum()<<" ";
00288 //  G4cout << aHadron.GetTotalMomentum()<<G4endl;
00289 }

void G4NeutronHPAngular::SetNeutron ( const G4ReactionProduct aNeutron  )  [inline]

Definition at line 70 of file G4NeutronHPAngular.hh.

Referenced by G4FissionLibrary::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), G4NeutronHPInelasticCompFS::InitDistributionInitialState(), G4NeutronHPFSFissionFS::SetNeutron(), and G4NeutronHPFissionBaseFS::SetNeutron().

00070 { theNeutron = aNeutron; }

void G4NeutronHPAngular::SetTarget ( const G4ReactionProduct aTarget  )  [inline]

Definition at line 68 of file G4NeutronHPAngular.hh.

Referenced by G4FissionLibrary::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), G4NeutronHPInelasticCompFS::InitDistributionInitialState(), G4NeutronHPFSFissionFS::SetTarget(), and G4NeutronHPFissionBaseFS::SetTarget().

00068 { theTarget = aTarget; }


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