G4NeutronHPElasticFS Class Reference

#include <G4NeutronHPElasticFS.hh>

Inheritance diagram for G4NeutronHPElasticFS:

G4NeutronHPFinalState

Public Member Functions

 G4NeutronHPElasticFS ()
 ~G4NeutronHPElasticFS ()
void Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType)
G4HadFinalState * ApplyYourself (const G4HadProjectile &theTrack)
G4NeutronHPFinalState * New ()

Detailed Description

Definition at line 45 of file G4NeutronHPElasticFS.hh.

Constructor & Destructor Documentation

G4NeutronHPElasticFS::G4NeutronHPElasticFS (  )  [inline]

Definition at line 49 of file G4NeutronHPElasticFS.hh.

References G4NeutronHPFinalState::hasXsec.

00050   {
00051     hasXsec = false; 
00052     theCoefficients = 0;
00053     theProbArray = 0;
00054   }

G4NeutronHPElasticFS::~G4NeutronHPElasticFS (  )  [inline]

Definition at line 55 of file G4NeutronHPElasticFS.hh.

00056   {
00057     if(theCoefficients!=0) delete theCoefficients;
00058     if(theProbArray!=0) delete theProbArray;
00059   }

Member Function Documentation

G4HadFinalState * G4NeutronHPElasticFS::ApplyYourself ( const G4HadProjectile &  theTrack  )  [virtual]

Reimplemented from G4NeutronHPFinalState.

Definition at line 182 of file G4NeutronHPElasticFS.cc.

