source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/cross_sections/src/G4QANuMuNuclearCrossSection.cc @ 1340

Last change on this file since 1340 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: G4QANuMuNuclearCrossSection.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: G4QANuMuNuclearCrossSection for (anu_mu,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 "G4QANuMuNuclearCrossSection.hh"
49
50// Initialization of the
51G4bool    G4QANuMuNuclearCrossSection::onlyCS=true;//Flag to calculate only CS (not QE)
52G4double  G4QANuMuNuclearCrossSection::lastSig=0.;//Last calculated total cross section
53G4double  G4QANuMuNuclearCrossSection::lastQEL=0.;//Last calculated quasi-el. cross section
54G4int     G4QANuMuNuclearCrossSection::lastL=0;   //Last used in cross section TheLastBin
55G4double  G4QANuMuNuclearCrossSection::lastE=0.;  //Last used in cross section TheEnergy
56G4double* G4QANuMuNuclearCrossSection::lastEN=0;  //Pointer to the Energy Scale of TX & QE
57G4double* G4QANuMuNuclearCrossSection::lastTX=0;  //Pointer to the LastArray of TX function
58G4double* G4QANuMuNuclearCrossSection::lastQE=0;  //Pointer to the LastArray of QE function
59G4int     G4QANuMuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
60G4int     G4QANuMuNuclearCrossSection::lastN=0;   // The last N of calculated nucleus
61G4int     G4QANuMuNuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
62G4double  G4QANuMuNuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
63G4double  G4QANuMuNuclearCrossSection::lastTH=0.; // Last threshold momentum
64G4double  G4QANuMuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
65G4int     G4QANuMuNuclearCrossSection::lastI=0;   // The last position in the DAMDB
66std::vector<G4double*>* G4QANuMuNuclearCrossSection::TX = new std::vector<G4double*>;
67std::vector<G4double*>* G4QANuMuNuclearCrossSection::QE = new std::vector<G4double*>;
68
69// Returns Pointer to the G4VQCrossSection class
70G4VQCrossSection* G4QANuMuNuclearCrossSection::GetPointer()
71{
72  static G4QANuMuNuclearCrossSection theCrossSection;//**Static body of the Cross Section**
73  return &theCrossSection;
74}
75
76G4QANuMuNuclearCrossSection::~G4QANuMuNuclearCrossSection()
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 G4QANuMuNuclearCrossSection::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<<"G4QAMNCS::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 debug
109    G4cout<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"---G4QAMNCrossSec::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<<"G4QAMNCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lstI="<<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<<"G4QAMNCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
190#endif
191        if(pEn>lastTH)
192        {
193#ifdef pdebug
194          G4cout<<"G4QAMNCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
195#endif
196          lastTH=pEn;
197        }
198      }
199#ifdef pdebug
200      G4cout<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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<<"G4QAMNCS::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 G4QANuMuNuclearCrossSection::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 G4QANuMuNuclearCrossSection::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 antiNuEnergy 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 G4QANuMuNuclearCrossSection::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    .077498,.247583,.329691,.386384,.429087,.462699,.489899,.512316,.530996,.546614,
403    .559616,.570292,.578840,.585395,.590053,.593083,.594197,.593614,.591396,.587611,
404    .582335,.575653,.567667,.558490,.548417,.537270,.525352,.512825,.499857,.486620,
405    .473283,.460014,.446970,.434294,.422116,.410656,.399782,.389665,.380349,.371860,
406    .364207,.357387,.351388,.346192,.341778,.338122,.335198,.332980,.331439,.330544,
407    .330263,.330558,.331391,.332718,.334494,.336667,.339182,.341697,.344470,.348125,
408    .351322,.354481,.357507,.359239};
409  static const G4double QELX[nE]={0.,
410    .008683,.028739,.039700,.048327,.055820,.062693,.069235,.075631,.082010,.088463,
411    .095059,.101851,.108883,.116192,.123814,.131826,.140185,.148962,.158197,.167933,
412    .178221,.189119,.200700,.213045,.226326,.240454,.255277,.270612,.286388,.302608,
413    .319318,.336582,.354468,.373031,.392427,.412445,.433146,.454448,.476222,.498289,
414    .520430,.542558,.564130,.585003,.604928,.623680,.641266,.657255,.671704,.684586,
415    .696111,.706028,.714553,.721951,.728085,.733182,.737348,.740958,.743716,.746059,
416    .747806,.749129,.750331,.751100};
417
418  // --------------------------------
419  G4int first=0;
420  if(z<0.)
421  {
422    first=1;
423    z=-z;
424  }
425  if(z<1 || z>92)             // neutron & plutonium are forbidden
426  {
427    G4cout<<"**G4QANuMuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
428    return -1;
429  }
430  for(G4int k=0; k<nE; k++)
431  {
432    G4double a=n+z;
433    G4double za=z+a;
434    G4double dz=z+z;
435    G4double da=a+a;
436    G4double ta=da+a;
437    if(first) e[k]=nuEn[k];       // Energy of neutrino E (first bin k=0 can be modified)
438    t[k]=TOTX[k]*nuEn[k]*(za+za)/ta+QELX[k]*(dz+dz-da)/ta; // TotalCrossSection
439    q[k]=QELX[k]*dz/a;                                     // QuasiElasticCrossSection
440  }
441  return first;
442}
443
444// Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic
445G4double G4QANuMuNuclearCrossSection::GetQEL_ExchangeQ2()
446{
447  static const G4double mmu=.105658369;// Mass of muon in GeV
448  static const G4double mmu2=mmu*mmu;  // Squared Mass of muon in GeV^2
449  static const double hmmu2=mmu2/2;    // .5*m_mu^2 in GeV^2
450  static const double MN=.931494043;   // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
451  static const double MN2=MN*MN;       // M_N^2 in GeV^2
452  static const G4double power=-3.5;    // direct power for the magic variable
453  static const G4double pconv=1./power;// conversion power for the magic variable
454  static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
455  static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
456  static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
457  // Reversed table
458  static const G4double Xl[nQ2]={5.20224e-16,
459 .006125,.0137008,.0218166,.0302652,.0389497,.0478144,.0568228,.0659497,.0751768,.0844898,
460 .093878, .103332, .112844, .122410, .132023, .141680, .151376, .161109, .170875, .180672,
461 .190499, .200352, .210230, .220131, .230055, .239999, .249963, .259945, .269944, .279960,
462 .289992, .300039, .310099, .320173, .330260, .340359, .350470, .360592, .370724, .380867,
463 .391019, .401181, .411352, .421531, .431719, .441915, .452118, .462329, .472547, .482771,
464 .493003, .503240, .513484, .523734, .533989, .544250, .554517, .564788, .575065, .585346,
465 .595632, .605923, .616218, .626517, .636820, .647127, .657438, .667753, .678072, .688394,
466 .698719, .709048, .719380, .729715, .740053, .750394, .760738, .771085, .781434, .791786,
467 .802140, .812497, .822857, .833219, .843582, .853949, .864317, .874687, .885060, .895434,
468 .905810, .916188, .926568, .936950, .947333, .957719, .968105, .978493, .988883, .999275};
469   // Direct table
470  static const G4double Xmax=Xl[lQ2];
471  static const G4double Xmin=Xl[0];
472  static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
473  static const G4double inl[nQ2]={0,
474 1.52225, 2.77846, 3.96651, 5.11612, 6.23990, 7.34467, 8.43466, 9.51272, 10.5809, 11.6406,
475 12.6932, 13.7394, 14.7801, 15.8158, 16.8471, 17.8743, 18.8979, 19.9181, 20.9353, 21.9496,
476 22.9614, 23.9707, 24.9777, 25.9826, 26.9855, 27.9866, 28.9860, 29.9837, 30.9798, 31.9745,
477 32.9678, 33.9598, 34.9505, 35.9400, 36.9284, 37.9158, 38.9021, 39.8874, 40.8718, 41.8553,
478 42.8379, 43.8197, 44.8007, 45.7810, 46.7605, 47.7393, 48.7174, 49.6950, 50.6718, 51.6481,
479 52.6238, 53.5990, 54.5736, 55.5476, 56.5212, 57.4943, 58.4670, 59.4391, 60.4109, 61.3822,
480 62.3531, 63.3236, 64.2937, 65.2635, 66.2329, 67.2019, 68.1707, 69.1390, 70.1071, 71.0748,
481 72.0423, 73.0095, 73.9763, 74.9429, 75.9093, 76.8754, 77.8412, 78.8068, 79.7721, 80.7373,
482 81.7022, 82.6668, 83.6313, 84.5956, 85.5596, 86.5235, 87.4872, 88.4507, 89.4140, 90.3771,
483 91.3401, 92.3029, 93.2656, 94.2281, 95.1904, 96.1526, 97.1147, 98.0766, 99.0384, 100.000};
484  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
485  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
486  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
487  G4double ME=Enu*MN;                 // M*E
488  G4double dME=ME+ME;                 // 2*M*E
489  G4double dEMN=(dEnu+MN)*ME;
490  G4double MEm=ME-hmmu2;
491  G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
492  G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
493  G4double ymax=(E2M+sqE)/dEMN;
494  G4double ymin=(E2M-sqE)/dEMN;
495  G4double rmin=1.-ymin;
496  G4double rhm2E=hmmu2/Enu2;
497  G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
498  G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
499  G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
500  G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
501  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
502  G4double rXi=(Xmi-Xmin)/dX;
503  G4int    iXi=static_cast<int>(rXi);
504  if(iXi<0) iXi=0;
505  if(iXi>bQ2) iXi=bQ2;
506  G4double dXi=rXi-iXi;
507  G4double bnti=inl[iXi];
508  G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
509  //
510  G4double rXa=(Xma-Xmin)/dX;
511  G4int    iXa=static_cast<int>(rXa);
512  if(iXa<0) iXa=0;
513  if(iXa>bQ2) iXa=bQ2;
514  G4double dXa=rXa-iXa;
515  G4double bnta=inl[iXa];
516  G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
517  // *** Find X using the reversed table ***
518  G4double intx=inti+(inta-inti)*G4UniformRand();
519  G4int    intc=static_cast<int>(intx);
520  if(intc<0) intc=0;
521  if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
522  G4double dint=intx-intc;
523  G4double mX=Xl[intc];
524  G4double X=mX+dint*(Xl[intc+1]-mX);
525  G4double Q2=std::pow(X,pconv)-1.;
526  return Q2*GeV*GeV;
527}
528
529// Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic
530G4double G4QANuMuNuclearCrossSection::GetNQE_ExchangeQ2()
531{
532  static const double mpi=.13957018;    // charged pi meson mass in GeV
533  static const G4double mmu=.105658369; // Mass of muon in GeV
534  static const G4double mmu2=mmu*mmu;   // Squared Mass of muon in GeV^2
535  static const double hmmu2=mmu2/2;     // .5*m_mu^2 in GeV^2
536  static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
537  static const double MN2=MN*MN;        // M_N^2 in GeV^2
538  static const double dMN=MN+MN;        // 2*M_N in GeV
539  static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
540  static const G4int power=7;           // direct power for the magic variable
541  static const G4double pconv=1./power; // conversion power for the magic variable
542  static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
543  static const G4int lX=nX-1;           // index of the last in the Xl table
544  static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
545  static const G4int nE=20;             // #Of point in the El table (in GeV^2)
546  static const G4int bE=nE-1;           // index of the last in the El table
547  static const G4int pE=bE-1;           // index of the before last in the El table
548  // Reversed table
549  static const G4double X0[nX]={5.21412e-05,
550 .437860, .681908, .891529, 1.08434, 1.26751, 1.44494, 1.61915, 1.79198, 1.96493, 2.13937,
551 2.31664, 2.49816, 2.68559, 2.88097, 3.08705, 3.30774, 3.54917, 3.82233, 4.15131, 4.62182};
552  static const G4double X1[nX]={.00102591,
553 1.00443, 1.55828, 2.03126, 2.46406, 2.87311, 3.26723, 3.65199, 4.03134, 4.40835, 4.78561,
554 5.16549, 5.55031, 5.94252, 6.34484, 6.76049, 7.19349, 7.64917, 8.13502, 8.66246, 9.25086};
555  static const G4double X2[nX]={.0120304,
556 2.59903, 3.98637, 5.15131, 6.20159, 7.18024, 8.10986, 9.00426, 9.87265, 10.7217, 11.5564,
557 12.3808, 13.1983, 14.0116, 14.8234, 15.6359, 16.4515, 17.2723, 18.1006, 18.9386, 19.7892};
558  static const G4double X3[nX]={.060124,
559 5.73857, 8.62595, 10.9849, 13.0644, 14.9636, 16.7340, 18.4066, 20.0019, 21.5342, 23.0142,
560 24.4497, 25.8471, 27.2114, 28.5467, 29.8564, 31.1434, 32.4102, 33.6589, 34.8912, 36.1095};
561  static const G4double X4[nX]={.0992363,
562 8.23746, 12.1036, 15.1740, 17.8231, 20.1992, 22.3792, 24.4092, 26.3198, 28.1320, 29.8615,
563 31.5200, 33.1169, 34.6594, 36.1536, 37.6044, 39.0160, 40.3920, 41.7353, 43.0485, 44.3354};
564  static const G4double X5[nX]={.0561127,
565 7.33661, 10.5694, 13.0778, 15.2061, 17.0893, 18.7973, 20.3717, 21.8400, 23.2211, 24.5291,
566 25.7745, 26.9655, 28.1087, 29.2094, 30.2721, 31.3003, 32.2972, 33.2656, 34.2076, 35.1265};
567  static const G4double X6[nX]={.0145859,
568 4.81774, 6.83565, 8.37399, 9.66291, 10.7920, 11.8075, 12.7366, 13.5975, 14.4025, 15.1608,
569 15.8791, 16.5628, 17.2162, 17.8427, 18.4451, 19.0259, 19.5869, 20.1300, 20.6566, 21.1706};
570  static const G4double X7[nX]={.00241155,
571 2.87095, 4.02492, 4.89243, 5.61207, 6.23747, 6.79613, 7.30433, 7.77270, 8.20858, 8.61732,
572 9.00296, 9.36863, 9.71682, 10.0495, 10.3684, 10.6749, 10.9701, 11.2550, 11.5306, 11.7982};
573  static const G4double X8[nX]={.000316863,
574 1.76189, 2.44632, 2.95477, 3.37292, 3.73378, 4.05420, 4.34415, 4.61009, 4.85651, 5.08666,
575 5.30299, 5.50738, 5.70134, 5.88609, 6.06262, 6.23178, 6.39425, 6.55065, 6.70149, 6.84742};
576  static const G4double X9[nX]={3.73544e-05,
577 1.17106, 1.61289, 1.93763, 2.20259, 2.42976, 2.63034, 2.81094, 2.97582, 3.12796, 3.26949,
578 3.40202, 3.52680, 3.64482, 3.75687, 3.86360, 3.96557, 4.06323, 4.15697, 4.24713, 4.33413};
579  static const G4double XA[nX]={4.19131e-06,
580 .849573, 1.16208, 1.38955, 1.57379, 1.73079, 1.86867, 1.99221, 2.10451, 2.20770, 2.30332,
581 2.39252, 2.47622, 2.55511, 2.62977, 2.70066, 2.76818, 2.83265, 2.89437, 2.95355, 3.01051};
582  static const G4double XB[nX]={4.59981e-07,
583 .666131, .905836, 1.07880, 1.21796, 1.33587, 1.43890, 1.53080, 1.61399, 1.69011, 1.76040,
584 1.82573, 1.88682, 1.94421, 1.99834, 2.04959, 2.09824, 2.14457, 2.18878, 2.23107, 2.27162};
585  static const G4double XC[nX]={4.99861e-08,
586 .556280, .752730, .893387, 1.00587, 1.10070, 1.18317, 1.25643, 1.32247, 1.38269, 1.43809,
587 1.48941, 1.53724, 1.58203, 1.62416, 1.66391, 1.70155, 1.73728, 1.77128, 1.80371, 1.83473};
588  static const G4double XD[nX]={5.40832e-09,
589 .488069, .657650, .778236, .874148, .954621, 1.02432, 1.08599, 1.14138, 1.19172, 1.23787,
590 1.28049, 1.32008, 1.35705, 1.39172, 1.42434, 1.45514, 1.48429, 1.51197, 1.53829, 1.56339};
591  static const G4double XE[nX]={5.84029e-10,
592 .445057, .597434, .705099, .790298, .861468, .922865, .976982, 1.02542, 1.06930, 1.10939,
593 1.14630, 1.18050, 1.21233, 1.24208, 1.27001, 1.29630, 1.32113, 1.34462, 1.36691, 1.38812};
594  static const G4double XF[nX]={6.30137e-11,
595 .418735, .560003, .659168, .737230, .802138, .857898, .906854, .950515, .989915, 1.02580,
596 1.05873, 1.08913, 1.11734, 1.14364, 1.16824, 1.19133, 1.21306, 1.23358, 1.25298, 1.27139};
597  static const G4double XG[nX]={6.79627e-12,
598 .405286, .539651, .633227, .706417, .766929, .818642, .863824, .903931, .939963, .972639,
599 1.00250, 1.02995, 1.05532, 1.07887, 1.10082, 1.12134, 1.14058, 1.15867, 1.17572, 1.19183};
600  static const G4double XH[nX]={7.32882e-13,
601 .404391, .535199, .625259, .695036, .752243, .800752, .842823, .879906, .912994, .942802,
602 .969862, .994583, 1.01729, 1.03823, 1.05763, 1.07566, 1.09246, 1.10816, 1.12286, 1.13667};
603  static const G4double XI[nX]={7.90251e-14,
604 .418084, .548382, .636489, .703728, .758106, .803630, .842633, .876608, .906576, .933269,
605 .957233, .978886, .998556, 1.01651, 1.03295, 1.04807, 1.06201, 1.07489, 1.08683, 1.09792};
606  static const G4double XJ[nX]={8.52083e-15,
607 .447299, .579635, .666780, .731788, .783268, .825512, .861013, .891356, .917626, .940597,
608 .960842, .978802, .994820, 1.00917, 1.02208, 1.03373, 1.04427, 1.05383, 1.06253, 1.07046};
609  // Direct table
610  static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
611                        X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
612  static const G4double dX[nE]={
613    (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
614    (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
615    (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
616    (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
617    (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
618  static const G4double* Xl[nE]=
619                             {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
620  static const G4double I0[nX]={0,
621 .354631, 1.08972, 2.05138, 3.16564, 4.38343, 5.66828, 6.99127, 8.32858, 9.65998, 10.9680,
622 12.2371, 13.4536, 14.6050, 15.6802, 16.6686, 17.5609, 18.3482, 19.0221, 19.5752, 20.0000};
623  static const G4double I1[nX]={0,
624 .281625, .877354, 1.67084, 2.60566, 3.64420, 4.75838, 5.92589, 7.12829, 8.34989, 9.57708,
625 10.7978, 12.0014, 13.1781, 14.3190, 15.4162, 16.4620, 17.4496, 18.3724, 19.2245, 20.0000};
626  static const G4double I2[nX]={0,
627 .201909, .642991, 1.24946, 1.98463, 2.82370, 3.74802, 4.74263, 5.79509, 6.89474, 8.03228,
628 9.19947, 10.3889, 11.5938, 12.8082, 14.0262, 15.2427, 16.4527, 17.6518, 18.8356, 20.0000};
629  static const G4double I3[nX]={0,
630 .140937, .461189, .920216, 1.49706, 2.17728, 2.94985, 3.80580, 4.73758, 5.73867, 6.80331,
631 7.92637, 9.10316, 10.3294, 11.6013, 12.9150, 14.2672, 15.6548, 17.0746, 18.5239, 20.0000};
632  static const G4double I4[nX]={0,
633 .099161, .337358, .694560, 1.16037, 1.72761, 2.39078, 3.14540, 3.98768, 4.91433, 5.92245,
634 7.00942, 8.17287, 9.41060, 10.7206, 12.1010, 13.5500, 15.0659, 16.6472, 18.2924, 20.0000};
635  static const G4double I5[nX]={0,
636 .071131, .255084, .543312, .932025, 1.41892, 2.00243, 2.68144, 3.45512, 4.32283, 5.28411,
637 6.33859, 7.48602, 8.72621, 10.0590, 11.4844, 13.0023, 14.6128, 16.3158, 18.1115, 20.0000};
638  static const G4double I6[nX]={0,
639 .053692, .202354, .443946, .778765, 1.20774, 1.73208, 2.35319, 3.07256, 3.89177, 4.81249,
640 5.83641, 6.96528, 8.20092, 9.54516, 10.9999, 12.5670, 14.2486, 16.0466, 17.9630, 20.0000};
641  static const G4double I7[nX]={0,
642 .043065, .168099, .376879, .672273, 1.05738, 1.53543, 2.10973, 2.78364, 3.56065, 4.44429,
643 5.43819, 6.54610, 7.77186, 9.11940, 10.5928, 12.1963, 13.9342, 15.8110, 17.8313, 20.0000};
644  static const G4double I8[nX]={0,
645 .036051, .143997, .327877, .592202, .941572, 1.38068, 1.91433, 2.54746, 3.28517, 4.13277,
646 5.09574, 6.17984, 7.39106, 8.73568, 10.2203, 11.8519, 13.6377, 15.5854, 17.7033, 20.0000};
647  static const G4double I9[nX]={0,
648 .030977, .125727, .289605, .528146, .846967, 1.25183, 1.74871, 2.34384, 3.04376, 3.85535,
649 4.78594, 5.84329, 7.03567, 8.37194, 9.86163, 11.5150, 13.3430, 15.3576, 17.5719, 20.0000};
650  static const G4double IA[nX]={0,
651 .027129, .111420, .258935, .475812, .768320, 1.14297, 1.60661, 2.16648, 2.83034, 3.60650,
652 4.50394, 5.53238, 6.70244, 8.02569, 9.51488, 11.1841, 13.0488, 15.1264, 17.4362, 20.0000};
653  static const G4double IB[nX]={0,
654 .024170, .100153, .234345, .433198, .703363, 1.05184, 1.48607, 2.01409, 2.64459, 3.38708,
655 4.25198, 5.25084, 6.39647, 7.70319, 9.18708, 10.8663, 12.7617, 14.8968, 17.2990, 20.0000};
656  static const G4double IC[nX]={0,
657 .021877, .091263, .214670, .398677, .650133, .976322, 1.38510, 1.88504, 2.48555, 3.19709,
658 4.03129, 5.00127, 6.12184, 7.40989, 8.88482, 10.5690, 12.4888, 14.6748, 17.1638, 20.0000};
659  static const G4double ID[nX]={0,
660 .020062, .084127, .198702, .370384, .606100, .913288, 1.30006, 1.77535, 2.34912, 3.03253,
661 3.83822, 4.78063, 5.87634, 7.14459, 8.60791, 10.2929, 12.2315, 14.4621, 17.0320, 20.0000};
662  static const G4double IE[nX]={0,
663 .018547, .078104, .185102, .346090, .567998, .858331, 1.22535, 1.67824, 2.22735, 2.88443,
664 3.66294, 4.57845, 5.64911, 6.89637, 8.34578, 10.0282, 11.9812, 14.2519, 16.8993, 20.0000};
665  static const G4double IF[nX]={0,
666 .017143, .072466, .172271, .323007, .531545, .805393, 1.15288, 1.58338, 2.10754, 2.73758,
667 3.48769, 4.37450, 5.41770, 6.64092, 8.07288, 9.74894, 11.7135, 14.0232, 16.7522, 20.0000};
668  static const G4double IG[nX]={0,
669 .015618, .066285, .158094, .297316, .490692, .745653, 1.07053, 1.47479, 1.96931, 2.56677,
670 3.28205, 4.13289, 5.14068, 6.33158, 7.73808, 9.40133, 11.3745, 13.7279, 16.5577, 20.0000};
671  static const G4double IH[nX]={0,
672 .013702, .058434, .139923, .264115, .437466, .667179, .961433, 1.32965, 1.78283, 2.33399,
673 2.99871, 3.79596, 4.74916, 5.88771, 7.24937, 8.88367, 10.8576, 13.2646, 16.2417, 20.0000};
674  static const G4double II[nX]={0,
675 .011264, .048311, .116235, .220381, .366634, .561656, .813132, 1.13008, 1.52322, 2.00554,
676 2.59296, 3.30542, 4.16834, 5.21490, 6.48964, 8.05434, 9.99835, 12.4580, 15.6567, 20.0000};
677  static const G4double IJ[nX]={0,
678 .008628, .037206, .089928, .171242, .286114, .440251, .640343, .894382, 1.21208, 1.60544,
679 2.08962, 2.68414, 3.41486, 4.31700, 5.44048, 6.85936, 8.69067, 11.1358, 14.5885, 20.0000};
680  static const G4double* Il[nE]=
681                             {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
682  static const G4double lE[nE]={
683-1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
684 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
685  static const G4double lEmi=lE[0];
686  static const G4double lEma=lE[nE-1];
687  static const G4double dlE=(lEma-lEmi)/bE;
688  //***************************************************************************************
689  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
690  G4double lEn=std::log(Enu);         // log(E) for interpolation
691  G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
692  G4int fE=static_cast<int>(rE);      // Left bin for interpolation
693  if(fE<0) fE=0;
694  if(fE>pE)fE=pE;
695  G4int    sE=fE+1;                   // Right bin for interpolation
696  G4double dE=rE-fE;                  // relative log shift from the left bin
697  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
698  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
699  G4double Emu=Enu-mmu;               // Free Energy of neutrino/anti-neutrino
700  G4double ME=Enu*MN;                 // M*E
701  G4double dME=ME+ME;                 // 2*M*E
702  G4double dEMN=(dEnu+MN)*ME;
703  G4double MEm=ME-hmmu2;
704  G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2);
705  G4double E2M=MN*Enu2-(Enu+MN)*hmmu2;
706  G4double ymax=(E2M+sqE)/dEMN;
707  G4double ymin=(E2M-sqE)/dEMN;
708  G4double rmin=1.-ymin;
709  G4double rhm2E=hmmu2/Enu2;
710  G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu)
711  G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
712  G4double Q2nq=Emu*dMN-mcV;
713  if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
714  // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
715  G4double Rmi=Q2mi/Q2ma;
716  G4double shift=.875/(1.+.2977/Enu/Enu)/std::pow(Enu,.78);
717  // --- E-interpolation must be done in a log scale ---
718  G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
719  G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
720  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
721  G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
722  G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
723  G4double rXi=(Xmi-iXmi)/idX;
724  G4int    iXi=static_cast<int>(rXi);
725  if(iXi<0) iXi=0;
726  if(iXi>bX) iXi=bX;
727  G4double dXi=rXi-iXi;
728  G4double bntil=Il[fE][iXi];
729  G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
730  G4double bntir=Il[sE][iXi];
731  G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
732  G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
733  //
734  G4double rXa=(Xma-iXmi)/idX;
735  G4int    iXa=static_cast<int>(rXa);
736  if(iXa<0) iXa=0;
737  if(iXa>bX) iXa=bX;
738  G4double dXa=rXa-iXa;
739  G4double bntal=Il[fE][iXa];
740  G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
741  G4double bntar=Il[sE][iXa];
742  G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
743  G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
744  //
745  // *** Find X using the reversed table ***
746  G4double intx=inti+(inta-inti)*G4UniformRand(); 
747  G4int    intc=static_cast<int>(intx);
748  if(intc<0) intc=0;
749  if(intc>bX) intc=bX;
750  G4double dint=intx-intc;
751  G4double mXl=Xl[fE][intc];
752  G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
753  G4double mXr=Xl[sE][intc];
754  G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
755  G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
756  G4double R=shift-std::pow(X,pconv);
757  G4double Q2=R*Q2ma;
758  return Q2*GeV*GeV;
759}
760
761// It returns a fraction of the direct interaction of the neutrino with quark-partons
762G4double G4QANuMuNuclearCrossSection::GetDirectPart(G4double Q2)
763{
764  G4double f=Q2/4.62;
765  G4double ff=f*f;
766  G4double r=ff*ff;
767  G4double s=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
768  //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
769  return 1.-s*(1.-s/2);
770}
771
772// #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
773G4double G4QANuMuNuclearCrossSection::GetNPartons(G4double Q2)
774{
775  return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
776}
777
778// This class can provide only virtual exchange pi- (a substitute for W- boson)
779G4int G4QANuMuNuclearCrossSection::GetExchangePDGCode() {return -211;}
Note: See TracBrowser for help on using the repository browser.