source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/cross_sections/src/G4QNuMuNuclearCrossSection.cc @ 1353

Last change on this file since 1353 was 1340, checked in by garnier, 14 years ago

update ti head
File size: 40.8 KB
Line 
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: G4QNuMuNuclearCrossSection.cc,v 1.1 2009/11/16 18:15:43 mkossov Exp $
28// GEANT4 tag $Name: hadr-chips-V09-03-08 $
29//
30//
31// G4 Physics class: G4QNuMuNuclearCrossSection for (nu,mu-)A cross sections
32// Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01
33// The last update: M.V. Kossov, CERN/ITEP (Moscow) 17-Oct-03
34//
35// ****************************************************************************************
36// ********** This CLASS is temporary moved from the photolepton_hadron directory *********
37// ******* DO NOT MAKE ANY CHANGE! With time it'll move back to photolepton...(M.K.) ******
38// ****************************************************************************************
39//=========================================================================================
40
41//#define debug
42//#define edebug
43//#define pdebug
44//#define ppdebug
45//#define tdebug
46//#define sdebug
47
48#include "G4QNuMuNuclearCrossSection.hh"
49
50// Initialization of the
51G4bool    G4QNuMuNuclearCrossSection::onlyCS=true;// Flag to calculate only CS (not QE)
52G4double  G4QNuMuNuclearCrossSection::lastSig=0.;// Last calculated total cross section
53G4double  G4QNuMuNuclearCrossSection::lastQEL=0.;// Last calculated quasi-el. cross section
54G4int     G4QNuMuNuclearCrossSection::lastL=0;   // Last used in cross section TheLastBin
55G4double  G4QNuMuNuclearCrossSection::lastE=0.;  // Last used in cross section TheEnergy
56G4double* G4QNuMuNuclearCrossSection::lastEN=0;  // Pointer to the Energy Scale of TX & QE
57G4double* G4QNuMuNuclearCrossSection::lastTX=0;  // Pointer to the LastArray of TX function
58G4double* G4QNuMuNuclearCrossSection::lastQE=0;  // Pointer to the LastArray of QE function
59G4int     G4QNuMuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
60G4int     G4QNuMuNuclearCrossSection::lastN=0;   // The last N of calculated nucleus
61G4int     G4QNuMuNuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
62G4double  G4QNuMuNuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
63G4double  G4QNuMuNuclearCrossSection::lastTH=0.; // Last threshold momentum
64G4double  G4QNuMuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
65G4int     G4QNuMuNuclearCrossSection::lastI=0;   // The last position in the DAMDB
66std::vector<G4double*>* G4QNuMuNuclearCrossSection::TX = new std::vector<G4double*>;
67std::vector<G4double*>* G4QNuMuNuclearCrossSection::QE = new std::vector<G4double*>;
68
69// Returns Pointer to the G4VQCrossSection class
70G4VQCrossSection* G4QNuMuNuclearCrossSection::GetPointer()
71{
72  static G4QNuMuNuclearCrossSection theCrossSection; //**Static body of the Cross Section**
73  return &theCrossSection;
74}
75
76G4QNuMuNuclearCrossSection::~G4QNuMuNuclearCrossSection()
77{
78  G4int lens=TX->size();
79  for(G4int i=0; i<lens; ++i) delete[] (*TX)[i];
80  delete TX;
81  G4int hens=QE->size();
82  for(G4int i=0; i<hens; ++i) delete[] (*QE)[i];
83  delete QE;
84}
85
86// The main member function giving the collision cross section (P is in IU, CS is in mb)
87// Make pMom in independent units ! (Now it is MeV)
88G4double G4QNuMuNuclearCrossSection::GetCrossSection(G4bool fCS, G4double pMom,
89                                                     G4int tgZ, G4int tgN, G4int pPDG)
90{
91  static G4int j;                      // A#0f records found in DB for this projectile
92  static std::vector <G4int>    colPDG;// Vector of the projectile PDG code
93  static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
94  static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
95  static std::vector <G4double> colP;  // Vector of last momenta for the reaction
96  static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
97  static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
98  // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
99  G4double pEn=pMom;
100#ifdef debug
101  G4cout<<"G4QNMNCS::GetCS:>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="<<tgN
102        <<"("<<lastN<<"),PDG="<<pPDG<<"("<<lastPDG<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="
103        <<colN.size()<<G4endl;
104  //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
105#endif
106  if(pPDG!=14)
107  {
108#ifdef pdebug
109    G4cout<<"G4QNMNCS::GetCS: *** Found pPDG="<<pPDG<<" ====> CS=0"<<G4endl;
110    //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
111#endif
112    return 0.;                         // projectile PDG=0 is a mistake (?!) @@
113  }
114  G4bool in=false;                     // By default the isotope must be found in the AMDB
115  if(tgN!=lastN || tgZ!=lastZ || pPDG!=lastPDG)// The nucleus was not the last used isotope
116  {
117    in = false;                        // By default the isotope haven't be found in AMDB 
118    lastP   = 0.;                      // New momentum history (nothing to compare with)
119    lastPDG = pPDG;                    // The last PDG of the projectile
120    lastN   = tgN;                     // The last N of the calculated nucleus
121    lastZ   = tgZ;                     // The last Z of the calculated nucleus
122    lastI   = colN.size();             // Size of the Associative Memory DB in the heap
123    j  = 0;                            // A#0f records found in DB for this projectile
124    if(lastI) for(G4int i=0; i<lastI; i++) if(colPDG[i]==pPDG) // The partType is found
125    {                                  // The nucleus with projPDG is found in AMDB
126      if(colN[i]==tgN && colZ[i]==tgZ)
127      {
128        lastI=i;
129        lastTH =colTH[i];                // Last THreshold (A-dependent)
130#ifdef pdebug
131        G4cout<<"G4QNMNCS::GetCS:*Found*P="<<pMom<<",Threshold="<<lastTH<<",j="<<j<<G4endl;
132        //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
133#endif
134        if(pEn<=lastTH)
135        {
136#ifdef pdebug
137          G4cout<<"G4QNMNCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",X=0"<<G4endl;
138          //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
139#endif
140          return 0.;                     // Energy is below the Threshold value
141        }
142        lastP  =colP [i];                // Last Momentum  (A-dependent)
143        lastCS =colCS[i];                // Last CrossSect (A-dependent)
144        if(std::fabs(lastP/pMom-1.)<tolerance)
145        {
146#ifdef pdebug
147          G4cout<<"G4QNMNCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
148          //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
149#endif
150          return lastCS*millibarn;     // Use theLastCS
151        }
152        in = true;                       // This is the case when the isotop is found in DB
153        // Momentum pMom is in IU ! @@ Units
154#ifdef pdebug
155        G4cout<<"G4QNMNCS::G:UpdaDB P="<<pMom<<",f="<<fCS<<",lI="<<lastI<<",j="<<j<<G4endl;
156#endif
157        lastCS=CalculateCrossSection(fCS,-1,j,lastPDG,lastZ,lastN,pMom); // read & update
158#ifdef pdebug
159        G4cout<<"G4QNMNCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
160        //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
161#endif
162        if(lastCS<=0. && pEn>lastTH)    // Correct the threshold
163        {
164#ifdef pdebug
165          G4cout<<"G4QNMNCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
166#endif
167          lastTH=pEn;
168        }
169        break;                           // Go out of the LOOP
170      }
171#ifdef pdebug
172      G4cout<<"---G4QNMNCrossSec::GetCrosSec:pPDG="<<pPDG<<",j="<<j<<",N="<<colN[i]
173            <<",Z["<<i<<"]="<<colZ[i]<<",cPDG="<<colPDG[i]<<G4endl;
174      //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
175#endif
176      j++;                             // Increment a#0f records found in DB for this pPDG
177    }
178    if(!in)                            // This nucleus has not been calculated previously
179    {
180#ifdef pdebug
181      G4cout<<"G4QNMNCS::GetCrSec: CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
182#endif
183      //!!The slave functions must provide cross-sections in millibarns (mb) !! (not in IU)
184      lastCS=CalculateCrossSection(fCS,0,j,lastPDG,lastZ,lastN,pMom); //calculate & create
185      if(lastCS<=0.)
186      {
187        lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
188#ifdef pdebug
189        G4cout<<"G4QNMNCrossSection::GetCrossSect:NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
190#endif
191        if(pEn>lastTH)
192        {
193#ifdef pdebug
194          G4cout<<"G4QNMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
195#endif
196          lastTH=pEn;
197        }
198      }
199#ifdef pdebug
200      G4cout<<"G4QNMNCS::GetCrosSec:New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
201      //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
202#endif
203      colN.push_back(tgN);
204      colZ.push_back(tgZ);
205      colPDG.push_back(pPDG);
206      colP.push_back(pMom);
207      colTH.push_back(lastTH);
208      colCS.push_back(lastCS);
209#ifdef pdebug
210      G4cout<<"G4QNMNCS::GetCS:1st,P="<<pMom<<"(MeV),X="<<lastCS*millibarn<<"(mb)"<<G4endl;
211      //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
212#endif
213      return lastCS*millibarn;
214    } // End of creation of the new set of parameters
215    else
216    {
217#ifdef pdebug
218      G4cout<<"G4QNMNCS::GetCS: Update lastI="<<lastI<<",j="<<j<<G4endl;
219#endif
220      colP[lastI]=pMom;
221      colPDG[lastI]=pPDG;
222      colCS[lastI]=lastCS;
223    }
224  } // End of parameters udate
225  else if(pEn<=lastTH)
226  {
227#ifdef pdebug
228    G4cout<<"G4QNMNCS::GetCS: Current T="<<pEn<<" < Threshold="<<lastTH<<", CS=0"<<G4endl;
229    //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
230#endif
231    return 0.;                         // Momentum is below the Threshold Value -> CS=0
232  }
233  else if(std::fabs(lastP/pMom-1.)<tolerance)
234  {
235#ifdef pdebug
236    G4cout<<"G4QNMNCS::GetCS:OldCur P="<<pMom<<"="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
237    //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
238#endif
239    return lastCS*millibarn;     // Use theLastCS
240  }
241  else
242  {
243#ifdef pdebug
244    G4cout<<"G4QNMNCS::GetCS:UpdaCur P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",j="<<j<<G4endl;
245#endif
246    lastCS=CalculateCrossSection(fCS,1,j,lastPDG,lastZ,lastN,pMom); // Only UpdateDB
247    lastP=pMom;
248  }
249#ifdef pdebug
250  G4cout<<"G4QNMNCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
251  //CalculateCrossSection(fCS,-27,j,lastPDG,lastZ,lastN,pMom); // DUMMY TEST
252#endif
253  return lastCS*millibarn;
254}
255
256// Gives the threshold energy = the same for all nuclei (@@ can be reduced for hevy nuclei)
257G4double G4QNuMuNuclearCrossSection::ThresholdEnergy(G4int Z, G4int N, G4int)
258{
259  //static const G4double mNeut = G4NucleiProperties::GetNuclearMass(1,0)/GeV;
260  //static const G4double mProt = G4NucleiProperties::GetNuclearMass(1,1)/GeV;
261  //static const G4double mDeut = G4NucleiProperties::GetNuclearMass(2,1)/GeV/2.;
262  static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
263  static const G4double dmN=mN+mN;    // Doubled nucleon mass (2*AtomicMassUnit, GeV)
264  static const G4double mmu=.105658369; // Mass of a muon in GeV
265  static const G4double mmu2=mmu*mmu;   // Squared mass of a muon in GeV^2
266  static const G4double thresh=mmu+mmu2/dmN; // Universal threshold in GeV
267  // ---------
268  //static const G4double infEn = 9.e27;
269  G4double dN=0.;
270  if(Z>0||N>0) dN=thresh*GeV; // @@ if upgraded, change it in a total cross section
271  //@@ "dN=mmu+mmu2/G4NucleiProperties::GetNuclearMass(<G4double>(Z+N),<G4double>(Z)/GeV"
272  return dN;
273}
274
275// The main member function giving the gamma-A cross section (E_kin in MeV, CS in mb)
276G4double G4QNuMuNuclearCrossSection::CalculateCrossSection(G4bool CS, G4int F, G4int I,
277                                        G4int, G4int targZ, G4int targN, G4double Momentum)
278{
279  static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2
280  static const G4int nE=65;   // !! If change this, change it in GetFunctions() (*.hh) !!
281  static const G4int mL=nE-1;
282  static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
283  static const G4double dmN=mN+mN;   // Doubled nucleon mass (2*AtomicMassUnit, GeV)
284  static const G4double mmu=.105658369; // Mass of a muon in GeV
285  static const G4double mmu2=mmu*mmu;// Squared mass of a muon in GeV^2
286  static const G4double EMi=mmu+mmu2/dmN; // Universal threshold of the reaction in GeV
287  static const G4double EMa=300.;    // Maximum tabulated Energy of nu_mu in GeV
288  // *** Begin of the Associative memory for acceleration of the cross section calculations
289  static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional)
290  static G4bool first=true;          // Flag of initialization of the energy axis
291  // *** End of Static Definitions (Associative Memory) ***
292  //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
293  //G4double TotEnergy2=Momentum;
294  onlyCS=CS;                         // Flag to calculate only CS (not TX & QE)
295  lastE=Momentum/GeV;                // Kinetic energy of the muon neutrino (in GeV!)
296  if (lastE<=EMi)                    // Energy is below the minimum energy in the table
297  {
298    lastE=0.;
299    lastSig=0.;
300    return 0.;
301  }
302  G4int Z=targZ;                     // New Z, which can change the sign
303  if(F<=0)                           // This isotope was not the last used isotop
304  {
305    if(F<0)                          // This isotope was found in DAMDB =========> RETRIEVE
306    {
307      lastTX =(*TX)[I];              // Pointer to the prepared TX function (same isotope)
308      lastQE =(*QE)[I];              // Pointer to the prepared QE function (same isotope)
309   }
310   else                              // This isotope wasn't calculated previously => CREATE
311   {
312      if(first)
313      {
314        lastEN = new G4double[nE];   // This must be done only once!
315        Z=-Z;                        // To explain GetFunctions that E-axis must be filled
316        first=false;                 // To make it only once
317      }
318      lastTX = new G4double[nE];     // Allocate memory for the new TX function
319      lastQE = new G4double[nE];     // Allocate memory for the new QE function
320      G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK)
321      if(res<0) G4cerr<<"*W*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
322      // *** The synchronization check ***
323      G4int sync=TX->size();
324      if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
325      TX->push_back(lastTX);
326      QE->push_back(lastQE);
327    } // End of creation of the new set of parameters
328  } // End of parameters udate
329  // ============================== NOW Calculate the Cross Section =====================
330  if (lastE<=EMi)                   // Check that the neutrinoEnergy is higher than ThreshE
331  {
332    lastE=0.;
333    lastSig=0.;
334    return 0.;
335  }
336  if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization
337  {
338    G4int chk=1;
339    G4int ran=mL/2;
340    G4int sep=ran;  // as a result = an index of the left edge of the interval
341    while(ran>=2)
342    {
343      G4int newran=ran/2;
344      if(lastE<=lastEN[sep]) sep-=newran;
345      else                   sep+=newran;
346      ran=newran;
347      chk=chk+chk; 
348    }
349    if(chk+chk!=mL) G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
350    G4double lowE=lastEN[sep];
351    G4double highE=lastEN[sep+1];
352    G4double lowTX=lastTX[sep];
353    if(lastE<lowE||sep>=mL||lastE>highE)
354      G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
355            <<", sep="<<sep<<", mL="<<mL<<G4endl;
356    lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
357    if(!onlyCS)                       // Skip the differential cross-section parameters
358    {
359      G4double lowQE=lastQE[sep];
360      lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
361#ifdef pdebug
362      G4cout<<"G4NuMuNuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
363#endif
364    }
365  }
366  else
367  {
368    lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
369    lastQEL=lastQE[mL];
370  }
371  if(lastQEL<0.) lastQEL = 0.;
372  if(lastSig<0.) lastSig = 0.;
373  // The cross-sections are expected to be in mb
374  lastSig*=mb38;
375  if(!onlyCS) lastQEL*=mb38;
376  return lastSig;
377}
378
379// Calculate the cros-section functions
380// ****************************************************************************************
381// *** This tables are the same for all lepto-nuclear reactions, only mass is different ***
382// ***@@ IT'S REASONABLE TO MAKE ADDiTIONAL VIRTUAL CLASS FOR LEPTO-NUCLEAR WITH THIS@@ ***
383// ****************************************************************************************
384G4int G4QNuMuNuclearCrossSection::GetFunctions(G4int z, G4int n,
385                                                     G4double* t, G4double* q, G4double* e)
386{
387  static const G4double mN=.931494043;// Nucleon mass (inside nucleus, AtomicMassUnit, GeV)
388  static const G4double dmN=mN+mN;    // Doubled nucleon mass (2*AtomicMassUnit, GeV)
389  static const G4double mmu=.105658369; // Mass of a muon in GeV
390  static const G4double mmu2=mmu*mmu;   // Squared mass of a muon in GeV^2
391  static const G4double thresh=mmu+mmu2/dmN; // Universal threshold in GeV
392  static const G4int nE=65; // !! If change this, change it in GetCrossSection() (*.cc) !!
393  static const G4double nuEn[nE]={thresh,
394    .112039,.116079,.120416,.125076,.130090,.135494,.141324,.147626,.154445,.161838,
395    .169864,.178594,.188105,.198485,.209836,.222272,.235923,.250941,.267497,.285789,
396    .306045,.328530,.353552,.381466,.412689,.447710,.487101,.531538,.581820,.638893,
397    .703886,.778147,.863293,.961275,1.07445,1.20567,1.35843,1.53701,1.74667,1.99390,
398    2.28679,2.63542,3.05245,3.55386,4.15990,4.89644,5.79665,6.90336,8.27224,9.97606,
399    12.1106,14.8029,18.2223,22.5968,28.2351,35.5587,45.1481,57.8086,74.6682,97.3201,
400    128.036,170.085,228.220,309.420};
401  static const G4double TOTX[nE]={0.,
402    .108618,.352160,.476083,.566575,.639014,.699871,.752634,.799407,.841524,.879844,
403    .914908,.947050,.976456,1.00321,1.02734,1.04881,1.06755,1.08349,1.09653,1.10657,
404    1.11355,1.11739,1.11806,1.11556,1.10992,1.10124,1.08964,1.07532,1.05851,1.03950,
405    1.01859,.996169,.972593,.948454,.923773,.899081,.874713,.850965,.828082,.806265,
406    .785659,.766367,.748450,.731936,.716824,.703098,.690723,.679652,.669829,.661187,
407    .653306,.646682,.640986,.636125,.631993,.628479,.625458,.622800,.620364,.616231,
408    .614986,.612563,.609807,.606511};
409  static const G4double QELX[nE]={0.,
410    .012170,.040879,.057328,.070865,.083129,.094828,.106366,.118013,.129970,.142392,
411    .155410,.169138,.183676,.199123,.215573,.233120,.251860,.271891,.293317,.316246,
412    .340796,.367096,.395292,.425547,.458036,.491832,.524989,.556457,.585692,.612377,
413    .636544,.657790,.676260,.692007,.705323,.716105,.724694,.731347,.736340,.740172,
414    .742783,.744584,.745804,.746829,.747479,.747995,.748436,.749047,.749497,.749925,
415    .750486,.750902,.751268,.751566,.752026,.752266,.752428,.752761,.752873,.753094,
416    .753161,.753164,.753340,.753321};
417  // --------------------------------
418  G4int first=0;
419  if(z<0.)
420  {
421    first=1;
422    z=-z;
423  }
424  if(z<1 || z>92)             // neutron and plutonium are forbidden
425  {
426    G4cout<<"***G4QNuMuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
427    return -1;
428  }
429  for(G4int k=0; k<nE; k++)
430  {
431    G4double a=n+z;
432    G4double na=n+a;
433    G4double dn=n+n;
434    G4double da=a+a;
435    G4double ta=da+a;
436    if(first) e[k]=nuEn[k];       // Energy of neutrino E (first bin k=0 can be modified)
437    t[k]=TOTX[k]*nuEn[k]*(na+na)/ta+QELX[k]*(dn+dn-da)/ta; // TotalCrossSection
438    q[k]=QELX[k]*dn/a;                                     // QuasiElasticCrossSection
439  }
440  return first;
441}
442
443// Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic
444G4double G4QNuMuNuclearCrossSection::GetQEL_ExchangeQ2()
445{
446  static const G4double mmu=.105658369;// Mass of muon in GeV
447  static const G4double mmu2=mmu*mmu;  // Squared Mass of muon in GeV^2
448  static const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
449  static const double MN=.931494043;   // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
450  static const double MN2=MN*MN;       // M_N^2 in GeV^2
451  static const G4double power=-3.5;    // direct power for the magic variable
452  static const G4double pconv=1./power;// conversion power for the magic variable
453  static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
454  static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
455  static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
456  // Reversed table
457  static const G4double Xl[nQ2]={1.87905e-10,
458 .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374,
459 .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069,
460 .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718,
461 .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169,
462 .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670,
463 .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720,
464 .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457,
465 .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366,
466 .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128,
467 .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539};
468  // Direct table
469  static const G4double Xmax=Xl[lQ2];
470  static const G4double Xmin=Xl[0];
471  static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
472  static const G4double inl[nQ2]={0,
473 1.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811,
474 15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039,
475 26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876,
476 36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950,
477 46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377,
478 55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162,
479 65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852,
480 74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767,
481 83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106,
482 92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 100.000};
483  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
484  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
485  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
486  G4double ME=Enu*MN;                 // M*E
487  G4double dME=ME+ME;                 // 2*M*E
488  G4double dEMN=(dEnu+MN)*ME;
489  G4double MEm=ME-hmmu2;
490  G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
491  G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
492  G4double ymax=(E2M+sqE)/dEMN;
493  G4double ymin=(E2M-sqE)/dEMN;
494  G4double rmin=1.-ymin;
495  G4double rhm2E=hmmu2/Enu2;
496  G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
497  G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
498  G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
499  G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
500  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
501  G4double rXi=(Xmi-Xmin)/dX;
502  G4int    iXi=static_cast<int>(rXi);
503  if(iXi<0) iXi=0;
504  if(iXi>bQ2) iXi=bQ2;
505  G4double dXi=rXi-iXi;
506  G4double bnti=inl[iXi];
507  G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
508  //
509  G4double rXa=(Xma-Xmin)/dX;
510  G4int    iXa=static_cast<int>(rXa);
511  if(iXa<0) iXa=0;
512  if(iXa>bQ2) iXa=bQ2;
513  G4double dXa=rXa-iXa;
514  G4double bnta=inl[iXa];
515  G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
516  // *** Find X using the reversed table ***
517  G4double intx=inti+(inta-inti)*G4UniformRand();
518  G4int    intc=static_cast<int>(intx);
519  if(intc<0) intc=0;
520  if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
521  G4double dint=intx-intc;
522  G4double mX=Xl[intc];
523  G4double X=mX+dint*(Xl[intc+1]-mX);
524  G4double Q2=std::pow(X,pconv)-1.;
525  return Q2*GeV*GeV;
526}
527
528// Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic
529G4double G4QNuMuNuclearCrossSection::GetNQE_ExchangeQ2()
530{
531  static const double mpi=.13957018;    // charged pi meson mass in GeV
532  static const G4double mmu=.105658369; // Mass of muon in GeV
533  static const G4double mmu2=mmu*mmu;   // Squared Mass of muon in GeV^2
534  static const double hmmu2=mmu2/2;     // .5*m_mu^2 in GeV^2
535  static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
536  static const double MN2=MN*MN;        // M_N^2 in GeV^2
537  static const double dMN=MN+MN;        // 2*M_N in GeV
538  static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
539  static const G4int power=7;           // direct power for the magic variable
540  static const G4double pconv=1./power; // conversion power for the magic variable
541  static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
542  static const G4int lX=nX-1;           // index of the last in the Xl table
543  static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
544  static const G4int nE=20;             // #Of point in the El table (in GeV^2)
545  static const G4int bE=nE-1;           // index of the last in the El table
546  static const G4int pE=bE-1;           // index of the before last in the El table
547  // Reversed table
548  static const G4double X0[nX]={6.14081e-05,
549 .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244,
550 2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419};
551  static const G4double X1[nX]={.00125268,
552 .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049,
553 4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330};
554  static const G4double X2[nX]={.015694,
555 1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728,
556 10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276};
557  static const G4double X3[nX]={.0866877,
558 4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302,
559 21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115};
560  static const G4double X4[nX]={.160483,
561 5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127,
562 29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265};
563  static const G4double X5[nX]={.0999307,
564 5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398,
565 26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649};
566  static const G4double X6[nX]={.0276367,
567 3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519,
568 16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148};
569  static const G4double X7[nX]={.00472383,
570 2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473,
571 9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057};
572  static const G4double X8[nX]={.000630783,
573 1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480,
574 5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037};
575  static const G4double X9[nX]={7.49179e-05,
576 .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823,
577 3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028};
578  static const G4double XA[nX]={8.43437e-06,
579 .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704,
580 2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823};
581  static const G4double XB[nX]={9.27028e-07,
582 .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059,
583 1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341};
584  static const G4double XC[nX]={1.00807e-07,
585 .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668,
586 1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014};
587  static const G4double XD[nX]={1.09102e-08,
588 .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878,
589 1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617};
590  static const G4double XE[nX]={1.17831e-09,
591 .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374,
592 .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571};
593  static const G4double XF[nX]={1.27141e-10,
594 .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046,
595 .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338};
596  static const G4double XG[nX]={1.3713e-11,
597 .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870,
598 .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379};
599  static const G4double XH[nX]={1.47877e-12,
600 .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421,
601 .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197};
602  static const G4double XI[nX]={1.59454e-13,
603 .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806,
604 .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872};
605  static const G4double XJ[nX]={1.71931e-14,
606 .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147,
607 .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813};
608  // Direct table
609  static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
610                        X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
611  static const G4double dX[nE]={
612    (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
613    (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
614    (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
615    (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
616    (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
617  static const G4double* Xl[nE]=
618                             {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
619  static const G4double I0[nX]={0,
620 .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555,
621 13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000};
622  static const G4double I1[nX]={0,
623 .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365,
624 13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000};
625  static const G4double I2[nX]={0,
626 .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586,
627 12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000};
628  static const G4double I3[nX]={0,
629 .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461,
630 12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000};
631  static const G4double I4[nX]={0,
632 .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515,
633 11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000};
634  static const G4double I5[nX]={0,
635 .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040,
636 10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000};
637  static const G4double I6[nX]={0,
638 .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243,
639 9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000};
640  static const G4double I7[nX]={0,
641 .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717,
642 8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000};
643  static const G4double I8[nX]={0,
644 .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415,
645 8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000};
646  static const G4double I9[nX]={0,
647 .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672,
648 8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000};
649  static const G4double IA[nX]={0,
650 .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265,
651 7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000};
652  static const G4double IB[nX]={0,
653 .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298,
654 7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000};
655  static const G4double IC[nX]={0,
656 .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696,
657 7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000};
658  static const G4double ID[nX]={0,
659 .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212,
660 7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000};
661  static const G4double IE[nX]={0,
662 .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949,
663 7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000};
664  static const G4double IF[nX]={0,
665 .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164,
666 7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000};
667  static const G4double IG[nX]={0,
668 .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249,
669 7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000};
670  static const G4double IH[nX]={0,
671 .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804,
672 7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000};
673  static const G4double II[nX]={0,
674 .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722,
675 7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000};
676  static const G4double IJ[nX]={0,
677 .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196,
678 6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 20.0000};
679  static const G4double* Il[nE]=
680                             {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
681  static const G4double lE[nE]={
682-1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
683 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
684  static const G4double lEmi=lE[0];
685  static const G4double lEma=lE[nE-1];
686  static const G4double dlE=(lEma-lEmi)/bE;
687  //***************************************************************************************
688  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
689  G4double lEn=std::log(Enu);         // log(E) for interpolation
690  G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
691  G4int fE=static_cast<int>(rE);      // Left bin for interpolation
692  if(fE<0) fE=0;
693  if(fE>pE)fE=pE;
694  G4int    sE=fE+1;                   // Right bin for interpolation
695  G4double dE=rE-fE;                  // relative log shift from the left bin
696  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
697  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
698  G4double Emu=Enu-mmu;               // Free Energy of neutrino/anti-neutrino
699  G4double ME=Enu*MN;                 // M*E
700  G4double dME=ME+ME;                 // 2*M*E
701  G4double dEMN=(dEnu+MN)*ME;
702  G4double MEm=ME-hmmu2;
703  G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
704  G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
705  G4double ymax=(E2M+sqE)/dEMN;
706  G4double ymin=(E2M-sqE)/dEMN;
707  G4double rmin=1.-ymin;
708  G4double rhm2E=hmmu2/Enu2;
709  G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
710  G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
711  G4double Q2nq=Emu*dMN-mcV;
712  if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
713  // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
714  G4double Rmi=Q2mi/Q2ma;
715  G4double shift=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu
716  // --- E-interpolation must be done in a log scale ---
717  G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
718  G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
719  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
720  G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
721  G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
722  G4double rXi=(Xmi-iXmi)/idX;
723  G4int    iXi=static_cast<int>(rXi);
724  if(iXi<0) iXi=0;
725  if(iXi>bX) iXi=bX;
726  G4double dXi=rXi-iXi;
727  G4double bntil=Il[fE][iXi];
728  G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
729  G4double bntir=Il[sE][iXi];
730  G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
731  G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
732  //
733  G4double rXa=(Xma-iXmi)/idX;
734  G4int    iXa=static_cast<int>(rXa);
735  if(iXa<0) iXa=0;
736  if(iXa>bX) iXa=bX;
737  G4double dXa=rXa-iXa;
738  G4double bntal=Il[fE][iXa];
739  G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
740  G4double bntar=Il[sE][iXa];
741  G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
742  G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
743  //
744  // *** Find X using the reversed table ***
745  G4double intx=inti+(inta-inti)*G4UniformRand(); 
746  G4int    intc=static_cast<int>(intx);
747  if(intc<0) intc=0;
748  if(intc>bX) intc=bX;
749  G4double dint=intx-intc;
750  G4double mXl=Xl[fE][intc];
751  G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
752  G4double mXr=Xl[sE][intc];
753  G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
754  G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
755  G4double R=shift-std::pow(X,pconv);
756  G4double Q2=R*Q2ma;
757  return Q2*GeV*GeV;
758}
759
760// It returns a fraction of the direct interaction of the neutrino with quark-partons
761G4double G4QNuMuNuclearCrossSection::GetDirectPart(G4double Q2)
762{
763  G4double f=Q2/4.62;
764  G4double ff=f*f;
765  G4double r=ff*ff;
766  G4double s=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
767  //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
768  return 1.-s*(1.-s/2);
769}
770
771// #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
772G4double G4QNuMuNuclearCrossSection::GetNPartons(G4double Q2)
773{
774  return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
775}
776
777// This class can provide only virtual exchange pi+ (a substitute for W+ boson)
778G4int G4QNuMuNuclearCrossSection::GetExchangePDGCode() {return 211;}
Note: See TracBrowser for help on using the repository browser.