| Geant4 Cross Reference |
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 //
28 //
29 // G4 Model: hadron diffraction elastic scattering with 4-momentum balance
30 //
31 // Class Description
32 // Final state production model for hadron-hadron elastic scattering
33 // in the framework of quark-diquark model with springy Pomeron.
34 // Projectiles are proton, neutron, pions, kaons.
35 // Targets are proton (and neutron).
36 // Class Description - End
37 //
38 // 02.05.14 V. Grichine - 1-st implementation
39 // 10.10.14 V. Grichine - change to combine with low mass diffraction
40
41 #ifndef G4hhElastic_h
42 #define G4hhElastic_h 1
43
44 #include "globals.hh"
45 #include <complex>
46 #include "G4Integrator.hh"
47
48 #include "G4HadronElastic.hh"
49 #include "G4HadProjectile.hh"
50 #include "G4Nucleus.hh"
51 #include "G4HadronNucleonXsc.hh"
52
53 #include "G4Exp.hh"
54 #include "G4Log.hh"
55
56
57 class G4ParticleDefinition;
58 class G4PhysicsTable;
59 class G4PhysicsLogVector;
60
61 class G4hhElastic : public G4HadronElastic
62 {
63 public:
64 // PL constructor
65 G4hhElastic();
66 // test constructor
67 G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile,
68 G4double plab );
69 // constructor used for low mass diffraction
70 G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile);
71
72 virtual ~G4hhElastic();
73
74 virtual G4bool IsApplicable(const G4HadProjectile &/*aTrack*/,
75 G4Nucleus & /*targetNucleus*/);
76
77
78 void Initialise();
79
80 void BuildTableT( G4ParticleDefinition* target, G4ParticleDefinition* projectile); // , G4double plab );
81 void BuildTableTest( G4ParticleDefinition* target, G4ParticleDefinition* projectile,
82 G4double plab );
83
84 G4double SampleInvariantT( const G4ParticleDefinition* p, G4double plab, G4int, G4int);
85 G4double SampleBisectionalT( const G4ParticleDefinition* p, G4double plab);
86
87 G4double SampleTest(G4double tMin ); // const G4ParticleDefinition* p, );
88
89 G4double GetTransfer( G4int iMomentum, G4int iTransfer, G4double position );
90
91 private:
92
93
94 G4ParticleDefinition* fTarget;
95 G4ParticleDefinition* fProjectile;
96 G4ParticleDefinition* theProton;
97 G4ParticleDefinition* theNeutron;
98 G4ParticleDefinition* thePionPlus;
99 G4ParticleDefinition* thePionMinus;
100
101
102 G4double lowEnergyRecoilLimit;
103 G4double lowEnergyLimitHE;
104 G4double lowEnergyLimitQ;
105 G4double lowestEnergyLimit;
106 G4double plabLowLimit;
107
108 G4int fEnergyBin;
109 G4int fBinT;
110
111 G4PhysicsLogVector* fEnergyVector;
112 G4PhysicsTable* fTableT;
113 std::vector<G4PhysicsTable*> fBankT;
114
115
116 // Gauss model parameters
117
118
119 G4double fMff2;
120 G4double fMQ;
121 G4double fMq;
122
123 G4double fMassTarg; // ~ A
124 G4double fMassProj; // ~ B
125 G4double fMassSum2;
126 G4double fMassDif2;
127
128
129 G4double fRA; // hadron A
130 G4double fRQ;
131 G4double fRq;
132 G4double fAlpha;
133 G4double fBeta;
134
135 G4double fRB; // hadron B
136 G4double fRG;
137 G4double fRg;
138 G4double fGamma;
139 G4double fDelta;
140
141 G4double fAlphaP;
142 G4double fLambdaFF;
143 G4double fLambda;
144 G4double fEta;
145 G4double fImCof;
146 G4double fCofF2;
147 G4double fCofF3;
148 G4double fRhoReIm;
149 G4double fExpSlope;
150 G4double fSo;
151
152 G4double fSigmaTot;
153 G4double fBq;
154 G4double fBQ;
155 G4double fBqQ;
156 G4double fOptRatio;
157 G4double fSpp;
158 G4double fPcms;
159 G4double fQcof; // q prime when integrate
160
161 public: // Gauss model methods
162
163 void SetParameters();
164 void SetSigmaTot(G4double stot){fSigmaTot = stot;};
165 void SetSpp(G4double spp){fSpp = spp;};
166 G4double GetSpp(){return fSpp;};
167 void SetParametersCMS(G4double plab);
168
169 G4double GetBq(){ return fBq;};
170 G4double GetBQ(){ return fBQ;};
171 G4double GetBqQ(){ return fBqQ;};
172 void SetBq(G4double b){fBq = b;};
173 void SetBQ(G4double b){fBQ = b;};
174 void SetBqQ(G4double b){fBqQ = b;};
175 G4double GetRhoReIm(){ return fRhoReIm;};
176
177 void CalculateBQ(G4double b);
178 void CalculateBqQ13(G4double b);
179 void CalculateBqQ12(G4double b);
180 void CalculateBqQ123(G4double b);
181 void SetRA(G4double rn, G4double pq, G4double pQ);
182 void SetRB(G4double rn, G4double pq, G4double pQ);
183
184 void SetAlphaP(G4double a){fAlphaP = a;};
185 void SetImCof(G4double a){fImCof = a;};
186 G4double GetImCof(){return fImCof;};
187 void SetLambda(G4double L){fLambda = L;};
188 void SetEta(G4double E){fEta = E;};
189 void SetCofF2(G4double f){fCofF2 = f;};
190 void SetCofF3(G4double f){fCofF3 = f;};
191 G4double GetCofF2(){return fCofF2;};
192 G4double GetCofF3(){return fCofF3;};
193
194 G4double GetRA(){ return fRA;};
195 G4double GetRq(){ return fRq;};
196 G4double GetRQ(){ return fRQ;};
197
198 G4double GetRB(){ return fRB;};
199 G4double GetRg(){ return fRg;};
200 G4double GetRG(){ return fRG;};
201
202 // FqQgG stuff
203
204 G4complex Pomeron();
205
206 G4complex Phi13();
207 G4complex Phi14();
208 G4complex Phi23();
209 G4complex Phi24();
210
211 G4complex GetF1qQgG(G4double qp);
212 G4double GetdsdtF1qQgG(G4double s, G4double q);
213 G4complex GetF2qQgG(G4double qp);
214 G4double GetdsdtF12qQgG(G4double s, G4double q);
215 G4complex GetF3qQgG(G4double qp);
216 G4double GetdsdtF123qQgG(G4double q); // sampling ds/dt
217 G4double GetdsdtF13qQG(G4double s, G4double q);
218
219
220 // F123 stuff
221
222 G4complex GetAqq();
223 G4complex GetAQQ();
224 G4complex GetAqQ();
225
226 G4double GetCofS1();
227 G4double GetCofS2();
228 G4double GetCofS3();
229 G4double GetOpticalRatio();
230
231 G4complex GetF1(G4double qp);
232 G4double GetdsdtF1(G4double s, G4double q);
233 G4complex GetF2(G4double qp);
234 G4double GetdsdtF12(G4double s, G4double q);
235 G4complex GetF3(G4double qp);
236 G4double GetdsdtF123(G4double q); // sampling ds/dt
237 G4double GetExpRatioF123(G4double s, G4double q);
238
239 // parameter arrays
240
241 private:
242
243 G4int fInTkin;
244 G4double fOldTkin;
245 static const G4double theNuclNuclData[19][6];
246 static const G4double thePiKaNuclData[8][6];
247 G4HadronNucleonXsc* fHadrNuclXsc;
248 };
249
250 //////////////////////////////////////////////////////////////////////
251 //////////////////////////////////////////////////////////////////////
252 ////////////////////////////////////////////////////////////////////////
253
254
255
256 inline G4bool G4hhElastic::IsApplicable(const G4HadProjectile & projectile,
257 G4Nucleus & nucleus)
258 {
259 if( ( projectile.GetDefinition() == G4Proton::Proton() ||
260 projectile.GetDefinition() == G4Neutron::Neutron() ||
261 projectile.GetDefinition() == G4PionPlus::PionPlus() ||
262 projectile.GetDefinition() == G4PionMinus::PionMinus() ||
263 projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
264 projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
265
266 nucleus.GetZ_asInt() < 2 ) return true;
267 else return false;
268 }
269
270
271 inline void G4hhElastic::SetParameters()
272 {
273 // masses
274
275 fMq = 0.36*CLHEP::GeV; // 0.441*GeV; // 0.36*GeV;
276 fMQ = 0.441*CLHEP::GeV;
277 fMff2 = 0.26*CLHEP::GeV*CLHEP::GeV; // 0.25*GeV*GeV; // 0.5*GeV*GeV;
278
279 fAlpha = 1./3.;
280 fBeta = 1. - fAlpha;
281
282 fGamma = 1./2.; // 1./3.; //
283 fDelta = 1. - fGamma; // 1./2.;
284
285 // radii and exp cof
286
287 fRA = 6.5/CLHEP::GeV; // 7.3/GeV; // 3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
288 fRq = 0.173*fRA; // 2.4/GeV;
289 fRQ = 0.316*fRA; // 1./GeV; // 2./GeV; // 1./GeV;
290 fRB = 6.5/CLHEP::GeV; // 7.3/GeV; // 3.25/GeV; // 7./GeV; // 2./GeV; // 1./GeV;
291 fRg = 0.173*fRA; // 2.4/GeV;
292 fRG = 0.173*fRA; // 1./GeV; // 2./GeV; // 1./GeV;
293
294 fAlphaP = 0.15/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
295 fLambda = 0.25*fRA*fRA; // 0.25
296 fEta = 0.25*fRB*fRB; // 0.25
297 fImCof = 6.5;
298 fCofF2 = 1.;
299 fCofF3 = 1.;
300
301 fBq = 0.02; // 0.21; // 1./3.;
302 fBQ = 1. + fBq - 2*std::sqrt(fBq); // 1 - fBq; // 2./3.;
303 fBqQ = std::sqrt(fBq*fBQ);
304
305 fLambdaFF = 1.5/CLHEP::GeV/CLHEP::GeV; // 0.15/GeV/GeV;
306 fSo = 1.*CLHEP::GeV*CLHEP::GeV;
307 fQcof = 0.009*CLHEP::GeV;
308 fExpSlope = 19.9/CLHEP::GeV/CLHEP::GeV;
309 }
310
311
312 ////////////////////////////////////////////////////////////////////////
313 //
314 // Set target and projectile masses and calculate mass sum and difference squared for Pcms
315
316 inline void G4hhElastic::SetParametersCMS(G4double plab)
317 {
318 G4int i;
319 G4double trMass = 900.*CLHEP::MeV, Tkin;
320 G4double sl, sh, ds, rAl, rAh, drA, rBl, rBh, drB, bql, bqh, dbq, bQl, bQh, dbQ, cIl, cIh, dcI;
321
322 Tkin = std::sqrt(fMassProj*fMassProj + plab*plab) - fMassProj;
323
324 G4DynamicParticle* theDynamicParticle = new G4DynamicParticle(fProjectile,
325 G4ParticleMomentum(0.,0.,1.),
326 Tkin);
327 fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, fTarget );
328
329 delete theDynamicParticle;
330
331 fSpp = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj);
332 fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
333
334 G4double sCMS = std::sqrt(fSpp);
335
336 if( fMassProj > trMass ) // p,n,pb on p
337 {
338 this->SetCofF2(1.);
339 this->SetCofF3(1.);
340 fGamma = 1./3.; // 1./3.; //
341 fDelta = 1. - fGamma; // 1./2.;
342
343 if( sCMS <= theNuclNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
344 {
345 this->SetRA(theNuclNuclData[0][1]/CLHEP::GeV,0.173,0.316);
346 this->SetRB(theNuclNuclData[0][2]/CLHEP::GeV,0.173,0.316);
347
348 this->SetBq(theNuclNuclData[0][3]);
349 this->SetBQ(theNuclNuclData[0][4]);
350 this->SetImCof(theNuclNuclData[0][5]);
351
352 this->SetLambda(0.25*this->GetRA()*this->GetRA());
353 this->SetEta(0.25*this->GetRB()*this->GetRB());
354 }
355 else if( sCMS >= theNuclNuclData[17][0]*CLHEP::GeV ) // high edge, as s=7000 ???
356 {
357 this->SetRA(theNuclNuclData[17][1]/CLHEP::GeV,0.173,0.316);
358 this->SetRB(theNuclNuclData[17][2]/CLHEP::GeV,0.173,0.316);
359
360 this->SetBq(theNuclNuclData[17][3]);
361 this->SetBQ(theNuclNuclData[17][4]);
362 this->SetImCof(theNuclNuclData[17][5]);
363
364 this->SetLambda(0.25*this->GetRA()*this->GetRA());
365 this->SetEta(0.25*this->GetRB()*this->GetRB());
366 }
367 else // in approximation between array points
368 {
369 for( i = 0; i < 19; i++ ) if( sCMS <= theNuclNuclData[i][0]*CLHEP::GeV ) break;
370 if( i == 0 ) i++;
371 if( i == 19 ) i--;
372
373 sl = theNuclNuclData[i-1][0]*CLHEP::GeV;
374 sh = theNuclNuclData[i][0]*CLHEP::GeV;
375 ds = (sCMS - sl)/(sh - sl);
376
377 rAl = theNuclNuclData[i-1][1]/CLHEP::GeV;
378 rAh = theNuclNuclData[i][1]/CLHEP::GeV;
379 drA = rAh - rAl;
380
381 rBl = theNuclNuclData[i-1][2]/CLHEP::GeV;
382 rBh = theNuclNuclData[i][2]/CLHEP::GeV;
383 drB = rBh - rBl;
384
385 bql = theNuclNuclData[i-1][3];
386 bqh = theNuclNuclData[i][3];
387 dbq = bqh - bql;
388
389 bQl = theNuclNuclData[i-1][4];
390 bQh = theNuclNuclData[i][4];
391 dbQ = bQh - bQl;
392
393 cIl = theNuclNuclData[i-1][5];
394 cIh = theNuclNuclData[i][5];
395 dcI = cIh - cIl;
396
397 this->SetRA(rAl+drA*ds,0.173,0.316);
398 this->SetRB(rBl+drB*ds,0.173,0.316);
399
400 this->SetBq(bql+dbq*ds);
401 this->SetBQ(bQl+dbQ*ds);
402 this->SetImCof(cIl+dcI*ds);
403
404 this->SetLambda(0.25*this->GetRA()*this->GetRA());
405 this->SetEta(0.25*this->GetRB()*this->GetRB());
406 }
407 }
408 else // pi, K
409 {
410 this->SetCofF2(1.);
411 this->SetCofF3(-1.);
412 fGamma = 1./2.; // 1./3.; //
413 fDelta = 1. - fGamma; // 1./2.;
414
415 if( sCMS <= thePiKaNuclData[0][0]*CLHEP::GeV ) // low edge, as s=2.76754
416 {
417 this->SetRA(thePiKaNuclData[0][1]/CLHEP::GeV,0.173,0.316);
418 this->SetRB(thePiKaNuclData[0][2]/CLHEP::GeV,0.173,0.173);
419
420 this->SetBq(thePiKaNuclData[0][3]);
421 this->SetBQ(thePiKaNuclData[0][4]);
422 this->SetImCof(thePiKaNuclData[0][5]);
423
424 this->SetLambda(0.25*this->GetRA()*this->GetRA());
425 this->SetEta(this->GetRB()*this->GetRB()/6.);
426 }
427 else if( sCMS >= thePiKaNuclData[7][0]*CLHEP::GeV ) // high edge, as s=7000 ???
428 {
429 this->SetRA(thePiKaNuclData[7][1]/CLHEP::GeV,0.173,0.316);
430 this->SetRB(thePiKaNuclData[7][2]/CLHEP::GeV,0.173,0.173);
431
432 this->SetBq(thePiKaNuclData[7][3]);
433 this->SetBQ(thePiKaNuclData[7][4]);
434 this->SetImCof(thePiKaNuclData[7][5]);
435
436 this->SetLambda(0.25*this->GetRA()*this->GetRA());
437 this->SetEta(this->GetRB()*this->GetRB()/6.);
438 }
439 else // in approximation between array points
440 {
441 for( i = 0; i < 8; i++ ) if( sCMS <= thePiKaNuclData[i][0]*CLHEP::GeV ) break;
442 if( i == 0 ) i++;
443 if( i == 8 ) i--;
444
445 sl = thePiKaNuclData[i-1][0]*CLHEP::GeV;
446 sh = thePiKaNuclData[i][0]*CLHEP::GeV;
447 ds = (sCMS - sl)/(sh - sl);
448
449 rAl = thePiKaNuclData[i-1][1]/CLHEP::GeV;
450 rAh = thePiKaNuclData[i][1]/CLHEP::GeV;
451 drA = rAh - rAl;
452
453 rBl = thePiKaNuclData[i-1][2]/CLHEP::GeV;
454 rBh = thePiKaNuclData[i][2]/CLHEP::GeV;
455 drB = rBh - rBl;
456
457 bql = thePiKaNuclData[i-1][3];
458 bqh = thePiKaNuclData[i][3];
459 dbq = bqh - bql;
460
461 bQl = thePiKaNuclData[i-1][4];
462 bQh = thePiKaNuclData[i][4];
463 dbQ = bQh - bQl;
464
465 cIl = thePiKaNuclData[i-1][5];
466 cIh = thePiKaNuclData[i][5];
467 dcI = cIh - cIl;
468
469 this->SetRA(rAl+drA*ds,0.173,0.316);
470 this->SetRB(rBl+drB*ds,0.173,0.173);
471
472 this->SetBq(bql+dbq*ds);
473 this->SetBQ(bQl+dbQ*ds);
474 this->SetImCof(cIl+dcI*ds);
475
476 this->SetLambda(0.25*this->GetRA()*this->GetRA());
477 this->SetEta(this->GetRB()*this->GetRB()/6.);
478 }
479 }
480 return;
481 }
482
483 /////////////////////////////////////////////////////
484 //
485 // RA for qQ
486
487 inline void G4hhElastic::SetRA(G4double rA, G4double pq, G4double pQ)
488 {
489 fRA = rA;
490 fRq = fRA*pq;
491 fRQ = fRA*pQ;
492 }
493
494 /////////////////////////////////////////////////////
495 //
496 // RB for gG
497
498 inline void G4hhElastic::SetRB(G4double rB, G4double pg, G4double pG)
499 {
500 fRB = rB;
501 fRg = fRB*pg;
502 fRG = fRB*pG;
503 }
504
505 ////////////////////////////////////////////////////
506 //
507 // Returns Pomeron parametrization with Im part modified, *= fImCof
508
509 inline G4complex G4hhElastic::Pomeron()
510 {
511 G4double re, im;
512
513 re = fAlphaP*G4Log(fSpp/fSo);
514 im = -0.5*fAlphaP*fImCof*CLHEP::pi;
515 return G4complex(re,im);
516 }
517
518 //////////////////////////////////////////////////
519
520 inline G4complex G4hhElastic::Phi13()
521 {
522 G4double re = (fRq*fRq + fRg*fRg)/16.;
523 G4complex result(re,0.);
524 result += Pomeron();
525 return result;
526 }
527
528 //////////////////////////////////////////////////
529
530 inline G4complex G4hhElastic::Phi14()
531 {
532 G4double re = (fRq*fRq + fRG*fRG)/16.;
533 G4complex result(re,0.);
534 result += Pomeron();
535 return result;
536 }
537
538 //////////////////////////////////////////////////
539
540 inline G4complex G4hhElastic::Phi23()
541 {
542 G4double re = (fRQ*fRQ + fRg*fRg)/16.;
543 G4complex result(re,0.);
544 result += Pomeron();
545 return result;
546 }
547
548 //////////////////////////////////////////////////
549
550 inline G4complex G4hhElastic::Phi24()
551 {
552 G4double re = (fRQ*fRQ + fRG*fRG)/16.;
553 G4complex result(re,0.);
554 result += Pomeron();
555 return result;
556 }
557
558 /////////////////////////////////////////////////////
559 //
560 // F1, case qQ-gG
561
562 inline G4complex G4hhElastic::GetF1qQgG(G4double t)
563 {
564 G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
565 G4double k = p/CLHEP::hbarc;
566
567 G4complex exp13 = fBq*std::exp(-(Phi13() + fBeta*fBeta*fLambda + fDelta*fDelta*fEta)*t);
568 G4complex exp14 = fBq*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
569 G4complex exp23 = fBQ*std::exp(-(Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta)*t);
570 G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);
571
572 G4complex res = exp13 + exp14 + exp23 + exp24;
573
574 res *= 0.25*k*fSigmaTot/CLHEP::pi;
575 res *= G4complex(0.,1.);
576
577 return res;
578 }
579
580 /////////////////////////////////////////////////////
581 //
582 //
583
584 inline G4double G4hhElastic::GetdsdtF13qQG(G4double spp, G4double t)
585 {
586 fSpp = spp;
587 G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
588 G4double k = p/CLHEP::hbarc;
589
590 G4complex exp14 = fBqQ*std::exp(-(Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta)*t);
591 G4complex exp24 = fBQ*std::exp(-(Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta)*t);
592
593 G4complex F1 = exp14 + exp24;
594
595 F1 *= 0.25*k*fSigmaTot/CLHEP::pi;
596 F1 *= G4complex(0.,1.);
597
598 // 1424
599
600 G4complex z1424 = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
601 z1424 /= Phi14() + Phi24() + fLambda;
602 z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
603
604
605 G4complex exp1424 = std::exp(-z1424*t);
606 exp1424 /= Phi14() + Phi24() + fLambda;
607
608 G4complex F3 = fBqQ*fBQ*exp1424;
609
610
611 F3 *= 0.25*k/CLHEP::pi;
612 F3 *= G4complex(0.,1.);
613 F3 *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
614
615 G4complex F13 = F1 - F3;
616
617 G4double dsdt = CLHEP::pi/p/p;
618 dsdt *= real(F13)*real(F13) + imag(F13)*imag(F13);
619
620 return dsdt;
621 }
622
623 //////////////////////////////////////////////////////////////
624 //
625 // dsigma/dt(s,t) F1qQgG
626
627 inline G4double G4hhElastic::GetdsdtF1qQgG(G4double spp, G4double t)
628 {
629 fSpp = spp;
630 G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
631
632 G4complex F1 = GetF1qQgG(t);
633
634 G4double dsdt = CLHEP::pi/p/p;
635 dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
636 return dsdt;
637 }
638
639 /////////////////////////////////////////////////////
640 //
641 //
642
643 inline G4complex G4hhElastic::GetF2qQgG(G4double t)
644 {
645 G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
646 G4double k = p/CLHEP::hbarc;
647
648 G4complex z1324 = -(Phi24() + fAlpha*fLambda + fGamma*fEta)*(Phi24() + fAlpha*fLambda + fGamma*fEta);
649 z1324 /= Phi13() + Phi24() + fLambda + fEta;
650 z1324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
651
652 G4complex exp1324 = std::exp(-z1324*t);
653 exp1324 /= Phi13() + Phi24() + fLambda + fEta;
654
655 G4complex z1423 = -(Phi23() + fAlpha*fLambda + fDelta*fEta)*(Phi24() + fAlpha*fLambda + fDelta*fEta);;
656 z1423 /= Phi14() + Phi23() + fLambda + fEta;
657 z1423 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;
658
659 G4complex exp1423 = std::exp(-z1423*t);
660 exp1423 /= Phi14() + Phi23() + fLambda + fEta;
661
662 G4complex res = exp1324 + exp1423;
663
664
665 res *= 0.25*k/CLHEP::pi;
666 res *= G4complex(0.,1.);
667 res *= fBq*fBQ*fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc); // or 4. ???
668
669 return res;
670 }
671
672 //////////////////////////////////////////////////////////////
673 //
674 // dsigma/dt(s,t) F12
675
676 inline G4double G4hhElastic::GetdsdtF12qQgG( G4double spp, G4double t)
677 {
678 fSpp = spp;
679 G4double p = std::sqrt((fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
680
681 G4complex F12 = GetF1qQgG(t) - GetF2qQgG(t);
682
683 G4double dsdt = CLHEP::pi/p/p;
684 dsdt *= real(F12)*real(F12) + imag(F12)*imag(F12);
685 return dsdt;
686 }
687
688 /////////////////////////////////////////////////////
689 //
690 //
691
692 inline G4complex G4hhElastic::GetF3qQgG(G4double t)
693 {
694 G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
695 G4double k = p/CLHEP::hbarc;
696
697 // 1314
698
699 G4complex z1314 = -(Phi14() + fGamma*fEta)*(Phi14() + fGamma*fEta);
700 z1314 /= Phi13() + Phi14() + fEta;
701 z1314 += Phi14() + fBeta*fBeta*fLambda + fGamma*fGamma*fEta;
702
703 G4complex exp1314 = std::exp(-z1314*t);
704 exp1314 /= Phi13() + Phi14() + fEta;
705
706 // 2324
707
708 G4complex z2324 = -(Phi24() + fGamma*fEta)*(Phi24() + fGamma*fEta);;
709 z2324 /= Phi24() + Phi23() + fEta;
710 z2324 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
711
712 G4complex exp2324 = std::exp(-z2324*t);
713 exp2324 /= Phi24() + Phi23() + fEta;
714
715 // 1323
716
717 G4complex z1323 = -(Phi23() + fAlpha*fLambda)*(Phi23() + fAlpha*fLambda);
718 z1323 /= Phi13() + Phi23() + fLambda;
719 z1323 += Phi23() + fAlpha*fAlpha*fLambda + fDelta*fDelta*fEta;
720
721 G4complex exp1323 = std::exp(-z1323*t);
722 exp1323 /= Phi13() + Phi23() + fLambda;
723
724 // 1424
725
726 G4complex z1424 = -(Phi24() + fAlpha*fLambda)*(Phi24() + fAlpha*fLambda);
727 z1424 /= Phi14() + Phi24() + fLambda;
728 z1424 += Phi24() + fAlpha*fAlpha*fLambda + fGamma*fGamma*fEta;
729
730 G4complex exp1424 = std::exp(-z1424*t);
731 exp1424 /= Phi14() + Phi24() + fLambda;
732
733 G4complex res = fBq*fBq*exp1314 + fBQ*fBQ*exp2324 + fBq*fBQ*exp1323 + fBq*fBQ*exp1424;
734
735 res *= 0.25*k/CLHEP::pi;
736 res *= G4complex(0.,1.);
737 res *= fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
738
739 return res;
740 }
741
742 //////////////////////////////////////////////////////////////
743 //
744 // dsigma/dt(s,t) F123 sampling ds/dt
745
746 inline G4double G4hhElastic::GetdsdtF123qQgG(G4double t)
747 {
748 G4double p = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp );
749
750 G4complex F123 = GetF1qQgG(t); // - fCofF2*GetF2qQgG(t) - fCofF3*GetF3qQgG(t);
751 F123 -= fCofF2*GetF2qQgG(t);
752 F123 -= fCofF3*GetF3qQgG(t);
753
754 G4double dsdt = CLHEP::pi/p/p;
755 dsdt *= real(F123)*real(F123) + imag(F123)*imag(F123);
756 return dsdt;
757 }
758
759 /////////////////////////////////////////////////////
760 //
761 // Set fBqQ at a given fBQ=b2 according to the optical theorem,qQ-G
762
763 inline void G4hhElastic::CalculateBqQ13(G4double b2)
764 {
765 fBQ = b2;
766
767 G4complex z1424 = G4complex(1./8./CLHEP::pi,0.);
768 z1424 /= Phi14() + Phi24() + fAlpha;
769 G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);
770
771 fBqQ = 1. - fBQ;
772 fBQ /= 1. - fSigmaTot*fBQ*c1424;
773
774 G4cout<<"fSigmaTot*fBQ*c1424 = "<<fSigmaTot*fBQ*c1424<<G4endl;
775
776 G4double ratio = fBqQ + fBQ - fSigmaTot*fBqQ*fBQ*c1424;
777 G4cout<<"ratio = "<<ratio<<G4endl;
778
779 return ;
780 }
781
782 /////////////////////////////////////////////////////
783 //
784 // Set fBQ at a given fBq=b according to the optical theorem, F1-F2
785
786 inline void G4hhElastic::CalculateBqQ12(G4double b1)
787 {
788 fBq = b1;
789
790 G4complex z1324 = G4complex(1./8./CLHEP::pi,0.);
791 z1324 /= Phi13() + Phi24() + fLambda + fEta;
792 G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);
793
794 G4complex z1423 = G4complex(1./8./CLHEP::pi,0.);
795 z1423 /= Phi14() + Phi23() + fLambda + fEta;
796 G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);
797
798 fBQ = 1. - 2.*fBq;
799 fBQ /= 2. - fSigmaTot*fBq*(c1324+1423);
800
801 G4double ratio = 2.*(fBq + fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423);
802 G4cout<<"ratio = "<<ratio<<G4endl;
803
804 return ;
805 }
806
807 /////////////////////////////////////////////////////
808 //
809 // Set fBQ at a given fBq=b according to the optical theorem, F1-F2-F3,
810 // simplified meson-barion case g=G=q
811
812 inline void G4hhElastic::CalculateBqQ123(G4double b1)
813 {
814 fBq = b1;
815
816 G4complex z1324 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
817 z1324 /= Phi13() + Phi24() + fLambda + fEta;
818 G4double c1324 = real(z1324)/(CLHEP::hbarc*CLHEP::hbarc);
819
820 G4complex z1423 = fCofF2*G4complex(1./8./CLHEP::pi,0.);
821 z1423 /= Phi14() + Phi23() + fLambda + fEta;
822 G4double c1423 = real(z1423)/(CLHEP::hbarc*CLHEP::hbarc);
823
824 G4complex z1314 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
825 z1314 /= Phi13() + Phi14() + fEta;
826 G4double c1314 = real(z1314)/(CLHEP::hbarc*CLHEP::hbarc);
827
828 G4complex z2324 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
829 z2324 /= Phi23() + Phi24() + fEta;
830 G4double c2324 = real(z2324)/(CLHEP::hbarc*CLHEP::hbarc);
831
832 G4complex z1323 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
833 z1323 /= Phi13() + Phi23() + fLambda;
834 G4double c1323 = real(z1323)/(CLHEP::hbarc*CLHEP::hbarc);
835
836 G4complex z1424 = fCofF3*G4complex(1./8./CLHEP::pi,0.);
837 z1424 /= Phi14() + Phi24() + fLambda;
838 G4double c1424 = real(z1424)/(CLHEP::hbarc*CLHEP::hbarc);
839
840 G4double A = fSigmaTot*c2324;
841 G4double B = fSigmaTot*fBq*(c1324 + c1423 + c1323 + c1424) - 2.;
842 G4double C = 1. + fSigmaTot*fBq*fBq*c1314 - 2*fBq;
843 G4cout<<"A = "<<A<<"; B = "<<B<<"; C = "<<C<<G4endl;
844 G4cout<<"determinant = "<<B*B-4.*A*C<<G4endl;
845
846 G4double x1 = ( -B - std::sqrt(B*B-4.*A*C) )/2./A;
847 G4double x2 = ( -B + std::sqrt(B*B-4.*A*C) )/2./A;
848 G4cout<<"x1 = "<<x1<<"; x2 = "<<x2<<G4endl;
849
850 if( B*B-4.*A*C < 1.e-6 ) fBQ = std::abs(-B/2./A);
851 else if ( B < 0.) fBQ = std::abs( ( -B - std::sqrt(B*B-4.*A*C) )/2./A);
852 else fBQ = std::abs( ( -B + std::sqrt(B*B-4.*A*C) )/2./A);
853
854 fOptRatio = 2*(fBq+fBQ) - fSigmaTot*fBq*fBQ*(c1324 + c1423 + c1323 + c1424);
855 fOptRatio -= fSigmaTot*fBq*fBq*c1314 + fSigmaTot*c2324*fBQ*fBQ;
856 G4cout<<"BqQ123, fOptRatio = "<<fOptRatio<<G4endl;
857
858 return ;
859 }
860
861
862
863 ///////////////////// F123 stuff hh-elastic, qQ-qQ ///////////////////
864 //////////////////////////////////////////////////////////////
865 /////////////////////////////////////////////////////////////////
866
867 inline G4complex G4hhElastic::GetAqq()
868 {
869 G4double re, im;
870
871 re = fRq*fRq/8. + fAlphaP*G4Log(fSpp/fSo) + 8.*fLambda/9.;
872 im = -0.5*fAlphaP*fImCof*CLHEP::pi;
873 return G4complex(re,im);
874 }
875
876 /////////////////////////////////////////////////////
877 //
878 //
879
880 inline G4complex G4hhElastic::GetAQQ()
881 {
882 G4double re, im;
883
884 re = fRQ*fRQ/8. + fAlphaP*G4Log(fSpp/fSo) + 2.*fLambda/9.;
885 im = -0.5*fAlphaP*fImCof*CLHEP::pi;
886 return G4complex(re,im);
887 }
888
889 /////////////////////////////////////////////////////
890 //
891 //
892
893 inline G4complex G4hhElastic::GetAqQ()
894 {
895 G4complex z = 0.5*( GetAqq() + GetAQQ() );
896 return z;
897 }
898
899 /////////////////////////////////////////////////////
900 //
901 //
902
903 inline G4double G4hhElastic::GetCofS1()
904 {
905 G4complex z = 1./( GetAqQ() + 4.*fLambda/9. );
906 G4double result = real(z);
907 result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
908 result *= fSigmaTot*fCofF2;
909 return result;
910 }
911
912 /////////////////////////////////////////////////////
913 //
914 //
915
916 inline G4double G4hhElastic::GetCofS2()
917 {
918 G4complex z = 1./( GetAqq() + GetAqQ() - 4.*fLambda/9. );
919 G4double result = real(z);
920 result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
921 result *= fSigmaTot*fCofF3;
922 return result;
923 }
924
925 /////////////////////////////////////////////////////
926 //
927 //
928
929 inline G4double G4hhElastic::GetCofS3()
930 {
931 G4complex z = 1./( GetAQQ() + GetAqQ() + 2.*fLambda/9. );
932 G4double result = real(z);
933 result /= 4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc;
934 result *= fSigmaTot*fCofF3;
935 return result;
936 }
937
938 /////////////////////////////////////////////////////
939 //
940 //
941
942 inline G4double G4hhElastic::GetOpticalRatio()
943 {
944 return fOptRatio;
945 // G4double sqrtBqBQ = std::sqrt(fBq*fBQ);
946 // G4double result = fBq + fBQ + 2.*sqrtBqBQ - 1.;
947 // result /= sqrtBqBQ*( GetCofS1()*sqrtBqBQ + GetCofS2()*fBq + GetCofS3()*fBQ );
948 // return result;
949 }
950
951
952
953 /////////////////////////////////////////////////////
954 //
955 // Set fBQ at a given fBq=b according to the optical theorem
956
957 inline void G4hhElastic::CalculateBQ(G4double b1)
958 {
959 fBq = b1;
960 G4double s1 = GetCofS1();
961 G4double s2 = GetCofS2();
962 G4double s3 = GetCofS3();
963 G4double sqrtBq = std::sqrt(fBq);
964
965 // cofs of the fBQ 3rd equation
966
967 G4double a = s3*sqrtBq;
968 G4double b = s1*fBq - 1.;
969 G4double c = (s2*fBq - 2.)*sqrtBq;
970 G4double d = 1. - fBq;
971
972 // cofs of the incomplete 3rd equation
973
974 G4double p = c/a;
975 p -= b*b/a/a/3.;
976 G4double q = d/a;
977 q -= b*c/a/a/3.;
978 q += 2*b*b*b/a/a/a/27.;
979
980 // cofs for the incomplete colutions
981
982 G4double D = p*p*p/3./3./3.;
983 D += q*q/2./2.;
984 G4complex A1 = G4complex(- q/2., std::sqrt(-D) );
985 G4complex A = std::pow(A1,1./3.);
986
987 G4complex B1 = G4complex(- q/2., -std::sqrt(-D) );
988 G4complex B = std::pow(B1,1./3.);
989
990 // roots of the incomplete 3rd equation
991
992 G4complex y1 = A + B;
993 G4complex y2 = -0.5*(A + B) + 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
994 G4complex y3 = -0.5*(A + B) - 0.5*std::sqrt(3.)*(A - B)*G4complex(0.,1.);
995
996 G4complex x1 = y1 - b/a/3.;
997 G4complex x2 = y2 - b/a/3.;
998 G4complex x3 = y3 - b/a/3.;
999
1000 G4cout<<"re_x1 = "<<real(x1)<<"; re_x2 = "<<real(x2)<<"; re_x3 = "<<real(x3)<<G4endl;
1001 G4cout<<"im_x1 = "<<imag(x1)<<"; im_x2 = "<<imag(x2)<<"; im_x3 = "<<imag(x3)<<G4endl;
1002
1003 G4double r1 = real(x1)*real(x1);
1004 G4double r2 = real(x2)*real(x2);
1005 G4double r3 = real(x3)*real(x3);
1006
1007 if( r1 <= 1. && r1 >= 0. ) fBQ = r1;
1008 else if( r2 <= 1. && r2 >= 0. ) fBQ = r2;
1009 else if( r3 <= 1. && r3 >= 0. ) fBQ = r3;
1010 else fBQ = 1.;
1011 // fBQ = real(x3)*real(x3);
1012 G4double sqrtBqBQ = std::sqrt(fBq*fBQ);
1013 fOptRatio = fBq + fBQ + 2.*sqrtBqBQ - 1.;
1014 fOptRatio /= sqrtBqBQ*( GetCofS1()*sqrtBqBQ + GetCofS2()*fBq + GetCofS3()*fBQ );
1015 G4cout<<"F123, fOptRatio = "<<fOptRatio<<G4endl;
1016
1017 return ;
1018 }
1019
1020 /////////////////////////////////////////////////////
1021 //
1022 //
1023
1024 inline G4complex G4hhElastic::GetF1(G4double t)
1025 {
1026 G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1027 G4double k = p/CLHEP::hbarc;
1028 G4complex exp1 = fBq*std::exp(-GetAqq()*t);
1029 G4complex exp2 = fBQ*std::exp(-GetAQQ()*t);
1030 G4complex exp3 = 2.*std::sqrt(fBq*fBQ)*std::exp(-GetAqQ()*t);
1031
1032 G4complex res = exp1 + exp2 + exp3;
1033 res *= 0.25*k*fSigmaTot/CLHEP::pi;
1034 res *= G4complex(0.,1.);
1035 return res;
1036 }
1037
1038 //////////////////////////////////////////////////////////////
1039 //
1040 // dsigma/dt(s,t) F1
1041
1042 inline G4double G4hhElastic::GetdsdtF1(G4double spp, G4double t)
1043 {
1044 fSpp = spp;
1045 G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1046 G4complex F1 = GetF1(t);
1047
1048 G4double dsdt = CLHEP::pi/p/p;
1049 dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
1050 return dsdt;
1051 }
1052
1053 /////////////////////////////////////////////////////
1054 //
1055 //
1056
1057 inline G4complex G4hhElastic::GetF2(G4double t)
1058 {
1059 G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1060 G4double k = p/CLHEP::hbarc;
1061 G4complex z1 = GetAqq()*GetAQQ() - 16.*fLambda*fLambda/81.;
1062 z1 /= 2.*(GetAqQ() + 4.*fLambda/9.);
1063 G4complex exp1 = std::exp(-z1*t);
1064
1065 G4complex z2 = 0.5*( GetAqQ() - 4.*fLambda/9.);
1066
1067 G4complex exp2 = std::exp(-z2*t);
1068
1069 G4complex res = exp1 + exp2;
1070
1071 G4complex z3 = GetAqQ() + 4.*fLambda/9.;
1072
1073 res *= 0.25*k/CLHEP::pi;
1074 res *= G4complex(0.,1.);
1075 res /= z3;
1076 res *= fBq*fBQ*fSigmaTot*fSigmaTot/(8.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
1077
1078 return res;
1079 }
1080
1081 //////////////////////////////////////////////////////////////
1082 //
1083 // dsigma/dt(s,t) F12
1084
1085 inline G4double G4hhElastic::GetdsdtF12(G4double spp, G4double t)
1086 {
1087 fSpp = spp;
1088 G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1089 G4complex F1 = GetF1(t) - GetF2(t);
1090
1091 G4double dsdt = CLHEP::pi/p/p;
1092 dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
1093 return dsdt;
1094 }
1095
1096 /////////////////////////////////////////////////////
1097 //
1098 //
1099
1100 inline G4complex G4hhElastic::GetF3(G4double t)
1101 {
1102 G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1103 G4double k = p/CLHEP::hbarc;
1104 G4complex z1 = GetAqq()*GetAqQ() - 4.*fLambda*fLambda/81.;
1105 z1 /= GetAqq() + GetAqQ() - 4.*fLambda/9.;
1106
1107 G4complex exp1 = std::exp(-z1*t)*fBq/(GetAqq() + GetAqQ() - 4.*fLambda/9.);
1108
1109 G4complex z2 = GetAqQ()*GetAQQ() - 1.*fLambda*fLambda/81.;
1110 z2 /= GetAQQ() + GetAqQ() + 2.*fLambda/9.;
1111
1112 G4complex exp2 = std::exp(-z2*t)*fBQ/(GetAQQ() + GetAqQ() + 2.*fLambda/9.);
1113
1114 G4complex res = exp1 + exp2;
1115
1116
1117 res *= 0.25*k/CLHEP::pi;
1118 res *= G4complex(0.,1.);
1119 res *= std::sqrt(fBq*fBQ)*fSigmaTot*fSigmaTot/(4.*CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
1120
1121 return res;
1122 }
1123
1124 //////////////////////////////////////////////////////////////
1125 //
1126 // dsigma/dt(s,t) F123, sampling ds/dt
1127
1128 inline G4double G4hhElastic::GetdsdtF123(G4double t)
1129 {
1130 G4double p = std::sqrt(0.25*fSpp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1131 G4complex F1 = GetF1(t);
1132 F1 -= fCofF2*GetF2(t);
1133 F1 -= fCofF3*GetF3(t);
1134 G4double dsdt = CLHEP::pi/p/p;
1135 dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
1136 return dsdt;
1137 }
1138
1139 //////////////////////////////////////////////////////////////
1140 //
1141 // dsigma/dt(s,t) F123
1142
1143 inline G4double G4hhElastic::GetExpRatioF123(G4double spp, G4double t)
1144 {
1145 fSpp = spp;
1146 G4double p = std::sqrt(0.25*spp - CLHEP::proton_mass_c2*CLHEP::proton_mass_c2);
1147
1148 // qQ-ds/dt
1149
1150 G4complex F1 = GetF1(t) - fCofF2*GetF2(t) - fCofF3*GetF3(t);
1151
1152 G4double dsdt = CLHEP::pi/p/p;
1153 dsdt *= real(F1)*real(F1) + imag(F1)*imag(F1);
1154
1155 // exponent ds/dt
1156
1157 G4complex F10 = GetF1(0.) - fCofF2*GetF2(0.) - fCofF3*GetF3(0.);
1158
1159 fRhoReIm = real(F10)/imag(F10);
1160
1161 G4double dsdt0 = CLHEP::pi/p/p;
1162 dsdt0 *= real(F10)*real(F10) + imag(F10)*imag(F10);
1163
1164 dsdt0 *= G4Exp(-fExpSlope*t);
1165
1166 G4double ratio = dsdt/dsdt0;
1167
1168 return ratio;
1169 }
1170
1171
1172 //
1173 //
1174 ////////////////////////////////////////////////////////////////////////
1175
1176 #endif
1177
1178
1179
1180
1181
1182
1183