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