00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "G4NeutronHPLabAngularEnergy.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "Randomize.hh"
00036 #include "G4Gamma.hh"
00037 #include "G4Electron.hh"
00038 #include "G4Positron.hh"
00039 #include "G4Neutron.hh"
00040 #include "G4Proton.hh"
00041 #include "G4Deuteron.hh"
00042 #include "G4Triton.hh"
00043 #include "G4He3.hh"
00044 #include "G4Alpha.hh"
00045
00046 void G4NeutronHPLabAngularEnergy::Init(std::ifstream & aDataFile)
00047 {
00048 aDataFile >> nEnergies;
00049 theManager.Init(aDataFile);
00050 theEnergies = new G4double[nEnergies];
00051 nCosTh = new G4int[nEnergies];
00052 theData = new G4NeutronHPVector * [nEnergies];
00053 theSecondManager = new G4InterpolationManager [nEnergies];
00054 for(G4int i=0; i<nEnergies; i++)
00055 {
00056 aDataFile >> theEnergies[i];
00057 theEnergies[i]*=eV;
00058 aDataFile >> nCosTh[i];
00059 theSecondManager[i].Init(aDataFile);
00060 theData[i] = new G4NeutronHPVector[nCosTh[i]];
00061 G4double label;
00062 for(G4int ii=0; ii<nCosTh[i]; ii++)
00063 {
00064 aDataFile >> label;
00065 theData[i][ii].SetLabel(label);
00066 theData[i][ii].Init(aDataFile, eV);
00067 }
00068 }
00069 }
00070
00071 G4ReactionProduct * G4NeutronHPLabAngularEnergy::Sample(G4double anEnergy, G4double massCode, G4double )
00072 {
00073 G4ReactionProduct * result = new G4ReactionProduct;
00074 G4int Z = static_cast<G4int>(massCode/1000);
00075 G4int A = static_cast<G4int>(massCode-1000*Z);
00076
00077 if(massCode==0)
00078 {
00079 result->SetDefinition(G4Gamma::Gamma());
00080 }
00081 else if(A==0)
00082 {
00083 result->SetDefinition(G4Electron::Electron());
00084 if(Z==1) result->SetDefinition(G4Positron::Positron());
00085 }
00086 else if(A==1)
00087 {
00088 result->SetDefinition(G4Neutron::Neutron());
00089 if(Z==1) result->SetDefinition(G4Proton::Proton());
00090 }
00091 else if(A==2)
00092 {
00093 result->SetDefinition(G4Deuteron::Deuteron());
00094 }
00095 else if(A==3)
00096 {
00097 result->SetDefinition(G4Triton::Triton());
00098 if(Z==2) result->SetDefinition(G4He3::He3());
00099 }
00100 else if(A==4)
00101 {
00102 result->SetDefinition(G4Alpha::Alpha());
00103 if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
00104 }
00105 else
00106 {
00107 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPLabAngularEnergy: Unknown ion case 2");
00108 }
00109
00110
00111 G4double cosTh, secEnergy;
00112 G4int i, it(0);
00113
00114 for(i=0; i<nEnergies; i++)
00115 {
00116 it = i;
00117 if ( anEnergy < theEnergies[i] ) break;
00118 }
00119
00120
00121 if ( it == 0 )
00122 {
00123 if(it==0) G4cout << "080808 Something unexpected is happen in G4NeutronHPLabAngularEnergy " << G4endl;
00124
00125 G4double * running = new G4double [nCosTh[it]];
00126 running[0]=0;
00127 for(i=0;i<nCosTh[it]; i++)
00128 {
00129 if(i!=0) running[i] = running[i-1];
00130 running[i]+=theData[it][i].GetIntegral();
00131 }
00132 G4double random = running[nCosTh[it]-1]*G4UniformRand();
00133 G4int ith(0);
00134 for(i=0;i<nCosTh[it]; i++)
00135 {
00136 ith = i;
00137 if(random<running[i]) break;
00138 }
00139
00140
00141 if ( ith == 0 )
00142 {
00143 cosTh = theData[it][ith].GetLabel();
00144 secEnergy = theData[it][ith].Sample();
00145 currentMeanEnergy = theData[it][ith].GetMeanX();
00146 }
00147 else
00148 {
00149
00150
00151
00152 G4double x1 = running [ ith-1 ];
00153 G4double x2 = running [ ith ];
00154 G4double x = random;
00155 G4double y1 = theData[it][ith-1].GetLabel();
00156 G4double y2 = theData[it][ith].GetLabel();
00157 cosTh = theInt.Interpolate(theSecondManager[it].GetInverseScheme(ith),
00158 x, x1, x2, y1, y2);
00159 G4NeutronHPVector theBuff1;
00160 theBuff1.SetInterpolationManager(theData[it][ith-1].GetInterpolationManager());
00161 G4NeutronHPVector theBuff2;
00162 theBuff2.SetInterpolationManager(theData[it][ith].GetInterpolationManager());
00163 x1=y1;
00164 x2=y2;
00165 G4double y, mu;
00166 for(i=0;i<theData[it][ith-1].GetVectorLength(); i++)
00167 {
00168 mu = theData[it][ith-1].GetX(i);
00169 y1 = theData[it][ith-1].GetY(i);
00170 y2 = theData[it][ith].GetY(mu);
00171
00172 y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
00173 cosTh, x1,x2,y1,y2);
00174 theBuff1.SetData(i, mu, y);
00175 }
00176 for(i=0;i<theData[it][ith].GetVectorLength(); i++)
00177 {
00178 mu = theData[it][ith].GetX(i);
00179 y1 = theData[it][ith-1].GetY(mu);
00180 y2 = theData[it][ith].GetY(i);
00181 y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
00182 cosTh, x1,x2,y1,y2);
00183 theBuff2.SetData(i, mu, y);
00184 }
00185 G4NeutronHPVector theStore;
00186 theStore.Merge(&theBuff1, &theBuff2);
00187 secEnergy = theStore.Sample();
00188 currentMeanEnergy = theStore.GetMeanX();
00189 }
00190 delete [] running;
00191 }
00192 else
00193 {
00194 G4double x, x1, x2, y1, y2, y, tmp, E;
00195
00196 G4NeutronHPVector run1;
00197 run1.SetY(0, 0.);
00198 for(i=0;i<nCosTh[it-1]; i++)
00199 {
00200 if(i!=0) run1.SetY(i, run1.GetY(i-1));
00201 run1.SetX(i, theData[it-1][i].GetLabel());
00202 run1.SetY(i, run1.GetY(i)+theData[it-1][i].GetIntegral());
00203 }
00204 G4NeutronHPVector run2;
00205 run2.SetY(0, 0.);
00206 for(i=0;i<nCosTh[it]; i++)
00207 {
00208 if(i!=0) run2.SetY(i, run2.GetY(i-1));
00209 run2.SetX(i, theData[it][i].GetLabel());
00210 run2.SetY(i, run2.GetY(i)+theData[it][i].GetIntegral());
00211 }
00212
00213 x = anEnergy;
00214 x1 = theEnergies[it-1];
00215 x2 = theEnergies[it];
00216 G4NeutronHPVector thBuff1;
00217 thBuff1.SetInterpolationManager(theSecondManager[it-1]);
00218 for(i=0; i<run1.GetVectorLength(); i++)
00219 {
00220 tmp = run1.GetX(i);
00221 y1 = run1.GetY(i);
00222 y2 = run2.GetY(tmp);
00223 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00224 thBuff1.SetData(i, tmp, y);
00225 }
00226 G4NeutronHPVector thBuff2;
00227 thBuff2.SetInterpolationManager(theSecondManager[it]);
00228 for(i=0; i<run2.GetVectorLength(); i++)
00229 {
00230 tmp = run2.GetX(i);
00231 y1 = run1.GetY(tmp);
00232 y2 = run2.GetY(i);
00233 y = theInt.Lin(x, x1,x2,y1,y2);
00234 thBuff2.SetData(i, tmp, y);
00235 }
00236 G4NeutronHPVector theThVec;
00237 theThVec.Merge(&thBuff1 ,&thBuff2);
00238 G4double random = (theThVec.GetY(theThVec.GetVectorLength()-1)
00239 -theThVec.GetY(0)) *G4UniformRand();
00240 G4int ith(0);
00241 for(i=1;i<theThVec.GetVectorLength(); i++)
00242 {
00243 ith = i;
00244 if(random<theThVec.GetY(i)-theThVec.GetY(0)) break;
00245 }
00246 {
00247
00248 G4double xx, xx1, xx2, yy1, yy2;
00249 xx = random;
00250 xx1 = theThVec.GetY(ith-1)-theThVec.GetY(0);
00251 xx2 = theThVec.GetY(ith)-theThVec.GetY(0);
00252 yy1 = theThVec.GetX(ith-1);
00253 yy2 = theThVec.GetX(ith);
00254 cosTh = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
00255 xx, xx1,xx2,yy1,yy2);
00256 }
00257 G4int i1(0), i2(0);
00258
00259
00260 for(i=0; i<nCosTh[it-1]; i++)
00261 {
00262 i1 = i;
00263 if(cosTh<theData[it-1][i].GetLabel()) break;
00264 }
00265
00266 x = cosTh;
00267 x1 = theData[it-1][i1-1].GetLabel();
00268 x2 = theData[it-1][i1].GetLabel();
00269 G4NeutronHPVector theBuff1a;
00270 theBuff1a.SetInterpolationManager(theData[it-1][i1-1].GetInterpolationManager());
00271 for(i=0;i<theData[it-1][i1-1].GetVectorLength(); i++)
00272 {
00273 E = theData[it-1][i1-1].GetX(i);
00274 y1 = theData[it-1][i1-1].GetY(i);
00275 y2 = theData[it-1][i1].GetY(E);
00276 y = theInt.Lin(x, x1,x2,y1,y2);
00277 theBuff1a.SetData(i, E, y);
00278 }
00279 G4NeutronHPVector theBuff2a;
00280 theBuff2a.SetInterpolationManager(theData[it-1][i1].GetInterpolationManager());
00281 for(i=0;i<theData[it-1][i1].GetVectorLength(); i++)
00282 {
00283 E = theData[it-1][i1].GetX(i);
00284 y1 = theData[it-1][i1-1].GetY(E);
00285 y2 = theData[it-1][i1].GetY(i);
00286 y = theInt.Lin(x, x1,x2,y1,y2);
00287 theBuff2a.SetData(i, E, y);
00288 }
00289 G4NeutronHPVector theStore1;
00290 theStore1.Merge(&theBuff1a, &theBuff2a);
00291
00292
00293
00294 for(i=0; i<nCosTh[it]; i++)
00295 {
00296 i2 = i;
00297 if(cosTh<theData[it][i2].GetLabel()) break;
00298 }
00299 x1 = theData[it][i2-1].GetLabel();
00300 x2 = theData[it][i2].GetLabel();
00301 G4NeutronHPVector theBuff1b;
00302 theBuff1b.SetInterpolationManager(theData[it][i2-1].GetInterpolationManager());
00303 for(i=0;i<theData[it][i2-1].GetVectorLength(); i++)
00304 {
00305 E = theData[it][i2-1].GetX(i);
00306 y1 = theData[it][i2-1].GetY(i);
00307 y2 = theData[it][i2].GetY(E);
00308 y = theInt.Lin(x, x1,x2,y1,y2);
00309 theBuff1b.SetData(i, E, y);
00310 }
00311 G4NeutronHPVector theBuff2b;
00312 theBuff2b.SetInterpolationManager(theData[it][i2].GetInterpolationManager());
00313
00314
00315 for(i=0;i<theData[it][i2].GetVectorLength(); i++)
00316 {
00317
00318
00319
00320 E = theData[it][i2].GetX(i);
00321 y1 = theData[it][i2-1].GetY(E);
00322 y2 = theData[it][i2].GetY(i);
00323 y = theInt.Lin(x, x1,x2,y1,y2);
00324 theBuff2b.SetData(i, E, y);
00325 }
00326 G4NeutronHPVector theStore2;
00327 theStore2.Merge(&theBuff1b, &theBuff2b);
00328
00329
00330 x = anEnergy;
00331 x1 = theEnergies[it-1];
00332 x2 = theEnergies[it];
00333 G4NeutronHPVector theOne1;
00334 theOne1.SetInterpolationManager(theStore1.GetInterpolationManager());
00335 for(i=0; i<theStore1.GetVectorLength(); i++)
00336 {
00337 E = theStore1.GetX(i);
00338 y1 = theStore1.GetY(i);
00339 y2 = theStore2.GetY(E);
00340 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00341 theOne1.SetData(i, E, y);
00342 }
00343 G4NeutronHPVector theOne2;
00344 theOne2.SetInterpolationManager(theStore2.GetInterpolationManager());
00345 for(i=0; i<theStore2.GetVectorLength(); i++)
00346 {
00347 E = theStore2.GetX(i);
00348 y1 = theStore1.GetY(E);
00349 y2 = theStore2.GetY(i);
00350 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00351 theOne2.SetData(i, E, y);
00352 }
00353 G4NeutronHPVector theOne;
00354 theOne.Merge(&theOne1, &theOne2);
00355
00356 secEnergy = theOne.Sample();
00357 currentMeanEnergy = theOne.GetMeanX();
00358 }
00359
00360
00361
00362 result->SetKineticEnergy(secEnergy);
00363
00364 G4double phi = twopi*G4UniformRand();
00365 G4double theta = std::acos(cosTh);
00366 G4double sinth = std::sin(theta);
00367 G4double mtot = result->GetTotalMomentum();
00368 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
00369 result->SetMomentum(tempVector);
00370
00371 return result;
00372 }