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
G4QANuANuNuclearCrossSection Class Reference

#include <G4QANuANuNuclearCrossSection.hh>

+ Inheritance diagram for G4QANuANuNuclearCrossSection:

Public Member Functions

 ~G4QANuANuNuclearCrossSection ()
 
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

 G4QANuANuNuclearCrossSection ()
 
- 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 G4QANuANuNuclearCrossSection.hh.

Constructor & Destructor Documentation

◆ G4QANuANuNuclearCrossSection()

G4QANuANuNuclearCrossSection::G4QANuANuNuclearCrossSection ( )
inlineprotected

Definition at line 53 of file G4QANuANuNuclearCrossSection.hh.

53{};

◆ ~G4QANuANuNuclearCrossSection()

G4QANuANuNuclearCrossSection::~G4QANuANuNuclearCrossSection ( )

Definition at line 78 of file G4QANuANuNuclearCrossSection.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 G4QANuANuNuclearCrossSection::CalculateCrossSection ( G4bool  CS,
G4int  F,
G4int  I,
G4int  PDG,
G4int  Z,
G4int  N,
G4double  Momentum 
)
virtual

Implements G4VQCrossSection.

Definition at line 262 of file G4QANuANuNuclearCrossSection.cc.

264{
265 static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2
266 static const G4int nE=65; // !! If change this, change it in GetFunctions() (*.hh) !!
267 static const G4int mL=nE-1;
268 static const G4double EMi=0.; // Universal threshold of the reaction in GeV
269 static const G4double EMa=300.; // Maximum tabulated Energy of nu_mu in GeV
270 // *** Begin of the Associative memory for acceleration of the cross section calculations
271 static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional)
272 static G4bool first=true; // Flag of initialization of the energy axis
273 // *** End of Static Definitions (Associative Memory) ***
274 //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
275 //G4double TotEnergy2=Momentum;
276 onlyCS=CS; // Flag to calculate only CS (not TX & QE)
277 lastE=Momentum/GeV; // Kinetic energy of the muon neutrino (in GeV!)
278 if (lastE<=EMi) // Energy is below the minimum energy in the table
279 {
280 lastE=0.;
281 lastSig=0.;
282 return 0.;
283 }
284 G4int Z=targZ; // New Z, which can change the sign
285 if(F<=0) // This isotope was not the last used isotop
286 {
287 if(F<0) // This isotope was found in DAMDB =-------=> RETRIEVE
288 {
289 lastTX =(*TX)[I]; // Pointer to the prepared TX function (same isotope)
290 lastQE =(*QE)[I]; // Pointer to the prepared QE function (same isotope)
291 }
292 else // This isotope wasn't calculated previously => CREATE
293 {
294 if(first)
295 {
296 lastEN = new G4double[nE]; // This must be done only once!
297 Z=-Z; // To explain GetFunctions that E-axis must be filled
298 first=false; // To make it only once
299 }
300 lastTX = new G4double[nE]; // Allocate memory for the new TX function
301 lastQE = new G4double[nE]; // Allocate memory for the new QE function
302 G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK)
303 if(res<0) G4cerr<<"*W*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
304 // *** The synchronization check ***
305 G4int sync=TX->size();
306 if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
307 TX->push_back(lastTX);
308 QE->push_back(lastQE);
309 } // End of creation of the new set of parameters
310 } // End of parameters udate
311 // =-----------------= NOW Calculate the Cross Section =---------------=
312 if (lastE<=EMi) // Check that antiNuEnergy is higher than ThreshE
313 {
314 lastE=0.;
315 lastSig=0.;
316 return 0.;
317 }
318 if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization
319 {
320 G4int chk=1;
321 G4int ran=mL/2;
322 G4int sep=ran; // as a result = an index of the left edge of the interval
323 while(ran>=2)
324 {
325 G4int newran=ran/2;
326 if(lastE<=lastEN[sep]) sep-=newran;
327 else sep+=newran;
328 ran=newran;
329 chk=chk+chk;
330 }
331 if(chk+chk!=mL) G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
332 G4double lowE=lastEN[sep];
333 G4double highE=lastEN[sep+1];
334 G4double lowTX=lastTX[sep];
335 if(lastE<lowE||sep>=mL||lastE>highE)
336 G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
337 <<", sep="<<sep<<", mL="<<mL<<G4endl;
338 lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
339 if(!onlyCS) // Skip the differential cross-section parameters
340 {
341 G4double lowQE=lastQE[sep];
342 lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
343#ifdef pdebug
344 G4cout<<"G4NuMuNuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
345#endif
346 }
347 }
348 else
349 {
350 lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
351 lastQEL=lastQE[mL];
352 }
353 if(lastQEL<0.) lastQEL = 0.;
354 if(lastSig<0.) lastSig = 0.;
355 // The cross-sections are expected to be in mb
356 lastSig*=mb38;
357 if(!onlyCS) lastQEL*=mb38;
358 return lastSig;
359}
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 G4QANuANuNuclearCrossSection::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 G4QANuANuNuclearCrossSection.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<<"G4QAMNCS::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 debug
111 G4cout<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"---G4QAMNCrossSec::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<<"G4QAMNCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lstI="<<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<<"G4QAMNCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
192#endif
193 if(pEn>lastTH)
194 {
195#ifdef pdebug
196 G4cout<<"G4QAMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
197#endif
198 lastTH=pEn;
199 }
200 }
201#ifdef pdebug
202 G4cout<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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}
G4QANuANuNuclearCrossSection::CalculateCrossSection
G4double CalculateCrossSection(G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum)
Definition: G4QANuANuNuclearCrossSection.cc:262
G4QANuANuNuclearCrossSection::ThresholdEnergy
G4double ThresholdEnergy(G4int Z, G4int N, G4int PDG=-14)
Definition: G4QANuANuNuclearCrossSection.cc:259
G4VQCrossSection::tolerance
static G4double tolerance
Definition: G4VQCrossSection.hh:175

