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