EvtD0ToKKpi0.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of models developed at BES collaboration
//      based on the EvtGen framework.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/BesCopyright
//      Copyright (A) 2006      Ping Rong-Gang @IHEP
//
// Module:  EvtD0ToKKpi0.cc
//
// Description: Model provided by user, see BAM-628
//
// Modification history:
//
//    Liaoyuan Dong  Aug. 11, 2022   Module created
//
//------------------------------------------------------------------------
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtVector3R.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtHelSys.hh"
#include "EvtGenModels/EvtD0ToKKpi0.hh"
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <string>
#include <complex>
#include <vector>
#include <math.h>
using namespace std;
#define PI acos(-1)
 
EvtD0ToKKpi0::~EvtD0ToKKpi0() {}
 
void EvtD0ToKKpi0::getName(std::string& model_name){
  model_name="D0ToKKpi0";     
}
 
EvtDecayBase* EvtD0ToKKpi0::clone(){
  return new EvtD0ToKKpi0;
}
 
void EvtD0ToKKpi0::init(){
  checkNArg(2);
  checkNDaug(3);
  charm   = getArg(0);
  tagmode = getArg(1);
  //std::cout << "Initializing EvtD0ToKKpi0: charm= "<<charm<<" tagmode= "<<tagmode<<std::endl;
 
 
    g_uv.clear();
    for(int i=0; i<4; i++){
        for(int j=0; j<4; j++){
            if(i!=j){
                g_uv.push_back(0.0);
            }else if(i<3){
                g_uv.push_back(-1.0);
            }else if(i==3){
                g_uv.push_back(1.0);
            }
        }
    }
 
    epsilon_uvmn.clear();
    for(int i=0; i<4; i++){
        for(int j=0; j<4; j++){
            for(int k=0; k<4; k++){
                for(int l=0; l<4; l++){
                    if(i==j || i==k || i==l || j==k || j==l || k==l){
                        epsilon_uvmn.push_back(0.0);
                    }else{
                        if(i==0 && j==1 && k==2 && l==3) epsilon_uvmn.push_back(1.0);
                        if(i==0 && j==1 && k==3 && l==2) epsilon_uvmn.push_back(-1.0);
                        if(i==0 && j==2 && k==1 && l==3) epsilon_uvmn.push_back(-1.0);
                        if(i==0 && j==2 && k==3 && l==1) epsilon_uvmn.push_back(1.0);
                        if(i==0 && j==3 && k==1 && l==2) epsilon_uvmn.push_back(1.0);
                        if(i==0 && j==3 && k==2 && l==1) epsilon_uvmn.push_back(-1.0);
 
                        if(i==1 && j==0 && k==2 && l==3) epsilon_uvmn.push_back(-1.0);
                        if(i==1 && j==0 && k==3 && l==2) epsilon_uvmn.push_back(1.0);
                        if(i==1 && j==2 && k==0 && l==3) epsilon_uvmn.push_back(1.0);
                        if(i==1 && j==2 && k==3 && l==0) epsilon_uvmn.push_back(-1.0);
                        if(i==1 && j==3 && k==0 && l==2) epsilon_uvmn.push_back(-1.0);
                        if(i==1 && j==3 && k==2 && l==0) epsilon_uvmn.push_back(1.0);
 
                        if(i==2 && j==0 && k==1 && l==3) epsilon_uvmn.push_back(1.0);
                        if(i==2 && j==0 && k==3 && l==1) epsilon_uvmn.push_back(-1.0);
                        if(i==2 && j==1 && k==0 && l==3) epsilon_uvmn.push_back(-1.0);
                        if(i==2 && j==1 && k==3 && l==0) epsilon_uvmn.push_back(1.0);
                        if(i==2 && j==3 && k==0 && l==1) epsilon_uvmn.push_back(1.0);
                        if(i==2 && j==3 && k==1 && l==0) epsilon_uvmn.push_back(-1.0);
 
                        if(i==3 && j==0 && k==1 && l==2) epsilon_uvmn.push_back(-1.0);
                        if(i==3 && j==0 && k==2 && l==1) epsilon_uvmn.push_back(1.0);
                        if(i==3 && j==1 && k==0 && l==2) epsilon_uvmn.push_back(1.0);
                        if(i==3 && j==1 && k==2 && l==0) epsilon_uvmn.push_back(-1.0);
                        if(i==3 && j==2 && k==0 && l==1) epsilon_uvmn.push_back(-1.0);
                        if(i==3 && j==2 && k==1 && l==0) epsilon_uvmn.push_back(1.0);
 
                    }
                }
            }
        }
    }
 
  _nd = 3;
    math_pi = 3.1415926;
    mass_Pion = 0.13957;
 
    rRes = 3.0*0.197321;
    rD   = 5.0*0.197321;
    m_Pi  = mass_Pion;
    m2_Pi = m_Pi*m_Pi;
    m_Pi0  = 0.134977;
    m2_Pi0 = m_Pi0*m_Pi0;
    m_Ka  = 0.493677;
    m2_Ka  = m_Ka*m_Ka;
 
    m0_Kst892 = 0.89167;
    w0_Kst892 = 0.0514;
 
    m0_Kst1410 = 1.414;
    w0_Kst1410 = 0.232;
 
    m0_phi1020 = 1.019461;
    w0_phi1020 = 0.004249;
 
    m0_f0980  = 0.965;
    g1_f0980  = 0.165;
    g2_f0980  = 4.210;
 
    m0_K01430 = 1.441;
    w0_K01430 = 0.193;
 
    a_lass    = 0.113;
    r_lass    = -33.8;
 
    fitpara.clear();
    fitpara.push_back(complex<double>(100,0));
    fitpara.push_back(complex<double>(66.2031*cos(2.5703),66.2031*sin(2.5703)));
    fitpara.push_back(complex<double>(50.1991*cos(3.34488),50.1991*sin(3.34488)));
    fitpara.push_back(complex<double>(236.463*cos(-0.541955),236.463*sin(-0.541955)));
    fitpara.push_back(complex<double>(9389.32*cos(2.1910213),9389.32*sin(2.1910213)));
 
 
  return;
}
 
void EvtD0ToKKpi0::initProbMax(){
  setProbMax(1432000000); // 1431906201.4
}
 
void EvtD0ToKKpi0::decay( EvtParticle *p ){
/*
   double maxprob = 0.0;
   for(int ir=0;ir<=60000000;ir++){
      p->initializePhaseSpace(getNDaug(),getDaugs());
      for(int i=0; i<_nd; i++){
          _p4Lab[i]=p->getDaug(i)->getP4Lab();
          _p4CM[i]=p->getDaug(i)->getP4();
      }
      double Prob = AmplitudeSquare(charm, tagmode);
      if(Prob>maxprob) {
          maxprob=Prob;
          std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
      }
   }
   std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
*/
  p->initializePhaseSpace(getNDaug(),getDaugs());
  for(int i=0; i<_nd; i++){
          _p4Lab[i]=p->getDaug(i)->getP4Lab();
          _p4CM[i]=p->getDaug(i)->getP4();
  }
  double prob = AmplitudeSquare(charm, tagmode);
  setProb(prob);
  return;
}
 
void EvtD0ToKKpi0::setInput(double* kap, double* kam, double* pi0){
    p4_Kap.clear(); p4_Kam.clear(); p4_Pi0.clear();
    p4_Kap.push_back(kap[0]); p4_Kam.push_back(kam[0]); p4_Pi0.push_back(pi0[0]); 
    p4_Kap.push_back(kap[1]); p4_Kam.push_back(kam[1]); p4_Pi0.push_back(pi0[1]); 
    p4_Kap.push_back(kap[2]); p4_Kam.push_back(kam[2]); p4_Pi0.push_back(pi0[2]); 
    p4_Kap.push_back(kap[3]); p4_Kam.push_back(kam[3]); p4_Pi0.push_back(pi0[3]); 
}
 