◆ GetDirectPart()

G4double G4QANuANuNuclearCrossSection::GetDirectPart ( G4double  Q2)
virtual

Reimplemented from G4VQCrossSection.

Definition at line 725 of file G4QANuANuNuclearCrossSection.cc.

726{
727 G4double f=Q2/4.62;
728 G4double ff=f*f;
729 G4double r=ff*ff;
730 G4double s_value=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
731 //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
732 return 1.-s_value*(1.-s_value/2);
733}

◆ GetExchangePDGCode()

G4int G4QANuANuNuclearCrossSection::GetExchangePDGCode ( )
virtual

Reimplemented from G4VQCrossSection.

Definition at line 742 of file G4QANuANuNuclearCrossSection.cc.

742{return 22;}

◆ GetLastQELCS()

G4double G4QANuANuNuclearCrossSection::GetLastQELCS ( )
inlinevirtual

Reimplemented from G4VQCrossSection.

Definition at line 81 of file G4QANuANuNuclearCrossSection.hh.

81{return lastQEL;}

◆ GetLastTOTCS()

G4double G4QANuANuNuclearCrossSection::GetLastTOTCS ( )
inlinevirtual

Reimplemented from G4VQCrossSection.

Definition at line 80 of file G4QANuANuNuclearCrossSection.hh.

80{return lastSig;}

◆ GetNPartons()

G4double G4QANuANuNuclearCrossSection::GetNPartons ( G4double  Q2)
virtual

Reimplemented from G4VQCrossSection.

Definition at line 736 of file G4QANuANuNuclearCrossSection.cc.

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

◆ GetNQE_ExchangeQ2()

G4double G4QANuANuNuclearCrossSection::GetNQE_ExchangeQ2 ( )
virtual

Reimplemented from G4VQCrossSection.

Definition at line 501 of file G4QANuANuNuclearCrossSection.cc.

