Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 24.04.20 V. Grichine 27 // 28 // (nu_e,anti_nu_e)-nucleus xsc 29 30 31 32 #include "G4ElNeutrinoNucleusTotXsc.hh" 33 #include "G4PhysicalConstants.hh" 34 #include "G4SystemOfUnits.hh" 35 #include "G4DynamicParticle.hh" 36 #include "G4ParticleTable.hh" 37 #include "G4IonTable.hh" 38 #include "G4HadTmpUtil.hh" 39 #include "G4NistManager.hh" 40 #include "G4Material.hh" 41 #include "G4Element.hh" 42 #include "G4Isotope.hh" 43 #include "G4ElementVector.hh" 44 45 #include "G4Electron.hh" 46 #include "G4Positron.hh" 47 48 using namespace std; 49 using namespace CLHEP; 50 51 G4ElNeutrinoNucleusTotXsc::G4ElNeutrinoNucleus 52 : G4VCrossSectionDataSet("NuElNuclTotXsc") 53 { 54 fCofXsc = 1.e-38*cm2/GeV; 55 56 // G4cout<<"fCofXsc = "<<fCofXsc*GeV/cm2<<" 57 58 // PDG2016: sin^2 theta Weinberg 59 60 fSin2tW = 0.23129; // 0.2312; 61 62 // 9 <-> 6, 5/9 or 5/6 ? 63 64 fCofS = 5.*fSin2tW*fSin2tW/9.; 65 fCofL = 1. - fSin2tW + fCofS; 66 67 // G4cout<<"fCosL = "<<fCofL<<", fCofS = "<< 68 69 fCutEnergy = 0.; // default value 70 71 fBiasingFactor = 1.; // default as physics 72 73 fIndex = 50; 74 75 fTotXsc = 0.; 76 fCcTotRatio = 0.75; // from nc/cc~0.33 ratio 77 fCcFactor = fNcFactor = 1.; 78 79 theElectron = G4Electron::Electron(); 80 thePositron = G4Positron::Positron(); 81 } 82 83 G4ElNeutrinoNucleusTotXsc::~G4ElNeutrinoNucleu 84 {} 85 86 ////////////////////////////////////////////// 87 88 89 G4bool 90 G4ElNeutrinoNucleusTotXsc::IsIsoApplicable( co 91 { 92 G4bool result = false; 93 G4String pName = aPart->GetDefinition()->Get 94 95 if( pName == "nu_e" || pName == "anti 96 { 97 result = true; 98 } 99 return result; 100 } 101 102 ////////////////////////////////////// 103 104 G4double G4ElNeutrinoNucleusTotXsc::GetElement 105 G4int Z, const G4Material* 106 { 107 G4int Zi(0); 108 size_t i(0), j(0); 109 const G4ElementVector* theElementVector = ma 110 111 for ( i = 0; i < theElementVector->size(); + 112 { 113 Zi = (*theElementVector)[i]->GetZasInt(); 114 if( Zi == Z ) break; 115 } 116 const G4Element* elm = (*theElementVector)[i 117 size_t nIso = elm->GetNumberOfIsotopes(); 118 G4double fact = 0.0; 119 G4double xsec = 0.0; 120 const G4Isotope* iso = nullptr; 121 const G4IsotopeVector* isoVector = elm->GetI 122 const G4double* abundVector = elm->GetRelati 123 124 for (j = 0; j<nIso; ++j) 125 { 126 iso = (*isoVector)[j]; 127 G4int A = iso->GetN(); 128 129 if( abundVector[j] > 0.0 && IsIsoApplicabl 130 { 131 fact += abundVector[j]; 132 xsec += abundVector[j]*GetIsoCrossSectio 133 } 134 } 135 if( fact > 0.0) { xsec /= fact; } 136 return xsec; 137 } 138 139 140 ////////////////////////////////////////////// 141 // 142 // 143 144 G4double G4ElNeutrinoNucleusTotXsc::GetIsoCros 145 const G4Isotope*, const G4Element* 146 { 147 fCcFactor = fNcFactor = 1.; 148 fCcTotRatio = 0.25; 149 150 G4double ccnuXsc, ccanuXsc, ncXsc, totXsc(0. 151 152 G4double energy = aPart->GetTotalEnergy(); 153 G4String pName = aPart->GetDefinition()->G 154 155 G4int index = GetEnergyIndex(energy); 156 157 if( index >= fIndex ) 158 { 159 G4double pm = proton_mass_c2; 160 G4double s2 = 2.*energy*pm+pm*pm; 161 G4double aa = 1.; 162 G4double bb = 1.085; 163 G4double mw = 80.385*GeV; 164 fCcFactor = bb/(1.+ aa*s2/mw/mw); 165 166 G4double mz = 91.1876*GeV; 167 fNcFactor = bb/(1.+ aa*s2/mz/mz); 168 } 169 ccnuXsc = GetNuElTotCsXsc(index, energy); 170 ccnuXsc *= fCcFactor; 171 ccanuXsc = GetANuElTotCsXsc(index, energy); 172 ccanuXsc *= fCcFactor; 173 174 if( pName == "nu_e") 175 { 176 ncXsc = fCofL*ccnuXsc + fCofS*ccanuXsc; 177 ncXsc *= fNcFactor/fCcFactor; 178 totXsc = ccnuXsc + ncXsc; 179 if( totXsc > 0.) fCcTotRatio = ccnuXsc/tot 180 } 181 else if( pName == "anti_nu_e") 182 { 183 ncXsc = fCofL*ccanuXsc + fCofS*ccnuXsc; 184 ncXsc *= fNcFactor/fCcFactor; 185 totXsc = ccanuXsc + ncXsc; 186 if( totXsc > 0.) fCcTotRatio = ccanuXsc/to 187 } 188 else return totXsc; 189 190 totXsc *= fCofXsc; 191 totXsc *= energy; 192 totXsc *= A; // incoherent sum over all is 193 194 totXsc *= fBiasingFactor; // biasing up, if 195 196 fTotXsc = totXsc; 197 198 return totXsc; 199 } 200 201 ////////////////////////////////////////////// 202 // 203 // Return index of nu/anu energy array corresp 204 205 G4int G4ElNeutrinoNucleusTotXsc::GetEnergyInde 206 { 207 G4int i, eIndex = 0; 208 209 for( i = 0; i < fIndex; i++) 210 { 211 if( energy <= fNuElEnergy[i]*GeV ) 212 { 213 eIndex = i; 214 break; 215 } 216 } 217 if( i >= fIndex-1 ) eIndex = fIndex-1; 218 // G4cout<<"eIndex = "<<eIndex<<G4endl; 219 return eIndex; 220 } 221 222 ////////////////////////////////////////////// 223 // 224 // nu_e xsc for index-1, index linear over ene 225 226 G4double G4ElNeutrinoNucleusTotXsc::GetNuElTot 227 { 228 G4double xsc(0.); 229 230 if( index <= 0 || energy < theElectron->GetP 231 else if (index >= fIndex) xsc = fNuElTotXsc[ 232 else 233 { 234 G4double x1 = fNuElEnergy[index-1]*GeV; 235 G4double x2 = fNuElEnergy[index]*GeV; 236 G4double y1 = fNuElTotXsc[index-1]; 237 G4double y2 = fNuElTotXsc[index]; 238 239 if(x1 >= x2) return fNuElTotXsc[index]; 240 else 241 { 242 G4double angle = (y2-y1)/(x2-x1); 243 xsc = y1 + (energy-x1)*angle; 244 } 245 } 246 return xsc; 247 } 248 249 ////////////////////////////////////////////// 250 // 251 // anu_e xsc for index-1, index linear over en 252 253 G4double G4ElNeutrinoNucleusTotXsc::GetANuElTo 254 { 255 G4double xsc(0.); 256 257 if( index <= 0 || energy < thePositron->GetP 258 else if (index >= fIndex) xsc = fANuElTotXsc 259 else 260 { 261 G4double x1 = fNuElEnergy[index-1]*GeV; 262 G4double x2 = fNuElEnergy[index]*GeV; 263 G4double y1 = fANuElTotXsc[index-1]; 264 G4double y2 = fANuElTotXsc[index]; 265 266 if( x1 >= x2 ) return fANuElTotXsc[index]; 267 else 268 { 269 G4double angle = (y2-y1)/(x2-x1); 270 xsc = y1 + (energy-x1)*angle; 271 } 272 } 273 return xsc; 274 } 275 276 ////////////////////////////////////////////// 277 // 278 // return fNuElTotXsc[index] if the index is i 279 280 G4double G4ElNeutrinoNucleusTotXsc::GetNuElTot 281 { 282 if( index >= 0 && index < fIndex) return fNu 283 else 284 { 285 G4cout<<"Improper index of fNuElTotXsc arr 286 return 0.; 287 } 288 } 289 290 ////////////////////////////////////////////// 291 // 292 // return fANuElTotXsc[index] if the index is 293 294 G4double G4ElNeutrinoNucleusTotXsc::GetANuElTo 295 { 296 if( index >= 0 && index < fIndex) return fAN 297 else 298 { 299 G4cout<<"Improper index of fANuElTotXsc ar 300 return 0.; 301 } 302 } 303 304 305 ////////////////////////////////////////////// 306 // 307 // E_nu in GeV 308 309 const G4double G4ElNeutrinoNucleusTotXsc::fNuE 310 { 311 0.000561138, 0.000735091, 0.000962969, 0.001 312 0.00216484, 0.00283594, 0.00371508, 0.004866 313 0.00835185, 0.0109409, 0.0143326, 0.0187757, 314 0.032221, 0.0422095, 0.0552945, 0.0724358, 0 315 0.124307, 0.162842, 0.213323, 0.279453, 0.36 316 0.47957, 0.628237, 0.82299, 1.07812, 1.41233 317 1.85016, 2.42371, 3.17505, 4.15932, 5.44871, 318 7.13781, 9.35053, 12.2492, 16.0464, 21.0208, 319 27.5373, 36.0739, 47.2568, 61.9064, 81.0973, 320 106.238, 139.171, 182.314, 238.832, 312.869 321 }; 322 323 ////////////////////////////////////////////// 324 // 325 // nu_e CC xsc_tot/E_nu, in 10^-38 cm2/GeV 326 327 const G4double G4ElNeutrinoNucleusTotXsc::fNuE 328 { 329 0.0026484, 0.00609503, 0.00939421, 0.0132163 330 0.0237692, 0.0312066, 0.0406632, 0.0526867, 331 0.0871913, 0.111359, 0.141458, 0.178584, 0.2 332 0.27822, 0.342461, 0.416865, 0.501361, 0.596 333 0.713623, 0.905749, 1.20718, 1.52521, 1.7528 334 1.82072, 1.67119, 1.50074, 1.3077, 1.14923, 335 1.0577, 0.977911, 0.918526, 0.792889, 0.7022 336 0.678615, 0.687099, 0.725167, 0.706795, 0.67 337 0.649791, 0.651328, 0.651934, 0.658062, 0.66 338 0.662534, 0.662601, 0.660261, 0.656724, 0.65 339 }; 340 341 342 343 ////////////////////////////////////////////// 344 // 345 // anu_e CC xsc_tot/E_nu, in 10^-38 cm2/GeV 346 347 const G4double G4ElNeutrinoNucleusTotXsc::fANu 348 { 349 0.00103385, 0.00237807, 0.00366358, 0.005151 350 0.00925859, 0.0121508, 0.0158252, 0.0204908, 351 0.0338304, 0.0431234, 0.0546346, 0.068735, 0 352 0.106025, 0.129614, 0.15643, 0.186063, 0.217 353 0.251065, 0.28525, 0.319171, 0.348995, 0.369 354 0.378165, 0.377353, 0.371224, 0.363257, 0.35 355 0.348618, 0.343082, 0.338825, 0.33574, 0.333 356 0.332504, 0.332052, 0.332187, 0.332781, 0.33 357 0.33489, 0.336213, 0.337608, 0.339008, 0.340 358 0.341606, 0.342706, 0.343628, 0.344305, 0.34 359 }; 360