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