BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
D0Topipipi0.cxx
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of models developed at BES collaboration
5// based on the EvtGen framework. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/BesCopyright
9// Copyright (A) 2006 Ping Rong-Gang @IHEP
10//
11// Module: EvtD0Topipipi0.cc
12//
13// Description: Model provided by user, see BAM-628
14//
15// Modification history:
16//
17// Liaoyuan Dong Aug. 11, 2022 Module created
18//
19//------------------------------------------------------------------------
20#include "QCMCFilterAlg/D0Topipipi0.h"
21#include <stdlib.h>
22#include <iostream>
23#include <fstream>
24#include <vector>
25#include <map>
26#include <string>
27#include <complex>
28#include <vector>
29#include <math.h>
30using namespace std;
31#define PI acos(-1)
32
33D0Topipipi0::~D0Topipipi0() {}
34
35void D0Topipipi0::init(){
36 //std::cout << "Initializing D0Topipipi0: charm= "<<charm<<" tagmode= "<<tagmode<<std::endl;
37
38
39 g_uv.clear();
40 for(int i=0; i<4; i++){
41 for(int j=0; j<4; j++){
42 if(i!=j){
43 g_uv.push_back(0.0);
44 }else if(i<3){
45 g_uv.push_back(-1.0);
46 }else if(i==3){
47 g_uv.push_back(1.0);
48 }
49 }
50 }
51
52 epsilon_uvmn.clear();
53 for(int i=0; i<4; i++){
54 for(int j=0; j<4; j++){
55 for(int k=0; k<4; k++){
56 for(int l=0; l<4; l++){
57 if(i==j || i==k || i==l || j==k || j==l || k==l){
58 epsilon_uvmn.push_back(0.0);
59 }else{
60 if(i==0 && j==1 && k==2 && l==3) epsilon_uvmn.push_back(1.0);
61 if(i==0 && j==1 && k==3 && l==2) epsilon_uvmn.push_back(-1.0);
62 if(i==0 && j==2 && k==1 && l==3) epsilon_uvmn.push_back(-1.0);
63 if(i==0 && j==2 && k==3 && l==1) epsilon_uvmn.push_back(1.0);
64 if(i==0 && j==3 && k==1 && l==2) epsilon_uvmn.push_back(1.0);
65 if(i==0 && j==3 && k==2 && l==1) epsilon_uvmn.push_back(-1.0);
66
67 if(i==1 && j==0 && k==2 && l==3) epsilon_uvmn.push_back(-1.0);
68 if(i==1 && j==0 && k==3 && l==2) epsilon_uvmn.push_back(1.0);
69 if(i==1 && j==2 && k==0 && l==3) epsilon_uvmn.push_back(1.0);
70 if(i==1 && j==2 && k==3 && l==0) epsilon_uvmn.push_back(-1.0);
71 if(i==1 && j==3 && k==0 && l==2) epsilon_uvmn.push_back(-1.0);
72 if(i==1 && j==3 && k==2 && l==0) epsilon_uvmn.push_back(1.0);
73
74 if(i==2 && j==0 && k==1 && l==3) epsilon_uvmn.push_back(1.0);
75 if(i==2 && j==0 && k==3 && l==1) epsilon_uvmn.push_back(-1.0);
76 if(i==2 && j==1 && k==0 && l==3) epsilon_uvmn.push_back(-1.0);
77 if(i==2 && j==1 && k==3 && l==0) epsilon_uvmn.push_back(1.0);
78 if(i==2 && j==3 && k==0 && l==1) epsilon_uvmn.push_back(1.0);
79 if(i==2 && j==3 && k==1 && l==0) epsilon_uvmn.push_back(-1.0);
80
81 if(i==3 && j==0 && k==1 && l==2) epsilon_uvmn.push_back(-1.0);
82 if(i==3 && j==0 && k==2 && l==1) epsilon_uvmn.push_back(1.0);
83 if(i==3 && j==1 && k==0 && l==2) epsilon_uvmn.push_back(1.0);
84 if(i==3 && j==1 && k==2 && l==0) epsilon_uvmn.push_back(-1.0);
85 if(i==3 && j==2 && k==0 && l==1) epsilon_uvmn.push_back(-1.0);
86 if(i==3 && j==2 && k==1 && l==0) epsilon_uvmn.push_back(1.0);
87
88 }
89 }
90 }
91 }
92 }
93
94 _nd = 3;
95 math_pi = 3.1415926;
96 mass_Pion = 0.13957;
97
98 rRes = 3.0*0.197321;
99 rD = 5.0*0.197321;
100 m_Pi = mass_Pion;
101 m2_Pi = m_Pi*m_Pi;
102 m_Pi0 = 0.134977;
103 m2_Pi0 = m_Pi0*m_Pi0;
104
105 m0_rho7700 = 0.77526;
106 w0_rho7700 = 0.1478;
107
108 m0_rho770p = 0.77511;
109 w0_rho770p = 0.1491;
110
111 m0_f21270 = 1.2755;
112 w0_f21270 = 0.1867;
113
114 s0_prod = -5.0;
115
116 fitpara.clear();
117 fitpara.push_back(complex<double>(100,0));
118 fitpara.push_back(complex<double>(69.8939*cos(3.14983),69.8939*sin(3.14983)));
119 fitpara.push_back(complex<double>(58.521*cos(-2.90685),58.521*sin(-2.90685)));
120 fitpara.push_back(complex<double>(483.035*cos(-0.679009),483.035*sin(-0.679009)));
121 fitpara.push_back(complex<double>(441.921*cos(-0.879847),441.921*sin(-0.879847)));
122 fitpara.push_back(complex<double>(1356.95*cos(-0.206653),1356.95*sin(-0.206653)));
123 fitpara.push_back(complex<double>(559.218*cos(0.501728),559.218*sin(0.501728)));
124 fitpara.push_back(complex<double>(3165.25*cos(3.39939),3165.25*sin(3.39939)));
125 fitpara.push_back(complex<double>(1422.93*cos(3.05347),1422.93*sin(3.05347)));
126 fitpara.push_back(complex<double>(2399.8*cos(2.24983),2399.8*sin(2.24983)));
127 fitpara.push_back(complex<double>(4601.19*cos(2.74388),4601.19*sin(2.74388)));
128 fitpara.push_back(complex<double>(1684.1*cos(1.99894),1684.1*sin(1.99894)));
129 fitpara.push_back(complex<double>(678.674*cos(-2.510691),678.674*sin(-2.510691)));
130 fitpara.push_back(complex<double>(2.19068*cos(0.991805),2.19068*sin(0.991805)));
131
132
133 return;
134}
135
136vector<double> D0Topipipi0::sum_tensor(vector<double> pa, vector<double> pb) {
137 if(pa.size()!=pb.size()){
138 cout<<"error sum tensor"<<endl;
139 exit(0);
140 }
141 vector<double> temp; temp.clear();
142 for(int i=0; i<pa.size(); i++){
143 double sum = pa[i] + pb[i];
144 temp.push_back(sum);
145 }
146 return temp;
147}
148
149double D0Topipipi0::contract_11_0(vector<double> pa, vector<double> pb){
150 if(pa.size()!=pb.size() || pa.size()!=4) {
151 cout<<"error contract 11->0"<<endl;
152 exit(0);
153 }
154 double temp = pa[3]*pb[3] - pa[0]*pb[0] - pa[1]*pb[1] - pa[2]*pb[2];
155 return temp;
156
157}
158
159vector<double> D0Topipipi0::contract_21_1(vector<double> pa, vector<double> pb){
160 if(pa.size()!=16 || pb.size()!=4) {
161 cout<<"error contract 21->1"<<endl;
162 exit(0);
163 }
164 vector<double> temp; temp.clear();
165 for(int i=0; i<4; i++){
166 double sum = 0;
167 for(int j=0; j<4; j++){
168 int idx = i*4+j;
169 sum += pa[idx]*pb[j]*g_uv[4*j+j];
170 }
171 temp.push_back(sum);
172 }
173 return temp;
174
175}
176
177double D0Topipipi0::contract_22_0(vector<double> pa, vector<double> pb){
178 if(pa.size()!=pb.size() || pa.size()!=16) {
179 cout<<"error contract 22->0"<<endl;
180 exit(0);
181 }
182 double temp = 0;
183 for(int i=0; i<4; i++){
184 for(int j=0; j<4; j++){
185 int idx = i*4+j;
186 temp += pa[idx]*pb[idx]*g_uv[4*i+i]*g_uv[4*j+j];
187 }
188 }
189 return temp;
190
191}
192
193vector<double> D0Topipipi0::contract_31_2(vector<double> pa, vector<double> pb){
194 if(pa.size()!=64 || pb.size()!=4) {
195 cout<<"error contract 31->2"<<endl;
196 exit(0);
197 }
198 vector<double> temp; temp.clear();
199 for(int i=0; i<16; i++){
200 double sum = 0;
201 for(int j=0; j<4; j++){
202 int idx = i*4+j;
203 sum += pa[idx]*pb[j]*g_uv[4*j+j];
204 }
205 temp.push_back(sum);
206 }
207 return temp;
208
209}
210
211vector<double> D0Topipipi0::contract_41_3(vector<double> pa, vector<double> pb){
212 if(pa.size()!=256|| pb.size()!=4) {
213 cout<<"error contract 41->3"<<endl;
214 exit(0);
215 }
216 vector<double> temp; temp.clear();
217 for(int i=0; i<64; i++){
218 double sum = 0;
219 for(int j=0; j<4; j++){
220 int idx = i*4+j;
221 sum += pa[idx]*pb[j]*g_uv[4*j+j];
222 }
223 temp.push_back(sum);
224 }
225 return temp;
226
227}
228
229vector<double> D0Topipipi0::contract_42_2(vector<double> pa, vector<double> pb){
230 if(pa.size()!=256|| pb.size()!=16) {
231 cout<<"error contract 42->2"<<endl;
232 exit(0);
233 }
234 vector<double> temp; temp.clear();
235 for(int i=0; i<16; i++){
236 double sum = 0;
237 for(int j=0; j<4; j++){
238 for(int k=0; k<4; k++){
239 int idxa = i*16+j*4+k;
240 int idxb = j*4+k;
241 sum += pa[idxa] * pb[idxb] * g_uv[4*j+j] * g_uv[4*k+k];
242 }
243 }
244 temp.push_back(sum);
245 }
246
247 return temp;
248
249}
250
251vector<double> D0Topipipi0::contract_22_2(vector<double> pa, vector<double> pb){
252 if(pa.size()!=16|| pb.size()!=16) {
253 cout<<"error contract 42->2"<<endl;
254 exit(0);
255 }
256 vector<double> temp; temp.clear();
257 for(int i=0; i<4; i++){
258 for(int j=0; j<4; j++){
259 double sum = 0;
260 for(int k=0; k<4; k++){
261 int idxa = i*4+k;
262 int idxb = j*4+k;
263 sum += pa[idxa] * pb[idxb] * g_uv[4*k+k];
264 }
265 temp.push_back(sum);
266 }
267 }
268
269 return temp;
270
271}
272
273/*
274 vector<double> D0Topipipi0::contract_11_2(vector<double> pa, vector<double> pb){
275 if(pa.size()!=pb.size() || pa.size()!=4) {
276 cout<<"error contract 11->2"<<endl;
277 exit(0);
278 }
279 vector<double> temp; temp.clear();
280 for(int i=0; i<4; i++){
281 for(int j=0; j<4; j++){
282 double prod = pa[i]*pb[j];
283 temp.push_back(prod);
284 }
285 }
286 return temp;
287 }
288
289 vector<double> D0Topipipi0::contract_22_4(vector<double> pa, vector<double> pb){
290 if(pa.size()!=pb.size() || pa.size()!=16) {
291 cout<<"error contract 22->4"<<endl;
292 exit(0);
293 }
294 vector<double> temp; temp.clear();
295 for(int i=0; i<16; i++){
296 for(int j=0; j<16; j++){
297 double prod = pa[i]*pb[j];
298 temp.push_back(prod);
299 }
300 }
301 return temp;
302 }
303 */
304
305//OrbitalTensors
306vector<double> D0Topipipi0::OrbitalTensors(vector<double> pa, vector<double> pb, vector<double> pc, double r, int rank)
307{
308
309 if(pa.size()!=4 || pb.size()!=4 || pc.size()!=4) {
310 cout<<"Error: pa, pb, pc"<<endl;
311 exit(0);
312 }
313 if(rank<0){
314 cout<<"Error: L<0 !!!"<<endl;
315 exit(0);
316 }
317
318 // relative momentum
319 vector<double> mr; mr.clear();
320
321 for(int i=0; i<4; i++){
322 double temp = pb[i] - pc[i];
323 mr.push_back(temp);
324 }
325
326 // "Masses" of particles
327 double msa = contract_11_0(pa, pa);
328 double msb = contract_11_0(pb, pb);
329 double msc = contract_11_0(pc, pc);
330
331 // Q^2_abc
332 double top = msa + msb - msc;
333 double Q2abc = top*top/(4.0*msa) - msb;
334
335 // Barrier factors
336 double Q_0 = 0.197321f/r;
337 double Q_02 = Q_0*Q_0;
338 double Q_04 = Q_02*Q_02;
339 // double Q_06 = Q_04*Q_02;
340 // double Q_08 = Q_04*Q_04;
341
342 double Q4abc = Q2abc*Q2abc;
343 // double Q6abc = Q4abc*Q2abc;
344 // double Q8abc = Q4abc*Q4abc;
345
346 double mB1 = sqrt(2.0f/(Q2abc + Q_02));
347 double mB2 = sqrt(13.0f/(Q4abc + (3.0f*Q_02)*Q2abc + 9.0f*Q_04));
348 // mB3 = &sqrt(277.0f/(Q6abc + (6.0f*Q_02)*Q4abc + (45.0f*Q_04)*Q2abc + 225.0f*Q_06));
349 // mB4 = &sqrt(12746.0f/(Q8abc + (10.0f*Q_02)*Q6abc + (135.0f*Q_04)*Q4abc + (1575.0f*Q_06)*Q2abc + 11025.0f*Q_08));
350
351 // Projection Operator 2-rank
352 vector<double> proj_uv; proj_uv.clear();
353 for(int i=0; i<4; i++){
354 for(int j=0; j<4; j++){
355 int idx = i*4 + j;
356 double temp = -g_uv[idx] + pa[i]*pa[j]/msa;
357 proj_uv.push_back(temp);
358 }
359 }
360
361 // Orbital Tensors
362 if(rank==0)
363 {
364 vector<double> t; t.clear();
365 t.push_back(1.0);
366 return t;
367
368 }else if(rank<3)
369 {
370 vector<double> t_u; t_u.clear();
371 vector<double> Bt_u; Bt_u.clear();
372 for(int i=0; i<4; i++){
373 double temp = 0;
374 for(int j=0; j<4; j++){
375 int idx = i*4 + j;
376 temp += -proj_uv[idx]*mr[j]*g_uv[j*4+j];
377 }
378 t_u.push_back(temp);
379 Bt_u.push_back(temp*mB1);
380 }
381 if(rank==1) return Bt_u;
382
383 double t_u2 = contract_11_0(t_u,t_u);
384
385 vector<double> Bt_uv; Bt_uv.clear();
386 for(int i=0; i<4; i++){
387 for(int j=0; j<4; j++){
388 int idx = 4*i + j;
389 double temp = t_u[i]*t_u[j] + (1.0/3.0)*proj_uv[idx]*t_u2;
390 Bt_uv.push_back(temp*mB2);
391 }
392 }
393 if(rank==2) return Bt_uv;
394
395 }else
396 {
397 cout<<"rank>2: please add it by yourself!!!"<<endl;
398 exit(0);
399 }
400
401}
402
403// projection Tensor
404vector<double> D0Topipipi0::ProjectionTensors(vector<double> pa, int rank)
405{
406
407 if(pa.size()!=4) {
408 cout<<"Error: pa"<<endl;
409 exit(0);
410 }
411 if(rank<0){
412 cout<<"Error: L<0 !!!"<<endl;
413 exit(0);
414 }
415
416 double msa = contract_11_0(pa, pa);
417
418 // Projection Operator 2-rank
419 vector<double> proj_uv; proj_uv.clear();
420 for(int i=0; i<4; i++){
421 for(int j=0; j<4; j++){
422 int idx = i*4 + j;
423 double temp = -g_uv[idx] + pa[i]*pa[j]/msa;
424 proj_uv.push_back(temp);
425 }
426 }
427
428 // Orbital Tensors
429 if(rank==0)
430 {
431 vector<double> t; t.clear();
432 t.push_back(1.0);
433 return t;
434
435 }else if(rank==1)
436 {
437 return proj_uv;
438 }else if(rank==2)
439 {
440 vector<double> proj_uvmn; proj_uvmn.clear();
441 for(int i=0; i<4; i++){
442 for(int j=0; j<4; j++){
443 for(int k=0; k<4; k++){
444 for(int l=0; l<4; l++){
445
446 int idx1_1 = 4*i + k;
447 int idx1_2 = 4*i + l;
448 int idx1_3 = 4*i + j;
449
450 int idx2_1 = 4*j + l;
451 int idx2_2 = 4*j + k;
452 int idx2_3 = 4*k + l;
453
454 double temp = (1.0/2.0)*(proj_uv[idx1_1]*proj_uv[idx2_1] + proj_uv[idx1_2]*proj_uv[idx2_2]) - (1.0/3.0)*proj_uv[idx1_3]*proj_uv[idx2_3];
455 proj_uvmn.push_back(temp);
456 }
457 }
458 }
459 }
460 return proj_uvmn;
461
462 }else
463 {
464 cout<<"rank>2: please add it by yourself!!!"<<endl;
465 exit(0);
466 }
467
468}
469double D0Topipipi0::fundecaymomentum(double mr2, double m1_2, double m2_2){
470 double mr = sqrt(mr2);
471 double poly = mr2*mr2 + m1_2*m1_2 + m2_2*m2_2 - 2*m1_2*mr2 -2*m2_2*mr2 -2*m1_2*m2_2;
472 double ret = sqrt(poly)/(2*mr);
473 if(poly<0)
474 //if(sqrt(m1_2) +sqrt(m2_2) > mr)
475 ret = 0.0f;
476 return ret;
477}
478
479complex<double> D0Topipipi0::breitwigner(double mx2, double mr, double wr)
480{
481 double output_x = 0;
482 double output_y = 0;
483
484 double mr2 = mr*mr;
485 double diff = mr2-mx2;
486 double denom = diff*diff + wr*wr*mr2;
487 if (wr<0) {
488 output_x = 0;
489 output_y = 0;
490 }
491 else if (wr<10) {
492 output_x = diff/denom;
493 output_y = wr*mr/denom;
494 }
495 else { /* phase space */
496 output_x = 1;
497 output_y = 0;
498 }
499 complex<double> output(output_x,output_y);
500 return output;
501
502}
503
504/* GS propagator */
505double D0Topipipi0::h(double m, double q){
506 double h = 2.0/math_pi*q/m*log((m+2.0*q)/(2.0*mass_Pion));
507 return h;
508}
509
510double D0Topipipi0::dh(double m0, double q0){
511 double dh = h(m0,q0)*(1.0/(8.0*q0*q0)-1.0/(2.0*m0*m0))+1.0/(2.0*math_pi*m0*m0);
512 return dh;
513}
514
515double D0Topipipi0::f(double m0, double sx, double q0, double q){
516 double m = sqrt(sx);
517 double f = m0*m0/(q0*q0*q0)*(q*q*(h(m,q)-h(m0,q0))+(m0*m0-sx)*q0*q0*dh(m0,q0));
518 return f;
519}
520
521double D0Topipipi0::d(double m0, double q0){
522 double d = 3.0/math_pi*mass_Pion*mass_Pion/(q0*q0)*log((m0+2.0*q0)/(2.0*mass_Pion)) + m0/(2.0*math_pi*q0) - (mass_Pion*mass_Pion*m0)/(math_pi*q0*q0*q0);
523 return d;
524}
525
526double D0Topipipi0::fundecaymomentum2(double mr2, double m1_2, double m2_2){
527 double mr = sqrt(mr2);
528 double poly = mr2*mr2 + m1_2*m1_2 + m2_2*m2_2 - 2*m1_2*mr2 -2*m2_2*mr2 -2*m1_2*m2_2;
529 double ret = poly/(4.0f*mr2);
530 if(poly<0)
531 ret = 0.0f;
532 return ret;
533}
534
535double D0Topipipi0::wid(double mass, double sa, double sb, double sc, double r, int l){
536 double widm = 1.0;
537 double sa0 = mass*mass;
538 double m = sqrt(sa);
539 double q = fundecaymomentum2(sa,sb,sc);
540 double q0 = fundecaymomentum2(sa0,sb,sc);
541 r = r/0.197321;
542 double z = q*r*r;
543 double z0 = q0*r*r;
544 double F = 0.0;
545 if(l == 0) F = 1.0;
546 if(l == 1) F = sqrt((1.0+z0)/(1.0+z));
547 if(l == 2) F = sqrt((9.0+3.0*z0+z0*z0)/(9.0+3.0*z+z*z));
548 if(l == 3) F = sqrt((225.0+45.0*z0+6.0*z0*z0+z0*z0*z0)/(225.0+45.0*z+6.0*z*z+z*z*z));
549 if(l == 4) F = sqrt((11025.0+1575.0*z0+135.0*z0*z0+10.0*z0*z0*z0+z0*z0*z0*z0)/(11025.0+1575.0*z+135.0*z*z+10.0*z*z*z+z*z*z*z));
550 double t = sqrt(q/q0);
551 //printf("sa0 = %f, sb = %f, sc = %f, q = %f, q0 = %0.15f, qq0 = %f, t = %f \n",sa0,sb,sc,q,q0,q/q0,t);
552 uint i=0;
553 for(i=0; i<(2*l+1); i++) {
554 widm *= t;
555 }
556 widm *= (mass/m*F*F);
557 return widm;
558}
559
560/* for rho0, use GS, rho0->pi+ pi-, only angular momentum 1*/
561complex<double> D0Topipipi0::GS(double mx2, double mr, double wr, double m1_2, double m2_2, double r, int l){
562
563 double mr2 = mr*mr;
564 double q = fundecaymomentum(mx2, m1_2, m2_2);
565 double q0 = fundecaymomentum(mr2, m1_2, m2_2);
566 double numer = 1.0+d(mr,q0)*wr/mr;
567 double denom_real = mr2-mx2+wr*f(mr,mx2,q0,q);
568 double denom_imag = mr*wr*wid(mr,mx2,m1_2,m2_2,r,l);//real-i*imag;
569
570 double denom = denom_real*denom_real+denom_imag*denom_imag;
571 double output_x = denom_real*numer/denom;
572 double output_y = denom_imag*numer/denom;
573
574 complex<double> output(output_x,output_y);
575 return output;
576}
577
578complex<double> D0Topipipi0::irho(double mr2, double m1_2, double m2_2){
579 double mr = sqrt(mr2);
580 double poly = mr2*mr2 + m1_2*m1_2 + m2_2*m2_2 - 2*m1_2*mr2 - 2*m2_2*mr2 - 2*m1_2*m2_2;
581 double ret_real = 0;
582 double ret_imag = 0;
583 if(poly>=0) {
584 ret_real = 2.0f*sqrt(poly)/(2.0f*mr2);
585 ret_imag = 0;
586 }
587 else {
588 ret_real = 0;
589 ret_imag = 2.0f*sqrt(-1.0f*poly)/(2.0f*mr2);
590 }
591 complex<double> ret(ret_real,ret_imag);
592 return ret;
593}
594
595complex<double> D0Topipipi0::Flatte(double mx2, double mr, double g1, double g2, double m1a, double m1b, double m2a, double m2b){
596
597 double mr2 = mr*mr;
598 double diff = mr2-mx2;
599 double g22 = g2*g1;
600 complex<double> ps1 = irho(mx2, m1a*m1a, m1b*m1b);
601 complex<double> ps2 = irho(mx2, m2a*m2a, m2b*m2b);
602 complex<double> iws = g1*ps1+g22*ps2; /*mass dependent width */
603
604 double denom_real = diff + iws.imag();
605 double denom_imag = iws.real();
606 double denom = denom_real*denom_real+denom_imag*denom_imag;
607
608 double output_x = denom_real/denom;
609 double output_y = denom_imag/denom;
610
611 complex<double> output(output_x,output_y);
612 return output;
613
614}
615
616complex<double> D0Topipipi0::RBW(double mx2, double mr, double wr, double m1_2, double m2_2, double r, int l){
617 double mx = sqrt(mx2);
618 double mr2 = mr*mr;
619 double denom_real = mr2-mx2;
620 double denom_imag = 0;
621 if(m1_2>0 && m2_2>0){
622 denom_imag = mr*wr*wid(mr,mx2,m1_2,m2_2,r,l);//real-i*imag;
623 }else{
624 denom_imag = mr*wr;
625 }
626
627 double denom = denom_real*denom_real+denom_imag*denom_imag;
628 double output_x = denom_real/denom;
629 double output_y = denom_imag/denom;
630
631 complex<double> output(output_x,output_y);
632 return output;
633}
634
635// PiPi-S wave K-Matrix
636double D0Topipipi0::rho22(double sc){
637 double rho[689] = { 3.70024e-18, 8.52763e-15, 1.87159e-13, 1.3311e-12, 5.61842e-12, 1.75224e-11, 4.48597e-11, 9.99162e-11, 2.00641e-10, 3.71995e-10,
638 6.47093e-10, 1.06886e-09, 1.69124e-09, 2.58031e-09, 3.8168e-09, 5.49601e-09, 7.72996e-09, 1.06509e-08, 1.44078e-08, 1.91741e-08,
639 2.51445e-08, 3.25345e-08, 4.15946e-08, 5.25949e-08, 6.58316e-08, 8.16443e-08, 1.00389e-07, 1.22455e-07, 1.48291e-07, 1.78348e-07,
640 2.1313e-07, 2.53192e-07, 2.99086e-07, 3.51462e-07, 4.10993e-07, 4.78349e-07, 5.54327e-07, 6.3972e-07, 7.35316e-07, 8.42099e-07,
641 9.61004e-07, 1.09295e-06, 1.2391e-06, 1.40051e-06, 1.57824e-06, 1.77367e-06, 1.98805e-06, 2.22257e-06, 2.47877e-06, 2.7581e-06,
642 3.06186e-06, 3.39182e-06, 3.74971e-06, 4.137e-06, 4.5555e-06, 5.00725e-06, 5.4939e-06, 6.01725e-06, 6.57992e-06, 7.18371e-06,
643 7.83044e-06, 8.52301e-06, 9.26342e-06, 1.00535e-05, 1.08967e-05, 1.17953e-05, 1.27514e-05, 1.37679e-05, 1.48482e-05, 1.59943e-05,
644 1.72088e-05, 1.84961e-05, 1.98586e-05, 2.12987e-05, 2.28207e-05, 2.44279e-05, 2.61228e-05, 2.79084e-05, 2.97906e-05, 3.17718e-05,
645 3.38544e-05, 3.60443e-05, 3.8345e-05, 4.07591e-05, 4.32903e-05, 4.59459e-05, 4.87285e-05, 5.16403e-05, 5.46887e-05, 5.7878e-05,
646 6.12111e-05, 6.46908e-05, 6.83274e-05, 7.21231e-05, 7.60817e-05, 8.0208e-05, 8.45102e-05, 8.89919e-05, 9.36544e-05, 9.85082e-05,
647 0.000103559, 0.000108812, 0.000114267, 0.000119938, 0.000125827, 0.00013194, 0.000138278, 0.000144857, 0.000151681, 0.000158752,
648 0.000166074, 0.000173663, 0.000181521, 0.000189652, 0.000198059, 0.000206761, 0.000215761, 0.000225063, 0.00023467, 0.000244599,
649 0.000254855, 0.00026544, 0.000276357, 0.000287629, 0.00029926, 0.000311253, 0.000323609, 0.000336351, 0.000349483, 0.000363009,
650 0.000376926, 0.000391264, 0.000406029, 0.000421225, 0.000436848, 0.000452921, 0.000469458, 0.000486461, 0.00050393, 0.00052187,
651 0.000540322, 0.000559278, 0.000578746, 0.00059872, 0.000619236, 0.0006403, 0.000661911, 0.000684074, 0.000706799, 0.000730127,
652 0.00075405, 0.000778569, 0.000803686, 0.000829443, 0.000855839, 0.000882879, 0.000910561, 0.000938898, 0.000967939, 0.000997674,
653 0.00102811, 0.00105923, 0.0010911, 0.0011237, 0.00115706, 0.00119117, 0.00122601, 0.00126168, 0.00129815, 0.00133543,
654 0.00137351, 0.00141242, 0.00145219, 0.00149283, 0.00153434, 0.0015767, 0.00161995, 0.00166415, 0.00170928, 0.00175534,
655 0.00180232, 0.00185028, 0.00189924, 0.00194919, 0.00200014, 0.00205207, 0.00210503, 0.0021591, 0.00221421, 0.0022704,
656 0.00232766, 0.00238602, 0.00244554, 0.00250619, 0.00256799, 0.0026309, 0.002695, 0.00276033, 0.00282689, 0.00289467,
657 0.00296367, 0.00303389, 0.00310543, 0.0031783, 0.00325244, 0.0033279, 0.0034046, 0.00348275, 0.00356229, 0.00364322,
658 0.00372555, 0.00380924, 0.00389438, 0.00398104, 0.00406914, 0.00415877, 0.00424985, 0.00434235, 0.00443651, 0.00453224,
659 0.00462954, 0.00472848, 0.00482894, 0.00493102, 0.00503483, 0.00514029, 0.00524749, 0.0053563, 0.00546675, 0.00557905,
660 0.0056931, 0.00580901, 0.0059267, 0.00604613, 0.00616735, 0.00629049, 0.00641557, 0.00654254, 0.00667142, 0.00680216,
661 0.00693472, 0.00706946, 0.00720621, 0.00734497, 0.0074858, 0.00762855, 0.00777338, 0.00792036, 0.00806957, 0.00822087,
662 0.00837426, 0.00852982, 0.0086875, 0.00884756, 0.00900991, 0.00917447, 0.00934137, 0.00951052, 0.00968194, 0.0098558,
663 0.010032, 0.0102108, 0.0103919, 0.0105754, 0.0107612, 0.0109496, 0.0111406, 0.0113343, 0.0115305, 0.0117293,
664 0.0119303, 0.0121343, 0.0123409, 0.0125502, 0.0127623, 0.0129771, 0.0131944, 0.0134145, 0.0136376, 0.0138636,
665 0.0140924, 0.0143241, 0.0145587, 0.0147959, 0.0150363, 0.0152797, 0.0155262, 0.0157758, 0.0160283, 0.0162838,
666 0.0165421, 0.016804, 0.0170691, 0.0173374, 0.0176087, 0.0178835, 0.0181612, 0.0184423, 0.0187269, 0.0190149,
667 0.0193063, 0.0196009, 0.0198991, 0.0202003, 0.0205052, 0.0208137, 0.0211259, 0.0214418, 0.0217611, 0.0220841,
668 0.0224105, 0.0227406, 0.0230746, 0.0234125, 0.0237542, 0.0240996, 0.0244486, 0.0248012, 0.025158, 0.0255188,
669 0.0258837, 0.0262527, 0.0266256, 0.0270025, 0.0273833, 0.027768, 0.0281572, 0.0285505, 0.0289483, 0.0293503,
670 0.0297564, 0.0301665, 0.0305808, 0.0309997, 0.0314231, 0.0318511, 0.0322835, 0.0327205, 0.0331616, 0.0336073,
671 0.0340576, 0.0345128, 0.0349727, 0.0354373, 0.0359066, 0.0363807, 0.0368589, 0.0373419, 0.0378302, 0.0383234,
672 0.0388218, 0.0393252, 0.0398336, 0.040347, 0.0408652, 0.041388, 0.0419165, 0.0424502, 0.0429893, 0.0435338,
673 0.0440833, 0.044638, 0.0451976, 0.0457627, 0.0463338, 0.0469103, 0.047492, 0.0480797, 0.0486729, 0.0492716,
674 0.0498757, 0.0504852, 0.0511009, 0.0517229, 0.0523503, 0.0529838, 0.0536231, 0.0542678, 0.054918, 0.0555743,
675 0.0562372, 0.0569065, 0.0575818, 0.0582634, 0.0589511, 0.0596454, 0.0603451, 0.061051, 0.0617635, 0.0624826,
676 0.0632084, 0.0639409, 0.06468, 0.0654254, 0.0661772, 0.0669346, 0.0676994, 0.0684714, 0.0692503, 0.0700354,
677 0.0708285, 0.0716277, 0.0724347, 0.0732479, 0.0740671, 0.0748947, 0.0757299, 0.0765715, 0.0774207, 0.0782771,
678 0.0791407, 0.0800119, 0.0808897, 0.0817743, 0.0826672, 0.0835684, 0.0844769, 0.0853938, 0.0863179, 0.0872493,
679 0.0881882, 0.0891349, 0.090089, 0.0910523, 0.0920236, 0.093002, 0.0939894, 0.094985, 0.0959887, 0.0970003,
680 0.0980191, 0.0990454, 0.100081, 0.101126, 0.10218, 0.103242, 0.104312, 0.105392, 0.10648, 0.107576,
681 0.10868, 0.109793, 0.110916, 0.112048, 0.113188, 0.114339, 0.115498, 0.116666, 0.117843, 0.119028,
682 0.120223, 0.121427, 0.122641, 0.123865, 0.125098, 0.126342, 0.127595, 0.128857, 0.130128, 0.131409,
683 0.132701, 0.134002, 0.135314, 0.136635, 0.137966, 0.139308, 0.14066, 0.142022, 0.143394, 0.144774,
684 0.146166, 0.14757, 0.148985, 0.15041, 0.151845, 0.153291, 0.154749, 0.156215, 0.157694, 0.159182,
685 0.160682, 0.162194, 0.163718, 0.165251, 0.166797, 0.168354, 0.169921, 0.1715, 0.17309, 0.17469,
686 0.176304, 0.177929, 0.179566, 0.181216, 0.182878, 0.184553, 0.186238, 0.187934, 0.189642, 0.191362,
687 0.193096, 0.194842, 0.196602, 0.198374, 0.200158, 0.201954, 0.203764, 0.205586, 0.207421, 0.209266,
688 0.211124, 0.212997, 0.214882, 0.216783, 0.218697, 0.220624, 0.222565, 0.224518, 0.226486, 0.228466,
689 0.230458, 0.232463, 0.234484, 0.23652, 0.238569, 0.240633, 0.242711, 0.244803, 0.246909, 0.249031,
690 0.251165, 0.253313, 0.255475, 0.257649, 0.259841, 0.262051, 0.264274, 0.266514, 0.268768, 0.271036,
691 0.273319, 0.275618, 0.277932, 0.280259, 0.282602, 0.28496, 0.287338, 0.28973, 0.292138, 0.294563,
692 0.297003, 0.299458, 0.30193, 0.304417, 0.306919, 0.309437, 0.311972, 0.314526, 0.317095, 0.319684,
693 0.322289, 0.324911, 0.327551, 0.330205, 0.332876, 0.335567, 0.338271, 0.340993, 0.343736, 0.346496,
694 0.349272, 0.352065, 0.354878, 0.35771, 0.360561, 0.363426, 0.366311, 0.369212, 0.372128, 0.375067,
695 0.378027, 0.381006, 0.384001, 0.387014, 0.39005, 0.393106, 0.396181, 0.399271, 0.402384, 0.405513,
696 0.408661, 0.41183, 0.41502, 0.418233, 0.421462, 0.424709, 0.42798, 0.43127, 0.434583, 0.437914,
697 0.441267, 0.444637, 0.448022, 0.451434, 0.454868, 0.458328, 0.461805, 0.465302, 0.468821, 0.472364,
698 0.475928, 0.47951, 0.483119, 0.486748, 0.490397, 0.494066, 0.497758, 0.501477, 0.505217, 0.508977,
699 0.512762, 0.516567, 0.520394, 0.524247, 0.528125, 0.532027, 0.535947, 0.53989, 0.543852, 0.547844,
700 0.551863, 0.555904, 0.559966, 0.56406, 0.568177, 0.572312, 0.576471, 0.580662, 0.584875, 0.58911,
701 0.593373, 0.597653, 0.601965, 0.606301, 0.610663, 0.615051, 0.619465, 0.623907, 0.62837, 0.632863,
702 0.637383, 0.641924, 0.646494, 0.651091, 0.655708, 0.660356, 0.665027, 0.669732, 0.674464, 0.679227,
703 0.684016, 0.688827, 0.693664, 0.698532, 0.703428, 0.708353, 0.713307, 0.718283, 0.72329, 0.728322,
704 0.733387, 0.738479, 0.743605, 0.748763, 0.753949, 0.759163, 0.764407, 0.769674, 0.774973, 0.780311,
705 0.78567, 0.791057, 0.796476, 0.801922, 0.8074, 0.812919, 0.818466, 0.824044 };
706
707 double m2 = 0.13957*0.13957;
708 double smin = (0.13957*4)*(0.13957*4);
709 double dh = 0.001;
710 int od = (sc - 0.312)/dh;
711 double sc_m = 0.312 + od*dh;
712 double rhouse = 0;
713 if(sc>=0.312 && sc<1){
714 rhouse = ((sc-sc_m)/dh)*(rho[od+1]-rho[od])+rho[od];
715 }else if(sc<0.312 && sc>=smin){
716 rhouse = ((sc - smin)/(0.312-smin))*rho[0];
717 }else if(sc>=1){
718 // rhouse = (1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc);
719 rhouse = sqrt(1-16*m2/sc);
720 }else{
721 rhouse = 0;
722 }
723 return rhouse;
724}
725
726complex<double> D0Topipipi0::rhoMTX(int i, int j, double s){
727
728 double rhoijx;
729 double rhoijy;
730 double mpi = 0.13957;
731 if(i==j && i==0 ){
732 double m2 = 0.13957*0.13957;
733 if((1-(4*m2)/s)>0){
734 rhoijx = sqrt(1.0f - (4*m2)/s);
735 rhoijy = 0;
736 }else{
737 rhoijy = sqrt((4*m2)/s - 1.0f);
738 rhoijx = 0;
739 }
740 }
741 if(i==j && i==1 ){
742 double m2 = 0.493677*0.493677;
743 if((1-(4*m2)/s)>0){
744 rhoijx = sqrt(1.0f - (4*m2)/s);
745 rhoijy = 0;
746 }else{
747 rhoijy = sqrt((4*m2)/s - 1.0f);
748 rhoijx = 0;
749 }
750 }
751 if(i==j && i==2){
752 rhoijx = rho22(s);
753 rhoijy = 0;
754 }
755 if(i==j && i==3){
756 double m2 = 0.547862*0.547862;
757 if((1-(4*m2)/s)>0){
758 rhoijx = sqrt(1.0f - (4*m2)/s);
759 rhoijy = 0;
760 }else{
761 rhoijy = sqrt((4*m2)/s - 1.0f);
762 rhoijx = 0;
763 }
764 }
765 if(i==j && i==4){
766 double m_1 = 0.547862;
767 double m_2 = 0.95778;
768 double mp2 = (m_1+m_2)*(m_1+m_2);
769 double mm2 = (m_1-m_2)*(m_1-m_2);
770 if((1-mp2/s)>0){
771 rhoijx = sqrt(1.0f - mp2/s);
772 rhoijy = 0;
773 }else{
774 rhoijy = sqrt(mp2/s - 1.0f);
775 rhoijx = 0;
776 }
777 }
778
779 if(i!=j){
780 rhoijx = 0;
781 rhoijy = 0;
782 }
783 complex<double> rhoij(rhoijx,rhoijy);
784 return rhoij;
785
786}
787
788/* Kij = (sum_a gigj/(ma^2-s) + fij(1-s0)/(s-s0))((s-0.5sampi2)(1-sa0)/(s-sa0)) */
789complex<double> D0Topipipi0::KMTX(int i, int j, double s){
790
791 double Kijx;
792 double Kijy;
793 double mpi = 0.13957;
794 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
795
796 double g1[5] = { 0.22889,-0.55377, 0.00000,-0.39899,-0.34639};
797 double g2[5] = { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503};
798 double g3[5] = { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681};
799 double g4[5] = { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984};
800 double g5[5] = { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358};
801
802 double f1[5] = { 0.23399, 0.15044,-0.20545, 0.32825, 0.35412};
803
804 double eps = 1e-11;
805
806 double down[5] = { 0,0,0,0,0};
807 double upreal[5] = { 0,0,0,0,0};
808 double upimag[5] = { 0,0,0,0,0};
809
810 for(int k=0; k<5; k++){
811
812 /* down[k] = (m[k]*m[k]-s)*(m[k]*m[k]-s)+eps*eps;
813 upreal[k] = (m[k]*m[k]-s)/down[k];
814 upimag[k] = -1.0f * eps/down[k]; */
815
816 double dm2 = m[k]*m[k]-s;
817 if(fabs(dm2)<eps && dm2<=0) dm2 = -eps;
818 if(fabs(dm2)<eps && dm2>0) dm2 = eps;
819 upreal[k] = 1.0f/dm2;
820 upimag[k] = 0;
821 }
822
823 double tmp1x = g1[i]*g1[j]*upreal[0] + g2[i]*g2[j]*upreal[1] + g3[i]*g3[j]*upreal[2] + g4[i]*g4[j]*upreal[3] + g5[i]*g5[j]*upreal[4];
824 double tmp1y = g1[i]*g1[j]*upimag[0] + g2[i]*g2[j]*upimag[1] + g3[i]*g3[j]*upimag[2] + g4[i]*g4[j]*upimag[3] + g5[i]*g5[j]*upimag[4];
825
826 double tmp2 = 0;
827 if(i==0){
828 tmp2 = f1[j]*(1+3.92637)/(s+3.92637);
829 }
830 if(j==0){
831 tmp2 = f1[i]*(1+3.92637)/(s+3.92637);
832 }
833 double tmp3 = (s-0.5*mpi*mpi)*(1+0.15)/(s+0.15);
834
835 Kijx = (tmp1x+tmp2)*tmp3;
836 Kijy = (tmp1y)*tmp3;
837
838 complex<double> Kij(Kijx,Kijy);
839 return Kij;
840}
841
842complex<double> D0Topipipi0::IMTX(int i, int j){
843
844 double Iijx;
845 double Iijy;
846 if(i==j){
847 Iijx = 1;
848 Iijy = 0;
849 }else{
850 Iijx = 0;
851 Iijy = 0;
852 }
853 complex<double> Iij(Iijx,Iijy);
854 return Iij;
855
856}
857
858/* F = I - i*K*rho */
859complex<double> D0Topipipi0::FMTX(double Kijx, double Kijy, double rhojjx, double rhojjy, int i, int j){
860
861 double Fijx;
862 double Fijy;
863
864 double tmpx = rhojjx*Kijx - rhojjy*Kijy;
865 double tmpy = rhojjx*Kijy + rhojjy*Kijx;
866
867 Fijx = IMTX(i,j).real() + tmpy;
868 Fijy = -tmpx;
869
870 complex<double> Fij(Fijx,Fijy);
871 return Fij;
872}
873
874/* inverse for Matrix (I - i*rho*K) using LUP */
875double D0Topipipi0::FINVMTX(double s, double *FINVx, double *FINVy){
876
877 int P[5] = { 0,1,2,3,4};
878
879 double Fx[5][5];
880 double Fy[5][5];
881
882 double Ux[5][5];
883 double Uy[5][5];
884 double Lx[5][5];
885 double Ly[5][5];
886
887 double UIx[5][5];
888 double UIy[5][5];
889 double LIx[5][5];
890 double LIy[5][5];
891
892 for(int k=0; k<5; k++){
893 double rhokkx = rhoMTX(k,k,s).real();
894 double rhokky = rhoMTX(k,k,s).imag();
895 Ux[k][k] = rhokkx;
896 Uy[k][k] = rhokky;
897 for(int l=k; l<5; l++){
898 double Kklx = KMTX(k,l,s).real();
899 double Kkly = KMTX(k,l,s).imag();
900 Lx[k][l] = Kklx;
901 Ly[k][l] = Kkly;
902 Lx[l][k] = Lx[k][l];
903 Ly[l][k] = Ly[k][l];
904 }
905 }
906
907 for(int k=0; k<5; k++){
908 for(int l=0; l<5; l++){
909 double Fklx = FMTX(Lx[k][l],Ly[k][l],Ux[l][l],Uy[l][l],k,l).real();
910 double Fkly = FMTX(Lx[k][l],Ly[k][l],Ux[l][l],Uy[l][l],k,l).imag();
911 Fx[k][l] = Fklx;
912 Fy[k][l] = Fkly;
913 }
914 }
915
916 for(int k=0; k<5; k++){
917 double tmprM = (Fx[k][k]*Fx[k][k]+Fy[k][k]*Fy[k][k]);
918 int tmpID = 0;
919 for(int l=k; l<5; l++){
920 double tmprF = (Fx[l][k]*Fx[l][k]+Fy[l][k]*Fy[l][k]);
921 if(tmprM<=tmprF){
922 tmprM = tmprF;
923 tmpID = l;
924 }
925 }
926
927 int tmpP = P[k];
928 P[k] = P[tmpID];
929 P[tmpID] = tmpP;
930
931 for(int l=0; l<5; l++){
932
933 double tmpFx = Fx[k][l];
934 double tmpFy = Fy[k][l];
935
936 Fx[k][l] = Fx[tmpID][l];
937 Fy[k][l] = Fy[tmpID][l];
938
939 Fx[tmpID][l] = tmpFx;
940 Fy[tmpID][l] = tmpFy;
941
942 }
943
944 for(int l=k+1; l<5; l++){
945 double rFkk = Fx[k][k]*Fx[k][k] + Fy[k][k]*Fy[k][k];
946 double Fxlk = Fx[l][k];
947 double Fylk = Fy[l][k];
948 double Fxkk = Fx[k][k];
949 double Fykk = Fy[k][k];
950 Fx[l][k] = (Fxlk*Fxkk + Fylk*Fykk)/rFkk;
951 Fy[l][k] = (Fylk*Fxkk - Fxlk*Fykk)/rFkk;
952 for(int m=k+1; m<5; m++){
953 Fx[l][m] = Fx[l][m] - (Fx[l][k]*Fx[k][m] - Fy[l][k]*Fy[k][m]);
954 Fy[l][m] = Fy[l][m] - (Fx[l][k]*Fy[k][m] + Fy[l][k]*Fx[k][m]);
955 }
956 }
957 }
958
959 for(int k=0; k<5; k++){
960 for(int l=0; l<5 ;l++){
961 if(k==l){
962 Lx[k][k] = 1;
963 Ly[k][k] = 0;
964 Ux[k][k] = Fx[k][k];
965 Uy[k][k] = Fy[k][k];
966 }
967 if(k>l){
968 Lx[k][l] = Fx[k][l];
969 Ly[k][l] = Fy[k][l];
970 Ux[k][l] = 0;
971 Uy[k][l] = 0;
972 }
973 if(k<l){
974 Ux[k][l] = Fx[k][l];
975 Uy[k][l] = Fy[k][l];
976 Lx[k][l] = 0;
977 Ly[k][l] = 0;
978 }
979 }
980 }
981
982 // calculate Inverse for L and U
983 for(int k=0; k<5; k++){
984
985 LIx[k][k] = 1;
986 LIy[k][k] = 0;
987
988 double rUkk = Ux[k][k]*Ux[k][k] + Uy[k][k]*Uy[k][k];
989 UIx[k][k] = Ux[k][k]/rUkk;
990 UIy[k][k] = -1.0f * Uy[k][k]/rUkk ;
991
992 for(int l=(k+1); l<5; l++){
993 LIx[k][l] = 0;
994 LIy[k][l] = 0;
995 UIx[l][k] = 0;
996 UIy[l][k] = 0;
997 }
998
999 for(int l=(k-1); l>=0; l--){ // U-1
1000 double sx = 0;
1001 double c_sx = 0;
1002 double sy = 0;
1003 double c_sy = 0;
1004 for(int m=l+1; m<=k; m++){
1005 sx = sx - c_sx;
1006 double sx_tmp = sx + Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k];
1007 c_sx = (sx_tmp - sx) - (Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k]);
1008 sx = sx_tmp;
1009
1010 sy = sy - c_sy;
1011 double sy_tmp = sy + Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k];
1012 c_sy = (sy_tmp - sy) - (Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k]);
1013 sy = sy_tmp;
1014 }
1015 UIx[l][k] = -1.0f * (UIx[l][l]*sx - UIy[l][l]*sy);
1016 UIy[l][k] = -1.0f * (UIy[l][l]*sx + UIx[l][l]*sy);
1017 }
1018
1019 for(int l=k+1; l<5; l++){ // L-1
1020 double sx = 0;
1021 double c_sx = 0;
1022 double sy = 0;
1023 double c_sy = 0;
1024 for(int m=k; m<l; m++){
1025 sx = sx - c_sx;
1026 double sx_tmp = sx + Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k];
1027 c_sx = (sx_tmp - sx) - (Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k]);
1028 sx = sx_tmp;
1029
1030 sy = sy - c_sy;
1031 double sy_tmp = sy + Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k];
1032 c_sy = (sy_tmp - sy) - (Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k]);
1033 sy = sy_tmp;
1034 }
1035 LIx[l][k] = -1.0f * sx;
1036 LIy[l][k] = -1.0f * sy;
1037 }
1038 }
1039
1040 for(int m=0; m<5; m++){
1041 double resX = 0;
1042 double c_resX = 0;
1043 double resY = 0;
1044 double c_resY = 0;
1045 for(int k=0; k<5; k++){
1046 for(int l=0; l<5; l++){
1047 double Plm = 0;
1048 if(P[l] == m) Plm = 1;
1049
1050 resX = resX - c_resX;
1051 double resX_tmp = resX + (UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm;
1052 c_resX = (resX_tmp - resX) - ((UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm);
1053 resX = resX_tmp;
1054
1055 resY = resY - c_resY;
1056 double resY_tmp = resY + (UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm;
1057 c_resY = (resY_tmp - resY) - ((UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm);
1058 resY = resY_tmp;
1059 }
1060 }
1061 FINVx[m] = resX;
1062 FINVy[m] = resY;
1063 }
1064
1065 return 1.0f;
1066
1067}
1068
1069complex<double> D0Topipipi0::PVTR(int ID, double s){
1070
1071 double VPix;
1072 double VPiy;
1073 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
1074
1075 double eps = 1e-11;
1076
1077 /* double down = (m[ID]*m[ID]-s)*(m[ID]*m[ID]-s)+eps*eps;
1078 double upreal = (m[ID]*m[ID]-s)/down;
1079 double upimag = -1.0f * eps/down; */
1080
1081 double dm2 = m[ID]*m[ID]-s;
1082
1083 if(fabs(dm2)<eps && dm2<=0) dm2 = -eps;
1084 if(fabs(dm2)<eps && dm2>0) dm2 = eps;
1085
1086 VPix = 1.0f/dm2;
1087 VPiy = 0;
1088
1089 complex<double> VPi(VPix,VPiy);
1090 return VPi;
1091}
1092
1093complex<double> D0Topipipi0::Fvector(double sa, double s0, int l){
1094
1095 double outputx = 0;
1096 double outputy = 0;
1097
1098 double FINVx[5] = {0,0,0,0,0};
1099 double FINVy[5] = {0,0,0,0,0};
1100
1101 double tmpFLAG = FINVMTX(sa,FINVx,FINVy);
1102
1103 if(l<5){
1104 double g[5][5] = {{ 0.22889,-0.55377, 0.00000,-0.39899,-0.34639},
1105 { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503},
1106 { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681},
1107 { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984},
1108 { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358}};
1109 double resx = 0;
1110 double c_resx = 0;
1111 double resy = 0;
1112 double c_resy = 0;
1113 double Plx = PVTR(l,sa).real();
1114 double Ply = PVTR(l,sa).imag();
1115 for(int j=0; j<5; j++){
1116 resx = resx - c_resx;
1117 double resx_tmp = resx + (FINVx[j]*g[l][j]*Plx - FINVy[j]*g[l][j]*Ply);
1118 c_resx = (resx_tmp - resx) - (FINVx[j]*g[l][j]*Plx - FINVy[j]*g[l][j]*Ply);
1119 resx = resx_tmp;
1120
1121 resy = resy - c_resy;
1122 double resy_tmp = resy + (FINVx[j]*g[l][j]*Ply + FINVy[j]*g[l][j]*Plx);
1123 c_resy = (resy_tmp - resy) - (FINVx[j]*g[l][j]*Ply + FINVy[j]*g[l][j]*Plx);
1124 resy = resy_tmp;
1125 }
1126 outputx = resx;
1127 outputy = resy;
1128 }else{
1129 int idx = l-5;
1130 double eps = 1e-11;
1131 double ds = sa-s0;
1132 if(fabs(ds)<eps && ds<=0) ds = -eps;
1133 if(fabs(ds)<eps && ds>0) ds = eps;
1134 double tmp = (1-s0)/ds;
1135 outputx = FINVx[idx]*tmp;
1136 outputy = FINVy[idx]*tmp;
1137 }
1138
1139 complex<double> output(outputx,outputy);
1140 return output;
1141}
1142
1143complex<double> D0Topipipi0::CalD0Amp(){
1144 return Amp(p4_Pip, p4_Pim, p4_Pi0);
1145}
1146complex<double> D0Topipipi0::CalDbAmp(){
1147
1148 vector<double> cpPip; cpPip.clear();
1149 vector<double> cpPim; cpPim.clear();
1150 vector<double> cpPi0; cpPi0.clear();
1151
1152 cpPip.push_back(-p4_Pim[0]); cpPim.push_back(-p4_Pip[0]); cpPi0.push_back(-p4_Pi0[0]);
1153 cpPip.push_back(-p4_Pim[1]); cpPim.push_back(-p4_Pip[1]); cpPi0.push_back(-p4_Pi0[1]);
1154 cpPip.push_back(-p4_Pim[2]); cpPim.push_back(-p4_Pip[2]); cpPi0.push_back(-p4_Pi0[2]);
1155 cpPip.push_back( p4_Pim[3]); cpPim.push_back( p4_Pip[3]); cpPi0.push_back( p4_Pi0[3]);
1156
1157 return (-1.0)*Amp(cpPip, cpPim, cpPi0);
1158}
1159
1160complex<double> D0Topipipi0::Amp(vector<double> Pip, vector<double> Pim, vector<double> Pi0)
1161{
1162
1163 if(fitpara.size()!=14) {
1164 cout<<"Error!!! number of para: "<<fitpara.size()<<endl;
1165 exit(0);
1166 }
1167
1168 vector<double> PipPim; PipPim.clear();
1169 vector<double> PipPi0; PipPi0.clear();
1170 vector<double> PimPi0; PimPi0.clear();
1171
1172 PipPim = sum_tensor(Pip, Pim);
1173 PipPi0 = sum_tensor(Pip, Pi0);
1174 PimPi0 = sum_tensor(Pim, Pi0);
1175
1176 vector<double> D0; D0.clear();
1177 D0 = sum_tensor(PipPim, Pi0);
1178
1179 double M2_PipPim = contract_11_0(PipPim, PipPim);
1180 double M2_PipPi0 = contract_11_0(PipPi0, PipPi0);
1181 double M2_PimPi0 = contract_11_0(PimPi0, PimPi0);
1182
1183 double M2_D0 = contract_11_0(D0, D0);
1184
1185 complex<double> GS_rho770_z = GS(M2_PipPim, m0_rho7700, w0_rho7700, m2_Pi, m2_Pi, rRes, 1);
1186 complex<double> GS_rho770_p = GS(M2_PipPi0, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1187 complex<double> GS_rho770_m = GS(M2_PimPi0, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1188
1189 complex<double> PiPiS_pm_0 = Fvector(M2_PipPim, s0_prod, 0);
1190 complex<double> PiPiS_pm_1 = Fvector(M2_PipPim, s0_prod, 1);
1191 complex<double> PiPiS_pm_2 = Fvector(M2_PipPim, s0_prod, 2);
1192 complex<double> PiPiS_pm_3 = Fvector(M2_PipPim, s0_prod, 3);
1193 complex<double> PiPiS_pm_4 = Fvector(M2_PipPim, s0_prod, 4);
1194 complex<double> PiPiS_pm_5 = Fvector(M2_PipPim, s0_prod, 5);
1195 complex<double> PiPiS_pm_6 = Fvector(M2_PipPim, s0_prod, 6);
1196 complex<double> PiPiS_pm_7 = Fvector(M2_PipPim, s0_prod, 7);
1197 complex<double> PiPiS_pm_8 = Fvector(M2_PipPim, s0_prod, 8);
1198 complex<double> PiPiS_pm_9 = Fvector(M2_PipPim, s0_prod, 9);
1199
1200 complex<double> RBW_f21270 = RBW(M2_PipPim, m0_f21270, w0_f21270, m2_Pi, m2_Pi, rRes, 2);
1201
1202 // X->PP Orbital
1203 vector<double> T1_PipPim; T1_PipPim.clear();
1204 vector<double> T1_PipPi0; T1_PipPi0.clear();
1205 vector<double> T1_PimPi0; T1_PimPi0.clear();
1206
1207 T1_PipPim = OrbitalTensors(PipPim, Pip, Pim, rRes, 1);
1208 T1_PipPi0 = OrbitalTensors(PipPi0, Pip, Pi0, rRes, 1);
1209 T1_PimPi0 = OrbitalTensors(PimPi0, Pim, Pi0, rRes, 1);
1210
1211 vector<double> T2_PipPim; T2_PipPim.clear();
1212
1213 T2_PipPim = OrbitalTensors(PipPim, Pip, Pim, rRes, 2);
1214
1215 // D->XP Orbital
1216 vector<double> T1_PipPimPi0; T1_PipPimPi0.clear();
1217 vector<double> T1_PipPi0Pim; T1_PipPi0Pim.clear();
1218 vector<double> T1_PimPi0Pip; T1_PimPi0Pip.clear();
1219
1220 T1_PipPimPi0 = OrbitalTensors(D0, PipPim, Pi0, rD, 1);
1221 T1_PipPi0Pim = OrbitalTensors(D0, PipPi0, Pim, rD, 1);
1222 T1_PimPi0Pip = OrbitalTensors(D0, PimPi0, Pip, rD, 1);
1223
1224 vector<double> T2_PipPimPi0; T2_PipPimPi0.clear();
1225
1226 T2_PipPimPi0 = OrbitalTensors(D0, PipPim, Pi0, rD, 2);
1227
1228 complex<double> amplitude(0,0);
1229
1230
1231 // D-> VP
1232 double SF_VpPm = contract_11_0(T1_PipPi0Pim, T1_PipPi0);
1233 amplitude += fitpara[0]*(SF_VpPm*GS_rho770_p);
1234
1235 double SF_VmPp = contract_11_0(T1_PimPi0Pip, T1_PimPi0);
1236 amplitude += fitpara[1]*(SF_VmPp*GS_rho770_m);
1237
1238 double SF_VzPz = contract_11_0(T1_PipPimPi0, T1_PipPim);
1239 amplitude += fitpara[2]*(SF_VzPz*GS_rho770_z);
1240
1241 // D-> SP
1242 amplitude += fitpara[3]*(PiPiS_pm_0);
1243 amplitude += fitpara[4]*(PiPiS_pm_1);
1244 amplitude += fitpara[5]*(PiPiS_pm_2);
1245 amplitude += fitpara[6]*(PiPiS_pm_3);
1246 amplitude += fitpara[7]*(PiPiS_pm_4);
1247 amplitude += fitpara[8]*(PiPiS_pm_5);
1248 amplitude += fitpara[9]*(PiPiS_pm_6);
1249 amplitude += fitpara[10]*(PiPiS_pm_7);
1250 amplitude += fitpara[11]*(PiPiS_pm_8);
1251 amplitude += fitpara[12]*(PiPiS_pm_9);
1252
1253 // D0 -> TP
1254 double SF_TzPz = contract_22_0(T2_PipPimPi0, T2_PipPim);
1255 amplitude += fitpara[13]*(SF_TzPz*RBW_f21270);
1256
1257 return amplitude;
1258
1259}
1260
1261int D0Topipipi0::CalAmp()
1262{
1263 m_AmpD0 = CalD0Amp();
1264 m_AmpDb = CalDbAmp();
1265 return 0;
1266}
1267
1268double D0Topipipi0::mag2(complex<double> x)
1269{
1270 double temp = x.real()*x.real() + x.imag()*x.imag();
1271 return temp;
1272}
1273
1274double D0Topipipi0::mag2_itf(complex<double> x, complex<double> y)
1275{
1276 double temp = 2*(x.real()*y.real() + x.imag()*y.imag());
1277 return temp;
1278}
1279
1280double D0Topipipi0::arg(complex<double> x)
1281{
1282 double temp = atan(x.imag()/x.real());
1283 if(x.real()<0) temp=temp+PI;
1284 return temp;
1285}
1286double D0Topipipi0::Get_strongPhase()
1287{
1288 double temp = arg(m_AmpD0) - arg(m_AmpDb);
1289 while (temp < -PI){
1290 temp += 2.0*PI;
1291 }
1292 while (temp > PI){
1293 temp -= 2.0*PI;
1294 }
1295 return temp;
1296}
1297
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
P
double P(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:87
mass
double mass
Definition CosmicGenerator.cxx:138
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
f1
TFile * f1
Definition DataBase/tau_mode.c:4
g1
TF1 * g1
Definition DataBase/tau_mode.c:66
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
ID
int ID[no]
Definition DocaEangle/fitDocaeangle.h:2
PI
#define PI
Definition EvtD0ToKKpi0.cc:41
eps
EvtTensor3C eps(const EvtVector3R &v)
Definition EvtTensor3C.cc:307
output
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in output
Definition FoamA.h:89
mpi
const double mpi
Definition Gam4pikp.cxx:47
s
XmlRpcServer s
Definition HelloServer.cpp:11
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
t
TTree * t
Definition binning.cxx:23
D0Topipipi0::~D0Topipipi0
virtual ~D0Topipipi0()
Definition D0Topipipi0.cxx:33
D0Topipipi0::Amp
complex< double > Amp(vector< double > Pip, vector< double > Pim, vector< double > Pi0)
Definition D0Topipipi0.cxx:1160
y
double y[1000]
Definition draw_charge_space_coarse_II.cxx:8
m2
double double * m2
Definition qcdloop1.h:75