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
00033
00034
00035
00036
00037
00038
00039
00040 #include "G4NeutronHPContAngularPar.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4NeutronHPLegendreStore.hh"
00044 #include "G4Gamma.hh"
00045 #include "G4Electron.hh"
00046 #include "G4Positron.hh"
00047 #include "G4Neutron.hh"
00048 #include "G4Proton.hh"
00049 #include "G4Deuteron.hh"
00050 #include "G4Triton.hh"
00051 #include "G4He3.hh"
00052 #include "G4Alpha.hh"
00053 #include "G4NeutronHPVector.hh"
00054 #include "G4NucleiProperties.hh"
00055 #include "G4NeutronHPKallbachMannSyst.hh"
00056 #include "G4ParticleTable.hh"
00057
00058 void G4NeutronHPContAngularPar::Init(std::ifstream & aDataFile)
00059 {
00060 aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
00061 theEnergy *= eV;
00062 theAngular = new G4NeutronHPList [nEnergies];
00063 for(G4int i=0; i<nEnergies; i++)
00064 {
00065 G4double sEnergy;
00066 aDataFile >> sEnergy;
00067 sEnergy*=eV;
00068 theAngular[i].SetLabel(sEnergy);
00069 theAngular[i].Init(aDataFile, nAngularParameters, 1.);
00070 }
00071 }
00072
00073 G4ReactionProduct *
00074 G4NeutronHPContAngularPar::Sample(G4double anEnergy, G4double massCode, G4double ,
00075 G4int angularRep, G4int )
00076 {
00077 G4ReactionProduct * result = new G4ReactionProduct;
00078 G4int Z = static_cast<G4int>(massCode/1000);
00079 G4int A = static_cast<G4int>(massCode-1000*Z);
00080 if(massCode==0)
00081 {
00082 result->SetDefinition(G4Gamma::Gamma());
00083 }
00084 else if(A==0)
00085 {
00086 result->SetDefinition(G4Electron::Electron());
00087 if(Z==1) result->SetDefinition(G4Positron::Positron());
00088 }
00089 else if(A==1)
00090 {
00091 result->SetDefinition(G4Neutron::Neutron());
00092 if(Z==1) result->SetDefinition(G4Proton::Proton());
00093 }
00094 else if(A==2)
00095 {
00096 result->SetDefinition(G4Deuteron::Deuteron());
00097 }
00098 else if(A==3)
00099 {
00100 result->SetDefinition(G4Triton::Triton());
00101 if(Z==2) result->SetDefinition(G4He3::He3());
00102 }
00103 else if(A==4)
00104 {
00105 result->SetDefinition(G4Alpha::Alpha());
00106 if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar: Unknown ion case 1");
00107 }
00108 else
00109 {
00110 result->SetDefinition(G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z));
00111 }
00112 G4int i(0);
00113 G4int it(0);
00114 G4double fsEnergy(0);
00115 G4double cosTh(0);
00116
00117 if( angularRep == 1 )
00118 {
00119
00120
00121
00122
00123 if ( nDiscreteEnergies != 0 )
00124 {
00125
00126
00127
00128 if ( fresh == true )
00129 {
00130
00131
00132 remaining_energy = std::max ( theAngular[0].GetLabel() , theAngular[nEnergies-1].GetLabel() );
00133 fresh = false;
00134 }
00135
00136
00137
00138 if ( nDiscreteEnergies == nEnergies )
00139 {
00140 remaining_energy = std::max ( remaining_energy , theAngular[nDiscreteEnergies-1].GetLabel() );
00141 }
00142 else
00143 {
00144
00145
00146 G4double cont_min=0.0;
00147 for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
00148 {
00149 cont_min = theAngular[j].GetLabel();
00150 if ( theAngular[j].GetValue(0) != 0.0 ) break;
00151 }
00152 remaining_energy = std::max ( remaining_energy , std::min ( theAngular[nDiscreteEnergies-1].GetLabel() , cont_min ) );
00153 }
00154
00155 G4double random = G4UniformRand();
00156
00157 G4double * running = new G4double[nEnergies+1];
00158 running[0] = 0.0;
00159
00160 for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
00161 {
00162 G4double delta = 0.0;
00163 if ( theAngular[j].GetLabel() <= remaining_energy ) delta = theAngular[i].GetValue(0);
00164 running[j+1] = running[j] + delta;
00165 }
00166 G4double tot_prob_DIS = running[ nDiscreteEnergies ];
00167
00168 for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
00169 {
00170 G4double delta = 0.0;
00171 G4double e_low = 0.0;
00172 G4double e_high = 0.0;
00173 if ( theAngular[j].GetLabel() <= remaining_energy ) delta = theAngular[j].GetValue(0);
00174
00175
00176
00177
00178
00179
00180 if ( theAngular[j].GetLabel() != 0 )
00181 {
00182 if ( j == nDiscreteEnergies ) {
00183 e_low = 0.0/eV;
00184 } else {
00185 e_low = theAngular[j-1].GetLabel()/eV;
00186 }
00187 e_high = theAngular[j].GetLabel()/eV;
00188 }
00189
00190
00191 if ( theAngular[j].GetLabel() == 0.0 ) {
00192 e_low = theAngular[j].GetLabel()/eV;
00193 if ( j != nEnergies-1 ) {
00194 e_high = theAngular[j+1].GetLabel()/eV;
00195 } else {
00196 e_high = theAngular[j].GetLabel()/eV;
00197 if ( theAngular[j].GetValue(0) != 0.0 ) {
00198 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar: Unexpected non zero value of theAngular[nEnergies-1].GetValue(0)");
00199 }
00200 }
00201 }
00202
00203 running[j+1] = running[j] + ( ( e_high - e_low ) * delta );
00204 }
00205 G4double tot_prob_CON = running[ nEnergies ] - running[ nDiscreteEnergies ];
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219 random *= (tot_prob_DIS + tot_prob_CON);
00220
00221 if ( random <= ( tot_prob_DIS / ( tot_prob_DIS + tot_prob_CON ) ) || nDiscreteEnergies == nEnergies )
00222 {
00223
00224 for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
00225 {
00226
00227 if ( random < running[ j+1 ] )
00228 {
00229 it = j;
00230 break;
00231 }
00232 }
00233 fsEnergy = theAngular[ it ].GetLabel();
00234
00235 G4NeutronHPLegendreStore theStore(1);
00236 theStore.Init(0,fsEnergy,nAngularParameters);
00237 for (G4int j=0;j<nAngularParameters;j++)
00238 {
00239 theStore.SetCoeff(0,j,theAngular[it].GetValue(j));
00240 }
00241
00242 cosTh = theStore.SampleMax(fsEnergy);
00243
00244 }
00245 else
00246 {
00247
00248 for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
00249 {
00250
00251 if ( random < running[ j ] )
00252 {
00253 it = j;
00254 break;
00255 }
00256 }
00257
00258 G4double x1 = running[it-1];
00259 G4double x2 = running[it];
00260
00261 G4double y1 = 0.0;
00262 if ( it != nDiscreteEnergies )
00263 y1 = theAngular[it-1].GetLabel();
00264 G4double y2 = theAngular[it].GetLabel();
00265
00266 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
00267 random,x1,x2,y1,y2);
00268
00269 G4NeutronHPLegendreStore theStore(2);
00270 theStore.Init(0,y1,nAngularParameters);
00271 theStore.Init(1,y2,nAngularParameters);
00272 theStore.SetManager(theManager);
00273 for (G4int j=0;j<nAngularParameters;j++)
00274 {
00275 G4int itt = it;
00276 if ( it == nDiscreteEnergies ) itt = it+1;
00277 if ( it == 0 )
00278 {
00279
00280
00281 itt = it+1;
00282 }
00283 theStore.SetCoeff(0,j,theAngular[itt-1].GetValue(j));
00284 theStore.SetCoeff(1,j,theAngular[itt].GetValue(j));
00285 }
00286
00287 cosTh = theStore.SampleMax(fsEnergy);
00288
00289
00290 }
00291
00292
00293 remaining_energy -= fsEnergy;
00294
00295
00296
00297 delete[] running;
00298
00299 }
00300 else
00301 {
00302
00303
00304
00305 if ( fresh == true )
00306 {
00307 remaining_energy = theAngular[ nEnergies-1 ].GetLabel();
00308 fresh = false;
00309 }
00310
00311 G4double random = G4UniformRand();
00312 G4double * running = new G4double[nEnergies];
00313 running[0]=0;
00314 G4double weighted = 0;
00315 for(i=1; i<nEnergies; i++)
00316 {
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330 running[i]=running[i-1];
00331 if ( remaining_energy >= theAngular[i].GetLabel() )
00332 {
00333 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
00334 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
00335 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
00336 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
00337 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
00338 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
00339 }
00340 }
00341
00342
00343 if ( nEnergies == 1 || running[nEnergies-1] == 0 )
00344 currentMeanEnergy = 0.0;
00345 else
00346 {
00347 currentMeanEnergy = weighted/running[nEnergies-1];
00348 }
00349
00350
00351 if ( nEnergies == 1 ) it = 0;
00352
00353
00354 if ( running[nEnergies-1] != 0 )
00355 {
00356 for ( i = 1 ; i < nEnergies ; i++ )
00357 {
00358 it = i;
00359 if ( random < running [ i ] / running [ nEnergies-1 ] ) break;
00360 }
00361 }
00362
00363
00364 if ( running [ nEnergies-1 ] == 0 ) it = 0;
00365
00366
00367 if (it<nDiscreteEnergies||it==0)
00368 {
00369 if(it == 0)
00370 {
00371 fsEnergy = theAngular[it].GetLabel();
00372 G4NeutronHPLegendreStore theStore(1);
00373 theStore.Init(0,fsEnergy,nAngularParameters);
00374 for(i=0;i<nAngularParameters;i++)
00375 {
00376 theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
00377 }
00378
00379 cosTh = theStore.SampleMax(fsEnergy);
00380 }
00381 else
00382 {
00383 G4double e1, e2;
00384 e1 = theAngular[it-1].GetLabel();
00385 e2 = theAngular[it].GetLabel();
00386 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
00387 random,
00388 running[it-1]/running[nEnergies-1],
00389 running[it]/running[nEnergies-1],
00390 e1, e2);
00391
00392 G4NeutronHPLegendreStore theStore(2);
00393 theStore.Init(0,e1,nAngularParameters);
00394 theStore.Init(1,e2,nAngularParameters);
00395 for(i=0;i<nAngularParameters;i++)
00396 {
00397 theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
00398 theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
00399 }
00400
00401 theStore.SetManager(theManager);
00402 cosTh = theStore.SampleMax(fsEnergy);
00403 }
00404 }
00405 else
00406 {
00407 G4double x1 = running[it-1]/running[nEnergies-1];
00408 G4double x2 = running[it]/running[nEnergies-1];
00409 G4double y1 = theAngular[it-1].GetLabel();
00410 G4double y2 = theAngular[it].GetLabel();
00411 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
00412 random,x1,x2,y1,y2);
00413 G4NeutronHPLegendreStore theStore(2);
00414 theStore.Init(0,y1,nAngularParameters);
00415 theStore.Init(1,y2,nAngularParameters);
00416 theStore.SetManager(theManager);
00417 for(i=0;i<nAngularParameters;i++)
00418 {
00419 theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
00420 theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
00421 }
00422
00423 cosTh = theStore.SampleMax(fsEnergy);
00424 }
00425 delete [] running;
00426
00427
00428 remaining_energy -= fsEnergy;
00429
00430 }
00431 }
00432 else if(angularRep==2)
00433 {
00434
00435 G4int j;
00436 G4double * running = new G4double[nEnergies];
00437 running[0]=0;
00438 G4double weighted = 0;
00439 for(j=1; j<nEnergies; j++)
00440 {
00441 if(j!=0) running[j]=running[j-1];
00442 running[j] += theInt.GetBinIntegral(theManager.GetScheme(j-1),
00443 theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
00444 theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
00445 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(j-1),
00446 theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
00447 theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
00448 }
00449
00450
00451
00452 if ( nEnergies == 1 )
00453 currentMeanEnergy = 0.0;
00454 else
00455 currentMeanEnergy = weighted/running[nEnergies-1];
00456
00457 G4int itt(0);
00458 G4double randkal = G4UniformRand();
00459
00460
00461 for(j=1; j<nEnergies; j++)
00462 {
00463 itt = j;
00464 if(randkal<running[j]/running[nEnergies-1]) break;
00465 }
00466
00467
00468 G4double x, x1,x2,y1,y2;
00469 if(itt==0) itt=1;
00470 x = randkal*running[nEnergies-1];
00471 x1 = running[itt-1];
00472 x2 = running[itt];
00473 G4double compoundFraction;
00474
00475 y1 = theAngular[itt-1].GetLabel();
00476 y2 = theAngular[itt].GetLabel();
00477 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(itt-1),
00478 x, x1,x2,y1,y2);
00479
00480 G4double cLow = theAngular[itt-1].GetValue(1);
00481 G4double cHigh = theAngular[itt].GetValue(1);
00482 compoundFraction = theInt.Interpolate(theManager.GetScheme(itt),
00483 fsEnergy, y1, y2, cLow,cHigh);
00484 delete [] running;
00485
00486
00487 G4double incidentEnergy = anEnergy;
00488 G4double incidentMass = G4Neutron::Neutron()->GetPDGMass();
00489 G4double productEnergy = fsEnergy;
00490 G4double productMass = result->GetMass();
00491 G4int targetZ = G4int(theTargetCode/1000);
00492 G4int targetA = G4int(theTargetCode-1000*targetZ);
00493
00494 if ( targetA == 0 )
00495 targetA = G4int ( theTarget->GetMass()/amu_c2 + 0.5 );
00496 G4double targetMass = theTarget->GetMass();
00497 G4int residualA = targetA+1-A;
00498 G4int residualZ = targetZ-Z;
00499 G4double residualMass = residualZ*G4Proton::Proton()->GetPDGMass();
00500 residualMass +=(residualA-residualZ)*G4Neutron::Neutron()->GetPDGMass();
00501 residualMass -= G4NucleiProperties::GetBindingEnergy( residualA , residualZ );
00502 G4NeutronHPKallbachMannSyst theKallbach(compoundFraction,
00503 incidentEnergy, incidentMass,
00504 productEnergy, productMass,
00505 residualMass, residualA, residualZ,
00506 targetMass, targetA, targetZ);
00507 cosTh = theKallbach.Sample(anEnergy);
00508 }
00509 else if(angularRep>10&&angularRep<16)
00510 {
00511 G4double random = G4UniformRand();
00512 G4double * running = new G4double[nEnergies];
00513 running[0]=0;
00514 G4double weighted = 0;
00515 for(i=1; i<nEnergies; i++)
00516 {
00517 if(i!=0) running[i]=running[i-1];
00518 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
00519 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
00520 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
00521 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
00522 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
00523 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
00524 }
00525
00526
00527 if ( nEnergies == 1 )
00528 currentMeanEnergy = 0.0;
00529 else
00530 currentMeanEnergy = weighted/running[nEnergies-1];
00531
00532
00533 if ( nEnergies == 1 ) it = 0;
00534
00535 for(i=1; i<nEnergies; i++)
00536 {
00537 it = i;
00538 if(random<running[i]/running[nEnergies-1]) break;
00539 }
00540 if(it<nDiscreteEnergies||it==0)
00541 {
00542 if(it==0)
00543 {
00544 fsEnergy = theAngular[0].GetLabel();
00545 G4NeutronHPVector theStore;
00546 G4int aCounter = 0;
00547 for(G4int j=1; j<nAngularParameters; j+=2)
00548 {
00549 theStore.SetX(aCounter, theAngular[0].GetValue(j));
00550 theStore.SetY(aCounter, theAngular[0].GetValue(j+1));
00551 aCounter++;
00552 }
00553 G4InterpolationManager aMan;
00554 aMan.Init(angularRep-10, nAngularParameters-1);
00555 theStore.SetInterpolationManager(aMan);
00556 cosTh = theStore.Sample();
00557 }
00558 else
00559 {
00560 fsEnergy = theAngular[it].GetLabel();
00561 G4NeutronHPVector theStore;
00562 G4InterpolationManager aMan;
00563 aMan.Init(angularRep-10, nAngularParameters-1);
00564 theStore.SetInterpolationManager(aMan);
00565 G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
00566 G4int aCounter = 0;
00567 for(G4int j=1; j<nAngularParameters; j+=2)
00568 {
00569 theStore.SetX(aCounter, theAngular[it].GetValue(j));
00570 theStore.SetY(aCounter, theInt.Interpolate(currentScheme,
00571 random,
00572 running[it-1]/running[nEnergies-1],
00573 running[it]/running[nEnergies-1],
00574 theAngular[it-1].GetValue(j+1),
00575 theAngular[it].GetValue(j+1)));
00576 aCounter++;
00577 }
00578 cosTh = theStore.Sample();
00579 }
00580 }
00581 else
00582 {
00583 G4double x1 = running[it-1]/running[nEnergies-1];
00584 G4double x2 = running[it]/running[nEnergies-1];
00585 G4double y1 = theAngular[it-1].GetLabel();
00586 G4double y2 = theAngular[it].GetLabel();
00587 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
00588 random,x1,x2,y1,y2);
00589 G4NeutronHPVector theBuff1;
00590 G4NeutronHPVector theBuff2;
00591 G4InterpolationManager aMan;
00592 aMan.Init(angularRep-10, nAngularParameters-1);
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604 {
00605 G4int j;
00606 for(i=0,j=1; i<nAngularParameters; i++,j+=2)
00607 {
00608 theBuff1.SetX(i, theAngular[it-1].GetValue(j));
00609 theBuff1.SetY(i, theAngular[it-1].GetValue(j+1));
00610 theBuff2.SetX(i, theAngular[it].GetValue(j));
00611 theBuff2.SetY(i, theAngular[it].GetValue(j+1));
00612 }
00613 }
00614 G4NeutronHPVector theStore;
00615 theStore.SetInterpolationManager(aMan);
00616 x1 = y1;
00617 x2 = y2;
00618 G4double x, y;
00619
00620 for(i=0;i<theBuff1.GetVectorLength(); i++)
00621 {
00622 x = theBuff1.GetX(i);
00623 y1 = theBuff1.GetY(i);
00624 y2 = theBuff2.GetY(i);
00625 y = theInt.Interpolate(theManager.GetScheme(it),
00626 fsEnergy, theAngular[it-1].GetLabel(),
00627 theAngular[it].GetLabel(), y1, y2);
00628 theStore.SetX(i, x);
00629 theStore.SetY(i, y);
00630 }
00631 cosTh = theStore.Sample();
00632 }
00633 delete [] running;
00634 }
00635 else
00636 {
00637 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar::Sample: Unknown angular representation");
00638 }
00639 result->SetKineticEnergy(fsEnergy);
00640 G4double phi = twopi*G4UniformRand();
00641 G4double theta = std::acos(cosTh);
00642 G4double sinth = std::sin(theta);
00643 G4double mtot = result->GetTotalMomentum();
00644 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
00645 result->SetMomentum(tempVector);
00646
00647 return result;
00648 }