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

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

update ti head
File size: 38.9 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: G4QNuNuNuclearCrossSection.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: G4QNuNuNuclearCrossSection for (nu,nu)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 "G4QNuNuNuclearCrossSection.hh"
49
50// Initialization of the
51G4bool    G4QNuNuNuclearCrossSection::onlyCS=true;// Flag to calculate only CS (not QE)
52G4double  G4QNuNuNuclearCrossSection::lastSig=0.;// Last calculated total cross section
53G4double  G4QNuNuNuclearCrossSection::lastQEL=0.;// Last calculated quasi-el. cross section
54G4int     G4QNuNuNuclearCrossSection::lastL=0;   // Last used in cross section TheLastBin
55G4double  G4QNuNuNuclearCrossSection::lastE=0.;  // Last used in cross section TheEnergy
56G4double* G4QNuNuNuclearCrossSection::lastEN=0;  // Pointer to the Energy Scale of TX & QE
57G4double* G4QNuNuNuclearCrossSection::lastTX=0;  // Pointer to the LastArray of TX function
58G4double* G4QNuNuNuclearCrossSection::lastQE=0;  // Pointer to the LastArray of QE function
59G4int     G4QNuNuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile
60G4int     G4QNuNuNuclearCrossSection::lastN=0;   // The last N of calculated nucleus
61G4int     G4QNuNuNuclearCrossSection::lastZ=0;   // The last Z of calculated nucleus
62G4double  G4QNuNuNuclearCrossSection::lastP=0.;  // Last used in cross section Momentum
63G4double  G4QNuNuNuclearCrossSection::lastTH=0.; // Last threshold momentum
64G4double  G4QNuNuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section
65G4int     G4QNuNuNuclearCrossSection::lastI=0;   // The last position in the DAMDB
66std::vector<G4double*>* G4QNuNuNuclearCrossSection::TX = new std::vector<G4double*>;
67std::vector<G4double*>* G4QNuNuNuclearCrossSection::QE = new std::vector<G4double*>;
68
69// Returns Pointer to the G4VQCrossSection class
70G4VQCrossSection* G4QNuNuNuclearCrossSection::GetPointer()
71{
72  static G4QNuNuNuclearCrossSection theCrossSection; //**Static body of the Cross Section**
73  return &theCrossSection;
74}
75
76G4QNuNuNuclearCrossSection::~G4QNuNuNuclearCrossSection()
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 G4QNuNuNuclearCrossSection::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 G4QNuNuNuclearCrossSection::ThresholdEnergy(G4int, G4int, G4int) {return 0.;}
258
259// The main member function giving the gamma-A cross section (E_kin in MeV, CS in mb)
260G4double G4QNuNuNuclearCrossSection::CalculateCrossSection(G4bool CS, G4int F, G4int I,
261                                        G4int, G4int targZ, G4int targN, G4double Momentum)
262{
263  static const G4double mb38=1.E-11; // Conversion 10^-38 cm^2 to mb=10^-27 cm^2
264  static const G4int nE=65;   // !! If change this, change it in GetFunctions() (*.hh) !!
265  static const G4int mL=nE-1;
266  static const G4double EMi=0.;      // Universal threshold of the reaction in GeV
267  static const G4double EMa=300.;    // Maximum tabulated Energy of nu_mu in GeV
268  // *** Begin of the Associative memory for acceleration of the cross section calculations
269  static std::vector <G4double> colH;//?? Vector of HighEnergyCoefficients (functional)
270  static G4bool first=true;          // Flag of initialization of the energy axis
271  // *** End of Static Definitions (Associative Memory) ***
272  //const G4double Energy = aPart->GetKineticEnergy()/MeV; // Energy of the Muon
273  //G4double TotEnergy2=Momentum;
274  onlyCS=CS;                         // Flag to calculate only CS (not TX & QE)
275  lastE=Momentum/GeV;                // Kinetic energy of the muon neutrino (in GeV!)
276  if (lastE<=EMi)                    // Energy is below the minimum energy in the table
277  {
278    lastE=0.;
279    lastSig=0.;
280    return 0.;
281  }
282  G4int Z=targZ;                     // New Z, which can change the sign
283  if(F<=0)                           // This isotope was not the last used isotop
284  {
285    if(F<0)                          // This isotope was found in DAMDB =========> RETRIEVE
286    {
287      lastTX =(*TX)[I];              // Pointer to the prepared TX function (same isotope)
288      lastQE =(*QE)[I];              // Pointer to the prepared QE function (same isotope)
289   }
290   else                              // This isotope wasn't calculated previously => CREATE
291   {
292      if(first)
293      {
294        lastEN = new G4double[nE];   // This must be done only once!
295        Z=-Z;                        // To explain GetFunctions that E-axis must be filled
296        first=false;                 // To make it only once
297      }
298      lastTX = new G4double[nE];     // Allocate memory for the new TX function
299      lastQE = new G4double[nE];     // Allocate memory for the new QE function
300      G4int res=GetFunctions(Z,targN,lastTX,lastQE,lastEN);//@@analize(0=first,-1=bad,1=OK)
301      if(res<0) G4cerr<<"*W*G4NuMuNuclearCS::CalcCrossSect:Bad Function Retrieve"<<G4endl;
302      // *** The synchronization check ***
303      G4int sync=TX->size();
304      if(sync!=I) G4cerr<<"***G4NuMuNuclearCS::CalcCrossSect:Sync.="<<sync<<"#"<<I<<G4endl;
305      TX->push_back(lastTX);
306      QE->push_back(lastQE);
307    } // End of creation of the new set of parameters
308  } // End of parameters udate
309  // ============================== NOW Calculate the Cross Section =====================
310  if (lastE<=EMi)                   // Check that the neutrinoEnergy is higher than ThreshE
311  {
312    lastE=0.;
313    lastSig=0.;
314    return 0.;
315  }
316  if(lastE<EMa) // Linear fit is made explicitly to fix the last bin for the randomization
317  {
318    G4int chk=1;
319    G4int ran=mL/2;
320    G4int sep=ran;  // as a result = an index of the left edge of the interval
321    while(ran>=2)
322    {
323      G4int newran=ran/2;
324      if(lastE<=lastEN[sep]) sep-=newran;
325      else                   sep+=newran;
326      ran=newran;
327      chk=chk+chk; 
328    }
329    if(chk+chk!=mL) G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Table! mL="<<mL<<G4endl;
330    G4double lowE=lastEN[sep];
331    G4double highE=lastEN[sep+1];
332    G4double lowTX=lastTX[sep];
333    if(lastE<lowE||sep>=mL||lastE>highE)
334      G4cerr<<"*Warn*G4NuMuNuclearCS::CalcCS:Bin! "<<lowE<<" < "<<lastE<<" < "<<highE
335            <<", sep="<<sep<<", mL="<<mL<<G4endl;
336    lastSig=lastE*(lastE-lowE)*(lastTX[sep+1]-lowTX)/(highE-lowE)+lowTX; // Recover *E
337    if(!onlyCS)                       // Skip the differential cross-section parameters
338    {
339      G4double lowQE=lastQE[sep];
340      lastQEL=(lastE-lowE)*(lastQE[sep+1]-lowQE)/(highE-lowE)+lowQE;
341#ifdef pdebug
342      G4cout<<"G4NuMuNuclearCS::CalcCS: T="<<lastSig<<",Q="<<lastQEL<<",E="<<lastE<<G4endl;
343#endif
344    }
345  }
346  else
347  {
348    lastSig=lastTX[mL]; // @@ No extrapolation, just a const, while it looks shrinking...
349    lastQEL=lastQE[mL];
350  }
351  if(lastQEL<0.) lastQEL = 0.;
352  if(lastSig<0.) lastSig = 0.;
353  // The cross-sections are expected to be in mb
354  lastSig*=mb38;
355  if(!onlyCS) lastQEL*=mb38;
356  return lastSig;
357}
358
359// Calculate the cros-section functions
360// ****************************************************************************************
361// *** This tables are the same for all lepto-nuclear reactions, only mass is different ***
362// ***@@ IT'S REASONABLE TO MAKE ADDiTIONAL VIRTUAL CLASS FOR LEPTO-NUCLEAR WITH THIS@@ ***
363// ****************************************************************************************
364G4int G4QNuNuNuclearCrossSection::GetFunctions(G4int z, G4int n,
365                                                     G4double* t, G4double* q, G4double* e)
366{
367  static const G4int nE=65; // !! If change this, change it in GetCrossSection() (*.cc) !!
368  static const G4double nuEn[nE]={0.,
369  1.00463e-5,1.05336e-5,1.10692e-5,1.16592e-5,1.23109e-5,1.30323e-5,1.38331e-5,1.47245e-5,
370  1.57194e-5,1.68335e-5,1.80848e-5,1.94948e-5,2.10894e-5,2.28991e-5,2.49608e-5,2.73189e-5,
371  3.00273e-5,3.31516e-5,3.67722e-5,4.09881e-5,4.59217e-5,5.17255e-5,5.85908e-5,6.67583e-5,
372  7.65338e-5,8.83078e-5,.000102583,.000120011,.000141441,.000167995,.000201160,.000242926,
373  .000295985,.000364008,.000452051,.000567152,.000719210,.000922307,.001196710,.001571930,
374  .002091530,.002820590,.003857810,.005354930,.007548840,.010815300,.015760100,.023376900,
375  .035325600,.054430800,.085595700,.137508000,.225898000,.379892000,.654712000,1.15767000,
376  2.10277000,3.92843000,7.55861000,14.9991000,30.7412000,65.1734000,143.155000,326.326000};
377  static const G4double TOTX[nE]={0.,
378  3.18319e-5,3.33759e-5,3.50729e-5,3.69425e-5,3.90071e-5,4.12928e-5,4.38300e-5,4.66541e-5,
379  4.98065e-5,5.33360e-5,5.73004e-5,6.17678e-5,6.68196e-5,7.25529e-5,7.90845e-5,8.65551e-5,
380  9.51355e-5,.000105033,.000116503,.000129858,.000145486,.000163871,.000185616,.000211484,
381  .000242444,.000279730,.000324932,.000380110,.000447953,.000531999,.000636946,.000769076,
382  .000936873,.001151900,.001430050,.001793410,.002272980,.002912690,.003775490,.004952530,
383  .006577150,.008846400,.012054100,.017404800,.027443300,.040285600,.058251500,.084054600,
384  .121309000,.173273000,.240200000,.312865000,.367478000,.379071000,.348461000,.301143000,
385  .259802000,.231449000,.214007000,.203827000,.198265000,.195186000,.192821000,.190555000};
386  static const G4double QELX[nE]={0.,
387  .319793e-9,.351567e-9,.388228e-9,.430721e-9,.480212e-9,.538141e-9,.606306e-9,.686956e-9,
388  .782930e-9,.897832e-9,1.03626e-9,1.20415e-9,1.40918e-9,1.66139e-9,1.97401e-9,2.36458e-9,
389  2.85666e-9,3.48202e-9,4.28408e-9,5.32264e-9,6.68098e-9,8.47630e-9,1.08754e-8,1.41183e-8,
390  1.85552e-8,2.47023e-8,3.33325e-8,4.56172e-8,6.33590e-8,8.93735e-8,1.28128e-7,1.86829e-7,
391  2.77301e-7,4.19300e-7,6.46455e-7,1.01714e-6,1.63475e-6,2.68640e-6,4.51816e-6,7.78503e-6,
392  1.37563e-5,2.49521e-5,4.65025e-5,9.32015e-5,.000207165,.000435700,.000918046,.001964940,
393  .004285300,.009431390,.020560100,.043021600,.083012600,.142280000,.194239000,.222213000,
394  .231903000,.234030000,.234691000,.235204000,.235575000,.235793000,.235902000,.235946000};
395  // --------------------------------
396  G4int first=0;
397  if(z<0.)
398  {
399    first=1;
400    z=-z;
401  }
402  if(z<1 || z>92)             // neutron and plutonium are forbidden
403  {
404    G4cout<<"***G4QNuNuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl;
405    return -1;
406  }
407  for(G4int k=0; k<nE; k++)
408  {
409    G4double a=n+z;
410    G4double na=n+a;
411    G4double dn=n+n;
412    G4double da=a+a;
413    G4double ta=da+a;
414    if(first) e[k]=nuEn[k];       // Energy of neutrino E (first bin k=0 can be modified)
415    t[k]=TOTX[k]*nuEn[k]*(na+na)/ta+QELX[k]*(dn+dn-da)/ta; // TotalCrossSection
416    q[k]=QELX[k]*dn/a;                                     // QuasiElasticCrossSection
417  }
418  return first;
419}
420
421// Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic
422G4double G4QNuNuNuclearCrossSection::GetQEL_ExchangeQ2()
423{
424  static const double MN=.931494043;   // Nucleon mass (inside nucleus, atomicMassUnit,GeV)
425  static const G4double power=-3.5;    // direct power for the magic variable
426  static const G4double pconv=1./power;// conversion power for the magic variable
427  static const G4int nQ2=101;          // #Of point in the Q2l table (in GeV^2)
428  static const G4int lQ2=nQ2-1;        // index of the last in the Q2l table
429  static const G4int bQ2=lQ2-1;        // index of the before last in the Q2 ltable
430  // Reversed table
431  static const G4double Xl[nQ2]={1.87905e-10,
432 .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374,
433 .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069,
434 .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718,
435 .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169,
436 .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670,
437 .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720,
438 .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457,
439 .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366,
440 .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128,
441 .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539};
442  // Direct table
443  static const G4double Xmax=Xl[lQ2];
444  static const G4double Xmin=Xl[0];
445  static const G4double dX=(Xmax-Xmin)/lQ2;  // step in X(Q2, GeV^2)
446  static const G4double inl[nQ2]={0,
447 1.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811,
448 15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039,
449 26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876,
450 36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950,
451 46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377,
452 55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162,
453 65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852,
454 74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767,
455 83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106,
456 92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 100.000};
457  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
458  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
459  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
460  G4double ME=Enu*MN;                 // M*E
461  G4double dME=ME+ME;                 // 2*M*E
462  G4double dEMN=(dEnu+MN)*ME;
463  G4double sqE=Enu*ME;
464  G4double E2M=MN*Enu2;
465  G4double ymax=(E2M+sqE)/dEMN;
466  G4double Q2mi=0; // Q2_min(E_nu)
467  G4double Q2ma=dME*ymax;                                                  // Q2_max(E_nu)
468  G4double Xma=std::pow((1.+Q2mi),power);  // X_max(E_nu)
469  G4double Xmi=std::pow((1.+Q2ma),power);  // X_min(E_nu)
470  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
471  G4double rXi=(Xmi-Xmin)/dX;
472  G4int    iXi=static_cast<int>(rXi);
473  if(iXi<0) iXi=0;
474  if(iXi>bQ2) iXi=bQ2;
475  G4double dXi=rXi-iXi;
476  G4double bnti=inl[iXi];
477  G4double inti=bnti+dXi*(inl[iXi+1]-bnti);
478  //
479  G4double rXa=(Xma-Xmin)/dX;
480  G4int    iXa=static_cast<int>(rXa);
481  if(iXa<0) iXa=0;
482  if(iXa>bQ2) iXa=bQ2;
483  G4double dXa=rXa-iXa;
484  G4double bnta=inl[iXa];
485  G4double inta=bnta+dXa*(inl[iXa+1]-bnta);
486  // *** Find X using the reversed table ***
487  G4double intx=inti+(inta-inti)*G4UniformRand();
488  G4int    intc=static_cast<int>(intx);
489  if(intc<0) intc=0;
490  if(intc>bQ2) intc=bQ2;         // If it is more than max, then the BAD extrapolation
491  G4double dint=intx-intc;
492  G4double mX=Xl[intc];
493  G4double X=mX+dint*(Xl[intc+1]-mX);
494  G4double Q2=std::pow(X,pconv)-1.;
495  return Q2*GeV*GeV;
496}
497
498// Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic
499G4double G4QNuNuNuclearCrossSection::GetNQE_ExchangeQ2()
500{
501  static const double mpi=.13957018;    // charged pi meson mass in GeV
502  static const double MN=.931494043;    // Nucleon mass (inside nucleus,atomicMassUnit,GeV)
503  static const double dMN=MN+MN;        // 2*M_N in GeV
504  static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic
505  static const G4int power=7;           // direct power for the magic variable
506  static const G4double pconv=1./power; // conversion power for the magic variable
507  static const G4int nX=21;             // #Of point in the Xl table (in GeV^2)
508  static const G4int lX=nX-1;           // index of the last in the Xl table
509  static const G4int bX=lX-1;           // @@ index of the before last in the Xl table
510  static const G4int nE=20;             // #Of point in the El table (in GeV^2)
511  static const G4int bE=nE-1;           // index of the last in the El table
512  static const G4int pE=bE-1;           // index of the before last in the El table
513  // Reversed table
514  static const G4double X0[nX]={6.14081e-05,
515 .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244,
516 2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419};
517  static const G4double X1[nX]={.00125268,
518 .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049,
519 4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330};
520  static const G4double X2[nX]={.015694,
521 1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728,
522 10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276};
523  static const G4double X3[nX]={.0866877,
524 4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302,
525 21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115};
526  static const G4double X4[nX]={.160483,
527 5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127,
528 29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265};
529  static const G4double X5[nX]={.0999307,
530 5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398,
531 26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649};
532  static const G4double X6[nX]={.0276367,
533 3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519,
534 16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148};
535  static const G4double X7[nX]={.00472383,
536 2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473,
537 9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057};
538  static const G4double X8[nX]={.000630783,
539 1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480,
540 5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037};
541  static const G4double X9[nX]={7.49179e-05,
542 .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823,
543 3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028};
544  static const G4double XA[nX]={8.43437e-06,
545 .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704,
546 2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823};
547  static const G4double XB[nX]={9.27028e-07,
548 .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059,
549 1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341};
550  static const G4double XC[nX]={1.00807e-07,
551 .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668,
552 1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014};
553  static const G4double XD[nX]={1.09102e-08,
554 .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878,
555 1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617};
556  static const G4double XE[nX]={1.17831e-09,
557 .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374,
558 .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571};
559  static const G4double XF[nX]={1.27141e-10,
560 .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046,
561 .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338};
562  static const G4double XG[nX]={1.3713e-11,
563 .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870,
564 .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379};
565  static const G4double XH[nX]={1.47877e-12,
566 .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421,
567 .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197};
568  static const G4double XI[nX]={1.59454e-13,
569 .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806,
570 .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872};
571  static const G4double XJ[nX]={1.71931e-14,
572 .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147,
573 .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813};
574  // Direct table
575  static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0],
576                        X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]};
577  static const G4double dX[nE]={
578    (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX,
579    (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX,
580    (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX,
581    (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX,
582    (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX};
583  static const G4double* Xl[nE]=
584                             {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ};
585  static const G4double I0[nX]={0,
586 .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555,
587 13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000};
588  static const G4double I1[nX]={0,
589 .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365,
590 13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000};
591  static const G4double I2[nX]={0,
592 .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586,
593 12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000};
594  static const G4double I3[nX]={0,
595 .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461,
596 12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000};
597  static const G4double I4[nX]={0,
598 .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515,
599 11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000};
600  static const G4double I5[nX]={0,
601 .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040,
602 10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000};
603  static const G4double I6[nX]={0,
604 .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243,
605 9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000};
606  static const G4double I7[nX]={0,
607 .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717,
608 8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000};
609  static const G4double I8[nX]={0,
610 .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415,
611 8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000};
612  static const G4double I9[nX]={0,
613 .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672,
614 8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000};
615  static const G4double IA[nX]={0,
616 .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265,
617 7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000};
618  static const G4double IB[nX]={0,
619 .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298,
620 7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000};
621  static const G4double IC[nX]={0,
622 .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696,
623 7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000};
624  static const G4double ID[nX]={0,
625 .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212,
626 7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000};
627  static const G4double IE[nX]={0,
628 .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949,
629 7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000};
630  static const G4double IF[nX]={0,
631 .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164,
632 7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000};
633  static const G4double IG[nX]={0,
634 .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249,
635 7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000};
636  static const G4double IH[nX]={0,
637 .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804,
638 7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000};
639  static const G4double II[nX]={0,
640 .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722,
641 7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000};
642  static const G4double IJ[nX]={0,
643 .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196,
644 6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 20.0000};
645  static const G4double* Il[nE]=
646                             {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ};
647  static const G4double lE[nE]={
648-1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292,
649 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218};
650  static const G4double lEmi=lE[0];
651  static const G4double lEma=lE[nE-1];
652  static const G4double dlE=(lEma-lEmi)/bE;
653  //***************************************************************************************
654  G4double Enu=lastE;                 // Get energy of the last calculated cross-section
655  G4double lEn=std::log(Enu);         // log(E) for interpolation
656  G4double rE=(lEn-lEmi)/dlE;         // Position of the energy
657  G4int fE=static_cast<int>(rE);      // Left bin for interpolation
658  if(fE<0) fE=0;
659  if(fE>pE)fE=pE;
660  G4int    sE=fE+1;                   // Right bin for interpolation
661  G4double dE=rE-fE;                  // relative log shift from the left bin
662  G4double dEnu=Enu+Enu;              // doubled energy of nu/anu
663  G4double Enu2=Enu*Enu;              // squared energy of nu/anu
664  G4double Emu=Enu;                   // Free Energy of neutrino/anti-neutrino
665  G4double ME=Enu*MN;                 // M*E
666  G4double dME=ME+ME;                 // 2*M*E
667  G4double dEMN=(dEnu+MN)*ME;
668  G4double sqE=Enu*ME;
669  G4double E2M=MN*Enu2;
670  G4double ymax=(E2M+sqE)/dEMN;
671  G4double Q2mi=0.;                   // Q2_min(E_nu)
672  G4double Q2ma=dME*ymax;             // Q2_max(E_nu)
673  G4double Q2nq=Emu*dMN-mcV;
674  if(Q2ma>Q2nq) Q2ma=Q2nq;            // Correction for Non Quasi Elastic
675  // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma ---
676  G4double Rmi=Q2mi/Q2ma;
677  G4double shift=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu
678  // --- E-interpolation must be done in a log scale ---
679  G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu)
680  G4double Xma=std::pow((shift-1.),power); // X_max(E_nu)
681  // Find the integral values integ(Xmi) & integ(Xma) using the direct table
682  G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step
683  G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum
684  G4double rXi=(Xmi-iXmi)/idX;
685  G4int    iXi=static_cast<int>(rXi);
686  if(iXi<0) iXi=0;
687  if(iXi>bX) iXi=bX;
688  G4double dXi=rXi-iXi;
689  G4double bntil=Il[fE][iXi];
690  G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil);
691  G4double bntir=Il[sE][iXi];
692  G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir);
693  G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral
694  //
695  G4double rXa=(Xma-iXmi)/idX;
696  G4int    iXa=static_cast<int>(rXa);
697  if(iXa<0) iXa=0;
698  if(iXa>bX) iXa=bX;
699  G4double dXa=rXa-iXa;
700  G4double bntal=Il[fE][iXa];
701  G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal);
702  G4double bntar=Il[sE][iXa];
703  G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar);
704  G4double inta=intal+dE*(intar-intal);// interpolated end of the integral
705  //
706  // *** Find X using the reversed table ***
707  G4double intx=inti+(inta-inti)*G4UniformRand(); 
708  G4int    intc=static_cast<int>(intx);
709  if(intc<0) intc=0;
710  if(intc>bX) intc=bX;
711  G4double dint=intx-intc;
712  G4double mXl=Xl[fE][intc];
713  G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl);
714  G4double mXr=Xl[sE][intc];
715  G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr);
716  G4double X=Xlb+dE*(Xrb-Xlb);        // interpolated X value
717  G4double R=shift-std::pow(X,pconv);
718  G4double Q2=R*Q2ma;
719  return Q2*GeV*GeV;
720}
721
722// It returns a fraction of the direct interaction of the neutrino with quark-partons
723G4double G4QNuNuNuclearCrossSection::GetDirectPart(G4double Q2)
724{
725  G4double f=Q2/4.62;
726  G4double ff=f*f;
727  G4double r=ff*ff;
728  G4double s=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3)));
729  //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par)
730  return 1.-s*(1.-s/2);
731}
732
733// #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut
734G4double G4QNuNuNuclearCrossSection::GetNPartons(G4double Q2)
735{
736  return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space
737}
738
739// This class can provide only virtual exchange gamma (a substitute for Z0 boson)
740G4int G4QNuNuNuclearCrossSection::GetExchangePDGCode() {return 22;}
Note: See TracBrowser for help on using the repository browser.