G4QTauNuclearCrossSection Class Reference

#include <G4QTauNuclearCrossSection.hh>

Inheritance diagram for G4QTauNuclearCrossSection:

G4VQCrossSection

Public Member Functions

 ~G4QTauNuclearCrossSection ()
G4double ThresholdEnergy (G4int Z, G4int N, G4int PDG=15)
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

 G4QTauNuclearCrossSection ()

Detailed Description

Definition at line 57 of file G4QTauNuclearCrossSection.hh.


Constructor & Destructor Documentation

G4QTauNuclearCrossSection::G4QTauNuclearCrossSection (  )  [inline, protected]

Definition at line 61 of file G4QTauNuclearCrossSection.hh.

00061 {};

G4QTauNuclearCrossSection::~G4QTauNuclearCrossSection (  ) 

Definition at line 84 of file G4QTauNuclearCrossSection.cc.

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


Member Function Documentation

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

Implements G4VQCrossSection.

Definition at line 321 of file G4QTauNuclearCrossSection.cc.

References G4cerr, G4cout, and G4endl.

Referenced by GetCrossSection().

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

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

Reimplemented from G4VQCrossSection.

Definition at line 99 of file G4QTauNuclearCrossSection.cc.

References CalculateCrossSection(), G4cout, G4endl, ThresholdEnergy(), and G4VQCrossSection::tolerance.

00101 {
00102   static const G4double mtu=1777.;     // Mass of a tau lepton in MeV
00103   static const G4double mtu2=mtu*mtu;  // Squared Mass of a tau-lepton in MeV^2
00104   static G4int j;                      // A#0f records found in DB for this projectile
00105   static std::vector <G4int>    colPDG;// Vector of the projectile PDG code
00106   static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
00107   static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
00108   static std::vector <G4double> colP;  // Vector of last momenta for the reaction
00109   static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
00110   static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
00111   // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
00112   G4double pEn=std::sqrt(pMom*pMom+mtu2)-mtu; // ==> tau-/tau+ kinEnergy
00113 #ifdef pdebug
00114   G4cout<<"G4QTNCS::GetCS:>>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
00115         <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="
00116         <<colN.size()<<G4endl;
00117   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00118 #endif
00119   if(std::abs(pPDG)!=15)
00120   {
00121 #ifdef pdebug
00122     G4cout<<"G4QTNCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00123     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00124 #endif
00125     return 0.;                         // projectile PDG=0 is a mistake (?!) @@
00126   }
00127   G4bool in=false;                     // By default the isotope must be found in the AMDB
00128   if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
00129   {
00130     in = false;                        // By default the isotope haven't be found in AMDB  
00131     lastP   = 0.;                      // New momentum history (nothing to compare with)
00132     lastPDG = pPDG;                    // The last PDG of the projectile
00133     lastN   = tgN;                     // The last N of the calculated nucleus
00134     lastZ   = tgZ;                     // The last Z of the calculated nucleus
00135     lastI   = colN.size();             // Size of the Associative Memory DB in the heap
00136     j  = 0;                            // A#0f records found in DB for this projectile
00137     if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found
00138     {                                  // The nucleus with projPDG is found in AMDB
00139       if(colN[i]==tgN && colZ[i]==tgZ)
00140       {
00141         lastI=i;
00142         lastTH =colTH[i];                // Last THreshold (A-dependent)
00143 #ifdef pdebug
00144         G4cout<<"G4QTNCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
00145         //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00146 #endif
00147         if(pEn<=lastTH)
00148         {
00149 #ifdef pdebug
00150           G4cout<<"G4QTNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",CS=0"<<G4endl;
00151           //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00152 #endif
00153           return 0.;                     // Energy is below the Threshold value
00154         }
00155         lastP  =colP [i];                // Last Momentum  (A-dependent)
00156         lastCS =colCS[i];                // Last CrossSect (A-dependent)
00157         if(std::fabs(lastP/pMom-1.)<tolerance)
00158         {
00159 #ifdef pdebug
00160           G4cout<<"G4QTNCS::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<<"G4QTNCS::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<<"G4QTNCS::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<<"G4QTNCS::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<<"---G4QTNCrossSec::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<<"G4QTNCS::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<<"G4QTNCrossSection::GetCrossSect: NewThresh="<<lastTH<<", T="<<pEn<<G4endl;
00203 #endif
00204         if(pEn>lastTH)
00205         {
00206 #ifdef pdebug
00207           G4cout<<"G4QTNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00208 #endif
00209           lastTH=pEn;
00210         }
00211       }
00212 #ifdef pdebug
00213       G4cout<<"G4QTNCS::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<<"G4QTNCS::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<<"G4QTNCS::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<<"G4QTNCS::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   {
00248 #ifdef pdebug
00249     G4cout<<"G4QTNCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<", CS="<<lastCS*millibarn<<G4endl;
00250     //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00251 #endif
00252     return lastCS*millibarn;     // Use theLastCS
00253   }
00254   else
00255   {
00256 #ifdef pdebug
00257     G4cout<<"G4QTNCS::GetCS:UpdatCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
00258 #endif
00259     lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
00260     lastP=pMom;
00261   }
00262 #ifdef pdebug
00263   G4cout<<"G4QTNCS::GetCroSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00264   //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
00265 #endif
00266   return lastCS*millibarn;
00267 }

G4double G4QTauNuclearCrossSection::GetExchangeEnergy (  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 2622 of file G4QTauNuclearCrossSection.cc.

References G4cerr, G4cout, G4endl, and G4UniformRand.

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

G4int G4QTauNuclearCrossSection::GetExchangePDGCode (  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 2791 of file G4QTauNuclearCrossSection.cc.

02791 {return 22;}

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

Reimplemented from G4VQCrossSection.

Definition at line 2737 of file G4QTauNuclearCrossSection.cc.

References G4cerr, and G4UniformRand.

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

G4VQCrossSection * G4QTauNuclearCrossSection::GetPointer (  )  [static]

Definition at line 78 of file G4QTauNuclearCrossSection.cc.

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

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

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

Reimplemented from G4VQCrossSection.

Definition at line 2793 of file G4QTauNuclearCrossSection.cc.

References G4cerr.

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

G4double G4QTauNuclearCrossSection::ThresholdEnergy ( G4int  Z,
G4int  N,
G4int  PDG = 15 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 275 of file G4QTauNuclearCrossSection.cc.

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

Referenced by GetCrossSection().

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


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