source: trunk/source/processes/hadronic/models/coherent_elastic/src/G4ElasticHadrNucleusHE.cc@ 1200

Last change on this file since 1200 was 1196, checked in by garnier, 16 years ago

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