vector<double> EvtD0ToKKpi0::sum_tensor(vector<double> pa, vector<double> pb) {
    if(pa.size()!=pb.size()){
        cout<<"error sum tensor"<<endl;
        exit(0);
    }
    vector<double> temp; temp.clear();
    for(int i=0; i<pa.size(); i++){
        double sum = pa[i] + pb[i];
        temp.push_back(sum);
    }
    return temp;
}
 
double EvtD0ToKKpi0::contract_11_0(vector<double> pa, vector<double> pb){
    if(pa.size()!=pb.size() || pa.size()!=4) {
        cout<<"error contract 11->0"<<endl;
        exit(0);
    }
    double temp = pa[3]*pb[3] - pa[0]*pb[0] - pa[1]*pb[1] - pa[2]*pb[2];
    return temp;
 
}
 
vector<double> EvtD0ToKKpi0::contract_21_1(vector<double> pa, vector<double> pb){
    if(pa.size()!=16 || pb.size()!=4) {
        cout<<"error contract 21->1"<<endl;
        exit(0);
    }
    vector<double> temp; temp.clear();
    for(int i=0; i<4; i++){
        double sum = 0;
        for(int j=0; j<4; j++){
            int idx   = i*4+j;
            sum += pa[idx]*pb[j]*g_uv[4*j+j];
        }
        temp.push_back(sum);
    }
    return temp;
 
}
 
double EvtD0ToKKpi0::contract_22_0(vector<double> pa, vector<double> pb){
    if(pa.size()!=pb.size() || pa.size()!=16) {
        cout<<"error contract 22->0"<<endl;
        exit(0);
    }
    double temp = 0;
    for(int i=0; i<4; i++){
        for(int j=0; j<4; j++){
            int idx    = i*4+j;
            temp += pa[idx]*pb[idx]*g_uv[4*i+i]*g_uv[4*j+j];
        }
    }
    return temp;
 
}
 
vector<double> EvtD0ToKKpi0::contract_31_2(vector<double> pa, vector<double> pb){
    if(pa.size()!=64 || pb.size()!=4) {
        cout<<"error contract 31->2"<<endl;
        exit(0);
    }
    vector<double> temp; temp.clear();
    for(int i=0; i<16; i++){
        double sum = 0;
        for(int j=0; j<4; j++){
            int idx   = i*4+j;
            sum += pa[idx]*pb[j]*g_uv[4*j+j];
        }
        temp.push_back(sum);
    }
    return temp;
 
}
 
vector<double> EvtD0ToKKpi0::contract_41_3(vector<double> pa, vector<double> pb){
    if(pa.size()!=256|| pb.size()!=4) {
        cout<<"error contract 41->3"<<endl;
        exit(0);
    }
    vector<double> temp; temp.clear();
    for(int i=0; i<64; i++){
        double sum = 0;
        for(int j=0; j<4; j++){
            int idx   = i*4+j;
            sum += pa[idx]*pb[j]*g_uv[4*j+j];
        }
        temp.push_back(sum);
    }
    return temp;
 
}
 
vector<double> EvtD0ToKKpi0::contract_42_2(vector<double> pa, vector<double> pb){
    if(pa.size()!=256|| pb.size()!=16) {
        cout<<"error contract 42->2"<<endl;
        exit(0);
    }
    vector<double> temp; temp.clear();
    for(int i=0; i<16; i++){
        double sum = 0;
        for(int j=0; j<4; j++){
            for(int k=0; k<4; k++){
                int idxa = i*16+j*4+k;
                int idxb = j*4+k;
                sum += pa[idxa] * pb[idxb] * g_uv[4*j+j] * g_uv[4*k+k];
            }
        }
        temp.push_back(sum);
    }
 
    return temp;
 
}
 
vector<double> EvtD0ToKKpi0::contract_22_2(vector<double> pa, vector<double> pb){
    if(pa.size()!=16|| pb.size()!=16) {
        cout<<"error contract 42->2"<<endl;
        exit(0);
    }
    vector<double> temp; temp.clear();
    for(int i=0; i<4; i++){
        for(int j=0; j<4; j++){
            double sum = 0;
            for(int k=0; k<4; k++){
                int idxa = i*4+k;
                int idxb = j*4+k;
                sum += pa[idxa] * pb[idxb] * g_uv[4*k+k];
            }
            temp.push_back(sum);
        }
    }
 
    return temp;
 
}
 
/*
   vector<double> EvtD0ToKKpi0::contract_11_2(vector<double> pa, vector<double> pb){
   if(pa.size()!=pb.size() || pa.size()!=4) {
   cout<<"error contract 11->2"<<endl;
   exit(0);
   }
   vector<double> temp; temp.clear();
   for(int i=0; i<4; i++){
   for(int j=0; j<4; j++){
   double prod = pa[i]*pb[j];
   temp.push_back(prod);
   }
   }
   return temp;
   }
 
   vector<double> EvtD0ToKKpi0::contract_22_4(vector<double> pa, vector<double> pb){
   if(pa.size()!=pb.size() || pa.size()!=16) {
   cout<<"error contract 22->4"<<endl;
   exit(0);
   }
   vector<double> temp; temp.clear();
   for(int i=0; i<16; i++){
   for(int j=0; j<16; j++){
   double prod = pa[i]*pb[j];
   temp.push_back(prod);
   }
   }
   return temp;
   }
   */
 
//OrbitalTensors
vector<double> EvtD0ToKKpi0::OrbitalTensors(vector<double> pa, vector<double> pb, vector<double> pc, double r, int rank)
{
 
    if(pa.size()!=4 || pb.size()!=4 || pc.size()!=4) {
        cout<<"Error: pa, pb, pc"<<endl;
        exit(0);
    }
    if(rank<0){ 
        cout<<"Error: L<0 !!!"<<endl;
        exit(0);
    }
 
    // relative momentum
    vector<double> mr; mr.clear();
 
    for(int i=0; i<4; i++){
        double temp = pb[i] - pc[i];
        mr.push_back(temp);
    }
 
    // "Masses" of particles
    double msa = contract_11_0(pa, pa);
    double msb = contract_11_0(pb, pb);
    double msc = contract_11_0(pc, pc);
 
    // Q^2_abc
    double top = msa + msb - msc;
    double Q2abc = top*top/(4.0*msa) - msb;
 
    // Barrier factors
    double Q_0 = 0.197321f/r;
    double Q_02 = Q_0*Q_0;
    double Q_04 = Q_02*Q_02;
    //  double Q_06 = Q_04*Q_02;
    //  double Q_08 = Q_04*Q_04;
 
    double Q4abc = Q2abc*Q2abc;
    //  double Q6abc = Q4abc*Q2abc;
    //  double Q8abc = Q4abc*Q4abc;
 
    double mB1 = sqrt(2.0f/(Q2abc + Q_02));
    double mB2 = sqrt(13.0f/(Q4abc + (3.0f*Q_02)*Q2abc + 9.0f*Q_04));
    //  mB3 = &sqrt(277.0f/(Q6abc + (6.0f*Q_02)*Q4abc + (45.0f*Q_04)*Q2abc + 225.0f*Q_06));
    //  mB4 = &sqrt(12746.0f/(Q8abc + (10.0f*Q_02)*Q6abc + (135.0f*Q_04)*Q4abc + (1575.0f*Q_06)*Q2abc + 11025.0f*Q_08));
 
    // Projection Operator 2-rank
    vector<double> proj_uv; proj_uv.clear();
    for(int i=0; i<4; i++){
        for(int j=0; j<4; j++){
            int idx = i*4 + j;
            double temp = -g_uv[idx] + pa[i]*pa[j]/msa;
            proj_uv.push_back(temp);
        }
    }
 
    // Orbital Tensors
    if(rank==0)
    {
        vector<double> t; t.clear();
        t.push_back(1.0);
        return t;
 
    }else if(rank<3)
    {
        vector<double> t_u; t_u.clear();
        vector<double> Bt_u; Bt_u.clear();
        for(int i=0; i<4; i++){
            double temp = 0;
            for(int j=0; j<4; j++){
                int idx = i*4 + j;
                temp += -proj_uv[idx]*mr[j]*g_uv[j*4+j];
            }
            t_u.push_back(temp);
            Bt_u.push_back(temp*mB1);
        }
        if(rank==1) return Bt_u;
 
        double t_u2 = contract_11_0(t_u,t_u);
 
        vector<double> Bt_uv; Bt_uv.clear();
        for(int i=0; i<4; i++){
            for(int j=0; j<4; j++){
                int idx = 4*i + j;
                double temp = t_u[i]*t_u[j] + (1.0/3.0)*proj_uv[idx]*t_u2;
                Bt_uv.push_back(temp*mB2);
            }
        }
        if(rank==2) return Bt_uv;
 
    }else
    {
        cout<<"rank>2: please add it by yourself!!!"<<endl;
        exit(0);
    }
 
}
 
