source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QNuMuNuclearCrossSection.cc @ 1197

Last change on this file since 1197 was 1055, checked in by garnier, 15 years ago

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