Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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