// projection Tensor
vector<double> EvtD0ToKKpi0::ProjectionTensors(vector<double> pa, int rank)
{
 
    if(pa.size()!=4) {
        cout<<"Error: pa"<<endl;
        exit(0);
    }
    if(rank<0){ 
        cout<<"Error: L<0 !!!"<<endl;
        exit(0);
    }
 
    double msa = contract_11_0(pa, pa);
 
    // Projection Operator 2-rank
    vector<double> proj_uv; proj_uv.clear();
    for(int i=0; i<4; i++){
        for(int j=0; j<4; j++){
            int idx = i*4 + j;
            double temp = -g_uv[idx] + pa[i]*pa[j]/msa;
            proj_uv.push_back(temp);
        }
    }
 
    // Orbital Tensors
    if(rank==0)
    {
        vector<double> t; t.clear();
        t.push_back(1.0);
        return t;
 
    }else if(rank==1)
    {
        return  proj_uv;
    }else if(rank==2)
    {
        vector<double> proj_uvmn; proj_uvmn.clear();
        for(int i=0; i<4; i++){
            for(int j=0; j<4; j++){
                for(int k=0; k<4; k++){
                    for(int l=0; l<4; l++){
 
                        int idx1_1 = 4*i + k;
                        int idx1_2 = 4*i + l;
                        int idx1_3 = 4*i + j;
 
                        int idx2_1 = 4*j + l;
                        int idx2_2 = 4*j + k;
                        int idx2_3 = 4*k + l;
 
                        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];
                        proj_uvmn.push_back(temp);
                    }
                }
            }
        }
        return proj_uvmn;
 
    }else
    {
        cout<<"rank>2: please add it by yourself!!!"<<endl;
        exit(0);
    }
 
}
double EvtD0ToKKpi0::fundecaymomentum(double mr2, double m1_2, double m2_2){
    double mr = sqrt(mr2);
    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;
    double ret = sqrt(poly)/(2*mr);
    if(poly<0)
        //if(sqrt(m1_2) +sqrt(m2_2) > mr)
        ret = 0.0f;
    return ret;
}
 
complex<double> EvtD0ToKKpi0::breitwigner(double mx2, double mr, double wr)
{   
    double output_x = 0;
    double output_y = 0;
 
    double mr2 = mr*mr;
    double diff = mr2-mx2;
    double denom = diff*diff + wr*wr*mr2;
    if (wr<0) {
        output_x = 0;
        output_y = 0;
    }
    else if (wr<10) {
        output_x = diff/denom;
        output_y = wr*mr/denom;
    }
    else { /* phase space */
        output_x = 1;
        output_y = 0;
    }
    complex<double>  output(output_x,output_y);
    return output;
 
}
 
/* GS propagator */
double EvtD0ToKKpi0::h(double m, double q){
    double h = 2.0/math_pi*q/m*log((m+2.0*q)/(2.0*mass_Pion));
    return h;
} 
 
double EvtD0ToKKpi0::dh(double m0, double q0){
    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);
    return dh;
} 
 
double EvtD0ToKKpi0::f(double m0, double sx, double q0, double q){
    double m = sqrt(sx);
    double f = m0*m0/(q0*q0*q0)*(q*q*(h(m,q)-h(m0,q0))+(m0*m0-sx)*q0*q0*dh(m0,q0));
    return f;
}
 
double EvtD0ToKKpi0::d(double m0, double q0){
    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);
    return d;
}
 
double EvtD0ToKKpi0::fundecaymomentum2(double mr2, double m1_2, double m2_2){
    double mr = sqrt(mr2);
    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;
    double ret = poly/(4.0f*mr2);
    if(poly<0)
        ret = 0.0f;
    return ret;
}
 
double EvtD0ToKKpi0::wid(double mass, double sa, double sb, double sc, double r, int l){
    double widm = 1.0;
    double sa0 = mass*mass;
    double m = sqrt(sa);
    double q = fundecaymomentum2(sa,sb,sc);
    double q0 = fundecaymomentum2(sa0,sb,sc);
    r = r/0.197321;
    double z = q*r*r;
    double z0 = q0*r*r;
    double F = 0.0;
    if(l == 0) F = 1.0;
    if(l == 1) F = sqrt((1.0+z0)/(1.0+z));
    if(l == 2) F = sqrt((9.0+3.0*z0+z0*z0)/(9.0+3.0*z+z*z));
    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));
    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));
    double t = sqrt(q/q0);
    //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);
    uint i=0;
    for(i=0; i<(2*l+1); i++) {
        widm *= t;
    }
    widm *= (mass/m*F*F);
    return widm; 
}
 
/* for rho0, use GS, rho0->pi+ pi-, only angular momentum 1*/
complex<double> EvtD0ToKKpi0::GS(double mx2, double mr, double wr, double m1_2, double m2_2, double r, int l){
 
    double mr2 = mr*mr;
    double q = fundecaymomentum(mx2, m1_2, m2_2);
    double q0 = fundecaymomentum(mr2, m1_2, m2_2);
    double numer = 1.0+d(mr,q0)*wr/mr;
    double denom_real = mr2-mx2+wr*f(mr,mx2,q0,q);
    double denom_imag = mr*wr*wid(mr,mx2,m1_2,m2_2,r,l);//real-i*imag;
 
    double denom = denom_real*denom_real+denom_imag*denom_imag;
    double output_x = denom_real*numer/denom;
    double output_y = denom_imag*numer/denom;
 
    complex<double>  output(output_x,output_y);
    return output;
}
 
complex<double> EvtD0ToKKpi0::irho(double mr2, double m1_2, double m2_2){
        double mr = sqrt(mr2);
    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;
        double ret_real = 0;
        double ret_imag = 0;
        if(poly>=0) {
            ret_real = 2.0f*sqrt(poly)/(2.0f*mr2);
            ret_imag = 0;
        }
        else {
            ret_real = 0;
            ret_imag = 2.0f*sqrt(-1.0f*poly)/(2.0f*mr2);
        }
    complex<double>  ret(ret_real,ret_imag);
        return ret;
}
 
