Geant4: G4NeutronHPLabAngularEnergy Class Reference

G4NeutronHPLabAngularEnergy::G4NeutronHPLabAngularEnergy ( ) [inline]

Definition at line 46 of file G4NeutronHPLabAngularEnergy.hh.

00047   {
00048     theEnergies = 0;
00049     theData = 0;
00050     nCosTh = 0;
00051     theSecondManager = 0;
00052   }

G4NeutronHPLabAngularEnergy::~G4NeutronHPLabAngularEnergy ( ) [inline]

Definition at line 53 of file G4NeutronHPLabAngularEnergy.hh.

00054   {
00055     if(theEnergies != 0) delete [] theEnergies;
00056     if(nCosTh != 0) delete [] nCosTh;
00057     if(theData != 0) 
00058     {
00059       for(G4int i=0; i<nEnergies; i++)
00060         delete [] theData[i];
00061       delete [] theData;
00062     }
00063     if(theSecondManager != 0) delete [] theSecondManager;
00064   }

void G4NeutronHPLabAngularEnergy::Init ( std::ifstream & aDataFile ) [virtual]

Implements G4VNeutronHPEnergyAngular.

Definition at line 46 of file G4NeutronHPLabAngularEnergy.cc.

References G4NeutronHPVector::Init(), G4InterpolationManager::Init(), and G4NeutronHPVector::SetLabel().

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 }

G4double G4NeutronHPLabAngularEnergy::MeanEnergyOfThisInteraction ( ) [inline, virtual]

Implements G4VNeutronHPEnergyAngular.

Definition at line 70 of file G4NeutronHPLabAngularEnergy.hh.

00071   {
00072     return currentMeanEnergy;
00073   }

G4ReactionProduct * G4NeutronHPLabAngularEnergy::Sample	(	G4double	anEnergy,
		G4double	massCode,
		G4double	mass
	)			`[virtual]`

Implements G4VNeutronHPEnergyAngular.

Definition at line 71 of file G4NeutronHPLabAngularEnergy.cc.

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    // get theta, E
00111    G4double cosTh, secEnergy;
00112    G4int i, it(0);
00113    // find the energy bin
00114    for(i=0; i<nEnergies; i++)
00115    {
00116      it = i;
00117      if ( anEnergy < theEnergies[i] ) break;
00118    }
00119    //080808
00120    //if ( it == 0 || it == nEnergies-1 ) // it marks the energy bin
00121    if ( it == 0 ) // it marks the energy bin
00122    {
00123 if(it==0) G4cout << "080808 Something unexpected is happen in G4NeutronHPLabAngularEnergy " << G4endl;
00124      // integrate the prob for each costh, and select theta.
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(); // Does interpolated integral.
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      //080807
00140      //if ( ith == 0 || ith == nCosTh[it]-1 ) //ith marks the angluar bin
00141      if ( ith == 0 ) //ith marks the angluar bin
00142      {
00143         cosTh = theData[it][ith].GetLabel();
00144         secEnergy = theData[it][ith].Sample();
00145         currentMeanEnergy = theData[it][ith].GetMeanX();
00146      }
00147      else
00148      {
00149        //080808
00150        //G4double x1 = theData[it][ith-1].GetIntegral();
00151        //G4double x2 = theData[it][ith].GetIntegral();
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 // this is the small big else.
00193    {
00194      G4double x, x1, x2, y1, y2, y, tmp, E;
00195      // integrate the prob for each costh, and select theta.
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      // get the distributions for the correct neutron energy
00213      x = anEnergy;
00214      x1 = theEnergies[it-1];
00215      x2 = theEnergies[it];
00216      G4NeutronHPVector thBuff1; // to be interpolated as run1.
00217      thBuff1.SetInterpolationManager(theSecondManager[it-1]);
00218      for(i=0; i<run1.GetVectorLength(); i++)
00219      {
00220        tmp = run1.GetX(i); //theta
00221        y1 = run1.GetY(i); // integral
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); //theta
00231        y1 = run1.GetY(tmp); // integral
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); // takes care of interpolation
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        // calculate theta
00248        G4double xx, xx1, xx2, yy1, yy2;
00249        xx =  random;
00250        xx1 = theThVec.GetY(ith-1)-theThVec.GetY(0); // integrals
00251        xx2 = theThVec.GetY(ith)-theThVec.GetY(0);
00252        yy1 = theThVec.GetX(ith-1); // std::cos(theta)
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      // get the indixes of the vectors close to theta for low energy
00259      // first it-1 !!!! i.e. low in energy
00260      for(i=0; i<nCosTh[it-1]; i++)
00261      {
00262        i1 = i;
00263        if(cosTh<theData[it-1][i].GetLabel()) break;
00264      }
00265      // now get the prob at this energy for the right theta value
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); // wrong E, right theta.
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); // wrong E, right theta.
00288      }
00289      G4NeutronHPVector theStore1;
00290      theStore1.Merge(&theBuff1a, &theBuff2a); // wrong E, right theta, complete binning
00291 
00292      // get the indixes of the vectors close to theta for high energy
00293      // then it !!!! i.e. high in energy
00294      for(i=0; i<nCosTh[it]; i++)
00295      {
00296        i2 = i;
00297        if(cosTh<theData[it][i2].GetLabel()) break;
00298      }                  // sonderfaelle mit i1 oder i2 head on fehlen. @@@@@
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); // wrong E, right theta.
00310      }
00311      G4NeutronHPVector theBuff2b;
00312      theBuff2b.SetInterpolationManager(theData[it][i2].GetInterpolationManager());
00313      //080808  i1 -> i2
00314      //for(i=0;i<theData[it][i1].GetVectorLength(); i++)
00315      for(i=0;i<theData[it][i2].GetVectorLength(); i++)
00316      {
00317        //E = theData[it][i1].GetX(i);
00318        //y1 = theData[it][i1-1].GetY(E);
00319        //y2 = theData[it][i1].GetY(i);
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); // wrong E, right theta.
00325      }
00326      G4NeutronHPVector theStore2;
00327      theStore2.Merge(&theBuff1b, &theBuff2b); // wrong E, right theta, complete binning
00328      // now get to the right energy.
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); // both correct
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); // both correct
00352      }
00353      G4NeutronHPVector theOne;
00354      theOne.Merge(&theOne1, &theOne2); // both correct, complete binning
00355 
00356      secEnergy = theOne.Sample();
00357      currentMeanEnergy = theOne.GetMeanX();
00358    }
00359 
00360 // now do random direction in phi, and fill the result.
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 }


Public Member Functions
	G4NeutronHPLabAngularEnergy ()
	~G4NeutronHPLabAngularEnergy ()
void	Init (std::ifstream &aDataFile)
G4ReactionProduct *	Sample (G4double anEnergy, G4double massCode, G4double mass)
G4double	MeanEnergyOfThisInteraction ()

G4NeutronHPLabAngularEnergy Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation