Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // >> 27 // $Id: G4ElasticHadrNucleusHE.cc,v 1.79 2008/01/14 10:39:13 vnivanch Exp $ >> 28 // GEANT4 tag $Name: geant4-09-02 $ >> 29 // >> 30 // 26 // The generator of high energy hadron-nucleu 31 // The generator of high energy hadron-nucleus elastic scattering 27 // The hadron kinetic energy T > 1 GeV 32 // The hadron kinetic energy T > 1 GeV 28 // N.Starkov 2003. << 33 // N. Starkov 2003. 29 // 34 // >> 35 // 19.05.04 Variant for G4 6.1: The 'ApplyYourself' was changed 30 // 19.11.05 The HE elastic scattering on prot 36 // 19.11.05 The HE elastic scattering on proton is added (N.Starkov) 31 // 16.11.06 The low energy boundary is shifte 37 // 16.11.06 The low energy boundary is shifted to T = 400 MeV (N.Starkov) 32 // 23.11.06 General cleanup, ONQ0=3, use poin 38 // 23.11.06 General cleanup, ONQ0=3, use pointer instead of particle name (VI) 33 // 02.05.07 Scale sampled t as p^2 (VI) 39 // 02.05.07 Scale sampled t as p^2 (VI) 34 // 15.05.07 Redesign and cleanup (V.Ivanchenk 40 // 15.05.07 Redesign and cleanup (V.Ivanchenko) 35 // 17.05.07 cleanup (V.Grichine) 41 // 17.05.07 cleanup (V.Grichine) 36 // 19.04.12 Fixed reproducibility violation ( << 37 // 12.06.12 Fixed warnings of shadowed variab << 38 // 42 // 39 43 40 #include "G4ElasticHadrNucleusHE.hh" 44 #include "G4ElasticHadrNucleusHE.hh" 41 #include "G4PhysicalConstants.hh" << 42 #include "G4SystemOfUnits.hh" << 43 #include "Randomize.hh" 45 #include "Randomize.hh" 44 #include "G4ios.hh" 46 #include "G4ios.hh" 45 #include "G4ParticleTable.hh" 47 #include "G4ParticleTable.hh" 46 #include "G4NucleiProperties.hh" << 47 #include "G4IonTable.hh" 48 #include "G4IonTable.hh" 48 #include "G4Proton.hh" 49 #include "G4Proton.hh" 49 #include "G4PionPlus.hh" << 50 #include "G4PionMinus.hh" << 51 #include "G4NistManager.hh" 50 #include "G4NistManager.hh" 52 #include "G4ProductionCutsTable.hh" << 53 #include "G4MaterialCutsCouple.hh" << 54 #include "G4Material.hh" << 55 #include "G4Element.hh" << 56 #include "G4Log.hh" << 57 #include "G4Exp.hh" << 58 51 59 using namespace std; 52 using namespace std; 60 53 61 const G4int G4ElasticHadrNucleusHE::fHadronCod << 62 {211,-211,2112,2212,321,-321,130,310,311,-311, << 63 3122,3222,3112,3212,3312,3322,3334, << 64 -2212,-2112,-3122,-3222,-3112,-3212,-3312,-33 << 65 << 66 const G4int G4ElasticHadrNucleusHE::fHadronTyp << 67 {2,3,6,0,4,5,4,4,4,5, << 68 0,0,0,0,0,0,0, << 69 1,7,1,1,1,1,1,1,1}; << 70 << 71 const G4int G4ElasticHadrNucleusHE::fHadronTyp << 72 {3,4,1,0,5,6,5,5,5,6, << 73 0,0,0,0,0,0,0, << 74 2,2,2,2,2,2,2,2,2}; << 75 << 76 G4double G4ElasticHadrNucleusHE::fLineF[] = { << 77 G4double G4ElasticHadrNucleusHE::fEnergy[] = { << 78 G4double G4ElasticHadrNucleusHE::fLowEdgeEnerg << 79 G4double G4ElasticHadrNucleusHE::fBinom[240][2 << 80 << 81 G4ElasticData* << 82 G4ElasticHadrNucleusHE::fElasticData[NHADRONS] << 83 << 84 #ifdef G4MULTITHREADED << 85 G4Mutex G4ElasticHadrNucleusHE::elasticMutex << 86 #endif << 87 << 88 G4bool G4ElasticHadrNucleusHE::fStoreToFile = << 89 G4bool G4ElasticHadrNucleusHE::fRetrieveFromFi << 90 << 91 const G4double invGeV = 1.0/CLHEP::GeV; << 92 const G4double MbToGeV2 = 2.568; << 93 const G4double GeV2 = CLHEP::GeV*CLHEP:: << 94 const G4double invGeV2 = 1.0/GeV2; << 95 const G4double protonM = CLHEP::proton_mass << 96 const G4double protonM2 = protonM*protonM; << 97 << 98 ////////////////////////////////////////////// 54 /////////////////////////////////////////////////////////////// >> 55 // >> 56 // 99 57 100 G4ElasticData::G4ElasticData(const G4ParticleD 58 G4ElasticData::G4ElasticData(const G4ParticleDefinition* p, 101 G4int Z, G4int A, const G4double* e << 59 G4int Z, G4double AWeight, G4double* eGeV) 102 { 60 { 103 G4double massGeV = p->GetPDGMass()*invGeV; << 61 hadr = p; 104 G4double mass2GeV2= massGeV*massGeV; << 62 massGeV = p->GetPDGMass()/GeV; >> 63 mass2GeV2= massGeV*massGeV; >> 64 AtomicWeight = G4int(AWeight + 0.5); 105 65 106 DefineNucleusParameters(A); << 66 DefineNucleusParameters(AWeight); 107 G4double limitQ2 = 35./(R1*R1); // (GeV << 67 >> 68 limitQ2 = 35./(R1*R1); // (GeV/c)^2 >> 69 >> 70 G4double dQ2 = limitQ2/(ONQ2 - 1.); >> 71 >> 72 TableQ2[0] = 0.0; >> 73 >> 74 for(G4int ii = 1; ii < ONQ2; ii++) >> 75 { >> 76 TableQ2[ii] = TableQ2[ii-1]+dQ2; >> 77 } 108 78 109 massA = G4NucleiProperties::GetNuclearMass( << 79 massA = AWeight*amu_c2/GeV; 110 massA2 = massA*massA; 80 massA2 = massA*massA; 111 /* << 81 112 G4cout << " G4ElasticData for " << p->GetPar << 82 for(G4int kk = 0; kk < NENERGY; kk++) 113 << " Z= " << Z << " A= " << A << " R1= " << << 114 << " R2= " << R2 << G4endl; << 115 */ << 116 for(G4int kk = 0; kk<NENERGY; ++kk) << 117 { 83 { 118 G4double elab = e[kk] + massGeV; << 84 dnkE[kk] = 0; 119 G4double plab2= e[kk]*(e[kk] + 2.0*massGeV << 85 G4double elab = eGeV[kk] + massGeV; 120 G4double Q2m = 4.0*plab2*massA2/(mass2GeV << 86 G4double plab2= eGeV[kk]*(eGeV[kk] + 2.0*massGeV); >> 87 G4double Q2m = 4.0*plab2*massA2/(mass2GeV2 + massA2 + 2.*massA2*elab); 121 88 122 if(Z == 1 && p == G4Proton::Proton()) { Q2 << 89 if(Z == 1 && p == G4Proton::Proton()) Q2m *= 0.5; 123 90 124 maxQ2[kk] = Q2m; << 91 maxQ2[kk] = std::min(limitQ2, Q2m); 125 /* << 92 TableCrossSec[ONQ2*kk] = 0.0; 126 G4cout << " Ekin= " << e[kk] << " Q2m= " < << 127 << " limitQ2= " << limitQ2 << G4endl; << 128 */ << 129 } << 130 93 131 dQ2 = limitQ2/(G4double)(ONQ2-2); << 94 // G4cout<<" kk eGeV[kk] "<<kk<<" "<<eGeV[kk]<<G4endl; >> 95 } 132 } 96 } 133 97 134 ////////////////////////////////////////////// 98 ///////////////////////////////////////////////////////////////////////// >> 99 // >> 100 // 135 101 136 void G4ElasticData::DefineNucleusParameters(G4 << 102 void G4ElasticData::DefineNucleusParameters(G4double A) 137 { 103 { 138 switch (A) { << 104 switch (AtomicWeight) >> 105 { 139 case 207: 106 case 207: 140 case 208: 107 case 208: 141 R1 = 20.5; 108 R1 = 20.5; >> 109 // R1 = 17.5; >> 110 // R1 = 21.3; 142 R2 = 15.74; 111 R2 = 15.74; >> 112 // R2 = 10.74; >> 113 143 Pnucl = 0.4; 114 Pnucl = 0.4; 144 Aeff = 0.7; 115 Aeff = 0.7; 145 break; 116 break; 146 case 237: 117 case 237: 147 case 238: 118 case 238: 148 R1 = 21.7; 119 R1 = 21.7; 149 R2 = 16.5; 120 R2 = 16.5; 150 Pnucl = 0.4; 121 Pnucl = 0.4; 151 Aeff = 0.7; 122 Aeff = 0.7; 152 break; 123 break; 153 case 90: 124 case 90: 154 case 91: 125 case 91: 155 R1 = 16.5; << 126 R1 = 16.5*1.0; 156 R2 = 11.62; 127 R2 = 11.62; 157 Pnucl = 0.4; 128 Pnucl = 0.4; 158 Aeff = 0.7; 129 Aeff = 0.7; 159 break; 130 break; 160 case 58: 131 case 58: 161 case 59: 132 case 59: 162 R1 = 15.75; << 133 R1 = 15.0*1.05; 163 R2 = 9.9; 134 R2 = 9.9; 164 Pnucl = 0.45; 135 Pnucl = 0.45; 165 Aeff = 0.85; 136 Aeff = 0.85; 166 break; 137 break; 167 case 48: 138 case 48: 168 case 47: 139 case 47: 169 R1 = 14.0; 140 R1 = 14.0; 170 R2 = 9.26; 141 R2 = 9.26; 171 Pnucl = 0.31; 142 Pnucl = 0.31; 172 Aeff = 0.75; 143 Aeff = 0.75; 173 break; 144 break; 174 case 40: 145 case 40: 175 case 41: 146 case 41: 176 R1 = 13.3; 147 R1 = 13.3; 177 R2 = 9.26; 148 R2 = 9.26; 178 Pnucl = 0.31; 149 Pnucl = 0.31; 179 Aeff = 0.75; 150 Aeff = 0.75; 180 break; 151 break; 181 case 28: 152 case 28: 182 case 29: 153 case 29: 183 R1 = 12.0; 154 R1 = 12.0; 184 R2 = 7.64; 155 R2 = 7.64; 185 Pnucl = 0.253; 156 Pnucl = 0.253; 186 Aeff = 0.8; 157 Aeff = 0.8; 187 break; 158 break; 188 case 16: 159 case 16: 189 R1 = 10.50; 160 R1 = 10.50; 190 R2 = 5.5; 161 R2 = 5.5; 191 Pnucl = 0.7; 162 Pnucl = 0.7; 192 Aeff = 0.98; 163 Aeff = 0.98; 193 break; 164 break; 194 case 12: 165 case 12: 195 R1 = 9.3936; 166 R1 = 9.3936; 196 R2 = 4.63; 167 R2 = 4.63; 197 Pnucl = 0.7; 168 Pnucl = 0.7; >> 169 // Pnucl = 0.5397; 198 Aeff = 1.0; 170 Aeff = 1.0; 199 break; 171 break; 200 case 11: 172 case 11: 201 R1 = 9.0; 173 R1 = 9.0; 202 R2 = 5.42; 174 R2 = 5.42; 203 Pnucl = 0.19; 175 Pnucl = 0.19; >> 176 // Pnucl = 0.5397; 204 Aeff = 0.9; 177 Aeff = 0.9; 205 break; 178 break; 206 case 9: 179 case 9: 207 R1 = 9.9; 180 R1 = 9.9; 208 R2 = 6.5; 181 R2 = 6.5; 209 Pnucl = 0.690; 182 Pnucl = 0.690; 210 Aeff = 0.95; 183 Aeff = 0.95; 211 break; 184 break; 212 case 4: 185 case 4: 213 R1 = 5.3; 186 R1 = 5.3; 214 R2 = 3.7; 187 R2 = 3.7; 215 Pnucl = 0.4; 188 Pnucl = 0.4; 216 Aeff = 0.75; 189 Aeff = 0.75; 217 break; 190 break; 218 case 1: 191 case 1: 219 R1 = 4.5; 192 R1 = 4.5; 220 R2 = 2.3; 193 R2 = 2.3; 221 Pnucl = 0.177; 194 Pnucl = 0.177; 222 Aeff = 0.9; 195 Aeff = 0.9; 223 break; 196 break; 224 default: 197 default: 225 R1 = 4.45*G4Exp(G4Log((G4double)(A - << 198 R1 = 4.45*std::pow(A - 1.,0.309)*0.9; 226 R2 = 2.3 *G4Exp(G4Log((G4double)A)* 0 << 199 R2 = 2.3 *std::pow(A, 0.36); >> 200 >> 201 if(A < 100 && A > 3) Pnucl = 0.176 + 0.00275*A; >> 202 else Pnucl = 0.4; >> 203 >> 204 //G4cout<<" Deault: A= "<<A<<" R1 R2 Aeff Pnucl "<<R1<<" "<<R2<<" " >> 205 // <<Aeff<<" "<<Pnucl<<G4endl; 227 206 228 if(A < 100 && A > 3) { Pnucl = 0.176 + 0 << 207 if(A >= 100) Aeff = 0.7; 229 else { Pnucl = 0.4; } << 208 else if(A < 100 && A > 75) Aeff = 1.5 - 0.008*A; 230 //G4cout<<" Deault: A= "<<A<<" R1 R2 Ae << 209 else Aeff = 0.9; 231 // <<Aeff<<" "<<Pnucl<<G4endl; << 232 << 233 if(A >= 100) { Aeff = 0.7; << 234 else if(A < 100 && A > 75) { Aeff = 1.5 << 235 else { Aeff = 0.9; << 236 break; 210 break; 237 } << 211 } 238 //G4cout<<" Result: A= "<<A<<" R1 R2 Aeff P << 212 //G4cout<<" Result: A= "<<A<<" R1 R2 Aeff Pnucl "<<R1<<" "<<R2<<" " 239 // <<Aeff<<" "<<Pnucl<<G4endl; << 213 // <<Aeff<<" "<<Pnucl<<G4endl; 240 } 214 } 241 215 242 ////////////////////////////////////////////// 216 //////////////////////////////////////////////////////////////////// >> 217 // >> 218 // The constructor for the generating of events >> 219 // 243 220 244 G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE << 221 G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE() 245 : G4HadronElastic(name), fDirectory(nullptr) << 222 :G4HadronicInteraction("G4ElasticHadrNucleusHE") 246 { 223 { 247 dQ2 = hMass = hMass2 = hLabMomentum = hLabMo << 248 = R1 = R2 = Pnucl = Aeff = HadrTot = HadrS << 249 = DDSect3 = ConstU = Slope1 = Slope2 = Coe << 250 = Slope0 = Coeff0 = aAIm = aDIm = Dtot11 = << 251 iHadrCode = iHadron = iHadron1 = 0; << 252 << 253 verboseLevel = 0; 224 verboseLevel = 0; 254 ekinLowLimit = 400.0*CLHEP::MeV; << 225 plabLowLimit = 20.0*MeV; 255 << 226 lowestEnergyLimit = 0.0; 256 BoundaryP[0]=9.0; BoundaryTG[0]=5.0;Boundary << 257 BoundaryP[1]=20.0;BoundaryTG[1]=1.5;Boundary << 258 BoundaryP[2]=5.0; BoundaryTG[2]=1.0;Boundary << 259 BoundaryP[3]=8.0; BoundaryTG[3]=3.0;Boundary << 260 BoundaryP[4]=7.0; BoundaryTG[4]=3.0;Boundary << 261 BoundaryP[5]=5.0; BoundaryTG[5]=2.0;Boundary << 262 BoundaryP[6]=5.0; BoundaryTG[6]=1.5;Boundary << 263 227 >> 228 MbToGeV2 = 2.568; >> 229 sqMbToGeV = 1.602; >> 230 Fm2ToGeV2 = 25.68; >> 231 GeV2 = GeV*GeV; >> 232 protonM = proton_mass_c2/GeV; >> 233 protonM2 = protonM*protonM; >> 234 >> 235 BoundaryP[0]=9.0;BoundaryTG[0]=5.0;BoundaryTL[0]=0.; >> 236 BoundaryP[1]=20.0;BoundaryTG[1]=1.5;BoundaryTL[1]=0.; >> 237 BoundaryP[2]=5.0; BoundaryTG[2]=1.0;BoundaryTL[2]=1.5; >> 238 BoundaryP[3]=8.0; BoundaryTG[3]=3.0;BoundaryTL[3]=0.; >> 239 BoundaryP[4]=7.0; BoundaryTG[4]=3.0;BoundaryTL[4]=0.; >> 240 BoundaryP[5]=5.0; BoundaryTG[5]=2.0;BoundaryTL[5]=0.; >> 241 BoundaryP[6]=5.0; BoundaryTG[6]=1.5;BoundaryTL[6]=3.0; >> 242 >> 243 Binom(); >> 244 // energy in GeV >> 245 Energy[0] = 0.4; >> 246 Energy[1] = 0.6; >> 247 Energy[2] = 0.8; >> 248 LowEdgeEnergy[0] = 0.0; >> 249 LowEdgeEnergy[1] = 0.5; >> 250 LowEdgeEnergy[2] = 0.7; >> 251 G4double e = 1.0; >> 252 G4double f = std::pow(10.,0.1); >> 253 for(G4int i=3; i<NENERGY; i++) { >> 254 Energy[i] = e; >> 255 LowEdgeEnergy[i] = e/f; >> 256 e *= f*f; >> 257 } 264 nistManager = G4NistManager::Instance(); 258 nistManager = G4NistManager::Instance(); 265 259 266 if(fEnergy[0] == 0.0) { << 260 // PDG code for hadrons 267 #ifdef G4MULTITHREADED << 261 G4int ic[NHADRONS] = {211,-211,2112,2212,321,-321,130,310,311,-311, 268 G4MUTEXLOCK(&elasticMutex); << 262 3122,3222,3112,3212,3312,3322,3334, 269 if(fEnergy[0] == 0.0) { << 263 -2212,-2112,-3122,-3222,-3112,-3212,-3312,-3322,-3334}; 270 #endif << 264 // internal index 271 isMaster = true; << 265 G4int id[NHADRONS] = {2,3,6,0,4,5,4,4,4,5, 272 Binom(); << 266 0,0,0,0,0,0,0, 273 // energy in GeV << 267 1,7,1,1,1,1,1,1,1}; 274 fEnergy[0] = 0.4; << 268 275 fEnergy[1] = 0.6; << 269 G4int id1[NHADRONS] = {3,4,1,0,5,6,5,5,5,6, 276 fEnergy[2] = 0.8; << 270 0,0,0,0,0,0,0, 277 fEnergy[3] = 1.0; << 271 2,2,2,2,2,2,2,2,2}; 278 fLowEdgeEnergy[0] = 0.0; << 279 fLowEdgeEnergy[1] = 0.5; << 280 fLowEdgeEnergy[2] = 0.7; << 281 fLowEdgeEnergy[3] = 0.9; << 282 G4double f = G4Exp(G4Log(10.)*0.1); << 283 G4double e = f*f; << 284 for(G4int i=4; i<NENERGY; ++i) { << 285 fEnergy[i] = e; << 286 fLowEdgeEnergy[i] = e/f; << 287 e *= f*f; << 288 } << 289 if(verboseLevel > 0) { << 290 G4cout << "### G4ElasticHadrNucleusHE: energ << 291 for(G4int i=0; i<NENERGY; ++i) { << 292 G4cout << " " << i << " " << fLowEdgeEn << 293 << " " << fEnergy[i] << G4endl; << 294 } << 295 } << 296 #ifdef G4MULTITHREADED << 297 } << 298 G4MUTEXUNLOCK(&elasticMutex); << 299 #endif << 300 } << 301 } << 302 272 303 ////////////////////////////////////////////// << 273 for(G4int j=0; j<NHADRONS; j++) >> 274 { >> 275 HadronCode[j] = ic[j]; >> 276 HadronType[j] = id[j]; >> 277 HadronType1[j] = id1[j]; 304 278 305 void G4ElasticHadrNucleusHE::ModelDescription( << 279 for(G4int k = 0; k < 93; k++) SetOfElasticData[j][k] = 0; 306 { << 280 } 307 outFile << "G4ElasticHadrNucleusHE is a hadr << 308 << "model developed by N. Starkov which us << 309 << "parameterization to calculate the fina << 310 << "for all hadrons with incident momentum << 311 } 281 } 312 282 313 ////////////////////////////////////////////// 283 /////////////////////////////////////////////////////////////////// >> 284 // >> 285 // 314 286 315 G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusH 287 G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusHE() 316 { 288 { 317 if(isMaster) { << 289 for(G4int j = 0; j < NHADRONS; j++) 318 for(G4int j = 0; j < NHADRONS; ++j) { << 290 { 319 for(G4int k = 0; k < ZMAX; ++k) { << 291 for(G4int k = 0; k < 93; k++) 320 G4ElasticData* ptr = fElasticData[j][k]; << 292 { 321 if(ptr) { << 293 if(SetOfElasticData[j][k]) delete SetOfElasticData[j][k]; 322 delete ptr; << 323 fElasticData[j][k] = nullptr; << 324 for(G4int l = j+1; l < NHADRONS; ++l) { << 325 if(ptr == fElasticData[l][k]) { fElastic << 326 } << 327 } << 328 } << 329 } 294 } 330 delete fDirectory; << 295 } 331 fDirectory = nullptr; << 332 } << 333 } 296 } 334 297 335 ////////////////////////////////////////////// << 298 //////////////////////////////////////////////////////////////////// >> 299 // >> 300 // 336 301 337 void G4ElasticHadrNucleusHE::InitialiseModel() << 302 G4HadFinalState * G4ElasticHadrNucleusHE::ApplyYourself( >> 303 const G4HadProjectile &aTrack, >> 304 G4Nucleus &targetNucleus) 338 { 305 { 339 if(!isMaster) { return; } << 306 theParticleChange.Clear(); 340 G4ProductionCutsTable* theCoupleTable= << 307 341 G4ProductionCutsTable::GetProductionCutsTa << 308 const G4HadProjectile* aParticle = &aTrack; 342 G4int numOfCouples = (G4int)theCoupleTable-> << 309 G4double ekin = aParticle->GetKineticEnergy(); >> 310 >> 311 if( ekin <= lowestEnergyLimit ) >> 312 { >> 313 theParticleChange.SetEnergyChange(ekin); >> 314 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); >> 315 return &theParticleChange; >> 316 } >> 317 >> 318 G4double aTarget = targetNucleus.GetN(); >> 319 G4double zTarget = targetNucleus.GetZ(); >> 320 >> 321 G4double plab = aParticle->GetTotalMomentum(); >> 322 >> 323 if (verboseLevel >1) >> 324 { >> 325 G4cout << "G4ElasticHadrNucleusHE: Incident particle plab=" >> 326 << plab/GeV << " GeV/c " >> 327 << " ekin(MeV) = " << ekin/MeV << " " >> 328 << aParticle->GetDefinition()->GetParticleName() << G4endl; >> 329 } >> 330 // Scattered particle referred to axis of incident particle >> 331 >> 332 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); >> 333 G4double m1 = theParticle->GetPDGMass(); >> 334 >> 335 G4int Z = static_cast<G4int>(zTarget+0.5); >> 336 G4int A = static_cast<G4int>(aTarget+0.5); >> 337 G4int projPDG = theParticle->GetPDGEncoding(); >> 338 >> 339 if (verboseLevel>1) >> 340 { >> 341 G4cout << "G4ElasticHadrNucleusHE for " << theParticle->GetParticleName() >> 342 << " PDGcode= " << projPDG << " on nucleus Z= " << Z >> 343 << " A= " << A >> 344 << G4endl; >> 345 } >> 346 G4ParticleDefinition * theDef = 0; >> 347 >> 348 if (Z == 1 && A == 1) theDef = G4Proton::Proton(); >> 349 else if (Z == 1 && A == 2) theDef = G4Deuteron::Deuteron(); >> 350 else if (Z == 1 && A == 3) theDef = G4Triton::Triton(); >> 351 else if (Z == 2 && A == 3) theDef = G4He3::He3(); >> 352 else if (Z == 2 && A == 4) theDef = G4Alpha::Alpha(); >> 353 else theDef = G4ParticleTable:: >> 354 GetParticleTable()->FindIon(Z,A,0,Z); >> 355 >> 356 G4double m2 = theDef->GetPDGMass(); >> 357 G4LorentzVector lv1 = aParticle->Get4Momentum(); >> 358 G4LorentzVector lv(0.0,0.0,0.0,m2); >> 359 lv += lv1; >> 360 >> 361 G4ThreeVector bst = lv.boostVector(); >> 362 lv1.boost(-bst); >> 363 >> 364 G4ThreeVector p1 = lv1.vect(); >> 365 G4double ptot = p1.mag(); >> 366 G4double tmax = 4.0*ptot*ptot; >> 367 G4double t = 0.0; >> 368 >> 369 // Choose generator >> 370 G4bool swave = false; >> 371 >> 372 // S-wave for very low energy >> 373 if(plab < plabLowLimit) swave = true; >> 374 >> 375 // normal sampling >> 376 if(!swave) { >> 377 t = SampleT(theParticle,plab,Z,A); >> 378 if(t > tmax) swave = true; >> 379 } >> 380 >> 381 if(swave) t = G4UniformRand()*tmax; >> 382 >> 383 // Sampling in CM system >> 384 G4double phi = G4UniformRand()*twopi; >> 385 G4double cost = 1. - 2.0*t/tmax; >> 386 G4double sint; >> 387 >> 388 if( cost >= 1.0 ) >> 389 { >> 390 cost = 1.0; >> 391 sint = 0.0; >> 392 } >> 393 else if( cost <= -1.0) >> 394 { >> 395 cost = -1.0; >> 396 sint = 0.0; >> 397 } >> 398 else >> 399 { >> 400 sint = std::sqrt((1.0-cost)*(1.0+cost)); >> 401 } >> 402 if (verboseLevel>1) >> 403 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl; >> 404 >> 405 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); >> 406 v1 *= ptot; >> 407 G4LorentzVector nlv1( v1.x(), v1.y(), v1.z(), std::sqrt(ptot*ptot + m1*m1)); >> 408 >> 409 nlv1.boost(bst); >> 410 >> 411 G4double eFinal = nlv1.e() - m1; >> 412 >> 413 if (verboseLevel > 1) >> 414 { >> 415 G4cout << "Scattered: " >> 416 << nlv1<<" m= " << m1 << " ekin(MeV)= " << eFinal >> 417 << " Proj: 4-mom " << lv1 >> 418 <<G4endl; >> 419 } >> 420 if( eFinal < 0.0 ) >> 421 { >> 422 G4cout << "G4ElasticHadrNucleusHE WARNING ekin= " << eFinal >> 423 << " after scattering of " >> 424 << aParticle->GetDefinition()->GetParticleName() >> 425 << " p(GeV/c)= " << plab >> 426 << " on " << theDef->GetParticleName() >> 427 << G4endl; >> 428 eFinal = 0.0; >> 429 nlv1.setE(m1); >> 430 } >> 431 >> 432 theParticleChange.SetMomentumChange( nlv1.vect().unit() ); >> 433 theParticleChange.SetEnergyChange(eFinal); 343 434 344 for(G4int i=0; i<2; ++i) { << 435 G4LorentzVector nlv0 = lv - nlv1; 345 const G4ParticleDefinition* p = G4PionPlus << 436 G4double erec = nlv0.e() - m2; 346 if(1 == i) { p = G4PionMinus::PionMinus(); << 437 347 iHadrCode = fHadronCode[i]; << 438 if (verboseLevel > 1) 348 iHadron = fHadronType[i]; << 439 { 349 iHadron1 = fHadronType1[i]; << 440 G4cout << "Recoil: " 350 hMass = p->GetPDGMass()*invGeV; << 441 << nlv0<<" m= " << m2 << " ekin(MeV)= " << erec 351 hMass2 = hMass*hMass; << 442 <<G4endl; 352 for(G4int j=0; j<numOfCouples; ++j) { << 443 } 353 auto mat = theCoupleTable->GetMaterialCu << 444 if(erec < 0.0) 354 auto elmVec = mat->GetElementVector(); << 445 { 355 std::size_t numOfElem = mat->GetNumberOf << 446 G4cout << "G4ElasticHadrNucleusHE WARNING Erecoil(MeV)= " << erec 356 for(std::size_t k=0; k<numOfElem; ++k) { << 447 << " after scattering of " 357 G4int Z = std::min((*elmVec)[k]->GetZa << 448 << aParticle->GetDefinition()->GetParticleName() 358 if(!fElasticData[i][Z]) { << 449 << " p(GeV/c)= " << plab 359 if(1 == i && Z > 1) { << 450 << " on " << theDef->GetParticleName() 360 fElasticData[1][Z] = fElasticData[ << 451 << G4endl; 361 } else { << 452 nlv0.setE(m2); 362 FillData(p, i, Z); << 363 } << 364 } << 365 } << 366 } << 367 } 453 } >> 454 G4DynamicParticle * aSec = new G4DynamicParticle(theDef, nlv0); >> 455 theParticleChange.AddSecondary(aSec); >> 456 >> 457 return &theParticleChange; 368 } 458 } 369 459 370 ////////////////////////////////////////////// << 460 ////////////////////////////////////////////////////////////////////////// >> 461 // >> 462 // 371 463 372 G4double << 464 G4double G4ElasticHadrNucleusHE:: 373 G4ElasticHadrNucleusHE::SampleInvariantT(const << 465 SampleT( const G4ParticleDefinition* p, 374 G4double inLabMom, << 466 G4double inLabMom, G4int Z, G4int N) 375 G4int iZ, G4int A) << 376 { 467 { 377 G4double mass = p->GetPDGMass(); << 468 G4double plab = inLabMom/GeV; // (GeV/c) 378 G4double kine = sqrt(inLabMom*inLabMom + mas << 469 G4double Q2 = 0; 379 if(kine <= ekinLowLimit) { << 470 380 return G4HadronElastic::SampleInvariantT(p << 381 } << 382 G4int Z = std::min(iZ,ZMAX-1); << 383 G4double Q2 = 0.0; << 384 iHadrCode = p->GetPDGEncoding(); 471 iHadrCode = p->GetPDGEncoding(); 385 472 386 // below computations in GeV/c << 473 NumbN = N; 387 hMass = mass*invGeV; << 388 hMass2 = hMass*hMass; << 389 G4double plab = inLabMom*invGeV; << 390 G4double tmax = pLocalTmax*invGeV2; << 391 474 392 if(verboseLevel > 1) { << 475 if(verboseLevel > 1) 393 G4cout<< "G4ElasticHadrNucleusHE::SampleT: << 476 { >> 477 G4cout<< " G4ElasticHadrNucleusHE::SampleT: " 394 << " for " << p->GetParticleName() 478 << " for " << p->GetParticleName() 395 << " at Z= " << Z << " A= " << A << 479 << " at Z= " << Z << " A= " << N 396 << " plab(GeV)= " << plab 480 << " plab(GeV)= " << plab 397 << " hadrCode= " << iHadrCode << 398 << G4endl; 481 << G4endl; 399 } 482 } 400 iHadron = -1; << 401 G4int idx; 483 G4int idx; 402 for(idx=0; idx<NHADRONS; ++idx) { << 484 403 if(iHadrCode == fHadronCode[idx]) { << 485 for( idx = 0 ; idx < NHADRONS; idx++) 404 iHadron = fHadronType[idx]; << 486 { 405 iHadron1 = fHadronType1[idx]; << 487 if(iHadrCode == HadronCode[idx]) break; 406 break; << 407 } << 408 } 488 } >> 489 409 // Hadron is not in the list 490 // Hadron is not in the list 410 if(0 > iHadron) { return 0.0; } << 411 491 412 if(Z==1) { << 492 if( idx >= NHADRONS ) return Q2; 413 Q2 = HadronProtonQ2(plab, tmax); << 414 493 415 if (verboseLevel>1) { << 494 iHadron = HadronType[idx]; 416 G4cout<<" Proton : Q2 "<<Q2<<G4endl; << 495 iHadrCode = HadronCode[idx]; 417 } << 496 418 } else { << 497 if(Z==1) 419 const G4ElasticData* ElD1 = fElasticData[i << 498 { >> 499 hMass = p->GetPDGMass()/GeV; >> 500 hMass2 = hMass*hMass; >> 501 >> 502 G4double T = sqrt(plab*plab+hMass2)-hMass; 420 503 421 // Construct elastic data << 504 if(T > 0.4) Q2 = HadronProtonQ2(p, plab); 422 if(!ElD1) { << 505 423 FillData(p, idx, Z); << 506 if (verboseLevel>1) 424 ElD1 = fElasticData[idx][Z]; << 507 G4cout<<" Proton : Q2 "<<Q2<<G4endl; 425 if(!ElD1) { return 0.0; } << 426 } 508 } >> 509 else >> 510 { >> 511 G4ElasticData* ElD1 = SetOfElasticData[idx][Z]; 427 512 428 // sample scattering << 513 // Construct elastic data 429 Q2 = HadronNucleusQ2_2(ElD1, plab, tmax); << 514 if(!ElD1) >> 515 { >> 516 G4double AWeight = nistManager->GetAtomicMassAmu(Z); >> 517 ElD1 = new G4ElasticData(p, Z, AWeight, Energy); >> 518 SetOfElasticData[idx][Z] = ElD1; >> 519 >> 520 if(verboseLevel > 1) >> 521 { >> 522 G4cout<< " G4ElasticHadrNucleusHE::SampleT: new record " << idx >> 523 << " for " << p->GetParticleName() << " Z= " << Z >> 524 << G4endl; >> 525 } >> 526 } >> 527 hMass = ElD1->massGeV; >> 528 hMass2 = ElD1->mass2GeV2; >> 529 G4double M = ElD1->massA; >> 530 G4double M2 = ElD1->massA2; >> 531 G4double plab2 = plab*plab; >> 532 G4double Q2max = 4.*plab2*M2/ >> 533 (hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2)); >> 534 >> 535 // sample scattering >> 536 G4double T = sqrt(plab2+hMass2)-hMass; 430 537 431 if(verboseLevel > 1) { << 538 if(T > 0.4) Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max); 432 G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< " << 539 433 << Q2/tmax <<G4endl; << 540 if(verboseLevel > 1) >> 541 G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< " t/tmax= " << Q2/Q2max <<G4endl; 434 } 542 } 435 } << 543 return Q2*GeV2; 436 return Q2*GeV2; << 437 } 544 } 438 545 439 ////////////////////////////////////////////// << 546 ////////////////////////////////////////////////////////////////////////// >> 547 // >> 548 // 440 549 441 void G4ElasticHadrNucleusHE::FillData(const G4 << 550 G4double G4ElasticHadrNucleusHE:: 442 G4int id << 551 HadronNucleusQ2_2(G4ElasticData* pElD, G4int Z, >> 552 G4double plab, G4double tmax) 443 { 553 { 444 #ifdef G4MULTITHREADED << 554 G4double LineFq2[ONQ2]; 445 G4MUTEXLOCK(&elasticMutex); << 446 if(!fElasticData[idx][Z]) { << 447 #endif << 448 G4int A = G4lrint(nistManager->GetAtomicMa << 449 G4ElasticData* pElD = new G4ElasticData(p, << 450 if(fRetrieveFromFile) { << 451 std::ostringstream ss; << 452 InFileName(ss, p, Z); << 453 std::ifstream infile(ss.str(), std::ios: << 454 for(G4int i=0; i<NENERGY; ++i) { << 455 if(ReadLine(infile, pElD->fCumProb[i])) { << 456 continue; << 457 } else { << 458 fRetrieveFromFile = false; << 459 break; << 460 } << 461 } << 462 infile.close(); << 463 } << 464 R1 = pElD->R1; << 465 R2 = pElD->R2; << 466 Aeff = pElD->Aeff; << 467 Pnucl = pElD->Pnucl; << 468 dQ2 = pElD->dQ2; << 469 if(verboseLevel > 0) { << 470 G4cout<<"### FillData for " << p->GetPar << 471 << " Z= " << Z << " idx= " << idx << " i << 472 <<" iHadron1= " << iHadron1 << " iHadrCo << 473 <<"\n R1= " << R1 << " R2= " << << 474 <<" Pnucl= " << Pnucl << G4endl; << 475 } << 476 555 477 if(!fRetrieveFromFile) { << 556 G4double Rand = G4UniformRand(); 478 for(G4int i=0; i<NENERGY; ++i) { << 479 G4double T = fEnergy[i]; << 480 hLabMomentum2 = T*(T + 2.*hMass); << 481 hLabMomentum = std::sqrt(hLabMomentum2); << 482 HadrEnergy = hMass + T; << 483 DefineHadronValues(Z); << 484 Q2max = pElD->maxQ2[i]; << 485 << 486 G4int length = FillFq2(A); << 487 (pElD->fCumProb[i]).reserve(length); << 488 G4double norm = 1.0/fLineF[length-1]; << 489 << 490 if(verboseLevel > 0) { << 491 G4cout << "### i= " << i << " Z= " << Z << << 492 << " length= " << length << " Q2max= " << << 493 } << 494 557 495 (pElD->fCumProb[i]).push_back(0.0); << 558 G4int iNumbQ2 = 0; 496 for(G4int ii=1; ii<length-1; ++ii) { << 559 G4double Q2 = 0.0; 497 (pElD->fCumProb[i]).push_back(fLineF[ii]*n << 498 if(verboseLevel > 2) { << 499 G4cout << " ii= " << ii << " val= " << 500 << (pElD->fCumProb[i])[ii] << G4endl; << 501 } << 502 } << 503 (pElD->fCumProb[i]).push_back(1.0); << 504 } << 505 } << 506 560 507 if(fStoreToFile) { << 561 G4double ptot2 = plab*plab; 508 std::ostringstream ss; << 562 G4double ekin = std::sqrt(hMass2 + ptot2) - hMass; 509 OutFileName(ss, p, Z); << 563 510 std::ofstream fileout(ss.str()); << 564 if(verboseLevel > 1) 511 for(G4int i=0; i<NENERGY; ++i) { << 565 G4cout<<"Q2_2: ekin plab "<<ekin<<" "<<plab<<" tmax "<<tmax<<G4endl; 512 WriteLine(fileout, pElD->fCumProb[i]); << 513 } << 514 fileout.close(); << 515 } << 516 << 517 if(verboseLevel > 0) { << 518 G4cout << " G4ElasticHadrNucleusHE::Fill << 519 << " for " << p->GetParticleName() << " << 520 << " A= " << A << G4endl; << 521 } << 522 fElasticData[idx][Z] = pElD; << 523 566 524 #ifdef G4MULTITHREADED << 567 // Find closest energy bin >> 568 G4int NumbOnE; >> 569 for( NumbOnE = 0; NumbOnE < NENERGY-1; NumbOnE++ ) >> 570 { >> 571 if( ekin <= LowEdgeEnergy[NumbOnE+1] ) break; 525 } 572 } 526 G4MUTEXUNLOCK(&elasticMutex); << 573 G4double* dNumbQ2 = pElD->TableQ2; 527 #endif << 528 } << 529 574 530 ////////////////////////////////////////////// << 575 G4int index = NumbOnE*ONQ2; 531 576 532 void G4ElasticHadrNucleusHE::InterpolateHN(G4i << 577 // Select kinematics for node energy 533 const G4double C0P[], const << 578 G4double T = Energy[NumbOnE]; 534 const G4double B0P[], const << 579 hLabMomentum2 = T*(T + 2.*hMass); 535 { << 580 G4double Q2max = pElD->maxQ2[NumbOnE]; 536 G4int i; << 581 G4int length = pElD->dnkE[NumbOnE]; 537 582 538 for(i=1; i<n; ++i) { if(hLabMomentum <= EnP[ << 583 // Build vector 539 if(i == n) { i = n - 1; } << 584 if(length == 0) >> 585 { >> 586 R1 = pElD->R1; >> 587 R2 = pElD->R2; >> 588 Aeff = pElD->Aeff; >> 589 Pnucl = pElD->Pnucl; >> 590 hLabMomentum = std::sqrt(hLabMomentum2); >> 591 >> 592 DefineHadronValues(Z); 540 593 541 Coeff0 = LineInterpol(EnP[i], EnP[i-1], C0P[ << 594 if(verboseLevel >0) 542 Coeff1 = LineInterpol(EnP[i], EnP[i-1], C1P[ << 595 { 543 Slope0 = LineInterpol(EnP[i], EnP[i-1], B0P[ << 596 G4cout<<"1 plab T "<<plab<<" "<<T<<" sigTot B ReIm " 544 Slope1 = LineInterpol(EnP[i], EnP[i-1], B1P[ << 597 <<HadrTot<<" "<<HadrSlope<<" "<<HadrReIm<<G4endl; >> 598 G4cout<<" R1 R2 Aeff p "<<R1<<" "<<R2<<" "<<Aeff<<" " >> 599 <<Pnucl<<G4endl; >> 600 } 545 601 546 // G4cout<<" InterpolHN: n i "<<n<<" "<<i< << 602 pElD->CrossSecMaxQ2[NumbOnE] = 1.0; 547 // <<hLabMomentum<<G4endl; << 548 } << 549 603 550 ////////////////////////////////////////////// << 604 if(verboseLevel > 1) >> 605 G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE >> 606 << " length= " << length >> 607 << " Q2max= " << Q2max >> 608 << " ekin= " << ekin <<G4endl; >> 609 >> 610 pElD->TableCrossSec[index] = 0; 551 611 552 G4double << 553 G4ElasticHadrNucleusHE::HadronNucleusQ2_2(cons << 554 G4do << 555 { << 556 G4double ekin = std::sqrt(hMass2 + plab*pla << 557 612 558 if(verboseLevel > 1) { << 613 dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0]; 559 G4cout<<"Q2_2: ekin(GeV)= " << ekin << " << 560 <<" tmax(GeV2)= " << tmax <<G4endl; << 561 } << 562 // Find closest energy bin << 563 G4int idx; << 564 for(idx=0; idx<NENERGY-1; ++idx) { << 565 if(ekin <= fLowEdgeEnergy[idx+1]) { break; << 566 } << 567 //G4cout << " idx= " << idx << G4endl; << 568 614 569 // Select kinematics for node energy << 615 GetHeavyFq2(NumbN, LineFq2); // %%%%%%%%%%%%%%%%%%%%%%%%% 570 R1 = pElD->R1; << 571 dQ2 = pElD->dQ2; << 572 Q2max = pElD->maxQ2[idx]; << 573 G4int length = (G4int)(pElD->fCumProb[idx]). << 574 616 575 G4double Rand = G4UniformRand(); << 617 for(G4int ii=0; ii<ONQ2; ii++) >> 618 { >> 619 //if(verboseLevel > 2) >> 620 // G4cout<<" ii LineFq2 "<<ii<<" "<<LineFq2[ii]/LineFq2[ONQ2-1] >> 621 // <<" dF(q2) "<<LineFq2[ii]-LineFq2[ii-1]<<G4endl; 576 622 577 G4int iNumbQ2 = 0; << 623 pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1]; 578 for(iNumbQ2=1; iNumbQ2<length; ++iNumbQ2) { << 624 } 579 if(Rand <= (pElD->fCumProb[idx])[iNumbQ2]) << 625 580 } << 626 pElD->dnkE[NumbOnE] = ONQ2; 581 iNumbQ2 = std::min(iNumbQ2, length - 1); << 627 length = ONQ2; 582 G4double Q2 = GetQ2_2(iNumbQ2, length, pElD- << 628 } 583 Q2 = std::min(Q2, Q2max); << 584 Q2 *= tmax/Q2max; << 585 629 586 if(verboseLevel > 1) { << 630 G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]); >> 631 >> 632 for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ ) >> 633 { >> 634 if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break; >> 635 } >> 636 Q2 = GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand); >> 637 >> 638 if(tmax < Q2max) Q2 *= tmax/Q2max; >> 639 >> 640 if(verboseLevel > 1) 587 G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<" 641 G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<" iNumbQ2= " << iNumbQ2 588 << " rand= " << Rand << " Q2max= " << Q2ma << 642 << " rand= " << Rand << G4endl; 589 << " tmax= " << tmax << G4endl; << 643 590 } << 591 return Q2; 644 return Q2; 592 } 645 } 593 646 594 ////////////////////////////////////////////// 647 /////////////////////////////////////////////////////////////////////// 595 // 648 // 596 // The randomization of one dimensional array 649 // The randomization of one dimensional array 597 // 650 // 598 << 651 G4double G4ElasticHadrNucleusHE::GetQ2_2(G4int kk, G4double * Q, 599 G4double G4ElasticHadrNucleusHE::GetQ2_2(G4int << 652 G4double * F, G4double ranUni) 600 const std::vector<G4double>& F, << 601 G4dou << 602 { 653 { 603 //G4cout << "GetQ2_2 kk= " << kk << " kmax= << 654 G4double ranQ2; 604 // << F.size() << " rand= " << ranUni << << 605 if(kk == kmax-1) { << 606 G4double X1 = dQ2*kk; << 607 G4double F1 = F[kk-1]; << 608 G4double X2 = Q2max; << 609 G4double xx = R1*(X2 - X1); << 610 xx = (xx > 20.) ? 0.0 : G4Exp(-xx); << 611 G4double Y = X1 - G4Log(1.0 - (ranUni - F1 << 612 return Y; << 613 } << 614 G4double F1, F2, F3, X1, X2, X3; << 615 655 616 if(kk == 1 || kk == 0) { << 656 G4double F1 = F[kk-2]; 617 F1 = F[0]; << 657 G4double F2 = F[kk-1]; 618 F2 = F[1]; << 658 G4double F3 = F[kk]; 619 F3 = F[2]; << 659 G4double X1 = Q[kk-2]; 620 X1 = 0.0; << 660 G4double X2 = Q[kk-1]; 621 X2 = dQ2; << 661 G4double X3 = Q[kk]; 622 X3 = dQ2*2; << 662 623 } else { << 663 if(verboseLevel > 2) 624 F1 = F[kk-2]; << 625 F2 = F[kk-1]; << 626 F3 = F[kk]; << 627 X1 = dQ2*(kk-2); << 628 X2 = dQ2*(kk-1); << 629 X3 = dQ2*kk; << 630 } << 631 if(verboseLevel > 1) { << 632 G4cout << "GetQ2_2 kk= " << kk << " X2= " 664 G4cout << "GetQ2_2 kk= " << kk << " X2= " << X2 << " X3= " << X3 633 << " F2= " << F2 << " F3= " << F3 << " Rn 665 << " F2= " << F2 << " F3= " << F3 << " Rndm= " << ranUni << G4endl; >> 666 >> 667 if(kk == 1 || kk == 0) >> 668 { >> 669 F1 = F[0]; >> 670 F2 = F[1]; >> 671 F3 = F[2]; >> 672 X1 = Q[0]; >> 673 X2 = Q[1]; >> 674 X3 = Q[2]; 634 } 675 } 635 676 636 G4double F12 = F1*F1; 677 G4double F12 = F1*F1; 637 G4double F22 = F2*F2; 678 G4double F22 = F2*F2; 638 G4double F32 = F3*F3; 679 G4double F32 = F3*F3; 639 680 640 G4double D0 = F12*F2+F1*F32+F3*F22-F32*F2-F 681 G4double D0 = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3; 641 682 642 if(verboseLevel > 2) { << 683 if(verboseLevel > 2) 643 G4cout << " X1= " << X1 << " F1= " < 684 G4cout << " X1= " << X1 << " F1= " << F1 << " D0= " 644 << D0 << G4endl; 685 << D0 << G4endl; 645 } << 686 646 G4double Y; << 687 if(std::abs(D0) < 0.00000001) 647 if(std::abs(D0) < 1.e-9) { << 688 { 648 Y = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2) << 689 ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2); 649 } else { << 690 } 650 G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F << 691 else 651 G4double DB = X2*F12+X1*F32+X3*F22-X2*F32- << 692 { 652 G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F2 << 693 G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1; >> 694 G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*F22; >> 695 G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22 653 -X1*F2*F32-X2*F3*F12-X3*F1*F22; 696 -X1*F2*F32-X2*F3*F12-X3*F1*F22; 654 Y = (DA*ranUni*ranUni + DB*ranUni + DC)/D0 << 697 ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0; 655 } << 698 } 656 return Y; << 699 return ranQ2; // MeV^2 657 } 700 } 658 701 659 ////////////////////////////////////////////// 702 //////////////////////////////////////////////////////////////////////// 660 << 703 // 661 G4int G4ElasticHadrNucleusHE::FillFq2(G4int A) << 704 // >> 705 G4double G4ElasticHadrNucleusHE::GetHeavyFq2(G4int Nucleus, G4double* LineF) 662 { 706 { >> 707 G4int ii, jj, aSimp; 663 G4double curQ2, curSec; 708 G4double curQ2, curSec; 664 G4double curSum = 0.0; 709 G4double curSum = 0.0; 665 G4double totSum = 0.0; 710 G4double totSum = 0.0; 666 711 667 G4double ddQ2 = dQ2*0.1; << 712 G4double ddQ2 = dQ2/20; 668 G4double Q2l = 0.0; << 713 G4double Q2l = 0; 669 714 670 G4int ii = 0; << 715 LineF[0] = 0; 671 for(ii=1; ii<ONQ2-1; ++ii) { << 716 for(ii = 1; ii<ONQ2; ii++) 672 curSum = curSec = 0.0; << 717 { 673 << 718 curSum = 0; 674 for(G4int jj=0; jj<10; ++jj) { << 719 aSimp = 4; 675 curQ2 = Q2l+(jj + 0.5)*ddQ2; << 720 676 if(curQ2 >= Q2max) { break; } << 721 for(jj = 0; jj<20; jj++) 677 curSec = HadrNucDifferCrSec(A, curQ2); << 722 { 678 curSum += curSec; << 723 curQ2 = Q2l+jj*ddQ2; 679 } << 724 680 G4double del = (curQ2 >= Q2max) ? Q2max - << 725 curSec = HadrNucDifferCrSec(Nucleus, curQ2); 681 Q2l += del; << 726 curSum += curSec*aSimp; 682 curSum *= del*0.1; << 727 683 totSum += curSum; << 728 if(aSimp > 3) aSimp = 2; 684 fLineF[ii] = totSum; << 729 else aSimp = 4; 685 if (verboseLevel>2) { << 730 686 G4cout<<ii << ". FillFq2: A= " << A << " << 731 if(jj == 0 && verboseLevel>2) 687 <<dQ2<<" Tot= "<<totSum << " dTot " <<cu << 732 G4cout<<" Q2 "<<curQ2<<" AIm "<<aAIm<<" DIm "<<aDIm 688 <<" curSec= " <<curSec<<G4endl; << 733 <<" Diff "<<curSec<<" totSum "<<totSum<<G4endl; 689 } << 734 } 690 if(totSum*1.e-4 > curSum || Q2l >= Q2max) << 735 691 } << 736 Q2l += dQ2; 692 ii = std::min(ii, ONQ2-2); << 737 curSum *= ddQ2/2.3; // $$$$$$$$$$$$$$$$$$$$$$$ 693 curQ2 = Q2l; << 738 totSum += curSum; 694 G4double xx = R1*(Q2max - curQ2); << 739 695 if(xx > 0.0) { << 740 LineF[ii] = totSum; 696 xx = (xx > 20.) ? 0.0 : G4Exp(-xx); << 741 697 curSec = HadrNucDifferCrSec(A, curQ2); << 742 if (verboseLevel>2) 698 totSum += curSec*(1.0 - xx)/R1; << 743 G4cout<<" GetHeavy: Q2 dQ2 totSum "<<Q2l<<" "<<dQ2<<" "<<totSum 699 } << 744 <<" curSec " 700 fLineF[ii + 1] = totSum; << 745 <<curSec<<" totSum "<< totSum<<" DTot " 701 if (verboseLevel>1) { << 746 <<curSum<<G4endl; 702 G4cout << "### FillFq2 done curQ2= " << cu << 747 } 703 << " sumG= " << fLineF[ONQ2-2] << " << 748 return totSum; 704 << " Nbins= " << ii + 1 << G4endl; << 705 } << 706 return ii + 2; << 707 } 749 } 708 750 709 ////////////////////////////////////////////// 751 //////////////////////////////////////////////////////////////////////// >> 752 // >> 753 // 710 754 711 G4double G4ElasticHadrNucleusHE::GetLightFq2(G << 755 G4double G4ElasticHadrNucleusHE::GetLightFq2(G4int Z, G4int Nucleus, >> 756 G4double Q2) 712 { 757 { 713 // Scattering off proton << 758 // Scattering of proton 714 if(Z == 1) 759 if(Z == 1) 715 { 760 { 716 G4double SqrQ2 = std::sqrt(Q2); 761 G4double SqrQ2 = std::sqrt(Q2); 717 G4double valueConstU = 2.*(hMass2 + proton << 762 G4double ConstU = 2.*(hMass2 + protonM2) - Q2; 718 763 719 G4double y = (1.-Coeff1-Coeff0)/HadrSlope* << 764 G4double y = (1.-Coeff1-Coeff0)/HadrSlope*(1.-std::exp(-HadrSlope*Q2)) 720 + Coeff0*(1.-G4Exp(-Slope0*Q2)) << 765 + Coeff0*(1.-std::exp(-Slope0*Q2)) 721 + Coeff2/Slope2*G4Exp(Slope2*valueConstU << 766 + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1.) 722 + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1 << 767 + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2)); 723 768 724 return y; 769 return y; 725 } 770 } 726 771 727 // The preparing of probability function 772 // The preparing of probability function 728 773 729 G4double prec = A > 208 ? 1.0e-7 : 1.0e-6; << 774 G4double prec = Nucleus > 208 ? 1.0e-7 : 1.0e-6; 730 775 731 G4double Stot = HadrTot*MbToGeV2; 776 G4double Stot = HadrTot*MbToGeV2; // Gev^-2 732 G4double Bhad = HadrSlope; // 777 G4double Bhad = HadrSlope; // GeV^-2 733 G4double Asq = 1+HadrReIm*HadrReIm; 778 G4double Asq = 1+HadrReIm*HadrReIm; 734 G4double Rho2 = std::sqrt(Asq); 779 G4double Rho2 = std::sqrt(Asq); 735 780 736 if(verboseLevel >1) { << 781 // if(verboseLevel >1) 737 G4cout<<" Fq2 Before for i Tot B Im "<<Had 782 G4cout<<" Fq2 Before for i Tot B Im "<<HadrTot<<" "<<HadrSlope<<" " 738 <<HadrReIm<<G4endl; 783 <<HadrReIm<<G4endl; 739 } << 784 740 if(verboseLevel > 1) { 785 if(verboseLevel > 1) { 741 G4cout << "GetFq2: Stot= " << Stot << " Bh 786 G4cout << "GetFq2: Stot= " << Stot << " Bhad= " << Bhad 742 <<" Im "<<HadrReIm 787 <<" Im "<<HadrReIm 743 << " Asq= " << Asq << G4endl; 788 << " Asq= " << Asq << G4endl; 744 G4cout << "R1= " << R1 << " R2= " << R2 << 789 G4cout << "R1= " << R1 << " R2= " << R2 << " Pnucl= " << Pnucl <<G4endl; 745 } 790 } 746 G4double R12 = R1*R1; 791 G4double R12 = R1*R1; 747 G4double R22 = R2*R2; 792 G4double R22 = R2*R2; 748 G4double R12B = R12+2*Bhad; 793 G4double R12B = R12+2*Bhad; 749 G4double R22B = R22+2*Bhad; 794 G4double R22B = R22+2*Bhad; 750 795 751 G4double Norm = (R12*R1-Pnucl*R22*R2) 796 G4double Norm = (R12*R1-Pnucl*R22*R2); // HP->Aeff; 752 797 753 G4double R13 = R12*R1/R12B; 798 G4double R13 = R12*R1/R12B; 754 G4double R23 = Pnucl*R22*R2/R22B; 799 G4double R23 = Pnucl*R22*R2/R22B; 755 G4double Unucl = Stot/twopi*R13/Norm; << 800 G4double Unucl = Stot/twopi/Norm*R13; 756 G4double UnucRho2 = -Unucl*Rho2; 801 G4double UnucRho2 = -Unucl*Rho2; 757 802 758 G4double FiH = std::asin(HadrReIm/Rh 803 G4double FiH = std::asin(HadrReIm/Rho2); 759 G4double NN2 = R23/R13; 804 G4double NN2 = R23/R13; 760 805 761 if(verboseLevel > 2) { << 806 if(verboseLevel > 2) 762 G4cout << "UnucRho2= " << UnucRho2 << " Fi 807 G4cout << "UnucRho2= " << UnucRho2 << " FiH= " << FiH << " NN2= " << NN2 763 << " Norm= " << Norm << G4endl; 808 << " Norm= " << Norm << G4endl; 764 } << 765 G4double Prod0 = 0.; << 766 G4double N1 = -1.0; << 767 809 768 for(G4int i1 = 1; i1<= A; ++i1) ////++++++++ << 810 G4double dddd; >> 811 >> 812 G4double Prod0 = 0; >> 813 G4double N1 = -1.0; >> 814 G4double Tot0 = 0; >> 815 G4double exp1; >> 816 >> 817 G4double Prod3 ; >> 818 G4double exp2 ; >> 819 G4double N4, N5, N2, Prod1, Prod2; >> 820 G4int i1, i2, m1, m2; >> 821 >> 822 for(i1 = 1; i1<= Nucleus; i1++) ////++++++++++ i1 769 { 823 { 770 N1 *= (-Unucl*Rho2*(A-i1+1)/(G4double)i1 << 824 N1 = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2; 771 G4double Prod1 = 0.; << 825 Prod1 = 0; 772 G4double N2 = -1.; << 826 Tot0 = 0; >> 827 N2 = -1; 773 828 774 for(G4int i2 = 1; i2<=A; ++i2) ////+++++ << 829 for(i2 = 1; i2<=Nucleus; i2++) ////+++++++++ i2 775 { 830 { 776 N2 *= (-Unucl*Rho2*(A-i2+1)/(G4doubl << 831 N2 = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2; 777 G4double Prod2 = 0; << 832 Prod2 = 0; 778 G4double N5 = -1/NN2; << 833 N5 = -1/NN2; 779 for(G4int j2=0; j2<= i2; ++j2) ////+++++++ << 834 for(m2=0; m2<= i2; m2++) ////+++++++++ m2 780 { 835 { 781 G4double Prod3 = 0; << 836 Prod3 = 0; 782 G4double exp2 = 1./((G4double)j << 837 exp2 = 1/(m2/R22B+(i2-m2)/R12B); 783 N5 *= (-NN2); << 838 N5 = -N5*NN2; 784 G4double N4 = -1./NN2; << 839 N4 = -1/NN2; 785 for(G4int j1=0; j1<=i1; ++j1) ////++++ << 840 for(m1=0; m1<=i1; m1++) ////++++++++ m1 786 { 841 { 787 G4double exp1 = 1./((G4double)j1/R22B+( << 842 exp1 = 1/(m1/R22B+(i1-m1)/R12B); 788 G4double dddd = 0.25*(exp1+exp2); << 843 dddd = exp1+exp2; 789 N4 *= (-NN2); << 844 N4 = -N4*NN2; 790 Prod3 += << 845 Prod3 = Prod3+N4*exp1*exp2* 791 N4*exp1*exp2*(1.-G4Exp(-Q2 << 846 (1-std::exp(-Q2*dddd/4))/dddd*4*SetBinom[i1][m1]; 792 } // j1 << 847 } // m1 793 Prod2 += Prod3*N5*GetBinomCof(i2,j2); << 848 Prod2 = Prod2 +Prod3*N5*SetBinom[i2][m2]; 794 } / << 849 } // m2 795 Prod1 += Prod2*N2*std::cos(FiH*(i1-i2)); << 850 Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2)); 796 851 797 if (std::abs(Prod2*N2/Prod1)<prec) break; << 852 if (std::fabs(Prod2*N2/Prod1)<prec) break; 798 } 853 } // i2 799 Prod0 += Prod1*N1; << 854 Prod0 = Prod0 + Prod1*N1; 800 if(std::abs(N1*Prod1/Prod0) < prec) brea << 855 if(std::fabs(N1*Prod1/Prod0) < prec) break; >> 856 801 } 857 } // i1 802 858 803 const G4double fact = 0.25*CLHEP::pi/MbToGeV << 859 Prod0 *= 0.25*pi/MbToGeV2; // This is in mb 804 Prod0 *= fact; // This is in mb << 805 860 806 if(verboseLevel>1) { << 861 if(verboseLevel>1) 807 G4cout << "GetLightFq2 Z= " << Z << " A= " << 862 G4cout << "GetLightFq2 Z= " << Z << " A= " << Nucleus 808 <<" Q2= " << Q2 << " Res= " << Prod0 << G 863 <<" Q2= " << Q2 << " Res= " << Prod0 << G4endl; 809 } << 810 return Prod0; 864 return Prod0; 811 } 865 } 812 << 866 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 813 ////////////////////////////////////////////// << 867 G4double G4ElasticHadrNucleusHE:: 814 << 868 HadrNucDifferCrSec(G4int Nucleus, G4double aQ2) 815 G4double << 816 G4ElasticHadrNucleusHE::HadrNucDifferCrSec(G4i << 817 { 869 { 818 // ------ All external kinematical variabl << 870 // ------ All external kinematical variables are in MeV ------- 819 // ------ but internal in GeV !!! << 871 // ------ but internal in GeV !!!! ------ >> 872 >> 873 G4double theQ2 = aQ2; ///GeV/GeV; 820 874 821 // Scattering of proton 875 // Scattering of proton 822 if(A == 1) << 876 if(Nucleus == 1) 823 { 877 { 824 G4double SqrQ2 = std::sqrt(aQ2); 878 G4double SqrQ2 = std::sqrt(aQ2); 825 G4double valueConstU = hMass2 + protonM2-2 << 879 G4double ConstU = hMass2 + protonM2-2*protonM*HadrEnergy - aQ2; 826 << 880 827 BoundaryTL[0] = Q2max; << 881 G4double MaxT = 4*MomentumCM*MomentumCM; 828 BoundaryTL[1] = Q2max; << 882 829 BoundaryTL[3] = Q2max; << 883 BoundaryTL[0] = MaxT; 830 BoundaryTL[4] = Q2max; << 884 BoundaryTL[1] = MaxT; 831 BoundaryTL[5] = Q2max; << 885 BoundaryTL[3] = MaxT; 832 << 886 BoundaryTL[4] = MaxT; 833 G4double dSigPodT = HadrTot*HadrTot*(1+Had << 887 BoundaryTL[5] = MaxT; 834 ( Coeff1*G4Exp(-Slope1*SqrQ2)+ << 888 835 Coeff2*G4Exp( Slope2*(valueConstU)+aQ2)+ << 889 G4double dSigPodT; 836 (1-Coeff1-Coeff0)*G4Exp(-HadrSlope*aQ2)+ << 890 837 Coeff0*G4Exp(-Slope0*aQ2) )*2.568/(16*pi); << 891 dSigPodT = HadrTot*HadrTot*(1+HadrReIm*HadrReIm)* >> 892 ( >> 893 Coeff1*std::exp(-Slope1*SqrQ2)+ >> 894 Coeff2*std::exp( Slope2*(ConstU)+aQ2)+ >> 895 (1-Coeff1-Coeff0)*std::exp(-HadrSlope*aQ2)+ >> 896 +Coeff0*std::exp(-Slope0*aQ2) >> 897 // +0.1*(1-std::fabs(CosTh)) >> 898 )/16/3.1416*2.568; 838 899 839 return dSigPodT; 900 return dSigPodT; 840 } 901 } 841 902 842 G4double Stot = HadrTot*MbToGeV2; << 903 G4double Stot = HadrTot*MbToGeV2; 843 G4double Bhad = HadrSlope; << 904 G4double Bhad = HadrSlope; 844 G4double Asq = 1+HadrReIm*HadrReIm; << 905 G4double Asq = 1+HadrReIm*HadrReIm; 845 G4double Rho2 = std::sqrt(Asq); << 906 G4double Rho2 = std::sqrt(Asq); 846 G4double R12 = R1*R1; << 907 G4double Pnuclp = 0.001; 847 G4double R22 = R2*R2; << 908 Pnuclp = Pnucl; 848 G4double R12B = R12+2*Bhad; << 909 G4double R12 = R1*R1; 849 G4double R22B = R22+2*Bhad; << 910 G4double R22 = R2*R2; 850 G4double R12Ap = R12+20; << 911 G4double R12B = R12+2*Bhad; 851 G4double R22Ap = R22+20; << 912 G4double R22B = R22+2*Bhad; 852 G4double R13Ap = R12*R1/R12Ap; << 913 G4double R12Ap = R12+20; 853 G4double R23Ap = R22*R2*Pnucl/R22Ap; << 914 G4double R22Ap = R22+20; 854 G4double R23dR13 = R23Ap/R13Ap; << 915 G4double R13Ap = R12*R1/R12Ap; 855 G4double R12Apd = 2/R12Ap; << 916 G4double R23Ap = R22*R2/R22Ap*Pnuclp; 856 G4double R22Apd = 2/R22Ap; << 917 G4double R23dR13 = R23Ap/R13Ap; 857 G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd); << 918 G4double R12Apd = 2/R12Ap; 858 << 919 G4double R22Apd = 2/R22Ap; 859 G4double DDSec1p = (DDSect2+DDSect3*G4Log(0 << 920 G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd); 860 G4double DDSec2p = (DDSect2+DDSect3*G4Log(0 << 921 861 std::sqrt((R12+R2 << 922 G4double DDSec1p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R1/4)); 862 G4double DDSec3p = (DDSect2+DDSect3*G4Log(0 << 923 G4double DDSec2p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/ >> 924 std::sqrt((R12+R22)/2)/4)); >> 925 G4double DDSec3p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R2/4)); >> 926 >> 927 G4double Norm = (R12*R1-Pnucl*R22*R2)*Aeff; >> 928 G4double Normp = (R12*R1-Pnuclp*R22*R2)*Aeff; >> 929 G4double R13 = R12*R1/R12B; >> 930 G4double R23 = Pnucl*R22*R2/R22B; >> 931 G4double Unucl = Stot/twopi/Norm*R13; >> 932 G4double UnuclScr = Stot/twopi/Normp*R13Ap; >> 933 G4double SinFi = HadrReIm/Rho2; >> 934 G4double FiH = std::asin(SinFi); >> 935 G4double N = -1; >> 936 G4double N2 = R23/R13; >> 937 >> 938 G4double ImElasticAmpl0 = 0; >> 939 G4double ReElasticAmpl0 = 0; >> 940 >> 941 G4double exp1; >> 942 G4double N4; >> 943 G4double Prod1, Tot1=0, medTot, DTot1, DmedTot; >> 944 G4int i; 863 945 864 G4double Norm = (R12*R1-Pnucl*R22*R2) << 946 for( i=1; i<=Nucleus; i++) 865 G4double R13 = R12*R1/R12B; << 947 { 866 G4double R23 = Pnucl*R22*R2/R22B; << 948 N = -N*Unucl*(Nucleus-i+1)/i*Rho2; 867 G4double Unucl = Stot/(twopi*Norm)*R13 << 949 N4 = 1; 868 G4double UnuclScr = Stot/(twopi*Norm)*R13 << 950 Prod1 = std::exp(-theQ2/i*R12B/4)/i*R12B; 869 G4double SinFi = HadrReIm/Rho2; << 951 medTot = R12B/i; 870 G4double FiH = std::asin(SinFi); << 952 871 G4double N = -1; << 953 for(G4int l=1; l<=i; l++) 872 G4double N2 = R23/R13; << 954 { 873 << 955 exp1 = l/R22B+(i-l)/R12B; 874 G4double ImElasticAmpl0 = 0; << 956 N4 = -N4*(i-l+1)/l*N2; 875 G4double ReElasticAmpl0 = 0; << 957 Prod1 = Prod1+N4/exp1*std::exp(-theQ2/exp1/4); 876 G4double exp1; << 958 medTot = medTot+N4/exp1; 877 << 959 } // end l 878 for(G4int i=1; i<=A; ++i) { << 960 879 N *= (-Unucl*Rho2*(A-i+1)/(G4double)i); << 961 ReElasticAmpl0 = ReElasticAmpl0+Prod1*N*std::sin(FiH*i); 880 G4double N4 = 1; << 962 ImElasticAmpl0 = ImElasticAmpl0+Prod1*N*std::cos(FiH*i); 881 G4double medTot = R12B/(G4double)i; << 963 Tot1 = Tot1+medTot*N*std::cos(FiH*i); 882 G4double Prod1 = G4Exp(-aQ2*R12B/(G4doubl << 964 if(std::fabs(Prod1*N/ImElasticAmpl0) < 0.000001) break; 883 << 965 } // i 884 for(G4int l=1; l<=i; ++l) { << 966 885 exp1 = l/R22B+(i-l)/R12B; << 967 ImElasticAmpl0 = ImElasticAmpl0*pi/2.568; // The amplitude in mB 886 N4 *= (-N2*(i-l+1)/(G4double)l); << 968 ReElasticAmpl0 = ReElasticAmpl0*pi/2.568; // The amplitude in mB 887 G4double expn4 = N4/exp1; << 969 Tot1 = Tot1*twopi/2.568; 888 Prod1 += expn4*G4Exp(-aQ2/(exp1*4)); << 970 889 medTot += expn4; << 971 G4double C1 = R13Ap*R13Ap/2*DDSec1p; 890 } // end l << 972 G4double C2 = 2*R23Ap*R13Ap/2*DDSec2p; 891 << 973 G4double C3 = R23Ap*R23Ap/2*DDSec3p; 892 G4double dcos = N*std::cos(FiH*i); << 974 893 ReElasticAmpl0 += Prod1*N*std::sin(FiH*i) << 975 G4double N1p = 1; 894 ImElasticAmpl0 += Prod1*dcos; << 976 895 if(std::abs(Prod1*N/ImElasticAmpl0) < 0.00 << 977 G4double Din1 = 0.5; 896 } // i << 978 897 << 979 Din1 = 0.5*(C1*std::exp(-theQ2/8*R12Ap)/2*R12Ap- 898 static const G4double pi25 = CLHEP::pi/2.568 << 980 C2/R12ApdR22Ap*std::exp(-theQ2/4/R12ApdR22Ap)+ 899 ImElasticAmpl0 *= pi25; // The amplitude i << 981 C3*R22Ap/2*std::exp(-theQ2/8*R22Ap)); 900 ReElasticAmpl0 *= pi25; // The amplitude i << 982 901 << 983 DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2); 902 G4double C1 = R13Ap*R13Ap*0.5*DDSec1p; << 984 903 G4double C2 = 2*R23Ap*R13Ap*0.5*DDSec2p; << 985 G4double exp1p; 904 G4double C3 = R23Ap*R23Ap*0.5*DDSec3p; << 986 G4double exp2p; 905 << 987 G4double exp3p; 906 G4double N1p = 1; << 988 G4double N2p; 907 G4double Din1 = 0.5*(C1*G4Exp(-aQ2/8*R12Ap)/ << 989 G4double Din2, BinCoeff; 908 C2/R12ApdR22Ap*G4Exp(-aQ2/(4*R12Apd << 909 C3*R22Ap/2*G4Exp(-aQ2/8*R22Ap)); << 910 << 911 G4double DTot1 = 0.5*(C1*0.5*R12Ap-C2/R12Apd << 912 << 913 for(G4int i=1; i<= A-2; ++i) { << 914 N1p *= (-UnuclScr*Rho2*(A-i-1)/(G4double)i << 915 G4double N2p = 1; << 916 G4double Din2 = 0; << 917 G4double DmedTot = 0; << 918 G4double BinCoeff = 1.0; << 919 for(G4int l=0; l<=i; ++l) { << 920 if(l > 0) { BinCoeff *= (i-l+1)/(G4doubl << 921 << 922 exp1 = l/R22B+(i-l)/R12B; << 923 G4double exp1p = exp1+R12Apd; << 924 G4double exp2p = exp1+R12ApdR22Ap; << 925 G4double exp3p = exp1+R22Apd; << 926 << 927 Din2 += N2p*BinCoeff*(C1/exp1p*G4Exp(-aQ << 928 C2/exp2p*G4Exp(-aQ2/(4*exp2p))+ << 929 C3/exp3p*G4Exp(-aQ2/(4*exp3p))); << 930 << 931 DmedTot += N2p*BinCoeff*(C1/exp1p-C2/exp << 932 << 933 N2p *= -R23dR13; << 934 } // l << 935 << 936 G4double dcos = N1p*std::cos(FiH*i)/(G4dou << 937 Din1 += Din2*dcos; << 938 DTot1 += DmedTot*dcos; << 939 << 940 if(std::abs(Din2*N1p/Din1) < 0.000001) bre << 941 } // i << 942 G4double gg = (G4double)(A*(A-1)*4)/(Norm*No << 943 990 944 Din1 *= (-gg); << 991 BinCoeff = 1; 945 DTot1 *= 5*gg; << 946 992 947 // ---------------- dSigma/d|-t|, mb/(GeV << 993 for( i = 1; i<= Nucleus-2; i++) >> 994 { >> 995 N1p = -N1p*UnuclScr*(Nucleus-i-1)/i*Rho2; >> 996 N2p = 1; >> 997 Din2 = 0; >> 998 DmedTot = 0; >> 999 for(G4int l = 0; l<=i; l++) >> 1000 { >> 1001 if(l == 0) BinCoeff = 1; >> 1002 else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l; >> 1003 >> 1004 exp1 = l/R22B+(i-l)/R12B; >> 1005 exp1p = exp1+R12Apd; >> 1006 exp2p = exp1+R12ApdR22Ap; >> 1007 exp3p = exp1+R22Apd; >> 1008 >> 1009 Din2 = Din2 + N2p*BinCoeff* >> 1010 (C1/exp1p*std::exp(-theQ2/4/exp1p)- >> 1011 C2/exp2p*std::exp(-theQ2/4/exp2p)+ >> 1012 C3/exp3p*std::exp(-theQ2/4/exp3p)); 948 1013 949 G4double DiffCrSec2 = (ReElasticAmpl0*ReElas << 1014 DmedTot = DmedTot + N2p*BinCoeff* 950 (ImElasticAmpl0+Din1)* << 1015 (C1/exp1p-C2/exp2p+C3/exp3p); 951 (ImElasticAmpl0+Din1))/twopi; << 1016 >> 1017 N2p = -N2p*R23dR13; >> 1018 } // l >> 1019 >> 1020 Din1 = Din1+Din2*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i); >> 1021 DTot1 = DTot1+DmedTot*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i); >> 1022 >> 1023 if(std::fabs(Din2*N1p/Din1) < 0.000001) break; >> 1024 } // i 952 1025 953 Dtot11 = DTot1; << 1026 Din1 = -Din1*Nucleus*(Nucleus-1) 954 aAIm = ImElasticAmpl0; << 1027 /2/pi/Normp/2/pi/Normp*16*pi*pi; 955 aDIm = Din1; << 956 1028 957 return DiffCrSec2; // dSig/d|-t|, mb/(GeV << 1029 DTot1 = DTot1*Nucleus*(Nucleus-1) 958 } << 1030 /2/pi/Normp/2/pi/Normp*16*pi*pi; >> 1031 >> 1032 DTot1 *= 5; // $$$$$$$$$$$$$$$$$$$$$$$$ >> 1033 // Din1 *= 0.2; // %%%%%%%%%%%%%%%%%%%%%%% proton >> 1034 // Din1 *= 0.05; // %%%%%%%%%%%%%%%%%%%%%%% pi+ >> 1035 // ---------------- dSigma/d|-t|, mb/(GeV/c)^-2 ----------------- >> 1036 >> 1037 G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+ >> 1038 (ImElasticAmpl0+Din1)* >> 1039 (ImElasticAmpl0+Din1))*2/4/pi; >> 1040 >> 1041 Tot1 = Tot1-DTot1; >> 1042 // Tott1 = Tot1*1.0; >> 1043 Dtot11 = DTot1; >> 1044 aAIm = ImElasticAmpl0; >> 1045 aDIm = Din1; >> 1046 >> 1047 return DiffCrSec2*1.0; // dSig/d|-t|, mb/(GeV/c)^-2 >> 1048 } // function >> 1049 // ############################################## 959 1050 960 ////////////////////////////////////////////// 1051 //////////////////////////////////////////////////////////////// >> 1052 // >> 1053 // 961 1054 962 void G4ElasticHadrNucleusHE::DefineHadronValue << 1055 void G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z) 963 { 1056 { >> 1057 HadrEnergy = std::sqrt(hMass2 + hLabMomentum2); >> 1058 964 G4double sHadr = 2.*HadrEnergy*protonM+proto 1059 G4double sHadr = 2.*HadrEnergy*protonM+protonM2+hMass2; 965 G4double sqrS = std::sqrt(sHadr); 1060 G4double sqrS = std::sqrt(sHadr); >> 1061 G4double Ecm = 0.5*(sHadr-hMass2+protonM2)/sqrS; >> 1062 MomentumCM = std::sqrt(Ecm*Ecm-protonM2); 966 1063 967 if(verboseLevel>2) { << 1064 if(verboseLevel>2) 968 G4cout << "GetHadrValues: Z= " << Z << " i << 1065 G4cout << "GetHadrVall.: Z= " << Z << " iHadr= " << iHadron 969 << " E(GeV)= " << HadrEnergy << " sqrS= " 1066 << " E(GeV)= " << HadrEnergy << " sqrS= " << sqrS 970 << " plab= " << hLabMomentum 1067 << " plab= " << hLabMomentum 971 <<" E - m "<<HadrEnergy - hMass<< G4end 1068 <<" E - m "<<HadrEnergy - hMass<< G4endl; 972 } << 1069 973 G4double TotN = 0.0; 1070 G4double TotN = 0.0; 974 G4double logE = G4Log(HadrEnergy); << 1071 G4double logE = std::log(HadrEnergy); 975 G4double logS = G4Log(sHadr); << 1072 G4double logS = std::log(sHadr); 976 TotP = 0.0; 1073 TotP = 0.0; 977 1074 978 switch (iHadron) { << 1075 switch (iHadron) 979 case 0: // proton, neutron << 1076 { 980 case 6: << 1077 case 0: // proton, neutron 981 << 1078 case 6: 982 if(hLabMomentum > 10) { << 1079 983 TotP = TotN = 7.5*logE - 40.12525 + 103* << 1080 if(hLabMomentum > 10) 984 << 1081 TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165); // mb 985 } else { << 1082 986 // ================== neutron ======== << 1083 else 987 << 1084 { 988 if( hLabMomentum > 1.4 ) { << 1085 // ================== neutron ================ 989 TotN = 33.3+15.2*(hLabMomentum2-1.35)/ << 1086 990 (G4Exp(G4Log(hLabMomentum)*2.37)+0.95); << 1087 //// if(iHadrCode == 2112) >> 1088 >> 1089 >> 1090 if( hLabMomentum > 1.4 ) >> 1091 TotN = 33.3+15.2*(hLabMomentum2-1.35)/ >> 1092 (std::pow(hLabMomentum,2.37)+0.95); 991 1093 992 } else if(hLabMomentum > 0.8) { << 1094 else if(hLabMomentum > 0.8) 993 G4double A0 = logE + 0.0513; << 1095 { 994 TotN = 33.0 + 25.5*A0*A0; << 1096 G4double A0 = logE + 0.0513; 995 } else { << 1097 TotN = 33.0 + 25.5*A0*A0; 996 G4double A0 = logE - 0.2634; // log(1.3) << 1098 } 997 TotN = 33.0 + 30.*A0*A0*A0*A0; << 1099 else 998 } << 1100 { 999 // ================= proton ========= << 1101 G4double A0 = logE - 0.2634; // log(1.3) >> 1102 TotN = 33.0 + 30.*A0*A0*A0*A0; >> 1103 } >> 1104 // ================= proton =============== >> 1105 // else if(iHadrCode == 2212) >> 1106 { >> 1107 if(hLabMomentum >= 1.05) >> 1108 { >> 1109 TotP = 39.0+75.*(hLabMomentum-1.2)/ >> 1110 (hLabMomentum2*hLabMomentum+0.15); >> 1111 } >> 1112 >> 1113 else if(hLabMomentum >= 0.7) >> 1114 { >> 1115 G4double A0 = logE + 0.3147; >> 1116 TotP = 23.0 + 40.*A0*A0; >> 1117 } >> 1118 else >> 1119 { >> 1120 TotP = 23.+50.*std::pow(std::log(0.73/hLabMomentum),3.5); >> 1121 } >> 1122 } >> 1123 } 1000 1124 1001 if(hLabMomentum >= 1.05) { << 1125 // HadrTot = 0.5*(82*TotP+126*TotN)/104; // $$$$$$$$$$$$$$$$$$ 1002 TotP = 39.0+75.*(hLabMomentum-1.2)/(hLabMom << 1126 HadrTot = 0.5*(TotP+TotN); 1003 } else if(hLabMomentum >= 0.7) { << 1127 // ................................................... 1004 G4double A0 = logE + 0.3147; << 1128 // Proton slope 1005 TotP = 23.0 + 40.*A0*A0; << 1129 if(hLabMomentum >= 2.) HadrSlope = 5.44 + 0.88*logS; 1006 } else { << 1130 1007 TotP = 23.+50.*G4Exp(G4Log(G4Log(0.73/hLabM << 1131 else if(hLabMomentum >= 0.5) HadrSlope = 3.73*hLabMomentum-0.37; 1008 } << 1132 1009 } << 1133 else HadrSlope = 1.5; 1010 HadrTot = 0.5*(TotP+TotN); << 1134 1011 // ..................................... << 1135 // ................................................... 1012 // Proton slope << 1136 if(hLabMomentum >= 1.2) 1013 if(hLabMomentum >= 2.) { HadrSlope << 1137 HadrReIm = 0.13*(logS - 5.8579332)*std::pow(sHadr,-0.18); 1014 else if(hLabMomentum >= 0.5) { HadrSlope << 1138 1015 else { HadrSlope << 1139 else if(hLabMomentum >= 0.6) 1016 << 1140 HadrReIm = -75.5*(std::pow(hLabMomentum,0.25)-0.95)/ 1017 // ..................................... << 1141 (std::pow(3*hLabMomentum,2.2)+1); 1018 if(hLabMomentum >= 1.2) { << 1142 1019 HadrReIm = 0.13*(logS - 5.8579332)*G4E << 1143 else 1020 } else if(hLabMomentum >= 0.6) { << 1144 HadrReIm = 15.5*hLabMomentum/(27*hLabMomentum2*hLabMomentum+2); 1021 HadrReIm = -75.5*(G4Exp(G4Log(hLabMomen << 1145 // ................................................... 1022 (G4Exp(G4Log(3*hLabMomentum)*2.2)+1); << 1146 DDSect2 = 2.2; //mb*GeV-2 1023 } else { << 1147 DDSect3 = 0.6; //mb*GeV-2 1024 HadrReIm = 15.5*hLabMomentum/(27*hLabMo << 1148 // ================== lambda ================== 1025 } << 1149 if( iHadrCode == 3122) 1026 // ..................................... << 1150 { 1027 DDSect2 = 2.2; << 1151 HadrTot *= 0.88; 1028 DDSect3 = 0.6; << 1152 HadrSlope *=0.85; 1029 // ================== lambda ========== << 1153 } 1030 if( iHadrCode == 3122) { << 1031 HadrTot *= 0.88; << 1032 HadrSlope *=0.85; << 1033 // ================== sigma + ======= 1154 // ================== sigma + ================== 1034 } else if( iHadrCode == 3222) { << 1155 else if( iHadrCode == 3222) 1035 HadrTot *=0.81; << 1156 { 1036 HadrSlope *=0.85; << 1157 HadrTot *=0.81; >> 1158 HadrSlope *=0.85; >> 1159 } 1037 // ================== sigma 0,- ===== 1160 // ================== sigma 0,- ================== 1038 } else if(iHadrCode == 3112 || iHadrCode << 1161 else if(iHadrCode == 3112 || iHadrCode == 3212 ) 1039 HadrTot *=0.88; << 1162 { 1040 HadrSlope *=0.85; << 1163 HadrTot *=0.88; >> 1164 HadrSlope *=0.85; >> 1165 } 1041 // =================== xi ========== 1166 // =================== xi ================= 1042 } else if( iHadrCode == 3312 || iHadrCode << 1167 else if( iHadrCode == 3312 || iHadrCode == 3322 ) 1043 HadrTot *=0.77; << 1168 { 1044 HadrSlope *=0.75; << 1169 HadrTot *=0.77; >> 1170 HadrSlope *=0.75; >> 1171 } 1045 // ================= omega ========= 1172 // ================= omega ================= 1046 } else if( iHadrCode == 3334) { << 1173 else if( iHadrCode == 3334) 1047 HadrTot *=0.78; << 1174 { 1048 HadrSlope *=0.7; << 1175 HadrTot *=0.78; 1049 } << 1176 HadrSlope *=0.7; 1050 break; << 1177 } 1051 // ===================================== << 1178 1052 case 1: // antiproton << 1179 break; 1053 case 7: // antineutron << 1180 // =========================================================== 1054 << 1181 case 1: // antiproton 1055 HadrTot = 5.2+5.2*logE + 123.2/sqrS; << 1182 case 7: // antineutron 1056 HadrSlope = 8.32+0.57*logS; << 1183 1057 << 1184 HadrTot = 5.2+5.2*logE + 123.2/sqrS; // mb 1058 if( HadrEnergy < 1000 ) { << 1185 HadrSlope = 8.32+0.57*logS; //(GeV/c)^-2 1059 HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14. << 1186 1060 } else { << 1187 if( HadrEnergy < 1000 ) 1061 HadrReIm = 0.6*(logS - 5.8579332)*G4Ex << 1188 HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8); 1062 } << 1189 else 1063 DDSect2 = 11; << 1190 HadrReIm = 0.6*(logS - 5.8579332)*std::pow(sHadr,-0.25); 1064 DDSect3 = 3; << 1191 1065 // ================== lambda ========== << 1192 DDSect2 = 11; //mb*(GeV/c)^-2 1066 if( iHadrCode == -3122) { << 1193 DDSect3 = 3; //mb*(GeV/c)^-2 1067 HadrTot *= 0.88; << 1194 // ================== lambda ================== 1068 HadrSlope *=0.85; << 1195 if( iHadrCode == -3122) >> 1196 { >> 1197 HadrTot *= 0.88; >> 1198 HadrSlope *=0.85; >> 1199 } 1069 // ================== sigma + ======= 1200 // ================== sigma + ================== 1070 } else if( iHadrCode == -3222) { << 1201 else if( iHadrCode == -3222) 1071 HadrTot *=0.81; << 1202 { 1072 HadrSlope *=0.85; << 1203 HadrTot *=0.81; >> 1204 HadrSlope *=0.85; >> 1205 } 1073 // ================== sigma 0,- ===== 1206 // ================== sigma 0,- ================== 1074 } else if(iHadrCode == -3112 || iHadrCode << 1207 else if(iHadrCode == -3112 || iHadrCode == -3212 ) 1075 HadrTot *=0.88; << 1208 { 1076 HadrSlope *=0.85; << 1209 HadrTot *=0.88; 1077 // =================== xi ============ << 1210 HadrSlope *=0.85; 1078 } else if( iHadrCode == -3312 || iHadrCod << 1211 } 1079 HadrTot *=0.77; << 1212 // =================== xi ================= 1080 HadrSlope *=0.75; << 1213 else if( iHadrCode == -3312 || iHadrCode == -3322 ) >> 1214 { >> 1215 HadrTot *=0.77; >> 1216 HadrSlope *=0.75; >> 1217 } 1081 // ================= omega ========= 1218 // ================= omega ================= 1082 } else if( iHadrCode == -3334) { << 1219 else if( iHadrCode == -3334) 1083 HadrTot *=0.78; << 1220 { 1084 HadrSlope *=0.7; << 1221 HadrTot *=0.78; 1085 } << 1222 HadrSlope *=0.7; 1086 break; << 1223 } 1087 // ------------------------------------- << 1088 case 2: // pi plus, pi minus << 1089 case 3: << 1090 << 1091 if(hLabMomentum >= 3.5) { << 1092 TotP = 10.6+2.*logE + 25.*G4Exp(-logE*0 << 1093 // =================================== << 1094 } else if(hLabMomentum >= 1.15) { << 1095 G4double x = (hLabMomentum - 2.55)/0.55 << 1096 G4double y = (hLabMomentum - 1.47)/0.22 << 1097 TotP = 3.2*G4Exp(-x*x) + 12.*G4Exp(-y*y << 1098 // =================================== << 1099 } else if(hLabMomentum >= 0.4) { << 1100 TotP = 88*(logE+0.2877)*(logE+0.2877)+ << 1101 // ===================================== << 1102 } else { << 1103 G4double x = (hLabMomentum - 0.29)/0.08 << 1104 TotP = 20. + 180.*G4Exp(-x*x); << 1105 } << 1106 // ------------------------------------- << 1107 1224 1108 if(hLabMomentum >= 3.0 ) { << 1225 break; 1109 TotN = 10.6 + 2.*logE + 30.*G4Exp(-logE << 1226 // ------------------------------------------- 1110 } else if(hLabMomentum >= 1.3) { << 1227 case 2: // pi plus, pi minus 1111 G4double x = (hLabMomentum - 2.1)/0.4; << 1228 case 3: 1112 G4double y = (hLabMomentum - 1.4)/0.12; << 1229 1113 TotN = 36.1+0.079 - 4.313*logE + 3.*G4E << 1230 if(hLabMomentum >= 3.5) 1114 } else if(hLabMomentum >= 0.65) { << 1231 TotP = 10.6+2.*logE + 25.*std::pow(HadrEnergy,-0.43); // mb 1115 G4double x = (hLabMomentum - 0.72)/0.06 << 1232 // ========================================= 1116 G4double y = (hLabMomentum - 1.015)/0.0 << 1233 else if(hLabMomentum >= 1.15) 1117 TotN = 36.1 + 10.*G4Exp(-x*x) + 24*G4Ex << 1234 { 1118 } else if(hLabMomentum >= 0.37) { << 1235 G4double x = (hLabMomentum - 2.55)/0.55; 1119 G4double x = G4Log(hLabMomentum/0.48); << 1236 G4double y = (hLabMomentum - 1.47)/0.225; 1120 TotN = 26. + 110.*x*x; << 1237 TotP = 3.2*std::exp(-x*x) + 12.*std::exp(-y*y) + 27.5; 1121 } else { << 1238 } 1122 G4double x = (hLabMomentum - 0.29)/0.07 << 1239 // ========================================= 1123 TotN = 28.0 + 40.*G4Exp(-x*x); << 1240 else if(hLabMomentum >= 0.4) 1124 } << 1241 { 1125 HadrTot = (TotP+TotN)*0.5; << 1242 TotP = 88*(logE+0.2877)*(logE+0.2877)+14.0; 1126 // ..................................... << 1243 } 1127 HadrSlope = 7.28+0.245*logS; // Ge << 1244 // ========================================= 1128 HadrReIm = 0.2*(logS - 4.6051702)*G4Exp( << 1245 else 1129 << 1246 { 1130 DDSect2 = 0.7; << 1247 G4double x = (hLabMomentum - 0.29)/0.085; 1131 DDSect3 = 0.27; << 1248 TotP = 20. + 180.*std::exp(-x*x); 1132 << 1249 } 1133 break; << 1250 // ------------------------------------------- 1134 // ===================================== << 1251 1135 case 4: // K plus << 1252 if(hLabMomentum >= 3.0 ) 1136 << 1253 TotN = 10.6 + 2.*logE + 30.*std::pow(HadrEnergy,-0.43); // mb 1137 HadrTot = 10.6+1.8*logE + 9.0*G4Exp(-lo << 1254 1138 if(HadrEnergy>100) { HadrSlope = 15.0; } << 1255 else if(hLabMomentum >= 1.3) 1139 else { HadrSlope = 1.0+1.76*logS - 2.84/s << 1256 { 1140 << 1257 G4double x = (hLabMomentum - 2.1)/0.4; 1141 HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*G4 << 1258 G4double y = (hLabMomentum - 1.4)/0.12; 1142 DDSect2 = 0.7; << 1259 TotN = 36.1+0.079 - 4.313*logE + 3.*std::exp(-x*x) + 1143 DDSect3 = 0.21; << 1260 1.5*std::exp(-y*y); 1144 break; << 1261 } 1145 // ===================================== << 1262 else if(hLabMomentum >= 0.65) 1146 case 5: // K minus << 1263 { 1147 << 1264 G4double x = (hLabMomentum - 0.72)/0.06; 1148 HadrTot = 10+1.8*logE + 25./sqrS; // mb << 1265 G4double y = (hLabMomentum - 1.015)/0.075; 1149 HadrSlope = 6.98+0.127*logS; // G << 1266 TotN = 36.1 + 10.*std::exp(-x*x) + 24*std::exp(-y*y); 1150 HadrReIm = 0.4*(sHadr-20)*(sHadr-20)*G4E << 1267 } 1151 DDSect2 = 0.7; << 1268 else if(hLabMomentum >= 0.37) 1152 DDSect3 = 0.27; << 1269 { 1153 break; << 1270 G4double x = std::log(hLabMomentum/0.48); >> 1271 TotN = 26. + 110.*x*x; >> 1272 } >> 1273 else >> 1274 { >> 1275 G4double x = (hLabMomentum - 0.29)/0.07; >> 1276 TotN = 28.0 + 40.*std::exp(-x*x); >> 1277 } >> 1278 HadrTot = (TotP+TotN)/2; >> 1279 // ........................................ >> 1280 HadrSlope = 7.28+0.245*logS; // GeV-2 >> 1281 HadrReIm = 0.2*(logS - 4.6051702)*std::pow(sHadr,-0.15); >> 1282 >> 1283 DDSect2 = 0.7; //mb*GeV-2 >> 1284 DDSect3 = 0.27; //mb*GeV-2 >> 1285 >> 1286 break; >> 1287 // ========================================================== >> 1288 case 4: // K plus >> 1289 >> 1290 HadrTot = 10.6+1.8*logE + 9.0*std::pow(HadrEnergy,-0.55); // mb >> 1291 if(HadrEnergy>100) HadrSlope = 15.0; >> 1292 else HadrSlope = 1.0+1.76*logS - 2.84/sqrS; // GeV-2 >> 1293 >> 1294 HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1); >> 1295 DDSect2 = 0.7; //mb*GeV-2 >> 1296 DDSect3 = 0.21; //mb*GeV-2 >> 1297 break; >> 1298 // ========================================================= >> 1299 case 5: // K minus >> 1300 >> 1301 HadrTot = 10+1.8*logE + 25./sqrS; // mb >> 1302 HadrSlope = 6.98+0.127*logS; // GeV-2 >> 1303 HadrReIm = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1); >> 1304 DDSect2 = 0.7; //mb*GeV-2 >> 1305 DDSect3 = 0.27; //mb*GeV-2 >> 1306 break; 1154 } 1307 } 1155 // ======================================= << 1308 // ========================================================= 1156 if(verboseLevel>2) { << 1309 if(verboseLevel>2) 1157 G4cout << "HadrTot= " << HadrTot << " Had 1310 G4cout << "HadrTot= " << HadrTot << " HadrSlope= " << HadrSlope 1158 << " HadrReIm= " << HadrReIm << " DDSect 1311 << " HadrReIm= " << HadrReIm << " DDSect2= " << DDSect2 1159 << " DDSect3= " << DDSect3 << G4endl; 1312 << " DDSect3= " << DDSect3 << G4endl; 1160 } << 1313 1161 if(Z != 1) return; 1314 if(Z != 1) return; 1162 1315 1163 // Scattering of protons 1316 // Scattering of protons 1164 1317 1165 Coeff0 = Coeff1 = Coeff2 = 0.0; 1318 Coeff0 = Coeff1 = Coeff2 = 0.0; 1166 Slope0 = Slope1 = 1.0; 1319 Slope0 = Slope1 = 1.0; 1167 Slope2 = 5.0; 1320 Slope2 = 5.0; 1168 1321 1169 // data for iHadron=0 1322 // data for iHadron=0 1170 static const G4double EnP0[6]={1.5,3.0,5.0, 1323 static const G4double EnP0[6]={1.5,3.0,5.0,9.0,14.0,19.0}; 1171 static const G4double C0P0[6]={0.15,0.02,0. 1324 static const G4double C0P0[6]={0.15,0.02,0.06,0.08,0.0003,0.0002}; 1172 static const G4double C1P0[6]={0.05,0.02,0. 1325 static const G4double C1P0[6]={0.05,0.02,0.03,0.025,0.0,0.0}; 1173 static const G4double B0P0[6]={1.5,2.5,3.0, 1326 static const G4double B0P0[6]={1.5,2.5,3.0,4.5,1.4,1.25}; 1174 static const G4double B1P0[6]={5.0,1.0,3.5, 1327 static const G4double B1P0[6]={5.0,1.0,3.5,4.0,4.8,4.8}; 1175 1328 1176 // data for iHadron=6,7 1329 // data for iHadron=6,7 1177 static const G4double EnN[5]={1.5,5.0,10.0, 1330 static const G4double EnN[5]={1.5,5.0,10.0,14.0,20.0}; 1178 static const G4double C0N[5]={0.0,0.0,0.02, 1331 static const G4double C0N[5]={0.0,0.0,0.02,0.02,0.01}; 1179 static const G4double C1N[5]={0.06,0.008,0. 1332 static const G4double C1N[5]={0.06,0.008,0.0015,0.001,0.0003}; 1180 static const G4double B0N[5]={1.5,2.5,3.8,3 1333 static const G4double B0N[5]={1.5,2.5,3.8,3.8,3.5}; 1181 static const G4double B1N[5]={1.5,2.2,3.6,4 1334 static const G4double B1N[5]={1.5,2.2,3.6,4.5,4.8}; 1182 1335 1183 // data for iHadron=1 1336 // data for iHadron=1 1184 static const G4double EnP[2]={1.5,4.0}; 1337 static const G4double EnP[2]={1.5,4.0}; 1185 static const G4double C0P[2]={0.001,0.0005} 1338 static const G4double C0P[2]={0.001,0.0005}; 1186 static const G4double C1P[2]={0.003,0.001}; 1339 static const G4double C1P[2]={0.003,0.001}; 1187 static const G4double B0P[2]={2.5,4.5}; 1340 static const G4double B0P[2]={2.5,4.5}; 1188 static const G4double B1P[2]={1.0,4.0}; 1341 static const G4double B1P[2]={1.0,4.0}; 1189 1342 1190 // data for iHadron=2 1343 // data for iHadron=2 1191 static const G4double EnPP[4]={1.0,2.0,3.0, 1344 static const G4double EnPP[4]={1.0,2.0,3.0,4.0}; 1192 static const G4double C0PP[4]={0.0,0.0,0.0, 1345 static const G4double C0PP[4]={0.0,0.0,0.0,0.0}; 1193 static const G4double C1PP[4]={0.15,0.08,0. 1346 static const G4double C1PP[4]={0.15,0.08,0.02,0.01}; 1194 static const G4double B0PP[4]={1.5,2.8,3.8, 1347 static const G4double B0PP[4]={1.5,2.8,3.8,3.8}; 1195 static const G4double B1PP[4]={0.8,1.6,3.6, 1348 static const G4double B1PP[4]={0.8,1.6,3.6,4.6}; 1196 1349 1197 // data for iHadron=3 1350 // data for iHadron=3 1198 static const G4double EnPPN[4]={1.0,2.0,3.0 1351 static const G4double EnPPN[4]={1.0,2.0,3.0,4.0}; 1199 static const G4double C0PPN[4]={0.0,0.0,0.0 1352 static const G4double C0PPN[4]={0.0,0.0,0.0,0.0}; 1200 static const G4double C1PPN[4]={0.0,0.0,0.0 1353 static const G4double C1PPN[4]={0.0,0.0,0.0,0.0}; 1201 static const G4double B0PPN[4]={1.5,2.8,3.8 1354 static const G4double B0PPN[4]={1.5,2.8,3.8,3.8}; 1202 static const G4double B1PPN[4]={0.8,1.6,3.6 1355 static const G4double B1PPN[4]={0.8,1.6,3.6,4.6}; 1203 1356 1204 // data for iHadron=4 1357 // data for iHadron=4 1205 static const G4double EnK[4]={1.4,2.33,3.0, 1358 static const G4double EnK[4]={1.4,2.33,3.0,5.0}; 1206 static const G4double C0K[4]={0.0,0.0,0.0,0 1359 static const G4double C0K[4]={0.0,0.0,0.0,0.0}; 1207 static const G4double C1K[4]={0.01,0.007,0. 1360 static const G4double C1K[4]={0.01,0.007,0.005,0.003}; 1208 static const G4double B0K[4]={1.5,2.0,3.8,3 1361 static const G4double B0K[4]={1.5,2.0,3.8,3.8}; 1209 static const G4double B1K[4]={1.6,1.6,1.6,1 1362 static const G4double B1K[4]={1.6,1.6,1.6,1.6}; 1210 1363 1211 // data for iHadron=5 1364 // data for iHadron=5 1212 static const G4double EnKM[2]={1.4,4.0}; 1365 static const G4double EnKM[2]={1.4,4.0}; 1213 static const G4double C0KM[2]={0.006,0.002} 1366 static const G4double C0KM[2]={0.006,0.002}; 1214 static const G4double C1KM[2]={0.00,0.00}; 1367 static const G4double C1KM[2]={0.00,0.00}; 1215 static const G4double B0KM[2]={2.5,3.5}; 1368 static const G4double B0KM[2]={2.5,3.5}; 1216 static const G4double B1KM[2]={1.6,1.6}; 1369 static const G4double B1KM[2]={1.6,1.6}; 1217 1370 1218 switch(iHadron) { << 1371 switch(iHadron) 1219 case 0: << 1372 { >> 1373 case 0 : 1220 1374 1221 if(hLabMomentum <BoundaryP[0]) { << 1375 if(hLabMomentum <BoundaryP[0]) 1222 InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P << 1376 InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P0); 1223 } << 1224 Coeff2 = 0.8/hLabMomentum2; << 1225 break; << 1226 1377 1227 case 6: << 1378 Coeff2 = 0.8/hLabMomentum/hLabMomentum; >> 1379 break; 1228 1380 1229 if(hLabMomentum < BoundaryP[1]) { << 1381 case 6 : 1230 InterpolateHN(5,EnN,C0N,C1N,B0N,B1N); << 1231 } << 1232 Coeff2 = 0.8/hLabMomentum2; << 1233 break; << 1234 1382 1235 case 1: << 1383 if(hLabMomentum < BoundaryP[1]) 1236 case 7: << 1384 InterpolateHN(5,EnN,C0N,C1N,B0N,B1N); 1237 if(hLabMomentum < BoundaryP[2]) { << 1385 1238 InterpolateHN(2,EnP,C0P,C1P,B0P,B1P); << 1386 Coeff2 = 0.8/hLabMomentum/hLabMomentum; 1239 } << 1387 break; 1240 break; << 1241 1388 1242 case 2: << 1389 case 1 : >> 1390 case 7 : >> 1391 if(hLabMomentum < BoundaryP[2]) >> 1392 InterpolateHN(2,EnP,C0P,C1P,B0P,B1P); >> 1393 break; 1243 1394 1244 if(hLabMomentum < BoundaryP[3]) { << 1395 case 2 : 1245 InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1P << 1246 } << 1247 Coeff2 = 0.02/hLabMomentum; << 1248 break; << 1249 1396 1250 case 3: << 1397 if(hLabMomentum < BoundaryP[3]) >> 1398 InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1PP); 1251 1399 1252 if(hLabMomentum < BoundaryP[4]) { << 1400 Coeff2 = 0.02/hLabMomentum; 1253 InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN << 1401 break; 1254 } << 1402 1255 Coeff2 = 0.02/hLabMomentum; << 1403 case 3 : 1256 break; << 1404 >> 1405 if(hLabMomentum < BoundaryP[4]) >> 1406 InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN,B1PPN); >> 1407 >> 1408 Coeff2 = 0.02/hLabMomentum; >> 1409 break; 1257 1410 1258 case 4: << 1411 case 4 : 1259 1412 1260 if(hLabMomentum < BoundaryP[5]) { << 1413 if(hLabMomentum < BoundaryP[5]) 1261 InterpolateHN(4,EnK,C0K,C1K,B0K,B1K); << 1414 InterpolateHN(4,EnK,C0K,C1K,B0K,B1K); 1262 } << 1415 1263 if(hLabMomentum < 1) { Coeff2 = 0.34; } << 1416 if(hLabMomentum < 1) Coeff2 = 0.34; 1264 else { Coeff2 = 0.34/(hLabMomentum2*hLab << 1417 else Coeff2 = 0.34/hLabMomentum2/hLabMomentum; 1265 break; << 1418 break; 1266 << 1419 1267 case 5: << 1420 case 5 : 1268 if(hLabMomentum < BoundaryP[6]) { << 1421 if(hLabMomentum < BoundaryP[6]) 1269 InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1K << 1422 InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1KM); >> 1423 >> 1424 if(hLabMomentum < 1) Coeff2 = 0.01; >> 1425 else Coeff2 = 0.01/hLabMomentum2/hLabMomentum; >> 1426 break; 1270 } 1427 } 1271 if(hLabMomentum < 1) { Coeff2 = 0.01; } << 1272 else { Coeff2 = 0.01/(hLabMomentum2*hLab << 1273 break; << 1274 } << 1275 1428 1276 if(verboseLevel > 2) { << 1429 if(verboseLevel > 2) 1277 G4cout<<" HadrVal : Plasb "<<hLabMoment 1430 G4cout<<" HadrVal : Plasb "<<hLabMomentum 1278 <<" iHadron "<<iHadron<<" HadrTot "<< 1431 <<" iHadron "<<iHadron<<" HadrTot "<<HadrTot<<G4endl; 1279 } << 1280 } 1432 } 1281 1433 1282 ///////////////////////////////////////////// << 1434 // ===================================================== >> 1435 void G4ElasticHadrNucleusHE:: >> 1436 GetKinematics(const G4ParticleDefinition * aHadron, >> 1437 G4double MomentumH) >> 1438 { >> 1439 if (verboseLevel>1) >> 1440 G4cout<<"1 GetKin.: HadronName MomentumH " >> 1441 <<aHadron->GetParticleName()<<" "<<MomentumH<<G4endl; >> 1442 >> 1443 DefineHadronValues(1); >> 1444 >> 1445 G4double Sh = 2.0*protonM*HadrEnergy+protonM2+hMass2; // GeV >> 1446 >> 1447 ConstU = 2*protonM2+2*hMass2-Sh; >> 1448 >> 1449 G4double MaxT = 4*MomentumCM*MomentumCM; >> 1450 >> 1451 BoundaryTL[0] = MaxT; //2.0; >> 1452 BoundaryTL[1] = MaxT; >> 1453 BoundaryTL[3] = MaxT; >> 1454 BoundaryTL[4] = MaxT; >> 1455 BoundaryTL[5] = MaxT; >> 1456 >> 1457 G4int NumberH=0; >> 1458 >> 1459 while(iHadrCode!=HadronCode[NumberH]) NumberH++; >> 1460 >> 1461 NumberH = HadronType1[NumberH]; 1283 1462 >> 1463 if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH]; >> 1464 else MaxTR = BoundaryTG[NumberH]; >> 1465 >> 1466 if (verboseLevel>1) >> 1467 G4cout<<"3 GetKin. : NumberH "<<NumberH >> 1468 <<" Bound.P[NumberH] "<<BoundaryP[NumberH] >> 1469 <<" Bound.TL[NumberH] "<<BoundaryTL[NumberH] >> 1470 <<" Bound.TG[NumberH] "<<BoundaryTG[NumberH] >> 1471 <<" MaxT MaxTR "<<MaxT<<" "<<MaxTR<<G4endl; >> 1472 >> 1473 // GetParametersHP(aHadron, MomentumH); >> 1474 } >> 1475 // ============================================================ 1284 G4double G4ElasticHadrNucleusHE::GetFt(G4doub 1476 G4double G4ElasticHadrNucleusHE::GetFt(G4double Q2) 1285 { 1477 { 1286 G4double Fdistr=0; 1478 G4double Fdistr=0; 1287 G4double SqrQ2 = std::sqrt(Q2); 1479 G4double SqrQ2 = std::sqrt(Q2); 1288 1480 1289 Fdistr = (1-Coeff1-Coeff0) / HadrSlope*(1-G << 1481 Fdistr = (1-Coeff1-Coeff0) //-0.0*Coeff2*std::exp(ConstU)) 1290 + Coeff0*(1-G4Exp(-Slope0*Q2)) << 1482 /HadrSlope*(1-std::exp(-HadrSlope*Q2)) 1291 + Coeff2/Slope2*G4Exp(Slope2*ConstU)*(G4E << 1483 + Coeff0*(1-std::exp(-Slope0*Q2)) 1292 + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+Sqr << 1484 + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1) >> 1485 + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2)); 1293 1486 1294 if (verboseLevel>1) { << 1487 if (verboseLevel>1) 1295 G4cout<<"Old: Coeff0 Coeff1 Coeff2 "<<Co 1488 G4cout<<"Old: Coeff0 Coeff1 Coeff2 "<<Coeff0<<" " 1296 <<Coeff1<<" "<<Coeff2<<" Slope Sl 1489 <<Coeff1<<" "<<Coeff2<<" Slope Slope0 Slope1 Slope2 " 1297 <<HadrSlope<<" "<<Slope0<<" "<<Sl 1490 <<HadrSlope<<" "<<Slope0<<" "<<Slope1<<" "<<Slope2 1298 <<" Fdistr "<<Fdistr<<G4endl; 1491 <<" Fdistr "<<Fdistr<<G4endl; 1299 } << 1300 return Fdistr; 1492 return Fdistr; 1301 } 1493 } 1302 << 1494 // +++++++++++++++++++++++++++++++++++++++ 1303 ///////////////////////////////////////////// << 1495 G4double G4ElasticHadrNucleusHE::GetQ2(G4double Ran) 1304 << 1305 G4double << 1306 G4ElasticHadrNucleusHE::HadronProtonQ2(G4doub << 1307 { 1496 { 1308 hLabMomentum = plab; << 1497 G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta; 1309 hLabMomentum2 = hLabMomentum*hLabMomentum; << 1498 G4double Q2=0; 1310 HadrEnergy = std::sqrt(hMass2 + hLabMome << 1311 DefineHadronValues(1); << 1312 1499 1313 G4double Sh = 2.0*protonM*HadrEnergy+proton << 1500 FmaxT = GetFt(MaxTR); 1314 ConstU = 2*protonM2+2*hMass2-Sh; << 1501 delta = GetDistrFun(DDD0)-Ran; 1315 1502 1316 BoundaryTL[0] = tmax; << 1503 while(std::fabs(delta) > 0.0001) 1317 BoundaryTL[1] = tmax; << 1318 BoundaryTL[3] = tmax; << 1319 BoundaryTL[4] = tmax; << 1320 BoundaryTL[5] = tmax; << 1321 << 1322 G4double MaxTR = (plab < BoundaryP[iHadron1 << 1323 BoundaryTL[iHadron1] : BoundaryTG[iHadron << 1324 << 1325 if (verboseLevel>1) { << 1326 G4cout<<"3 GetKin. : iHadron1 "<<iHadro << 1327 <<" Bound.P[iHadron1] "<<BoundaryP[iHadr << 1328 <<" Bound.TL[iHadron1] "<<BoundaryTL[iHa << 1329 <<" Bound.TG[iHadron1] "<<BoundaryTG[iHa << 1330 <<" MaxT MaxTR "<<tmax<<" "<<MaxTR<<G4e << 1331 } << 1332 << 1333 G4double rand = G4UniformRand(); << 1334 << 1335 G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=Max << 1336 << 1337 G4double norm = 1.0/GetFt(MaxTR); << 1338 G4double delta = GetFt(DDD0)*norm - rand; << 1339 << 1340 static const G4int maxNumberOfLoops = 10000 << 1341 G4int loopCounter = -1; << 1342 while ( (std::abs(delta) > 0.0001) && << 1343 ++loopCounter < maxNumberOfLoops ) << 1344 { 1504 { 1345 if(delta>0) 1505 if(delta>0) 1346 { 1506 { 1347 DDD2 = DDD0; 1507 DDD2 = DDD0; 1348 DDD0 = (DDD0+DDD1)*0.5; 1508 DDD0 = (DDD0+DDD1)*0.5; 1349 } 1509 } 1350 else if(delta<0.0) << 1510 else if(delta<0) 1351 { 1511 { 1352 DDD1 = DDD0; 1512 DDD1 = DDD0; 1353 DDD0 = (DDD0+DDD2)*0.5; 1513 DDD0 = (DDD0+DDD2)*0.5; 1354 } 1514 } 1355 delta = GetFt(DDD0)*norm - rand; << 1515 delta = GetDistrFun(DDD0)-Ran; 1356 } 1516 } 1357 return (loopCounter >= maxNumberOfLoops) ? << 1358 } << 1359 1517 1360 ///////////////////////////////////////////// << 1518 Q2 = DDD0; 1361 1519 1362 void G4ElasticHadrNucleusHE::Binom() << 1520 return Q2; 1363 { << 1364 for(G4int N = 0; N < 240; ++N) { << 1365 G4double J = 1.0; << 1366 for(G4int M = 0; M <= N; ++M) { << 1367 G4double Fact1 = 1.0; << 1368 if (N > 0 && N > M && M > 0 ) { << 1369 J *= (G4double)(N-M+1)/(G4double)M; << 1370 Fact1 = J; << 1371 } << 1372 fBinom[N][M] = Fact1; << 1373 } << 1374 } << 1375 } 1521 } >> 1522 // ++++++++++++++++++++++++++++++++++++++++++ >> 1523 G4double G4ElasticHadrNucleusHE:: >> 1524 HadronProtonQ2(const G4ParticleDefinition * p, >> 1525 G4double inLabMom) >> 1526 { 1376 1527 1377 ///////////////////////////////////////////// << 1528 hMass = p->GetPDGMass()/GeV; >> 1529 hMass2 = hMass*hMass; >> 1530 hLabMomentum = inLabMom; >> 1531 hLabMomentum2 = hLabMomentum*hLabMomentum; >> 1532 HadrEnergy = sqrt(hLabMomentum2+hMass2); 1378 1533 1379 void << 1534 G4double Rand = G4UniformRand(); 1380 G4ElasticHadrNucleusHE::InFileName(std::ostri << 1381 const G4ParticleDefinition* p, G4i << 1382 { << 1383 if(!fDirectory) { << 1384 fDirectory = G4FindDataDir("G4LEDATA"); << 1385 if (fDirectory) { << 1386 ss << fDirectory << "/"; << 1387 } << 1388 } << 1389 OutFileName(ss, p, Z); << 1390 } << 1391 1535 1392 ///////////////////////////////////////////// << 1536 GetKinematics(p, inLabMom); 1393 1537 1394 void << 1538 G4double Q2 = GetQ2(Rand); 1395 G4ElasticHadrNucleusHE::OutFileName(std::ostr << 1539 1396 const G4ParticleDefinition* p, G4 << 1540 return Q2; 1397 { << 1398 ss << "hedata/" << p->GetParticleName() << << 1399 } 1541 } 1400 1542 1401 ///////////////////////////////////////////// << 1543 // =========================================== 1402 1544 1403 G4bool G4ElasticHadrNucleusHE::ReadLine(std:: << 1545 /////////////////////////////////////////////////////////////////// 1404 std::vector<G4double>& v) << 1546 // >> 1547 // >> 1548 >> 1549 void G4ElasticHadrNucleusHE::Binom() 1405 { 1550 { 1406 G4int n(0); << 1551 G4int N, M; 1407 infile >> n; << 1552 G4double Fact1, J; 1408 if (infile.fail()) { return false; } << 1409 if(n > 0) { << 1410 v.reserve(n); << 1411 G4double x(0.0); << 1412 for(G4int i=0; i<n; ++i) { << 1413 infile >> x; << 1414 if (infile.fail()) { return false; } << 1415 v.emplace_back(x); << 1416 } << 1417 } << 1418 return true; << 1419 } << 1420 1553 1421 ///////////////////////////////////////////// << 1554 for(N = 0; N < 240; N++) >> 1555 { >> 1556 J = 1; 1422 1557 1423 void G4ElasticHadrNucleusHE::WriteLine(std::o << 1558 for( M = 0; M <= N; M++ ) 1424 std::vector<G4double>& v) << 1559 { 1425 { << 1560 Fact1 = 1; 1426 std::size_t n = v.size(); << 1561 1427 outfile << n << G4endl; << 1562 if ( ( N > 0 ) && ( N > M ) && M > 0 ) 1428 if(n > 0) { << 1563 { 1429 for(std::size_t i=0; i<n; ++i) { << 1564 J *= G4double(N-M+1)/G4double(M); 1430 outfile << v[i] << " "; << 1565 Fact1 = J; 1431 } << 1566 } 1432 outfile << G4endl; << 1567 SetBinom[N][M] = Fact1; >> 1568 } 1433 } 1569 } >> 1570 return; 1434 } 1571 } 1435 1572 >> 1573 >> 1574 // >> 1575 // 1436 ///////////////////////////////////////////// 1576 /////////////////////////////////////////////////////////// >> 1577 1437 1578