502{
503 static const double mpi=.13957018; // charged pi meson mass in GeV
504 static const double MN=.931494043; // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
505 static const double dMN=MN+MN; // 2*M_N in GeV
506 static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
507 static const G4int power=7; // direct power for the magic variable
508 static const G4double pconv=1./power; // conversion power for the magic variable
509 static const G4int nX=21; // #Of point in the Xl table (in GeV^2)
510 static const G4int lX=nX-1; // index of the last in the Xl table
511 static const G4int bX=lX-1; // @@ index of the before last in the Xl table
512 static const G4int nE=20; // #Of point in the El table (in GeV^2)
513 static const G4int bE=nE-1; // index of the last in the El table
514 static const G4int pE=bE-1; // index of the before last in the El table
515 // Reversed table
516 static const G4double X0[nX]={5.21412e-05,
517 .437860, .681908, .891529, 1.08434, 1.26751, 1.44494, 1.61915, 1.79198, 1.96493, 2.13937,
518 2.31664, 2.49816, 2.68559, 2.88097, 3.08705, 3.30774, 3.54917, 3.82233, 4.15131, 4.62182};
519 static const G4double X1[nX]={.00102591,
520 1.00443, 1.55828, 2.03126, 2.46406, 2.87311, 3.26723, 3.65199, 4.03134, 4.40835, 4.78561,
521 5.16549, 5.55031, 5.94252, 6.34484, 6.76049, 7.19349, 7.64917, 8.13502, 8.66246, 9.25086};
522 static const G4double X2[nX]={.0120304,
523 2.59903, 3.98637, 5.15131, 6.20159, 7.18024, 8.10986, 9.00426, 9.87265, 10.7217, 11.5564,
524 12.3808, 13.1983, 14.0116, 14.8234, 15.6359, 16.4515, 17.2723, 18.1006, 18.9386, 19.7892};
525 static const G4double X3[nX]={.060124,
526 5.73857, 8.62595, 10.9849, 13.0644, 14.9636, 16.7340, 18.4066, 20.0019, 21.5342, 23.0142,
527 24.4497, 25.8471, 27.2114, 28.5467, 29.8564, 31.1434, 32.4102, 33.6589, 34.8912, 36.1095};
528 static const G4double X4[nX]={.0992363,
529 8.23746, 12.1036, 15.1740, 17.8231, 20.1992, 22.3792, 24.4092, 26.3198, 28.1320, 29.8615,
530 31.5200, 33.1169, 34.6594, 36.1536, 37.6044, 39.0160, 40.3920, 41.7353, 43.0485, 44.3354};
531 static const G4double X5[nX]={.0561127,
532 7.33661, 10.5694, 13.0778, 15.2061, 17.0893, 18.7973, 20.3717, 21.8400, 23.2211, 24.5291,
533 25.7745, 26.9655, 28.1087, 29.2094, 30.2721, 31.3003, 32.2972, 33.2656, 34.2076, 35.1265};
534 static const G4double X6[nX]={.0145859,
535 4.81774, 6.83565, 8.37399, 9.66291, 10.7920, 11.8075, 12.7366, 13.5975, 14.4025, 15.1608,
536 15.8791, 16.5628, 17.2162, 17.8427, 18.4451, 19.0259, 19.5869, 20.1300, 20.6566, 21.1706};
537 static const G4double X7[nX]={.00241155,
538 2.87095, 4.02492, 4.89243, 5.61207, 6.23747, 6.79613, 7.30433, 7.77270, 8.20858, 8.61732,
539 9.00296, 9.36863, 9.71682, 10.0495, 10.3684, 10.6749, 10.9701, 11.2550, 11.5306, 11.7982};
540 static const G4double X8[nX]={.000316863,
541 1.76189, 2.44632, 2.95477, 3.37292, 3.73378, 4.05420, 4.34415, 4.61009, 4.85651, 5.08666,
542 5.30299, 5.50738, 5.70134, 5.88609, 6.06262, 6.23178, 6.39425, 6.55065, 6.70149, 6.84742};
543 static const G4double X9[nX]={3.73544e-05,
544 1.17106, 1.61289, 1.93763, 2.20259, 2.42976, 2.63034, 2.81094, 2.97582, 3.12796, 3.26949,
545 3.40202, 3.52680, 3.64482, 3.75687, 3.86360, 3.96557, 4.06323, 4.15697, 4.24713, 4.33413};
546 static const G4double XA[nX]={4.19131e-06,
547 .849573, 1.16208, 1.38955, 1.57379, 1.73079, 1.86867, 1.99221, 2.10451, 2.20770, 2.30332,
548 2.39252, 2.47622, 2.55511, 2.62977, 2.70066, 2.76818, 2.83265, 2.89437, 2.95355, 3.01051};
549 static const G4double XB[nX]={4.59981e-07,
550 .666131, .905836, 1.07880, 1.21796, 1.33587, 1.43890, 1.53080, 1.61399, 1.69011, 1.76040,
551 1.82573, 1.88682, 1.94421, 1.99834, 2.04959, 2.09824, 2.14457, 2.18878, 2.23107, 2.27162};
552 static const G4double XC[nX]={4.99861e-08,
553 .556280, .752730, .893387, 1.00587, 1.10070, 1.18317, 1.25643, 1.32247, 1.38269, 1.43809,
554 1.48941, 1.53724, 1.58203, 1.62416, 1.66391, 1.70155, 1.73728, 1.77128, 1.80371, 1.83473};
555 static const G4double XD[nX]={5.40832e-09,
556 .488069, .657650, .778236, .874148, .954621, 1.02432, 1.08599, 1.14138, 1.19172, 1.23787,
557 1.28049, 1.32008, 1.35705, 1.39172, 1.42434, 1.45514, 1.48429, 1.51197, 1.53829, 1.56339};
558 static const G4double XE[nX]={5.84029e-10,
559 .445057, .597434, .705099, .790298, .861468, .922865, .976982, 1.02542, 1.06930, 1.10939,
560 1.14630, 1.18050, 1.21233, 1.24208, 1.27001, 1.29630, 1.32113, 1.34462, 1.36691, 1.38812};
561 static const G4double XF[nX]={6.30137e-11,
562 .418735, .560003, .659168, .737230, .802138, .857898, .906854, .950515, .989915, 1.02580,
563 1.05873, 1.08913, 1.11734, 1.14364, 1.16824, 1.19133, 1.21306, 1.23358, 1.25298, 1.27139};
564 static const G4double XG[nX]={6.79627e-12,
565 .405286, .539651, .633227, .706417, .766929, .818642, .863824, .903931, .939963, .972639,
566 1.00250, 1.02995, 1.05532, 1.07887, 1.10082, 1.12134, 1.14058, 1.15867, 1.17572, 1.19183};
567 static const G4double XH[nX]={7.32882e-13,
568 .404391, .535199, .625259, .695036, .752243, .800752, .842823, .879906, .912994, .942802,
569 .969862, .994583, 1.01729, 1.03823, 1.05763, 1.07566, 1.09246, 1.10816, 1.12286, 1.13667};
570 static const G4double XI[nX]={7.90251e-14,
571 .418084, .548382, .636489, .703728, .758106, .803630, .842633, .876608, .906576, .933269,
572 .957233, .978886, .998556, 1.01651, 1.03295, 1.04807, 1.06201, 1.07489, 1.08683, 1.09792};
573 static const G4double XJ[nX]={8.52083e-15,
574 .447299, .579635, .666780, .731788, .783268, .825512, .861013, .891356, .917626, .940597,
575 .960842, .978802, .994820, 1.00917, 1.02208, 1.03373, 1.04427, 1.05383, 1.06253, 1.07046};
576 // Direct table
577 static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
578 X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
579 static const G4double dX[nE]={
580 (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
581 (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
582 (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
583 (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
584 (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
585 static const G4double* Xl[nE]=
586 {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
587 static const G4double I0[nX]={0,
588 .354631, 1.08972, 2.05138, 3.16564, 4.38343, 5.66828, 6.99127, 8.32858, 9.65998, 10.9680,
589 12.2371, 13.4536, 14.6050, 15.6802, 16.6686, 17.5609, 18.3482, 19.0221, 19.5752, 20.0000};
590 static const G4double I1[nX]={0,
591 .281625, .877354, 1.67084, 2.60566, 3.64420, 4.75838, 5.92589, 7.12829, 8.34989, 9.57708,
592 10.7978, 12.0014, 13.1781, 14.3190, 15.4162, 16.4620, 17.4496, 18.3724, 19.2245, 20.0000};
593 static const G4double I2[nX]={0,
594 .201909, .642991, 1.24946, 1.98463, 2.82370, 3.74802, 4.74263, 5.79509, 6.89474, 8.03228,
595 9.19947, 10.3889, 11.5938, 12.8082, 14.0262, 15.2427, 16.4527, 17.6518, 18.8356, 20.0000};
596 static const G4double I3[nX]={0,
597 .140937, .461189, .920216, 1.49706, 2.17728, 2.94985, 3.80580, 4.73758, 5.73867, 6.80331,
598 7.92637, 9.10316, 10.3294, 11.6013, 12.9150, 14.2672, 15.6548, 17.0746, 18.5239, 20.0000};
599 static const G4double I4[nX]={0,
600 .099161, .337358, .694560, 1.16037, 1.72761, 2.39078, 3.14540, 3.98768, 4.91433, 5.92245,
601 7.00942, 8.17287, 9.41060, 10.7206, 12.1010, 13.5500, 15.0659, 16.6472, 18.2924, 20.0000};
602 static const G4double I5[nX]={0,
603 .071131, .255084, .543312, .932025, 1.41892, 2.00243, 2.68144, 3.45512, 4.32283, 5.28411,
604 6.33859, 7.48602, 8.72621, 10.0590, 11.4844, 13.0023, 14.6128, 16.3158, 18.1115, 20.0000};
605 static const G4double I6[nX]={0,
606 .053692, .202354, .443946, .778765, 1.20774, 1.73208, 2.35319, 3.07256, 3.89177, 4.81249,
607 5.83641, 6.96528, 8.20092, 9.54516, 10.9999, 12.5670, 14.2486, 16.0466, 17.9630, 20.0000};
608 static const G4double I7[nX]={0,
609 .043065, .168099, .376879, .672273, 1.05738, 1.53543, 2.10973, 2.78364, 3.56065, 4.44429,
610 5.43819, 6.54610, 7.77186, 9.11940, 10.5928, 12.1963, 13.9342, 15.8110, 17.8313, 20.0000};
611 static const G4double I8[nX]={0,
612 .036051, .143997, .327877, .592202, .941572, 1.38068, 1.91433, 2.54746, 3.28517, 4.13277,
613 5.09574, 6.17984, 7.39106, 8.73568, 10.2203, 11.8519, 13.6377, 15.5854, 17.7033, 20.0000};
614 static const G4double I9[nX]={0,
615 .030977, .125727, .289605, .528146, .846967, 1.25183, 1.74871, 2.34384, 3.04376, 3.85535,
616 4.78594, 5.84329, 7.03567, 8.37194, 9.86163, 11.5150, 13.3430, 15.3576, 17.5719, 20.0000};
617 static const G4double IA[nX]={0,
618 .027129, .111420, .258935, .475812, .768320, 1.14297, 1.60661, 2.16648, 2.83034, 3.60650,
619 4.50394, 5.53238, 6.70244, 8.02569, 9.51488, 11.1841, 13.0488, 15.1264, 17.4362, 20.0000};
620 static const G4double IB[nX]={0,
621 .024170, .100153, .234345, .433198, .703363, 1.05184, 1.48607, 2.01409, 2.64459, 3.38708,
622 4.25198, 5.25084, 6.39647, 7.70319, 9.18708, 10.8663, 12.7617, 14.8968, 17.2990, 20.0000};
623 static const G4double IC[nX]={0,
624 .021877, .091263, .214670, .398677, .650133, .976322, 1.38510, 1.88504, 2.48555, 3.19709,
625 4.03129, 5.00127, 6.12184, 7.40989, 8.88482, 10.5690, 12.4888, 14.6748, 17.1638, 20.0000};
626 static const G4double ID[nX]={0,
627 .020062, .084127, .198702, .370384, .606100, .913288, 1.30006, 1.77535, 2.34912, 3.03253,
628 3.83822, 4.78063, 5.87634, 7.14459, 8.60791, 10.2929, 12.2315, 14.4621, 17.0320, 20.0000};
629 static const G4double IE[nX]={0,
630 .018547, .078104, .185102, .346090, .567998, .858331, 1.22535, 1.67824, 2.22735, 2.88443,
631 3.66294, 4.57845, 5.64911, 6.89637, 8.34578, 10.0282, 11.9812, 14.2519, 16.8993, 20.0000};
632 static const G4double IF[nX]={0,
633 .017143, .072466, .172271, .323007, .531545, .805393, 1.15288, 1.58338, 2.10754, 2.73758,
634 3.48769, 4.37450, 5.41770, 6.64092, 8.07288, 9.74894, 11.7135, 14.0232, 16.7522, 20.0000};
635 static const G4double IG[nX]={0,
636 .015618, .066285, .158094, .297316, .490692, .745653, 1.07053, 1.47479, 1.96931, 2.56677,
637 3.28205, 4.13289, 5.14068, 6.33158, 7.73808, 9.40133, 11.3745, 13.7279, 16.5577, 20.0000};
638 static const G4double IH[nX]={0,
639 .013702, .058434, .139923, .264115, .437466, .667179, .961433, 1.32965, 1.78283, 2.33399,
640 2.99871, 3.79596, 4.74916, 5.88771, 7.24937, 8.88367, 10.8576, 13.2646, 16.2417, 20.0000};
641 static const G4double II[nX]={0,
642 .011264, .048311, .116235, .220381, .366634, .561656, .813132, 1.13008, 1.52322, 2.00554,
643 2.59296, 3.30542, 4.16834, 5.21490, 6.48964, 8.05434, 9.99835, 12.4580, 15.6567, 20.0000};
644 static const G4double IJ[nX]={0,
645 .008628, .037206, .089928, .171242, .286114, .440251, .640343, .894382, 1.21208, 1.60544,
646 2.08962, 2.68414, 3.41486, 4.31700, 5.44048, 6.85936, 8.69067, 11.1358, 14.5885, 20.0000};
647 static const G4double* Il[nE]=
648 {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
649 static const G4double lE[nE]={
650-1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
651 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
652 static const G4double lEmi=lE[0];
653 static const G4double lEma=lE[nE-1];
654 static const G4double dlE=(lEma-lEmi)/bE;
655 //***************************************************************************************
656 G4double Enu=lastE; // Get energy of the last calculated cross-section
657 G4double lEn=std::log(Enu); // log(E) for interpolation
658 G4double rE=(lEn-lEmi)/dlE; // Position of the energy
659 G4int fE=static_cast<int>(rE); // Left bin for interpolation
660 if(fE<0) fE=0;
661 if(fE>pE)fE=pE;
662 G4int sE=fE+1; // Right bin for interpolation
663 G4double dE=rE-fE; // relative log shift from the left bin
664 G4double dEnu=Enu+Enu; // doubled energy of nu/anu
665 G4double Enu2=Enu*Enu; // squared energy of nu/anu
666 G4double Emu=Enu; // Free Energy of neutrino/anti-neutrino
667 G4double ME=Enu*MN; // M*E
668 G4double dME=ME+ME; // 2*M*E
669 G4double dEMN=(dEnu+MN)*ME;
670 G4double sqE=Enu*ME;
671 G4double E2M=MN*Enu2;
672 G4double ymax=(E2M+sqE)/dEMN;
673 G4double Q2mi=0.; // Q2_min(E_nu)
674 G4double Q2ma=dME*ymax; // Q2_max(E_nu)
675 G4double Q2nq=Emu*dMN-mcV;
676 if(Q2ma>Q2nq) Q2ma=Q2nq; // Correction for Non Quasi Elastic
677 // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
678 G4double Rmi=Q2mi/Q2ma;
679 G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78);
680 // --- E-interpolation must be done in a log scale ---
681 G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
682 G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
683 // Find the integral values integ(Xmi) & integ(Xma) using the direct table
684 G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
685 G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
686 G4double rXi=(Xmi-iXmi)/idX;
687 G4int iXi=static_cast<int>(rXi);
688 if(iXi<0) iXi=0;
689 if(iXi>bX) iXi=bX;
690 G4double dXi=rXi-iXi;
691 G4double bntil=Il[fE][iXi];
692 G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
693 G4double bntir=Il[sE][iXi];
694 G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
695 G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
696 //
697 G4double rXa=(Xma-iXmi)/idX;
698 G4int iXa=static_cast<int>(rXa);
699 if(iXa<0) iXa=0;
700 if(iXa>bX) iXa=bX;
701 G4double dXa=rXa-iXa;
702 G4double bntal=Il[fE][iXa];
703 G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
704 G4double bntar=Il[sE][iXa];
705 G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
706 G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
707 //
708 // *** Find X using the reversed table ***
709 G4double intx=inti+(inta-inti)*G4UniformRand();
710 G4int intc=static_cast<int>(intx);
711 if(intc<0) intc=0;
712 if(intc>bX) intc=bX;
713 G4double dint=intx-intc;
714 G4double mXl=Xl[fE][intc];
715 G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
716 G4double mXr=Xl[sE][intc];
717 G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
718 G4double X=Xlb+dE*(Xrb-Xlb); // interpolated X value
719 G4double R=shift-std::pow(X,pconv);
720 G4double Q2=R*Q2ma;
721 return Q2*GeV*GeV;
722}
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53

◆ GetPointer()

G4VQCrossSection * G4QANuANuNuclearCrossSection::GetPointer ( )
static

Definition at line 72 of file G4QANuANuNuclearCrossSection.cc.

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

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

◆ GetQEL_ExchangeQ2()

G4double G4QANuANuNuclearCrossSection::GetQEL_ExchangeQ2 ( )
virtual

Reimplemented from G4VQCrossSection.

Definition at line 424 of file G4QANuANuNuclearCrossSection.cc.

425{
426 static const double MN=.931494043; // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
427 static const G4double power=-3.5; // direct power for the magic variable
428 static const G4double pconv=1./power;// conversion power for the magic variable
429 static const G4int nQ2=101; // #Of point in the Q2l table (in GeV^2)
430 static const G4int lQ2=nQ2-1; // index of the last in the Q2l table
431 static const G4int bQ2=lQ2-1; // index of the before last in the Q2 ltable
432 // Reversed table
433 static const G4double Xl[nQ2]={5.20224e-16,
434 .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898,
435 .093878, .103332, .112844, .122410, .132023, .141680, .151376, .161109, .170875, .180672,
436 .190499, .200352, .210230, .220131, .230055, .239999, .249963, .259945, .269944, .279960,
437 .289992, .300039, .310099, .320173, .330260, .340359, .350470, .360592, .370724, .380867,
438 .391019, .401181, .411352, .421531, .431719, .441915, .452118, .462329, .472547, .482771,
439 .493003, .503240, .513484, .523734, .533989, .544250, .554517, .564788, .575065, .585346,
440 .595632, .605923, .616218, .626517, .636820, .647127, .657438, .667753, .678072, .688394,
441 .698719, .709048, .719380, .729715, .740053, .750394, .760738, .771085, .781434, .791786,
442 .802140, .812497, .822857, .833219, .843582, .853949, .864317, .874687, .885060, .895434,
443 .905810, .916188, .926568, .936950, .947333, .957719, .968105, .978493, .988883, .999275};
444 // Direct table
445 static const G4double Xmax=Xl[lQ2];
446 static const G4double Xmin=Xl[0];
447 static const G4double dX=(Xmax-Xmin)/lQ2; // step in X(Q2, GeV^2)
448 static const G4double inl[nQ2]={0,
449 1.52225, 2.77846, 3.96651, 5.11612, 6.23990, 7.34467, 8.43466, 9.51272, 10.5809, 11.6406,
450 12.6932, 13.7394, 14.7801, 15.8158, 16.8471, 17.8743, 18.8979, 19.9181, 20.9353, 21.9496,
451 22.9614, 23.9707, 24.9777, 25.9826, 26.9855, 27.9866, 28.9860, 29.9837, 30.9798, 31.9745,
452 32.9678, 33.9598, 34.9505, 35.9400, 36.9284, 37.9158, 38.9021, 39.8874, 40.8718, 41.8553,
453 42.8379, 43.8197, 44.8007, 45.7810, 46.7605, 47.7393, 48.7174, 49.6950, 50.6718, 51.6481,
454 52.6238, 53.5990, 54.5736, 55.5476, 56.5212, 57.4943, 58.4670, 59.4391, 60.4109, 61.3822,
455 62.3531, 63.3236, 64.2937, 65.2635, 66.2329, 67.2019, 68.1707, 69.1390, 70.1071, 71.0748,
456 72.0423, 73.0095, 73.9763, 74.9429, 75.9093, 76.8754, 77.8412, 78.8068, 79.7721, 80.7373,
457 81.7022, 82.6668, 83.6313, 84.5956, 85.5596, 86.5235, 87.4872, 88.4507, 89.4140, 90.3771,
458 91.3401, 92.3029, 93.2656, 94.2281, 95.1904, 96.1526, 97.1147, 98.0766, 99.0384, 100.000};
459 G4double Enu=lastE; // Get energy of the last calculated cross-section
460 G4double dEnu=Enu+Enu; // doubled energy of nu/anu
461 G4double Enu2=Enu*Enu; // squared energy of nu/anu
462 G4double ME=Enu*MN; // M*E
463 G4double dME=ME+ME; // 2*M*E
464 G4double dEMN=(dEnu+MN)*ME;
465 G4double sqE=Enu*ME;
466 G4double E2M=MN*Enu2;
467 G4double ymax=(E2M+sqE)/dEMN;
468 G4double Q2mi=0.; // Q2_min(E_nu)
469 G4double Q2ma=dME*ymax; // Q2_max(E_nu)
470 G4double Xma=std::pow((1.+Q2mi),power); // X_max(E_nu)
471 G4double Xmi=std::pow((1.+Q2ma),power); // X_min(E_nu)
472 // Find the integral values integ(Xmi) & integ(Xma) using the direct table
473 G4double rXi=(Xmi-Xmin)/dX;
474 G4int iXi=static_cast<int>(rXi);
475 if(iXi<0) iXi=0;
476 if(iXi>bQ2) iXi=bQ2;
477 G4double dXi=rXi-iXi;
478 G4double bnti=inl[iXi];
479 G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
480 //
481 G4double rXa=(Xma-Xmin)/dX;
482 G4int iXa=static_cast<int>(rXa);
483 if(iXa<0) iXa=0;
484 if(iXa>bQ2) iXa=bQ2;
485 G4double dXa=rXa-iXa;
486 G4double bnta=inl[iXa];
487 G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
488 // *** Find X using the reversed table ***
489 G4double intx=inti+(inta-inti)*G4UniformRand();
490 G4int intc=static_cast<int>(intx);
491 if(intc<0) intc=0;
492 if(intc>bQ2) intc=bQ2; // If it is more than max, then the BAD extrapolation
493 G4double dint=intx-intc;
494 G4double mX=Xl[intc];
495 G4double X=mX+dint*(Xl[intc+1]-mX);
496 G4double Q2=std::pow(X,pconv)-1.;
497 return Q2*GeV*GeV;
498}

◆ ThresholdEnergy()

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

Reimplemented from G4VQCrossSection.

Definition at line 259 of file G4QANuANuNuclearCrossSection.cc.

259{return 0.;}

Referenced by GetCrossSection().

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