complex<double> EvtD0ToKKpi0::Flatte(double mx2, double mr, double g1, double g2, double m1a, double m1b, double m2a, double m2b){
 
    double mr2 = mr*mr;
    double diff = mr2-mx2;
    double g22 = g2*g1;
    complex<double> ps1 = irho(mx2, m1a*m1a, m1b*m1b);
    complex<double> ps2 = irho(mx2, m2a*m2a, m2b*m2b);
    complex<double> iws = g1*ps1+g22*ps2; /*mass dependent width */
 
    double denom_real = diff + iws.imag();
    double denom_imag = iws.real();
    double denom = denom_real*denom_real+denom_imag*denom_imag;
 
    double output_x = denom_real/denom;
    double output_y = denom_imag/denom;
  
    complex<double>  output(output_x,output_y);
    return output;
 
}
 
complex<double> EvtD0ToKKpi0::RBW(double mx2, double mr, double wr, double m1_2, double m2_2, double r, int l){
    double mx  = sqrt(mx2);
    double mr2 = mr*mr;
    double denom_real = mr2-mx2;
    double denom_imag = 0;
    if(m1_2>0 && m2_2>0){
        denom_imag = mr*wr*wid(mr,mx2,m1_2,m2_2,r,l);//real-i*imag;
    }else{
        denom_imag = mr*wr;
    }
 
    double denom = denom_real*denom_real+denom_imag*denom_imag;
    double output_x = denom_real/denom;
    double output_y = denom_imag/denom;
 
    complex<double>  output(output_x,output_y);
    return output;
}
 
complex<double> EvtD0ToKKpi0::KPiSLass(double mx2, double mr, double wr, double a, double r, double R, double phi_R, double phi_F){
    
    double mx  = sqrt(mx2);
    double mr2 = mr*mr;
    double q2abc  = fundecaymomentum2(mx2, m2_Ka, m2_Pi0);
    double q02abc = fundecaymomentum2(mr2, m2_Ka, m2_Pi0);
    
    double qr2 = q2abc/q02abc;
    if(qr2<0) qr2 = 0;
    double wid = wr*sqrt(qr2)*(mr/mx);
 
    double alphaR = atan(mr*wid/(mr2 - mx2));
    if((mr2 - mx2)<0) alphaR += math_pi;
 
    double alphaF = atan(2.0*a*sqrt(q2abc)/(2.0 + a*r*q2abc));
    if((2.0*a*sqrt(q2abc)/(2.0 + a*r*q2abc))<0) alphaF += math_pi;
 
    double delta_R = phi_R + alphaR;
    double delta_F = phi_F + alphaF;
 
    double output_x = R*sin(delta_R)*cos(delta_R+2.0*delta_F) + sin(delta_F)*cos(delta_F);
    double output_y = R*sin(delta_R)*sin(delta_R+2.0*delta_F) + sin(delta_F)*sin(delta_F);
 
    complex<double>  output(output_x,output_y);
    return output;
}
 
