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