Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Protected Member Functions | List of all members
G4QNuMuNuclearCrossSection Class Reference

#include <G4QNuMuNuclearCrossSection.hh>

+ Inheritance diagram for G4QNuMuNuclearCrossSection:

Public Member Functions

 ~G4QNuMuNuclearCrossSection ()
 
G4double ThresholdEnergy (G4int Z, G4int N, G4int PDG=14)
 
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 GetDirectPart (G4double Q2)
 
G4double GetNPartons (G4double Q2)
 
G4double GetQEL_ExchangeQ2 ()
 
G4double GetNQE_ExchangeQ2 ()
 
G4double GetLastTOTCS ()
 
G4double GetLastQELCS ()
 
- Public Member Functions inherited from G4VQCrossSection
virtual ~G4VQCrossSection ()
 
virtual G4double GetCrossSection (G4bool, G4double, G4int, G4int, G4int pPDG=0)
 
virtual G4double ThresholdEnergy (G4int Z, G4int N, G4int PDG=0)
 
virtual G4double CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int PDG, G4int tgZ, G4int tgN, G4double pMom)=0
 
virtual G4double GetLastTOTCS ()
 
virtual G4double GetLastQELCS ()
 
virtual G4double GetDirectPart (G4double Q2)
 
virtual G4double GetNPartons (G4double Q2)
 
virtual G4double GetExchangeEnergy ()
 
virtual G4double GetExchangeT (G4int tZ, G4int tN, G4int pPDG)
 
virtual G4double GetSlope (G4int tZ, G4int tN, G4int pPDG)
 
virtual G4double GetHMaxT ()
 
virtual G4double GetExchangeQ2 (G4double nu=0)
 
virtual G4double GetVirtualFactor (G4double nu, G4double Q2)
 
virtual G4double GetQEL_ExchangeQ2 ()
 
virtual G4double GetNQE_ExchangeQ2 ()
 
virtual G4int GetExchangePDGCode ()
 

Static Public Member Functions

static G4VQCrossSection * GetPointer ()
 
- Static Public Member Functions inherited from G4VQCrossSection
static void setTolerance (G4double tol)
 

Protected Member Functions

 G4QNuMuNuclearCrossSection ()
 
- Protected Member Functions inherited from G4VQCrossSection
 G4VQCrossSection ()
 
G4double LinearFit (G4double X, G4int N, G4double *XN, G4double *YN)
 
G4double EquLinearFit (G4double X, G4int N, G4double X0, G4double DX, G4double *Y)
 

Additional Inherited Members

- Static Protected Attributes inherited from G4VQCrossSection
static G4double tolerance =.001
 

Detailed Description

Definition at line 49 of file G4QNuMuNuclearCrossSection.hh.

Constructor & Destructor Documentation

◆ G4QNuMuNuclearCrossSection()

G4QNuMuNuclearCrossSection::G4QNuMuNuclearCrossSection ( )
inlineprotected

Definition at line 53 of file G4QNuMuNuclearCrossSection.hh.

53{};

◆ ~G4QNuMuNuclearCrossSection()

G4QNuMuNuclearCrossSection::~G4QNuMuNuclearCrossSection ( )

Definition at line 78 of file G4QNuMuNuclearCrossSection.cc.

79{
80 G4int lens=TX->size();
81 for(G4int i=0; i<lens; ++i) delete[] (*TX)[i];
82 delete TX;
83 G4int hens=QE->size();
84 for(G4int i=0; i<hens; ++i) delete[] (*QE)[i];
85 delete QE;
86}
G4int
int G4int
Definition: G4Types.hh:66

Member Function Documentation

◆ CalculateCrossSection()

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

Implements G4VQCrossSection.

Definition at line 278 of file G4QNuMuNuclearCrossSection.cc.