// PiPi-S wave K-Matrix
double EvtD0ToKKpi0::rho22(double sc){
    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, 
        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, 
        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, 
        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, 
        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, 
        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, 
        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, 
        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, 
        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, 
        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, 
        0.000103559, 0.000108812, 0.000114267, 0.000119938, 0.000125827, 0.00013194, 0.000138278, 0.000144857, 0.000151681, 0.000158752, 
        0.000166074, 0.000173663, 0.000181521, 0.000189652, 0.000198059, 0.000206761, 0.000215761, 0.000225063, 0.00023467, 0.000244599, 
        0.000254855, 0.00026544, 0.000276357, 0.000287629, 0.00029926, 0.000311253, 0.000323609, 0.000336351, 0.000349483, 0.000363009, 
        0.000376926, 0.000391264, 0.000406029, 0.000421225, 0.000436848, 0.000452921, 0.000469458, 0.000486461, 0.00050393, 0.00052187, 
        0.000540322, 0.000559278, 0.000578746, 0.00059872, 0.000619236, 0.0006403, 0.000661911, 0.000684074, 0.000706799, 0.000730127, 
        0.00075405, 0.000778569, 0.000803686, 0.000829443, 0.000855839, 0.000882879, 0.000910561, 0.000938898, 0.000967939, 0.000997674, 
        0.00102811, 0.00105923, 0.0010911, 0.0011237, 0.00115706, 0.00119117, 0.00122601, 0.00126168, 0.00129815, 0.00133543, 
        0.00137351, 0.00141242, 0.00145219, 0.00149283, 0.00153434, 0.0015767, 0.00161995, 0.00166415, 0.00170928, 0.00175534, 
        0.00180232, 0.00185028, 0.00189924, 0.00194919, 0.00200014, 0.00205207, 0.00210503, 0.0021591, 0.00221421, 0.0022704, 
        0.00232766, 0.00238602, 0.00244554, 0.00250619, 0.00256799, 0.0026309, 0.002695, 0.00276033, 0.00282689, 0.00289467, 
        0.00296367, 0.00303389, 0.00310543, 0.0031783, 0.00325244, 0.0033279, 0.0034046, 0.00348275, 0.00356229, 0.00364322, 
        0.00372555, 0.00380924, 0.00389438, 0.00398104, 0.00406914, 0.00415877, 0.00424985, 0.00434235, 0.00443651, 0.00453224, 
        0.00462954, 0.00472848, 0.00482894, 0.00493102, 0.00503483, 0.00514029, 0.00524749, 0.0053563, 0.00546675, 0.00557905, 
        0.0056931, 0.00580901, 0.0059267, 0.00604613, 0.00616735, 0.00629049, 0.00641557, 0.00654254, 0.00667142, 0.00680216, 
        0.00693472, 0.00706946, 0.00720621, 0.00734497, 0.0074858, 0.00762855, 0.00777338, 0.00792036, 0.00806957, 0.00822087, 
        0.00837426, 0.00852982, 0.0086875, 0.00884756, 0.00900991, 0.00917447, 0.00934137, 0.00951052, 0.00968194, 0.0098558, 
        0.010032, 0.0102108, 0.0103919, 0.0105754, 0.0107612, 0.0109496, 0.0111406, 0.0113343, 0.0115305, 0.0117293, 
        0.0119303, 0.0121343, 0.0123409, 0.0125502, 0.0127623, 0.0129771, 0.0131944, 0.0134145, 0.0136376, 0.0138636, 
        0.0140924, 0.0143241, 0.0145587, 0.0147959, 0.0150363, 0.0152797, 0.0155262, 0.0157758, 0.0160283, 0.0162838, 
        0.0165421, 0.016804, 0.0170691, 0.0173374, 0.0176087, 0.0178835, 0.0181612, 0.0184423, 0.0187269, 0.0190149, 
        0.0193063, 0.0196009, 0.0198991, 0.0202003, 0.0205052, 0.0208137, 0.0211259, 0.0214418, 0.0217611, 0.0220841, 
        0.0224105, 0.0227406, 0.0230746, 0.0234125, 0.0237542, 0.0240996, 0.0244486, 0.0248012, 0.025158, 0.0255188, 
        0.0258837, 0.0262527, 0.0266256, 0.0270025, 0.0273833, 0.027768, 0.0281572, 0.0285505, 0.0289483, 0.0293503, 
        0.0297564, 0.0301665, 0.0305808, 0.0309997, 0.0314231, 0.0318511, 0.0322835, 0.0327205, 0.0331616, 0.0336073, 
        0.0340576, 0.0345128, 0.0349727, 0.0354373, 0.0359066, 0.0363807, 0.0368589, 0.0373419, 0.0378302, 0.0383234, 
        0.0388218, 0.0393252, 0.0398336, 0.040347, 0.0408652, 0.041388, 0.0419165, 0.0424502, 0.0429893, 0.0435338, 
        0.0440833, 0.044638, 0.0451976, 0.0457627, 0.0463338, 0.0469103, 0.047492, 0.0480797, 0.0486729, 0.0492716, 
        0.0498757, 0.0504852, 0.0511009, 0.0517229, 0.0523503, 0.0529838, 0.0536231, 0.0542678, 0.054918, 0.0555743, 
        0.0562372, 0.0569065, 0.0575818, 0.0582634, 0.0589511, 0.0596454, 0.0603451, 0.061051, 0.0617635, 0.0624826, 
        0.0632084, 0.0639409, 0.06468, 0.0654254, 0.0661772, 0.0669346, 0.0676994, 0.0684714, 0.0692503, 0.0700354, 
        0.0708285, 0.0716277, 0.0724347, 0.0732479, 0.0740671, 0.0748947, 0.0757299, 0.0765715, 0.0774207, 0.0782771, 
        0.0791407, 0.0800119, 0.0808897, 0.0817743, 0.0826672, 0.0835684, 0.0844769, 0.0853938, 0.0863179, 0.0872493, 
        0.0881882, 0.0891349, 0.090089, 0.0910523, 0.0920236, 0.093002, 0.0939894, 0.094985, 0.0959887, 0.0970003, 
        0.0980191, 0.0990454, 0.100081, 0.101126, 0.10218, 0.103242, 0.104312, 0.105392, 0.10648, 0.107576, 
        0.10868, 0.109793, 0.110916, 0.112048, 0.113188, 0.114339, 0.115498, 0.116666, 0.117843, 0.119028, 
        0.120223, 0.121427, 0.122641, 0.123865, 0.125098, 0.126342, 0.127595, 0.128857, 0.130128, 0.131409, 
        0.132701, 0.134002, 0.135314, 0.136635, 0.137966, 0.139308, 0.14066, 0.142022, 0.143394, 0.144774, 
        0.146166, 0.14757, 0.148985, 0.15041, 0.151845, 0.153291, 0.154749, 0.156215, 0.157694, 0.159182, 
        0.160682, 0.162194, 0.163718, 0.165251, 0.166797, 0.168354, 0.169921, 0.1715, 0.17309, 0.17469, 
        0.176304, 0.177929, 0.179566, 0.181216, 0.182878, 0.184553, 0.186238, 0.187934, 0.189642, 0.191362, 
        0.193096, 0.194842, 0.196602, 0.198374, 0.200158, 0.201954, 0.203764, 0.205586, 0.207421, 0.209266, 
        0.211124, 0.212997, 0.214882, 0.216783, 0.218697, 0.220624, 0.222565, 0.224518, 0.226486, 0.228466, 
        0.230458, 0.232463, 0.234484, 0.23652, 0.238569, 0.240633, 0.242711, 0.244803, 0.246909, 0.249031, 
        0.251165, 0.253313, 0.255475, 0.257649, 0.259841, 0.262051, 0.264274, 0.266514, 0.268768, 0.271036, 
        0.273319, 0.275618, 0.277932, 0.280259, 0.282602, 0.28496, 0.287338, 0.28973, 0.292138, 0.294563, 
        0.297003, 0.299458, 0.30193, 0.304417, 0.306919, 0.309437, 0.311972, 0.314526, 0.317095, 0.319684, 
        0.322289, 0.324911, 0.327551, 0.330205, 0.332876, 0.335567, 0.338271, 0.340993, 0.343736, 0.346496, 
        0.349272, 0.352065, 0.354878, 0.35771, 0.360561, 0.363426, 0.366311, 0.369212, 0.372128, 0.375067, 
        0.378027, 0.381006, 0.384001, 0.387014, 0.39005, 0.393106, 0.396181, 0.399271, 0.402384, 0.405513, 
        0.408661, 0.41183, 0.41502, 0.418233, 0.421462, 0.424709, 0.42798, 0.43127, 0.434583, 0.437914, 
        0.441267, 0.444637, 0.448022, 0.451434, 0.454868, 0.458328, 0.461805, 0.465302, 0.468821, 0.472364, 
        0.475928, 0.47951, 0.483119, 0.486748, 0.490397, 0.494066, 0.497758, 0.501477, 0.505217, 0.508977, 
        0.512762, 0.516567, 0.520394, 0.524247, 0.528125, 0.532027, 0.535947, 0.53989, 0.543852, 0.547844, 
        0.551863, 0.555904, 0.559966, 0.56406, 0.568177, 0.572312, 0.576471, 0.580662, 0.584875, 0.58911, 
        0.593373, 0.597653, 0.601965, 0.606301, 0.610663, 0.615051, 0.619465, 0.623907, 0.62837, 0.632863, 
        0.637383, 0.641924, 0.646494, 0.651091, 0.655708, 0.660356, 0.665027, 0.669732, 0.674464, 0.679227, 
        0.684016, 0.688827, 0.693664, 0.698532, 0.703428, 0.708353, 0.713307, 0.718283, 0.72329, 0.728322, 
        0.733387, 0.738479, 0.743605, 0.748763, 0.753949, 0.759163, 0.764407, 0.769674, 0.774973, 0.780311, 
        0.78567, 0.791057, 0.796476, 0.801922, 0.8074, 0.812919, 0.818466, 0.824044 };
 
    double m2 = 0.13957*0.13957;
    double smin = (0.13957*4)*(0.13957*4);
    double dh = 0.001;
    int od  = (sc - 0.312)/dh;
    double sc_m = 0.312 + od*dh;
    double rhouse = 0;
    if(sc>=0.312 && sc<1){
        rhouse = ((sc-sc_m)/dh)*(rho[od+1]-rho[od])+rho[od];
    }else if(sc<0.312 && sc>=smin){
        rhouse = ((sc - smin)/(0.312-smin))*rho[0];
    }else if(sc>=1){
        //      rhouse = (1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc);
        rhouse = sqrt(1-16*m2/sc);
    }else{
        rhouse = 0;
    }
    return rhouse;
}
 
