00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "G4ElasticHadrNucleusHE.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "Randomize.hh"
00049 #include "G4ios.hh"
00050 #include "G4ParticleTable.hh"
00051 #include "G4IonTable.hh"
00052 #include "G4Proton.hh"
00053 #include "G4NistManager.hh"
00054
00055 #include <cmath>
00056
00057 using namespace std;
00058
00060
00061
00062
00063 G4ElasticData::G4ElasticData(const G4ParticleDefinition* p,
00064 G4int Z, G4double AWeight, G4double* eGeV)
00065 {
00066 hadr = p;
00067 massGeV = p->GetPDGMass()/GeV;
00068 mass2GeV2= massGeV*massGeV;
00069 AtomicWeight = G4int(AWeight + 0.5);
00070
00071 DefineNucleusParameters(AWeight);
00072
00073 limitQ2 = 35./(R1*R1);
00074
00075 G4double dQ2 = limitQ2/(ONQ2 - 1.);
00076
00077 TableQ2[0] = 0.0;
00078
00079 for(G4int ii = 1; ii < ONQ2; ii++)
00080 {
00081 TableQ2[ii] = TableQ2[ii-1]+dQ2;
00082 }
00083
00084 massA = AWeight*amu_c2/GeV;
00085 massA2 = massA*massA;
00086
00087 for(G4int kk = 0; kk < NENERGY; kk++)
00088 {
00089 dnkE[kk] = 0;
00090 G4double elab = eGeV[kk] + massGeV;
00091 G4double plab2= eGeV[kk]*(eGeV[kk] + 2.0*massGeV);
00092 G4double Q2m = 4.0*plab2*massA2/(mass2GeV2 + massA2 + 2.*massA2*elab);
00093
00094 if(Z == 1 && p == G4Proton::Proton()) Q2m *= 0.5;
00095
00096 maxQ2[kk] = std::min(limitQ2, Q2m);
00097 TableCrossSec[ONQ2*kk] = 0.0;
00098
00099
00100 }
00101 }
00102
00104
00105
00106
00107 void G4ElasticData::DefineNucleusParameters(G4double A)
00108 {
00109 switch (AtomicWeight)
00110 {
00111 case 207:
00112 case 208:
00113 R1 = 20.5;
00114
00115
00116 R2 = 15.74;
00117
00118
00119 Pnucl = 0.4;
00120 Aeff = 0.7;
00121 break;
00122 case 237:
00123 case 238:
00124 R1 = 21.7;
00125 R2 = 16.5;
00126 Pnucl = 0.4;
00127 Aeff = 0.7;
00128 break;
00129 case 90:
00130 case 91:
00131 R1 = 16.5*1.0;
00132 R2 = 11.62;
00133 Pnucl = 0.4;
00134 Aeff = 0.7;
00135 break;
00136 case 58:
00137 case 59:
00138 R1 = 15.0*1.05;
00139 R2 = 9.9;
00140 Pnucl = 0.45;
00141 Aeff = 0.85;
00142 break;
00143 case 48:
00144 case 47:
00145 R1 = 14.0;
00146 R2 = 9.26;
00147 Pnucl = 0.31;
00148 Aeff = 0.75;
00149 break;
00150 case 40:
00151 case 41:
00152 R1 = 13.3;
00153 R2 = 9.26;
00154 Pnucl = 0.31;
00155 Aeff = 0.75;
00156 break;
00157 case 28:
00158 case 29:
00159 R1 = 12.0;
00160 R2 = 7.64;
00161 Pnucl = 0.253;
00162 Aeff = 0.8;
00163 break;
00164 case 16:
00165 R1 = 10.50;
00166 R2 = 5.5;
00167 Pnucl = 0.7;
00168 Aeff = 0.98;
00169 break;
00170 case 12:
00171 R1 = 9.3936;
00172 R2 = 4.63;
00173 Pnucl = 0.7;
00174
00175 Aeff = 1.0;
00176 break;
00177 case 11:
00178 R1 = 9.0;
00179 R2 = 5.42;
00180 Pnucl = 0.19;
00181
00182 Aeff = 0.9;
00183 break;
00184 case 9:
00185 R1 = 9.9;
00186 R2 = 6.5;
00187 Pnucl = 0.690;
00188 Aeff = 0.95;
00189 break;
00190 case 4:
00191 R1 = 5.3;
00192 R2 = 3.7;
00193 Pnucl = 0.4;
00194 Aeff = 0.75;
00195 break;
00196 case 1:
00197 R1 = 4.5;
00198 R2 = 2.3;
00199 Pnucl = 0.177;
00200 Aeff = 0.9;
00201 break;
00202 default:
00203 R1 = 4.45*std::pow(A - 1.,0.309)*0.9;
00204 R2 = 2.3 *std::pow(A, 0.36);
00205
00206 if(A < 100 && A > 3) Pnucl = 0.176 + 0.00275*A;
00207 else Pnucl = 0.4;
00208
00209
00210
00211
00212 if(A >= 100) Aeff = 0.7;
00213 else if(A < 100 && A > 75) Aeff = 1.5 - 0.008*A;
00214 else Aeff = 0.9;
00215 break;
00216 }
00217
00218
00219 }
00220
00222
00223
00224
00225
00226 G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE(const G4String& name)
00227 : G4HadronElastic(name)
00228 {
00229 dQ2 = hMass = hMass2 = hLabMomentum = hLabMomentum2 = MomentumCM = HadrEnergy
00230 = R1 = R2 = Pnucl = Aeff = HadrTot = HadrSlope = HadrReIm = TotP = DDSect2
00231 = DDSect3 = ConstU = FmaxT = Slope1 = Slope2 = Coeff1 = Coeff2 = MaxTR
00232 = Slope0 = Coeff0 = aAIm = aDIm = Dtot11 = 0.0;
00233 NumbN = iHadrCode = iHadron = 0;
00234
00235 verboseLevel = 0;
00236 plabLowLimit = 20.0*MeV;
00237 lowestEnergyLimit = 0.0;
00238
00239
00240 MbToGeV2 = 2.568;
00241 sqMbToGeV = 1.602;
00242 Fm2ToGeV2 = 25.68;
00243 GeV2 = GeV*GeV;
00244 protonM = proton_mass_c2/GeV;
00245 protonM2 = protonM*protonM;
00246
00247 BoundaryP[0]=9.0;BoundaryTG[0]=5.0;BoundaryTL[0]=0.;
00248 BoundaryP[1]=20.0;BoundaryTG[1]=1.5;BoundaryTL[1]=0.;
00249 BoundaryP[2]=5.0; BoundaryTG[2]=1.0;BoundaryTL[2]=1.5;
00250 BoundaryP[3]=8.0; BoundaryTG[3]=3.0;BoundaryTL[3]=0.;
00251 BoundaryP[4]=7.0; BoundaryTG[4]=3.0;BoundaryTL[4]=0.;
00252 BoundaryP[5]=5.0; BoundaryTG[5]=2.0;BoundaryTL[5]=0.;
00253 BoundaryP[6]=5.0; BoundaryTG[6]=1.5;BoundaryTL[6]=3.0;
00254
00255 Binom();
00256
00257 Energy[0] = 0.4;
00258 Energy[1] = 0.6;
00259 Energy[2] = 0.8;
00260 LowEdgeEnergy[0] = 0.0;
00261 LowEdgeEnergy[1] = 0.5;
00262 LowEdgeEnergy[2] = 0.7;
00263 G4double e = 1.0;
00264 G4double f = std::pow(10.,0.1);
00265 for(G4int i=3; i<NENERGY; i++) {
00266 Energy[i] = e;
00267 LowEdgeEnergy[i] = e/f;
00268 e *= f*f;
00269 }
00270 nistManager = G4NistManager::Instance();
00271
00272
00273 G4int ic[NHADRONS] = {211,-211,2112,2212,321,-321,130,310,311,-311,
00274 3122,3222,3112,3212,3312,3322,3334,
00275 -2212,-2112,-3122,-3222,-3112,-3212,-3312,-3322,-3334};
00276
00277 G4int id[NHADRONS] = {2,3,6,0,4,5,4,4,4,5,
00278 0,0,0,0,0,0,0,
00279 1,7,1,1,1,1,1,1,1};
00280
00281 G4int id1[NHADRONS] = {3,4,1,0,5,6,5,5,5,6,
00282 0,0,0,0,0,0,0,
00283 2,2,2,2,2,2,2,2,2};
00284
00285 for(G4int j=0; j<NHADRONS; j++)
00286 {
00287 HadronCode[j] = ic[j];
00288 HadronType[j] = id[j];
00289 HadronType1[j] = id1[j];
00290
00291 for(G4int k = 0; k < 93; k++) { SetOfElasticData[j][k] = 0; }
00292 }
00293 }
00294
00295
00296 void G4ElasticHadrNucleusHE::Description() const
00297 {
00298 char* dirName = getenv("G4PhysListDocDir");
00299 if (dirName) {
00300 std::ofstream outFile;
00301 G4String outFileName = GetModelName() + ".html";
00302 G4String pathName = G4String(dirName) + "/" + outFileName;
00303 outFile.open(pathName);
00304 outFile << "<html>\n";
00305 outFile << "<head>\n";
00306
00307 outFile << "<title>Description of G4ElasticHadrNucleusHE Model</title>\n";
00308 outFile << "</head>\n";
00309 outFile << "<body>\n";
00310
00311 outFile << "G4ElasticHadrNucleusHE is a hadron-nucleus elastic scattering\n"
00312 << "model developed by N. Starkov which uses a Glauber model\n"
00313 << "parameterization to calculate the final state. It is valid\n"
00314 << "for all hadrons with incident energies above 1 GeV.\n";
00315
00316 outFile << "</body>\n";
00317 outFile << "</html>\n";
00318 outFile.close();
00319 }
00320 }
00321
00322
00324
00325
00326
00327 G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusHE()
00328 {
00329 for(G4int j = 0; j < NHADRONS; j++)
00330 {
00331 for(G4int k = 0; k < 93; k++)
00332 {
00333 if(SetOfElasticData[j][k]) delete SetOfElasticData[j][k];
00334 }
00335 }
00336 }
00337
00339
00340
00341
00342 G4double
00343 G4ElasticHadrNucleusHE::SampleInvariantT(const G4ParticleDefinition* p,
00344 G4double inLabMom,
00345 G4int Z, G4int N)
00346 {
00347 G4double plab = inLabMom/GeV;
00348 G4double Q2 = 0;
00349
00350 iHadrCode = p->GetPDGEncoding();
00351
00352 NumbN = N;
00353
00354 if(verboseLevel > 1)
00355 {
00356 G4cout<< " G4ElasticHadrNucleusHE::SampleT: "
00357 << " for " << p->GetParticleName()
00358 << " at Z= " << Z << " A= " << N
00359 << " plab(GeV)= " << plab
00360 << G4endl;
00361 }
00362 G4int idx;
00363
00364 for( idx = 0 ; idx < NHADRONS; idx++)
00365 {
00366 if(iHadrCode == HadronCode[idx]) break;
00367 }
00368
00369
00370
00371 if( idx >= NHADRONS ) return Q2;
00372
00373 iHadron = HadronType[idx];
00374 iHadrCode = HadronCode[idx];
00375
00376 if(Z==1)
00377 {
00378 hMass = p->GetPDGMass()/GeV;
00379 hMass2 = hMass*hMass;
00380
00381 G4double T = sqrt(plab*plab+hMass2)-hMass;
00382
00383 if(T > 0.4) Q2 = HadronProtonQ2(p, plab);
00384
00385 if (verboseLevel>1)
00386 G4cout<<" Proton : Q2 "<<Q2<<G4endl;
00387 }
00388 else
00389 {
00390 G4ElasticData* ElD1 = SetOfElasticData[idx][Z];
00391
00392
00393 if(!ElD1)
00394 {
00395 G4double AWeight = nistManager->GetAtomicMassAmu(Z);
00396 ElD1 = new G4ElasticData(p, Z, AWeight, Energy);
00397 SetOfElasticData[idx][Z] = ElD1;
00398
00399 if(verboseLevel > 1)
00400 {
00401 G4cout<< " G4ElasticHadrNucleusHE::SampleT: new record " << idx
00402 << " for " << p->GetParticleName() << " Z= " << Z
00403 << G4endl;
00404 }
00405 }
00406 hMass = ElD1->massGeV;
00407 hMass2 = ElD1->mass2GeV2;
00408 G4double M = ElD1->massA;
00409 G4double M2 = ElD1->massA2;
00410 G4double plab2 = plab*plab;
00411 G4double Q2max = 4.*plab2*M2/
00412 (hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2));
00413
00414
00415 G4double T = sqrt(plab2+hMass2)-hMass;
00416
00417 if(T > 0.4) Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max);
00418
00419 if(verboseLevel > 1)
00420 G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< " t/tmax= " << Q2/Q2max <<G4endl;
00421 }
00422 return Q2*GeV2;
00423 }
00424
00426
00427
00428
00429 G4double
00430 G4ElasticHadrNucleusHE::SampleT(const G4ParticleDefinition* p,
00431 G4double inLabMom,
00432 G4int Z, G4int N)
00433 {
00434 return SampleInvariantT(p, inLabMom, Z, N);
00435 }
00436
00438
00439
00440
00441 G4double G4ElasticHadrNucleusHE::
00442 HadronNucleusQ2_2(G4ElasticData* pElD, G4int Z,
00443 G4double plab, G4double tmax)
00444 {
00445 G4double LineFq2[ONQ2];
00446
00447 G4double Rand = G4UniformRand();
00448
00449 G4int iNumbQ2 = 0;
00450 G4double Q2 = 0.0;
00451
00452 G4double ptot2 = plab*plab;
00453 G4double ekin = std::sqrt(hMass2 + ptot2) - hMass;
00454
00455 if(verboseLevel > 1)
00456 G4cout<<"Q2_2: ekin plab "<<ekin<<" "<<plab<<" tmax "<<tmax<<G4endl;
00457
00458
00459 G4int NumbOnE;
00460 for( NumbOnE = 0; NumbOnE < NENERGY-1; NumbOnE++ )
00461 {
00462 if( ekin <= LowEdgeEnergy[NumbOnE+1] ) break;
00463 }
00464 G4double* dNumbQ2 = pElD->TableQ2;
00465
00466 G4int index = NumbOnE*ONQ2;
00467
00468
00469 G4double T = Energy[NumbOnE];
00470 hLabMomentum2 = T*(T + 2.*hMass);
00471 G4double Q2max = pElD->maxQ2[NumbOnE];
00472 G4int length = pElD->dnkE[NumbOnE];
00473
00474
00475 if(length == 0)
00476 {
00477 R1 = pElD->R1;
00478 R2 = pElD->R2;
00479 Aeff = pElD->Aeff;
00480 Pnucl = pElD->Pnucl;
00481 hLabMomentum = std::sqrt(hLabMomentum2);
00482
00483 DefineHadronValues(Z);
00484
00485 if(verboseLevel >0)
00486 {
00487 G4cout<<"1 plab T "<<plab<<" "<<T<<" sigTot B ReIm "
00488 <<HadrTot<<" "<<HadrSlope<<" "<<HadrReIm<<G4endl;
00489 G4cout<<" R1 R2 Aeff p "<<R1<<" "<<R2<<" "<<Aeff<<" "
00490 <<Pnucl<<G4endl;
00491 }
00492
00493
00494
00495 if(verboseLevel > 1)
00496 G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE
00497 << " length= " << length
00498 << " Q2max= " << Q2max
00499 << " ekin= " << ekin <<G4endl;
00500
00501 pElD->TableCrossSec[index] = 0;
00502
00503
00504 dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0];
00505
00506 GetHeavyFq2(Z, NumbN, LineFq2);
00507
00508 for(G4int ii=0; ii<ONQ2; ii++)
00509 {
00510
00511
00512
00513
00514 pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1];
00515 }
00516
00517 pElD->dnkE[NumbOnE] = ONQ2;
00518 length = ONQ2;
00519 }
00520
00521 G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]);
00522
00523 for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ )
00524 {
00525 if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break;
00526 }
00527 Q2 = GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand);
00528
00529 if(tmax < Q2max) Q2 *= tmax/Q2max;
00530
00531 if(verboseLevel > 1)
00532 G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<" iNumbQ2= " << iNumbQ2
00533 << " rand= " << Rand << G4endl;
00534
00535 return Q2;
00536 }
00537
00539
00540
00541
00542 G4double G4ElasticHadrNucleusHE::GetQ2_2(G4int kk, G4double * Q,
00543 G4double * F, G4double ranUni)
00544 {
00545 G4double ranQ2;
00546
00547 G4double F1 = F[kk-2];
00548 G4double F2 = F[kk-1];
00549 G4double F3 = F[kk];
00550 G4double X1 = Q[kk-2];
00551 G4double X2 = Q[kk-1];
00552 G4double X3 = Q[kk];
00553
00554 if(verboseLevel > 2)
00555 G4cout << "GetQ2_2 kk= " << kk << " X2= " << X2 << " X3= " << X3
00556 << " F2= " << F2 << " F3= " << F3 << " Rndm= " << ranUni << G4endl;
00557
00558 if(kk == 1 || kk == 0)
00559 {
00560 F1 = F[0];
00561 F2 = F[1];
00562 F3 = F[2];
00563 X1 = Q[0];
00564 X2 = Q[1];
00565 X3 = Q[2];
00566 }
00567
00568 G4double F12 = F1*F1;
00569 G4double F22 = F2*F2;
00570 G4double F32 = F3*F3;
00571
00572 G4double D0 = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3;
00573
00574 if(verboseLevel > 2)
00575 G4cout << " X1= " << X1 << " F1= " << F1 << " D0= "
00576 << D0 << G4endl;
00577
00578 if(std::abs(D0) < 0.00000001)
00579 {
00580 ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2);
00581 }
00582 else
00583 {
00584 G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1;
00585 G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*F22;
00586 G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22
00587 -X1*F2*F32-X2*F3*F12-X3*F1*F22;
00588 ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0;
00589 }
00590 return ranQ2;
00591 }
00592
00594
00595
00596 G4double G4ElasticHadrNucleusHE::GetHeavyFq2(G4int Z, G4int Nucleus, G4double* LineF)
00597 {
00598 G4int ii, jj, aSimp;
00599 G4double curQ2, curSec;
00600 G4double curSum = 0.0;
00601 G4double totSum = 0.0;
00602
00603 G4double ddQ2 = dQ2/20;
00604 G4double Q2l = 0;
00605
00606 LineF[0] = 0;
00607 for(ii = 1; ii<ONQ2; ii++)
00608 {
00609 curSum = 0;
00610 aSimp = 4;
00611
00612 for(jj = 0; jj<20; jj++)
00613 {
00614 curQ2 = Q2l+jj*ddQ2;
00615
00616 curSec = HadrNucDifferCrSec(Z, Nucleus, curQ2);
00617 curSum += curSec*aSimp;
00618
00619 if(aSimp > 3) aSimp = 2;
00620 else aSimp = 4;
00621
00622 if(jj == 0 && verboseLevel>2)
00623 G4cout<<" Q2 "<<curQ2<<" AIm "<<aAIm<<" DIm "<<aDIm
00624 <<" Diff "<<curSec<<" totSum "<<totSum<<G4endl;
00625 }
00626
00627 Q2l += dQ2;
00628 curSum *= ddQ2/2.3;
00629 totSum += curSum;
00630
00631 LineF[ii] = totSum;
00632
00633 if (verboseLevel>2)
00634 G4cout<<" GetHeavy: Q2 dQ2 totSum "<<Q2l<<" "<<dQ2<<" "<<totSum
00635 <<" curSec "
00636 <<curSec<<" totSum "<< totSum<<" DTot "
00637 <<curSum<<G4endl;
00638 }
00639 return totSum;
00640 }
00641
00643
00644
00645
00646 G4double G4ElasticHadrNucleusHE::GetLightFq2(G4int Z, G4int Nucleus,
00647 G4double Q2)
00648 {
00649
00650 if(Z == 1)
00651 {
00652 G4double SqrQ2 = std::sqrt(Q2);
00653 G4double valueConstU = 2.*(hMass2 + protonM2) - Q2;
00654
00655 G4double y = (1.-Coeff1-Coeff0)/HadrSlope*(1.-std::exp(-HadrSlope*Q2))
00656 + Coeff0*(1.-std::exp(-Slope0*Q2))
00657 + Coeff2/Slope2*std::exp(Slope2*valueConstU)*(std::exp(Slope2*Q2)-1.)
00658 + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
00659
00660 return y;
00661 }
00662
00663
00664
00665 G4double prec = Nucleus > 208 ? 1.0e-7 : 1.0e-6;
00666
00667 G4double Stot = HadrTot*MbToGeV2;
00668 G4double Bhad = HadrSlope;
00669 G4double Asq = 1+HadrReIm*HadrReIm;
00670 G4double Rho2 = std::sqrt(Asq);
00671
00672
00673 G4cout<<" Fq2 Before for i Tot B Im "<<HadrTot<<" "<<HadrSlope<<" "
00674 <<HadrReIm<<G4endl;
00675
00676 if(verboseLevel > 1) {
00677 G4cout << "GetFq2: Stot= " << Stot << " Bhad= " << Bhad
00678 <<" Im "<<HadrReIm
00679 << " Asq= " << Asq << G4endl;
00680 G4cout << "R1= " << R1 << " R2= " << R2 << " Pnucl= " << Pnucl <<G4endl;
00681 }
00682 G4double R12 = R1*R1;
00683 G4double R22 = R2*R2;
00684 G4double R12B = R12+2*Bhad;
00685 G4double R22B = R22+2*Bhad;
00686
00687 G4double Norm = (R12*R1-Pnucl*R22*R2);
00688
00689 G4double R13 = R12*R1/R12B;
00690 G4double R23 = Pnucl*R22*R2/R22B;
00691 G4double Unucl = Stot/twopi/Norm*R13;
00692 G4double UnucRho2 = -Unucl*Rho2;
00693
00694 G4double FiH = std::asin(HadrReIm/Rho2);
00695 G4double NN2 = R23/R13;
00696
00697 if(verboseLevel > 2)
00698 G4cout << "UnucRho2= " << UnucRho2 << " FiH= " << FiH << " NN2= " << NN2
00699 << " Norm= " << Norm << G4endl;
00700
00701 G4double dddd;
00702
00703 G4double Prod0 = 0;
00704 G4double N1 = -1.0;
00705
00706 G4double exp1;
00707
00708 G4double Prod3 ;
00709 G4double exp2 ;
00710 G4double N4, N5, N2, Prod1, Prod2;
00711 G4int i1, i2, j1, j2;
00712
00713 for(i1 = 1; i1<= Nucleus; i1++)
00714 {
00715 N1 = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2;
00716 Prod1 = 0;
00717
00718 N2 = -1;
00719
00720 for(i2 = 1; i2<=Nucleus; i2++)
00721 {
00722 N2 = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2;
00723 Prod2 = 0;
00724 N5 = -1/NN2;
00725 for(j2=0; j2<= i2; j2++)
00726 {
00727 Prod3 = 0;
00728 exp2 = 1/(j2/R22B+(i2-j2)/R12B);
00729 N5 = -N5*NN2;
00730 N4 = -1/NN2;
00731 for(j1=0; j1<=i1; j1++)
00732 {
00733 exp1 = 1/(j1/R22B+(i1-j1)/R12B);
00734 dddd = exp1+exp2;
00735 N4 = -N4*NN2;
00736 Prod3 = Prod3+N4*exp1*exp2*
00737 (1-std::exp(-Q2*dddd/4))/dddd*4*SetBinom[i1][j1];
00738 }
00739 Prod2 = Prod2 +Prod3*N5*SetBinom[i2][j2];
00740 }
00741 Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2));
00742
00743 if (std::fabs(Prod2*N2/Prod1)<prec) break;
00744 }
00745 Prod0 = Prod0 + Prod1*N1;
00746 if(std::fabs(N1*Prod1/Prod0) < prec) break;
00747
00748 }
00749
00750 Prod0 *= 0.25*pi/MbToGeV2;
00751
00752 if(verboseLevel>1)
00753 G4cout << "GetLightFq2 Z= " << Z << " A= " << Nucleus
00754 <<" Q2= " << Q2 << " Res= " << Prod0 << G4endl;
00755 return Prod0;
00756 }
00757
00758 G4double G4ElasticHadrNucleusHE::
00759 HadrNucDifferCrSec(G4int Z, G4int , G4double aQ2)
00760 {
00761
00762
00763
00764 G4int NWeight = int( nistManager->GetAtomicMassAmu(Z) + 0.5 );
00765
00766 G4double theQ2 = aQ2;
00767
00768
00769 if(NWeight == 1)
00770 {
00771 G4double SqrQ2 = std::sqrt(aQ2);
00772 G4double valueConstU = hMass2 + protonM2-2*protonM*HadrEnergy - aQ2;
00773
00774 G4double MaxT = 4*MomentumCM*MomentumCM;
00775
00776 BoundaryTL[0] = MaxT;
00777 BoundaryTL[1] = MaxT;
00778 BoundaryTL[3] = MaxT;
00779 BoundaryTL[4] = MaxT;
00780 BoundaryTL[5] = MaxT;
00781
00782 G4double dSigPodT;
00783
00784 dSigPodT = HadrTot*HadrTot*(1+HadrReIm*HadrReIm)*
00785 (
00786 Coeff1*std::exp(-Slope1*SqrQ2)+
00787 Coeff2*std::exp( Slope2*(valueConstU)+aQ2)+
00788 (1-Coeff1-Coeff0)*std::exp(-HadrSlope*aQ2)+
00789 +Coeff0*std::exp(-Slope0*aQ2)
00790
00791 )/16/3.1416*2.568;
00792
00793 return dSigPodT;
00794 }
00795
00796 G4double Stot = HadrTot*MbToGeV2;
00797 G4double Bhad = HadrSlope;
00798 G4double Asq = 1+HadrReIm*HadrReIm;
00799 G4double Rho2 = std::sqrt(Asq);
00800 G4double Pnuclp = 0.001;
00801 Pnuclp = Pnucl;
00802 G4double R12 = R1*R1;
00803 G4double R22 = R2*R2;
00804 G4double R12B = R12+2*Bhad;
00805 G4double R22B = R22+2*Bhad;
00806 G4double R12Ap = R12+20;
00807 G4double R22Ap = R22+20;
00808 G4double R13Ap = R12*R1/R12Ap;
00809 G4double R23Ap = R22*R2/R22Ap*Pnuclp;
00810 G4double R23dR13 = R23Ap/R13Ap;
00811 G4double R12Apd = 2/R12Ap;
00812 G4double R22Apd = 2/R22Ap;
00813 G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd);
00814
00815 G4double DDSec1p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R1/4));
00816 G4double DDSec2p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/
00817 std::sqrt((R12+R22)/2)/4));
00818 G4double DDSec3p = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R2/4));
00819
00820 G4double Norm = (R12*R1-Pnucl*R22*R2)*Aeff;
00821 G4double Normp = (R12*R1-Pnuclp*R22*R2)*Aeff;
00822 G4double R13 = R12*R1/R12B;
00823 G4double R23 = Pnucl*R22*R2/R22B;
00824 G4double Unucl = Stot/twopi/Norm*R13;
00825 G4double UnuclScr = Stot/twopi/Normp*R13Ap;
00826 G4double SinFi = HadrReIm/Rho2;
00827 G4double FiH = std::asin(SinFi);
00828 G4double N = -1;
00829 G4double N2 = R23/R13;
00830
00831 G4double ImElasticAmpl0 = 0;
00832 G4double ReElasticAmpl0 = 0;
00833
00834 G4double exp1;
00835 G4double N4;
00836 G4double Prod1, Tot1=0, medTot, DTot1, DmedTot;
00837 G4int i;
00838
00839 for( i=1; i<=NWeight; i++)
00840 {
00841 N = -N*Unucl*(NWeight-i+1)/i*Rho2;
00842 N4 = 1;
00843 Prod1 = std::exp(-theQ2/i*R12B/4)/i*R12B;
00844 medTot = R12B/i;
00845
00846 for(G4int l=1; l<=i; l++)
00847 {
00848 exp1 = l/R22B+(i-l)/R12B;
00849 N4 = -N4*(i-l+1)/l*N2;
00850 Prod1 = Prod1+N4/exp1*std::exp(-theQ2/exp1/4);
00851 medTot = medTot+N4/exp1;
00852 }
00853
00854 ReElasticAmpl0 = ReElasticAmpl0+Prod1*N*std::sin(FiH*i);
00855 ImElasticAmpl0 = ImElasticAmpl0+Prod1*N*std::cos(FiH*i);
00856 Tot1 = Tot1+medTot*N*std::cos(FiH*i);
00857 if(std::fabs(Prod1*N/ImElasticAmpl0) < 0.000001) break;
00858 }
00859
00860 ImElasticAmpl0 = ImElasticAmpl0*pi/2.568;
00861 ReElasticAmpl0 = ReElasticAmpl0*pi/2.568;
00862 Tot1 = Tot1*twopi/2.568;
00863
00864 G4double C1 = R13Ap*R13Ap/2*DDSec1p;
00865 G4double C2 = 2*R23Ap*R13Ap/2*DDSec2p;
00866 G4double C3 = R23Ap*R23Ap/2*DDSec3p;
00867
00868 G4double N1p = 1;
00869
00870 G4double Din1 = 0.5;
00871
00872 Din1 = 0.5*(C1*std::exp(-theQ2/8*R12Ap)/2*R12Ap-
00873 C2/R12ApdR22Ap*std::exp(-theQ2/4/R12ApdR22Ap)+
00874 C3*R22Ap/2*std::exp(-theQ2/8*R22Ap));
00875
00876 DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2);
00877
00878 G4double exp1p;
00879 G4double exp2p;
00880 G4double exp3p;
00881 G4double N2p;
00882 G4double Din2, BinCoeff;
00883
00884 BinCoeff = 1;
00885
00886 for( i = 1; i<= NWeight-2; i++)
00887 {
00888 N1p = -N1p*UnuclScr*(NWeight-i-1)/i*Rho2;
00889 N2p = 1;
00890 Din2 = 0;
00891 DmedTot = 0;
00892 for(G4int l = 0; l<=i; l++)
00893 {
00894 if(l == 0) BinCoeff = 1;
00895 else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l;
00896
00897 exp1 = l/R22B+(i-l)/R12B;
00898 exp1p = exp1+R12Apd;
00899 exp2p = exp1+R12ApdR22Ap;
00900 exp3p = exp1+R22Apd;
00901
00902 Din2 = Din2 + N2p*BinCoeff*
00903 (C1/exp1p*std::exp(-theQ2/4/exp1p)-
00904 C2/exp2p*std::exp(-theQ2/4/exp2p)+
00905 C3/exp3p*std::exp(-theQ2/4/exp3p));
00906
00907 DmedTot = DmedTot + N2p*BinCoeff*
00908 (C1/exp1p-C2/exp2p+C3/exp3p);
00909
00910 N2p = -N2p*R23dR13;
00911 }
00912
00913 Din1 = Din1+Din2*N1p/(i+2)/(i+1)*std::cos(FiH*i);
00914 DTot1 = DTot1+DmedTot*N1p/(i+2)/(i+1)*std::cos(FiH*i);
00915
00916 if(std::fabs(Din2*N1p/Din1) < 0.000001) break;
00917 }
00918
00919 Din1 = -Din1*NWeight*(NWeight-1)
00920 /2/pi/Normp/2/pi/Normp*16*pi*pi;
00921
00922 DTot1 = DTot1*NWeight*(NWeight-1)
00923 /2/pi/Normp/2/pi/Normp*16*pi*pi;
00924
00925 DTot1 *= 5;
00926
00927
00928
00929
00930 G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
00931 (ImElasticAmpl0+Din1)*
00932 (ImElasticAmpl0+Din1))*2/4/pi;
00933
00934 Tot1 = Tot1-DTot1;
00935
00936 Dtot11 = DTot1;
00937 aAIm = ImElasticAmpl0;
00938 aDIm = Din1;
00939
00940 return DiffCrSec2*1.0;
00941 }
00942
00943
00945
00946
00947
00948 void G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
00949 {
00950 HadrEnergy = std::sqrt(hMass2 + hLabMomentum2);
00951
00952 G4double sHadr = 2.*HadrEnergy*protonM+protonM2+hMass2;
00953 G4double sqrS = std::sqrt(sHadr);
00954 G4double Ecm = 0.5*(sHadr-hMass2+protonM2)/sqrS;
00955 MomentumCM = std::sqrt(Ecm*Ecm-protonM2);
00956
00957 if(verboseLevel>2)
00958 G4cout << "GetHadrVall.: Z= " << Z << " iHadr= " << iHadron
00959 << " E(GeV)= " << HadrEnergy << " sqrS= " << sqrS
00960 << " plab= " << hLabMomentum
00961 <<" E - m "<<HadrEnergy - hMass<< G4endl;
00962
00963 G4double TotN = 0.0;
00964 G4double logE = std::log(HadrEnergy);
00965 G4double logS = std::log(sHadr);
00966 TotP = 0.0;
00967
00968 switch (iHadron)
00969 {
00970 case 0:
00971 case 6:
00972
00973 if(hLabMomentum > 10)
00974 TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165);
00975
00976 else
00977 {
00978
00979
00981
00982
00983 if( hLabMomentum > 1.4 )
00984 TotN = 33.3+15.2*(hLabMomentum2-1.35)/
00985 (std::pow(hLabMomentum,2.37)+0.95);
00986
00987 else if(hLabMomentum > 0.8)
00988 {
00989 G4double A0 = logE + 0.0513;
00990 TotN = 33.0 + 25.5*A0*A0;
00991 }
00992 else
00993 {
00994 G4double A0 = logE - 0.2634;
00995 TotN = 33.0 + 30.*A0*A0*A0*A0;
00996 }
00997
00998
00999 {
01000 if(hLabMomentum >= 1.05)
01001 {
01002 TotP = 39.0+75.*(hLabMomentum-1.2)/
01003 (hLabMomentum2*hLabMomentum+0.15);
01004 }
01005
01006 else if(hLabMomentum >= 0.7)
01007 {
01008 G4double A0 = logE + 0.3147;
01009 TotP = 23.0 + 40.*A0*A0;
01010 }
01011 else
01012 {
01013 TotP = 23.+50.*std::pow(std::log(0.73/hLabMomentum),3.5);
01014 }
01015 }
01016 }
01017
01018
01019 HadrTot = 0.5*(TotP+TotN);
01020
01021
01022 if(hLabMomentum >= 2.) HadrSlope = 5.44 + 0.88*logS;
01023
01024 else if(hLabMomentum >= 0.5) HadrSlope = 3.73*hLabMomentum-0.37;
01025
01026 else HadrSlope = 1.5;
01027
01028
01029 if(hLabMomentum >= 1.2)
01030 HadrReIm = 0.13*(logS - 5.8579332)*std::pow(sHadr,-0.18);
01031
01032 else if(hLabMomentum >= 0.6)
01033 HadrReIm = -75.5*(std::pow(hLabMomentum,0.25)-0.95)/
01034 (std::pow(3*hLabMomentum,2.2)+1);
01035
01036 else
01037 HadrReIm = 15.5*hLabMomentum/(27*hLabMomentum2*hLabMomentum+2);
01038
01039 DDSect2 = 2.2;
01040 DDSect3 = 0.6;
01041
01042 if( iHadrCode == 3122)
01043 {
01044 HadrTot *= 0.88;
01045 HadrSlope *=0.85;
01046 }
01047
01048 else if( iHadrCode == 3222)
01049 {
01050 HadrTot *=0.81;
01051 HadrSlope *=0.85;
01052 }
01053
01054 else if(iHadrCode == 3112 || iHadrCode == 3212 )
01055 {
01056 HadrTot *=0.88;
01057 HadrSlope *=0.85;
01058 }
01059
01060 else if( iHadrCode == 3312 || iHadrCode == 3322 )
01061 {
01062 HadrTot *=0.77;
01063 HadrSlope *=0.75;
01064 }
01065
01066 else if( iHadrCode == 3334)
01067 {
01068 HadrTot *=0.78;
01069 HadrSlope *=0.7;
01070 }
01071
01072 break;
01073
01074 case 1:
01075 case 7:
01076
01077 HadrTot = 5.2+5.2*logE + 123.2/sqrS;
01078 HadrSlope = 8.32+0.57*logS;
01079
01080 if( HadrEnergy < 1000 )
01081 HadrReIm = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8);
01082 else
01083 HadrReIm = 0.6*(logS - 5.8579332)*std::pow(sHadr,-0.25);
01084
01085 DDSect2 = 11;
01086 DDSect3 = 3;
01087
01088 if( iHadrCode == -3122)
01089 {
01090 HadrTot *= 0.88;
01091 HadrSlope *=0.85;
01092 }
01093
01094 else if( iHadrCode == -3222)
01095 {
01096 HadrTot *=0.81;
01097 HadrSlope *=0.85;
01098 }
01099
01100 else if(iHadrCode == -3112 || iHadrCode == -3212 )
01101 {
01102 HadrTot *=0.88;
01103 HadrSlope *=0.85;
01104 }
01105
01106 else if( iHadrCode == -3312 || iHadrCode == -3322 )
01107 {
01108 HadrTot *=0.77;
01109 HadrSlope *=0.75;
01110 }
01111
01112 else if( iHadrCode == -3334)
01113 {
01114 HadrTot *=0.78;
01115 HadrSlope *=0.7;
01116 }
01117
01118 break;
01119
01120 case 2:
01121 case 3:
01122
01123 if(hLabMomentum >= 3.5)
01124 TotP = 10.6+2.*logE + 25.*std::pow(HadrEnergy,-0.43);
01125
01126 else if(hLabMomentum >= 1.15)
01127 {
01128 G4double x = (hLabMomentum - 2.55)/0.55;
01129 G4double y = (hLabMomentum - 1.47)/0.225;
01130 TotP = 3.2*std::exp(-x*x) + 12.*std::exp(-y*y) + 27.5;
01131 }
01132
01133 else if(hLabMomentum >= 0.4)
01134 {
01135 TotP = 88*(logE+0.2877)*(logE+0.2877)+14.0;
01136 }
01137
01138 else
01139 {
01140 G4double x = (hLabMomentum - 0.29)/0.085;
01141 TotP = 20. + 180.*std::exp(-x*x);
01142 }
01143
01144
01145 if(hLabMomentum >= 3.0 )
01146 TotN = 10.6 + 2.*logE + 30.*std::pow(HadrEnergy,-0.43);
01147
01148 else if(hLabMomentum >= 1.3)
01149 {
01150 G4double x = (hLabMomentum - 2.1)/0.4;
01151 G4double y = (hLabMomentum - 1.4)/0.12;
01152 TotN = 36.1+0.079 - 4.313*logE + 3.*std::exp(-x*x) +
01153 1.5*std::exp(-y*y);
01154 }
01155 else if(hLabMomentum >= 0.65)
01156 {
01157 G4double x = (hLabMomentum - 0.72)/0.06;
01158 G4double y = (hLabMomentum - 1.015)/0.075;
01159 TotN = 36.1 + 10.*std::exp(-x*x) + 24*std::exp(-y*y);
01160 }
01161 else if(hLabMomentum >= 0.37)
01162 {
01163 G4double x = std::log(hLabMomentum/0.48);
01164 TotN = 26. + 110.*x*x;
01165 }
01166 else
01167 {
01168 G4double x = (hLabMomentum - 0.29)/0.07;
01169 TotN = 28.0 + 40.*std::exp(-x*x);
01170 }
01171 HadrTot = (TotP+TotN)/2;
01172
01173 HadrSlope = 7.28+0.245*logS;
01174 HadrReIm = 0.2*(logS - 4.6051702)*std::pow(sHadr,-0.15);
01175
01176 DDSect2 = 0.7;
01177 DDSect3 = 0.27;
01178
01179 break;
01180
01181 case 4:
01182
01183 HadrTot = 10.6+1.8*logE + 9.0*std::pow(HadrEnergy,-0.55);
01184 if(HadrEnergy>100) HadrSlope = 15.0;
01185 else HadrSlope = 1.0+1.76*logS - 2.84/sqrS;
01186
01187 HadrReIm = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1);
01188 DDSect2 = 0.7;
01189 DDSect3 = 0.21;
01190 break;
01191
01192 case 5:
01193
01194 HadrTot = 10+1.8*logE + 25./sqrS;
01195 HadrSlope = 6.98+0.127*logS;
01196 HadrReIm = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1);
01197 DDSect2 = 0.7;
01198 DDSect3 = 0.27;
01199 break;
01200 }
01201
01202 if(verboseLevel>2)
01203 G4cout << "HadrTot= " << HadrTot << " HadrSlope= " << HadrSlope
01204 << " HadrReIm= " << HadrReIm << " DDSect2= " << DDSect2
01205 << " DDSect3= " << DDSect3 << G4endl;
01206
01207 if(Z != 1) return;
01208
01209
01210
01211 Coeff0 = Coeff1 = Coeff2 = 0.0;
01212 Slope0 = Slope1 = 1.0;
01213 Slope2 = 5.0;
01214
01215
01216 static const G4double EnP0[6]={1.5,3.0,5.0,9.0,14.0,19.0};
01217 static const G4double C0P0[6]={0.15,0.02,0.06,0.08,0.0003,0.0002};
01218 static const G4double C1P0[6]={0.05,0.02,0.03,0.025,0.0,0.0};
01219 static const G4double B0P0[6]={1.5,2.5,3.0,4.5,1.4,1.25};
01220 static const G4double B1P0[6]={5.0,1.0,3.5,4.0,4.8,4.8};
01221
01222
01223 static const G4double EnN[5]={1.5,5.0,10.0,14.0,20.0};
01224 static const G4double C0N[5]={0.0,0.0,0.02,0.02,0.01};
01225 static const G4double C1N[5]={0.06,0.008,0.0015,0.001,0.0003};
01226 static const G4double B0N[5]={1.5,2.5,3.8,3.8,3.5};
01227 static const G4double B1N[5]={1.5,2.2,3.6,4.5,4.8};
01228
01229
01230 static const G4double EnP[2]={1.5,4.0};
01231 static const G4double C0P[2]={0.001,0.0005};
01232 static const G4double C1P[2]={0.003,0.001};
01233 static const G4double B0P[2]={2.5,4.5};
01234 static const G4double B1P[2]={1.0,4.0};
01235
01236
01237 static const G4double EnPP[4]={1.0,2.0,3.0,4.0};
01238 static const G4double C0PP[4]={0.0,0.0,0.0,0.0};
01239 static const G4double C1PP[4]={0.15,0.08,0.02,0.01};
01240 static const G4double B0PP[4]={1.5,2.8,3.8,3.8};
01241 static const G4double B1PP[4]={0.8,1.6,3.6,4.6};
01242
01243
01244 static const G4double EnPPN[4]={1.0,2.0,3.0,4.0};
01245 static const G4double C0PPN[4]={0.0,0.0,0.0,0.0};
01246 static const G4double C1PPN[4]={0.0,0.0,0.0,0.0};
01247 static const G4double B0PPN[4]={1.5,2.8,3.8,3.8};
01248 static const G4double B1PPN[4]={0.8,1.6,3.6,4.6};
01249
01250
01251 static const G4double EnK[4]={1.4,2.33,3.0,5.0};
01252 static const G4double C0K[4]={0.0,0.0,0.0,0.0};
01253 static const G4double C1K[4]={0.01,0.007,0.005,0.003};
01254 static const G4double B0K[4]={1.5,2.0,3.8,3.8};
01255 static const G4double B1K[4]={1.6,1.6,1.6,1.6};
01256
01257
01258 static const G4double EnKM[2]={1.4,4.0};
01259 static const G4double C0KM[2]={0.006,0.002};
01260 static const G4double C1KM[2]={0.00,0.00};
01261 static const G4double B0KM[2]={2.5,3.5};
01262 static const G4double B1KM[2]={1.6,1.6};
01263
01264 switch(iHadron)
01265 {
01266 case 0 :
01267
01268 if(hLabMomentum <BoundaryP[0])
01269 InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P0);
01270
01271 Coeff2 = 0.8/hLabMomentum/hLabMomentum;
01272 break;
01273
01274 case 6 :
01275
01276 if(hLabMomentum < BoundaryP[1])
01277 InterpolateHN(5,EnN,C0N,C1N,B0N,B1N);
01278
01279 Coeff2 = 0.8/hLabMomentum/hLabMomentum;
01280 break;
01281
01282 case 1 :
01283 case 7 :
01284 if(hLabMomentum < BoundaryP[2])
01285 InterpolateHN(2,EnP,C0P,C1P,B0P,B1P);
01286 break;
01287
01288 case 2 :
01289
01290 if(hLabMomentum < BoundaryP[3])
01291 InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1PP);
01292
01293 Coeff2 = 0.02/hLabMomentum;
01294 break;
01295
01296 case 3 :
01297
01298 if(hLabMomentum < BoundaryP[4])
01299 InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN,B1PPN);
01300
01301 Coeff2 = 0.02/hLabMomentum;
01302 break;
01303
01304 case 4 :
01305
01306 if(hLabMomentum < BoundaryP[5])
01307 InterpolateHN(4,EnK,C0K,C1K,B0K,B1K);
01308
01309 if(hLabMomentum < 1) Coeff2 = 0.34;
01310 else Coeff2 = 0.34/hLabMomentum2/hLabMomentum;
01311 break;
01312
01313 case 5 :
01314 if(hLabMomentum < BoundaryP[6])
01315 InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1KM);
01316
01317 if(hLabMomentum < 1) Coeff2 = 0.01;
01318 else Coeff2 = 0.01/hLabMomentum2/hLabMomentum;
01319 break;
01320 }
01321
01322 if(verboseLevel > 2)
01323 G4cout<<" HadrVal : Plasb "<<hLabMomentum
01324 <<" iHadron "<<iHadron<<" HadrTot "<<HadrTot<<G4endl;
01325 }
01326
01327
01328 void G4ElasticHadrNucleusHE::
01329 GetKinematics(const G4ParticleDefinition * aHadron,
01330 G4double MomentumH)
01331 {
01332 if (verboseLevel>1)
01333 G4cout<<"1 GetKin.: HadronName MomentumH "
01334 <<aHadron->GetParticleName()<<" "<<MomentumH<<G4endl;
01335
01336 DefineHadronValues(1);
01337
01338 G4double Sh = 2.0*protonM*HadrEnergy+protonM2+hMass2;
01339
01340 ConstU = 2*protonM2+2*hMass2-Sh;
01341
01342 G4double MaxT = 4*MomentumCM*MomentumCM;
01343
01344 BoundaryTL[0] = MaxT;
01345 BoundaryTL[1] = MaxT;
01346 BoundaryTL[3] = MaxT;
01347 BoundaryTL[4] = MaxT;
01348 BoundaryTL[5] = MaxT;
01349
01350 G4int NumberH=0;
01351
01352 while(iHadrCode!=HadronCode[NumberH]) NumberH++;
01353
01354 NumberH = HadronType1[NumberH];
01355
01356 if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH];
01357 else MaxTR = BoundaryTG[NumberH];
01358
01359 if (verboseLevel>1)
01360 G4cout<<"3 GetKin. : NumberH "<<NumberH
01361 <<" Bound.P[NumberH] "<<BoundaryP[NumberH]
01362 <<" Bound.TL[NumberH] "<<BoundaryTL[NumberH]
01363 <<" Bound.TG[NumberH] "<<BoundaryTG[NumberH]
01364 <<" MaxT MaxTR "<<MaxT<<" "<<MaxTR<<G4endl;
01365
01366
01367 }
01368
01369 G4double G4ElasticHadrNucleusHE::GetFt(G4double Q2)
01370 {
01371 G4double Fdistr=0;
01372 G4double SqrQ2 = std::sqrt(Q2);
01373
01374 Fdistr = (1-Coeff1-Coeff0)
01375 /HadrSlope*(1-std::exp(-HadrSlope*Q2))
01376 + Coeff0*(1-std::exp(-Slope0*Q2))
01377 + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1)
01378 + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
01379
01380 if (verboseLevel>1)
01381 G4cout<<"Old: Coeff0 Coeff1 Coeff2 "<<Coeff0<<" "
01382 <<Coeff1<<" "<<Coeff2<<" Slope Slope0 Slope1 Slope2 "
01383 <<HadrSlope<<" "<<Slope0<<" "<<Slope1<<" "<<Slope2
01384 <<" Fdistr "<<Fdistr<<G4endl;
01385 return Fdistr;
01386 }
01387
01388 G4double G4ElasticHadrNucleusHE::GetQ2(G4double Ran)
01389 {
01390 G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta;
01391 G4double Q2=0;
01392
01393 FmaxT = GetFt(MaxTR);
01394 delta = GetDistrFun(DDD0)-Ran;
01395
01396 while(std::fabs(delta) > 0.0001)
01397 {
01398 if(delta>0)
01399 {
01400 DDD2 = DDD0;
01401 DDD0 = (DDD0+DDD1)*0.5;
01402 }
01403 else if(delta<0)
01404 {
01405 DDD1 = DDD0;
01406 DDD0 = (DDD0+DDD2)*0.5;
01407 }
01408 delta = GetDistrFun(DDD0)-Ran;
01409 }
01410
01411 Q2 = DDD0;
01412
01413 return Q2;
01414 }
01415
01416 G4double G4ElasticHadrNucleusHE::
01417 HadronProtonQ2(const G4ParticleDefinition * p,
01418 G4double inLabMom)
01419 {
01420
01421 hMass = p->GetPDGMass()/GeV;
01422 hMass2 = hMass*hMass;
01423 hLabMomentum = inLabMom;
01424 hLabMomentum2 = hLabMomentum*hLabMomentum;
01425 HadrEnergy = sqrt(hLabMomentum2+hMass2);
01426
01427 G4double Rand = G4UniformRand();
01428
01429 GetKinematics(p, inLabMom);
01430
01431 G4double Q2 = GetQ2(Rand);
01432
01433 return Q2;
01434 }
01435
01436
01437
01439
01440
01441
01442 void G4ElasticHadrNucleusHE::Binom()
01443 {
01444 G4int N, M;
01445 G4double Fact1, J;
01446
01447 for(N = 0; N < 240; N++)
01448 {
01449 J = 1;
01450
01451 for( M = 0; M <= N; M++ )
01452 {
01453 Fact1 = 1;
01454
01455 if ( ( N > 0 ) && ( N > M ) && M > 0 )
01456 {
01457 J *= G4double(N-M+1)/G4double(M);
01458 Fact1 = J;
01459 }
01460 SetBinom[N][M] = Fact1;
01461 }
01462 }
01463 return;
01464 }
01465
01466
01467
01468
01470