G4QMuonNuclearCrossSection Class Reference

#include <G4QMuonNuclearCrossSection.hh>

Inheritance diagram for G4QMuonNuclearCrossSection:

G4VQCrossSection

Public Member Functions

 ~G4QMuonNuclearCrossSection ()
G4double ThresholdEnergy (G4int Z, G4int N, G4int PDG=13)
virtual G4double GetCrossSection (G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=0)
G4double CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum)
G4int GetExchangePDGCode ()
G4double GetExchangeEnergy ()
G4double GetVirtualFactor (G4double nu, G4double Q2)
G4double GetExchangeQ2 (G4double nu)

Static Public Member Functions

static G4VQCrossSectionGetPointer ()

Protected Member Functions

 G4QMuonNuclearCrossSection ()

Detailed Description

Definition at line 56 of file G4QMuonNuclearCrossSection.hh.


Constructor & Destructor Documentation

G4QMuonNuclearCrossSection::G4QMuonNuclearCrossSection (  )  [inline, protected]

Definition at line 60 of file G4QMuonNuclearCrossSection.hh.

00060 {};

G4QMuonNuclearCrossSection::~G4QMuonNuclearCrossSection (  ) 

Definition at line 83 of file G4QMuonNuclearCrossSection.cc.

00084 {
00085   G4int lens=J1->size();
00086   for(G4int i=0; i<lens; ++i) delete[] (*J1)[i];
00087   delete J1;
00088   G4int hens=J2->size();
00089   for(G4int i=0; i<hens; ++i) delete[] (*J2)[i];
00090   delete J2;
00091   G4int pens=J3->size();
00092   for(G4int i=0; i<pens; ++i) delete[] (*J3)[i];
00093   delete J3;
00094 }


Member Function Documentation

G4double G4QMuonNuclearCrossSection::CalculateCrossSection ( G4bool  CS,
G4int  F,
G4int  I,
G4int  PDG,
G4int  Z,
G4int  N,
G4double  Momentum 
) [virtual]

Implements G4VQCrossSection.

Definition at line 322 of file G4QMuonNuclearCrossSection.cc.

References G4cerr, G4cout, and G4endl.

Referenced by GetCrossSection().

00324 {
00325   static const G4int nE=336;    // !! If change this, change it in GetFunctions() (*.hh) !!
00326   static const G4int mL=nE-1;
00327   static const G4double EMi=2.0612;          // Minimum tabulated KineticEnergy of Muon
00328   static const G4double EMa=50000.;          // Maximum tabulated Energy of the Muon 
00329   static const G4double lEMi=std::log(EMi);  // Minimum tabulatedLogarithmKinEnergy of Muon
00330   static const G4double lEMa=std::log(EMa);  // Maximum tabulatedLogarithmKinEnergy of Muon
00331   static const G4double dlnE=(lEMa-lEMi)/mL; // Logarithmic table step for MuonKinEnergy
00332   static const G4double alop=1./137.036/3.14159265; //coeffitient for  E>50000 calculations
00333   static const G4double mmu=105.65839;       // Mass of the muon in MeV
00334   static const G4double mmu2=mmu*mmu;        // Squared Mass of muon in MeV^2
00335   static const G4double lmmu=std::log(mmu);       // Log of the muon mass
00336   // *** Begin of the Associative memory for acceleration of the cross section calculations
00337   static std::vector <G4int> colF;     // Vector of LastStartPosition in Ji-function tables
00338   static std::vector <G4double> colH;  // Vector of HighEnergyCoefficients (functional)
00339 #ifdef pdebug
00340   G4cout<<"G4QMuonNucCrossSection::CalculateCrossSection: ***Called*** "<<J3->size();
00341   if(J3->size()) G4cout<<", p="<<(*J3)[0];
00342   G4cout<<G4endl;
00343 #endif
00344   // *** End of Static Definitions (Associative Memory) ***
00345   //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
00346   onlyCS=CS;                           // Flag to calculate only CS (not Si/Bi)
00347   G4double TotEnergy2=Momentum*Momentum+mmu2;
00348   G4double TotEnergy=std::sqrt(TotEnergy2); // Total energy of the muon
00349   lastE=TotEnergy-mmu;               // Kinetic energy of the muon
00350 #ifdef pdebug
00351   G4cout<<"G4QMuonNucCS::CalcCS: P="<<Momentum<<", F="<<F<<", I="<<I<<", Z="<<targZ
00352         <<", N="<<targN<<", onlyCS="<<CS<<",E="<<lastE<<",th="<<EMi<<G4endl;
00353 #endif
00354   G4double A=targN+targZ;            // New A (can be different from targetAtomicNumber)
00355   if(F<=0)                           // >>>This isotope was not the last used isotop<<<
00356   {
00357     if(F<0)                          // This isotope was found in DAMDB =------=> RETRIEVE
00358     {                                // ...........................................=------=
00359       if (lastE<=EMi)                // Energy is below the minimum energy in the table
00360       {
00361         lastE=0.;
00362         lastG=0.;
00363         lastSig=0.;
00364 #ifdef pdebug
00365         G4cout<<"--> G4QMuonNucCS::CalcCS: Old CS=0  as lastE="<<lastE<<" < "<<EMi<<G4endl;
00366 #endif
00367         return 0.;
00368       }
00369       lastJ1 =(*J1)[I];              // Pointer to the prepared J1 function
00370       lastJ2 =(*J2)[I];              // Pointer to the prepared J2 function
00371       lastJ3 =(*J3)[I];              // Pointer to the prepared J3 function
00372       lastF  =colF[I];               // Last ZeroPosition in the J-functions
00373       lastH  =colH[I];               // Last High Energy Coefficient (A-dependent)
00374     }
00375     else                             // This isotope wasn't calculated previously => CREATE
00376     {                                // ............................................=-----=
00377       lastJ1 = new G4double[nE];     // Allocate memory for the new J1 function
00378       lastJ2 = new G4double[nE];     // Allocate memory for the new J2 function
00379       lastJ3 = new G4double[nE];     // Allocate memory for the new J3 function
00380       lastF   = GetFunctions(A,lastJ1,lastJ2,lastJ3);//newZeroPos and J-functions filling
00381       lastH   = alop*A*(1.-.072*std::log(A)); // like lastSP of G4PhotonuclearCrossSection
00382 #ifdef pdebug
00383       G4cout<<"==>G4QMuNCS::CalcCS:lJ1="<<lastJ1<<",lJ2="<<lastJ2<<",lJ3="<<lastJ3<<G4endl;
00384 #endif
00385       // *** The synchronization check ***
00386       G4int sync=J1->size();
00387       if(sync!=I) G4cerr<<"***G4QMuonNuclearCS::CalcCrS: PDG=13, S="<<sync<<"#"<<I<<G4endl;
00388       J1->push_back(lastJ1);
00389       J2->push_back(lastJ2);
00390       J3->push_back(lastJ3);
00391       colF.push_back(lastF);
00392       colH.push_back(lastH);
00393     } // End of creation of the new set of parameters
00394   } // End of parameters udate
00395   // =-------------------= NOW Calculate the Cross Section =--------------------------=
00396   if (lastE<=lastTH || lastE<=EMi)   // Check that muKiE is higher than ThreshE
00397   {
00398     lastE=0.;
00399     lastG=0.;
00400     lastSig=0.;
00401 #ifdef pdebug
00402     G4cout<<"--> G4QMuonNucCS::CalcCS:CS=0 as T="<<lastE<<"<"<<EMi<<" || "<<lastTH<<G4endl;
00403 #endif
00404     return 0.;
00405   }
00406   G4double lE=std::log(lastE);       // log(muE) (it is necessary for the fit)
00407   lastG=lE-lmmu;                     // Gamma of the muon (used to recover log(muE))
00408   G4double dlg1=lastG+lastG-1.;
00409   G4double lgoe=lastG/lastE;
00410   if(lE<lEMa) // Log fit is made explicitly to fix the last bin for the randomization
00411   {
00412     G4double shift=(lE-lEMi)/dlnE;
00413     G4int    blast=static_cast<int>(shift);
00414 #ifdef pdebug
00415     G4cout<<"-->G4QMuonNuclearCS::CalcCrossSect:LOGfit b="<<blast<<",max="<<mL<<",lJ1="
00416           <<lastJ1<<",lJ2="<<lastJ2<<",lJ3="<<lastJ3<<G4endl;
00417 #endif
00418     if(blast<0)   blast=0;
00419     if(blast>=mL) blast=mL-1;
00420     shift-=blast;
00421     lastL=blast+1;
00422     G4double YNi=dlg1*lastJ1[blast]-lgoe*(lastJ2[blast]+lastJ2[blast]-lastJ3[blast]/lastE);
00423     G4double YNj=dlg1*lastJ1[lastL]-lgoe*(lastJ2[lastL]+lastJ2[lastL]-lastJ3[lastL]/lastE);
00424     lastSig = YNi+shift*(YNj-YNi);
00425     if(lastSig>YNj)lastSig=YNj;
00426 #ifdef pdebug
00427     G4cout<<"G4QMuNucCS::CalcCS:S="<<lastSig<<",E="<<lE<<",Yi="<<YNi<<",Yj="<<YNj<<",M="
00428           <<lEMa<<G4endl;
00429     G4cout<<"G4QMuNucCS::CalcCS:s="<<shift<<",Jb="<<lastJ1[blast]<<",J="<<lastJ1[lastL]
00430           <<",b="<<blast<<G4endl;
00431 #endif
00432   }
00433   else
00434   {
00435 #ifdef pdebug
00436     G4cout<<"->G4QMuonNucCS::CCS:LOGex="<<lastJ1<<",lJ2="<<lastJ2<<",lJ3="<<lastJ3<<G4endl;
00437 #endif
00438     lastL=mL;
00439     G4double term1=lastJ1[mL]+lastH*HighEnergyJ1(lE);
00440     G4double term2=lastJ2[mL]+lastH*HighEnergyJ2(lE);
00441     G4double term3=lastJ3[mL]+lastH*HighEnergyJ3(lE);
00442     lastSig=dlg1*term1-lgoe*(term2+term2-term3/lastE);
00443 #ifdef pdebug
00444     G4cout<<"G4QMuNucCS::CalculateCrossSection:S="<<lastSig<<",lE="<<lE<<",J1="
00445           <<lastH*HighEnergyJ1(lE)<<",Pm="<<lastJ1[mL]<<",Fm="<<lastJ2[mL]<<",Fh="
00446           <<lastH*HighEnergyJ2(lE)<<",EM="<<lEMa<<G4endl;
00447 #endif
00448   }
00449   if(lastSig<0.) lastSig = 0.;
00450   return lastSig;
00451 }

G4double G4QMuonNuclearCrossSection::GetCrossSection ( G4bool  fCS,
G4double  pMom,
G4int  tgZ,
G4int  tgN,
G4int  pPDG = 0 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 98 of file G4QMuonNuclearCrossSection.cc.

References CalculateCrossSection(), G4cout, G4endl, and ThresholdEnergy().

00100 {
00101   static const G4double mmu=105.65839; // Mass of the muon in MeV
00102   static const G4double mmu2=mmu*mmu;  // Squared Mass of muon in MeV^2
00103   static G4int j;                      // A#0f records found in DB for this projectile
00104   static std::vector <G4int>    colPDG;// Vector of the projectile PDG code
00105   static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
00106   static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
00107   static std::vector <G4double> colP;  // Vector of last momenta for the reaction
00108   static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
00109   static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
00110   // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
00111   G4double pEn=std::sqrt(pMom*pMom+mmu2)-mmu; // ==> mu-/mu+ kinEnergy
00112 #ifdef pdebug
00113   G4cout<<"G4QMNCS::GetCS:>>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
00114         <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="
00115         <<colN.size()<<G4endl;
00116   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00117 #endif
00118   if(std::abs(pPDG)!=13)
00119   {
00120 #ifdef pdebug
00121     G4cout<<"G4QMNCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00122     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00123 #endif
00124     return 0.;                         // projectile PDG=0 is a mistake (?!) @@
00125   }
00126   G4bool in=false;                     // By default the isotope must be found in the AMDB
00127   if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
00128   {
00129     in = false;                        // By default the isotope haven't be found in AMDB  
00130     lastP   = 0.;                      // New momentum history (nothing to compare with)
00131     lastPDG = pPDG;                    // The last PDG of the projectile
00132     lastN   = tgN;                     // The last N of the calculated nucleus
00133     lastZ   = tgZ;                     // The last Z of the calculated nucleus
00134     lastI   = colN.size();             // Size of the Associative Memory DB in the heap
00135     j  = 0;                            // A#0f records found in DB for this projectile
00136     if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found
00137     {                                  // The nucleus with projPDG is found in AMDB
00138       if(colN[i]==tgN && colZ[i]==tgZ)
00139       {
00140         lastI=i;
00141         lastTH =colTH[i];                // Last THreshold (A-dependent)
00142 #ifdef pdebug
00143         G4cout<<"G4QMNCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00144         //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00145 #endif
00146         if(pEn<=lastTH)
00147         {
00148 #ifdef pdebug
00149           G4cout<<"G4QMNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",CS=0"<<G4endl;
00150           //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00151 #endif
00152           return 0.;                     // Energy is below the Threshold value
00153         }
00154         lastP  =colP [i];                // Last Momentum  (A-dependent)
00155         lastCS =colCS[i];                // Last CrossSect (A-dependent)
00156  //        if(std::fabs(lastP/pMom-1.)<tolerance)
00157         if(lastP==pMom)                  // VI do not use tolerance
00158         {
00159 #ifdef pdebug
00160           G4cout<<"G4QMNCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00161 #endif
00162           CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // Update param's only
00163           return lastCS*millibarn;     // Use theLastCS
00164         }
00165         in = true;                       // This is the case when the isotop is found in DB
00166         // Momentum pMom is in IU ! @@ Units
00167 #ifdef pdebug
00168         G4cout<<"G4QMNCS::G:UpdatDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
00169 #endif
00170         lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update
00171 #ifdef pdebug
00172         G4cout<<"G4QMNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00173         //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00174 #endif
00175         if(lastCS<=0. && pEn>lastTH)    // Correct the threshold
00176         {
00177 #ifdef pdebug
00178           G4cout<<"G4QMNCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00179 #endif
00180           lastTH=pEn;
00181         }
00182         break;                           // Go out of the LOOP
00183       }
00184 #ifdef pdebug
00185       G4cout<<"---G4QMNCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
00186             <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl;
00187       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00188 #endif
00189       j++;                             // Increment a#0f records found in DB for this pPDG
00190     }
00191     if(!in)                            // This nucleus has not been calculated previously
00192     {
00193 #ifdef pdebug
00194       G4cout<<"G4QMNCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
00195 #endif
00197       lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create
00198       if(lastCS<=0.)
00199       {
00200         lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
00201 #ifdef pdebug
00202         G4cout<<"G4QMNCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
00203 #endif
00204         if(pEn>lastTH)
00205         {
00206 #ifdef pdebug
00207           G4cout<<"G4QMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00208 #endif
00209           lastTH=pEn;
00210         }
00211       }
00212 #ifdef pdebug
00213       G4cout<<"G4QMNCS::GetCrosSec: New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00214       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00215 #endif
00216       colN.push_back(tgN);
00217       colZ.push_back(tgZ);
00218       colPDG.push_back(pPDG);
00219       colP.push_back(pMom);
00220       colTH.push_back(lastTH);
00221       colCS.push_back(lastCS);
00222 #ifdef pdebug
00223       G4cout<<"G4QMNCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00224       //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00225 #endif
00226       return lastCS*millibarn;
00227     } // End of creation of the new set of parameters
00228     else
00229     {
00230 #ifdef pdebug
00231       G4cout<<"G4QMNCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
00232 #endif
00233       colP[lastI]=pMom;
00234       colPDG[lastI]=pPDG;
00235       colCS[lastI]=lastCS;
00236     }
00237   } // End of parameters udate
00238   else if(pEn<=lastTH)
00239   {
00240 #ifdef pdebug
00241     G4cout<<"G4QMNCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
00242     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00243 #endif
00244     return 0.;                         // Momentum is below the Threshold Value -> CS=0
00245   }
00246   //  else if(std::fabs(lastP/pMom-1.)<tolerance)
00247   else if(lastP==pMom)                 // VI do not use tolerance
00248   {
00249 #ifdef pdebug
00250     G4cout<<"G4QMNCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", CS="<<lastCS*millibarn<<G4endl;
00251     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00252 #endif
00253     return lastCS*millibarn;     // Use theLastCS
00254   }
00255   else
00256   {
00257 #ifdef pdebug
00258     G4cout<<"G4QMNCS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00259 #endif
00260     lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
00261     lastP=pMom;
00262   }
00263 #ifdef pdebug
00264   G4cout<<"G4QMNCS::GetCroSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00265   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00266 #endif
00267   return lastCS*millibarn;
00268 }

G4double G4QMuonNuclearCrossSection::GetExchangeEnergy (  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 2618 of file G4QMuonNuclearCrossSection.cc.

References G4cerr, G4cout, G4endl, and G4UniformRand.

02619 {
02620   // @@ All constants are copy of that from GetCrossSection funct. => Make them general.
02621   static const G4int nE=336; // !!  If change this, change it in GetFunctions() (*.hh) !!
02622   static const G4int mL=nE-1;
02623   static const G4double EMi=2.0612;          // Minimum Energy
02624   static const G4double EMa=50000.;          // Maximum Energy
02625   static const G4double lEMi=std::log(EMi);       // Minimum logarithmic Energy
02626   static const G4double lEMa=std::log(EMa);       // Maximum logarithmic Energy
02627   static const G4double dlnE=(lEMa-lEMi)/mL; // Logarithmic step in Energy
02628   static const G4double mmu=105.65839;       // Mass of muon in MeV
02629   static const G4double lmmu=std::log(mmu);       // Log of muon mass
02630   G4double phLE=0.;                     // Prototype of the log(nu=E_gamma)
02631   G4double Y[nE];                       // Prepare the array for randomization
02632 #ifdef debug
02633   G4cout<<"G4QMuonNuclCrossSect::GetExchanEn:B="<<lastF<<",l="<<lastL<<",1="<<lastJ1[lastL]
02634         <<",2="<<lastJ2[lastL]<<",3="<<lastJ3[lastL]<<",S="<<lastSig<<",E="<<lastE<<G4endl;
02635 #endif
02636   G4double lastLE=lastG+lmmu;           // recover log(eE) from the gamma (lastG)
02637   G4double dlg1=lastG+lastG-1.;
02638   G4double lgoe=lastG/lastE;
02639   for(G4int i=lastF;i<=lastL;i++)
02640     Y[i]=dlg1*lastJ1[i]-lgoe*(lastJ2[i]+lastJ2[i]-lastJ3[i]/lastE);
02641   G4double ris=lastSig*G4UniformRand(); // Sig can be > Y[lastL=mL], then it is funct. reg.
02642 #ifdef debug
02643   G4cout<<"G4QMuonNuclearCrossSection::GetExchangeEnergy: "<<ris<<",Y="<<Y[lastL]<<G4endl;
02644 #endif
02645   if(ris<Y[lastL])                      // Search in the table
02646   {
02647     G4int j=lastF;
02648     G4double Yj=Y[j];                   // It mast be 0 (some times just very small)
02649     while (ris>Yj && j<lastL)           // Associative search
02650     {
02651       j++;
02652       Yj=Y[j];                          // High value
02653     }
02654     G4int j1=j-1;
02655     G4double Yi=Y[j1];                  // Low value
02656     phLE=lEMi+(j1+(ris-Yi)/(Yj-Yi))*dlnE;
02657 #ifdef debug
02658     G4cout<<"G4QMuNuclearCS::E="<<phLE<<",l="<<lEMi<<",j="<<j<<",ris="<<ris<<",Yi="<<Yi
02659           <<",Y="<<Yj<<G4endl;
02660 #endif
02661   }
02662   else                                  // Search with the function
02663   {
02664     if(lastL<mL)
02665       G4cerr<<"**G4QMuonNucCS::GetExEn:L="<<lastL<<",S="<<lastSig<<",Y="<<Y[lastL]<<G4endl;
02666     G4double f=(ris-Y[lastL])/lastH;    // ScaledResidualValue of the cross-sec. integral
02667 #ifdef pdebug
02668     G4cout<<"G4QMuNucCS::GetExEn:HighEnergy f="<<f<<",ris="<<ris<<",lastH="<<lastH<<G4endl;
02669 #endif
02670     phLE=SolveTheEquation(f);           // Solve equation to find theLog(phE) (comp lastLE)
02671 #ifdef pdebug
02672     G4cout<<"G4QMuonNuclearCrossSection::GetExchangeEnergy:HighEnergy lphE="<<phLE<<G4endl;
02673 #endif
02674   }
02675   if(phLE>lastLE)
02676   {
02677     G4cerr<<"***G4QMuNuclearCS::GetExEnergy: N="<<lastN<<",Z="<<lastZ<<",lpE"<<phLE<<">leE"
02678           <<lastLE<<",Sig="<<lastSig<<",rndSig="<<ris<<",Beg="<<lastF<<",End="<<lastL
02679           <<",Y="<<Y[lastL]<<G4endl;
02680     if(lastLE<7.2) phLE=std::log(std::exp(lastLE)-mmu);
02681     else phLE=7.;
02682   }
02683   return std::exp(phLE);
02684 }

G4int G4QMuonNuclearCrossSection::GetExchangePDGCode (  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 2787 of file G4QMuonNuclearCrossSection.cc.

02787 {return 22;}

G4double G4QMuonNuclearCrossSection::GetExchangeQ2 ( G4double  nu  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 2733 of file G4QMuonNuclearCrossSection.cc.

References G4cerr, and G4UniformRand.

02734 {
02735   static const G4double mmu=105.65839;  // Mass of muon in MeV
02736   static const G4double mmu2=mmu*mmu;   // Squared Mass of muon in MeV
02737   G4double y=nu/lastE;                  // Part of energy carried by the equivalent pfoton
02738   if(y>=1.-1./(lastG+lastG)) return 0.; // The region where the method does not work
02739   G4double y2=y*y;                      // Squared photonic part of energy
02740   G4double ye=1.-y;                     // Part of energy carried by the secondary muon
02741   G4double Qi2=mmu2*y2/ye;              // Minimum Q2
02742   G4double Qa2=4*lastE*lastE*ye;        // Maximum Q2
02743   G4double iar=Qi2/Qa2;                 // Q2min/Q2max ratio
02744   G4double Dy=ye+.5*y2;                 // D(y) function
02745   G4double Py=ye/Dy;                    // P(y) function
02746   G4double ePy=1.-std::exp(Py);              // 1-exp(P(y)) part
02747   G4double Uy=Py*(1.-iar);              // U(y) function
02748   G4double Fy=(ye+ye)*(1.+ye)*iar/y2;   // F(y) function
02749   G4double fr=iar/(1.-ePy*iar);         // Q-fraction
02750   if(Fy<=-fr)
02751   {
02752 #ifdef edebug
02753     G4cerr<<"***G4QMuonNucCS::GetExchQ2: Fy="<<Fy<<"+fr="<<fr<<" <0"<<",iar="<<iar<<G4endl;
02754 #endif
02755     return 0.;
02756   }    
02757   G4double LyQa2=std::log(Fy+fr);            // L(y,Q2max) function
02758   G4bool cond=true;
02759   G4int maxTry=3;
02760   G4int cntTry=0;
02761   G4double Q2=Qi2;
02762   while(cond&&cntTry<maxTry)            // The loop to avoid x>1.
02763   {
02764     G4double R=G4UniformRand();         // Random number (0,1)
02765     Q2=Qi2*(ePy+1./(std::exp(R*LyQa2-(1.-R)*Uy)-Fy));
02766     cntTry++;
02767     cond = Q2>1878.*nu;
02768   }
02769   if(Q2<Qi2)
02770   {
02771 #ifdef edebug
02772     G4cerr<<"***G4QMuonNuclearCrossSect::GetExchangeQ2: Q2="<<Q2<<" < Q2min="<<Qi2<<G4endl;
02773 #endif
02774     return Qi2;
02775   }  
02776   if(Q2>Qa2)
02777   {
02778 #ifdef edebug
02779     G4cerr<<"***G4QMuonNuclearCrossSect::GetExchangeQ2: Q2="<<Q2<<" > Q2max="<<Qi2<<G4endl;
02780 #endif
02781     return Qa2;
02782   }  
02783   return Q2;
02784 }

G4VQCrossSection * G4QMuonNuclearCrossSection::GetPointer (  )  [static]

Definition at line 77 of file G4QMuonNuclearCrossSection.cc.

Referenced by G4QInelastic::GetMeanFreePath(), G4QInelastic::PostStepDoIt(), and G4QAtomicElectronScattering::PostStepDoIt().

00078 {
00079   static G4QMuonNuclearCrossSection theCrossSection; //**Static body of the Cross Section**
00080   return &theCrossSection;
00081 }

G4double G4QMuonNuclearCrossSection::GetVirtualFactor ( G4double  nu,
G4double  Q2 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 2789 of file G4QMuonNuclearCrossSection.cc.

References G4cerr.

02790 {
02791   static const G4double dM=938.27+939.57;// Mean double nucleon mass = m_n+m_p (no binding)
02792   static const G4double Q0=843.;         // Coefficient of the dipole nucleonic form-factor
02793   static const G4double Q02=Q0*Q0;       // Squared coefficient of theDipoleNuclFormFactor
02794   static const G4double blK0=std::log(185.);  // Coefficient of the b-function
02795   static const G4double bp=0.85;          // Power of the b-function
02796   static const G4double clK0=std::log(1390.);  // Coefficient of the c-function
02797   static const G4double cp=3.;            // Power of the c-function
02798   //G4double x=Q2/dM/nu;                  // Direct x definition
02799   G4double K=nu-Q2/dM;                    // K=nu*(1-x)
02800   if(K<0.)
02801   {
02802 #ifdef edebug
02803     G4cerr<<"**G4QMuonNucCS::GetVirtFact:K="<<K<<",nu="<<nu<<",Q2="<<Q2<<",d="<<dM<<G4endl;
02804 #endif
02805     return 0.;
02806   }
02807   G4double lK=std::log(K);                // ln(K)
02808   G4double x=1.-K/nu;                     // This definitin saves one div.
02809   G4double GD=1.+Q2/Q02;                  // Reversed nucleonic form-factor
02810   G4double b=std::exp(bp*(lK-blK0));      // b-factor
02811   G4double c=std::exp(cp*(lK-clK0));      // c-factor
02812   G4double r=.5*std::log(Q2+nu*nu)-lK;    // r=.5*log((Q^2+nu^2)/K^2)
02813   G4double ef=std::exp(r*(b-c*r*r));      // exponential factor
02814   return (1.-x)*ef/GD/GD;
02815 }

G4double G4QMuonNuclearCrossSection::ThresholdEnergy ( G4int  Z,
G4int  N,
G4int  PDG = 13 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 276 of file G4QMuonNuclearCrossSection.cc.

References G4cout, G4endl, G4NucleiProperties::GetNuclearMass(), and G4NucleiProperties::IsInStableTable().

Referenced by GetCrossSection().

00277 {
00278   // CHIPS - Direct GEANT
00279   //static const G4double mNeut = G4QPDGCode(2112).GetMass();
00280   //static const G4double mProt = G4QPDGCode(2212).GetMass();
00281   static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1,0)/MeV;
00282   static const G4double mProt = G4NucleiProperties::GetNuclearMass(1,1)/MeV;
00283   static const G4double mAlph = G4NucleiProperties::GetNuclearMass(4,2)/MeV;
00284   // ---------
00285   static const G4double infEn = 9.e27;
00286 
00287   G4int A=Z+N;
00288   if(A<1) return infEn;
00289   else if(A==1) return 160.; // Pi0 threshold in MeV for the proton:T>m+(m^2+2lm)/2M,l=m_mu
00290   // CHIPS - Direct GEANT
00291   //G4double mT= G4QPDGCode(111).GetNuclMass(Z,N,0);
00292   G4double mT= 0.;
00293   if(Z&&G4NucleiProperties::IsInStableTable(A,Z))
00294                             mT = G4NucleiProperties::GetNuclearMass(A,Z)/MeV;
00295   else return 0.;       // If it is not in the Table of Stable Nuclei, then the Threshold=0
00296   // --------- Splitting thresholds
00297   G4double mP= infEn;
00298   if(A>1&&Z&&G4NucleiProperties::IsInStableTable(A-1,Z-1))
00299           mP = G4NucleiProperties::GetNuclearMass(A-1,Z-1)/MeV;// ResNucMass for a proton
00300 
00301   G4double mN= infEn;
00302   if(A>1&&G4NucleiProperties::IsInStableTable(A-1,Z))
00303           mN = G4NucleiProperties::GetNuclearMass(A-1,Z)/MeV;  // ResNucMass for a neutron
00304 
00305   G4double mA= infEn;
00306   if(A>4&&Z>1&&G4NucleiProperties::IsInStableTable(A-4,Z-2))
00307           mA=G4NucleiProperties::GetNuclearMass(A-4,Z-2)/MeV;  // ResNucMass for an alpha
00308 
00309   G4double dP= mP +mProt - mT;
00310   G4double dN= mN +mNeut - mT;
00311   G4double dA= mA +mAlph - mT;
00312 #ifdef pdebug
00313   G4cout<<"G4QMuonNucCS::ThreshEn: mP="<<mP<<",dP="<<dP<<",mN="<<mN<<",dN="<<dN<<",mA="
00314         <<mA<<",dA="<<dA<<",mT="<<mT<<",A="<<A<<",Z="<<Z<<G4endl;
00315 #endif
00316   if(dP<dN)dN=dP;
00317   if(dA<dN)dN=dA;
00318   return dN;
00319 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:11 2013 for Geant4 by  doxygen 1.4.7