280{
281 static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2
282 static const G4int nE=65; // !! If change this, change it in GetFunctions() (*.hh) !!
283 static const G4int mL=nE-1;
284 static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
285 static const G4double dmN=mN+mN; // Doubled nucleon mass (2*AtomicMassUnit, GeV)
286 static const G4double mmu=.105658369; // Mass of a muon in GeV
287 static const G4double mmu2=mmu*mmu;// Squared mass of a muon in GeV^2
288 static const G4double EMi=mmu+mmu2/dmN; // Universal threshold of the reaction in GeV
289 static const G4double EMa=300.; // Maximum tabulated Energy of nu_mu in GeV
290 // *** Begin of the Associative memory for acceleration of the cross section calculations
291 static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional)
292 static G4bool first=true; // Flag of initialization of the energy axis
293 // *** End of Static Definitions (Associative Memory) ***
294 //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
295 //G4double TotEnergy2=Momentum;
296 onlyCS=CS; // Flag to calculate only CS (not TX & QE)
297 lastE=Momentum/GeV; // Kinetic energy of the muon neutrino (in GeV!)
298 if (lastE<=EMi) // Energy is below the minimum energy in the table
299 {
300 lastE=0.;
301 lastSig=0.;
302 return 0.;
303 }
304 G4int Z=targZ; // New Z, which can change the sign
305 if(F<=0) // This isotope was not the last used isotop
306 {
307 if(F<0) // This isotope was found in DAMDB =-------=> RETRIEVE
308 {
309 lastTX =(*TX)[I]; // Pointer to the prepared TX function (same isotope)
310 lastQE =(*QE)[I]; // Pointer to the prepared QE function (same isotope)
311 }
312 else // This isotope wasn't calculated previously => CREATE
313 {
314 if(first)
315 {
316 lastEN = new G4double[nE]; // This must be done only once!
317 Z=-Z; // To explain GetFunctions that E-axis must be filled
318 first=false; // To make it only once
319 }
320 lastTX = new G4double[nE]; // Allocate memory for the new TX function
321 lastQE = new G4double[nE]; // Allocate memory for the new QE function
322 G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK)
323 if(res<0) G4cerr<<"*W*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
324 // *** The synchronization check ***
325 G4int sync=TX->size();
326 if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
327 TX->push_back(lastTX);
328 QE->push_back(lastQE);
329 } // End of creation of the new set of parameters
330 } // End of parameters udate
331 // =--------------------= NOW Calculate the Cross Section =-------------------=
332 if (lastE<=EMi) // Check that the neutrinoEnergy is higher than ThreshE
333 {
334 lastE=0.;
335 lastSig=0.;
336 return 0.;
337 }
338 if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization
339 {
340 G4int chk=1;
341 G4int ran=mL/2;
342 G4int sep=ran; // as a result = an index of the left edge of the interval
343 while(ran>=2)
344 {
345 G4int newran=ran/2;
346 if(lastE<=lastEN[sep]) sep-=newran;
347 else sep+=newran;
348 ran=newran;
349 chk=chk+chk;
350 }
351 if(chk+chk!=mL) G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
352 G4double lowE=lastEN[sep];
353 G4double highE=lastEN[sep+1];
354 G4double lowTX=lastTX[sep];
355 if(lastE<lowE||sep>=mL||lastE>highE)
356 G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
357 <<", sep="<<sep<<", mL="<<mL<<G4endl;
358 lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
359 if(!onlyCS) // Skip the differential cross-section parameters
360 {
361 G4double lowQE=lastQE[sep];
362 lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
363#ifdef pdebug
364 G4cout<<"G4NuMuNuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
365#endif
366 }
367 }
368 else
369 {
370 lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
371 lastQEL=lastQE[mL];
372 }
373 if(lastQEL<0.) lastQEL = 0.;
374 if(lastSig<0.) lastSig = 0.;
375 // The cross-sections are expected to be in mb
376 lastSig*=mb38;
377 if(!onlyCS) lastQEL*=mb38;
378 return lastSig;
379}
G4double
double G4double
Definition: G4Types.hh:64
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4cout
G4DLLIMPORT std::ostream G4cout

Referenced by GetCrossSection().

◆ GetCrossSection()

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

!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)

Reimplemented from G4VQCrossSection.

Definition at line 90 of file G4QNuMuNuclearCrossSection.cc.

