Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4QANuANuNuclearCrossSection.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id$
28//
29//
30// G4 Physics class: G4QANuANuNuclearCrossSection for gamma+A cross sections
31// Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01
32// The last update: M.V. Kossov, CERN/ITEP (Moscow) 17-Oct-03
33//
34// ****************************************************************************************
35// ********** This CLASS is temporary moved from the photolepton_hadron directory *********
36// ******* DO NOT MAKE ANY CHANGE! With time it'll move back to photolepton...(M.K.) ******
37// ****************************************************************************************
38// ----------------------------------------------------------------------------------------
39// Short description: Provides the (anti-nu,anti-nu)A cross-section (Created by M. Kossov)
40// ----------------------------------------------------------------------------------------
41
42//#define debug
43//#define edebug
44//#define pdebug
45//#define ppdebug
46//#define tdebug
47//#define sdebug
48
49#include "G4QANuANuNuclearCrossSection.hh"
50#include "G4SystemOfUnits.hh"
51
52// Initialization of the
53G4bool G4QANuANuNuclearCrossSection::onlyCS=true;//Flag to calculate only CS (not QE)
54G4double G4QANuANuNuclearCrossSection::lastSig=0.;//Last calculated total cross section
55G4double G4QANuANuNuclearCrossSection::lastQEL=0.;//Last calculated quasi-el cross section
56G4int G4QANuANuNuclearCrossSection::lastL=0; //Last used in cross section TheLastBin
57G4double G4QANuANuNuclearCrossSection::lastE=0.; //Last used in cross section TheEnergy
58G4double* G4QANuANuNuclearCrossSection::lastEN=0; //Pointer to theEnergy Scale of TX & QE
59G4double* G4QANuANuNuclearCrossSection::lastTX=0; //Pointer to theLastArray of TX function
60G4double* G4QANuANuNuclearCrossSection::lastQE=0; //Pointer to theLastArray of QE function
61G4int G4QANuANuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
62G4int G4QANuANuNuclearCrossSection::lastN=0; // The last N of calculated nucleus
63G4int G4QANuANuNuclearCrossSection::lastZ=0; // The last Z of calculated nucleus
64G4double G4QANuANuNuclearCrossSection::lastP=0.; // Last used in cross section Momentum
65G4double G4QANuANuNuclearCrossSection::lastTH=0.; // Last threshold momentum
66G4double G4QANuANuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
67G4int G4QANuANuNuclearCrossSection::lastI=0; // The last position in the DAMDB
68std::vector<G4double*>* G4QANuANuNuclearCrossSection::TX = new std::vector<G4double*>;
69std::vector<G4double*>* G4QANuANuNuclearCrossSection::QE = new std::vector<G4double*>;
70
71// Returns Pointer to the G4VQCrossSection class
72G4VQCrossSection* G4QANuANuNuclearCrossSection::GetPointer()
73{
74 static G4QANuANuNuclearCrossSection theCrossSection; //*Static body of the Cross Section*
75 return &theCrossSection;
76}
77
78G4QANuANuNuclearCrossSection::~G4QANuANuNuclearCrossSection()
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}
87
88// The main member function giving the collision cross section (P is in IU, CS is in mb)
89// Make pMom in independent units ! (Now it is MeV)
90G4double G4QANuANuNuclearCrossSection::GetCrossSection(G4bool fCS, G4double pMom,
91 G4int tgZ, G4int tgN, G4int pPDG)
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}
257
258// Gives the threshold energy = the same for all nuclei (@@ can be reduced for hevy nuclei)
259G4double G4QANuANuNuclearCrossSection::ThresholdEnergy(G4int, G4int, G4int) {return 0.;}
260
261// The main member function giving the gamma-A cross section (E_kin in MeV, CS in mb)
262G4double G4QANuANuNuclearCrossSection::CalculateCrossSection(G4bool CS, G4int F, G4int I,
263 G4int , G4int targZ, G4int targN, G4double Momentum)
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}
360
361// Calculate the cros-section functions
362// ****************************************************************************************
363// *** This tables are the same for all lepto-nuclear reactions, only mass is different ***
364// ***@@ IT'S REASONABLE TO MAKE ADDiTIONAL VIRTUAL CLASS FOR LEPTO-NUCLEAR WITH THIS@@ ***
365// ****************************************************************************************
366G4int G4QANuANuNuclearCrossSection::GetFunctions (G4int z, G4int n,
367 G4double* t, G4double* q, G4double* e)
368{
369 static const G4int nE=65; // !! If change this, change it in GetCrossSection() (*.cc) !!
370 static const G4double nuEn[nE]={0.,
371 1.00463e-5,1.05336e-5,1.10692e-5,1.16592e-5,1.23109e-5,1.30323e-5,1.38331e-5,1.47245e-5,
372 1.57194e-5,1.68335e-5,1.80848e-5,1.94948e-5,2.10894e-5,2.28991e-5,2.49608e-5,2.73189e-5,
373 3.00273e-5,3.31516e-5,3.67722e-5,4.09881e-5,4.59217e-5,5.17255e-5,5.85908e-5,6.67583e-5,
374 7.65338e-5,8.83078e-5,.000102583,.000120011,.000141441,.000167995,.000201160,.000242926,
375 .000295985,.000364008,.000452051,.000567152,.000719210,.000922307,.001196710,.001571930,
376 .002091530,.002820590,.003857810,.005354930,.007548840,.010815300,.015760100,.023376900,
377 .035325600,.054430800,.085595700,.137508000,.225898000,.379892000,.654712000,1.15767000,
378 2.10277000,3.92843000,7.55861000,14.9991000,30.7412000,65.1734000,143.155000,326.326000};
379 static const G4double TOTX[nE]={0.,
380 3.63538e-5,3.81165e-5,4.00539e-5,4.21884e-5,4.45454e-5,4.71548e-5,5.00511e-5,5.32747e-5,
381 5.68730e-5,6.09016e-5,6.54261e-5,7.05246e-5,7.62894e-5,8.28315e-5,9.02838e-5,9.88066e-5,
382 .000108594,.000119884,.000132964,.000148191,.000166007,.000186959,.000211736,.000241202,
383 .000276453,.000318890,.000370311,.000433042,.000510116,.000605516,.000724514,.000874144,
384 .001063870,.001306520,.001619660,.002027520,.002563780,.003275710,.004230080,.005521890,
385 .007286920,.009719890,.013099700,.018695100,.029208400,.042095000,.059253700,.082373900,
386 .113071000,.151041000,.191803000,.224208000,.234187000,.217774000,.187139000,.157818000,
387 .137494000,.125872000,.120462000,.119148000,.120418000,.123027000,.126408000,.129071000};
388 static const G4double QELX[nE]={0.,
389 .365220e-9,.401502e-9,.443364e-9,.491885e-9,.548393e-9,.614536e-9,.692362e-9,.784441e-9,
390 .894012e-9,1.02519e-9,1.18322e-9,1.37487e-9,1.60890e-9,1.89677e-9,2.25355e-9,2.69928e-9,
391 3.26079e-9,3.97433e-9,4.88937e-9,6.07407e-9,7.62331e-9,9.67058e-9,1.24058e-8,1.61022e-8,
392 2.11580e-8,2.81605e-8,3.79876e-8,5.19696e-8,7.21515e-8,1.01724e-7,1.45743e-7,2.12353e-7,
393 3.14890e-7,4.75585e-7,7.32171e-7,1.14991e-6,1.84390e-6,3.02121e-6,5.06217e-6,8.68002e-6,
394 1.52408e-5,2.74159e-5,5.05363e-5,.000100111,.000220489,.000455269,.000933841,.001925650,
395 .003994300,.008221270,.016417600,.030830400,.052902400,.082519200,.115560000,.149598000,
396 .184112000,.215102000,.238253000,.252949000,.261267000,.265626000,.267782000,.268791000};
397 // --------------------------------
398 G4int first=0;
399 if(z<0.)
400 {
401 first=1;
402 z=-z;
403 }
404 if(z<1 || z>92) // neutron & plutonium are forbidden
405 {
406 G4cout<<"*G4QANuANuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
407 return -1;
408 }
409 for(G4int k=0; k<nE; k++)
410 {
411 G4double a=n+z;
412 G4double za=z+a;
413 G4double dz=z+z;
414 G4double da=a+a;
415 G4double ta=da+a;
416 if(first) e[k]=nuEn[k]; // Energy of neutrino E (first bin k=0 can be modified)
417 t[k]=TOTX[k]*nuEn[k]*(za+za)/ta+QELX[k]*(dz+dz-da)/ta; // TotalCrossSection
418 q[k]=QELX[k]*dz/a; // QuasiElasticCrossSection
419 }
420 return first;
421}
422
423// Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic
424G4double G4QANuANuNuclearCrossSection::GetQEL_ExchangeQ2()
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}
499
500// Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic
501G4double G4QANuANuNuclearCrossSection::GetNQE_ExchangeQ2()
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}
723
724// It returns a fraction of the direct interaction of the neutrino with quark-partons
725G4double G4QANuANuNuclearCrossSection::GetDirectPart(G4double Q2)
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}
734
735// #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
736G4double G4QANuANuNuclearCrossSection::GetNPartons(G4double Q2)
737{
738 return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
739}
740
741// This class can provide only virtual exchange by gamma (a substitute for Z0 boson)
742G4int G4QANuANuNuclearCrossSection::GetExchangePDGCode() {return 22;}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4QANuANuNuclearCrossSection::CalculateCrossSection
G4double CalculateCrossSection(G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N, G4double Momentum)
Definition: G4QANuANuNuclearCrossSection.cc:262
G4QANuANuNuclearCrossSection::GetPointer
static G4VQCrossSection * GetPointer()
Definition: G4QANuANuNuclearCrossSection.cc:72
G4QANuANuNuclearCrossSection::ThresholdEnergy
G4double ThresholdEnergy(G4int Z, G4int N, G4int PDG=-14)
Definition: G4QANuANuNuclearCrossSection.cc:259
G4QANuANuNuclearCrossSection::GetCrossSection
virtual G4double GetCrossSection(G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=0)
Definition: G4QANuANuNuclearCrossSection.cc:90
G4VQCrossSection::tolerance
static G4double tolerance
Definition: G4VQCrossSection.hh:175