References G4HadFinalState::AddSecondary(), G4Alpha::Alpha(), G4InuclParticleNames::ap, G4HadFinalState::Clear(), G4Deuteron::Deuteron(), G4cout, G4endl, G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleus::GetBiasedThermalNucleus(), G4HadProjectile::GetDefinition(), G4ReactionProduct::GetKineticEnergy(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4HadProjectile::GetMaterial(), G4ReactionProduct::GetMomentum(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4Material::GetTemperature(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4He3::He3(), G4ReactionProduct::Lorentz(), G4Proton::Proton(), G4NeutronHPPartial::Sample(), G4NeutronHPLegendreStore::SampleElastic(), G4DynamicParticle::SetDefinition(), G4HadFinalState::SetEnergyChange(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4DynamicParticle::SetMomentum(), G4ReactionProduct::SetMomentum(), G4HadFinalState::SetMomentumChange(), G4HadFinalState::SetStatusChange(), G4ReactionProduct::SetTotalEnergy(), suspend, G4NeutronHPFinalState::theBaseA, G4NeutronHPFinalState::theBaseZ, G4NeutronHPFinalState::theResult, and G4Triton::Triton().

00183   {  
00184 //    G4cout << "G4NeutronHPElasticFS::ApplyYourself+"<<G4endl;
00185     theResult.Clear();
00186     G4double eKinetic = theTrack.GetKineticEnergy();
00187     const G4HadProjectile *incidentParticle = &theTrack;
00188     G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
00189     theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
00190     theNeutron.SetKineticEnergy( eKinetic );
00191 //    G4cout << "G4NeutronHPElasticFS::ApplyYourself++"<<eKinetic<<" "<<G4endl;
00192 //    G4cout << "CMSVALUES 0 "<<theNeutron.GetTotalMomentum()<<G4endl;
00193     
00194     G4ReactionProduct theTarget; 
00195     G4Nucleus aNucleus;
00196     G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
00197     theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
00198 //     G4cout << "Nucleus-test"<<" "<<targetMass<<" ";
00199 //     G4cout << theTarget.GetMomentum().x()<<" ";
00200 //     G4cout << theTarget.GetMomentum().y()<<" ";
00201 //     G4cout << theTarget.GetMomentum().z()<<G4endl;
00202     
00203     // neutron and target defined as reaction products.
00204 
00205 // prepare lorentz-transformation to Lab.
00206 
00207     G4ThreeVector the3Neutron = theNeutron.GetMomentum();
00208     G4double nEnergy = theNeutron.GetTotalEnergy();
00209     G4ThreeVector the3Target = theTarget.GetMomentum();
00210 //    cout << "@@@" << the3Target<<G4endl;
00211     G4double tEnergy = theTarget.GetTotalEnergy();
00212     G4ReactionProduct theCMS;
00213     G4double totE = nEnergy+tEnergy;
00214     G4ThreeVector the3CMS = the3Target+the3Neutron;
00215     theCMS.SetMomentum(the3CMS);
00216     G4double cmsMom = std::sqrt(the3CMS*the3CMS);
00217     G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
00218     theCMS.SetMass(sqrts);
00219     theCMS.SetTotalEnergy(totE);
00220     
00221     // data come as fcn of n-energy in nuclear rest frame
00222     G4ReactionProduct boosted;
00223     boosted.Lorentz(theNeutron, theTarget);
00224     eKinetic = boosted.GetKineticEnergy(); // get kinetic energy for scattering
00225     G4double cosTh = -2;
00226     if(repFlag == 1)
00227     {
00228       cosTh = theCoefficients->SampleElastic(eKinetic);
00229     }
00230     
00231     else if (repFlag==2)
00232     {
00233       cosTh = theProbArray->Sample(eKinetic);
00234     }
00235     else if (repFlag==3)
00236     {
00237        if ( eKinetic <= tE_of_repFlag3 )
00238        {
00239           cosTh = theCoefficients->SampleElastic(eKinetic);
00240        }
00241        else
00242        {
00243           cosTh = theProbArray->Sample(eKinetic);
00244        }
00245     }
00246     else if (repFlag==0)
00247     {
00248       cosTh = 2.*G4UniformRand()-1.;
00249     }
00250     else
00251     {
00252       G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
00253       throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::Init -- unusable number for repFlag");
00254     }
00255     if(cosTh<-1.1) { return 0; }
00256     G4double phi = twopi*G4UniformRand();
00257     G4double theta = std::acos(cosTh);
00258     G4double sinth = std::sin(theta);
00259     if (frameFlag == 1) // final state data given in target rest frame.
00260     {
00261       // we have the scattering angle, now we need the energy, then do the
00262       // boosting.
00263       // relativistic elastic scattering energy angular correlation:
00264       theNeutron.Lorentz(theNeutron, theTarget);
00265       G4double e0 = theNeutron.GetTotalEnergy();
00266       G4double p0 = theNeutron.GetTotalMomentum();
00267       G4double mN = theNeutron.GetMass();
00268       G4double mT = theTarget.GetMass();
00269       G4double eE = e0+mT;
00270       G4double ap = (mT+eE)*(mT-eE) + (p0+mN)*(p0-mN);
00271       G4double a = 4*(eE+p0*cosTh)*(eE-p0*cosTh);
00272       G4double b = 4*ap*p0*cosTh;
00273       G4double c = (2.*eE*mN-ap)*(2.*eE*mN+ap);
00274       G4double en = (-b+std::sqrt(b*b - 4*a*c) )/(2*a);
00275       G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00276       theNeutron.SetMomentum(tempVector);
00277       theNeutron.SetTotalEnergy(std::sqrt(en*en+theNeutron.GetMass()*theNeutron.GetMass()));
00278       // first to lab     
00279       theNeutron.Lorentz(theNeutron, -1.*theTarget);
00280       // now to CMS
00281       theNeutron.Lorentz(theNeutron, theCMS);
00282       theTarget.SetMomentum(-theNeutron.GetMomentum());
00283       theTarget.SetTotalEnergy(theNeutron.GetTotalEnergy());
00284       // and back to lab
00285       theNeutron.Lorentz(theNeutron, -1.*theCMS);
00286       theTarget.Lorentz(theTarget, -1.*theCMS);      
00287 //111005 Protection for not producing 0 kinetic energy target
00288       if ( theNeutron.GetKineticEnergy() <= 0 ) theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00289       if ( theTarget.GetKineticEnergy() <= 0 ) theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00290     }
00291     else if (frameFlag == 2) // CMS
00292     {
00293       theNeutron.Lorentz(theNeutron, theCMS);
00294       theTarget.Lorentz(theTarget, theCMS);
00295       G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
00296       G4ThreeVector cmsMom_tmp=theNeutron.GetMomentum(); // for neutron direction in CMS
00297       G4double cms_theta=cmsMom_tmp.theta();
00298       G4double cms_phi=cmsMom_tmp.phi();
00299       G4ThreeVector tempVector;
00300       tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
00301                       +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
00302                       -std::sin(theta)*std::sin(phi)*std::sin(cms_phi)  );
00303       tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
00304                       +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
00305                       +std::sin(theta)*std::sin(phi)*std::cos(cms_phi)  );
00306       tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
00307                       -std::sin(theta)*std::cos(phi)*std::sin(cms_theta)  );
00308       tempVector *= en;
00309       theNeutron.SetMomentum(tempVector);
00310       theTarget.SetMomentum(-tempVector);
00311       G4double tP = theTarget.GetTotalMomentum();
00312       G4double tM = theTarget.GetMass();
00313       theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
00314 
00315 /*
00316 For debug purpose. 
00317 Same transformation G4ReactionProduct.Lorentz() by 4vectors
00318 {
00319 G4LorentzVector n4p = G4LorentzVector ( theNeutron.GetMomentum() , theNeutron.GetKineticEnergy() + theNeutron.GetMass() );    
00320 G4cout << "before " << ( n4p.e() - n4p.m() ) / eV<< G4endl;
00321 G4LorentzVector cm4p = G4LorentzVector ( theCMS.GetMomentum() , theCMS.GetKineticEnergy() + theCMS.GetMass() );    
00322 n4p.boost( cm4p.boostVector() );
00323 G4cout << cm4p/eV << G4endl;
00324 G4cout << "after " <<  ( n4p.e() - n4p.m() ) / eV<< G4endl;
00325 }
00326 */
00327 
00328       theNeutron.Lorentz(theNeutron, -1.*theCMS);
00329 //080904 Add Protection for very low energy (1e-6eV) scattering 
00330       if ( theNeutron.GetKineticEnergy() <= 0 )
00331       {
00332          //theNeutron.SetMomentum( G4ThreeVector(0) ); 
00333          //theNeutron.SetTotalEnergy ( theNeutron.GetMass() );
00334 //110822 Protection for not producing 0 kinetic energy neutron
00335          theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00336       }
00337 
00338       theTarget.Lorentz(theTarget, -1.*theCMS);
00339 //080904 Add Protection for very low energy (1e-6eV) scattering 
00340       if ( theTarget.GetKineticEnergy() < 0 )
00341       {
00342          //theTarget.SetMomentum( G4ThreeVector(0) ); 
00343          //theTarget.SetTotalEnergy ( theTarget.GetMass()  );
00344 //110822 Protection for not producing 0 kinetic energy target
00345          theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00346       }
00347     }
00348     else
00349     {
00350       G4cout <<"Value of frameFlag (1=LAB, 2=CMS): "<<frameFlag;
00351       throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::ApplyYourSelf frameflag incorrect");
00352     }
00353     // now all in Lab
00354 // nun den recoil generieren...und energy change, momentum change angeben.
00355     theResult.SetEnergyChange(theNeutron.GetKineticEnergy());
00356     theResult.SetMomentumChange(theNeutron.GetMomentum().unit());
00357     G4DynamicParticle* theRecoil = new G4DynamicParticle;
00358     if(targetMass<4.5)
00359     {
00360       if(targetMass<1)
00361       {
00362         // proton
00363         theRecoil->SetDefinition(G4Proton::Proton());
00364       }
00365       else if(targetMass<2 )
00366       {
00367         // deuteron
00368         theRecoil->SetDefinition(G4Deuteron::Deuteron());
00369       }
00370       else if(targetMass<2.999 )
00371       {
00372         // 3He 
00373         theRecoil->SetDefinition(G4He3::He3());
00374       }
00375       else if(targetMass<3 )
00376       {
00377         // Triton
00378         theRecoil->SetDefinition(G4Triton::Triton());
00379       }
00380       else
00381       {
00382         // alpha
00383         theRecoil->SetDefinition(G4Alpha::Alpha());
00384       }
00385     }
00386     else
00387     {
00388       theRecoil->SetDefinition(G4ParticleTable::GetParticleTable()
00389                                ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ)));
00390     }
00391     theRecoil->SetMomentum(theTarget.GetMomentum());
00392     theResult.AddSecondary(theRecoil);
00393 //    G4cout << "G4NeutronHPElasticFS::ApplyYourself 10+"<<G4endl;
00394     // postpone the tracking of the primary neutron
00395      theResult.SetStatusChange(suspend);
00396     return &theResult;
00397   }