complex<double> EvtD0ToKKpi0::rhoMTX(int i, int j, double s){
 
    double rhoijx;  
    double rhoijy;  
    double mpi = 0.13957;
    if(i==j && i==0 ){
        double m2 = 0.13957*0.13957;
        if((1-(4*m2)/s)>0){
            rhoijx = sqrt(1.0f - (4*m2)/s);
            rhoijy = 0;
        }else{
            rhoijy = sqrt((4*m2)/s - 1.0f);
            rhoijx = 0;
        }
    }
    if(i==j && i==1 ){
        double m2 = 0.493677*0.493677;
        if((1-(4*m2)/s)>0){
            rhoijx = sqrt(1.0f - (4*m2)/s);
            rhoijy = 0;
        }else{
            rhoijy = sqrt((4*m2)/s - 1.0f);
            rhoijx = 0;
        }
    }
    if(i==j && i==2){
        rhoijx = rho22(s);
        rhoijy = 0;
    }
    if(i==j && i==3){
        double m2 = 0.547862*0.547862;
        if((1-(4*m2)/s)>0){
            rhoijx = sqrt(1.0f - (4*m2)/s);
            rhoijy = 0;
        }else{
            rhoijy = sqrt((4*m2)/s - 1.0f);
            rhoijx = 0;
        }
    }
    if(i==j && i==4){
        double m_1 = 0.547862;
        double m_2 = 0.95778;
        double mp2 = (m_1+m_2)*(m_1+m_2);
        double mm2 = (m_1-m_2)*(m_1-m_2);
        if((1-mp2/s)>0){
            rhoijx = sqrt(1.0f - mp2/s);
            rhoijy = 0;
        }else{
            rhoijy = sqrt(mp2/s - 1.0f);
            rhoijx = 0;
        }
    }
 
    if(i!=j){
        rhoijx = 0;
        rhoijy = 0;
    }
    complex<double> rhoij(rhoijx,rhoijy);
    return rhoij;
 
}
/* Kij = (sum_a gigj/(ma^2-s) + fij(1-s0)/(s-s0))((s-0.5sampi2)(1-sa0)/(s-sa0))  */
complex<double> EvtD0ToKKpi0::KMTX(int i, int j, double s){
 
    double Kijx;
    double Kijy;
    double mpi = 0.13957;
    double m[5]  = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
 
    double g1[5] = { 0.22889,-0.55377, 0.00000,-0.39899,-0.34639};
    double g2[5] = { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503};
    double g3[5] = { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681};
    double g4[5] = { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984};
    double g5[5] = { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358};
 
    double f1[5] = { 0.23399, 0.15044,-0.20545, 0.32825, 0.35412};
 
    double eps   = 1e-11; 
 
    double down[5]   = { 0,0,0,0,0};
    double upreal[5] = { 0,0,0,0,0};
    double upimag[5] = { 0,0,0,0,0};
 
    for(int k=0; k<5; k++){
 
        /*    down[k] = (m[k]*m[k]-s)*(m[k]*m[k]-s)+eps*eps;
              upreal[k] = (m[k]*m[k]-s)/down[k];
              upimag[k] = -1.0f * eps/down[k];  */
 
        double dm2 = m[k]*m[k]-s;
        if(fabs(dm2)<eps && dm2<=0) dm2 = -eps;
        if(fabs(dm2)<eps && dm2>0)  dm2 = eps;
        upreal[k] = 1.0f/dm2;
        upimag[k] = 0;  
    }
 
    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];
    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];
 
    double tmp2  = 0;
    if(i==0){
        tmp2 = f1[j]*(1+3.92637)/(s+3.92637);
    }
    if(j==0){
        tmp2 = f1[i]*(1+3.92637)/(s+3.92637);
    } 
    double tmp3 = (s-0.5*mpi*mpi)*(1+0.15)/(s+0.15);
 
    Kijx = (tmp1x+tmp2)*tmp3;
    Kijy = (tmp1y)*tmp3;
 
    complex<double> Kij(Kijx,Kijy);
    return Kij;
}
 
complex<double> EvtD0ToKKpi0::IMTX(int i, int j){
 
   double Iijx;
   double Iijy;
   if(i==j){
      Iijx = 1;
      Iijy = 0;
   }else{
      Iijx = 0;
      Iijy = 0;
   }
   complex<double> Iij(Iijx,Iijy);
   return Iij;
 
}
 
/* F = I - i*K*rho */
complex<double> EvtD0ToKKpi0::FMTX(double Kijx, double Kijy, double rhojjx, double rhojjy, int i, int j){
 
   double Fijx;
   double Fijy;
 
   double tmpx = rhojjx*Kijx - rhojjy*Kijy;
   double tmpy = rhojjx*Kijy + rhojjy*Kijx;
 
   Fijx = IMTX(i,j).real() + tmpy;
   Fijy = -tmpx;
 
   complex<double> Fij(Fijx,Fijy);
   return Fij;
}
 