92{
93 static G4int j; // A#0f records found in DB for this projectile
94 static std::vector <G4int> colPDG;// Vector of the projectile PDG code
95 static std::vector <G4int> colN; // Vector of N for calculated nuclei (isotops)
96 static std::vector <G4int> colZ; // Vector of Z for calculated nuclei (isotops)
97 static std::vector <G4double> colP; // Vector of last momenta for the reaction
98 static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
99 static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
100 // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
101 G4double pEn=pMom;
102#ifdef debug
103 G4cout<<"G4QNMNCS::GetCS:>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
104 <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="
105 <<colN.size()<<G4endl;
106 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
107#endif
108 if(pPDG!=14)
109 {
110#ifdef pdebug
111 G4cout<<"G4QNMNCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
112 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
113#endif
114 return 0.; // projectile PDG=0 is a mistake (?!) @@
115 }
116 G4bool in=false; // By default the isotope must be found in the AMDB
117 if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
118 {
119 in = false; // By default the isotope haven't be found in AMDB
120 lastP = 0.; // New momentum history (nothing to compare with)
121 lastPDG = pPDG; // The last PDG of the projectile
122 lastN = tgN; // The last N of the calculated nucleus
123 lastZ = tgZ; // The last Z of the calculated nucleus
124 lastI = colN.size(); // Size of the Associative Memory DB in the heap
125 j = 0; // A#0f records found in DB for this projectile
126 if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found
127 { // The nucleus with projPDG is found in AMDB
128 if(colN[i]==tgN && colZ[i]==tgZ)
129 {
130 lastI=i;
131 lastTH =colTH[i]; // Last THreshold (A-dependent)
132#ifdef pdebug
133 G4cout<<"G4QNMNCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
134 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
135#endif
136 if(pEn<=lastTH)
137 {
138#ifdef pdebug
139 G4cout<<"G4QNMNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",X=0"<<G4endl;
140 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
141#endif
142 return 0.; // Energy is below the Threshold value
143 }
144 lastP =colP [i]; // Last Momentum (A-dependent)
145 lastCS =colCS[i]; // Last CrossSect (A-dependent)
146 if(std::fabs(lastP/pMom-1.)<tolerance)
147 {
148#ifdef pdebug
149 G4cout<<"G4QNMNCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
150 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
151#endif
152 return lastCS*millibarn; // Use theLastCS
153 }
154 in = true; // This is the case when the isotop is found in DB
155 // Momentum pMom is in IU ! @@ Units
156#ifdef pdebug
157 G4cout<<"G4QNMNCS::G:UpdaDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
158#endif
159 lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update
160#ifdef pdebug
161 G4cout<<"G4QNMNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
162 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
163#endif
164 if(lastCS<=0. && pEn>lastTH) // Correct the threshold
165 {
166#ifdef pdebug
167 G4cout<<"G4QNMNCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
168#endif
169 lastTH=pEn;
170 }
171 break; // Go out of the LOOP
172 }
173#ifdef pdebug
174 G4cout<<"---G4QNMNCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
175 <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl;
176 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
177#endif
178 j++; // Increment a#0f records found in DB for this pPDG
179 }
180 if(!in) // This nucleus has not been calculated previously
181 {
182#ifdef pdebug
183 G4cout<<"G4QNMNCS::GetCrSec: CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
184#endif
185 //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
186 lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create
187 if(lastCS<=0.)
188 {
189 lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
190#ifdef pdebug
191 G4cout<<"G4QNMNCrossSection::GetCrossSect:NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
192#endif
193 if(pEn>lastTH)
194 {
195#ifdef pdebug
196 G4cout<<"G4QNMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
197#endif
198 lastTH=pEn;
199 }
200 }
201#ifdef pdebug
202 G4cout<<"G4QNMNCS::GetCrosSec:New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
203 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
204#endif
205 colN.push_back(tgN);
206 colZ.push_back(tgZ);
207 colPDG.push_back(pPDG);
208 colP.push_back(pMom);
209 colTH.push_back(lastTH);
210 colCS.push_back(lastCS);
211#ifdef pdebug
212 G4cout<<"G4QNMNCS::GetCS:1st,P="<<pMom<<"(MeV),X="<<lastCS*millibarn<<"(mb)"<<G4endl;
213 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
214#endif
215 return lastCS*millibarn;
216 } // End of creation of the new set of parameters
217 else
218 {
219#ifdef pdebug
220 G4cout<<"G4QNMNCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
221#endif
222 colP[lastI]=pMom;
223 colPDG[lastI]=pPDG;
224 colCS[lastI]=lastCS;
225 }
226 } // End of parameters udate
227 else if(pEn<=lastTH)
228 {
229#ifdef pdebug
230 G4cout<<"G4QNMNCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
231 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
232#endif
233 return 0.; // Momentum is below the Threshold Value -> CS=0
234 }
235 else if(std::fabs(lastP/pMom-1.)<tolerance)
236 {
237#ifdef pdebug
238 G4cout<<"G4QNMNCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
239 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
240#endif
241 return lastCS*millibarn; // Use theLastCS
242 }
243 else
244 {
245#ifdef pdebug
246 G4cout<<"G4QNMNCS::GetCS:UpdaCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
247#endif
248 lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
249 lastP=pMom;
250 }
251#ifdef pdebug
252 G4cout<<"G4QNMNCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
253 //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
254#endif
255 return lastCS*millibarn;
256}
G4QNuMuNuclearCrossSection::ThresholdEnergy
G4double ThresholdEnergy(G4int Z, G4int N, G4int PDG=14)
Definition: G4QNuMuNuclearCrossSection.cc:259
G4QNuMuNuclearCrossSection::CalculateCrossSection
G4double CalculateCrossSection(G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum)
Definition: G4QNuMuNuclearCrossSection.cc:278
G4VQCrossSection::tolerance
static G4double tolerance
Definition: G4VQCrossSection.hh:175

