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