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

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

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