◆ GetDirectPart()

G4double G4QNuMuNuclearCrossSection::GetDirectPart ( G4double  Q2)
virtual

Reimplemented from G4VQCrossSection.

Definition at line 763 of file G4QNuMuNuclearCrossSection.cc.

764{
765 G4double f=Q2/4.62;
766 G4double ff=f*f;
767 G4double r=ff*ff;
768 G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
769 //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
770 return 1.-s_value*(1.-s_value/2);
771}

◆ GetExchangePDGCode()

G4int G4QNuMuNuclearCrossSection::GetExchangePDGCode ( )
virtual

Reimplemented from G4VQCrossSection.

Definition at line 780 of file G4QNuMuNuclearCrossSection.cc.

780{return 211;}

◆ GetLastQELCS()

G4double G4QNuMuNuclearCrossSection::GetLastQELCS ( )
inlinevirtual

Reimplemented from G4VQCrossSection.

Definition at line 82 of file G4QNuMuNuclearCrossSection.hh.

82{return lastQEL;}

◆ GetLastTOTCS()

G4double G4QNuMuNuclearCrossSection::GetLastTOTCS ( )
inlinevirtual

Reimplemented from G4VQCrossSection.

Definition at line 81 of file G4QNuMuNuclearCrossSection.hh.

81{return lastSig;}

◆ GetNPartons()

G4double G4QNuMuNuclearCrossSection::GetNPartons ( G4double  Q2)
virtual

Reimplemented from G4VQCrossSection.

Definition at line 774 of file G4QNuMuNuclearCrossSection.cc.

775{
776 return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
777}

◆ GetNQE_ExchangeQ2()

G4double G4QNuMuNuclearCrossSection::GetNQE_ExchangeQ2 ( )
virtual

Reimplemented from G4VQCrossSection.

Definition at line 531 of file G4QNuMuNuclearCrossSection.cc.

