[819] | 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.12 2007/11/01 16:09:38 mkossov Exp $ |
---|
[962] | 28 | // GEANT4 tag $Name: geant4-09-02-ref-02 $ |
---|
[819] | 29 | // |
---|
| 30 | // |
---|
| 31 | // G4 Physics class: G4QNuMuNuclearCrossSection for gamma+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 |
---|
| 51 | G4bool G4QNuMuNuclearCrossSection::onlyCS=true;// Flag to calculate only CS (not QE) |
---|
| 52 | G4double G4QNuMuNuclearCrossSection::lastSig=0.;// Last calculated total cross section |
---|
| 53 | G4double G4QNuMuNuclearCrossSection::lastQEL=0.;// Last calculated quasi-el. cross section |
---|
| 54 | G4int G4QNuMuNuclearCrossSection::lastL=0; // Last used in cross section TheLastBin |
---|
| 55 | G4double G4QNuMuNuclearCrossSection::lastE=0.; // Last used in cross section TheEnergy |
---|
| 56 | G4double* G4QNuMuNuclearCrossSection::lastEN=0; // Pointer to the Energy Scale of TX & QE |
---|
| 57 | G4double* G4QNuMuNuclearCrossSection::lastTX=0; // Pointer to the LastArray of TX function |
---|
| 58 | G4double* G4QNuMuNuclearCrossSection::lastQE=0; // Pointer to the LastArray of QE function |
---|
| 59 | G4int G4QNuMuNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile |
---|
| 60 | G4int G4QNuMuNuclearCrossSection::lastN=0; // The last N of calculated nucleus |
---|
| 61 | G4int G4QNuMuNuclearCrossSection::lastZ=0; // The last Z of calculated nucleus |
---|
| 62 | G4double G4QNuMuNuclearCrossSection::lastP=0.; // Last used in cross section Momentum |
---|
| 63 | G4double G4QNuMuNuclearCrossSection::lastTH=0.; // Last threshold momentum |
---|
| 64 | G4double G4QNuMuNuclearCrossSection::lastCS=0.; // Last value of the Cross Section |
---|
| 65 | G4int G4QNuMuNuclearCrossSection::lastI=0; // The last position in the DAMDB |
---|
| 66 | |
---|
| 67 | // Returns Pointer to the G4VQCrossSection class |
---|
| 68 | G4VQCrossSection* 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) |
---|
| 76 | G4double 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) |
---|
| 245 | G4double 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) |
---|
| 264 | G4double 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 | // **************************************************************************************** |
---|
| 374 | G4int 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 | |
---|
| 400 | static const G4double QELX[nE]={0., |
---|
| 401 | .012170,.040879,.057328,.070865,.083129,.094828,.106366,.118013,.129970,.142392, |
---|
| 402 | .155410,.169138,.183676,.199123,.215573,.233120,.251860,.271891,.293317,.316246, |
---|
| 403 | .340796,.367096,.395292,.425547,.458036,.491832,.524989,.556457,.585692,.612377, |
---|
| 404 | .636544,.657790,.676260,.692007,.705323,.716105,.724694,.731347,.736340,.740172, |
---|
| 405 | .742783,.744584,.745804,.746829,.747479,.747995,.748436,.749047,.749497,.749925, |
---|
| 406 | .750486,.750902,.751268,.751566,.752026,.752266,.752428,.752761,.752873,.753094, |
---|
| 407 | .753161,.753164,.753340,.753321}; |
---|
| 408 | // -------------------------------- |
---|
| 409 | G4int first=0; |
---|
| 410 | if(z<0.) |
---|
| 411 | { |
---|
| 412 | first=1; |
---|
| 413 | z=-z; |
---|
| 414 | } |
---|
| 415 | if(z<1 || z>92) // neutron and plutonium are forbidden |
---|
| 416 | { |
---|
| 417 | G4cout<<"***G4QNuMuNuclearCrossSection::GetFunctions:Z="<<z<<".No CS returned"<<G4endl; |
---|
| 418 | return -1; |
---|
| 419 | } |
---|
| 420 | for(G4int k=0; k<nE; k++) |
---|
| 421 | { |
---|
| 422 | G4double a=n+z; |
---|
| 423 | G4double na=n+a; |
---|
| 424 | G4double dn=n+n; |
---|
| 425 | G4double da=a+a; |
---|
| 426 | G4double ta=da+a; |
---|
| 427 | if(first) e[k]=nuEn[k]; // Energy of neutrino E (first bin k=0 can be modified) |
---|
| 428 | t[k]=TOTX[k]*nuEn[k]*(na+na)/ta+QELX[k]*(dn+dn-da)/ta; // TotalCrossSection |
---|
| 429 | q[k]=QELX[k]*dn/a; // QuasiElasticCrossSection |
---|
| 430 | } |
---|
| 431 | return first; |
---|
| 432 | } |
---|
| 433 | |
---|
| 434 | // Randomize Q2 from neutrino to the scattered muon when the scattering is quasi-elastic |
---|
| 435 | G4double G4QNuMuNuclearCrossSection::GetQEL_ExchangeQ2() |
---|
| 436 | { |
---|
| 437 | static const G4double mmu=.105658369;// Mass of muon in GeV |
---|
| 438 | static const G4double mmu2=mmu*mmu; // Squared Mass of muon in GeV^2 |
---|
| 439 | static const double hmmu2=mmu2/2; // .5*m_mu^2 in GeV^2 |
---|
| 440 | static const double MN=.931494043; // Nucleon mass (inside nucleus, atomicMassUnit,GeV) |
---|
| 441 | static const double MN2=MN*MN; // M_N^2 in GeV^2 |
---|
| 442 | static const G4double power=-3.5; // direct power for the magic variable |
---|
| 443 | static const G4double pconv=1./power;// conversion power for the magic variable |
---|
| 444 | static const G4int nQ2=101; // #Of point in the Q2l table (in GeV^2) |
---|
| 445 | static const G4int lQ2=nQ2-1; // index of the last in the Q2l table |
---|
| 446 | static const G4int bQ2=lQ2-1; // index of the before last in the Q2 ltable |
---|
| 447 | // Reversed table |
---|
| 448 | static const G4double Xl[nQ2]={1.87905e-10, |
---|
| 449 | .005231, .010602, .016192, .022038, .028146, .034513, .041130, .047986, .055071, .062374, |
---|
| 450 | .069883, .077587, .085475, .093539, .101766, .110150, .118680, .127348, .136147, .145069, |
---|
| 451 | .154107, .163255, .172506, .181855, .191296, .200825, .210435, .220124, .229886, .239718, |
---|
| 452 | .249617, .259578, .269598, .279675, .289805, .299986, .310215, .320490, .330808, .341169, |
---|
| 453 | .351568, .362006, .372479, .382987, .393527, .404099, .414700, .425330, .435987, .446670, |
---|
| 454 | .457379, .468111, .478866, .489643, .500441, .511260, .522097, .532954, .543828, .554720, |
---|
| 455 | .565628, .576553, .587492, .598447, .609416, .620398, .631394, .642403, .653424, .664457, |
---|
| 456 | .675502, .686557, .697624, .708701, .719788, .730886, .741992, .753108, .764233, .775366, |
---|
| 457 | .786508, .797658, .808816, .819982, .831155, .842336, .853524, .864718, .875920, .887128, |
---|
| 458 | .898342, .909563, .920790, .932023, .943261, .954506, .965755, .977011, .988271, .999539}; |
---|
| 459 | // Direct table |
---|
| 460 | static const G4double Xmax=Xl[lQ2]; |
---|
| 461 | static const G4double Xmin=Xl[0]; |
---|
| 462 | static const G4double dX=(Xmax-Xmin)/lQ2; // step in X(Q2, GeV^2) |
---|
| 463 | static const G4double inl[nQ2]={0, |
---|
| 464 | 1.88843, 3.65455, 5.29282, 6.82878, 8.28390, 9.67403, 11.0109, 12.3034, 13.5583, 14.7811, |
---|
| 465 | 15.9760, 17.1466, 18.2958, 19.4260, 20.5392, 21.6372, 22.7215, 23.7933, 24.8538, 25.9039, |
---|
| 466 | 26.9446, 27.9766, 29.0006, 30.0171, 31.0268, 32.0301, 33.0274, 34.0192, 35.0058, 35.9876, |
---|
| 467 | 36.9649, 37.9379, 38.9069, 39.8721, 40.8337, 41.7920, 42.7471, 43.6992, 44.6484, 45.5950, |
---|
| 468 | 46.5390, 47.4805, 48.4197, 49.3567, 50.2916, 51.2245, 52.1554, 53.0846, 54.0120, 54.9377, |
---|
| 469 | 55.8617, 56.7843, 57.7054, 58.6250, 59.5433, 60.4603, 61.3761, 62.2906, 63.2040, 64.1162, |
---|
| 470 | 65.0274, 65.9375, 66.8467, 67.7548, 68.6621, 69.5684, 70.4738, 71.3784, 72.2822, 73.1852, |
---|
| 471 | 74.0875, 74.9889, 75.8897, 76.7898, 77.6892, 78.5879, 79.4860, 80.3835, 81.2804, 82.1767, |
---|
| 472 | 83.0724, 83.9676, 84.8622, 85.7563, 86.6499, 87.5430, 88.4356, 89.3277, 90.2194, 91.1106, |
---|
| 473 | 92.0013, 92.8917, 93.7816, 94.6711, 95.5602, 96.4489, 97.3372, 98.2252, 99.1128, 100.000}; |
---|
| 474 | G4double Enu=lastE; // Get energy of the last calculated cross-section |
---|
| 475 | G4double dEnu=Enu+Enu; // doubled energy of nu/anu |
---|
| 476 | G4double Enu2=Enu*Enu; // squared energy of nu/anu |
---|
| 477 | G4double ME=Enu*MN; // M*E |
---|
| 478 | G4double dME=ME+ME; // 2*M*E |
---|
| 479 | G4double dEMN=(dEnu+MN)*ME; |
---|
| 480 | G4double MEm=ME-hmmu2; |
---|
| 481 | G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2); |
---|
| 482 | G4double E2M=MN*Enu2-(Enu+MN)*hmmu2; |
---|
| 483 | G4double ymax=(E2M+sqE)/dEMN; |
---|
| 484 | G4double ymin=(E2M-sqE)/dEMN; |
---|
| 485 | G4double rmin=1.-ymin; |
---|
| 486 | G4double rhm2E=hmmu2/Enu2; |
---|
| 487 | G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu) |
---|
| 488 | G4double Q2ma=dME*ymax; // Q2_max(E_nu) |
---|
| 489 | G4double Xma=std::pow((1.+Q2mi),power); // X_max(E_nu) |
---|
| 490 | G4double Xmi=std::pow((1.+Q2ma),power); // X_min(E_nu) |
---|
| 491 | // Find the integral values integ(Xmi) & integ(Xma) using the direct table |
---|
| 492 | G4double rXi=(Xmi-Xmin)/dX; |
---|
| 493 | G4int iXi=static_cast<int>(rXi); |
---|
| 494 | if(iXi<0) iXi=0; |
---|
| 495 | if(iXi>bQ2) iXi=bQ2; |
---|
| 496 | G4double dXi=rXi-iXi; |
---|
| 497 | G4double bnti=inl[iXi]; |
---|
| 498 | G4double inti=bnti+dXi*(inl[iXi+1]-bnti); |
---|
| 499 | // |
---|
| 500 | G4double rXa=(Xma-Xmin)/dX; |
---|
| 501 | G4int iXa=static_cast<int>(rXa); |
---|
| 502 | if(iXa<0) iXa=0; |
---|
| 503 | if(iXa>bQ2) iXa=bQ2; |
---|
| 504 | G4double dXa=rXa-iXa; |
---|
| 505 | G4double bnta=inl[iXa]; |
---|
| 506 | G4double inta=bnta+dXa*(inl[iXa+1]-bnta); |
---|
| 507 | // *** Find X using the reversed table *** |
---|
| 508 | G4double intx=inti+(inta-inti)*G4UniformRand(); |
---|
| 509 | G4int intc=static_cast<int>(intx); |
---|
| 510 | if(intc<0) intc=0; |
---|
| 511 | if(intc>bQ2) intc=bQ2; // If it is more than max, then the BAD extrapolation |
---|
| 512 | G4double dint=intx-intc; |
---|
| 513 | G4double mX=Xl[intc]; |
---|
| 514 | G4double X=mX+dint*(Xl[intc+1]-mX); |
---|
| 515 | G4double Q2=std::pow(X,pconv)-1.; |
---|
| 516 | return Q2*GeV*GeV; |
---|
| 517 | } |
---|
| 518 | |
---|
| 519 | // Randomize Q2 from neutrino to the scattered muon when the scattering is not quasiElastic |
---|
| 520 | G4double G4QNuMuNuclearCrossSection::GetNQE_ExchangeQ2() |
---|
| 521 | { |
---|
| 522 | static const double mpi=.13957018; // charged pi meson mass in GeV |
---|
| 523 | static const G4double mmu=.105658369; // Mass of muon in GeV |
---|
| 524 | static const G4double mmu2=mmu*mmu; // Squared Mass of muon in GeV^2 |
---|
| 525 | static const double hmmu2=mmu2/2; // .5*m_mu^2 in GeV^2 |
---|
| 526 | static const double MN=.931494043; // Nucleon mass (inside nucleus,atomicMassUnit,GeV) |
---|
| 527 | static const double MN2=MN*MN; // M_N^2 in GeV^2 |
---|
| 528 | static const double dMN=MN+MN; // 2*M_N in GeV |
---|
| 529 | static const double mcV=(dMN+mpi)*mpi;// constant of W>M+mc cut for Quasi-Elastic |
---|
| 530 | static const G4int power=7; // direct power for the magic variable |
---|
| 531 | static const G4double pconv=1./power; // conversion power for the magic variable |
---|
| 532 | static const G4int nX=21; // #Of point in the Xl table (in GeV^2) |
---|
| 533 | static const G4int lX=nX-1; // index of the last in the Xl table |
---|
| 534 | static const G4int bX=lX-1; // @@ index of the before last in the Xl table |
---|
| 535 | static const G4int nE=20; // #Of point in the El table (in GeV^2) |
---|
| 536 | static const G4int bE=nE-1; // index of the last in the El table |
---|
| 537 | static const G4int pE=bE-1; // index of the before last in the El table |
---|
| 538 | // Reversed table |
---|
| 539 | static const G4double X0[nX]={6.14081e-05, |
---|
| 540 | .413394, .644455, .843199, 1.02623, 1.20032, 1.36916, 1.53516, 1.70008, 1.86539, 2.03244, |
---|
| 541 | 2.20256, 2.37723, 2.55818, 2.74762, 2.94857, 3.16550, 3.40582, 3.68379, 4.03589, 4.77419}; |
---|
| 542 | static const G4double X1[nX]={.00125268, |
---|
| 543 | .861178, 1.34230, 1.75605, 2.13704, 2.49936, 2.85072, 3.19611, 3.53921, 3.88308, 4.23049, |
---|
| 544 | 4.58423, 4.94735, 5.32342, 5.71700, 6.13428, 6.58447, 7.08267, 7.65782, 8.38299, 9.77330}; |
---|
| 545 | static const G4double X2[nX]={.015694, |
---|
| 546 | 1.97690, 3.07976, 4.02770, 4.90021, 5.72963, 6.53363, 7.32363, 8.10805, 8.89384, 9.68728, |
---|
| 547 | 10.4947, 11.3228, 12.1797, 13.0753, 14.0234, 15.0439, 16.1692, 17.4599, 19.0626, 21.7276}; |
---|
| 548 | static const G4double X3[nX]={.0866877, |
---|
| 549 | 4.03498, 6.27651, 8.20056, 9.96931, 11.6487, 13.2747, 14.8704, 16.4526, 18.0351, 19.6302, |
---|
| 550 | 21.2501, 22.9075, 24.6174, 26.3979, 28.2730, 30.2770, 32.4631, 34.9243, 37.8590, 41.9115}; |
---|
| 551 | static const G4double X4[nX]={.160483, |
---|
| 552 | 5.73111, 8.88884, 11.5893, 14.0636, 16.4054, 18.6651, 20.8749, 23.0578, 25.2318, 27.4127, |
---|
| 553 | 29.6152, 31.8540, 34.1452, 36.5074, 38.9635, 41.5435, 44.2892, 47.2638, 50.5732, 54.4265}; |
---|
| 554 | static const G4double X5[nX]={.0999307, |
---|
| 555 | 5.25720, 8.11389, 10.5375, 12.7425, 14.8152, 16.8015, 18.7296, 20.6194, 22.4855, 24.3398, |
---|
| 556 | 26.1924, 28.0527, 29.9295, 31.8320, 33.7699, 35.7541, 37.7975, 39.9158, 42.1290, 44.4649}; |
---|
| 557 | static const G4double X6[nX]={.0276367, |
---|
| 558 | 3.53378, 5.41553, 6.99413, 8.41629, 9.74057, 10.9978, 12.2066, 13.3796, 14.5257, 15.6519, |
---|
| 559 | 16.7636, 17.8651, 18.9603, 20.0527, 21.1453, 22.2411, 23.3430, 24.4538, 25.5765, 26.7148}; |
---|
| 560 | static const G4double X7[nX]={.00472383, |
---|
| 561 | 2.08253, 3.16946, 4.07178, 4.87742, 5.62140, 6.32202, 6.99034, 7.63368, 8.25720, 8.86473, |
---|
| 562 | 9.45921, 10.0430, 10.6179, 11.1856, 11.7475, 12.3046, 12.8581, 13.4089, 13.9577, 14.5057}; |
---|
| 563 | static const G4double X8[nX]={.000630783, |
---|
| 564 | 1.22723, 1.85845, 2.37862, 2.84022, 3.26412, 3.66122, 4.03811, 4.39910, 4.74725, 5.08480, |
---|
| 565 | 5.41346, 5.73457, 6.04921, 6.35828, 6.66250, 6.96250, 7.25884, 7.55197, 7.84232, 8.13037}; |
---|
| 566 | static const G4double X9[nX]={7.49179e-05, |
---|
| 567 | .772574, 1.16623, 1.48914, 1.77460, 2.03586, 2.27983, 2.51069, 2.73118, 2.94322, 3.14823, |
---|
| 568 | 3.34728, 3.54123, 3.73075, 3.91638, 4.09860, 4.27779, 4.45428, 4.62835, 4.80025, 4.97028}; |
---|
| 569 | static const G4double XA[nX]={8.43437e-06, |
---|
| 570 | .530035, .798454, 1.01797, 1.21156, 1.38836, 1.55313, 1.70876, 1.85712, 1.99956, 2.13704, |
---|
| 571 | 2.27031, 2.39994, 2.52640, 2.65007, 2.77127, 2.89026, 3.00726, 3.12248, 3.23607, 3.34823}; |
---|
| 572 | static const G4double XB[nX]={9.27028e-07, |
---|
| 573 | .395058, .594211, .756726, .899794, 1.03025, 1.15167, 1.26619, 1.37523, 1.47979, 1.58059, |
---|
| 574 | 1.67819, 1.77302, 1.86543, 1.95571, 2.04408, 2.13074, 2.21587, 2.29960, 2.38206, 2.46341}; |
---|
| 575 | static const G4double XC[nX]={1.00807e-07, |
---|
| 576 | .316195, .474948, .604251, .717911, .821417, .917635, 1.00829, 1.09452, 1.17712, 1.25668, |
---|
| 577 | 1.33364, 1.40835, 1.48108, 1.55207, 1.62150, 1.68954, 1.75631, 1.82193, 1.88650, 1.95014}; |
---|
| 578 | static const G4double XD[nX]={1.09102e-08, |
---|
| 579 | .268227, .402318, .511324, .606997, .694011, .774803, .850843, .923097, .992243, 1.05878, |
---|
| 580 | 1.12309, 1.18546, 1.24613, 1.30530, 1.36313, 1.41974, 1.47526, 1.52978, 1.58338, 1.63617}; |
---|
| 581 | static const G4double XE[nX]={1.17831e-09, |
---|
| 582 | .238351, .356890, .453036, .537277, .613780, .684719, .751405, .814699, .875208, .933374, |
---|
| 583 | .989535, 1.04396, 1.09685, 1.14838, 1.19870, 1.24792, 1.29615, 1.34347, 1.38996, 1.43571}; |
---|
| 584 | static const G4double XF[nX]={1.27141e-10, |
---|
| 585 | .219778, .328346, .416158, .492931, .562525, .626955, .687434, .744761, .799494, .852046, |
---|
| 586 | .902729, .951786, .999414, 1.04577, 1.09099, 1.13518, 1.17844, 1.22084, 1.26246, 1.30338}; |
---|
| 587 | static const G4double XG[nX]={1.3713e-11, |
---|
| 588 | .208748, .310948, .393310, .465121, .530069, .590078, .646306, .699515, .750239, .798870, |
---|
| 589 | .845707, .890982, .934882, .977559, 1.01914, 1.05973, 1.09941, 1.13827, 1.17637, 1.21379}; |
---|
| 590 | static const G4double XH[nX]={1.47877e-12, |
---|
| 591 | .203089, .301345, .380162, .448646, .510409, .567335, .620557, .670820, .718647, .764421, |
---|
| 592 | .808434, .850914, .892042, .931967, .970812, 1.00868, 1.04566, 1.08182, 1.11724, 1.15197}; |
---|
| 593 | static const G4double XI[nX]={1.59454e-13, |
---|
| 594 | .201466, .297453, .374007, .440245, .499779, .554489, .605506, .653573, .699213, .742806, |
---|
| 595 | .784643, .824952, .863912, .901672, .938353, .974060, 1.00888, 1.04288, 1.07614, 1.10872}; |
---|
| 596 | static const G4double XJ[nX]={1.71931e-14, |
---|
| 597 | .202988, .297870, .373025, .437731, .495658, .548713, .598041, .644395, .688302, .730147, |
---|
| 598 | .770224, .808762, .845943, .881916, .916805, .950713, .983728, 1.01592, 1.04737, 1.07813}; |
---|
| 599 | // Direct table |
---|
| 600 | static const G4double Xmin[nE]={X0[0],X1[0],X2[0],X3[0],X4[0],X5[0],X6[0],X7[0],X8[0], |
---|
| 601 | X9[0],XA[0],XB[0],XC[0],XD[0],XE[0],XF[0],XG[0],XH[0],XI[0],XJ[0]}; |
---|
| 602 | static const G4double dX[nE]={ |
---|
| 603 | (X0[lX]-X0[0])/lX, (X1[lX]-X1[0])/lX, (X2[lX]-X2[0])/lX, (X3[lX]-X3[0])/lX, |
---|
| 604 | (X4[lX]-X4[0])/lX, (X5[lX]-X5[0])/lX, (X6[lX]-X6[0])/lX, (X7[lX]-X7[0])/lX, |
---|
| 605 | (X8[lX]-X8[0])/lX, (X9[lX]-X9[0])/lX, (XA[lX]-XA[0])/lX, (XB[lX]-XB[0])/lX, |
---|
| 606 | (XC[lX]-XC[0])/lX, (XD[lX]-XD[0])/lX, (XE[lX]-XE[0])/lX, (XF[lX]-XF[0])/lX, |
---|
| 607 | (XG[lX]-XG[0])/lX, (XH[lX]-XH[0])/lX, (XI[lX]-XI[0])/lX, (XJ[lX]-XJ[0])/lX}; |
---|
| 608 | static const G4double* Xl[nE]= |
---|
| 609 | {X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,XA,XB,XC,XD,XE,XF,XG,XH,XI,XJ}; |
---|
| 610 | static const G4double I0[nX]={0, |
---|
| 611 | .411893, 1.25559, 2.34836, 3.60264, 4.96046, 6.37874, 7.82342, 9.26643, 10.6840, 12.0555, |
---|
| 612 | 13.3628, 14.5898, 15.7219, 16.7458, 17.6495, 18.4217, 19.0523, 19.5314, 19.8501, 20.0000}; |
---|
| 613 | static const G4double I1[nX]={0, |
---|
| 614 | .401573, 1.22364, 2.28998, 3.51592, 4.84533, 6.23651, 7.65645, 9.07796, 10.4780, 11.8365, |
---|
| 615 | 13.1360, 14.3608, 15.4967, 16.5309, 17.4516, 18.2481, 18.9102, 19.4286, 19.7946, 20.0000}; |
---|
| 616 | static const G4double I2[nX]={0, |
---|
| 617 | .387599, 1.17339, 2.19424, 3.37090, 4.65066, 5.99429, 7.37071, 8.75427, 10.1232, 11.4586, |
---|
| 618 | 12.7440, 13.9644, 15.1065, 16.1582, 17.1083, 17.9465, 18.6634, 19.2501, 19.6982, 20.0000}; |
---|
| 619 | static const G4double I3[nX]={0, |
---|
| 620 | .366444, 1.09391, 2.04109, 3.13769, 4.33668, 5.60291, 6.90843, 8.23014, 9.54840, 10.8461, |
---|
| 621 | 12.1083, 13.3216, 14.4737, 15.5536, 16.5512, 17.4573, 18.2630, 18.9603, 19.5417, 20.0000}; |
---|
| 622 | static const G4double I4[nX]={0, |
---|
| 623 | .321962, .959681, 1.79769, 2.77753, 3.85979, 5.01487, 6.21916, 7.45307, 8.69991, 9.94515, |
---|
| 624 | 11.1759, 12.3808, 13.5493, 14.6720, 15.7402, 16.7458, 17.6813, 18.5398, 19.3148, 20.0000}; |
---|
| 625 | static const G4double I5[nX]={0, |
---|
| 626 | .257215, .786302, 1.49611, 2.34049, 3.28823, 4.31581, 5.40439, 6.53832, 7.70422, 8.89040, |
---|
| 627 | 10.0865, 11.2833, 12.4723, 13.6459, 14.7969, 15.9189, 17.0058, 18.0517, 19.0515, 20.0000}; |
---|
| 628 | static const G4double I6[nX]={0, |
---|
| 629 | .201608, .638914, 1.24035, 1.97000, 2.80354, 3.72260, 4.71247, 5.76086, 6.85724, 7.99243, |
---|
| 630 | 9.15826, 10.3474, 11.5532, 12.7695, 13.9907, 15.2117, 16.4275, 17.6337, 18.8258, 20.0000}; |
---|
| 631 | static const G4double I7[nX]={0, |
---|
| 632 | .168110, .547208, 1.07889, 1.73403, 2.49292, 3.34065, 4.26525, 5.25674, 6.30654, 7.40717, |
---|
| 633 | 8.55196, 9.73492, 10.9506, 12.1940, 13.4606, 14.7460, 16.0462, 17.3576, 18.6767, 20.0000}; |
---|
| 634 | static const G4double I8[nX]={0, |
---|
| 635 | .150652, .497557, .990048, 1.60296, 2.31924, 3.12602, 4.01295, 4.97139, 5.99395, 7.07415, |
---|
| 636 | 8.20621, 9.38495, 10.6057, 11.8641, 13.1561, 14.4781, 15.8267, 17.1985, 18.5906, 20.0000}; |
---|
| 637 | static const G4double I9[nX]={0, |
---|
| 638 | .141449, .470633, .941304, 1.53053, 2.22280, 3.00639, 3.87189, 4.81146, 5.81837, 6.88672, |
---|
| 639 | 8.01128, 9.18734, 10.4106, 11.6772, 12.9835, 14.3261, 15.7019, 17.1080, 18.5415, 20.0000}; |
---|
| 640 | static const G4double IA[nX]={0, |
---|
| 641 | .136048, .454593, .912075, 1.48693, 2.16457, 2.93400, 3.78639, 4.71437, 5.71163, 6.77265, |
---|
| 642 | 7.89252, 9.06683, 10.2916, 11.5631, 12.8780, 14.2331, .625500, 17.0525, 18.5115, 20.0000}; |
---|
| 643 | static const G4double IB[nX]={0, |
---|
| 644 | .132316, .443455, .891741, 1.45656, 2.12399, 2.88352, 3.72674, 4.64660, 5.63711, 6.69298, |
---|
| 645 | 7.80955, 8.98262, 10.2084, 11.4833, 12.8042, 14.1681, 15.5721, 17.0137, 18.4905, 20.0000}; |
---|
| 646 | static const G4double IC[nX]={0, |
---|
| 647 | .129197, .434161, .874795, 1.43128, 2.09024, 2.84158, 3.67721, 4.59038, 5.57531, 6.62696, |
---|
| 648 | 7.74084, 8.91291, 10.1395, 11.4173, 12.7432, 14.1143, 15.5280, 16.9817, 18.4731, 20.0000}; |
---|
| 649 | static const G4double ID[nX]={0, |
---|
| 650 | .126079, .424911, .857980, 1.40626, 2.05689, 2.80020, 3.62840, 4.53504, 5.51456, 6.56212, |
---|
| 651 | 7.67342, 8.84458, 10.0721, 11.3527, 12.6836, 14.0618, 15.4849, 16.9504, 18.4562, 20.0000}; |
---|
| 652 | static const G4double IE[nX]={0, |
---|
| 653 | .122530, .414424, .838964, 1.37801, 2.01931, 2.75363, 3.57356, 4.47293, 5.44644, 6.48949, |
---|
| 654 | 7.59795, 8.76815, 9.99673, 11.2806, 12.6170, 14.0032, 15.4369, 16.9156, 18.4374, 20.0000}; |
---|
| 655 | static const G4double IF[nX]={0, |
---|
| 656 | .118199, .401651, .815838, 1.34370, 1.97370, 2.69716, 3.50710, 4.39771, 5.36401, 6.40164, |
---|
| 657 | 7.50673, 8.67581, 9.90572, 11.1936, 12.5367, 13.9326, 15.3790, 16.8737, 18.4146, 20.0000}; |
---|
| 658 | static const G4double IG[nX]={0, |
---|
| 659 | .112809, .385761, .787075, 1.30103, 1.91700, 2.62697, 3.42451, 4.30424, 5.26158, 6.29249, |
---|
| 660 | 7.39341, 8.56112, 9.79269, 11.0855, 12.4369, 13.8449, 15.3071, 16.8216, 18.3865, 20.0000}; |
---|
| 661 | static const G4double IH[nX]={0, |
---|
| 662 | .106206, .366267, .751753, 1.24859, 1.84728, 2.54062, 3.32285, .189160, 5.13543, 6.15804, |
---|
| 663 | 7.25377, 8.41975, 9.65334, 10.9521, 12.3139, 13.7367, 15.2184, 16.7573, 18.3517, 20.0000}; |
---|
| 664 | static const G4double II[nX]={0, |
---|
| 665 | .098419, .343194, .709850, 1.18628, 1.76430, 2.43772, 3.20159, 4.05176, 4.98467, 5.99722, |
---|
| 666 | 7.08663, 8.25043, 9.48633, 10.7923, 12.1663, 13.6067, 15.1118, 16.6800, 18.3099, 20.0000}; |
---|
| 667 | static const G4double IJ[nX]={0, |
---|
| 668 | .089681, .317135, .662319, 1.11536, 1.66960, 2.32002, 3.06260, 3.89397, 4.81126, 5.81196, |
---|
| 669 | 6.89382, 8.05483, 9.29317, 10.6072, 11.9952, 13.4560, 14.9881, 16.5902, 18.2612, 20.0000}; |
---|
| 670 | static const G4double* Il[nE]= |
---|
| 671 | {I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,IA,IB,IC,ID,IE,IF,IG,IH,II,IJ}; |
---|
| 672 | static const G4double lE[nE]={ |
---|
| 673 | -1.98842,-1.58049,-1.17256,-.764638,-.356711, .051215, .459141, .867068, 1.27499, 1.68292, |
---|
| 674 | 2.09085, 2.49877, 2.90670, 3.31463, 3.72255, 4.13048, 4.53840, 4.94633, 5.35426, 5.76218}; |
---|
| 675 | static const G4double lEmi=lE[0]; |
---|
| 676 | static const G4double lEma=lE[nE-1]; |
---|
| 677 | static const G4double dlE=(lEma-lEmi)/bE; |
---|
| 678 | //*************************************************************************************** |
---|
| 679 | G4double Enu=lastE; // Get energy of the last calculated cross-section |
---|
| 680 | G4double lEn=std::log(Enu); // log(E) for interpolation |
---|
| 681 | G4double rE=(lEn-lEmi)/dlE; // Position of the energy |
---|
| 682 | G4int fE=static_cast<int>(rE); // Left bin for interpolation |
---|
| 683 | if(fE<0) fE=0; |
---|
| 684 | if(fE>pE)fE=pE; |
---|
| 685 | G4int sE=fE+1; // Right bin for interpolation |
---|
| 686 | G4double dE=rE-fE; // relative log shift from the left bin |
---|
| 687 | G4double dEnu=Enu+Enu; // doubled energy of nu/anu |
---|
| 688 | G4double Enu2=Enu*Enu; // squared energy of nu/anu |
---|
| 689 | G4double Emu=Enu-mmu; // Free Energy of neutrino/anti-neutrino |
---|
| 690 | G4double ME=Enu*MN; // M*E |
---|
| 691 | G4double dME=ME+ME; // 2*M*E |
---|
| 692 | G4double dEMN=(dEnu+MN)*ME; |
---|
| 693 | G4double MEm=ME-hmmu2; |
---|
| 694 | G4double sqE=Enu*std::sqrt(MEm*MEm-mmu2*MN2); |
---|
| 695 | G4double E2M=MN*Enu2-(Enu+MN)*hmmu2; |
---|
| 696 | G4double ymax=(E2M+sqE)/dEMN; |
---|
| 697 | G4double ymin=(E2M-sqE)/dEMN; |
---|
| 698 | G4double rmin=1.-ymin; |
---|
| 699 | G4double rhm2E=hmmu2/Enu2; |
---|
| 700 | G4double Q2mi=(Enu2+Enu2)*(rmin-rhm2E-std::sqrt(rmin*rmin-rhm2E-rhm2E)); // Q2_min(E_nu) |
---|
| 701 | G4double Q2ma=dME*ymax; // Q2_max(E_nu) |
---|
| 702 | G4double Q2nq=Emu*dMN-mcV; |
---|
| 703 | if(Q2ma>Q2nq) Q2ma=Q2nq; // Correction for Non Quasi Elastic |
---|
| 704 | // --- now r_min=Q2mi/Q2ma and r_max=1.; when r is randomized -> Q2=r*Q2ma --- |
---|
| 705 | G4double Rmi=Q2mi/Q2ma; |
---|
| 706 | G4double shift=1.+.9673/(1.+.323/Enu/Enu)/std::pow(Enu,.78); //@@ different for anti-nu |
---|
| 707 | // --- E-interpolation must be done in a log scale --- |
---|
| 708 | G4double Xmi=std::pow((shift-Rmi),power);// X_min(E_nu) |
---|
| 709 | G4double Xma=std::pow((shift-1.),power); // X_max(E_nu) |
---|
| 710 | // Find the integral values integ(Xmi) & integ(Xma) using the direct table |
---|
| 711 | G4double idX=dX[fE]+dE*(dX[sE]-dX[fE]); // interpolated X step |
---|
| 712 | G4double iXmi=Xmin[fE]+dE*(Xmin[sE]-Xmin[fE]); // interpolated X minimum |
---|
| 713 | G4double rXi=(Xmi-iXmi)/idX; |
---|
| 714 | G4int iXi=static_cast<int>(rXi); |
---|
| 715 | if(iXi<0) iXi=0; |
---|
| 716 | if(iXi>bX) iXi=bX; |
---|
| 717 | G4double dXi=rXi-iXi; |
---|
| 718 | G4double bntil=Il[fE][iXi]; |
---|
| 719 | G4double intil=bntil+dXi*(Il[fE][iXi+1]-bntil); |
---|
| 720 | G4double bntir=Il[sE][iXi]; |
---|
| 721 | G4double intir=bntir+dXi*(Il[sE][iXi+1]-bntir); |
---|
| 722 | G4double inti=intil+dE*(intir-intil);// interpolated begin of the integral |
---|
| 723 | // |
---|
| 724 | G4double rXa=(Xma-iXmi)/idX; |
---|
| 725 | G4int iXa=static_cast<int>(rXa); |
---|
| 726 | if(iXa<0) iXa=0; |
---|
| 727 | if(iXa>bX) iXa=bX; |
---|
| 728 | G4double dXa=rXa-iXa; |
---|
| 729 | G4double bntal=Il[fE][iXa]; |
---|
| 730 | G4double intal=bntal+dXa*(Il[fE][iXa+1]-bntal); |
---|
| 731 | G4double bntar=Il[sE][iXa]; |
---|
| 732 | G4double intar=bntar+dXa*(Il[sE][iXa+1]-bntar); |
---|
| 733 | G4double inta=intal+dE*(intar-intal);// interpolated end of the integral |
---|
| 734 | // |
---|
| 735 | // *** Find X using the reversed table *** |
---|
| 736 | G4double intx=inti+(inta-inti)*G4UniformRand(); |
---|
| 737 | G4int intc=static_cast<int>(intx); |
---|
| 738 | if(intc<0) intc=0; |
---|
| 739 | if(intc>bX) intc=bX; |
---|
| 740 | G4double dint=intx-intc; |
---|
| 741 | G4double mXl=Xl[fE][intc]; |
---|
| 742 | G4double Xlb=mXl+dint*(Xl[fE][intc+1]-mXl); |
---|
| 743 | G4double mXr=Xl[sE][intc]; |
---|
| 744 | G4double Xrb=mXr+dint*(Xl[sE][intc+1]-mXr); |
---|
| 745 | G4double X=Xlb+dE*(Xrb-Xlb); // interpolated X value |
---|
| 746 | G4double R=shift-std::pow(X,pconv); |
---|
| 747 | G4double Q2=R*Q2ma; |
---|
| 748 | return Q2*GeV*GeV; |
---|
| 749 | } |
---|
| 750 | |
---|
| 751 | // It returns a fraction of the direct interaction of the neutrino with quark-partons |
---|
| 752 | G4double G4QNuMuNuclearCrossSection::GetDirectPart(G4double Q2) |
---|
| 753 | { |
---|
| 754 | G4double f=Q2/4.62; |
---|
| 755 | G4double ff=f*f; |
---|
| 756 | G4double r=ff*ff; |
---|
| 757 | G4double s=std::pow((1.+.6/Q2),(-1.-(1.+r)/(12.5+r/.3))); |
---|
| 758 | //@@ It is the same for nu/anu, but for nu it is a bit less, and for anu a bit more (par) |
---|
| 759 | return 1.-s*(1.-s/2); |
---|
| 760 | } |
---|
| 761 | |
---|
| 762 | // #of quark-partons in the nonperturbative phase space is the same for neut and anti-neut |
---|
| 763 | G4double G4QNuMuNuclearCrossSection::GetNPartons(G4double Q2) |
---|
| 764 | { |
---|
| 765 | return 3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space |
---|
| 766 | } |
---|
| 767 | |
---|
| 768 | // This class can provide only virtual exchange pi+ (a substitute for W+ boson) |
---|
| 769 | G4int G4QNuMuNuclearCrossSection::GetExchangePDGCode() {return 211;} |
---|