[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 | // G4 Tools program: Glauber Hadron-Nuclear Cross Sections |
---|
| 28 | // ..................................................... |
---|
| 29 | // Created: M.V. Kossov, CERN/ITEP(Moscow), 20-Jan-2005 |
---|
| 30 | // |
---|
| 31 | //===================================================================== |
---|
| 32 | |
---|
| 33 | //#define debug |
---|
| 34 | //#define bestest |
---|
| 35 | //#define integrcs |
---|
| 36 | #define eldiffcs |
---|
| 37 | |
---|
| 38 | #include "globals.hh" |
---|
| 39 | #include <iostream> |
---|
| 40 | #include <fstream> |
---|
| 41 | #include <vector> |
---|
| 42 | #include "G4ios.hh" |
---|
| 43 | #include "G4QBesIKJY.hh" |
---|
| 44 | #include "G4QPDGCode.hh" |
---|
| 45 | #include "G4NucleiPropertiesTable.hh" |
---|
| 46 | #include "G4NucleiProperties.hh" |
---|
| 47 | //#include <CLHEP/GenericFunctions/Bessel.hh> |
---|
| 48 | |
---|
| 49 | // Integrated cross-sections |
---|
| 50 | G4double TotalCrSec, InelCrSec, ProdCrSec, ElasticCrSec, QuasyElasticCrSec; |
---|
| 51 | // Parameters of hadron-nucleon interaction |
---|
| 52 | G4double HadrTot, HadrSlope, HadrReIm, DDSect2, DDSect3, Tot00; |
---|
| 53 | // CHIPS parameters of hadron-proton & hadron-neutron interaction (Re/Im is common) |
---|
| 54 | G4double HadrTElN, HadrTElP, HadrSlpN, HadrSlpP, HadrPexN, HadrPexP; |
---|
| 55 | |
---|
| 56 | // CHIPS Parameters (only for protons, nn is necessary for neutrons np+nn) |
---|
| 57 | void CHIPSParameters(G4int iPDG, G4double HadrMoment) // HadrMoment is in MeV |
---|
| 58 | { |
---|
| 59 | G4double p=HadrMoment/GeV; // Momentum in GeV |
---|
| 60 | G4double p2=p*p; // Squared momentum in GeV^2 |
---|
| 61 | G4double sp=std::sqrt(p); |
---|
| 62 | G4double p2s=p2*sp; |
---|
| 63 | G4double p3=p2*p; |
---|
| 64 | G4double p4=p2*p2; |
---|
| 65 | G4double p5=p4*p; |
---|
| 66 | G4double lp=std::log(p); |
---|
| 67 | G4double d2=(lp-5.)*(lp-5.); |
---|
| 68 | G4double p8=p4*p4; |
---|
| 69 | if (iPDG==2212) |
---|
| 70 | { |
---|
| 71 | HadrTElN=12./(p2s+.05*p+.0001/sp)+.35/p+(6.75+.14*d2+19./p)/(1.+.6/p4); // pn (mb) |
---|
| 72 | HadrTElP=2.91/(p2s+.0024)+(6.75+.14*d2+23./p)/(1.+.4/p4); // pp (mb) |
---|
| 73 | HadrSlpN=(7.2+4.32/(p8+.012*p3))/(1.+2.5/p4); // pn (GeV^-2) |
---|
| 74 | HadrSlpP=8.*std::pow(p,.055)/(1.+3.64/p4); // pp (GeV^-2) |
---|
| 75 | HadrPexN=(6.75+.14*d2+13./p)/(1.+.14/p4)+.6/(p4+.00013); // pn (mb/sr) |
---|
| 76 | HadrPexP=(74.+3.*d2)/(1.+3.4/p5)+(.2/p2+17.*p)/(p4+.001*sp); // pn (mb/sr) |
---|
| 77 | HadrReIm=-.54+.0954*lp+.975/(p2+.09865/p); // no unit (the same for pp & pn) |
---|
| 78 | //G4cout<<"CHIPSPar:p="<<p<<",N="<<HadrTotN<<",P="<<HadrTotN<<",RI="<<HadrReIm<<G4endl; |
---|
| 79 | } |
---|
| 80 | else G4cout<<"CHIPSPar: Only proton projectile (2212) is implemented PDG="<<iPDG<<G4endl; |
---|
| 81 | } |
---|
| 82 | |
---|
| 83 | // Dubna Parameters |
---|
| 84 | G4int CalculateParameters(G4int iPDG, G4double HadrMoment) // HadrMoment is in MeV |
---|
| 85 | { |
---|
| 86 | static const G4double mN=.938; // Dubna's mass of the nucleon in nuclei in GeV |
---|
| 87 | //static const G4double mN=.9315; // AtomicMassUnit (mass of the nucleon in nuclei) in GeV |
---|
| 88 | static const G4double dmN=mN+mN; // Doubled mass of the nucleon in nuclei in GeV |
---|
| 89 | static const G4double mN2=mN*mN; // Squared mass of the nucleon in nuclei in GeV^2 |
---|
| 90 | G4int iHadron=-1; |
---|
| 91 | if(iPDG==2212||iPDG==2112||iPDG==3122||iPDG==3222||iPDG==3112|| |
---|
| 92 | iPDG==3212||iPDG==3312||iPDG==3322||iPDG==3334) iHadron = 0; |
---|
| 93 | else if(iPDG==-2212||iPDG==-2112||iPDG==-3122||iPDG==-3222|| |
---|
| 94 | iPDG==-3112||iPDG==-3212||iPDG==-3312||iPDG==-3322|| |
---|
| 95 | iPDG==-3334) iHadron = 1; // Anti-nucleons, Anti-hyperons |
---|
| 96 | else if(iPDG== 211) iHadron = 2; // Pi+ |
---|
| 97 | else if(iPDG==-211) iHadron = 3; // Pi- |
---|
| 98 | else if(iPDG== 321) iHadron = 4; // K+ |
---|
| 99 | else if(iPDG==-321) iHadron = 5; // K- @@ What about K0? |
---|
| 100 | else |
---|
| 101 | { |
---|
| 102 | G4cout<<"CalculateParameters: PDG= "<<iPDG<<" is not supported"<<G4endl; |
---|
| 103 | return iHadron; |
---|
| 104 | } |
---|
| 105 | G4double mHadr = G4QPDGCode(iPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
| 106 | G4double mHadr2 = mHadr*mHadr; // Squared Hadron Mass in GeV |
---|
| 107 | HadrMoment/= 1000.; // Momentum in GeV (input parameter) |
---|
| 108 | G4double HadrMoment2= HadrMoment*HadrMoment; // Squared Momentum in GeV^2 |
---|
| 109 | G4double HadrEnergy2= mHadr2+HadrMoment2; // Tot energy in GeV |
---|
| 110 | G4double HadrEnergy = std::sqrt(HadrEnergy2); // Tot energy in GeV (global parameter) |
---|
| 111 | G4double sHadr = HadrEnergy*dmN+mN2+mHadr2; // s in GeV^2 |
---|
| 112 | G4double sqrS = std::sqrt(sHadr); // W in GeV |
---|
| 113 | //G4cout<<"GHAD: E="<<HadrEnergy<<",dM="<<dmN<<",sN="<<mN2<<",sH="<<mHadr2<<G4endl; |
---|
| 114 | switch (iHadron) |
---|
| 115 | { |
---|
| 116 | case 0: // =====>>> proton or hyperons (the same?) |
---|
| 117 | //G4double Delta=1; // LowEnergy corr |
---|
| 118 | //if(HadrEnergy<40.) Delta = 0.916+0.0021*HadrEnergy; |
---|
| 119 | HadrTot = 5.2+5.2*std::log(HadrEnergy)+51*std::pow(HadrEnergy,-0.35); // mb |
---|
| 120 | HadrSlope = 6.44+0.88*std::log(sHadr)-1; // GeV^-2 |
---|
| 121 | HadrReIm = 0.13*std::log(sHadr/350)*std::pow(sHadr,-0.18); // no unit |
---|
| 122 | //G4cout<<"GHAD:S="<<sHadr<<",T="<<HadrTot<<",R="<<HadrReIm<<",B="<<HadrSlope<<G4endl; |
---|
| 123 | DDSect2 = 11.; // mb*GeV^-2 |
---|
| 124 | DDSect3 = 3.; // mb*GeV^-2 |
---|
| 125 | // Hyperon corrections |
---|
| 126 | if(iPDG==3122||iPDG==3222||iPDG==3112||iPDG==3212) //Lambda,Sig0+- |
---|
| 127 | { |
---|
| 128 | HadrTot *= 0.80; |
---|
| 129 | HadrSlope *= 0.85; |
---|
| 130 | } |
---|
| 131 | else if(iPDG==3312||iPDG==3322) // ---> Xi0- |
---|
| 132 | { |
---|
| 133 | HadrTot *= 0.70; |
---|
| 134 | HadrSlope *= 0.75; |
---|
| 135 | } |
---|
| 136 | else if(iPDG==3334) // ---> Omega- |
---|
| 137 | { |
---|
| 138 | HadrTot *= 0.60; |
---|
| 139 | HadrSlope *= 0.65; |
---|
| 140 | } |
---|
| 141 | break; |
---|
| 142 | case 1: // =====>>> antiproton and antihyperons |
---|
| 143 | HadrTot = 5.2+5.2*std::log(HadrEnergy)+123.2*std::pow(HadrEnergy,-0.5); // mb |
---|
| 144 | HadrSlope = 8.32+0.57*std::log(sHadr); // GeV^-2 |
---|
| 145 | if(HadrEnergy<1000.) HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8); |
---|
| 146 | else HadrReIm = 0.6*std::log(sHadr/350)*std::pow(sHadr,-0.25); |
---|
| 147 | DDSect2 = 11.; //mb*GeV-2 @@ the same as in case 0 |
---|
| 148 | DDSect3 = 3.; //mb*GeV-2 @@ the same as in case 0 |
---|
| 149 | // Hyperon corrections |
---|
| 150 | if(iPDG==-3122||iPDG==-3222||iPDG==-3112||iPDG==-3212)//AntiLam,Sig |
---|
| 151 | { |
---|
| 152 | HadrTot *=0.75; |
---|
| 153 | HadrSlope *=0.85; |
---|
| 154 | } |
---|
| 155 | else if(iPDG==-3312||iPDG==-3322) // =====>>> AntiXi0- |
---|
| 156 | { |
---|
| 157 | HadrTot *=0.65; |
---|
| 158 | HadrSlope *=0.75; |
---|
| 159 | } |
---|
| 160 | if(iPDG==-3334) // ======>>>AntiOmega- |
---|
| 161 | { |
---|
| 162 | HadrTot *=0.55; |
---|
| 163 | HadrSlope *=0.65; |
---|
| 164 | } |
---|
| 165 | break; |
---|
| 166 | case 2: // =====>>> pi plus |
---|
| 167 | if(HadrMoment>2.) HadrTot = 10.6+2.*std::log(HadrEnergy)+25*std::pow(HadrEnergy,-0.43);// mb |
---|
| 168 | else HadrTot = 40-50*(HadrMoment-1.5)*(HadrMoment-1.7); |
---|
| 169 | HadrSlope = 7.28+0.245*std::log(sHadr); // GeV^-2 |
---|
| 170 | HadrReIm = 0.2*std::log(sHadr/100)*std::pow(sHadr,-0.15); // no dim |
---|
| 171 | DDSect2 = 4.6; // mb*GeV^-2 |
---|
| 172 | DDSect3 = 1.33; // mb*GeV^-2 |
---|
| 173 | break; |
---|
| 174 | case 3: // =====>>> pi minus |
---|
| 175 | HadrTot = 10.6+2*std::log(HadrEnergy)+30*std::pow(HadrEnergy,-0.43); // mb @@ E->inf? |
---|
| 176 | if(HadrMoment<1.399) HadrTot = HadrTot+21.0/0.4*(1.4-HadrMoment); // A huge jump ? |
---|
| 177 | HadrSlope = 7.28+0.245*std::log(sHadr); // GeV-2 |
---|
| 178 | HadrReIm = 0.2*std::log(sHadr/100)*std::pow(sHadr,-0.15); |
---|
| 179 | DDSect2 = 4.6; //mb*GeV-2 @@ the same as in case 2 |
---|
| 180 | DDSect3 = 1.33; //mb*GeV-2 |
---|
| 181 | break; |
---|
| 182 | case 4: // =====>>> K plus |
---|
| 183 | HadrTot = 10.6+1.8*std::log(HadrEnergy)+9.*std::pow(HadrEnergy,-0.55); // mb |
---|
| 184 | if(HadrEnergy>100) HadrSlope = 15.; // Jump at 100 GeV (?) |
---|
| 185 | else HadrSlope = 1.0+1.76*std::log(sHadr)-2.84*std::pow(sHadr,-0.5); // GeV-2 |
---|
| 186 | //else HadrSlope = 5.28+1.76*std::log(sHadr)-2.84*std::pow(sHadr,-0.5); // GeV-2 |
---|
| 187 | HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1); |
---|
| 188 | DDSect2 = 3.5; //mb*GeV-2 |
---|
| 189 | DDSect3 = 1.03; //mb*GeV-2 |
---|
| 190 | break; |
---|
| 191 | case 5: // ======>>> K minus |
---|
| 192 | HadrTot = 10+1.8*std::log(HadrEnergy)+25*std::pow(HadrEnergy,-0.5); // mb |
---|
| 193 | HadrSlope = 6.98+0.127*std::log(sHadr); // GeV-2 |
---|
| 194 | //if(HadrEnergy<8) HadrReIm = 0.7; |
---|
| 195 | HadrReIm = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1); |
---|
| 196 | DDSect2 = 3.5; //mb*GeV-2 |
---|
| 197 | DDSect3 = 1.03; //mb*GeV-2 |
---|
| 198 | break; |
---|
| 199 | } |
---|
| 200 | return iHadron; |
---|
| 201 | } |
---|
| 202 | |
---|
| 203 | void CalculateIntegralCS(G4int Anuc, G4double HadrEnergy) |
---|
| 204 | { |
---|
| 205 | static const G4double eps = 0.0001; // Accuracy of calculations |
---|
| 206 | static const G4double mb2G2 = 2.568; // Transformation from mb to GeV^-2 |
---|
| 207 | static const G4double incrf = 10; // Increase factor of radius for intergation |
---|
| 208 | static const G4double m2G10 = incrf*mb2G2; // Big conversion factor |
---|
| 209 | static const G4double dR0 = 1.06+1.06; // Doubled nominal R0 |
---|
| 210 | G4double An = Anuc; // A float |
---|
| 211 | G4double Stot = HadrTot*mb2G2; // Total hadron-nucleon CS in GeV-2 |
---|
| 212 | G4double Bhad = HadrSlope; // Diffraction cone nH In GeV-2 |
---|
| 213 | G4double Asq = 1+HadrReIm*HadrReIm; // |M|^2/(ImM)^2 |
---|
| 214 | |
---|
| 215 | G4double R0 = std::sqrt(0.99); // in fermi (strange value? 1.07?) |
---|
| 216 | |
---|
| 217 | if (Anuc == 58) R0 = std::sqrt(0.6); // Personal correction for Fe |
---|
| 218 | if (Anuc >20) R0 = std::sqrt((35.34+0.5*Anuc)/(40.97+Anuc)); |
---|
| 219 | if (Anuc == 16) R0 = std::sqrt(0.75); // Personal correction for O16 |
---|
| 220 | if (Anuc >10) R0 = std::sqrt(0.84); // The same R0 for A=11-20 (except O16) |
---|
| 221 | if (Anuc == 4) R0 *= 1.2; // Personal correction for Helium |
---|
| 222 | |
---|
| 223 | G4double Rnucl = R0*std::pow(static_cast<double>(Anuc),.3333333); // in fermi |
---|
| 224 | G4double Rnuc2 = Rnucl*Rnucl*m2G10; // in GeV^-2 (10 -> to make it >> R) |
---|
| 225 | G4double RB = Rnuc2+Bhad; // Effective increase of R2 (Convolution?) |
---|
| 226 | G4double R2B = RB+Bhad; // Over-convolution (?) |
---|
| 227 | G4double Delta = Stot/R2B/twopi; // Integration step for Total & Inelastic CS |
---|
| 228 | G4double Delt2 = Delta+Delta; // Two delta |
---|
| 229 | G4double Delt3 = Stot*Asq*R2B/RB/Bhad/16/pi; // Integration step for Production CS |
---|
| 230 | |
---|
| 231 | G4double Tot0 = 0.; // =SUM[(-Delta)^(i-1)*(C^A_i)/i] |
---|
| 232 | G4double Inel0 = 0.; // =SUM[(-Delta)^(i-1)*(2A)!/(2A-i)!/i!/i] |
---|
| 233 | G4double N = -1./Delta; |
---|
| 234 | G4double N1 = N; |
---|
| 235 | G4double N3 = -1./Delt2; |
---|
| 236 | G4double Prod0 = 0.; |
---|
| 237 | G4int Ai = Anuc+1; // Working value to avoid multiplication |
---|
| 238 | G4double RBi = 0.; // Working value to avoid multiplication |
---|
| 239 | for (G4int i=1; i<= Anuc; i++) // ==> Loop over nucleons |
---|
| 240 | { |
---|
| 241 | Ai--; // Working value to avoid multiplication |
---|
| 242 | RBi += RB; // Working value to avoid multiplication |
---|
| 243 | G4double Di=Delta/i; // Working value to avoid repetition |
---|
| 244 | N *= -Di*Ai; // (-Delta)^(i-1)*A!/(A-i)!/i! (C^A_i) |
---|
| 245 | N1 *= -Di*(Anuc+Ai); // (-Delta)^(i-1)*(2A)!/(2A-i)!/i! |
---|
| 246 | N3 *= -Delt2*Ai/i; // (-2*Delt)^(i-1)*A!/(A-i)!/i! (C^A_i) |
---|
| 247 | Tot0 += N/i; // =SUM[(-Delta)^(i-1)*(C^A_i)/i] |
---|
| 248 | Inel0 += N1/i; // =SUM[(-Delta)^(i-1)*(2A)!/(2A-i)!/i!/i] |
---|
| 249 | //G4double N2 = 1.; |
---|
| 250 | G4double N4 = 1.; |
---|
| 251 | //G4double Inel1 = 0.; |
---|
| 252 | //G4double Inel2 = 0.; |
---|
| 253 | G4double Prod1 = 0.; |
---|
| 254 | G4double Bhl = 0.; |
---|
| 255 | for (G4int l=0; l<= i; l++) // ==> Loop over nucleon-nucleon pairs (0?) |
---|
| 256 | { |
---|
| 257 | //Inel1 += N2/(i+l); // SUM[(-Delta)^(l-1)*i!(i-1)!/(l-1)!/l!/(i+l)!] |
---|
| 258 | Prod1 += N4*RB/(RBi+Bhl); // |
---|
| 259 | //N2 *= -Delt*(i-l)/(l+1); // (-Delta)^l*i!/l!/(l+1)! |
---|
| 260 | N4 *= -Delt3*(i-l)/(l+1); |
---|
| 261 | Bhl += Bhad; |
---|
| 262 | } |
---|
| 263 | //Inel2 += Inel1*N3; |
---|
| 264 | Prod0 += Prod1*N3; |
---|
| 265 | if(std::fabs(N1/i/Inel0) < eps) break; |
---|
| 266 | } |
---|
| 267 | |
---|
| 268 | Tot0 *= HadrTot; // Multiply by the hN cross section |
---|
| 269 | Inel0 *= HadrTot/2; // Multiply by a half of hN cross section |
---|
| 270 | //Inel2 *= HadrTot; // Multiply by the hN cross section |
---|
| 271 | Prod0 *= HadrTot; // Multiply by the hN cross section |
---|
| 272 | Tot00 = Tot0; // Remember because it will be changed |
---|
| 273 | G4double ak = Rnuc2*twopi/Stot; // cross-section ratio for integration |
---|
| 274 | G4double Ank = An/ak; // Working value to avoid repetition |
---|
| 275 | G4double rak = 1./ak; // Working value to avoid repetition |
---|
| 276 | G4double ak1 = 1.-rak/4; // Working value to avoid repetition |
---|
| 277 | G4double Ank1 = Ank*ak1; // Working value to avoid repetition |
---|
| 278 | G4double DDSect1 = (DDSect2+DDSect3*std::log(dR0*HadrEnergy/Rnucl/std::sqrt(m2G10)/4)); |
---|
| 279 | G4double Dtot = 8*pi*ak/HadrTot*(1.-(1.+Ank)*std::exp(-Ank))*DDSect1/mb2G2;//TotCSCor |
---|
| 280 | G4double bk = (1.-rak)/Stot/ak1; // Working value to avoid power(,2) |
---|
| 281 | G4double bd = bk*bk*DDSect1*(1.-(1.+Ank1)*std::exp(-Ank1))*Rnuc2; // Working |
---|
| 282 | G4double Dprod = bd*4*pi2*mb2G2; // Correction for the Production & Inel CS |
---|
| 283 | TotalCrSec = Tot0-Dtot; // Total Cross Section (primary) |
---|
| 284 | InelCrSec = Inel0-Dprod; // Inelastic Cross Section (primary) |
---|
| 285 | //InelCrSec1 = Inel2; // Direct calculation without correction |
---|
| 286 | ProdCrSec = Prod0-Dprod; // Production Cross Section (primary) |
---|
| 287 | ElasticCrSec = TotalCrSec-InelCrSec; // Elastic Cross Section (secondary) |
---|
| 288 | QuasyElasticCrSec = InelCrSec-ProdCrSec; // Quasy-Elastic Cross Section (secondary) |
---|
| 289 | } |
---|
| 290 | |
---|
| 291 | #ifdef eldiffcs |
---|
| 292 | // Functions for calculation of differential elastic cross-sections |
---|
| 293 | // Static globals |
---|
| 294 | static const G4double rAmax = 2.5; // Factor of maximum radius |
---|
| 295 | static const G4int NptB = 500; // A#of steps in the impact parameter integration |
---|
| 296 | static const G4int nFact = 253; // Maximum A (for Factorials and Factors) |
---|
| 297 | |
---|
| 298 | // Globals |
---|
| 299 | G4double r0 = 1.1; // Wood-Saxon's factor |
---|
| 300 | G4double Mnoj[nFact], Factorials[nFact]; // Factors and factorials |
---|
| 301 | G4double R1, R2, Pnucl, Aeff; // Nuclear Parameters |
---|
| 302 | G4double InCoh; // Incoherent quasi-elastic cross section |
---|
| 303 | |
---|
| 304 | G4double binom(G4int N, G4int M) |
---|
| 305 | { |
---|
| 306 | if(N>170 && N>M && M) |
---|
| 307 | { |
---|
| 308 | G4double fN=N; |
---|
| 309 | G4double fM=M; |
---|
| 310 | G4double fD=N-M; |
---|
| 311 | G4double man=N*std::log(fN)-M*std::log(fM)-fD*std::log(fD)-1.; |
---|
| 312 | //G4cout<<"G4GCS::binom: N="<<N<<", M="<<M<<", C="<<man<<", E="<<std::exp(man)<<G4endl; |
---|
| 313 | return std::exp(man); |
---|
| 314 | } |
---|
| 315 | else if(N>1 && N>M && M) return Factorials[N]/Factorials[M]/Factorials[N-M]; |
---|
| 316 | return 1.; |
---|
| 317 | } |
---|
| 318 | |
---|
| 319 | void Initialize() |
---|
| 320 | { |
---|
| 321 | //static const G4double sat = 3.03; // Stirling saturation of Dubna |
---|
| 322 | G4double Val = Factorials[0] = Mnoj[0] = 1.; |
---|
| 323 | for (G4int i=1; i<nFact; i++) |
---|
| 324 | { |
---|
| 325 | G4double Si = 0.; |
---|
| 326 | Val *= i; // Calculation of the factorial values |
---|
| 327 | Factorials[i] = Val; |
---|
| 328 | for (G4int l = 0; l<=i; l++) |
---|
| 329 | { |
---|
| 330 | G4double Sm = 1.; // sum of reversed binom coefficients (?) |
---|
| 331 | for (G4int m=1; m<=i-l; m++) Sm += 1./binom(i-l,m); |
---|
| 332 | //if(i>169)G4cout<<"G4GCS::Init:S="<<Sm<<", b("<<i<<","<<l<<")="<<binom(i,l)<<G4endl; |
---|
| 333 | Si += Sm/binom(i,l); |
---|
| 334 | } // End of l LOOP |
---|
| 335 | Mnoj[i] = Si; |
---|
| 336 | //G4double d=Si-sat; |
---|
| 337 | //G4double r=d/sat; |
---|
| 338 | //if(i>100)G4cout<<"G4GCS::Init:i="<<i<<":"<<Si<<"-"<<sat<<" = "<<d<<", r="<<r<<G4endl; |
---|
| 339 | } // End of i LOOP |
---|
| 340 | } // End of Initialization |
---|
| 341 | |
---|
| 342 | // Dubna paranmeters of nuclei |
---|
| 343 | void CalculateNuclearParameters(G4int Anuc) // R1,R2,Pnuc,Aeff |
---|
| 344 | {// ============================================== |
---|
| 345 | G4double A=Anuc; |
---|
| 346 | // Personal corrections fore some nuclei |
---|
| 347 | if(Anuc == 4) // Personal correction for He (!) |
---|
| 348 | { |
---|
| 349 | R1 = 5.5; |
---|
| 350 | R2 = 3.7; |
---|
| 351 | Pnucl = 0.4; |
---|
| 352 | Aeff = 0.87; |
---|
| 353 | } |
---|
| 354 | else if(Anuc == 9) // Personal correction for Be (!) |
---|
| 355 | { |
---|
| 356 | R1 = 9.0; |
---|
| 357 | R2 = 7.0; |
---|
| 358 | Pnucl = 0.190; |
---|
| 359 | Aeff = 0.9; |
---|
| 360 | } |
---|
| 361 | //else if(Anuc == 11) // Personal correction for B (!) |
---|
| 362 | //{ |
---|
| 363 | // R1 = 10.8; |
---|
| 364 | // R2 = 7.5; |
---|
| 365 | // Pnucl = 0.85; |
---|
| 366 | // Aeff = 1.2; |
---|
| 367 | //} |
---|
| 368 | //else if(Anuc == 12) // Personal correction for C (!) |
---|
| 369 | //{ |
---|
| 370 | // R1 = 9.336; |
---|
| 371 | // R2 = 5.63; |
---|
| 372 | // Pnucl = 0.197; |
---|
| 373 | // Aeff = 1.0; |
---|
| 374 | //} |
---|
| 375 | else if(Anuc == 16) // Personal correction for O (!) |
---|
| 376 | { |
---|
| 377 | R1 = 10.50; |
---|
| 378 | R2 = 5.5; |
---|
| 379 | Pnucl = 0.7; |
---|
| 380 | Aeff = 0.98; |
---|
| 381 | //R1 = 11.3; |
---|
| 382 | //R2 = 2.5; |
---|
| 383 | //Pnucl = 0.75; |
---|
| 384 | //Aeff = 0.9; |
---|
| 385 | } |
---|
| 386 | else if(Anuc == 58) // Personal correction for Ni (!) |
---|
| 387 | { |
---|
| 388 | R1 = 15.0; |
---|
| 389 | R2 = 9.9; |
---|
| 390 | Pnucl = 0.45; |
---|
| 391 | Aeff = 0.85; |
---|
| 392 | } |
---|
| 393 | else if(Anuc == 90) // Personal correction for Zr (!) |
---|
| 394 | { |
---|
| 395 | //R1 = 16.5; |
---|
| 396 | //R2 = 11.62; |
---|
| 397 | //Pnucl = 0.4; |
---|
| 398 | //Aeff = 0.9; |
---|
| 399 | R1 = 16.5; |
---|
| 400 | R2 = 11.62; |
---|
| 401 | Pnucl = 0.4; |
---|
| 402 | Aeff = 0.7; |
---|
| 403 | } |
---|
| 404 | else if(Anuc == 208) // Personal correction for Pb (!) |
---|
| 405 | { |
---|
| 406 | // R1 = 20.73; R2 = 15.74. |
---|
| 407 | //R1 = 4.1408*std::pow(A,0.3018); |
---|
| 408 | //R2 = 3.806*std::pow(A-10.068,0.2685); |
---|
| 409 | //Pnucl = 0.9; |
---|
| 410 | //Aeff = 1.1; |
---|
| 411 | R1 = 19.5; |
---|
| 412 | R2 = 15.74; |
---|
| 413 | Pnucl = 0.4; |
---|
| 414 | Aeff = 0.7; |
---|
| 415 | } |
---|
| 416 | else |
---|
| 417 | { |
---|
| 418 | R1 = 4.45*std::pow(A-1.,.309); // First diffractional radius |
---|
| 419 | if(Anuc == 28) R1 *= .955; // Personal correction for Si (?!) |
---|
| 420 | R2 = 2.3*std::pow(A,.36); // Second diffraction radius |
---|
| 421 | Pnucl = 0.176+(.00167+.00000869*Anuc)*Anuc; // Momentum of mean Fermi motion |
---|
| 422 | Aeff = 0.9; // Effective screaning |
---|
| 423 | } |
---|
| 424 | } |
---|
| 425 | |
---|
| 426 | // Independent transition method calculation. Has problems for A>95 (short cut). |
---|
| 427 | G4double DiffElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) // All in MeV |
---|
| 428 | {// ================================================================= |
---|
| 429 | static const G4double eps = 0.000001; // Accuracy of calculations |
---|
| 430 | static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 |
---|
| 431 | static const G4double piMG = pi/mb2G2; // PiTrans from GeV^-2 to mb |
---|
| 432 | static const G4double dR0 = 1.06+1.06; // Doubled nominal R0 |
---|
| 433 | G4double hMom = HadrMom/1000.; // Momentum (GeV, inputParam) |
---|
| 434 | G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
| 435 | G4double MassH2 = MassH*MassH; |
---|
| 436 | G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2); // Tot energy in GeV |
---|
| 437 | //G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl; |
---|
| 438 | if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3"); |
---|
| 439 | if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!"); |
---|
| 440 | //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)"); |
---|
| 441 | CalculateParameters(hPDG, HadrMom); |
---|
| 442 | CalculateIntegralCS(A, HadrEnergy); |
---|
| 443 | CalculateNuclearParameters(A); |
---|
| 444 | G4double Q2 = aQ2/1000/1000; // in GeV^2 |
---|
| 445 | //MassN = A*0.938; // @@ This is bad! |
---|
| 446 | G4double MassN = A*0.9315; // @@ Atomic Unit is bad too |
---|
| 447 | G4double MassN2 = MassN*MassN; |
---|
| 448 | G4double S = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s |
---|
| 449 | G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // CM energy of a hadron |
---|
| 450 | G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum |
---|
| 451 | |
---|
| 452 | G4double Stot = HadrTot*mb2G2; // in GeV^-2 |
---|
| 453 | G4double Bhad = HadrSlope; // in GeV^-2 |
---|
| 454 | G4double Asq = 1+HadrReIm*HadrReIm; // |M|^2/(ImM)^2 |
---|
| 455 | G4double Rho2 = std::sqrt(Asq); // M/ImM |
---|
| 456 | G4double R12 = R1*R1; |
---|
| 457 | G4double R22 = R2*R2; |
---|
| 458 | G4double R12B = R12+Bhad+Bhad; |
---|
| 459 | G4double R22B = R22+Bhad+Bhad; |
---|
| 460 | G4double R12Ap = R12+20; |
---|
| 461 | G4double R22Ap = R22+20; |
---|
| 462 | G4double R13Ap = R12*R1/R12Ap; |
---|
| 463 | G4double R23Ap = R22*R2/R22Ap*Pnucl; |
---|
| 464 | G4double R23dR13 = R23Ap/R13Ap; |
---|
| 465 | G4double R12Apd = 2./R12Ap; |
---|
| 466 | G4double R22Apd = 2./R22Ap; |
---|
| 467 | G4double R122Apd = (R12Apd+R22Apd)/2; |
---|
| 468 | G4double dR0H = dR0*HadrEnergy/4; |
---|
| 469 | G4double DDSec1p = DDSect2+DDSect3*std::log(dR0H/R1); |
---|
| 470 | G4double DDSec2p = DDSect2+DDSect3*std::log(dR0H/std::sqrt((R12+R22)/2)); |
---|
| 471 | G4double DDSec3p = DDSect2+DDSect3*std::log(dR0H/R2); |
---|
| 472 | G4double Norm = (R12*R1-Pnucl*R22*R2)*Aeff; // Some questionable norming |
---|
| 473 | G4double R13 = R12*R1/R12B; |
---|
| 474 | G4double R23 = Pnucl*R22*R2/R22B; |
---|
| 475 | G4double norFac = Stot/twopi/Norm; // totCS (in GeV^-2) is used |
---|
| 476 | G4double Unucl = norFac*R13; |
---|
| 477 | G4double UnuclScr = norFac*R13Ap; |
---|
| 478 | G4double SinFi = HadrReIm/Rho2; |
---|
| 479 | G4double FiH = std::asin(SinFi); |
---|
| 480 | G4double N = -1.; |
---|
| 481 | G4double N2 = R23/R13; |
---|
| 482 | G4double ImElA0 = 0.; |
---|
| 483 | G4double ReElA0 = 0.; |
---|
| 484 | G4double Tot1 = 0.; |
---|
| 485 | for(G4int i=1; i<=A; i++) // @@ Make separately for n and p |
---|
| 486 | { |
---|
| 487 | N *= -Unucl*(A-i+1)*Rho2/i; |
---|
| 488 | G4double N4 = 1.; |
---|
| 489 | G4double Prod1 = std::exp(-Q2*R12B/i/4)*R12B/i; |
---|
| 490 | G4double medTot = R12B/i; |
---|
| 491 | for(G4int l=1; l<=i; l++) |
---|
| 492 | { |
---|
| 493 | G4double exp1 = l/R22B+(i-l)/R12B; |
---|
| 494 | N4 *= -(i-l+1)*N2/l; |
---|
| 495 | G4double Nexp = N4/exp1; |
---|
| 496 | Prod1 += Nexp*std::exp(-Q2/exp1/4); |
---|
| 497 | medTot += Nexp; |
---|
| 498 | } |
---|
| 499 | G4double FiHi = FiH*i; |
---|
| 500 | G4double nCos = N*std::cos(FiHi); |
---|
| 501 | ReElA0 += Prod1*N*std::sin(FiHi); |
---|
| 502 | ImElA0 += Prod1*nCos; |
---|
| 503 | Tot1 += medTot*nCos; |
---|
| 504 | if(std::abs(Prod1*N/ImElA0) < eps) break; |
---|
| 505 | } |
---|
| 506 | ImElA0 *= piMG; // The amplitude in mb |
---|
| 507 | ReElA0 *= piMG; // The amplitude in mb |
---|
| 508 | Tot1 *= piMG+piMG; |
---|
| 509 | G4double N1p = 1.; |
---|
| 510 | G4double R13Ap1 = DDSec1p*R13Ap*R13Ap/2; |
---|
| 511 | G4double R22Ap1 = DDSec2p*R13Ap*R23Ap; |
---|
| 512 | G4double R23Ap1 = DDSec3p*R23Ap*R23Ap*R22Ap/2; |
---|
| 513 | G4double R13Ap2 = R13Ap1*R12Ap/4; |
---|
| 514 | G4double R22Ap2 = R22Ap1/R122Apd/2; |
---|
| 515 | G4double R23Ap2 = R23Ap1*R22Ap/4; |
---|
| 516 | G4double Q28 = Q2/8; |
---|
| 517 | G4double Q24 = Q28+Q28; |
---|
| 518 | G4double Din1 = R13Ap2*std::exp(-Q28*R12Ap)-R22Ap2*std::exp(-(Q28+Q28)/R122Apd)+ |
---|
| 519 | R23Ap2*std::exp(-Q28*R22Ap); // i=0 start value |
---|
| 520 | G4double DTot1 = R13Ap2-R22Ap2+R23Ap2; // i=0 start value |
---|
| 521 | if (A<96) |
---|
| 522 | { |
---|
| 523 | for(G4int i=1; i<A-1; i++) |
---|
| 524 | { |
---|
| 525 | N1p *= -UnuclScr*(A-i-1)/i*Rho2; |
---|
| 526 | G4double N2p = 1.; |
---|
| 527 | G4double Din2 = 0.; |
---|
| 528 | G4double DmedTot = 0.; |
---|
| 529 | G4double BinCoeff = 1.; // Start for Binomial Coefficient |
---|
| 530 | for(G4int l = 0; l<=i; l++) |
---|
| 531 | { |
---|
| 532 | if(l) BinCoeff *= (i-l+1)/l; |
---|
| 533 | G4double exp1 = l/R22B+(i-l)/R12B; |
---|
| 534 | G4double exp1p = exp1+R12Apd; |
---|
| 535 | G4double exp2p = exp1+R122Apd; |
---|
| 536 | G4double exp3p = exp1+R22Apd; |
---|
| 537 | G4double NBinC = N2p*BinCoeff; |
---|
| 538 | Din2 += NBinC*(R13Ap1/exp1p*std::exp(-Q24/exp1p)- |
---|
| 539 | R22Ap1/exp2p*std::exp(-Q24/exp2p)+ |
---|
| 540 | R23Ap1/exp3p*std::exp(-Q24/exp3p)); |
---|
| 541 | DmedTot += NBinC * (R13Ap1/exp1p - R22Ap1/exp2p + R23Ap1/exp3p); |
---|
| 542 | N2p *= -R23dR13; |
---|
| 543 | } |
---|
| 544 | G4double comFac = N1p*Mnoj[i]/(i+2)/(i+1)*std::cos(FiH*i); |
---|
| 545 | Din1 += Din2*comFac; |
---|
| 546 | DTot1 += DmedTot*comFac; |
---|
| 547 | if(std::abs(Din2*N1p/Din1) < eps) break; |
---|
| 548 | } |
---|
| 549 | G4double cFac = A*(A-1)*4/Norm/Norm; |
---|
| 550 | Din1 *= -cFac; |
---|
| 551 | DTot1 *= cFac; |
---|
| 552 | } |
---|
| 553 | //else Din1=0.; |
---|
| 554 | |
---|
| 555 | G4double Corr0 = Tot00/Tot1; |
---|
| 556 | ImElA0 *= Corr0; |
---|
| 557 | |
---|
| 558 | return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*CMMom*CMMom*mb2G2/4/pi2; |
---|
| 559 | } // End of DiffElasticCS |
---|
| 560 | |
---|
| 561 | G4double CHIPSDiffElCS(G4int hPDG, G4double HadrMom, G4int Z, G4int A, G4double aQ2) // MeV |
---|
| 562 | {// ================================================================= |
---|
| 563 | static const G4double eps = 0.000001; // Accuracy of calculations |
---|
| 564 | static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 |
---|
| 565 | static const G4double piMG = pi/mb2G2; // PiTrans from GeV^-2 to mb |
---|
| 566 | G4double hMom = HadrMom/1000.; // Momentum (GeV, inputParam) |
---|
| 567 | G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
| 568 | G4double MassH2 = MassH*MassH; |
---|
| 569 | G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2); // Tot energy in GeV |
---|
| 570 | //G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl; |
---|
| 571 | if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3"); |
---|
| 572 | if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!"); |
---|
| 573 | //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)"); |
---|
| 574 | CalculateParameters(hPDG, HadrMom); |
---|
| 575 | CalculateIntegralCS(A, HadrEnergy); |
---|
| 576 | CalculateNuclearParameters(A); |
---|
| 577 | G4double Q2 = aQ2/1000/1000; // in GeV^2 |
---|
| 578 | //G4doubleMassN = A*0.938; // @@ This is bad! |
---|
| 579 | G4double MassN = A*0.9315; // @@ Atomic Unit is bad too |
---|
| 580 | if(G4NucleiPropertiesTable::IsInTable(Z,A)) |
---|
| 581 | MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.; // Geant4 NuclearMass |
---|
| 582 | G4double MassN2 = MassN*MassN; |
---|
| 583 | G4double S = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s |
---|
| 584 | G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // CM energy of a hadron |
---|
| 585 | G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum |
---|
| 586 | G4double Stot = HadrTot*mb2G2; // in GeV^-2 |
---|
| 587 | G4double Bhad = HadrSlope; // in GeV^-2 |
---|
| 588 | G4double Asq = 1+HadrReIm*HadrReIm; // |M|^2/(ImM)^2 |
---|
| 589 | G4double Rho2 = std::sqrt(Asq); // M/ImM |
---|
| 590 | G4double R12 = R1*R1; |
---|
| 591 | G4double R22 = R2*R2; |
---|
| 592 | G4double R12B = R12+Bhad+Bhad; // Slope is used |
---|
| 593 | G4double R22B = R22+Bhad+Bhad; // Slope is used |
---|
| 594 | G4double Norm = (R12*R1-Pnucl*R22*R2)*Aeff; // Some questionable norming |
---|
| 595 | G4double R13 = R12*R1/R12B; // Slope is used |
---|
| 596 | G4double R23 = Pnucl*R22*R2/R22B; // Slope is used |
---|
| 597 | G4double norFac = Stot/twopi/Norm; // totCS (in GeV^-2) is used |
---|
| 598 | G4double Unucl = norFac*R13; |
---|
| 599 | G4double SinFi = HadrReIm/Rho2; |
---|
| 600 | G4double FiH = std::asin(SinFi); |
---|
| 601 | G4double N = -1.; |
---|
| 602 | G4double N2 = R23/R13; // Slope is used |
---|
| 603 | G4double ImElA0 = 0.; |
---|
| 604 | G4double ReElA0 = 0.; |
---|
| 605 | G4double Tot1 = 0.; |
---|
| 606 | for(G4int i=1; i<=A; i++) |
---|
| 607 | { |
---|
| 608 | N *= -Unucl*(A-i+1)*Rho2/i; |
---|
| 609 | G4double N4 = 1.; |
---|
| 610 | G4double Prod1 = std::exp(-Q2*R12B/i/4)*R12B/i; // Slope is used |
---|
| 611 | G4double medTot = R12B/i; // Slope is used |
---|
| 612 | for(G4int l=1; l<=i; l++) |
---|
| 613 | { |
---|
| 614 | G4double exp1 = l/R22B+(i-l)/R12B; // Slope is used |
---|
| 615 | N4 *= -(i-l+1)*N2/l; // Slope is used |
---|
| 616 | G4double Nexp = N4/exp1; |
---|
| 617 | Prod1 += Nexp*std::exp(-Q2/exp1/4); |
---|
| 618 | medTot += Nexp; |
---|
| 619 | } |
---|
| 620 | G4double FiHi = FiH*i; |
---|
| 621 | G4double nCos = N*std::cos(FiHi); |
---|
| 622 | ReElA0 += Prod1*N*std::sin(FiHi); |
---|
| 623 | ImElA0 += Prod1*nCos; |
---|
| 624 | Tot1 += medTot*nCos; |
---|
| 625 | if(std::abs(Prod1*N/ImElA0) < eps) break; |
---|
| 626 | } |
---|
| 627 | ImElA0 *= piMG; // The amplitude in mb |
---|
| 628 | ReElA0 *= piMG; // The amplitude in mb |
---|
| 629 | Tot1 *= piMG+piMG; |
---|
| 630 | G4double Corr0 = Tot00/Tot1; |
---|
| 631 | ImElA0 *= Corr0; |
---|
| 632 | |
---|
| 633 | return (ReElA0*ReElA0+ImElA0*ImElA0)*CMMom*CMMom*mb2G2/4/pi2; |
---|
| 634 | } // End of DiffElasticCS |
---|
| 635 | |
---|
| 636 | G4double CHIPS_Tb(G4int A, G4double b) // T(b) in fm-2 |
---|
| 637 | { |
---|
| 638 | static G4int Am=0; // Associative memory value A |
---|
| 639 | static G4double B=0.; // Effective edge |
---|
| 640 | static G4double C=0.; // Normalization factor |
---|
| 641 | static G4double D=0.; // Slope parameter |
---|
| 642 | if(A!=Am) |
---|
| 643 | { |
---|
| 644 | B=.0008*A*A; // no units |
---|
| 645 | D=.42*std::pow(A,-.26); // fm^-2 |
---|
| 646 | C=A*D/pi/std::log(1.+B); // fm^-2 |
---|
| 647 | } |
---|
| 648 | G4double E=B*std::exp(-D*b*b); |
---|
| 649 | return C*E/(1.+E); // T(b) in fm^-2 |
---|
| 650 | } |
---|
| 651 | |
---|
| 652 | G4double Thickness(G4int A, G4double b, G4double R) // T(b) in fm^-2 |
---|
| 653 | {// ========================================== // rAmax=2.5 is a GlobStatConstPar |
---|
| 654 | static const G4double bTh = .545; // Surface thicknes |
---|
| 655 | static const G4double bTh2 = bTh*bTh; // SquaredSurface thicknes |
---|
| 656 | G4double b2 = b*b; // Squared impact parameter |
---|
| 657 | G4double dr = rAmax*R/(NptB-1); // Step of integration in z |
---|
| 658 | G4double R2 = R*R; // Working parameter |
---|
| 659 | G4double Norm = .75/pi/R2/R/(1.+bTh2*pi2/R2); // Corrected Volume (?) |
---|
| 660 | G4double r = -dr; // Pre value for the radius |
---|
| 661 | G4double SumZ = 0.; // Prototype of the result |
---|
| 662 | //G4double SumN = 0.; // is not used |
---|
| 663 | for(G4int k=0; k<NptB; k++) |
---|
| 664 | { |
---|
| 665 | r += dr; |
---|
| 666 | SumZ += 1./(1.+std::exp((std::sqrt(b2+r*r)-R)/bTh)); // Woods-Saxon density |
---|
| 667 | //SumN += r*r/(1.+std::exp((r-rAfm)/bTh)); // is not used |
---|
| 668 | } |
---|
| 669 | //SumN *= (twopi+twopi)*dr*Norm; // is not used |
---|
| 670 | return SumZ*(A+A)*dr*Norm; // T(b) in fm^-2 |
---|
| 671 | } // End of Thickness |
---|
| 672 | |
---|
| 673 | G4double CoherentDifElasticCS(G4int hPDG, G4double HadrMom, G4int A, G4double aQ2) //in MeV |
---|
| 674 | {// ========================================================================= |
---|
| 675 | //static const G4double eps = 0.000001; // CalculationAccuracy@@NotUsed |
---|
| 676 | //static const G4double mN = .9315; // Atomic Unit GeV |
---|
| 677 | static const G4double mN = .938; // Atomic Unit GeV |
---|
| 678 | static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 |
---|
| 679 | static const G4double incrf = 10; // Increase factor of radius for intergation |
---|
| 680 | static const G4double m2G10 = incrf*mb2G2; // Big conversion factor |
---|
| 681 | static const G4double s2G10 = std::sqrt(m2G10); // sqrt Big conversion factor |
---|
| 682 | G4double ReIntegrand[NptB],ImIntegrand[NptB],Thick[NptB]; // Calculated arrays |
---|
| 683 | G4QBesIKJY QI0(BessI0); // I0 Bessel function |
---|
| 684 | G4QBesIKJY QJ0(BessJ0); // J0 Bessel function |
---|
| 685 | G4double hMom = HadrMom/1000.; // Momentum (GeV, inputParam) |
---|
| 686 | G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
| 687 | G4double MassH2 = MassH*MassH; |
---|
| 688 | G4double HadrEnergy = std::sqrt(hMom*hMom+MassH2); // Tot energy in GeV |
---|
| 689 | if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3"); |
---|
| 690 | if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!"); |
---|
| 691 | //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)"); |
---|
| 692 | G4double Re = std::pow(A, 0.33333333); |
---|
| 693 | G4double Re2 = Re*Re; |
---|
| 694 | CalculateParameters(hPDG, HadrMom); |
---|
| 695 | G4double Q2 = aQ2/1000/1000; // GeV |
---|
| 696 | // Individual corrections for nuclei which had data |
---|
| 697 | if (A==208) r0 = 1.125; |
---|
| 698 | else if(A==90) r0 = 1.12; |
---|
| 699 | else if(A==64) r0 = 1.1; |
---|
| 700 | else if(A==58) r0 = 1.09; |
---|
| 701 | else if(A==48) r0 = 1.07; |
---|
| 702 | else if(A==40) r0 = 1.15; |
---|
| 703 | else if(A==28) r0 = 0.93; |
---|
| 704 | else if(A==16) r0 = 0.92; |
---|
| 705 | else if(A==12) r0 = 0.80; |
---|
| 706 | else r0 = 1.16*(1.-1.16/Re2); // For other nuclei which have not been tested |
---|
| 707 | G4double MassN = mN; // A * AtomicUnit(GeV) |
---|
| 708 | G4double MassN2 = MassN*MassN; |
---|
| 709 | G4double S = (MassN+MassN)*HadrEnergy+MassN2+MassH2; // Mondelstam S |
---|
| 710 | G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // Hadron CM Energy |
---|
| 711 | G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum |
---|
| 712 | G4double rAfm = r0*Re; // Nuclear radius |
---|
| 713 | G4double rAGeV = rAfm*s2G10; // Big integration radius |
---|
| 714 | G4double stepB = rAmax*rAGeV/(NptB-1); // dr step of integration |
---|
| 715 | G4double hTotG2 = HadrTot*mb2G2; |
---|
| 716 | G4double ReSum = 0.; |
---|
| 717 | G4double ImSum = 0.; |
---|
| 718 | G4double ValB = -stepB; |
---|
| 719 | for(G4int i=0; i<NptB; i++) |
---|
| 720 | { |
---|
| 721 | ValB += stepB; // An incident parameter |
---|
| 722 | G4double ValB2 = ValB*ValB; // A working value |
---|
| 723 | G4double IPH = ValB/HadrSlope; // A working value |
---|
| 724 | G4double dHS = HadrSlope+HadrSlope; // A working value |
---|
| 725 | G4double Integ = 0.; |
---|
| 726 | G4double ValS = 0.; |
---|
| 727 | for(G4int j=1; j<NptB; j++) |
---|
| 728 | { |
---|
| 729 | ValS += stepB; // Impact parameter GeV^-1 |
---|
| 730 | if(!i) Thick[j] = Thickness(A,ValS/s2G10,rAfm)/m2G10; // Calculate only once |
---|
| 731 | G4double FunS = IPH*ValS; |
---|
| 732 | if(FunS > 320.) break; |
---|
| 733 | Integ += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Thick[j]; |
---|
| 734 | } |
---|
| 735 | G4double InExp = -hTotG2*Integ*stepB/dHS; |
---|
| 736 | G4double expB = std::exp(InExp); |
---|
| 737 | G4double HRIE = HadrReIm*InExp; |
---|
| 738 | ReIntegrand[i] = (1.-expB*std::cos(HRIE)); |
---|
| 739 | ImIntegrand[i] = expB*std::sin(HRIE); |
---|
| 740 | } |
---|
| 741 | InCoh = 0.; // incohirent (quasi-elastic) |
---|
| 742 | ValB = -stepB; |
---|
| 743 | for(G4int k=0; k<NptB; k++) |
---|
| 744 | { |
---|
| 745 | ValB += stepB; |
---|
| 746 | InCoh += Thick[k]*ValB*std::exp(-hTotG2*Thick[k]); |
---|
| 747 | G4double J0qb = QJ0(std::sqrt(Q2)*ValB)*ValB; |
---|
| 748 | ReSum += J0qb*ReIntegrand[k]; |
---|
| 749 | ImSum += J0qb*ImIntegrand[k]; |
---|
| 750 | } |
---|
| 751 | //G4cout<<"GHAD:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2=" |
---|
| 752 | // <<Q2<<",m="<<mb2G2<<G4endl; |
---|
| 753 | InCoh *= stepB*hTotG2*hTotG2*(1.+HadrReIm*HadrReIm)*std::exp(-HadrSlope*Q2)/8/mb2G2; |
---|
| 754 | //G4cout<<"GHAD:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl; |
---|
| 755 | return (ReSum*ReSum+ImSum*ImSum)*m2G10*CMMom*CMMom*stepB*stepB/twopi; |
---|
| 756 | } |
---|
| 757 | |
---|
| 758 | G4double CHIPSDifElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2) |
---|
| 759 | {// ============================================================================ |
---|
| 760 | static const G4int Npb = 500; // A#of intergation points |
---|
| 761 | //static const G4double mN = .9315; // Atomic Unit GeV |
---|
| 762 | static const G4double mN = .938; // Atomic Unit GeV |
---|
| 763 | static const G4double hc2 = .3893793; // Transform from GeV^-2 to mb |
---|
| 764 | static const G4double mb2G2 = 1./hc2; // Transform from mb to GeV^-2 |
---|
| 765 | //static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 |
---|
| 766 | static const G4double f22mb = 10; // Transform from fermi^2 to mb |
---|
| 767 | static const G4double f22G2 = f22mb*mb2G2; // Transform from fm2 to GeV^-2 |
---|
| 768 | static const G4double f2Gm1 = std::sqrt(f22G2); // Transform from fm to GeV^-1 |
---|
| 769 | G4double Re = std::pow(A,.33333333); // A-dep coefficient |
---|
| 770 | G4double Lim = 50.*Re; // Integration accuracy limit |
---|
| 771 | G4double Tb[Npb]; // Calculated T(b) array |
---|
| 772 | G4QBesIKJY QI0(BessI0); // I0 Bessel function |
---|
| 773 | G4QBesIKJY QJ0(BessJ0); // J0 Bessel function |
---|
| 774 | hMom = hMom/1000.; // Momentum (GeV, inputParam) |
---|
| 775 | G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
| 776 | G4double MassH2 = MassH*MassH; // Squared mass of the hadron |
---|
| 777 | G4double hEnergy = std::sqrt(hMom*hMom+MassH2); // Tot energy in GeV |
---|
| 778 | G4double Q2 = aQ2/1000000.; // -t in GeV |
---|
| 779 | G4double MassN = mN; // AtomicUnit(GeV)[prototype] |
---|
| 780 | G4double MassN2 = MassN*MassN; // Squared mass of the target |
---|
| 781 | G4double S = (MassN+MassN)*hEnergy+MassN2+MassH2;// Mondelstam s |
---|
| 782 | G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // Hadron CM Energy |
---|
| 783 | G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum (to norm CS) |
---|
| 784 | // @@ Temporary only for nucleons |
---|
| 785 | //G4cout<<"CHPS:E="<<hEnergy<<",dM="<<2*MassN<<",sN="<<MassN2<<",sH="<<MassH2<<G4endl; |
---|
| 786 | G4double shNTot = 5.2+5.2*std::log(hEnergy)+51*std::pow(hEnergy,-.35); // SighN in mb |
---|
| 787 | G4double shNSl = 5.44+.88*std::log(S); // B-slope in GeV^-2 |
---|
| 788 | G4double shNReIm = .13*std::log(S/350)*std::pow(S,-.18); // Re/Im_hN in no unit |
---|
| 789 | //G4cout<<"CHPS: s="<<S<<",T="<<shNTot<<",R="<<shNReIm<<",B="<<shNSl<<G4endl; |
---|
| 790 | // @@ End of temporary ^^^^^^^ |
---|
| 791 | G4double dHS = shNSl+shNSl; // Working: doubled B-slope |
---|
| 792 | //G4double rAfm = 0.; |
---|
| 793 | //if (A==208) rAfm = 1.125*Re; |
---|
| 794 | //else if(A==90) rAfm = 1.12*Re; |
---|
| 795 | //else if(A==64) rAfm = 1.1*Re; |
---|
| 796 | //else if(A==58) rAfm = 1.09*Re; |
---|
| 797 | //else if(A==48) rAfm = 1.07*Re; |
---|
| 798 | //else if(A==40) rAfm = 1.15*Re; |
---|
| 799 | //else if(A==28) rAfm = 0.93*Re; |
---|
| 800 | //else if(A==16) rAfm = 0.92*Re; |
---|
| 801 | //else if(A==12) rAfm = 0.80*Re; |
---|
| 802 | //else rAfm = 1.16*(Re-1.16/Re); // For other nuclei |
---|
| 803 | //G4double stepB = 2.5*rAfm*f2Gm1/(Npb-1); // in GeV^-1, step of integral |
---|
| 804 | G4double stepB = (Re+Re+2.7)*f2Gm1/(Npb-1); // in GeV^-1, step of integral |
---|
| 805 | G4double hTotG2 = shNTot*mb2G2; // sigma_hN in GeV^-2 |
---|
| 806 | G4double ReSum = 0.; // Integration of RePart of Amp |
---|
| 807 | G4double ImSum = 0.; // Integration of ImPart of Amp |
---|
| 808 | G4double ValB = -stepB; |
---|
| 809 | for(G4int i=0; i<Npb; i++) |
---|
| 810 | { |
---|
| 811 | ValB += stepB; // An incident parameter |
---|
| 812 | G4double ValB2 = ValB*ValB; // A working value |
---|
| 813 | G4double IPH = ValB/shNSl; // A working value |
---|
| 814 | G4double Integ = 0.; // Integral over ImpactParam. |
---|
| 815 | G4double ValS = 0.; // Prototype of ImpactParameter |
---|
| 816 | for(G4int j=1; j<Npb; j++) |
---|
| 817 | { |
---|
| 818 | ValS += stepB; // back to fm // Impact parameter GeV^-1 |
---|
| 819 | if(!i) Tb[j] = CHIPS_Tb(A,ValS/f2Gm1)/f22G2; // GeV^2, calculate only once |
---|
| 820 | //if(!i) Tb[j] = Thickness(A,ValS/f2Gm1,rAfm)/f22G2; // Calculate T(b) only once |
---|
| 821 | G4double FunS = IPH*ValS; // Working product |
---|
| 822 | if(FunS > Lim) break; // (?) |
---|
| 823 | Integ += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Tb[j]; |
---|
| 824 | } |
---|
| 825 | G4double InExp = -hTotG2*Integ*stepB/dHS; // Working product |
---|
| 826 | G4double expB = std::exp(InExp); // Workung sqrt |
---|
| 827 | G4double HRIE = shNReIm*InExp; // Phase shift |
---|
| 828 | G4double J0qb = QJ0(std::sqrt(Q2)*ValB)*ValB; |
---|
| 829 | ReSum += J0qb*(1.-expB*std::cos(HRIE)); |
---|
| 830 | ImSum += J0qb*expB*std::sin(HRIE); |
---|
| 831 | } |
---|
| 832 | //G4cout<<"CHPS:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl; |
---|
| 833 | //return (ReSum*ReSum+ImSum*ImSum)*mb2G2*CMMom*CMMom*stepB*stepB/twopi; |
---|
| 834 | return (ReSum*ReSum+ImSum*ImSum)*f22G2*CMMom*CMMom*stepB*stepB/twopi; |
---|
| 835 | } |
---|
| 836 | G4double CHIPSDifQuasiElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2) |
---|
| 837 | {// ================================================================================= |
---|
| 838 | static const G4int Npb = 500; // A#of intergation points |
---|
| 839 | //static const G4double mN = .9315; // Atomic Unit GeV |
---|
| 840 | static const G4double mN = .938; // Mass of proton GeV |
---|
| 841 | static const G4double hc2 = .3893793; // Transform from GeV^-2 to mb |
---|
| 842 | static const G4double mb2G2 = 1./hc2; // Transform from mb to GeV^-2 |
---|
| 843 | //static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 |
---|
| 844 | static const G4double f22mb = 10; // Transform from fermi^2 to mb |
---|
| 845 | static const G4double f22G2 = f22mb*mb2G2; // Transform from fm2 to GeV^-2 |
---|
| 846 | static const G4double f2Gm1 = std::sqrt(f22G2); // Transform from fm to GeV^-1 |
---|
| 847 | hMom = hMom/1000.; // Momentum (GeV, inputParam) |
---|
| 848 | G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV |
---|
| 849 | G4double MassH2 = MassH*MassH; // Squared mass of the hadron |
---|
| 850 | G4double hEnergy = std::sqrt(hMom*hMom+MassH2); // Tot energy in GeV |
---|
| 851 | G4double Q2 = aQ2/1000000.; // -t in GeV |
---|
| 852 | G4double MassN = mN; // A*AtomicUnit(GeV)[prototype] |
---|
| 853 | //if(G4NucleiPropertiesTable::IsInTable(Z,A)) |
---|
| 854 | // MassN=G4NucleiProperties::GetNuclearMass(A,Z)/A/1000.; // Geant4 NuclearMass/A |
---|
| 855 | G4double MassN2 = MassN*MassN; // Squared mass of the target |
---|
| 856 | G4double S = (MassN+MassN)*hEnergy+MassN2+MassH2;// Mondelstam s |
---|
| 857 | // @@ Temporary only for nucleons |
---|
| 858 | G4double shNTot = 5.2+5.2*std::log(hEnergy)+51*std::pow(hEnergy,-.35); // SighN in mb |
---|
| 859 | G4double shNSl = 6.44+.88*std::log(S); // B-slope in GeV^-2 |
---|
| 860 | G4double shNReIm = .13*std::log(S/350)*std::pow(S,-.18); // Re/Im_hN in no unit |
---|
| 861 | // @@ End of temporary ^^^^^^^ |
---|
| 862 | G4double Re = std::pow(A,.33333333); // A-dep coefficient |
---|
| 863 | G4double stepB = (Re+Re+2.7)*f2Gm1/(Npb-1); // in GeV^-1, step of integral |
---|
| 864 | G4double InQE = 0.; // Quasielastic integral |
---|
| 865 | G4double ValB = -stepB; // Impact parameter |
---|
| 866 | G4double hTotG2 = shNTot*mb2G2; // sigma_hN in GeV^-2 |
---|
| 867 | for(G4int k=0; k<Npb; k++) |
---|
| 868 | { |
---|
| 869 | ValB += stepB; |
---|
| 870 | //G4double Tb = Thickness(A,ValB/f2Gm1,rAfm)/f22G2; // Calculate T(b) only once |
---|
| 871 | G4double Tb = CHIPS_Tb(A,ValB/f2Gm1)/f22G2; // GeV^2, calculate only once |
---|
| 872 | InQE += Tb*ValB*std::exp(-hTotG2*Tb); |
---|
| 873 | } |
---|
| 874 | //G4cout<<"CHPS:st="<<stepB<<",hT="<<hTotG2<<",HRI="<<HadrReIm<<",HS="<<HadrSlope<<",Q2=" |
---|
| 875 | // <<Q2<<",m="<<mb2G2<<G4endl; |
---|
| 876 | return InQE*stepB*hTotG2*hTotG2*(1.+shNReIm*shNReIm)*std::exp(-shNSl*Q2)/8/mb2G2; |
---|
| 877 | } |
---|
| 878 | #endif |
---|
| 879 | |
---|
| 880 | int main() |
---|
| 881 | { |
---|
| 882 | #ifdef bestest |
---|
| 883 | // *** Test of CHIPS implementation of Bessel functions (accuracy 1.E-8 and better) *** |
---|
| 884 | const G4int nj=21; |
---|
| 885 | const G4int ny=16; |
---|
| 886 | const G4int ni=20; |
---|
| 887 | const G4int nk=20; |
---|
| 888 | const G4double jX[nj]= |
---|
| 889 | {-5.,-4.,-3.,-2.,-1.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.}; |
---|
| 890 | const G4double yX[ny]={.1,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.}; |
---|
| 891 | const G4double iX[ni]= |
---|
| 892 | {0.,.2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,6.,8.,10.}; |
---|
| 893 | const G4double kX[nk]= |
---|
| 894 | {.1,.2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,6.,8.,10.}; |
---|
| 895 | const G4double vBJ0[nj]= |
---|
| 896 | {-.1775968,-.3971498,-.2600520,0.2238908,0.7651976,1.0000000,0.7651977, |
---|
| 897 | 0.2238908,-.2600520,-.3971498,-.1775968,0.1506453,0.3000793,0.1716508, |
---|
| 898 | -.0903336,-.2459358,-.1711903,0.0476893,0.2069261,0.1710735,-.0142245}; |
---|
| 899 | const G4double vBJ1[nj]= |
---|
| 900 | {0.3275791,0.0660433,-.3390590,-.5767248,-.4400506,0.0000000,0.4400506, |
---|
| 901 | 0.5767248,0.3390590,-.0660433,-.3275791,-.2766839,-.0046828,0.2346364, |
---|
| 902 | 0.2453118,0.0434728,-.1767853,-.2234471,-.0703181,0.1333752,0.2051040}; |
---|
| 903 | const G4double vBY0[ny]= |
---|
| 904 | {-1.534239,0.0882570,.51037567,.37685001,-.0169407,-.3085176,-.2881947,-.0259497, |
---|
| 905 | 0.2235215,0.2499367,0.0556712,-.1688473,-.2252373,-.0782079,0.1271926,0.2054743}; |
---|
| 906 | const G4double vBY1[ny]= |
---|
| 907 | {-6.458951,-.7812128,-.1070324,0.3246744,0.3979257,0.1478631,-.1750103,-.3026672, |
---|
| 908 | -.1580605,0.1043146,0.2490154,0.1637055,-.0570992,-.2100814,-.1666448,0.0210736}; |
---|
| 909 | const G4double vBI0[ni]={1.0000000,1.0100250,1.0404018,1.0920453,1.1665149, |
---|
| 910 | 1.2660658,1.3937256,1.5533951,1.7499807,1.9895593, |
---|
| 911 | 2.2795852,3.2898391,4.8807925,7.3782035,11.301922, |
---|
| 912 | 17.481172,27.239871,67.234406,427.56411,2815.7167}; |
---|
| 913 | const G4double vBI1[ni]={.00000000,.10050083,.20402675,.31370403,.43286480, |
---|
| 914 | .56515912,.71467794,.88609197,1.0848107,1.3171674, |
---|
| 915 | 1.5906369,2.5167163,3.9533700,6.2058350,9.7594652, |
---|
| 916 | 15.389221,24.335643,61.341937,399.87313,2670.9883}; |
---|
| 917 | const G4double vBK0[nk]={2.4270690,1.7527038,1.1145291,.77752208,.56534710, |
---|
| 918 | .42102445,.31850821,.24365506,.18795475,.14593140, |
---|
| 919 | .11389387,.062347553,.0347395,.019598897,.011159676, |
---|
| 920 | .0063998572,.0036910983,.0012439943,.00014647071,.000017780062}; |
---|
| 921 | const G4double vBK1[nk]={9.8538451,4.7759725,2.1843544,1.3028349,.86178163, |
---|
| 922 | .60190724,.43459241,.32083589,.24063392,.18262309, |
---|
| 923 | .13986588,.073890816,.040156431,.022239393,.012483499, |
---|
| 924 | .0070780949,.0040446134,.0013439197,.00015536921,.000018648773}; |
---|
| 925 | //Bessel J0('J',0); // CLHEP J/Y is possible insead of 'J' |
---|
| 926 | //Bessel J1('J',1); |
---|
| 927 | G4QBesIKJY QI0(BessI0); // CHIPS I/K/J/Y 0/1 Bessel functions |
---|
| 928 | G4QBesIKJY QI1(BessI1); |
---|
| 929 | G4QBesIKJY QJ0(BessJ0); |
---|
| 930 | G4QBesIKJY QJ1(BessJ1); |
---|
| 931 | G4QBesIKJY QK0(BessK0); |
---|
| 932 | G4QBesIKJY QK1(BessK1); |
---|
| 933 | G4QBesIKJY QY0(BessY0); |
---|
| 934 | G4QBesIKJY QY1(BessY1); |
---|
| 935 | G4cout<<"G4GCS: ***J0*** Test: ================================================"<<G4endl; |
---|
| 936 | for(G4int j0=0; j0<nj; j0++) |
---|
| 937 | { |
---|
| 938 | G4double F=QJ0(jX[j0]); |
---|
| 939 | G4double d=F-vBJ0[j0]; |
---|
| 940 | G4double r=std::abs(d/F); |
---|
| 941 | G4cout<<"G4GCS: x="<<jX[j0]<<", J0="<<F<<" - "<<vBJ0[j0]<<" = "<<d<<", r="<<r<<G4endl; |
---|
| 942 | } |
---|
| 943 | G4cout<<"G4GCS: ***J1*** Test: ================================================"<<G4endl; |
---|
| 944 | for(G4int j1=0; j1<nj; j1++) |
---|
| 945 | { |
---|
| 946 | G4double F=QJ1(jX[j1]); |
---|
| 947 | G4double d=F-vBJ1[j1]; |
---|
| 948 | G4double r=std::abs(d/F); |
---|
| 949 | G4cout<<"G4GCS: x="<<jX[j1]<<", J1="<<F<<" - "<<vBJ1[j1]<<" = "<<d<<", r="<<r<<G4endl; |
---|
| 950 | } |
---|
| 951 | G4cout<<"G4GCS: ***Y0*** Test: ================================================"<<G4endl; |
---|
| 952 | for(G4int y0=0; y0<ny; y0++) |
---|
| 953 | { |
---|
| 954 | G4double F=QY0(yX[y0]); |
---|
| 955 | G4double d=F-vBY0[y0]; |
---|
| 956 | G4double r=std::abs(d/F); |
---|
| 957 | G4cout<<"G4GCS: x="<<yX[y0]<<", Y0="<<F<<" - "<<vBY0[y0]<<" = "<<d<<", r="<<r<<G4endl; |
---|
| 958 | } |
---|
| 959 | G4cout<<"G4GCS: ***Y1*** Test: ================================================"<<G4endl; |
---|
| 960 | for(G4int y1=0; y1<ny; y1++) |
---|
| 961 | { |
---|
| 962 | G4double F=QY1(yX[y1]); |
---|
| 963 | G4double d=F-vBY1[y1]; |
---|
| 964 | G4double r=std::abs(d/F); |
---|
| 965 | G4cout<<"G4GCS: x="<<yX[y1]<<", Y1="<<F<<" - "<<vBY1[y1]<<" = "<<d<<", r="<<r<<G4endl; |
---|
| 966 | } |
---|
| 967 | G4cout<<"G4GCS: ***I0*** Test: ================================================"<<G4endl; |
---|
| 968 | for(G4int i0=0; i0<ni; i0++) |
---|
| 969 | { |
---|
| 970 | G4double F=QI0(iX[i0]); |
---|
| 971 | G4double d=F-vBI0[i0]; |
---|
| 972 | G4double r=std::abs(d/F); |
---|
| 973 | G4cout<<"G4GCS: x="<<iX[i0]<<", I0="<<F<<" - "<<vBI0[i0]<<" = "<<d<<", r="<<r<<G4endl; |
---|
| 974 | } |
---|
| 975 | G4cout<<"G4GCS: ***I1*** Test: ================================================"<<G4endl; |
---|
| 976 | for(G4int i1=0; i1<ni; i1++) |
---|
| 977 | { |
---|
| 978 | G4double F=QI1(iX[i1]); |
---|
| 979 | G4double d=F-vBI1[i1]; |
---|
| 980 | G4double r=std::abs(d/F); |
---|
| 981 | G4cout<<"G4GCS: x="<<iX[i1]<<", I1="<<F<<" - "<<vBI1[i1]<<" = "<<d<<", r="<<r<<G4endl; |
---|
| 982 | } |
---|
| 983 | G4cout<<"G4GCS: ***K0*** Test: ================================================"<<G4endl; |
---|
| 984 | for(G4int k0=0; k0<nk; k0++) |
---|
| 985 | { |
---|
| 986 | G4double F=QK0(kX[k0]); |
---|
| 987 | G4double d=F-vBK0[k0]; |
---|
| 988 | G4double r=std::abs(d/F); |
---|
| 989 | G4cout<<"G4GCS: x="<<kX[k0]<<", I1="<<F<<" - "<<vBK0[k0]<<" = "<<d<<", r="<<r<<G4endl; |
---|
| 990 | } |
---|
| 991 | G4cout<<"G4GCS: ***K1*** Test: ================================================"<<G4endl; |
---|
| 992 | for(G4int k1=0; k1<nk; k1++) |
---|
| 993 | { |
---|
| 994 | G4double F=QK1(kX[k1]); |
---|
| 995 | G4double d=F-vBK1[k1]; |
---|
| 996 | G4double r=std::abs(d/F); |
---|
| 997 | G4cout<<"G4GCS: x="<<kX[k1]<<", I1="<<F<<" - "<<vBK1[k1]<<" = "<<d<<", r="<<r<<G4endl; |
---|
| 998 | } |
---|
| 999 | G4cout<<"End .................................................................."<<G4endl; |
---|
| 1000 | #endif |
---|
| 1001 | |
---|
| 1002 | // Test of integrated cross sections |
---|
| 1003 | //const G4int na=12; |
---|
| 1004 | const G4int na=1; |
---|
| 1005 | // |
---|
| 1006 | const G4int np=7; |
---|
| 1007 | const G4int nm=1; |
---|
| 1008 | //// He Be C O Al Ti Ni Cu Sn Ta Pb U |
---|
| 1009 | //const G4int A[na]={4,9,12,16,27,48,58,64,120,181,207,238}; // A's of target nuclei |
---|
| 1010 | // He Al Pb |
---|
| 1011 | const G4int A[na]={208}; // A's of target nuclei |
---|
| 1012 | // p n pi+ pi- K+ K- antip |
---|
| 1013 | const G4int pdg[np]={2212,2112,211,-211,321,-321,-2212}; // projectiles |
---|
| 1014 | const G4double mom[nm]={1090.}; // momentum in MeV/c |
---|
| 1015 | #ifdef integrc |
---|
| 1016 | for(G4int ip=0; ip<np; ip++) |
---|
| 1017 | { |
---|
| 1018 | G4int PDG=pdg[ip]; |
---|
| 1019 | for(G4int im=0; im<nm; im++) |
---|
| 1020 | { |
---|
| 1021 | CalculateParameters(PDG, mom[im]); |
---|
| 1022 | G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<":A,Tot,Inel,Prod,El,QEl"<<G4endl; |
---|
| 1023 | for(G4int ia=0; ia<na; ia++) |
---|
| 1024 | { |
---|
| 1025 | CalculateIntegralCrossSections(A[ia]); |
---|
| 1026 | G4cout<<A[ia]<<" "<<TotalCrSec<<" "<<InelCrSec<<" "<<ProdCrSec |
---|
| 1027 | <<" "<<ElasticCrSec<<" "<<QuasyElasticCrSec<<G4endl; |
---|
| 1028 | } |
---|
| 1029 | } |
---|
| 1030 | } |
---|
| 1031 | #endif |
---|
| 1032 | |
---|
| 1033 | #ifdef eldiffcs |
---|
| 1034 | const G4double ms=.000001; |
---|
| 1035 | const G4int nt=200; |
---|
| 1036 | G4double t[nt]; // Q2=-t in Mev^2 |
---|
| 1037 | const G4double tMin=200.; |
---|
| 1038 | const G4double tMax=2000000.; |
---|
| 1039 | const G4double ltMi=std::log(tMin); |
---|
| 1040 | const G4double ltMa=std::log(tMax); |
---|
| 1041 | const G4double dlt=(ltMa-ltMi)/(nt-1); |
---|
| 1042 | G4double lt=ltMi-dlt; |
---|
| 1043 | for(G4int it=0; it<nt; it++) |
---|
| 1044 | { |
---|
| 1045 | lt+=dlt; |
---|
| 1046 | t[it]=std::exp(lt); |
---|
| 1047 | } |
---|
| 1048 | //// He Be C O Al Ti Ni Cu Sn Ta Pb U |
---|
| 1049 | //const G4int Z[na]={2,4, 6, 8,13,22,28,29, 50, 73, 82, 92}; // Z's of target nuclei |
---|
| 1050 | // He Al Pb |
---|
| 1051 | const G4int Z[na]={82}; // Z's of target nuclei |
---|
| 1052 | // Test of differential ellastic cross sections |
---|
| 1053 | Initialize(); |
---|
| 1054 | //for(G4int ip=0; ip<np; ip++) |
---|
| 1055 | for(G4int ip=0; ip<1; ip++) |
---|
| 1056 | { |
---|
| 1057 | G4int PDG=pdg[ip]; |
---|
| 1058 | for(G4int im=0; im<nm; im++) |
---|
| 1059 | { |
---|
| 1060 | CalculateParameters(PDG, mom[im]); |
---|
| 1061 | G4cout<<G4endl<<"G4GCS:PDG="<<PDG<<",P="<<mom[im]<<G4endl; |
---|
| 1062 | for(G4int ia=0; ia<na; ia++) |
---|
| 1063 | { |
---|
| 1064 | G4cout<<"G4GCS:-------------------A="<<A[ia]<<G4endl; |
---|
| 1065 | for(G4int it=0; it<nt; it++) |
---|
| 1066 | { |
---|
| 1067 | G4double Sig1 = CHIPSDiffElCS(PDG, mom[im], Z[ia], A[ia], t[it]); |
---|
| 1068 | //G4double Sig1 = DiffElasticCS(PDG, mom[im], A[ia], t[it]); |
---|
| 1069 | G4double Sig2 = CoherentDifElasticCS(PDG, mom[im], A[ia], t[it]); |
---|
| 1070 | G4double Sig3 = CHIPSDifElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]); |
---|
| 1071 | G4double CQEl = CHIPSDifQuasiElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]); |
---|
| 1072 | G4cout<<std::setw(12)<<ms*t[it]<<std::setw(12)<<Sig1<<std::setw(12)<<Sig2 |
---|
| 1073 | <<std::setw(12)<<InCoh<<std::setw(12)<<Sig3<<std::setw(12)<<CQEl<<G4endl; |
---|
| 1074 | } |
---|
| 1075 | } |
---|
| 1076 | } |
---|
| 1077 | } |
---|
| 1078 | #endif |
---|
| 1079 | return 0; |
---|
| 1080 | } |
---|