void G4NeutronHPElasticFS::Init ( G4double  A,
G4double  Z,
G4int  M,
G4String &  dirName,
G4String &  aFSType 
) [virtual]

Implements G4NeutronHPFinalState.

Definition at line 48 of file G4NeutronHPElasticFS.cc.

References G4cout, G4endl, G4NeutronHPDataUsed::GetName(), G4NeutronHPNames::GetName(), G4NeutronHPFinalState::hasAnyData, G4NeutronHPFinalState::hasFSData, G4NeutronHPFinalState::hasXsec, G4NeutronHPLegendreStore::Init(), G4NeutronHPPartial::InitInterpolation(), G4NeutronHPLegendreStore::InitInterpolation(), G4NeutronHPFinalState::SetAZMs(), G4NeutronHPLegendreStore::SetCoeff(), G4NeutronHPPartial::SetT(), G4NeutronHPLegendreStore::SetTemperature(), G4NeutronHPPartial::SetX(), G4NeutronHPPartial::SetY(), and G4NeutronHPFinalState::theNames.

00049   {
00050     G4String tString = "/FS";
00051     G4bool dbool;
00052     G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
00053     G4String filename = aFile.GetName();
00054     SetAZMs( A, Z, M, aFile ); 
00055     //theBaseA = aFile.GetA();
00056     //theBaseZ = aFile.GetZ();
00057     if(!dbool)
00058     {
00059       hasAnyData = false;
00060       hasFSData = false; 
00061       hasXsec = false;
00062       return;
00063     }
00064     std::ifstream theData(filename, std::ios::in);
00065     theData >> repFlag >> targetMass >> frameFlag;
00066     if(repFlag==1)
00067     {
00068       G4int nEnergy;
00069       theData >> nEnergy; 
00070       theCoefficients = new G4NeutronHPLegendreStore(nEnergy);
00071       theCoefficients->InitInterpolation(theData);
00072       G4double temp, energy;
00073       G4int tempdep, nLegendre;
00074       G4int i, ii;
00075       for (i=0; i<nEnergy; i++)
00076       {
00077         theData >> temp >> energy >> tempdep >> nLegendre;
00078         energy *=eV;
00079         theCoefficients->Init(i, energy, nLegendre);
00080         theCoefficients->SetTemperature(i, temp);
00081         G4double coeff=0;
00082         for(ii=0; ii<nLegendre; ii++)
00083         {
00084           // load legendre coefficients.
00085           theData >> coeff;
00086           theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
00087         }
00088       }
00089     }
00090     else if (repFlag==2)
00091     {
00092       G4int nEnergy;
00093       theData >> nEnergy;
00094       theProbArray = new G4NeutronHPPartial(nEnergy, nEnergy);
00095       theProbArray->InitInterpolation(theData);
00096       G4double temp, energy;
00097       G4int tempdep, nPoints;
00098       for(G4int i=0; i<nEnergy; i++)
00099       {
00100         theData >> temp >> energy >> tempdep >> nPoints;
00101         energy *= eV;
00102         theProbArray->InitInterpolation(i, theData);
00103         theProbArray->SetT(i, temp);
00104         theProbArray->SetX(i, energy);
00105         G4double prob, costh;
00106         for(G4int ii=0; ii<nPoints; ii++)
00107         {
00108           // fill probability arrays.
00109           theData >> costh >> prob;
00110           theProbArray->SetX(i, ii, costh);
00111           theProbArray->SetY(i, ii, prob);
00112         }
00113       }
00114     }
00115     else if ( repFlag==3 )
00116     {
00117        G4int nEnergy_Legendre;
00118        theData >> nEnergy_Legendre; 
00119        theCoefficients = new G4NeutronHPLegendreStore( nEnergy_Legendre );
00120        theCoefficients->InitInterpolation( theData );
00121        G4double temp, energy;
00122        G4int tempdep, nLegendre;
00123        //G4int i, ii;
00124        for ( G4int i = 0 ; i < nEnergy_Legendre ; i++ )
00125        {
00126           theData >> temp >> energy >> tempdep >> nLegendre;
00127           energy *=eV;
00128           theCoefficients->Init( i , energy , nLegendre );
00129           theCoefficients->SetTemperature( i , temp );
00130           G4double coeff = 0;
00131           for (G4int ii = 0 ; ii < nLegendre ; ii++ )
00132           {
00133              // load legendre coefficients.
00134              theData >> coeff;
00135              theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
00136           }
00137        } 
00138 
00139        tE_of_repFlag3 = energy; 
00140 
00141        G4int nEnergy_Prob;
00142        theData >> nEnergy_Prob;
00143        theProbArray = new G4NeutronHPPartial( nEnergy_Prob , nEnergy_Prob );
00144        theProbArray->InitInterpolation( theData );
00145        G4int nPoints;
00146        for ( G4int i=0 ; i < nEnergy_Prob ; i++ )
00147        {
00148           theData >> temp >> energy >> tempdep >> nPoints;
00149 
00150           energy *= eV;
00151 
00152 //        consistency check
00153           if ( i == 0 )
00154              //if ( energy != tE_of_repFlag3 ) //110620TK This is too tight for 32bit machines 
00155              if ( std::abs( energy - tE_of_repFlag3 ) / tE_of_repFlag3 > 1.0e-15 )
00156                 G4cout << "Warning Transition Energy of repFlag3 is not consistent." << G4endl; 
00157 
00158           theProbArray->InitInterpolation( i , theData );
00159           theProbArray->SetT( i , temp );
00160           theProbArray->SetX( i , energy );
00161           G4double prob, costh;
00162           for( G4int ii = 0 ; ii < nPoints ; ii++ )
00163           {
00164              // fill probability arrays.
00165              theData >> costh >> prob;
00166              theProbArray->SetX( i , ii , costh );
00167              theProbArray->SetY( i , ii , prob );
00168           }
00169        }
00170     }
00171     else if (repFlag==0)
00172     {
00173       theData >> frameFlag;
00174     }
00175     else
00176     {
00177       G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
00178       throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::Init -- unusable number for repFlag");
00179     }
00180     theData.close();
00181   }

G4NeutronHPFinalState* G4NeutronHPElasticFS::New (  )  [inline, virtual]

Implements G4NeutronHPFinalState.

Definition at line 62 of file G4NeutronHPElasticFS.hh.

00063   {
00064    G4NeutronHPElasticFS * theNew = new G4NeutronHPElasticFS;
00065    return theNew;
00066   }

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