G4QIonIonCrossSection.cc

Go to the documentation of this 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 //
00027 // The lust update: M.V. Kossov, CERN/ITEP(Moscow) 19-Aug-07
00028 // GEANT4 tag $Name: not supported by cvs2svn $
00029 //
00030 //
00031 // G4 Physics class: G4QIonIonCrossSection for gamma+A cross sections
00032 // Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Dec-03
00033 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Feb-04
00034 // --------------------------------------------------------------------------------
00035 // ****************************************************************************************
00036 // This Header is a part of the CHIPS physics package (author: M. Kosov)
00037 // ****************************************************************************************
00038 // Short description: CHIPS cross-sectons for Ion-Ion interactions
00039 // ---------------------------------------------------------------
00040 //
00041 //#define debug
00042 //#define pdebug
00043 //#define debug3
00044 //#define debugn
00045 //#define debugs
00046 
00047 #include "G4QIonIonCrossSection.hh"
00048 #include "G4SystemOfUnits.hh"
00049 
00050 // Initialization of the
00051 G4double* G4QIonIonCrossSection::lastLENI=0;// Pointer to the lastArray of LowEn Inelast CS
00052 G4double* G4QIonIonCrossSection::lastHENI=0;// Pointer to the lastArray of HighEn InelastCS
00053 G4double* G4QIonIonCrossSection::lastLENE=0;// Pointer to the lastArray of LowEn Elastic CS
00054 G4double* G4QIonIonCrossSection::lastHENE=0;// Pointer to the lastArray of HighEn ElasticCS
00055 G4int     G4QIonIonCrossSection::lastPDG=0; // The last PDG code of the projectile
00056 G4int     G4QIonIonCrossSection::lastN=0;   // The last N of calculated nucleus
00057 G4int     G4QIonIonCrossSection::lastZ=0;   // The last Z of calculated nucleus
00058 G4double  G4QIonIonCrossSection::lastP=0.;  // Last used in cross section Momentum
00059 G4double  G4QIonIonCrossSection::lastTH=0.; // Last threshold momentum
00060 G4double  G4QIonIonCrossSection::lastICS=0.;// Last value of the Inelastic Cross Section
00061 G4double  G4QIonIonCrossSection::lastECS=0.;// Last value of the Elastic Cross Section
00062 G4int     G4QIonIonCrossSection::lastI=0;   // The last position in the DAMDB
00063 
00064 // Returns Pointer to the G4VQCrossSection class
00065 G4VQCrossSection* G4QIonIonCrossSection::GetPointer()
00066 {
00067   static G4QIonIonCrossSection theCrossSection; //**Static body of Cross Section**
00068   return &theCrossSection;
00069 }
00070 
00071 // The main member function giving the collision cross section (P is in IU, CS is in mb)
00072 // Make pMom in independent units !(Now it is MeV): fCS=true->Inelastic, fCS=false->Elastic
00073 G4double G4QIonIonCrossSection::GetCrossSection(G4bool fCS, G4double pMom, G4int tZ,
00074                                                 G4int tN, G4int pPDG)
00075 {
00076   static G4int j;                      // A#0f records found in DB for this projectile
00077   static std::vector <G4int>    colPDG;// Vector of the projectile PDG code
00078   static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
00079   static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
00080   static std::vector <G4double> colP;  // Vector of last momenta for the reaction
00081   static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
00082   static std::vector <G4double> colICS;// Vector of last inelastic cross-sections
00083   static std::vector <G4double> colECS;// Vector of last elastic cross-sections
00084   // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
00085 #ifdef pdebug
00086   G4cout<<"G4QIICS::GetCS:>>> f="<<fCS<<", Z="<<tZ<<"("<<lastZ<<"), N="<<tN<<"("<<lastN
00087         <<"),PDG="<<pPDG<<"("<<lastPDG<<"), p="<<pMom<<"("<<lastTH<<")"<<",Sz="
00088         <<colN.size()<<G4endl;
00089 #endif
00090   if(!pPDG)
00091   {
00092 #ifdef pdebug
00093     G4cout<<"G4QIonIonCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00094 #endif
00095     return 0.;                         // projectile PDG=0 is a mistake (?!) @@
00096   }
00097   G4bool in=false;                     // By default the isotope must be found in the AMDB
00098   if(tN!=lastN || tZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
00099   {
00100     in = false;                        // By default the isotope haven't be found in AMDB  
00101     lastP   = 0.;                      // New momentum history (nothing to compare with)
00102     lastPDG = pPDG;                    // The last PDG of the projectile
00103     lastN   = tN;                      // The last N of the calculated nucleus
00104     lastZ   = tZ;                      // The last Z of the calculated nucleus
00105     lastI   = colN.size();             // Size of the Associative Memory DB in the heap
00106     j  = 0;                            // A#0f records found in DB for this projectile
00107 #ifdef pdebug
00108     G4cout<<"G4QIICS::GetCS:FindI="<<lastI<<",pPDG="<<pPDG<<",tN="<<tN<<",tZ="<<tZ<<G4endl;
00109 #endif
00110     if(lastI) for(G4int i=0; i<lastI; i++) // Loop over all DB
00111     {                                  // The nucleus with projPDG is found in AMDB
00112 #ifdef pdebug
00113       G4cout<<"G4QII::GCS:P="<<colPDG[i]<<",N="<<colN[i]<<",Z="<<colZ[i]<<",j="<<j<<G4endl;
00114 #endif
00115       if(colPDG[i]==pPDG && colN[i]==tN && colZ[i]==tZ)
00116       {
00117         lastI=i;
00118         lastTH =colTH[i];                // Last THreshold (A-dependent)
00119 #ifdef pdebug
00120         G4cout<<"G4QIICS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00121 #endif
00122         if(pMom<=lastTH)
00123         {
00124 #ifdef pdebug
00125           G4cout<<"G4QIICS::GetCS:Found P="<<pMom<<"<Threshold="<<lastTH<<"->XS=0"<<G4endl;
00126 #endif
00127           return 0.;                     // Energy is below the Threshold value
00128         }
00129         lastP  =colP [i];                // Last Momentum  (A-dependent)
00130         lastICS=colICS[i];               // Last Inelastic Cross-Section (A-dependent)
00131         lastECS=colECS[i];               // Last Elastic Cross-Section (A-dependent)
00132         if(std::fabs(lastP/pMom-1.)<tolerance)
00133         {
00134 #ifdef pdebug
00135           G4cout<<"G4QIonIonCS::GetCS:P="<<pMom<<",InXS="<<lastICS*millibarn<<",ElXS="
00136                 <<lastECS*millibarn<<G4endl;
00137 #endif
00138           CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only
00139           if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
00140           return         lastECS*millibarn;     // Use theLastElasticCS
00141         }
00142         in = true;                       // This is the case when the isotop is found in DB
00143         // Momentum pMom is in IU ! @@ Units
00144 #ifdef pdebug
00145         G4cout<<"G4QIICS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00146 #endif
00147         lastICS=CalculateCrossSection( true,-1,j,lastPDG,lastZ,lastN,pMom);// read & update
00148         lastECS=CalculateCrossSection(false,-1,j,lastPDG,lastZ,lastN,pMom);// read & update
00149 #ifdef pdebug
00150         G4cout<<"G4QIonIonCS::GetCS:=>New(inDB) InCS="<<lastICS<<",ElCS="<<lastECS<<G4endl;
00151 #endif
00152         if((lastICS<=0. || lastECS<=0.) && pMom>lastTH) // Correct the threshold
00153         {
00154 #ifdef pdebug
00155           G4cout<<"G4QIonIonCS::GetCS:New,T="<<pMom<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00156 #endif
00157           lastTH=pMom;
00158         }
00159         break;                           // Go out of the LOOP
00160       }
00161 #ifdef pdebug
00162       G4cout<<"--->G4QIonIonCrossSec::GetCrosSec: pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
00163             <<",Z["<<i<<"]="<<colZ[i]<<",PDG="<<colPDG[i]<<G4endl;
00164 #endif
00165       j++;                             // Increment a#0f records found in DB for this pPDG
00166     }
00167     if(!in)                            // This nucleus has not been calculated previously
00168     {
00169 #ifdef pdebug
00170       G4cout<<"G4QIICS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
00171 #endif
00173       lastICS=CalculateCrossSection(true ,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create
00174       lastECS=CalculateCrossSection(false,0,j,lastPDG,lastZ,lastN,pMom); //calculate&create
00175       if(lastICS<=0. || lastECS<=0.)
00176       {
00177         lastTH = ThresholdEnergy(tZ, tN); // Threshold Energy=Mom=0 which is now the last
00178 #ifdef pdebug
00179         G4cout<<"G4QIonIonCrossSect::GetCrossSect:NewThresh="<<lastTH<<",P="<<pMom<<G4endl;
00180 #endif
00181         if(pMom>lastTH)
00182         {
00183 #ifdef pdebug
00184           G4cout<<"G4QIonIonCS::GetCS:1-st,P="<<pMom<<">Thresh="<<lastTH<<"->XS=0"<<G4endl;
00185 #endif
00186           lastTH=pMom;
00187         }
00188       }
00189 #ifdef pdebug
00190       G4cout<<"G4QIICS::GetCS: *New* ICS="<<lastICS<<", ECS="<<lastECS<<",N="<<lastN<<",Z="
00191             <<lastZ<<G4endl;
00192 #endif
00193       colN.push_back(tN);
00194       colZ.push_back(tZ);
00195       colPDG.push_back(pPDG);
00196       colP.push_back(pMom);
00197       colTH.push_back(lastTH);
00198       colICS.push_back(lastICS);
00199       colECS.push_back(lastECS);
00200 #ifdef pdebug
00201       G4cout<<"G4QIICS::GetCS:*1st*, P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn
00202             <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
00203 #endif
00204       if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
00205       return         lastECS*millibarn;     // Use theLastElasticCS
00206     } // End of creation of the new set of parameters
00207     else
00208     {
00209 #ifdef pdebug
00210       G4cout<<"G4QIICS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00211 #endif
00212       colP[lastI]=pMom;
00213       colPDG[lastI]=pPDG;
00214       colICS[lastI]=lastICS;
00215       colECS[lastI]=lastECS;
00216     }
00217   } // End of parameters udate
00218   else if(pMom<=lastTH)
00219   {
00220 #ifdef pdebug
00221     G4cout<<"G4QIICS::GetCS: Current T="<<pMom<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00222 #endif
00223     return 0.;                         // Momentum is below the Threshold Value -> CS=0
00224   }
00225   else if(std::fabs(lastP/pMom-1.)<tolerance)
00226   {
00227 #ifdef pdebug
00228     G4cout<<"G4QIICS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", InCS="<<lastICS*millibarn
00229           <<", ElCS="<<lastECS*millibarn<<"(mb)"<<G4endl;
00230 #endif
00231     if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
00232     return         lastECS*millibarn;     // Use theLastElasticCS
00233   }
00234   else
00235   {
00236 #ifdef pdebug
00237     G4cout<<"G4QIICS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00238 #endif
00239     lastICS=CalculateCrossSection( true,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
00240     lastECS=CalculateCrossSection(false,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
00241     lastP=pMom;
00242   }
00243 #ifdef pdebug
00244   G4cout<<"G4QIICS::GetCroSec:*End*,P="<<pMom<<"(MeV), InCS="<<lastICS*millibarn<<", ElCS="
00245         <<lastECS*millibarn<<"(mb)"<<G4endl;
00246 #endif
00247     if(fCS) return lastICS*millibarn;     // Use theLastInelasticCS
00248     return         lastECS*millibarn;     // Use theLastElasticCS
00249 }
00250 
00251 // The main member function giving the A-A cross section (Momentum in MeV, CS in mb)
00252 G4double G4QIonIonCrossSection::CalculateCrossSection(G4bool XS,G4int F,G4int I,G4int pPDG,
00253                                                       G4int tZ,G4int tN, G4double TotMom)
00254 {
00255   //static const G4double third=1./3.; // power for A^P->R conversion [R=1.1*A^(1/3)]
00256   //static const G4double conv=38.; // coeff. R2->sig=c*(pR+tR)^2, c=pi*10(mb/fm^2)*1.21
00257   // If change the following, please change in ::GetFunctions:
00258   static const G4double THmin=0.;  // @@ start from threshold (?) minimum Energy Threshold
00259   static const G4double dP=10.;    // step for the LEN table
00260   static const G4int    nL=100;    // A#of LENesonance points in E (each MeV from 2 to 106)
00261   static const G4double Pmin=THmin+(nL-1)*dP; // minE for the HighE part
00262   static const G4double Pmax=300000.;   // maxE for the HighE part
00263   static const G4int    nH=100;         // A#of HResonance points in lnE
00264   static const G4double milP=std::log(Pmin); // Low logarithm energy for the HighE part
00265   static const G4double malP=std::log(Pmax); // High logarithm energy (each 2.75 percent)
00266   static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part
00267   //
00268   // Associative memory for acceleration
00269   static std::vector <G4double*> LENI;   // Vector of pointers: LowEnIneIonIonCrossSection
00270   static std::vector <G4double*> HENI;   // Vector of pointers: HighEnIneIonIonCrossSection
00271   static std::vector <G4double*> LENE;   // Vector of pointers: LowEnElaIonIonCrossSection
00272   static std::vector <G4double*> HENE;   // Vector of pointers: HighEnElaIonIonCrossSection
00273 #ifdef debug
00274   G4cout<<"G4QIonIonCrossSection::CalcCS: Z="<<tZ<<", N="<<tN<<", P="<<TotMom<<G4endl;
00275 #endif
00276   G4int dPDG=pPDG/10;                // 10SZZZAAA
00277   G4int zPDG=dPDG/1000;              // 10SZZZ (?)
00278   G4int zA=dPDG%1000;                // proj A
00279   G4int pZ=zPDG%1000;                // proj Z (?)
00280   G4int pN=zA-pZ;                    // proj N (?)
00281   G4double Momentum=TotMom/zA;       // Momentum per nucleon
00282   if (Momentum<THmin) return 0.;     // @@ This can be dangerouse for the heaviest nuc.!
00283   G4double sigma=0.;
00284   if(F&&I) sigma=0.;                 // @@ *!* Fake line *!* to use F & I !!!Temporary!!!
00285   G4double tA=tN+tZ;                 // Target weight
00286   G4double pA=zA;                    // Projectile weight
00287   if(F<=0)                           // This isotope was not the last used isotop
00288   {
00289     if(F<0 || !XS)                   // This isotope was found in DAMDB or Elast =>RETRIEVE
00290     {
00291       lastLENI=LENI[I];              // Pointer to Low Energy inelastic cross sections
00292       lastHENI=HENI[I];              // Pointer to High Energy inelastic cross sections
00293       lastLENE=LENE[I];              // Pointer to Low Energy inelastic cross sections
00294       lastHENE=HENE[I];              // Pointer to High Energy inelastic cross sections
00295     }
00296     else                             // This isotope wasn't calculated previously => CREATE
00297     {
00298       lastLENI = new G4double[nL];   // Allocate memory for the new LEN cross sections
00299       lastHENI = new G4double[nH];   // Allocate memory for the new HEN cross sections
00300       lastLENE = new G4double[nL];   // Allocate memory for the new LEN cross sections
00301       lastHENE = new G4double[nH];   // Allocate memory for the new HEN cross sections
00302       G4int er=GetFunctions(pZ,pN,tZ,tN,lastLENI,lastHENI,lastLENE,lastHENE);
00303       if(er<1) G4cerr<<"*W*G4QIonIonCroSec::CalcCrossSection: pA="<<tA<<",tA="<<tA<<G4endl;
00304 #ifdef debug
00305       G4cout<<"G4QIonIonCrossSection::CalcCS: GetFunctions er="<<er<<",pA="<<pA<<",tA="<<tA
00306             <<G4endl;
00307 #endif
00308       // *** The synchronization check ***
00309       G4int sync=LENI.size();
00310       if(sync!=I) G4cout<<"*W*G4IonIonCrossSec::CalcCrossSect:Sync="<<sync<<"#"<<I<<G4endl;
00311       LENI.push_back(lastLENI);      // added LEN Inelastic
00312       HENI.push_back(lastHENI);      // added HEN Inelastic
00313       LENE.push_back(lastLENE);      // added LEN Elastic
00314       HENE.push_back(lastHENE);      // added HEN Elastic
00315     } // End of creation of the new set of parameters
00316   } // End of parameters udate
00317   // =------------= NOW the Magic Formula =--------------------------=
00318   if (Momentum<lastTH) return 0.;    // It must be already checked in the interface class
00319   else if (Momentum<Pmin)            // LEN region (approximated in E, not in lnE)
00320   {
00321 #ifdef debug
00322     G4cout<<"G4QIICS::CalCS:p="<<pA<<",t="<<tA<<",n="<<nL<<",T="<<THmin<<",d="<<dP<<G4endl;
00323 #endif
00324     if(tA<1. || pA<1.)
00325     {
00326       G4cout<<"-Warning-G4QIICS::CalcCS: pA="<<pA<<" or tA="<<tA<<" aren't nuclei"<<G4endl;
00327       sigma=0.;
00328     }
00329     else
00330     {
00331       G4double dPp=dP*pA;
00332       if(XS) sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENI);
00333       else   sigma=EquLinearFit(Momentum,nL,THmin,dPp,lastLENE);
00334     }
00335 #ifdef debugn
00336     if(sigma<0.) G4cout<<"-Warning-G4QIICS::CalcCS:pA="<<pA<<",tA="<<tA<<",XS="<<XS<<",P="
00337                        <<Momentum<<", Th="<<THmin<<", dP="<<dP<<G4endl;
00338 #endif
00339   }
00340   else if (Momentum<Pmax*pA)                     // High Energy region
00341   {
00342     G4double lP=std::log(Momentum);
00343 #ifdef debug
00344     G4cout<<"G4QIonIonCS::CalcCS:before HEN nH="<<nH<<",iE="<<milP<<",dlP="<<dlP<<G4endl;
00345 #endif
00346     if(tA<1. || pA<1.)
00347     {
00348       G4cout<<"-Warning-G4QIICS::CalCS:pA="<<pA<<" or tA="<<tA<<" aren't composit"<<G4endl;
00349       sigma=0.;
00350     }
00351     else
00352     {
00353       G4double milPp=milP+std::log(pA);
00354       if(XS) sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENI);
00355       else   sigma=EquLinearFit(lP,nH,milPp,dlP,lastHENE);
00356     }
00357   }
00358   else                                      // UltraHighE region (not frequent)
00359   {
00360     std::pair<G4double, G4double> inelel = CalculateXS(pZ, pN, tZ, tN, Momentum);
00361     if(XS) sigma=inelel.first;
00362     else   sigma=inelel.second;
00363   }
00364 #ifdef debug
00365   G4cout<<"G4IonIonCrossSection::CalculateCrossSection: sigma="<<sigma<<G4endl;
00366 #endif
00367   if(sigma<0.) return 0.;
00368   return sigma;
00369 }
00370 
00371 // Linear fit for YN[N] tabulated (from X0 with fixed step DX) function to X point
00372 
00373 // Calculate the functions for the log(A)
00374 G4int G4QIonIonCrossSection::GetFunctions(G4int pZ,G4int pN,G4int tZ,G4int tN,G4double* li,
00375                                           G4double* hi, G4double* le, G4double* he)
00376 {
00377   // If change the following, please change in ::CalculateCrossSection:
00378   static const G4double THmin=0.;  // @@ start from threshold (?) minimum Energy Threshold
00379   static const G4double dP=10.;    // step for the LEN table
00380   static const G4int    nL=100;    // A#of LENesonance points in E (each MeV from 2 to 106)
00381   static const G4double Pmin=THmin+(nL-1)*dP;   // minE for the HighE part
00382   static const G4double Pmax=300000.;           // maxE for the HighE part
00383   static const G4int    nH=100;                 // A#of HResonance points in lnE
00384   static const G4double milP=std::log(Pmin);    // Low logarithm energy for the HighE part
00385   static const G4double malP=std::log(Pmax);    // High logarithm energy
00386   static const G4double dlP=(malP-milP)/(nH-1); // Step in log energy in the HighE part
00387   static const G4double lP=std::exp(dlP);       // Multiplication factor in the HighE part
00388   // If the cross section approximation formula is changed - replace from file.
00389   if(pZ<1 || pN<0 || tZ<1 || tN<0)
00390   {
00391     G4cout<<"-W-G4QIonIonCS::GetFunct:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl;
00392     return -1;
00393   }
00394   G4int pA=pN+pZ;
00395   G4double dPp=dP*pA;
00396   G4double Mom=THmin;
00397   for(G4int k=0; k<nL; k++)
00398   {
00399     std::pair<G4double,G4double> len = CalculateXS(pZ, pN, tZ, tN, Mom);
00400     li[k]=len.first;
00401     le[k]=len.second;
00402     Mom+=dPp;
00403   }
00404   G4double lMom=Pmin*pA;
00405   for(G4int j=0; j<nH; j++)
00406   {
00407     std::pair<G4double,G4double> len = CalculateXS(pZ, pN, pZ, pN, lMom);
00408     hi[j]=len.first;
00409     he[j]=len.second;
00410     lMom*=lP;
00411   }
00412 #ifdef debug
00413   G4cout<<"G4QIonIonCS::GetFunctions: pZ="<<pZ<<", tZ="<<tZ<<" pair is calculated"<<G4endl;
00414 #endif
00415   return 1;
00416 }
00417 
00418 // Momentum (Mom=p/A) is in MeV/c, first=InelasticXS, second=ElasticXS (mb)
00419 std::pair<G4double,G4double> G4QIonIonCrossSection::CalculateXS(G4int pZ,G4int pN,G4int tZ,
00420                                                                 G4int tN, G4double Mom)
00421 {
00422   static G4VQCrossSection* PElCSman = G4QProtonElasticCrossSection::GetPointer();
00423   static G4VQCrossSection* NElCSman = G4QNeutronElasticCrossSection::GetPointer();
00424   static G4VQCrossSection* InelPCSman = G4QProtonNuclearCrossSection::GetPointer();
00425   static G4VQCrossSection* InelNCSman = G4QNeutronNuclearCrossSection::GetPointer();
00426   G4double pA=pZ+pN;
00427   G4double tA=tZ+tN;
00428   if(pA<.9 || tA<.9 ||pA>239. || tA>239 || Mom < 0.) return std::make_pair(0.,0.);
00429   G4double inCS=0.;
00430   G4double elCS=0.;
00431   if(pA<1.1 )               // nucleon-ion interaction use NA(in,el)
00432   {
00433     if     (pZ == 1 && !pN) // proton-nuclear
00434     {
00435       inCS=InelPCSman->GetCrossSection(true, Mom, tZ, tN, 2212);
00436       elCS=PElCSman->GetCrossSection(true, Mom, tZ, tN, 2212);
00437     }
00438     else if(pN == 1 && !pZ) // neutron-nuclear
00439     {
00440       inCS=InelNCSman->GetCrossSection(true, Mom, tZ, tN, 2112);
00441       elCS=NElCSman->GetCrossSection(true, Mom, tZ, tN, 2112);
00442     }
00443     else G4cerr<<"-Warn-G4QIICS::CaCS:pZ="<<pZ<<",pN="<<pN<<",tZ="<<tZ<<",tN="<<tN<<G4endl;
00444   }
00445   else
00446   {
00447     G4double T=ThresholdMomentum(pZ, pN, tZ, tN); // @@ Can be cashed as lastTH (?)
00448     if(Mom<=T)
00449     {
00450       elCS=0.;
00451       inCS=0.;
00452     }
00453     else
00454     {
00455       G4double P2=Mom*Mom;
00456       G4double T2=T*T;
00457       G4double R=1.-T2/P2;                        // @@ Very rough threshold effect
00458       //G4double P4=P2*P2;
00459       //G4double P8=P4*P4;
00460       //G4double T4=T2*T2;
00461       //G4double tot=CalculateTotal(pA, tA, Mom)*P8/(P8+T4*T4); // @@ convert to IndepUnits
00462       G4double tot=R*CalculateTotal(pA, tA, Mom); // @@ convert to IndepUnits
00463       G4double rat=CalculateElTot(pA, tA, Mom);
00464       elCS=tot*rat;
00465       inCS=tot-elCS;
00466     }
00467   }
00468   return std::make_pair(inCS,elCS);
00469 }
00470 
00471 // Total Ion-ion cross-section (mb), Momentum (Mom) here is p/A (MeV/c=IU) (No Threshold)
00472 G4double G4QIonIonCrossSection::CalculateTotal(G4double pA, G4double tA, G4double Mom)
00473 {
00474   G4double y=std::log(Mom/1000.); // Log of momentum in GeV/c
00475   G4double ab=pA+tA;
00476   G4double al=std::log(ab);
00477   G4double ap=std::log(pA*tA);
00478   G4double e=std::pow(pA,0.1)+std::pow(tA,0.1);
00479   G4double d=e-1.55/std::pow(al,0.2);
00480   G4double f=4.;
00481   if(pA>4. && tA>4.) f=3.3+130./ab/ab+2.25/e;
00482   G4double c=(385.+11.*ab/(1.+.02*ab*al)+(.5*ab+500./al/al)/al)*std::pow(d,f);
00483   G4double r=y-3.-4./ap/ap;
00484 #ifdef pdebug
00485   G4cout<<"G4QIonIonCS::CalcTot:P="<<Mom<<", stot="<<c+d*r*r<<", d="<<d<<", r="<<r<<G4endl;
00486 #endif
00487   return c+d*r*r;
00488 }
00489 
00490 // Ratio elastic/Total, Momentum (Mom) here is p/A (MeV/c=IU)
00491 G4double G4QIonIonCrossSection::CalculateElTot(G4double pA, G4double tA, G4double Mom)
00492 {
00493   G4double y=std::log(Mom/1000.); // Log of momentum in GeV/c
00494   G4double ab=pA*tA;
00495   G4double ap=std::log(ab);
00496   G4double r=y-3.92-1.73/ap/ap;
00497   G4double d=.00166/(1.+.002*std::pow(ab,1.33333));
00498   G4double al=std::log(pA+tA);
00499   G4double e=1.+.235*(std::fabs(ap-5.78));
00500   if     (std::fabs(pA-2.)<.1 && std::fabs(tA-2.)<.1) e=2.18;
00501   else if(std::fabs(pA-2.)<.1 && std::fabs(tA-3.)<.1) e=1.90;
00502   else if(std::fabs(pA-3.)<.1 && std::fabs(tA-2.)<.1) e=1.90;
00503   else if(std::fabs(pA-2.)<.1 && std::fabs(tA-4.)<.1) e=1.65;
00504   else if(std::fabs(pA-4.)<.1 && std::fabs(tA-2.)<.1) e=1.65;
00505   else if(std::fabs(pA-3.)<.1 && std::fabs(tA-4.)<.1) e=1.32;
00506   else if(std::fabs(pA-4.)<.1 && std::fabs(tA-3.)<.1) e=1.32;
00507   else if(std::fabs(pA-4.)<.1 && std::fabs(tA-4.)<.1) e=1.;
00508   G4double f=.37+.0188*al;
00509   G4double g_value=std::log(std::pow(pA,0.35)+std::pow(tA,0.35));
00510   G4double h=g_value*g_value;
00511   G4double c=f/(1.+e/h/h);
00512 #ifdef pdebug
00513   G4cout<<"G4QIonIonCS::CalcElT:P="<<Mom<<",el/tot="<<c+d*r*r<<",d="<<d<<", r="<<r<<G4endl;
00514 #endif
00515   return c+d*r*r;
00516 }
00517 
00518 // Electromagnetic momentum/A-threshold (in MeV/c) 
00519 G4double G4QIonIonCrossSection::ThresholdMomentum(G4int pZ, G4int pN, G4int tZ, G4int tN)
00520 {
00521   static const G4double third=1./3.;
00522   static const G4double pM = G4QPDGCode(2212).GetMass(); // Proton mass in MeV
00523   static const G4double tpM= pM+pM;       // Doubled proton mass (MeV)
00524   if(pZ<.99 || pN<0. || tZ<.99 || tN<0.) return 0.;
00525   G4double tA=tZ+tN;
00526   G4double pA=pZ+pN;
00527   //G4double dE=1.263*tZ/(1.+std::pow(tA,third));
00528   G4double dE=pZ*tZ/(std::pow(pA,third)+std::pow(tA,third))/pA; // dE/pA (per projNucleon)
00529   return std::sqrt(dE*(tpM+dE));
00530 }

Generated on Mon May 27 17:49:36 2013 for Geant4 by  doxygen 1.4.7