532{
533 static const double mpi=.13957018; // charged pi meson mass in GeV
534 static const G4double mmu=.105658369; // Mass of muon in GeV
535 static const G4double mmu2=mmu*mmu; // Squared Mass of muon in GeV^2
536 static const double hmmu2=mmu2/2; // .5*m_mu^2 in GeV^2
537 static const double MN=.931494043; // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
538 static const double MN2=MN*MN; // M_N^2 in GeV^2
539 static const double dMN=MN+MN; // 2*M_N in GeV
540 static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
541 static const G4int power=7; // direct power for the magic variable
542 static const G4double pconv=1./power; // conversion power for the magic variable
543 static const G4int nX=21; // #Of point in the Xl table (in GeV^2)
544 static const G4int lX=nX-1; // index of the last in the Xl table
545 static const G4int bX=lX-1; // @@ index of the before last in the Xl table
546 static const G4int nE=20; // #Of point in the El table (in GeV^2)
547 static const G4int bE=nE-1; // index of the last in the El table
548 static const G4int pE=bE-1; // index of the before last in the El table
549 // Reversed table
550 static const G4double X0[nX]={6.14081e-05,
551 .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244,
552 2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419};
553 static const G4double X1[nX]={.00125268,
554 .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049,
555 4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330};
556 static const G4double X2[nX]={.015694,
557 1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728,
558 10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276};
559 static const G4double X3[nX]={.0866877,
560 4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302,
561 21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115};
562 static const G4double X4[nX]={.160483,
563 5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127,
564 29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265};
565 static const G4double X5[nX]={.0999307,
566 5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398,
567 26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649};
568 static const G4double X6[nX]={.0276367,
569 3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519,
570 16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148};
571 static const G4double X7[nX]={.00472383,
572 2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473,
573 9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057};
574 static const G4double X8[nX]={.000630783,
575 1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480,
576 5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037};
577 static const G4double X9[nX]={7.49179e-05,
578 .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823,
579 3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028};
580 static const G4double XA[nX]={8.43437e-06,
581 .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704,
582 2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823};
583 static const G4double XB[nX]={9.27028e-07,
584 .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059,
585 1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341};
586 static const G4double XC[nX]={1.00807e-07,
587 .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668,
588 1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014};
589 static const G4double XD[nX]={1.09102e-08,
590 .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878,
591 1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617};
592 static const G4double XE[nX]={1.17831e-09,
593 .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374,
594 .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571};
595 static const G4double XF[nX]={1.27141e-10,
596 .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046,
597 .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338};
598 static const G4double XG[nX]={1.3713e-11,
599 .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870,
600 .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379};
601 static const G4double XH[nX]={1.47877e-12,
602 .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421,
603 .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197};
604 static const G4double XI[nX]={1.59454e-13,
605 .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806,
606 .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872};
607 static const G4double XJ[nX]={1.71931e-14,
608 .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147,
609 .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813};
610 // Direct table
611 static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
612 X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
613 static const G4double dX[nE]={
614 (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
615 (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
616 (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
617 (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
618 (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
619 static const G4double* Xl[nE]=
620 {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
621 static const G4double I0[nX]={0,
622 .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555,
623 13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000};
624 static const G4double I1[nX]={0,
625 .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365,
626 13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000};
627 static const G4double I2[nX]={0,
628 .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586,
629 12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000};
630 static const G4double I3[nX]={0,
631 .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461,
632 12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000};
633 static const G4double I4[nX]={0,
634 .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515,
635 11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000};
636 static const G4double I5[nX]={0,
637 .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040,
638 10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000};
639 static const G4double I6[nX]={0,
640 .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243,
641 9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000};
642 static const G4double I7[nX]={0,
643 .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717,
644 8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000};
645 static const G4double I8[nX]={0,
646 .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415,
647 8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000};
648 static const G4double I9[nX]={0,
649 .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672,
650 8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000};
651 static const G4double IA[nX]={0,
652 .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265,
653 7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000};
654 static const G4double IB[nX]={0,
655 .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298,
656 7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000};
657 static const G4double IC[nX]={0,
658 .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696,
659 7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000};
660 static const G4double ID[nX]={0,
661 .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212,
662 7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000};
663 static const G4double IE[nX]={0,
664 .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949,
665 7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000};
666 static const G4double IF[nX]={0,
667 .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164,
668 7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000};
669 static const G4double IG[nX]={0,
670 .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249,
671 7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000};
672 static const G4double IH[nX]={0,
673 .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804,
674 7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000};
675 static const G4double II[nX]={0,
676 .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722,
677 7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000};
678 static const G4double IJ[nX]={0,
679 .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196,
680 6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 20.0000};
681 static const G4double* Il[nE]=
682 {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
683 static const G4double lE[nE]={
684-1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
685 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
686 static const G4double lEmi=lE[0];
687 static const G4double lEma=lE[nE-1];
688 static const G4double dlE=(lEma-lEmi)/bE;
689 //***************************************************************************************
690 G4double Enu=lastE; // Get energy of the last calculated cross-section
691 G4double lEn=std::log(Enu); // log(E) for interpolation
692 G4double rE=(lEn-lEmi)/dlE; // Position of the energy
693 G4int fE=static_cast<int>(rE); // Left bin for interpolation
694 if(fE<0) fE=0;
695 if(fE>pE)fE=pE;
696 G4int sE=fE+1; // Right bin for interpolation
697 G4double dE=rE-fE; // relative log shift from the left bin
698 G4double dEnu=Enu+Enu; // doubled energy of nu/anu
699 G4double Enu2=Enu*Enu; // squared energy of nu/anu
700 G4double Emu=Enu-mmu; // Free Energy of neutrino/anti-neutrino
701 G4double ME=Enu*MN; // M*E
702 G4double dME=ME+ME; // 2*M*E
703 G4double dEMN=(dEnu+MN)*ME;
704 G4double MEm=ME-hmmu2;
705 G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
706 G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
707 G4double ymax=(E2M+sqE)/dEMN;
708 G4double ymin=(E2M-sqE)/dEMN;
709 G4double rmin=1.-ymin;
710 G4double rhm2E=hmmu2/Enu2;
711 G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
712 G4double Q2ma=dME*ymax; // Q2_max(E_nu)
713 G4double Q2nq=Emu*dMN-mcV;
714 if(Q2ma>Q2nq) Q2ma=Q2nq; // Correction for Non Quasi Elastic
715 // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
716 G4double Rmi=Q2mi/Q2ma;
717 G4double shift=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu
718 // --- E-interpolation must be done in a log scale ---
719 G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
720 G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
721 // Find the integral values integ(Xmi) & integ(Xma) using the direct table
722 G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
723 G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
724 G4double rXi=(Xmi-iXmi)/idX;
725 G4int iXi=static_cast<int>(rXi);
726 if(iXi<0) iXi=0;
727 if(iXi>bX) iXi=bX;
728 G4double dXi=rXi-iXi;
729 G4double bntil=Il[fE][iXi];
730 G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
731 G4double bntir=Il[sE][iXi];
732 G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
733 G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
734 //
735 G4double rXa=(Xma-iXmi)/idX;
736 G4int iXa=static_cast<int>(rXa);
737 if(iXa<0) iXa=0;
738 if(iXa>bX) iXa=bX;
739 G4double dXa=rXa-iXa;
740 G4double bntal=Il[fE][iXa];
741 G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
742 G4double bntar=Il[sE][iXa];
743 G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
744 G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
745 //
746 // *** Find X using the reversed table ***
747 G4double intx=inti+(inta-inti)*G4UniformRand();
748 G4int intc=static_cast<int>(intx);
749 if(intc<0) intc=0;
750 if(intc>bX) intc=bX;
751 G4double dint=intx-intc;
752 G4double mXl=Xl[fE][intc];
753 G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
754 G4double mXr=Xl[sE][intc];
755 G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
756 G4double X=Xlb+dE*(Xrb-Xlb); // interpolated X value
757 G4double R=shift-std::pow(X,pconv);
758 G4double Q2=R*Q2ma;
759 return Q2*GeV*GeV;
760}
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53

◆ GetPointer()

G4VQCrossSection * G4QNuMuNuclearCrossSection::GetPointer ( )
static

Definition at line 72 of file G4QNuMuNuclearCrossSection.cc.

73{
74 static G4QNuMuNuclearCrossSection theCrossSection; //**Static body of the Cross Section**
75 return &theCrossSection;
76}

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

◆ GetQEL_ExchangeQ2()

G4double G4QNuMuNuclearCrossSection::GetQEL_ExchangeQ2 ( )
virtual

Reimplemented from G4VQCrossSection.

Definition at line 446 of file G4QNuMuNuclearCrossSection.cc.

447{
448 static const G4double mmu=.105658369;// Mass of muon in GeV
449 static const G4double mmu2=mmu*mmu; // Squared Mass of muon in GeV^2
450 static const double hmmu2=mmu2/2; // .5*m_mu^2 in GeV^2
451 static const double MN=.931494043; // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
452 static const double MN2=MN*MN; // M_N^2 in GeV^2
453 static const G4double power=-3.5; // direct power for the magic variable
454 static const G4double pconv=1./power;// conversion power for the magic variable
455 static const G4int nQ2=101; // #Of point in the Q2l table (in GeV^2)
456 static const G4int lQ2=nQ2-1; // index of the last in the Q2l table
457 static const G4int bQ2=lQ2-1; // index of the before last in the Q2 ltable
458 // Reversed table
459 static const G4double Xl[nQ2]={1.87905e-10,
460 .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374,
461 .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069,
462 .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718,
463 .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169,
464 .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670,
465 .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720,
466 .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457,
467 .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366,
468 .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128,
469 .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539};
470 // Direct table
471 static const G4double Xmax=Xl[lQ2];
472 static const G4double Xmin=Xl[0];
473 static const G4double dX=(Xmax-Xmin)/lQ2; // step in X(Q2, GeV^2)
474 static const G4double inl[nQ2]={0,
475 1.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811,
476 15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039,
477 26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876,
478 36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950,
479 46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377,
480 55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162,
481 65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852,
482 74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767,
483 83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106,
484 92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 100.000};
485 G4double Enu=lastE; // Get energy of the last calculated cross-section
486 G4double dEnu=Enu+Enu; // doubled energy of nu/anu
487 G4double Enu2=Enu*Enu; // squared energy of nu/anu
488 G4double ME=Enu*MN; // M*E
489 G4double dME=ME+ME; // 2*M*E
490 G4double dEMN=(dEnu+MN)*ME;
491 G4double MEm=ME-hmmu2;
492 G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
493 G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
494 G4double ymax=(E2M+sqE)/dEMN;
495 G4double ymin=(E2M-sqE)/dEMN;
496 G4double rmin=1.-ymin;
497 G4double rhm2E=hmmu2/Enu2;
498 G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
499 G4double Q2ma=dME*ymax; // Q2_max(E_nu)
500 G4double Xma=std::pow((1.+Q2mi),power); // X_max(E_nu)
501 G4double Xmi=std::pow((1.+Q2ma),power); // X_min(E_nu)
502 // Find the integral values integ(Xmi) & integ(Xma) using the direct table
503 G4double rXi=(Xmi-Xmin)/dX;
504 G4int iXi=static_cast<int>(rXi);
505 if(iXi<0) iXi=0;
506 if(iXi>bQ2) iXi=bQ2;
507 G4double dXi=rXi-iXi;
508 G4double bnti=inl[iXi];
509 G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
510 //
511 G4double rXa=(Xma-Xmin)/dX;
512 G4int iXa=static_cast<int>(rXa);
513 if(iXa<0) iXa=0;
514 if(iXa>bQ2) iXa=bQ2;
515 G4double dXa=rXa-iXa;
516 G4double bnta=inl[iXa];
517 G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
518 // *** Find X using the reversed table ***
519 G4double intx=inti+(inta-inti)*G4UniformRand();
520 G4int intc=static_cast<int>(intx);
521 if(intc<0) intc=0;
522 if(intc>bQ2) intc=bQ2; // If it is more than max, then the BAD extrapolation
523 G4double dint=intx-intc;
524 G4double mX=Xl[intc];
525 G4double X=mX+dint*(Xl[intc+1]-mX);
526 G4double Q2=std::pow(X,pconv)-1.;
527 return Q2*GeV*GeV;
528}

◆ ThresholdEnergy()

G4double G4QNuMuNuclearCrossSection::ThresholdEnergy ( G4int  Z,
G4int  N,
G4int  PDG = 14 
)
virtual

Reimplemented from G4VQCrossSection.

Definition at line 259 of file G4QNuMuNuclearCrossSection.cc.

260{
261 //static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1,0)/GeV;
262 //static const G4double mProt = G4NucleiProperties::GetNuclearMass(1,1)/GeV;
263 //static const G4double mDeut = G4NucleiProperties::GetNuclearMass(2,1)/GeV/2.;
264 static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
265 static const G4double dmN=mN+mN; // Doubled nucleon mass (2*AtomicMassUnit, GeV)
266 static const G4double mmu=.105658369; // Mass of a muon in GeV
267 static const G4double mmu2=mmu*mmu; // Squared mass of a muon in GeV^2
268 static const G4double thresh=mmu+mmu2/dmN; // Universal threshold in GeV
269 // ---------
270 //static const G4double infEn = 9.e27;
271 G4double dN=0.;
272 if(Z>0||N>0) dN=thresh*GeV; // @@ if upgraded, change it in a total cross section
273 //@@ "dN=mmu+mmu2/G4NucleiProperties::GetNuclearMass(<G4double>(Z+N),<G4double>(Z)/GeV"
274 return dN;
275}

Referenced by GetCrossSection().

The documentation for this class was generated from the following files: