Geant4: G4NeutronHPDiscreteTwoBody.cc Source File

00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // neutron_hp -- source file
00027 // J.P. Wellisch, Nov-1996
00028 // A prototype of the low energy neutron transport model.
00029 //
00030 //080612 Bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2,3
00031 //080709 Bug fix Sampling Legendre expansion by T. Koi   
00032 //101110 Bug fix in MF=6, LAW=2 case; contribution from E. Mendoza, D. Cano-Ott (CIEMAT)
00033 //
00034 #include "G4NeutronHPDiscreteTwoBody.hh"
00035 #include "G4PhysicalConstants.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 #include "G4NeutronHPVector.hh"
00046 #include "G4NeutronHPLegendreStore.hh"
00047 
00048 G4ReactionProduct * G4NeutronHPDiscreteTwoBody::Sample(G4double anEnergy, G4double massCode, G4double )
00049 { // Interpolation still only for the most used parts; rest to be Done @@@@@
00050    G4ReactionProduct * result = new G4ReactionProduct;
00051    G4int Z = static_cast<G4int>(massCode/1000);
00052    G4int A = static_cast<G4int>(massCode-1000*Z);
00053 
00054    if(massCode==0)
00055    {
00056      result->SetDefinition(G4Gamma::Gamma());
00057    }
00058    else if(A==0)
00059    {
00060      result->SetDefinition(G4Electron::Electron());     
00061      if(Z==1) result->SetDefinition(G4Positron::Positron());
00062    }
00063    else if(A==1)
00064    {
00065      result->SetDefinition(G4Neutron::Neutron());
00066      if(Z==1) result->SetDefinition(G4Proton::Proton());
00067    }
00068    else if(A==2)
00069    {
00070      result->SetDefinition(G4Deuteron::Deuteron());      
00071    }
00072    else if(A==3)
00073    {
00074      result->SetDefinition(G4Triton::Triton());  
00075      if(Z==2) result->SetDefinition(G4He3::He3());
00076    }
00077    else if(A==4)
00078    {
00079      result->SetDefinition(G4Alpha::Alpha());
00080      if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");    
00081    }
00082    else
00083    {
00084      throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPDiscreteTwoBody: Unknown ion case 2");
00085    }
00086    
00087 // get cosine(theta)
00088    G4int i(0), it(0);
00089    G4double cosTh(0);
00090    for(i=0; i<nEnergy; i++)
00091    {
00092      it = i;
00093      if(theCoeff[i].GetEnergy()>anEnergy) break;
00094    }
00095    if(it==0||it==nEnergy-1)
00096    {
00097      if(theCoeff[it].GetRepresentation()==0)
00098      {
00099 //TK Legendre expansion
00100        G4NeutronHPLegendreStore theStore(1);
00101        theStore.SetCoeff(0, theCoeff);
00102        theStore.SetManager(theManager);
00103        //cosTh = theStore.SampleMax(anEnergy);
00104        //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
00105        cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
00106      }
00107      else if(theCoeff[it].GetRepresentation()==12) // means LINLIN
00108      {
00109        G4NeutronHPVector theStore; 
00110        G4InterpolationManager aManager;
00111        aManager.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
00112        theStore.SetInterpolationManager(aManager);
00113        for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
00114        {
00115          //101110
00116          //theStore.SetX(i, theCoeff[it].GetCoeff(i));
00117          //theStore.SetY(i, theCoeff[it].GetCoeff(i));
00118          theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
00119          theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
00120          i++;
00121        }
00122        cosTh = theStore.Sample();
00123      }
00124      else if(theCoeff[it].GetRepresentation()==14) //this is LOGLIN
00125      {
00126        G4NeutronHPVector theStore;
00127        G4InterpolationManager aManager;
00128        aManager.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
00129        theStore.SetInterpolationManager(aManager);
00130        for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
00131        {
00132          //101110
00133          //theStore.SetX(i, theCoeff[it].GetCoeff(i));
00134          //theStore.SetY(i, theCoeff[it].GetCoeff(i));
00135          theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
00136          theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
00137          i++;
00138        }
00139        cosTh = theStore.Sample(); 
00140      }
00141      else
00142      {
00143        throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering");
00144      }
00145    }
00146    else
00147    {
00148      if(theCoeff[it].GetRepresentation() == theCoeff[it-1].GetRepresentation())
00149      {
00150        if(theCoeff[it].GetRepresentation()==0)
00151        {
00152 //TK Legendre expansion
00153          G4NeutronHPLegendreStore theStore(2);
00154          theStore.SetCoeff(0, &(theCoeff[it-1]));
00155          theStore.SetCoeff(1, &(theCoeff[it]));
00156          G4InterpolationManager aManager;
00157          aManager.Init(theManager.GetScheme(it), 2);
00158          theStore.SetManager(aManager);
00159          //cosTh = theStore.SampleMax(anEnergy);
00160 //080709 TKDB
00161          cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
00162        }
00163        else if(theCoeff[it].GetRepresentation()==12) // LINLIN
00164        {
00165          G4NeutronHPVector theBuff1;
00166          G4InterpolationManager aManager1;
00167          aManager1.Init(LINLIN, theCoeff[it-1].GetNumberOfPoly()/2);
00168          theBuff1.SetInterpolationManager(aManager1);
00169          for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i++)
00170          {
00171            //101110
00172            //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
00173            //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
00174            theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
00175            theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
00176            i++;
00177          }
00178          G4NeutronHPVector theBuff2;
00179          G4InterpolationManager aManager2;
00180          aManager2.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
00181          theBuff2.SetInterpolationManager(aManager2);
00182          for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
00183          {
00184            //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
00185            //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
00186            theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
00187            theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
00188            i++;
00189          }
00190 
00191          G4double x1 = theCoeff[it-1].GetEnergy();
00192          G4double x2 = theCoeff[it].GetEnergy();
00193          G4double x = anEnergy;
00194          G4double y1, y2, y, mu;
00195 
00196          G4NeutronHPVector theStore1;
00197          theStore1.SetInterpolationManager(aManager1);
00198          G4NeutronHPVector theStore2;
00199          theStore2.SetInterpolationManager(aManager2);
00200          G4NeutronHPVector theStore;
00201          
00202          // for fixed mu get p1, p2 and interpolate according to x
00203          for(i=0; i<theBuff1.GetVectorLength(); i++)
00204          {
00205            mu = theBuff1.GetX(i);
00206            y1 = theBuff1.GetY(i);
00207            y2 = theBuff2.GetY(mu);
00208            y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00209            theStore1.SetData(i, mu, y);
00210          }
00211          for(i=0; i<theBuff2.GetVectorLength(); i++)
00212          {
00213            mu = theBuff2.GetX(i);
00214            y1 = theBuff2.GetY(i);
00215            y2 = theBuff1.GetY(mu);
00216            y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00217            theStore2.SetData(i, mu, y);
00218          }
00219          theStore.Merge(&theStore1, &theStore2); // merge takes care of interpolationschemes
00220          cosTh = theStore.Sample();
00221        }
00222        else if(theCoeff[it].GetRepresentation()==14) //TK LOG_LIN
00223        {
00224          G4NeutronHPVector theBuff1;
00225          G4InterpolationManager aManager1;
00226          aManager1.Init(LOGLIN, theCoeff[it-1].GetNumberOfPoly()/2);
00227          theBuff1.SetInterpolationManager(aManager1);
00228          for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i++)
00229          {
00230            //101110
00231            //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
00232            //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
00233            theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
00234            theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
00235            i++;
00236          }
00237          
00238          G4NeutronHPVector theBuff2;
00239          G4InterpolationManager aManager2;
00240          aManager2.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
00241          theBuff2.SetInterpolationManager(aManager2);
00242          for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
00243          {
00244            //101110
00245            //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
00246            //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
00247            theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
00248            theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
00249            i++;
00250          }
00251 
00252          G4double x1 = theCoeff[it-1].GetEnergy();
00253          G4double x2 = theCoeff[it].GetEnergy();
00254          G4double x = anEnergy;
00255          G4double y1, y2, y, mu;
00256 
00257          G4NeutronHPVector theStore1;
00258          theStore1.SetInterpolationManager(aManager1);
00259          G4NeutronHPVector theStore2;
00260          theStore2.SetInterpolationManager(aManager2);
00261          G4NeutronHPVector theStore;
00262          
00263          // for fixed mu get p1, p2 and interpolate according to x
00264          for(i=0; i<theBuff1.GetVectorLength(); i++)
00265          {
00266            mu = theBuff1.GetX(i);
00267            y1 = theBuff1.GetY(i);
00268            y2 = theBuff2.GetY(mu);
00269            y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00270            theStore1.SetData(i, mu, y);
00271          }
00272          for(i=0; i<theBuff2.GetVectorLength(); i++)
00273          {
00274            mu = theBuff2.GetX(i);
00275            y1 = theBuff2.GetY(i);
00276            y2 = theBuff1.GetY(mu);
00277            y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
00278            theStore2.SetData(i, mu, y);
00279          }
00280          theStore.Merge(&theStore1, &theStore2); 
00281          cosTh = theStore.Sample(); 
00282        }
00283        else
00284        {
00285          throw G4HadronicException(__FILE__, __LINE__, "Two neighbouring distributions with different interpolation");
00286        }
00287      }
00288      else
00289      {
00290        throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering, case 2");
00291      }
00292    }
00293    
00294 // now get the energy from kinematics and Q-value.
00295 
00296    //G4double restEnergy = anEnergy+GetQValue();
00297    
00298 // assumed to be in CMS @@@@@@@@@@@@@@@@@
00299 
00300    //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2
00301    //G4double residualMass =   GetTarget()->GetMass() + GetNeutron()->GetMass()
00302    //                        - result->GetMass() - GetQValue();
00303    //G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
00304    G4double A1     =  GetTarget()->GetMass()/GetNeutron()->GetMass(); 
00305    G4double A1prim =  result->GetMass()/GetNeutron()->GetMass();
00306    G4double E1     =  (A1+1)*(A1+1)/A1/A1*anEnergy; 
00307    G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*E1+(1+A1)*GetQValue());
00308 
00309    result->SetKineticEnergy(kinE); // non relativistic @@
00310    G4double phi = twopi*G4UniformRand();
00311    G4double theta = std::acos(cosTh);
00312    G4double sinth = std::sin(theta);
00313    G4double mtot = result->GetTotalMomentum(); 
00314    G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
00315    result->SetMomentum(tempVector);
00316    
00317 // some garbage collection
00318    
00319 // return the result   
00320    return result;
00321 }