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