#include <G4QTauNuclearCrossSection.hh>
Inheritance diagram for G4QTauNuclearCrossSection:
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 G4VQCrossSection * | GetPointer () |
Protected Member Functions | |
G4QTauNuclearCrossSection () |
Definition at line 57 of file G4QTauNuclearCrossSection.hh.
G4QTauNuclearCrossSection::G4QTauNuclearCrossSection | ( | ) | [inline, protected] |
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 }
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 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 }
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 }
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 }