#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 () |
Definition at line 42 of file G4NeutronHPAngular.hh.
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] |
G4double G4NeutronHPAngular::GetTargetMass | ( | ) | [inline] |
Definition at line 72 of file G4NeutronHPAngular.hh.
Referenced by G4NeutronHPInelasticBaseFS::BaseApply().
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().
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().