Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
FemtoWeightGeneratorKisiel.cxx
1/*
2 * FemtoWeightGeneratorKisiel.cxx
3 *
4 * Created on: 30 sie 2021
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9#include "FemtoWeightGeneratorKisiel.h"
10
11#include "FemtoPair.h"
12#include "Package.h"
13#include "Parameter.h"
14
15#include <TMath.h>
16
17#include <cmath>
18#include <iostream>
19
20#define COULOMBSTEPS 170
21
22namespace Hal {
23 FemtoWeightGeneratorKisiel::FemtoWeightGeneratorKisiel() :
24 FemtoWeightGenerator(),
25 fPionac(0),
26 fOneoveracsq(0),
27 fTwopioverac(0),
28 fCoulqscpart(0),
29 fTwospin(0),
30 fWritegrps(0),
31 fPcount(0),
32 fEuler(0),
33 fF0(0),
34 fD0(0),
35 fRStarOutS(0),
36 fRStarSideS(0),
37 fRStarLongS(0),
38 fRStarS(0) {}
39
40 void FemtoWeightGeneratorKisiel::InitializeGamow() {
41 fTwospin = 0;
42
43 fEuler = 0.577215665;
44 fF0 = 7.77 / 0.197327;
45 fD0 = 2.77 / 0.197327;
46
47 fD0s.re = 0.0;
48 fD0s.im = 0.0;
49 fD0t.re = 0.0;
50 fD0t.im = 0.0;
51
52 fF0s.re = 0.0;
53 fF0s.im = 0.0;
54 fF0t.re = 0.0;
55 fF0t.im = 0.0;
56
57 static unsigned short int isospin = 0;
58 switch (fPairType) {
59 case Femto::EPairType::kPionPlusPionPlus: { // 0
60 fPionac = 387.5 / 0.197327;
61 } break;
62 case Femto::EPairType::kKaonPlusKaonPlus: { // 1
63 fPionac = 109.55 / 0.197327;
64 } break;
65 case Femto::EPairType::kPionPlusKaonPlus: { // 2
66 fPionac = 248.519 / 0.197327;
67 fWritegrps = 0;
68 } break;
69 case Femto::EPairType::kPionPlusKaonMinus: { // 3
70 fPionac = -248.519 / 0.197327;
71 fWritegrps = 0;
72 } break;
73 case Femto::EPairType::kPionPlusProton: { // 4
74 fPionac = 222.564 / 0.197327;
75 fWritegrps = 0;
76 } break;
77 case Femto::EPairType::kPionPlusAntiproton: { // 5
78 fPionac = -222.564 / 0.197327;
79 fWritegrps = 0;
80 } break;
81 case Femto::EPairType::kProtonSigmaPlus: { // 6
82 fPionac = 51.5553 / 0.197327;
83 fTwospin = 1;
84 // Rijken, Nagels, Yamamoto ESC08 model PTPS No. 185, 2010
85 fF0s.re = 3.85 / 0.197327;
86 fD0s.re = 3.40 / 0.197327;
87 fF0t.re = -0.62 / 0.197327;
88 fD0t.re = -2.13 / 0.197327;
89 } break;
90 case Femto::EPairType::kProtonAntiSigmaPlus: { // 7
91 fPionac = -51.5553 / 0.197327;
92
93 // from antiproton - lambda
94 fF0s.re = -1.20 / 0.197327;
95 fF0s.im = 0.37 / 0.197327;
96 fF0t.re = -1.20 / 0.197327;
97 fF0t.im = 0.37 / 0.197327;
98 } break;
99 case Femto::EPairType::kProtonLambda: { // 8
100 // ESC08
101 fF0s.re = 2.70 / 0.197327;
102 fD0s.re = 2.97 / 0.197327;
103 fF0t.re = 1.65 / 0.197327;
104 fD0t.re = 3.63 / 0.197327;
105 } break;
106 case Femto::EPairType::kProtonAntiLambda: { // 9
107 fPionac = 0.0;
108
109 fF0s.re = -1.20 / 0.197327;
110 fF0s.im = 0.37 / 0.197327;
111 fF0t.re = -1.20 / 0.197327;
112 fF0t.im = 0.37 / 0.197327;
113 } break;
114 case Femto::EPairType::kSigmaPlusSigmaPlus: { // 10
115 fPionac = 45.4709 / 0.197327;
116 fTwospin = 1;
117 } break;
118 case Femto::EPairType::kSigmaPlusAntiSigmaPlus: { // 11
119 fPionac = -45.4709 / 0.197327;
120 } break;
121 case Femto::EPairType::kProtonProton: { // 12
122 fPionac = 57.63975274 / 0.197327;
123 fTwospin = 1;
124
125 // d0s.re = 2.77 / 0.197327;
126 // f0s.re = 7.77 / 0.197327;
127 fD0t.re = 1.7 / 0.197327;
128 fF0t.re = -5.4 / 0.197327;
129
130 // ESC08
131 fF0s.re = 7.771 / 0.197327;
132 fD0s.re = 2.754 / 0.197327;
133 // f0t.re = ? / 0.197327; no data from ESC08
134 // d0t.re = ? / 0.197327;
135 } break;
136 case Femto::EPairType::kKaonPlusProton: { // 13
137 fPionac = 83.59432006 / 0.197327;
138 } break;
139 case Femto::EPairType::kKaonPlusAntiproton: { // 14
140 fPionac = -83.59432006 / 0.197327;
141 } break;
142 case Femto::EPairType::kPionPlusPionMinus: { // 15
143 fPionac = -387.5 / 0.197327;
144 } break;
145 /* case Femto::EPairType::kPionPlusPionPlus: {//16 ?? same
146 pionac = -387.5 / 0.197327;
147 } break;*/
148 case Femto::EPairType::kLambdaLambda: { // 17
149 fPionac = 0;
150
151 // ESC08
152 fF0s.re = 0.88 / 0.197327;
153 fD0s.re = 4.34 / 0.197327;
154 fF0t.re = 0.88 / 0.197327;
155 fD0t.re = 4.34 / 0.197327;
156 } break;
157 case Femto::EPairType::kProtonXiZero: { // 18
158 fPionac = 0;
159
160 // ESC08
161 if (isospin == 0) {
162 fF0s.re = 0.0;
163 fD0s.re = 0.0;
164 fF0t.re = -6.90 / 0.197327;
165 fD0t.re = 1.18 / 0.197327;
166 isospin = 1;
167 } else if (isospin == 1) {
168 fF0s.re = -0.58 / 0.197327;
169 fD0s.re = -2.71 / 0.197327;
170 fF0t.re = -3.49 / 0.197327;
171 fD0t.re = 0.60 / 0.197327;
172 isospin = 0;
173 }
174 } break;
175 case Femto::EPairType::kNeutronXiMinus: { // 19
176 fPionac = 0;
177
178 // ESC08
179 if (isospin == 0) {
180 fF0s.re = 0.0;
181 fD0s.re = 0.0;
182 fF0t.re = -6.90 / 0.197327;
183 fD0t.re = 1.18 / 0.197327;
184 isospin = 1;
185 } else if (isospin == 1) {
186 fF0s.re = -0.58 / 0.197327;
187 fD0s.re = -2.71 / 0.197327;
188 fF0t.re = -3.49 / 0.197327;
189 fD0t.re = 0.60 / 0.197327;
190 isospin = 0;
191 }
192 } break;
193 case Femto::EPairType::kProtonXIMInus: { // 20
194 fPionac = 49.2788901 / 0.197327;
195
196 // ESC08
197 if (isospin == 0) {
198 fF0s.re = 0.0;
199 fD0s.re = 0.0;
200 fF0t.re = -6.90 / 0.197327;
201 fD0t.re = 1.18 / 0.197327;
202 isospin = 1;
203 } else if (isospin == 1) {
204 fF0s.re = -0.58 / 0.197327;
205 fD0s.re = -2.71 / 0.197327;
206 fF0t.re = -3.49 / 0.197327;
207 fD0t.re = 0.60 / 0.197327;
208 isospin = 0;
209 }
210 } break;
211 case Femto::EPairType::kNeutronXiZero: { // 21
212 fPionac = 0;
213
214 // ESC08
215 if (isospin == 0) {
216 fF0s.re = 0.0;
217 fD0s.re = 0.0;
218 fF0t.re = -6.90 / 0.197327;
219 fD0t.re = 1.18 / 0.197327;
220 isospin = 1;
221 } else if (isospin == 1) {
222 fF0s.re = -0.58 / 0.197327;
223 fD0s.re = -2.71 / 0.197327;
224 fF0t.re = -3.49 / 0.197327;
225 fD0t.re = 0.60 / 0.197327;
226 isospin = 0;
227 }
228 } break;
229 case Femto::EPairType::kNeutronProton: { // 22
230 fPionac = 0;
231
232 // ESC08
233 fF0s.re = 23.735 / 0.197327;
234 fD0s.re = 2.694 / 0.197327;
235 fF0t.re = -5.428 / 0.197327;
236 fD0t.re = 1.753 / 0.197327;
237 } break;
238 case Femto::EPairType::kNeutronNeutron: { // 23
239 fPionac = 0;
240
241 // ESC08
242 fF0s.re = 16.51 / 0.197327;
243 fD0s.re = 2.85 / 0.197327;
244 } break;
245 case Femto::EPairType::kProtonSigmaZero: { // 24
246 fPionac = 0.0;
247
248 // ESC08
249 fF0s.re = 2.70 / 0.197327;
250 fD0s.re = 2.97 / 0.197327;
251 fF0t.re = 1.65 / 0.197327;
252 fD0t.re = 3.63 / 0.197327;
253 } break;
254 case Femto::EPairType::kSigmaZeroSigmaZero: { // 25
255 fPionac = 0;
256
257 // ESC08
258 fF0s.re = 0.88 / 0.197327;
259 fD0s.re = 4.34 / 0.197327;
260 fF0t.re = 0.88 / 0.197327;
261 fD0t.re = 4.34 / 0.197327;
262 } break;
263 case Femto::EPairType::kLambdaSigmaZero: { // 26
264 fPionac = 0;
265 // ESC08
266 fF0s.re = 0.88 / 0.197327;
267 fD0s.re = 4.34 / 0.197327;
268 fF0t.re = 0.88 / 0.197327;
269 fD0t.re = 4.34 / 0.197327;
270 } break;
271 default: fPionac = 0; break;
272 }
273 fOneoveracsq = 1.0 / (fPionac * fPionac);
274 fTwopioverac = 2.0 * TMath::Pi() / fPionac;
275 double tpaoverk;
276
277 for (int iter = 0; iter < 2000; iter++) {
278 tpaoverk = fTwopioverac / (iter * 0.0002 + 0.0001);
279 // gamov[iter] = tpaoverk * 1.0 / (exp(tpaoverk) - 1);
280 }
281 }
282
283 void FemtoWeightGeneratorKisiel::PairKinematics(FemtoPair* pair) {
284 Double_t tPx = pair->TruePx1() + pair->TruePx2();
285 Double_t tPy = pair->TruePy1() + pair->TruePy2();
286 Double_t tPz = pair->TruePz1() + pair->TruePz2();
287 Double_t tE1 = pair->TrueE1();
288 Double_t tE2 = pair->TrueE2();
289 Double_t tE = tE1 + tE2;
290 Double_t tPt = tPx * tPx + tPy * tPy;
291 Double_t tMt = tE * tE - tPz * tPz; // mCVK;
292 Double_t tM = sqrt(tMt - tPt);
293 tMt = sqrt(tMt);
294 tPt = sqrt(tPt);
295 // Double_t tBetat = tPt/tMt;
296
297 // Boost to LCMS
298 Double_t tBeta = tPz / tE;
299 Double_t tGamma = tE / tMt;
300 fKStarLong = tGamma * (pair->TruePz1() - tBeta * tE1);
301 Double_t tE1L = tGamma * (tE1 - tBeta * pair->TruePz1());
302
303 fKStarOut = (pair->TruePx1() * tPx + pair->TruePy1() * tPy) / tPt;
304 fKStarSide = (-pair->TruePx1() * tPy + pair->TruePy1() * tPx) / tPt;
305
306 fKStarOut = tMt / tM * (fKStarOut - tPt / tMt * tE1L);
307
308
309 Double_t tDX = pair->GetX1() - pair->GetX2();
310 Double_t tDY = pair->GetY1() - pair->GetY2();
311 Double_t tRLong = pair->GetZ1() - pair->GetZ2();
312 Double_t tDTime = pair->GetT1() - pair->GetT2();
313
314 Double_t tROut = (tDX * tPx + tDY * tPy) / tPt;
315 Double_t tRSide = (-tDX * tPy + tDY * tPx) / tPt;
316
317 fRStarSide = tRSide;
318 fRStarSideS = fRStarSide / 0.197327;
319
320 fRStarLong = tGamma * (tRLong - tBeta * tDTime);
321 Double_t tDTimePairLCMS = tGamma * (tDTime - tBeta * tRLong);
322
323 fRStarLongS = fRStarLong / 0.197327;
324 // 1/.1973
325 tBeta = tPt / tMt;
326 tGamma = tMt / tM;
327 fRStarOut = tGamma * (tROut - tBeta * tDTimePairLCMS);
328 fRStarOutS = fRStarOut / 0.197327;
329 // Double_t tDTimePairCMS = tGamma*(tDTimePairLCMS - tBeta* tROut);
330 fRStar = ::sqrt(fRStarOut * fRStarOut + fRStarSide * fRStarSide + fRStarLong * fRStarLong);
331 fKStar = ::sqrt(fKStarOut * fKStarOut + fKStarSide * fKStarSide + fKStarLong * fKStarLong);
332 if (fKStarOut < 0) fKStar = -fKStar;
333 fRStarS = fRStar / 0.197327;
334 }
335
336 Double_t FemtoWeightGeneratorKisiel::GenerateWeight(FemtoPair* pair) {
337 double coulombweight;
338 fPairType = Femto::PidToPairType(pair->GetPdg1(), pair->GetPdg2());
339
340 InitializeGamow();
341 // std::cout << "Czy sa identyczne? " << IsIdentical(localpairtype) << " z pary: "<< localpairtype << std::endl;
342
343 if (Femto::IsPairIdentical(fPairType) == kFALSE) {
344 if (fabs(fPionac) > 0.1) {
345
346 if (fabs(fD0s.re) > 0.00001)
347 coulombweight = GetCoulombStrong();
348 else
349 coulombweight = GetCoulomb();
350 } else {
351 coulombweight = GetStrong();
352 if (std::isnan(coulombweight) == 1) return 1;
353 }
354 } else {
355 if (fPairType == Femto::EPairType::kPionZeroPionZero) { return GetQuantum(); }
356 if (fabs(fPionac) > 0.1) {
357 if (fPairType != Femto::EPairType::kProtonProton)
358 coulombweight = GetQuantumCoulomb();
359 else
360 coulombweight = GetQuantumCoulombStrong();
361 } else {
362 coulombweight = GetQuantumStrong();
363 }
364 }
365
366
367 return coulombweight;
368 }
369
370 double FemtoWeightGeneratorKisiel::Gamow(double arg) const {
371 // if (arg<0.0001)
372 // return 0.0;
373 // if (arg > 0.4)
374 long double eta = fTwopioverac / arg;
375 return (eta) *1.0 / (exp(eta) - 1.0);
376 // int bin = arg / 0.0002;
377 // double b1 = bin*0.0002 + 0.0001;
378 // double b2 = bin*0.0002 + 0.0003;
379 // return ((gamov[bin+1] * (arg - b1) + gamov[bin] * (b2 - arg)) / 0.0002);
380 }
381
382 double FemtoWeightGeneratorKisiel::GetQuantumCoulombStrong() {
383 if (fRStarS < 0.0000000001) return 1.0;
384
385 if (fRStarS < 1.0 / 0.197327) return GetQuantumCoulomb();
386
387 double tKstRst = fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS;
388 long double kstar = fabs(fKStar);
389 long double rho = fRStarS * kstar;
390
391 int ccase = 0;
392 static int pcount = 0;
393 int wavesign = 1;
394
395 // Classical limit - if distance is larger than Coulomb radius,
396 // the interaction does not matter
397 // if (fabs(fRStarS) > fabs(pionac)) return (1.0 + wavesign*cos(2*tKstRst));
398
399 // Classical limit - in the case of large k* we go to
400 // classical coulomb interaction
401 long double testp = rho * (1.0 + tKstRst / (rho));
402 long double testm = rho * (1.0 - tKstRst / (rho));
403
404 dcomplex ffplus, ffminus;
405 if ((testp > 15.0) && (testm > 15.0)) {
406 double asym;
407 asym = (1.0 - 1.0 / (fRStarS * (1.0 - tKstRst / rho) * fPionac * kstar * kstar)) / Gamow(kstar);
408 // std::cout << "as1 " << asym << std::endl;
409 asym = sqrt(asym);
410 if (asym < 1.0)
411 ffminus.re = 1.0 + (asym - 1.0) * 2.0;
412 else
413 ffminus.re = 1.0 + (asym - 1.0) / 2.0;
414 ffminus.im = sqrt(asym * asym - ffminus.re * ffminus.re);
415
416 asym = (1.0 - 1.0 / (fRStarS * (1.0 + tKstRst / rho) * fPionac * kstar * kstar)) / Gamow(kstar);
417 // std::cout << "as2 " << asym << std::endl;
418 asym = sqrt(asym);
419 if (asym < 1.0)
420 ffplus.re = 1.0 + (asym - 1.0) * 2.0;
421 else
422 ffplus.re = 1.0 + (asym - 1.0) / 2.0;
423 ffplus.im = sqrt(asym * asym - ffplus.re * ffplus.re);
424
425 }
426
427 // Check for the classical limit in both functions separately
428 else if (((testp < 15.0) && (testm < 15.0))) // ||
429 {
430 // Calculate the F function
431 GetFFdouble(ffplus, ffminus);
432 ccase = 1;
433 } else if (testp < 15.0) {
434 double asym;
435 GetFFsingle(ffplus, 1);
436 GetFFsingle(ffminus, -1);
437 if ((fabs(ffminus.re) > 2.0) || fabs(ffminus.im) > 2.0) {
438 asym = (1.0 - 1.0 / (fRStarS * (1.0 - tKstRst / (rho) *fPionac * kstar * kstar))) / Gamow(kstar);
439 asym = sqrt(asym);
440 if (asym < 1.0)
441 ffminus.re = 1.0 + (asym - 1.0) * 2.0;
442 else
443 ffminus.re = 1.0 + (asym - 1.0) / 2.0;
444 ffminus.im = sqrt(asym * asym - ffminus.re * ffminus.re);
445 }
446 ccase = 2;
447 } else {
448 double asym;
449 GetFFsingle(ffminus, -1);
450 GetFFsingle(ffplus, 1);
451 if ((fabs(ffplus.re) > 2.0) || fabs(ffplus.im) > 2.0) {
452 asym = (1.0 - 1.0 / (fRStarS * (1.0 + tKstRst / (rho) *fPionac * kstar * kstar))) / Gamow(kstar);
453 asym = sqrt(asym);
454 if (asym < 1.0)
455 ffplus.re = 1.0 + (asym - 1.0) * 2.0;
456 else
457 ffplus.re = 1.0 + (asym - 1.0) / 2.0;
458 ffplus.im = sqrt(asym * asym - ffplus.re * ffplus.re);
459 }
460 ccase = 3;
461 }
462
463 long double eta = 1.0 / (kstar * fPionac);
464 long double hfun = GetH(eta);
465 dcomplex gtilde = GetG(eta, rho, hfun);
466 dcomplex gtilor = mult(gtilde, 1.0 / fRStarS);
467
468 dcomplex fcouls, fcoult;
469 Getfc(kstar, eta, hfun, fcouls, fcoult);
470
471 dcomplex fgs = mult(gtilor, fcouls);
472 long double fgmods = modl2(fgs);
473
474 dcomplex fgt = mult(gtilor, fcoult);
475 long double fgmodt = modl2(fgt);
476
477 dcomplex expikr;
478 expikr.re = cos(tKstRst);
479 expikr.im = -sin(tKstRst);
480 dcomplex expikrc = conj(expikr);
481 dcomplex ffplusc = conj(ffplus);
482 dcomplex ffminusc = conj(ffminus);
483
484 dcomplex expikr2 = mult(expikr, expikr);
485 dcomplex expikrc2 = conj(expikr2);
486 dcomplex sterm = mult(expikr2, mult(ffplus, ffminusc));
487 dcomplex tterm = mult(expikrc2, mult(ffminus, ffplusc));
488
489 dcomplex epfpc = mult(expikr, ffplus);
490 dcomplex emfmc = mult(expikrc, ffminus);
491
492 long double fcgefhs = (fgs.re * emfmc.re + fgs.im * emfmc.im);
493 long double fcgefgs = (fgs.re * epfpc.re + fgs.im * epfpc.im);
494
495 long double fcgefht = (fgt.re * emfmc.re + fgt.im * emfmc.im);
496 long double fcgefgt = (fgt.re * epfpc.re + fgt.im * epfpc.im);
497
498 long double smult = 1 + wavesign;
499
500 if (!std::isfinite(ffminus.re))
501 std::cout << "FFMinus Re not a number ! " << testp << " " << testm << " " << ccase << std::endl;
502
503 if (!std::isfinite(ffminus.im))
504 std::cout << "FFMinus Im not a number !" << testp << " " << testm << " " << ccase << std::endl;
505
506 if ((ffplus.re > 2.0) || (ffplus.re < -2.0))
507 std::cout << std::endl << "FFplus Re wild !" << ffplus.re << " case " << ccase << " " << testp << " " << testm << std::endl;
508
509 if ((ffplus.im > 2.0) || (ffplus.im < -2.0))
510 std::cout << "FFplus Im wild !" << ffplus.im << " case " << ccase << " " << testp << " " << testm << std::endl;
511
512 if ((ffminus.re > 2.0) || (ffminus.re < -2.0))
513 std::cout << std::endl << "FFminus Re wild !" << ffminus.re << " case " << ccase << std::endl;
514
515 if ((ffminus.im > 2.0) || (ffminus.im < -2.0))
516 std::cout << "FFminus Im wild !" << ffminus.im << " case " << ccase << std::endl;
517
518 // coulqscpart = 0.5 * Gamow(kstar) * (modl2(ffplus) + modl2(ffminus));
519
520 // return (0.5 * Gamow(kstar) *
521 // (modl2(ffplus) + wavesign*sterm.re + wavesign*tterm.re + modl2(ffminus)));
522 if (fTwospin == 1) {
523 wavesign = 1;
524 smult = 2;
525 long double singlet = (0.5 * Gamow(kstar)
526 * (2.0 * fgmods * smult + modl2(ffplus) + modl2(ffminus) + wavesign * sterm.re + wavesign * tterm.re
527 + smult * 2 * (fcgefhs + fcgefgs)));
528 wavesign = -1;
529 smult = 0;
530 long double triplet = (0.5 * Gamow(kstar)
531 * (2.0 * fgmodt * smult + modl2(ffplus) + modl2(ffminus) + wavesign * sterm.re + wavesign * tterm.re
532 + smult * 2 * (fcgefht + fcgefgt)));
533 // scmp = singlet;
534 // tcmp = triplet;
535
536 // ccmp = 0.5 * Gamow(kstar) * (modl2(ffplus) + modl2(ffminus));
537 // gcmp = fgmod;
538
539 return (0.25 * singlet + 0.75 * triplet);
540 // return triplet;
541 } else
542 return (0.5 * Gamow(kstar)
543 * (2.0 * fgmods * smult + modl2(ffplus) + modl2(ffminus) + wavesign * sterm.re + wavesign * tterm.re
544 + smult * 2 * (fcgefhs + fcgefgs)));
545 }
546
547 double FemtoWeightGeneratorKisiel::GetCoulombStrong() {
548 if (fRStarS < 0.0000000001) return 1.0;
549
550 if (fRStarS < 1.0 / 0.197327) return GetCoulomb();
551
552 double tKstRst = fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS;
553 long double kstar = fabs(fKStar);
554 long double rho = fRStarS * kstar;
555
556 // Classical limit - in the case of large k* we go to
557 // classical coulomb interaction
558 long double testp = rho * (1.0 + tKstRst / (rho));
559
560 dcomplex ffplus;
561 if (testp > 15.0) {
562 double asym;
563 asym = (1.0 - 1.0 / (fRStarS * (1.0 + tKstRst / rho) * fPionac * kstar * kstar)) / Gamow(kstar);
564 asym = sqrt(asym);
565 if (asym < 1.0)
566 ffplus.re = 1.0 + (asym - 1.0) * 2.0;
567 else
568 ffplus.re = 1.0 + (asym - 1.0) / 2.0;
569 ffplus.im = sqrt(asym * asym - ffplus.re * ffplus.re);
570 } else
571 GetFFsingle(ffplus, 1);
572
573
574 long double eta = 1.0 / (kstar * fPionac);
575 long double hfun = GetH(eta);
576 dcomplex gtilde = GetG(eta, rho, hfun);
577 dcomplex gtilor = mult(gtilde, 1.0 / fRStarS);
578
579 dcomplex fcouls, fcoult;
580 Getfc(kstar, eta, hfun, fcouls, fcoult);
581
582 dcomplex fgs = mult(gtilor, fcouls);
583 // long double fgmods = modl2(fgs);
584
585 dcomplex fgt = mult(gtilor, fcoult);
586 // long double fgmodt = modl2(fgt);
587
588 dcomplex expikr;
589 expikr.re = cos(tKstRst);
590 expikr.im = -sin(tKstRst);
591 // dcomplex expikrc = conj(expikr);
592 // dcomplex ffplusc = conj(ffplus);
593
594 // dcomplex expikr2 = mult(expikr, expikr);
595 // dcomplex expikrc2 = conj(expikr2);
596
597 dcomplex epfpc = mult(expikr, ffplus);
598 // dcomplex emfmc = mult(expikrc, ffminus);
599
600 // long double fcgefhs = (fgs.re*emfmc.re + fgs.im*emfmc.im);
601 // long double fcgefgs = (fgs.re*epfpc.re + fgs.im*epfpc.im);
602 // long double fcgefgt = (fgt.re*epfpc.re + fgt.im*epfpc.im);
603 // long double fcgefht = (fgt.re*emfmc.re + fgt.im*emfmc.im);
604
605 dcomplex fvfs;
606 fvfs.re = epfpc.re + fgs.re;
607 fvfs.im = epfpc.im + fgs.im;
608
609 dcomplex fvft;
610 fvft.re = epfpc.re + fgt.re;
611 fvft.im = epfpc.im + fgt.im;
612
613 // long double smult = 1+wavesign;
614
615 // return 0.5 * Gamow(kstar);
616 // return (0.5 * Gamow(kstar) *
617 // (2.0 * fgmods * smult + modl2(ffplus) + modl2(ffminus) +
618 // wavesign*sterm.re + wavesign*tterm.re +
619 // smult * 2 * (fcgefhs + fcgefgs)));
620
621 if (fTwospin == 1) {
622 // wavesign = 1;
623 // smult = 2;
624 // long double singlet = (0.5 * Gamow(kstar) *
625 // (2.0 * fgmods * smult + modl2(ffplus) + modl2(ffminus) +
626 // wavesign*sterm.re + wavesign*tterm.re +
627 // smult * 2 * (fcgefhs + fcgefgs)));
628 // wavesign = -1;
629 // smult = 0;
630 // long double triplet = (0.5 * Gamow(kstar) *
631 // (2.0 * fgmodt * smult + modl2(ffplus) + modl2(ffminus) +
632 // wavesign*sterm.re + wavesign*tterm.re +
633 // smult * 2 * (fcgefht + fcgefgt)));
634 // scmp = singlet;
635 // tcmp = triplet;
636
637 // ccmp = 0.5 * Gamow(kstar) * (modl2(ffplus) + modl2(ffminus));
638 // gcmp = fgmod;
639
640 long double singlet = Gamow(kstar) * modl2(fvfs);
641 long double triplet = Gamow(kstar) * modl2(fvft);
642
643 // std::cout << " s t " << singlet << " " << triplet << " - " << (0.25*singlet+0.75*triplet) << std::endl;
644
645 return (0.25 * singlet + 0.75 * triplet);
646 // return triplet;
647 } else
648 return Gamow(kstar) * modl2(fvfs);
649 // return (0.5 * Gamow(kstar) *
650 // (2.0 * fgmods * smult + modl2(ffplus) + modl2(ffminus) +
651 // wavesign*sterm.re + wavesign*tterm.re +
652 // smult * 2 * (fcgefhs + fcgefgs)));
653 }
654
655 double FemtoWeightGeneratorKisiel::GetStrong() {
656 double fr, fi;
657 double dr;
658 double ar, ai;
659 double ir, ii;
660
661 double sr, si;
662 double tr, ti;
663
664 double skr, ckr;
665 double srk, crk;
666
667 double d0_sr = fD0s.re;
668 double d0_si = fD0s.im;
669 double d0_tr = fD0t.re;
670 double d0_ti = fD0t.im;
671
672 double f0_sr = fF0s.re;
673 double f0_si = fF0s.im;
674 double f0_tr = fF0t.re;
675 double f0_ti = fF0t.im;
676
677 double kstar = fabs(fKStar);
678 double rstar = fRStarS;
679 double kstrst = fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS;
680 /*
681 if (pairtype == 8) {
682 f0_sr = 2.88 /0.197327;
683 f0_si = 0.0 /0.197327;
684 d0_sr = 2.92 /0.197327;
685 d0_si = 0.0 /0.197327;
686 f0_tr = 1.66 /0.197327;
687 f0_ti = 0.0 /0.197327;
688 d0_tr = 3.78 /0.197327;
689 d0_ti = 0.0 /0.197327;
690 }
691 else if (pairtype == 9) {
692 f0_sr = -1.20 /0.197327;
693 f0_si = 0.37 /0.197327;
694 d0_sr = 0.0 /0.197327;
695 d0_si = 0.0 /0.197327;
696 f0_tr = -1.20 /0.197327;
697 f0_ti = 0.37 /0.197327;
698 d0_tr = 0.0 /0.197327;
699 d0_ti = 0.0 /0.197327;
700 }
701 */
702 if (rstar < 1.0 / 0.197327) return 1.0;
703
704 fr = f0_sr / (f0_sr * f0_sr + f0_si * f0_si);
705 fi = -f0_si / (f0_sr * f0_sr + f0_si * f0_si);
706
707 dr = 0.5 * d0_sr * kstar * kstar;
708
709 ir = fr + dr;
710 ii = fi - kstar;
711
712 ar = ir / (ir * ir + ii * ii);
713 ai = -ii / (ir * ir + ii * ii);
714
715 srk = sin(kstar * rstar);
716 crk = cos(kstar * rstar);
717
718 sr = (ar * crk - ai * srk) / rstar;
719 si = (ai * crk + ar * srk) / rstar;
720
721 fr = f0_tr / (f0_tr * f0_tr + f0_ti * f0_ti);
722 fi = -f0_ti / (f0_tr * f0_tr + f0_ti * f0_ti);
723
724 dr = 0.5 * d0_tr * kstar * kstar;
725
726 ir = fr + dr;
727 ii = fi - kstar;
728
729 ar = ir / (ir * ir + ii * ii);
730 ai = -ii / (ir * ir + ii * ii);
731
732 tr = (ar * crk - ai * srk) / rstar;
733 ti = (ai * crk + ar * srk) / rstar;
734
735 skr = sin(kstrst);
736 ckr = cos(kstrst);
737
738 return (0.25 * ((ckr + sr) * (ckr + sr) + (skr + si) * (skr + si))
739 + 0.75 * ((ckr + tr) * (ckr + tr) + (skr + ti) * (skr + ti)));
740 }
741
742 double FemtoWeightGeneratorKisiel::GetQuantumStrong() {
743 double fr, fi;
744 double dr;
745 double ar, ai;
746 double ir, ii;
747
748 double sgrp, sgip, sgrm, sgim;
749 double trrp, trip, trrm, trim;
750
751 double skr, ckr;
752 double skm, ckm;
753 double srk, crk;
754 double srm, crm;
755
756 // double d0_sr;
757 // double d0_si;
758 // double d0_tr;
759 // double d0_ti;
760
761 // double f0_sr;
762 // double f0_si;
763 // double f0_tr;
764 // double f0_ti;
765
766 double kstar = fabs(fKStar);
767 double rstar = fRStarS;
768 double kstrst = fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS;
769
770 // std::cout << "k o s l " << kstar << " " << fKStarOut << " " << fKStarSide << " " << fKStarLong << std::endl;
771 // std::cout << "r o s l " << fRStarS << " " << fRStarOutS << " " << fRStarSideS << " " << fRStarLongS << std::endl;
772 // std::cout << "kr " << kstrst << std::endl;
773
774 double snr, sni, trr, tri;
775
776 double singlet, triplet;
777
778 // if (pairtype == 8) {
779 // f0_sr = 2.88 /0.197327;
780 // f0_si = 0.0 /0.197327;
781 // d0_sr = 2.92 /0.197327;
782 // d0_si = 0.0 /0.197327;
783 // f0_tr = 1.66 /0.197327;
784 // f0_ti = 0.0 /0.197327;
785 // d0_tr = 3.78 /0.197327;
786 // d0_ti = 0.0 /0.197327;
787 // }
788 // else if (pairtype == 9) {
789 // f0_sr = -1.20 /0.197327;
790 // f0_si = 0.37 /0.197327;
791 // d0_sr = 0.0 /0.197327;
792 // d0_si = 0.0 /0.197327;
793 // f0_tr = -1.20 /0.197327;
794 // f0_ti = 0.37 /0.197327;
795 // d0_tr = 0.0 /0.197327;
796 // d0_ti = 0.0 /0.197327;
797 // }
798
799 if (rstar < 1.0 / 0.197327) return 1.0;
800
801 fr = fF0s.re / (fF0s.re * fF0s.re + fF0s.im * fF0s.im);
802 fi = -fF0s.im / (fF0s.re * fF0s.re + fF0s.im * fF0s.im);
803
804 dr = 0.5 * fD0s.re * kstar * kstar;
805
806 ir = fr + dr;
807 ii = fi - kstar;
808
809 ar = ir / (ir * ir + ii * ii);
810 ai = -ii / (ir * ir + ii * ii);
811
812 srk = sin(kstar * rstar);
813 crk = cos(kstar * rstar);
814
815 sgrp = (ar * crk - ai * srk) / rstar;
816 sgip = (ai * crk + ar * srk) / rstar;
817
818 srm = sin(kstar * rstar);
819 crm = cos(kstar * rstar);
820
821 sgrm = (ar * crm - ai * srm) / rstar;
822 sgim = (ai * crm + ar * srm) / rstar;
823
824 skr = sin(kstrst);
825 ckr = cos(kstrst);
826
827 skm = sin(-kstrst);
828 ckm = cos(-kstrst);
829
830 snr = ckr + sgrp + ckm + sgrm;
831 sni = skr + sgip + skm + sgim;
832
833 singlet = 0.5 * (snr * snr + sni * sni);
834
835 fr = fF0t.re / (fF0t.re * fF0t.re + fF0t.im * fF0t.im);
836 fi = -fF0t.im / (fF0t.re * fF0t.re + fF0t.im * fF0t.im);
837
838 dr = 0.5 * fD0t.re * kstar * kstar;
839
840 ir = fr + dr;
841 ii = fi - kstar;
842
843 ar = ir / (ir * ir + ii * ii);
844 ai = -ii / (ir * ir + ii * ii);
845
846 trrp = (ar * crk - ai * srk) / rstar;
847 trip = (ai * crk + ar * srk) / rstar;
848
849 trrm = (ar * crm - ai * srm) / rstar;
850 trim = (ai * crm + ar * srm) / rstar;
851
852 trr = ckr + trrp - ckm - trrm;
853 tri = skr + trip - skm - trim;
854
855 triplet = 0.5 * (trr * trr + tri * tri);
856
857 // return (0.25 * ((ckr+sr)*(ckr+sr) + (skr+si)*(skr+si)) +
858 // 0.75 * ((ckr+tr)*(ckr+tr) + (skr+ti)*(skr+ti)));
859
860 // std::cout << "sng tri full " << singlet << " "<< triplet << " " << (0.25 * singlet + 0.75 * triplet) << std::endl;
861
862 return (0.25 * singlet + 0.75 * triplet);
863 }
864
865 double FemtoWeightGeneratorKisiel::GetQuantum() {
866 double quantumweight;
867
868 if (fTwospin == 0) {
869 quantumweight = 1.0 + TMath::Cos(2.0 * (fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS));
870 } else if (fTwospin == 1) {
871 if (fPcount == 3) {
872 quantumweight = 1.0 + TMath::Cos(2.0 * (fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS));
873 } else {
874 quantumweight = 1.0 - TMath::Cos(2.0 * (fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS));
875 }
876 fPcount++;
877 if (fPcount == 4) fPcount = 0;
878 }
879
880 return quantumweight;
881 }
882
883 FemtoWeightGeneratorKisiel::dcomplex
884 FemtoWeightGeneratorKisiel::GetG(long double eta, long double rho, long double hfun) const {
885 dcomplex gtemp;
886 long double bres, pres;
887 long double bmult;
888
889 Bfunpfun(eta, rho, bres, pres);
890 bmult = 2.0 * eta * rho * bres;
891
892 gtemp.re = pres + bmult * (log(fabs(2.0 * eta * rho)) + 2.0 * fEuler - 1.0 + hfun);
893 gtemp.im = bmult * Chiim(eta);
894
895 return gtemp;
896 }
897
898 void
899 FemtoWeightGeneratorKisiel::Getfc(long double kstar, long double eta, long double hfun, dcomplex& fcs, dcomplex& fct) const {
900 dcomplex ci;
901 dcomplex cia;
902
903 dcomplex fis;
904 dcomplex dis;
905 dcomplex fcinvs;
906
907 dcomplex fit;
908 dcomplex dit;
909 dcomplex fcinvt;
910
911 ci.re = hfun;
912 ci.im = Chiim(eta);
913 cia = mult(ci, 2.0 / fPionac);
914
915 fis = invr(fF0s);
916 dis = mult(fD0s, 0.5 * kstar * kstar);
917 fcinvs.re = fis.re + dis.re - cia.re;
918 fcinvs.im = fis.im + dis.im - cia.im;
919 fcs = invr(fcinvs);
920
921 fit = invr(fF0t);
922 dit = mult(fD0t, 0.5 * kstar * kstar);
923 fcinvt.re = fit.re + dit.re - cia.re;
924 fcinvt.im = fit.im + dit.im - cia.im;
925 fct = invr(fcinvt);
926 }
927
928 void FemtoWeightGeneratorKisiel::Bfunpfun(long double eta, long double rho, long double& bret, long double& pret) const {
929 long double b0 = 1;
930 long double b1 = eta * rho;
931 long double bsum = b0 + b1;
932 long double bnpu, bn, bnmu;
933 long double p0 = 1.0;
934 long double p1 = 0.0;
935 long double psum = p0 + p1;
936 long double pnpu, pn, pnmu;
937
938 if (rho > TMath::TwoPi() * 2.0) {
939 bret = sin(rho) / rho;
940 pret = cos(rho);
941 return;
942 }
943
944
945 bn = b1;
946 bnmu = b0;
947 pn = p1;
948 pnmu = p0;
949 for (int iter = 1; iter < 100000; iter++) {
950 bnpu = 2 * eta * rho * bn - rho * rho * bnmu;
951 bnpu /= (1.0 * iter + 1.0) * (1.0 * iter + 2.0);
952 bsum += bnpu;
953 // std::cout << "B E " << iter << " " << bnpu << std::endl;
954
955 pnpu = 2 * eta * rho * pn - rho * rho * pnmu - (2.0 * iter + 1.0) * 2.0 * eta * rho * bn;
956 pnpu /= (1.0 * iter) * (1.0 * iter + 1.0);
957 psum += pnpu;
958 // std::cout << "P E " << iter << " " << pnpu << std::endl;
959
960 bnmu = bn;
961 bn = bnpu;
962
963 pnmu = pn;
964 pn = pnpu;
965 if ((fabs(pnmu) + fabs(bnmu) + fabs(bnpu) + fabs(pnpu)) < 1.0e-20) {
966 // std::cout << "iter " << iter << std::endl;
967 break;
968 }
969 }
970
971
972 bret = bsum;
973 pret = psum;
974 }
975
976 long double FemtoWeightGeneratorKisiel::GetH(long double eta) const {
977 long double etasum = log(1.0 / eta) - fEuler;
978 long double series = 0.0;
979 long double x2inv = (eta * eta);
980 long double element;
981 long double save;
982 for (int iter = 1; iter < 1000000; iter++) {
983 element = ((1.0 * iter) * (1.0 * iter) + x2inv) * (1.0 * iter);
984 // std::cout << "Element " << iter << " is " << element << std::endl;
985 element = 1.0 / element;
986 if (iter == 1) save = 1.0e-10 * element;
987 // std::cout << "Element " << iter << " is " << element << std::endl;
988
989 series += element;
990 if (element < save) break;
991 }
992 series *= x2inv;
993 etasum += series;
994
995 return etasum;
996 }
997
998 double FemtoWeightGeneratorKisiel::GetCoulomb() {
999 double kstrst = fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS;
1000
1001 // Classical limit - if distance is larger than Coulomb radius,
1002 // the interaction does not matter
1003 if (fabs(fRStarS) > fabs(fPionac)) return (1.0);
1004 if (fabs(fRStarS) == 0.0) return (Gamow(fabs(fKStar)));
1005
1006 // Classical limit - in the case of large k*r* product we go to
1007 // classical coulomb interaction
1008 if (fabs(fKStar) * fRStarS * (1.0 + kstrst / (fRStarS * fabs(fKStar))) > 15.0)
1009 return (1.0 - 1.0 / (fRStarS * fPionac * fKStar * fKStar));
1010
1011 // Calculate the F function
1012 dcomplex ffplus;
1013 GetFFsingle(ffplus);
1014
1015 if (!std::isfinite(ffplus.re)) {
1016 std::cout << "FFPlus Re not a number !"
1017 << " " << std::endl;
1018 std::cout << fRStarS << " " << fKStar << " " << fRStarOutS << " " << fRStarSideS << " " << fRStarLongS << std::endl;
1019 }
1020
1021 if (!std::isfinite(ffplus.im))
1022 std::cout << "FFPlus Im not a number !"
1023 << " " << std::endl;
1024
1025
1026 return (Gamow(fabs(fKStar)) * modl2(ffplus));
1027 }
1028
1029 double FemtoWeightGeneratorKisiel::GetQuantumCoulomb() {
1030 if (fRStarS < 0.0000000001) return 1.0;
1031
1032 double kstrst = fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS;
1033 int ccase = 0;
1034 int wavesign = 1;
1035
1036 if (fTwospin == 1) {
1037 if (fPcount == 3)
1038 wavesign = 1;
1039 else
1040 wavesign = -1;
1041 fPcount++;
1042 if (fPcount == 4) fPcount = 0;
1043 }
1044
1045 // Classical limit - if distance is larger than Coulomb radius,
1046 // the interaction does not matter
1047 if (fabs(fRStarS) > fabs(fPionac)) return (1.0 + wavesign * cos(2 * kstrst));
1048
1049 // Classical limit - in the case of large k* we go to
1050 // classical coulomb interaction
1051 long double testp = fabs(fKStar) * fRStarS * (1.0 + kstrst / (fRStarS * fabs(fKStar)));
1052 long double testm = fabs(fKStar) * fRStarS * (1.0 - kstrst / (fRStarS * fabs(fKStar)));
1053
1054 if ((testp > 15.0) && (testm > 15.0)) {
1055 double fasymplus = (1.0 - 1.0 / ((fRStarS + kstrst) * fPionac * fKStar * fKStar));
1056 double fasymminus = (1.0 - 1.0 / ((fRStarS - kstrst) * fPionac * fKStar * fKStar));
1057 return 0.5 * ((fasymplus + fasymminus) * cos(2 * kstrst) + (2.0 * sqrt(fasymplus * fasymminus)));
1058 }
1059 // return (1.0 - 1.0/(fRStarS*pionac*fKStar*fKStar))*(1.0+wavesign*cos(2*kstrst));
1060
1061 dcomplex ffplus, ffminus;
1062 // Check for the classical limit in both functions separately
1063 if (((testp < 15.0) && (testm < 15.0))) // ||
1064 // ((testp< 15.0) && (testm> 15.0) && (fabs(testp-testm < 1.0))) ||
1065 // ((testp> 15.0) && (testm< 15.0) && (fabs(testp-testm < 1.0))))
1066 {
1067 // Calculate the F function
1068 GetFFdouble(ffplus, ffminus);
1069 ccase = 1;
1070 } else if (testp < 15.0) {
1071 double asym;
1072 GetFFsingle(ffplus);
1073 // GetFFsingle(&ffminus,-1);
1074 asym =
1075 (1.0 - 1.0 / (fRStarS * (1.0 - kstrst / (fRStarS * fabs(fKStar)) * fPionac * fKStar * fKStar))) / Gamow(fabs(fKStar));
1076 asym = sqrt(asym);
1077 if (asym < 1.0)
1078 ffminus.re = 1.0 + (asym - 1.0) * 2.0;
1079 else
1080 ffminus.re = 1.0 + (asym - 1.0) / 2.0;
1081 ffminus.im = sqrt(asym * asym - ffminus.re * ffminus.re);
1082 ccase = 2;
1083 } else {
1084 double asym;
1085 // GetFFsingle(&ffplus);
1086 GetFFsingle(ffminus, -1);
1087 asym =
1088 (1.0 - 1.0 / (fRStarS * (1.0 + kstrst / (fRStarS * fabs(fKStar)) * fPionac * fKStar * fKStar))) / Gamow(fabs(fKStar));
1089 asym = sqrt(asym);
1090 if (asym < 1.0)
1091 ffplus.re = 1.0 + (asym - 1.0) * 2.0;
1092 else
1093 ffplus.re = 1.0 + (asym - 1.0) / 2.0;
1094 ffplus.im = sqrt(asym * asym - ffplus.re * ffplus.re);
1095 ccase = 3;
1096 }
1097
1098 dcomplex expikr;
1099 expikr.re = cos(kstrst);
1100 expikr.im = sin(kstrst);
1101 dcomplex expikrc = conj(expikr);
1102 dcomplex ffplusc = conj(ffplus);
1103 dcomplex ffminusc = conj(ffminus);
1104
1105 dcomplex expikr2 = mult(expikr, expikr);
1106 dcomplex expikrc2 = conj(expikr2);
1107 dcomplex sterm = mult(expikr2, mult(ffplus, ffminusc));
1108 dcomplex tterm = mult(expikrc2, mult(ffminus, ffplusc));
1109
1110
1111 if (!std::isfinite(ffminus.re))
1112 std::cout << "FFMinus Re not a number !"
1113 << " " << ccase << std::endl;
1114
1115 if (!std::isfinite(ffminus.im))
1116 std::cout << "FFMinus Im not a number !"
1117 << " " << ccase << std::endl;
1118
1119 if ((ffplus.re > 2.0) || (ffplus.re < -2.0)) std::cout << std::endl << "FFplus Re wild !" << ffplus.re << std::endl;
1120
1121 if ((ffplus.im > 2.0) || (ffplus.im < -2.0)) std::cout << "FFplus Im wild !" << ffplus.im << std::endl;
1122
1123 if ((ffminus.re > 2.0) || (ffminus.re < -2.0))
1124 std::cout << std::endl << "FFminus Re wild !" << ffminus.re << " " << ccase << std::endl;
1125
1126 if ((ffminus.im > 2.0) || (ffminus.im < -2.0)) std::cout << "FFminus Im wild !" << ffminus.im << " " << ccase << std::endl;
1127
1128 fCoulqscpart = 0.5 * Gamow(fabs(fKStar)) * (modl2(ffplus) + modl2(ffminus));
1129
1130 return (0.5 * Gamow(fabs(fKStar)) * (modl2(ffplus) + wavesign * sterm.re + wavesign * tterm.re + modl2(ffminus)));
1131 }
1132
1133 void FemtoWeightGeneratorKisiel::GetFFdouble(dcomplex& ffp, dcomplex& ffm) const {
1134 long double comprep[COULOMBSTEPS];
1135 long double compimp[COULOMBSTEPS];
1136 long double comprem[COULOMBSTEPS];
1137 long double compimm[COULOMBSTEPS];
1138 long double eta, ksip, ksim;
1139 dcomplex alfa, zetp, zetm;
1140
1141 int nsteps;
1142
1143 long double kstar = fabs(fKStar);
1144 long double kstrst = fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS;
1145 long double coskr = kstrst / (fabs(fKStar) * fRStarS);
1146
1147 if ((kstar * fRStarS * (1 + coskr) < 5.0) && (kstar * fRStarS * (1 - coskr) < 5.0))
1148 nsteps = 25;
1149 else if ((kstar * fRStarS * (1 + coskr) < 10.0) && (kstar * fRStarS * (1 - coskr) < 10.0))
1150 nsteps = 45;
1151 else if ((kstar * fRStarS * (1 + coskr) < 15.0) && (kstar * fRStarS * (1 - coskr) < 15.0))
1152 nsteps = 110;
1153 else
1154 nsteps = 150;
1155
1156 eta = 1.0 / (kstar * fPionac);
1157 alfa.re = 0.0;
1158 alfa.im = -eta;
1159
1160 dcomplex fcomp, scompp, scompm;
1161 long double tcomp;
1162 dcomplex sump, summ;
1163 dcomplex fcmult;
1164
1165 long double rad = fRStarS;
1166
1167 ksip = kstar * rad * (1 + coskr);
1168 ksim = kstar * rad * (1 - coskr);
1169
1170 zetp.re = 0.0;
1171 zetp.im = ksip;
1172
1173 zetm.re = 0.0;
1174 zetm.im = ksim;
1175
1176 fcomp.re = 1.0;
1177 fcomp.im = 0.0;
1178 scompp.re = 1.0;
1179 scompp.im = 0.0;
1180 scompm.re = 1.0;
1181 scompm.im = 0.0;
1182 tcomp = 1.0;
1183
1184 for (int istep = 0; istep < nsteps; istep++) {
1185 sump = mult(fcomp, scompp);
1186 summ = mult(fcomp, scompm);
1187
1188 sump = mult(sump, 1.0 / (tcomp * tcomp));
1189 summ = mult(summ, 1.0 / (tcomp * tcomp));
1190
1191
1192 if (istep == 0) {
1193 comprep[istep] = sump.re;
1194 compimp[istep] = sump.im;
1195
1196 comprem[istep] = summ.re;
1197 compimm[istep] = summ.im;
1198 } else {
1199 comprep[istep] = comprep[istep - 1] + sump.re;
1200 compimp[istep] = compimp[istep - 1] + sump.im;
1201
1202 comprem[istep] = comprem[istep - 1] + summ.re;
1203 compimm[istep] = compimm[istep - 1] + summ.im;
1204 }
1205
1206 fcmult.re = alfa.re + istep;
1207 fcmult.im = alfa.im;
1208
1209 fcomp = mult(fcomp, fcmult);
1210 scompp = mult(scompp, zetp);
1211 scompm = mult(scompm, zetm);
1212 tcomp *= (istep + 1);
1213 }
1214
1215 ffp.re = comprep[nsteps - 1];
1216 ffp.im = compimp[nsteps - 1];
1217
1218 ffm.re = comprem[nsteps - 1];
1219 ffm.im = compimm[nsteps - 1];
1220 }
1221
1222 void FemtoWeightGeneratorKisiel::GetFFsingle(dcomplex& ffp, int sign) const {
1223 double comprep[COULOMBSTEPS];
1224 double compimp[COULOMBSTEPS];
1225 double eta, ksip;
1226 dcomplex alfa, zetp;
1227
1228 int nsteps;
1229
1230 double kstar = fabs(fKStar);
1231 double kstrst = fKStarOut * fRStarOutS + fKStarSide * fRStarSideS + fKStarLong * fRStarLongS;
1232 double coskr = sign * kstrst / (fabs(fKStar) * fRStarS);
1233
1234 if (kstar * fRStarS * (1 + coskr) > 15.0)
1235 nsteps = 170;
1236 else if (kstar * fRStarS * (1 + coskr) > 10.0)
1237 nsteps = 45;
1238 else if (kstar * fRStarS * (1 + coskr) > 5.0)
1239 nsteps = 35;
1240 else
1241 nsteps = 15;
1242
1243 eta = 1.0 / (kstar * fPionac);
1244 alfa.re = 0.0;
1245 alfa.im = -eta;
1246
1247 dcomplex fcomp, scompp;
1248 double tcomp;
1249 dcomplex sump;
1250 dcomplex fcmult;
1251
1252 double rad = fRStarS;
1253
1254 ksip = kstar * rad * (1 + coskr);
1255
1256 zetp.re = 0.0;
1257 zetp.im = ksip;
1258
1259 fcomp.re = 1.0;
1260 fcomp.im = 0.0;
1261 scompp.re = 1.0;
1262 scompp.im = 0.0;
1263 tcomp = 1.0;
1264
1265 for (int istep = 0; istep < nsteps; istep++) {
1266 sump = mult(fcomp, scompp);
1267
1268 sump = mult(sump, 1.0 / (tcomp * tcomp));
1269
1270 if (istep == 0) {
1271 comprep[istep] = sump.re;
1272 compimp[istep] = sump.im;
1273 } else {
1274 comprep[istep] = comprep[istep - 1] + sump.re;
1275 compimp[istep] = compimp[istep - 1] + sump.im;
1276 }
1277
1278 fcmult.re = alfa.re + istep;
1279 fcmult.im = alfa.im;
1280
1281 fcomp = mult(fcomp, fcmult);
1282 scompp = mult(scompp, zetp);
1283 tcomp *= (istep + 1);
1284
1285 if ((sump.re * sump.re + sump.im * sump.im) < 1.0e-14) {
1286 nsteps = istep;
1287 break;
1288 }
1289 }
1290
1291 ffp.re = comprep[nsteps - 1];
1292 ffp.im = compimp[nsteps - 1];
1293 }
1294
1295 FemtoWeightGeneratorKisiel::dcomplex FemtoWeightGeneratorKisiel::conj(const dcomplex& arg) const {
1296 dcomplex res;
1297
1298 res.re = arg.re;
1299 res.im = -arg.im;
1300
1301 return res;
1302 }
1303
1304 FemtoWeightGeneratorKisiel::dcomplex FemtoWeightGeneratorKisiel::mult(const dcomplex& arga, long double argb) const {
1305 dcomplex res;
1306
1307 res.re = arga.re * argb;
1308 res.im = arga.im * argb;
1309
1310 return res;
1311 }
1312
1313 FemtoWeightGeneratorKisiel::dcomplex FemtoWeightGeneratorKisiel::mult(const dcomplex& arga, const dcomplex& argb) const {
1314 dcomplex res;
1315
1316 res.re = arga.re * argb.re - arga.im * argb.im;
1317 res.im = arga.re * argb.im + argb.re * arga.im;
1318
1319 return res;
1320 }
1321
1322 Package* FemtoWeightGeneratorKisiel::Report() const {
1323 Package* pack = FemtoWeightGenerator::Report();
1324 pack->AddObject(new ParameterString("PairType", Femto::PairTypeToString(fPairType)));
1325 return pack;
1326 }
1327
1328 Bool_t FemtoWeightGeneratorKisiel::IsPairSupported(Femto::EPairType type) const {
1329 if (static_cast<int>(type) >= 100) return kTRUE;
1330 return kFALSE;
1331 }
1332
1333 FemtoWeightGeneratorKisiel::~FemtoWeightGeneratorKisiel() {}
1334} // namespace Hal
Hal::FemtoPair::GetPdg2
Int_t GetPdg2() const
Definition FemtoPair.h:302
Hal::FemtoPair::GetPdg1
Int_t GetPdg1() const
Definition FemtoPair.h:297