/* inverse for Matrix (I - i*rho*K) using LUP */
double EvtD0ToKKpi0::FINVMTX(double s, double *FINVx, double *FINVy){
 
    int P[5] = { 0,1,2,3,4};
 
    double Fx[5][5];
    double Fy[5][5];
 
    double Ux[5][5];
    double Uy[5][5];
    double Lx[5][5];
    double Ly[5][5];
 
    double UIx[5][5];
    double UIy[5][5];
    double LIx[5][5];
    double LIy[5][5];
 
    for(int k=0; k<5; k++){
        double rhokkx = rhoMTX(k,k,s).real();
        double rhokky = rhoMTX(k,k,s).imag();
        Ux[k][k] = rhokkx;
        Uy[k][k] = rhokky;
        for(int l=k; l<5; l++){
            double Kklx = KMTX(k,l,s).real();
            double Kkly = KMTX(k,l,s).imag();
            Lx[k][l] = Kklx;
            Ly[k][l] = Kkly;
            Lx[l][k] = Lx[k][l];
            Ly[l][k] = Ly[k][l];
        }
    }
 
    for(int k=0; k<5; k++){
        for(int l=0; l<5; l++){
            double Fklx = FMTX(Lx[k][l],Ly[k][l],Ux[l][l],Uy[l][l],k,l).real();
            double Fkly = FMTX(Lx[k][l],Ly[k][l],Ux[l][l],Uy[l][l],k,l).imag();
            Fx[k][l] = Fklx;
            Fy[k][l] = Fkly;
        }
    }
 
    for(int k=0; k<5; k++){
        double tmprM = (Fx[k][k]*Fx[k][k]+Fy[k][k]*Fy[k][k]);
        int tmpID = 0;
        for(int l=k; l<5; l++){
            double tmprF = (Fx[l][k]*Fx[l][k]+Fy[l][k]*Fy[l][k]);
            if(tmprM<=tmprF){
                tmprM = tmprF;
                tmpID = l;
            }
        }
 
        int tmpP = P[k];
        P[k]     = P[tmpID];
        P[tmpID] = tmpP;
 
        for(int l=0; l<5; l++){
 
            double tmpFx = Fx[k][l];
            double tmpFy = Fy[k][l];
 
            Fx[k][l] = Fx[tmpID][l];
            Fy[k][l] = Fy[tmpID][l];
 
            Fx[tmpID][l] = tmpFx;
            Fy[tmpID][l] = tmpFy;
 
        }
 
        for(int l=k+1; l<5; l++){
            double rFkk = Fx[k][k]*Fx[k][k] + Fy[k][k]*Fy[k][k];
            double Fxlk = Fx[l][k];
            double Fylk = Fy[l][k];
            double Fxkk = Fx[k][k];
            double Fykk = Fy[k][k];
            Fx[l][k] = (Fxlk*Fxkk + Fylk*Fykk)/rFkk;
            Fy[l][k] = (Fylk*Fxkk - Fxlk*Fykk)/rFkk;
            for(int m=k+1; m<5; m++){
                Fx[l][m] = Fx[l][m] - (Fx[l][k]*Fx[k][m] - Fy[l][k]*Fy[k][m]);
                Fy[l][m] = Fy[l][m] - (Fx[l][k]*Fy[k][m] + Fy[l][k]*Fx[k][m]);
            }
        }
    }
 
    for(int k=0; k<5; k++){
        for(int l=0; l<5 ;l++){
            if(k==l){
                Lx[k][k] = 1;
                Ly[k][k] = 0;
                Ux[k][k] = Fx[k][k];
                Uy[k][k] = Fy[k][k];
            }
            if(k>l){
                Lx[k][l] = Fx[k][l];
                Ly[k][l] = Fy[k][l];
                Ux[k][l] = 0;
                Uy[k][l] = 0;
            }
            if(k<l){
                Ux[k][l] = Fx[k][l];
                Uy[k][l] = Fy[k][l];
                Lx[k][l] = 0;
                Ly[k][l] = 0;
            }
        }
    }
 
    // calculate Inverse for L and U
    for(int k=0; k<5; k++){
 
        LIx[k][k] = 1;
        LIy[k][k] = 0;
 
        double rUkk = Ux[k][k]*Ux[k][k] + Uy[k][k]*Uy[k][k];
        UIx[k][k] = Ux[k][k]/rUkk;
        UIy[k][k] = -1.0f * Uy[k][k]/rUkk ;
 
        for(int l=(k+1); l<5; l++){
            LIx[k][l] = 0;
            LIy[k][l] = 0;
            UIx[l][k] = 0;
            UIy[l][k] = 0;
        }
 
        for(int l=(k-1); l>=0; l--){ // U-1
            double sx   = 0;
            double c_sx = 0;
            double sy   = 0;
            double c_sy = 0;
            for(int m=l+1; m<=k; m++){
                sx = sx - c_sx;
                double sx_tmp = sx + Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k]; 
                c_sx = (sx_tmp - sx) - (Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k]);
                sx = sx_tmp;
 
                sy = sy - c_sy;
                double sy_tmp = sy + Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k]; 
                c_sy = (sy_tmp - sy) - (Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k]);
                sy = sy_tmp;
            }
            UIx[l][k] = -1.0f * (UIx[l][l]*sx - UIy[l][l]*sy);
            UIy[l][k] = -1.0f * (UIy[l][l]*sx + UIx[l][l]*sy);
        }
 
        for(int l=k+1; l<5; l++){ // L-1
            double sx = 0;
            double c_sx = 0;
            double sy = 0;
            double c_sy = 0;
            for(int m=k; m<l; m++){
                sx = sx - c_sx;
                double sx_tmp = sx + Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k];
                c_sx = (sx_tmp - sx) - (Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k]);
                sx = sx_tmp;
 
                sy = sy - c_sy;
                double sy_tmp = sy + Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k];
                c_sy = (sy_tmp - sy) - (Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k]);
                sy = sy_tmp;
            }
            LIx[l][k] = -1.0f * sx;
            LIy[l][k] = -1.0f * sy;
        }
    }
 
    for(int m=0; m<5; m++){
        double resX = 0; 
        double c_resX = 0; 
        double resY = 0; 
        double c_resY = 0; 
        for(int k=0; k<5; k++){
            for(int l=0; l<5; l++){
                double Plm = 0;
                if(P[l] == m) Plm = 1;
 
                resX = resX - c_resX;
                double resX_tmp = resX + (UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm;
                c_resX = (resX_tmp - resX) - ((UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm);
                resX = resX_tmp;
 
                resY = resY - c_resY;
                double resY_tmp = resY + (UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm;
                c_resY = (resY_tmp - resY) - ((UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm);
                resY = resY_tmp;
            }
        }
        FINVx[m] = resX;
        FINVy[m] = resY;
    }
 
    return 1.0f;
 
}
 
complex<double> EvtD0ToKKpi0::PVTR(int ID, double s){
 
   double VPix;   
   double VPiy;   
   double m[5]  = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
 
   double eps   = 1e-11; 
 
   /*  double down   = (m[ID]*m[ID]-s)*(m[ID]*m[ID]-s)+eps*eps;
       double upreal = (m[ID]*m[ID]-s)/down;
       double upimag = -1.0f * eps/down;  */
 
   double dm2 = m[ID]*m[ID]-s;
 
   if(fabs(dm2)<eps && dm2<=0) dm2 = -eps;
   if(fabs(dm2)<eps && dm2>0)  dm2 = eps;
 
   VPix = 1.0f/dm2;
   VPiy = 0; 
 
   complex<double> VPi(VPix,VPiy);
   return VPi;
}
 
complex<double> EvtD0ToKKpi0::Fvector(double sa, double s0, int l){
 
   double outputx = 0;
   double outputy = 0;
 
   double FINVx[5] = {0,0,0,0,0};
   double FINVy[5] = {0,0,0,0,0};
 
   double tmpFLAG = FINVMTX(sa,FINVx,FINVy);
 
   if(l<5){
      double g[5][5] = {{ 0.22889,-0.55377, 0.00000,-0.39899,-0.34639},
         { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503},
         { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681},
         { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984},
         { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358}};
      double resx = 0;
      double c_resx = 0;
      double resy = 0;
      double c_resy = 0;
      double Plx = PVTR(l,sa).real();
      double Ply = PVTR(l,sa).imag();
      for(int j=0; j<5; j++){
         resx = resx - c_resx;
         double resx_tmp = resx + (FINVx[j]*g[l][j]*Plx - FINVy[j]*g[l][j]*Ply);
         c_resx = (resx_tmp - resx) - (FINVx[j]*g[l][j]*Plx - FINVy[j]*g[l][j]*Ply);
         resx = resx_tmp;
 
         resy = resy - c_resy;
         double resy_tmp = resy + (FINVx[j]*g[l][j]*Ply + FINVy[j]*g[l][j]*Plx);
         c_resy = (resy_tmp - resy) - (FINVx[j]*g[l][j]*Ply + FINVy[j]*g[l][j]*Plx);
         resy = resy_tmp;
      }
      outputx = resx;
      outputy = resy;
   }else{
      int idx = l-5;
      double eps = 1e-11; 
      double ds = sa-s0;
      if(fabs(ds)<eps && ds<=0) ds = -eps;
      if(fabs(ds)<eps && ds>0)  ds = eps;
      double tmp = (1-s0)/ds;
      outputx = FINVx[idx]*tmp;
      outputy = FINVy[idx]*tmp;
   }
 
   complex<double> output(outputx,outputy);
   return output;
}
 
complex<double> EvtD0ToKKpi0::CalD0Amp(){
   return Amp(p4_Kap, p4_Kam, p4_Pi0);
}
complex<double> EvtD0ToKKpi0::CalDbAmp(){
 
    vector<double> cpKap;  cpKap.clear();
    vector<double> cpKam;  cpKam.clear();
    vector<double> cpPi0;  cpPi0.clear();
 
    cpKap.push_back(-p4_Kam[0]); cpKam.push_back(-p4_Kap[0]); cpPi0.push_back(-p4_Pi0[0]); 
    cpKap.push_back(-p4_Kam[1]); cpKam.push_back(-p4_Kap[1]); cpPi0.push_back(-p4_Pi0[1]); 
    cpKap.push_back(-p4_Kam[2]); cpKam.push_back(-p4_Kap[2]); cpPi0.push_back(-p4_Pi0[2]); 
    cpKap.push_back( p4_Kam[3]); cpKam.push_back( p4_Kap[3]); cpPi0.push_back( p4_Pi0[3]); 
 
    return (-1.0)*Amp(cpKap, cpKam, cpPi0); 
 
}
 
complex<double> EvtD0ToKKpi0::Amp(vector<double> Kap, vector<double> Kam, vector<double> Pi0)
{
 
    if(fitpara.size()!=5) {
        cout<<"Error!!! number of para: "<<fitpara.size()<<endl;
        exit(0);
    }
 
    vector<double> KapKam; KapKam.clear();
    vector<double> KapPi0; KapPi0.clear();
    vector<double> KamPi0; KamPi0.clear();
 
    KapKam  = sum_tensor(Kap, Kam);
    KapPi0  = sum_tensor(Kap, Pi0);
    KamPi0  = sum_tensor(Kam, Pi0);
 
    vector<double> D0; D0.clear();
    D0 = sum_tensor(KapKam, Pi0);
 
    double M2_KapKam  = contract_11_0(KapKam, KapKam);
    double M2_KapPi0  = contract_11_0(KapPi0, KapPi0);
    double M2_KamPi0  = contract_11_0(KamPi0, KamPi0);
 
    double M2_D0 = contract_11_0(D0, D0);
 
    complex<double> RBW_Kst892_p = RBW(M2_KapPi0, m0_Kst892, w0_Kst892, m2_Ka, m2_Pi0, rRes, 1);
    complex<double> RBW_Kst892_m = RBW(M2_KamPi0, m0_Kst892, w0_Kst892, m2_Ka, m2_Pi0, rRes, 1);
    complex<double> RBW_Kst1410_m = RBW(M2_KamPi0, m0_Kst1410, w0_Kst1410, m2_Ka, m2_Pi0, rRes, 1);
 
    complex<double> RBW_phi1020  = RBW(M2_KapKam, m0_phi1020, w0_phi1020,-1, -1,-1,-1);
 
    complex<double> KpPi0S_Lass  = KPiSLass(M2_KapPi0, m0_K01430, w0_K01430, a_lass, r_lass,  0.187955,  1.82603, -0.10329550);
    complex<double> FT_f0980     = Flatte(M2_KapKam, m0_f0980, g1_f0980, g2_f0980, m_Pi, m_Pi, m_Ka, m_Ka);
 
    // X->PP Orbital
    vector<double> T1_KapKam; T1_KapKam.clear();
    vector<double> T1_KapPi0; T1_KapPi0.clear();
    vector<double> T1_KamPi0; T1_KamPi0.clear();
 
    T1_KapKam = OrbitalTensors(KapKam, Kap, Kam, rRes, 1);
    T1_KapPi0 = OrbitalTensors(KapPi0, Kap, Pi0, rRes, 1);
    T1_KamPi0 = OrbitalTensors(KamPi0, Kam, Pi0, rRes, 1);
 
    vector<double> T2_KapKam; T2_KapKam.clear();
 
    T2_KapKam = OrbitalTensors(KapKam, Kap, Kam, rRes, 2);
 
    // D->XP Orbital
    vector<double> T1_KapKamPi0; T1_KapKamPi0.clear();
    vector<double> T1_KapPi0Kam; T1_KapPi0Kam.clear();
    vector<double> T1_KamPi0Kap; T1_KamPi0Kap.clear();
 
    T1_KapKamPi0 = OrbitalTensors(D0, KapKam, Pi0, rD, 1);
    T1_KapPi0Kam = OrbitalTensors(D0, KapPi0, Kam, rD, 1);
    T1_KamPi0Kap = OrbitalTensors(D0, KamPi0, Kap, rD, 1);
 
    vector<double> T2_KapKamPi0; T2_KapKamPi0.clear();
 
    T2_KapKamPi0 = OrbitalTensors(D0, KapKam, Pi0, rD, 2);
 
    complex<double> amplitude(0,0);
    
 
    // D-> VP
    double SF_VpPm = contract_11_0(T1_KapPi0Kam, T1_KapPi0);
    amplitude += fitpara[0]*(SF_VpPm*RBW_Kst892_p);
 
    double SF_VmPp = contract_11_0(T1_KamPi0Kap, T1_KamPi0);
    amplitude += fitpara[1]*(SF_VmPp*RBW_Kst892_m);
 
    double SF_VzPz = contract_11_0(T1_KapKamPi0, T1_KapKam);
    amplitude += fitpara[2]*(SF_VzPz*RBW_phi1020);
 
    amplitude += fitpara[3]*(SF_VmPp*RBW_Kst1410_m);
 
    // D-> SP
    amplitude += fitpara[4]*KpPi0S_Lass;
 
 
    return amplitude;
 
}
 
int EvtD0ToKKpi0::CalAmp()
{
   m_AmpD0 = CalD0Amp();
   m_AmpDb = CalDbAmp();
   return 0;
}
 
double EvtD0ToKKpi0::mag2(complex<double> x)
{
   double temp = x.real()*x.real() + x.imag()*x.imag();
   return temp;
}
 
double EvtD0ToKKpi0::mag2_itf(complex<double> x, complex<double> y)
{
    double temp = 2*(x.real()*y.real() + x.imag()*y.imag());
    return temp;
}
 
double EvtD0ToKKpi0::arg(complex<double> x)
{
   double temp = atan(x.imag()/x.real());
   if(x.real()<0) temp=temp+PI;
   return temp;
}
double EvtD0ToKKpi0::Get_strongPhase()
{
   double temp = arg(m_AmpD0) - arg(m_AmpDb);
   while (temp < -PI){
      temp += 2.0*PI;
   }
   while (temp > PI){
      temp -= 2.0*PI;
   }  
   return temp;
}
 
double EvtD0ToKKpi0::AmplitudeSquare(int charm, int tagmode){
 
    EvtVector4R dp1=GetDaugMomLab(0),dp2=GetDaugMomLab(1),dp3=GetDaugMomLab(2); // D0 -> K+ K- pi0
    EvtVector4R mp = dp1 + dp2 + dp3;
 
  double emp = mp.get(0);
  EvtVector3R boostp3(-mp.get(1)/emp, -mp.get(2)/emp, -mp.get(3)/emp);
 
    EvtVector4R dp1bst = boostTo(dp1, boostp3);
    EvtVector4R dp2bst = boostTo(dp2, boostp3);
    EvtVector4R dp3bst = boostTo(dp3, boostp3);
 
    double p4kp[4], p4km[4], p4pi0[4];
    for(int i=0; i<3; i++){
        p4kp[i] =dp1bst.get(i+1);
        p4km[i] =dp2bst.get(i+1);
        p4pi0[i]=dp3bst.get(i+1);
    }
 
    p4kp[3] =dp1bst.get(0);
    p4km[3] =dp2bst.get(0);
    p4pi0[3]=dp3bst.get(0);
 
    setInput(p4kp, p4km, p4pi0);
    CalAmp();
 
    complex<double> ampD0, ampDb;
  if(charm>0){
     ampD0 = Get_AmpD0();
     ampDb = Get_AmpDb();
  }else{
     ampD0 = Get_AmpDb();
     ampDb = Get_AmpD0();
  }
 
  double ampsq = 1e-20;
  double r_tag = 0, R_tag = 0, delta_tag = 0;
   
  if (tagmode==1||tagmode==2||tagmode==3) {
     if(tagmode == 1){// Kpi
        r_tag = 0.0586;
        R_tag = 1;
        delta_tag = 192.1/180.0*3.1415926;
     }else if(tagmode == 2){// Kpipi0
        r_tag = 0.0441;
        R_tag = 0.79;
        delta_tag = 196.0/180.0*3.1415926;
     }else if(tagmode == 3){// K3pi
        r_tag = 0.0550;
        R_tag = 0.44;
        delta_tag = 161.0/180.0*3.1415926;
     }
     complex<double> qcf(r_tag*R_tag*cos(-delta_tag), r_tag*R_tag*sin(-delta_tag));
     complex<double> ampD0_part1 = ampD0 - qcf*ampDb;
     ampsq = mag2(ampD0_part1) + r_tag*r_tag*(1-R_tag*R_tag)*(mag2(ampDb));
  } else {
     ampsq = mag2(ampD0);
  }
 
  return (ampsq <= 0) ? 1e-20 : ampsq;
}