EvtD0To2pip2pim.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:  EvtD0To2pip2pim.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/EvtVector4R.hh"
#include "EvtGenBase/EvtVector3R.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtD0To2pip2pim.hh"
#include <stdlib.h>
#include <iostream>
#include <string>
#include <complex>
#include <vector>
#include <math.h>
#include "TMath.h"
using namespace std;
 
EvtD0To2pip2pim::~EvtD0To2pip2pim() {}
 
void EvtD0To2pip2pim::getName(std::string& model_name){
  model_name="D0To2pip2pim";     
}
 
EvtDecayBase* EvtD0To2pip2pim::clone(){
  return new EvtD0To2pip2pim;
}
 
void EvtD0To2pip2pim::init(){
  checkNArg(2);
  checkNDaug(4);
  charm   = getArg(0);
  tagmode = getArg(1);
  //std::cout << "Initializing EvtD0To2pip2pim: charm= "<<charm<<" tagmode= "<<tagmode<<std::endl;
 
  double mag[30], pha[30];
  mag[0]= 100.0;    pha[0]= 0.0;
  mag[1]= 7.95507;  pha[1]= -0.0687407;
  mag[2]= 37.5559;  pha[2]= -1.74946;
  mag[3]= 61.2172;  pha[3]= 2.98079;
  mag[4]= 187.79;   pha[4]= 2.64471;
  mag[5]= 385.474;  pha[5]= -0.137107;
  mag[6]= 0.330788; pha[6]= 0.268133;
  mag[7]= 127.158;  pha[7]= -2.47773;
  mag[8]= 339.914;  pha[8]= 2.22856;
  mag[9]= 0.320888; pha[9]= -2.6194;
  mag[10]=0.366283; pha[10]=-0.26867;
  mag[11]=86.0865;  pha[11]=-2.49649;
  mag[12]=6.1541;   pha[12]=-1.18299;
  mag[13]=56.6067;  pha[13]=0.142977;
  mag[14]=92.3073;  pha[14]=-2.15881;
  mag[15]=10.5885;  pha[15]=-3.03166;
  mag[16]=8.36765;  pha[16]=1.8417;
  mag[17]=6.56437;  pha[17]=-2.93087;
  mag[18]=15.7197;  pha[18]=0.96925;
  mag[19]=21.4195;  pha[19]=-1.23701;
  mag[20]=56.8867;  pha[20]=-0.385837;
  mag[21]=231.626;  pha[21]=2.14842;
  mag[22]=2938.45;  pha[22]=-0.693491;
  mag[23]=7252.7;   pha[23]=2.23659;
  mag[24]=5165.87;  pha[24]=0.913557;
  mag[25]=11508.6;  pha[25]=-1.07187;
  mag[26]=2461.86;  pha[26]=1.8709;
  mag[27]=8757.75;  pha[27]=2.40756;
  mag[28]=19.7413;  pha[28]=-1.0753;
  mag[29]=66.3826;  pha[29]=2.34666;
 
  fitpara.clear();
  for(int i=0; i<30; i++){
     complex<double> ctemp(mag[i]*cos(pha[i]),mag[i]*sin(pha[i]));
     fitpara.push_back(ctemp);
  }
 
  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 = 4;
  math_pi = 3.1415926f;
  mass_Pion = 0.13957f;
 
  rRes = 3.0*0.197321;
  rD   = 5.0*0.197321;
  m_Pi  = 0.139570;
  m2_Pi = m_Pi*m_Pi;
 
  m0_rho770 = 0.77526;
  w0_rho770 = 0.1478;
 
  m0_rho1450 = 1.465;
  w0_rho1450 = 0.400;
 
  m0_f21270 = 1.2755;
  w0_f21270 = 0.1867;
 
  m0_a11260 = 1.1337;
  g1_a11260 = 0.00335;
  g2_a11260 = 0.0;
 
  m0_pi1300 = 1.498;
  w0_pi1300 = 0.590;
 
  m0_a11420 = 1.411;
  w0_a11420 = 0.161;
 
  m0_a11640 = 1.655;
  w0_a11640 = 0.254;
 
  m0_a21320 = 1.3186;
  w0_a21320 = 0.105;
 
  m0_pi11400 = 1.354;
  w0_pi11400 = 0.330;
 
  s0_prod = -5.0;
 
  return;
}
 
void EvtD0To2pip2pim::initProbMax(){
  setProbMax(1999257515.2);
}
 
void EvtD0To2pip2pim::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 EvtD0To2pip2pim::setInput(double* pip1, double* pim1, double* pip2, double* pim2){
   m_Pip1.clear(); m_Pim1.clear(); m_Pip2.clear(); m_Pim2.clear();
   m_Pip1.push_back(pip1[0]); m_Pim1.push_back(pim1[0]); m_Pip2.push_back(pip2[0]); m_Pim2.push_back(pim2[0]);
   m_Pip1.push_back(pip1[1]); m_Pim1.push_back(pim1[1]); m_Pip2.push_back(pip2[1]); m_Pim2.push_back(pim2[1]);
   m_Pip1.push_back(pip1[2]); m_Pim1.push_back(pim1[2]); m_Pip2.push_back(pip2[2]); m_Pim2.push_back(pim2[2]);
   m_Pip1.push_back(pip1[3]); m_Pim1.push_back(pim1[3]); m_Pip2.push_back(pip2[3]); m_Pim2.push_back(pim2[3]);
}
 
vector<double> EvtD0To2pip2pim::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 EvtD0To2pip2pim::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> EvtD0To2pip2pim::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 EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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 EvtD0To2pip2pim::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> EvtD0To2pip2pim::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 EvtD0To2pip2pim::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 EvtD0To2pip2pim::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 EvtD0To2pip2pim::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 EvtD0To2pip2pim::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 EvtD0To2pip2pim::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 EvtD0To2pip2pim::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);
   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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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;
}
 
/* build a1260 propagator */
double EvtD0To2pip2pim::widT1260(int i, double g1, double g2){
 
  double wid1[300]  = { 0.00100302, 0.0069383, 0.0223132, 0.0504984, 0.093998, 0.154569, 0.233464, 0.331844, 0.450141, 0.589068, 
  0.748192, 0.928578, 1.13001, 1.35227, 1.59548, 1.86005, 2.14633, 2.45252, 2.78199, 3.13055, 
  3.50351, 3.89773, 4.31274, 4.75409, 5.21133, 5.69991, 6.20735, 6.74638, 7.30128, 7.8858, 
  8.50289, 9.14654, 9.82395, 10.5209, 11.2643, 12.0436, 12.8585, 13.692, 14.598, 15.5291, 
  16.5158, 17.5337, 18.6289, 19.7599, 20.9847, 22.2557, 23.5959, 25.0095, 26.5123, 28.0789, 
  29.7542, 31.5143, 33.3769, 35.3462, 37.3911, 39.5988, 41.874, 44.2815, 46.7975, 49.401, 
  52.0553, 54.7753, 57.5932, 60.4542, 63.3049, 66.0665, 68.8987, 71.6282, 74.2613, 76.8713, 
  79.3528, 81.722, 84.1212, 86.227, 88.4243, 90.3478, 92.2478, 94.1483, 95.8541, 97.5086, 
  99.0092, 100.48, 101.861, 103.153, 104.338, 105.576, 106.696, 107.647, 108.761, 109.725, 
  110.625, 111.529, 112.426, 113.01, 113.877, 114.647, 115.086, 115.856, 116.533, 117.076, 
  117.646, 118.25, 118.653, 119.023, 119.554, 119.958, 120.384, 121.036, 121.402, 121.686, 
  122.44, 122.592, 122.979, 123.39, 123.819, 123.957, 124.459, 124.681, 125.071, 125.405, 
  125.769, 125.978, 126.542, 126.817, 127.017, 127.292, 127.765, 127.989, 128.542, 128.66, 
  128.923, 129.094, 129.441, 129.716, 130.23, 130.506, 130.658, 131.12, 131.308, 131.579, 
  131.994, 132.28, 132.594, 132.79, 133.107, 133.589, 133.935, 134.242, 134.484, 134.765, 
  135.208, 135.58, 135.922, 136.236, 136.545, 136.949, 137.216, 137.503, 137.994, 138.35, 
  138.62, 138.912, 139.413, 139.831, 140.137, 140.478, 141, 141.3, 141.807, 142.291, 
  142.864, 143.315, 143.678, 144.215, 144.587, 145.122, 145.8, 145.885, 146.583, 147.226, 
  147.661, 148.187, 148.698, 149.227, 149.832, 150.548, 151.122, 151.674, 152.074, 152.666, 
  153.295, 153.899, 154.661, 155.364, 155.908, 156.495, 157.36, 157.719, 158.533, 159.287, 
  159.79, 160.654, 161.257, 161.93, 162.437, 163.468, 163.957, 164.631, 165.414, 166.203, 
  166.738, 167.61, 168.453, 169.101, 170.111, 170.333, 171.123, 171.958, 173.018, 173.663, 
  174.213, 175.241, 175.579, 176.435, 177.291, 178.071, 178.969, 179.635, 180.118, 181.078, 
  182.007, 182.73, 183.282, 184.161, 184.981, 185.695, 186.506, 187.16, 187.996, 188.439, 
  189.416, 190.104, 190.759, 191.786, 192.331, 193.318, 193.836, 194.981, 195.634, 196.231, 
  196.832, 197.835, 198.608, 199.273, 199.854, 200.695, 201.719, 202.105, 202.958, 203.707, 
  204.306, 205.319, 205.977, 206.875, 207.687, 208.352, 209.04, 209.352, 210.313, 211.322, 
  212.02, 212.458, 213.246, 214.331, 214.923, 215.466, 216.536, 217.346, 217.867, 218.463, 
  219.201, 219.88, 220.829, 221.461, 222.399, 223.068, 223.712, 224.174, 224.837, 225.838, 
  227.019, 227.171, 227.797, 228.663, 229.429, 230.323, 230.845, 231.574, 232.417, 232.677 };
  double wid2[300]  = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
      1.87136e-06, 1.50063e-05, 5.10425e-05, 0.000122121, 0.000240853, 0.000420318, 0.000675161, 0.0010173, 0.00146434, 0.00203321, 
      0.00273489, 0.0035927, 0.00462579, 0.00584255, 0.00727372, 0.00895462, 0.0108831, 0.013085, 0.0156197, 0.0184865, 
      0.0217078, 0.0253423, 0.0294103, 0.0339191, 0.0389837, 0.0446351, 0.0508312, 0.0577268, 0.0653189, 0.0737049, 
      0.0829819, 0.0930611, 0.104328, 0.116663, 0.130105, 0.144922, 0.16122, 0.179091, 0.198759, 0.220133, 
      0.243916, 0.269803, 0.298861, 0.330061, 0.365741, 0.40437, 0.447191, 0.49501, 0.548576, 0.606445, 
      0.674414, 0.748353, 0.831686, 0.929938, 1.03771, 1.16187, 1.30387, 1.47341, 1.65629, 1.88318, 
      2.14353, 2.44169, 2.79831, 3.2009, 3.65522, 4.16317, 4.69597, 5.2585, 5.85965, 6.44984, 
      7.04202, 7.60113, 8.14571, 8.73195, 9.24537, 9.75717, 10.2093, 10.6731, 11.1487, 11.5819, 
      12.0158, 12.4253, 12.8113, 13.2073, 13.5995, 13.9317, 14.312, 14.6595, 14.9511, 15.2668, 
      15.6092, 15.9349, 16.1873, 16.5049, 16.819, 17.0743, 17.3621, 17.6094, 17.8418, 18.0681, 
      18.3141, 18.5914, 18.8187, 19.0562, 19.2282, 19.4918, 19.7326, 19.9112, 20.134, 20.3386, 
      20.511, 20.6865, 20.8958, 21.0518, 21.2967, 21.44, 21.6361, 21.8012, 21.9523, 22.1736, 
      22.2615, 22.4207, 22.6056, 22.7198, 22.9299, 23.0605, 23.2959, 23.3808, 23.4961, 23.6793, 
      23.7843, 23.9697, 24.0689, 24.1919, 24.405, 24.3898, 24.6018, 24.7294, 24.789, 24.9978, 
      25.0626, 25.1728, 25.2809, 25.3579, 25.5444, 25.5995, 25.7644, 25.8397, 25.9229, 26.095, 
      26.1495, 26.2899, 26.3871, 26.54, 26.6603, 26.7008, 26.7836, 26.907, 26.9653, 26.9969, 
      27.1226, 27.226, 27.3543, 27.4686, 27.4887, 27.6163, 27.6986, 27.7506, 27.7884, 27.8662, 
      27.9886, 28.0573, 28.1238, 28.2612, 28.3209, 28.3457, 28.4392, 28.5086, 28.6399, 28.7603, 
      28.788, 28.8502, 28.9038, 28.9667, 28.975, 29.0032, 29.2681, 29.2392, 29.2572, 29.3364 };
 
   return wid1[i]*g1+wid2[i]*g2;
}
 
double EvtD0To2pip2pim::anywid1260(double sc, double g1, double g2){
 
   double smin = (0.13957*3)*(0.13957*3);
   double dh = 0.01;
   int od  = (sc - 0.18)/dh;
   double sc_m = 0.18 + od*dh;
   double widuse = 0;
   if(sc>=0.18 && sc<=3.17){
      widuse = ((sc-sc_m)/dh)*(widT1260(od+1,g1,g2)-widT1260(od,g1,g2))+widT1260(od,g1,g2);
   }else if(sc<0.18 && sc>smin){
      widuse = ((sc - smin)/(0.18-smin))*widT1260(0,g1,g2);
   }else if(sc>3.17){
      widuse = widT1260(299,g1,g2);
   }else{
      widuse = 0;
   }
   return widuse;
 
}
 
complex<double> EvtD0To2pip2pim::RBWa1260(double mx2, double mr, double g1, double g2){
 
   double mx  = sqrt(mx2);
   double mr2 = mr*mr;
   double wid0 = anywid1260(mx2,g1,g2);
 
   double denom_real = mr2-mx2;
   double denom_imag = mr*wid0;//real-i*imag;
 
   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;
 
}
 
/* build pi1300 propagator */
double EvtD0To2pip2pim::widT1300(int i){
 
   double wid1[300]  = { 0.0702928, 0.399073, 0.991742, 1.82025, 2.85953, 4.08606, 5.48082, 7.02683, 8.70496, 10.5007, 
   12.4053, 14.4026, 16.4831, 18.6423, 20.8642, 23.1544, 25.4896, 27.8703, 30.3015, 32.7861, 
   35.2622, 37.8173, 40.3819, 42.974, 45.5732, 48.2303, 50.8659, 53.5741, 56.28, 59.0242, 
   61.738, 64.5642, 67.377, 70.1605, 73.0155, 75.8849, 78.7611, 81.7366, 84.7156, 87.7527, 
   90.7217, 93.8402, 96.8516, 100.036, 103.168, 106.483, 109.772, 113.098, 116.491, 120.013, 
   123.618, 127.069, 130.983, 134.868, 138.605, 142.625, 147.007, 151.154, 155.625, 160.1, 
   164.776, 169.651, 174.646, 179.669, 185.084, 190.409, 196.147, 201.788, 207.901, 214.041, 
   220.327, 226.505, 233.334, 239.816, 246.878, 253.563, 260.393, 267.453, 274.5, 282.15, 
   289.014, 296.45, 303.808, 311.427, 318.649, 326.965, 334.298, 341.576, 349.715, 356.89, 
   365.029, 372.677, 379.882, 387.677, 395.178, 402.445, 410.353, 418.649, 424.994, 432.156, 
   440.002, 448.394, 454.382, 460.97, 468.446, 475.847, 481.956, 489.729, 496.094, 501.22, 
   509.278, 514.618, 521.06, 528.247, 534.246, 540.312, 547.316, 552.549, 559.193, 566.059, 
   572.882, 578.147, 585.118, 589.989, 596.717, 601.222, 607.749, 613.96, 621.107, 625.218, 
   630.396, 635.57, 641.175, 646.024, 651.984, 657.156, 661.385, 666.804, 672.088, 675.939, 
   681.207, 685.072, 690.63, 694.767, 699.469, 704.1, 709.445, 713.704, 716.909, 720.681, 
   726.12, 730.403, 733.553, 739.123, 742.156, 746.6, 750.027, 753.462, 757.426, 761.595, 
   764.336, 768.251, 772.371, 775.963, 778.886, 781.905, 784.798, 788.825, 792.372, 796.27, 
   800.361, 803.544, 806.544, 808.819, 812.146, 814.989, 819.234, 820.073, 824.067, 828.047, 
   830.277, 833.013, 835.374, 838.463, 840.82, 844.655, 846.391, 849.408, 851.659, 853.977, 
   856.409, 860.029, 862.128, 866.104, 866.864, 869.24, 872.133, 872.591, 876.528, 879.029, 
   880.786, 883.8, 886.065, 887.511, 890.301, 892.086, 894.429, 895.666, 897.961, 900.712, 
   901.559, 904.787, 906.882, 908.034, 911.366, 911.249, 914.274, 916.238, 918.105, 920.585, 
   920.473, 924.468, 923.888, 926.046, 928.648, 930.3, 931.861, 934.253, 934.081, 936.95, 
   938.319, 940.464, 940.539, 943.393, 944.729, 946.944, 947.712, 948.948, 951.026, 952.121, 
   954.114, 955.146, 956.206, 959.056, 960.316, 962.919, 961.946, 964.324, 966.134, 967.689, 
   968.612, 970.357, 972.302, 973.514, 976.512, 975.815, 979.043, 979.486, 981.285, 983.173, 
   983.96, 985.947, 987.447, 988.455, 991.739, 992.1, 993.045, 995.918, 997.377, 999.136, 
   1001.51, 1001.12, 1002.46, 1004.57, 1005.76, 1007.12, 1009.23, 1011.7, 1012.48, 1014.84, 
   1014.21, 1017.28, 1017.22, 1018.95, 1021.8, 1021.94, 1023.22, 1025.13, 1026.01, 1027.8, 
   1030.04, 1030.12, 1031.54, 1033.2, 1034.62, 1035.83, 1037.33, 1037.92, 1038.9, 1041.69 };
   return wid1[i];
}
 
double EvtD0To2pip2pim::anywid1300(double sc){
 
   double smin = (0.13957*3)*(0.13957*3);
   double dh = 0.01;
   int od  = (sc - 0.18)/dh;
   double sc_m = 0.18 + od*dh;
   double widuse = 0;
   if(sc>=0.18 && sc<=3.17){
      widuse = ((sc-sc_m)/dh)*(widT1300(od+1)-widT1300(od))+widT1300(od);
   }else if(sc<0.18 && sc>smin){
      widuse = ((sc - smin)/(0.18-smin))*widT1300(0);
   }else if(sc>3.17){
      widuse = widT1300(299);
   }else{
      widuse = 0;
   }
   return widuse;
}
 
complex<double> EvtD0To2pip2pim::RBWpi1300(double mx2, double mr, double wr){
 
   double mx  = sqrt(mx2);
   double mr2 = mr*mr;
   double g1 = wr/anywid1300(mr2);
   double wid0 = anywid1300(mx2)*g1;
 
   double denom_real = mr2-mx2;
   double denom_imag = mr*wid0;//real-i*imag;
 
   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;
 
}
 
/* build a1640 propagator */
double EvtD0To2pip2pim::widT1640(int i){
   double wid1[300]  = { 1.38316e-05, 0.000403892, 0.00181814, 0.0048161, 0.00982907, 0.0172548, 0.0273979, 0.040567, 0.0569061, 0.0768551, 
      0.100513, 0.128031, 0.159729, 0.195626, 0.236099, 0.280881, 0.330745, 0.386095, 0.446448, 0.511879, 
      0.583827, 0.66167, 0.745453, 0.835386, 0.934317, 1.0386, 1.1513, 1.26975, 1.39901, 1.53362, 
      1.68291, 1.84163, 2.0066, 2.18366, 2.37394, 2.57742, 2.7905, 3.02463, 3.27434, 3.53467, 
      3.80737, 4.10838, 4.41975, 4.76341, 5.12572, 5.51301, 5.91839, 6.36597, 6.8457, 7.33806, 
      7.87328, 8.45901, 9.08869, 9.74744, 10.464, 11.2096, 12.0103, 12.8556, 13.7563, 14.7352, 
      15.7336, 16.7432, 17.8117, 18.9327, 20.0186, 21.1632, 22.3549, 23.5172, 24.6518, 25.7808, 
      26.9103, 28.016, 29.1542, 30.0458, 31.0808, 32.1018, 33.0395, 33.9151, 34.8873, 35.7289, 
      36.5603, 37.2489, 38.023, 38.7983, 39.55, 40.2977, 40.8819, 41.4564, 42.1864, 42.7368, 
      43.3923, 43.8651, 44.4667, 44.8108, 45.3935, 45.9551, 46.2652, 46.8683, 47.1943, 47.6864, 
      48.1666, 48.5599, 48.8894, 49.1867, 49.6234, 49.9326, 50.4594, 50.6707, 51.005, 51.2612, 
      51.7638, 51.8946, 52.3176, 52.5107, 52.7378, 52.9418, 53.4019, 53.3571, 53.7937, 54.137, 
      54.2265, 54.3471, 54.6637, 54.897, 55.2174, 55.1577, 55.7098, 55.8616, 55.8862, 56.2106, 
      56.3357, 56.5165, 56.6819, 56.7906, 56.9814, 57.0507, 57.3059, 57.4898, 57.5848, 57.5792, 
      57.7696, 58.0302, 58.1915, 58.3319, 58.3892, 58.4671, 58.6736, 58.7872, 58.7949, 58.8366, 
      59.0247, 59.0881, 59.2675, 59.479, 59.6261, 59.6111, 59.6055, 59.7286, 59.8806, 60.0424, 
      60.1126, 60.0742, 60.2066, 60.2253, 60.565, 60.6557, 60.7359, 60.6405, 60.6429, 60.8521, 
      60.8098, 61.0699, 61.1678, 61.0329, 61.0522, 61.1792, 61.3671, 61.4394, 61.5152, 61.6122, 
      61.584, 61.711, 61.707, 61.7254, 61.816, 61.9248, 61.9748, 61.9498, 62.0014, 62.0634, 
      62.2929, 62.2349, 62.2101, 62.4434, 62.4281, 62.4166, 62.4905, 62.6055, 62.5097, 62.5994, 
      62.6637, 62.6794, 62.7068, 62.7908, 62.8135, 63.0085, 62.8848, 62.8159, 63.047, 62.8632, 
      63.1119, 63.0864, 63.1423, 63.2334, 63.0695, 63.2902, 63.3719, 63.1882, 63.2649, 63.3338, 
      63.4709, 63.4662, 63.3746, 63.623, 63.6402, 63.5632, 63.6611, 63.6012, 63.5904, 63.7467, 
      63.5535, 63.7792, 63.5213, 63.829, 63.8696, 63.8047, 63.9557, 63.9433, 63.9363, 63.9436, 
      63.9804, 64.0707, 64.0105, 63.96, 64.0437, 64.0235, 64.1795, 64.1377, 64.073, 64.2282, 
      64.2933, 64.4369, 64.3887, 64.2474, 64.2373, 64.3553, 64.425, 64.4401, 64.3197, 64.4212, 
      64.5787, 64.4919, 64.6878, 64.4998, 64.5788, 64.6628, 64.6658, 64.5072, 64.7227, 64.7327, 
      64.4472, 64.6792, 64.7801, 64.5715, 64.7263, 64.8505, 64.7488, 64.6448, 64.8962, 64.8815, 
      64.821, 64.902, 64.8944, 64.8959, 64.8957, 64.7882, 65.0725, 64.8787, 64.797, 65.1112, 
      65.1212, 65.157, 64.9412, 65.2601, 65.0662, 65.0093, 65.0899, 65.1035, 65.0865, 65.3276 };
   return wid1[i];
}
 
double EvtD0To2pip2pim::anywid1640(double sc){
 
   double smin = (0.13957*3)*(0.13957*3);
   double dh = 0.01;
   int od  = (sc - 0.18)/dh;
   double sc_m = 0.18 + od*dh;
   double widuse = 0;
   if(sc>=0.18 && sc<=3.17){
      widuse = ((sc-sc_m)/dh)*(widT1640(od+1)-widT1640(od))+widT1640(od);
   }else if(sc<0.18 && sc>smin){
      widuse = ((sc - smin)/(0.18-smin))*widT1640(0);
   }else if(sc>3.17){
      widuse = widT1640(299);
   }else{
      widuse = 0;
   }
   return widuse;
}
 
complex<double> EvtD0To2pip2pim::RBWa1640(double mx2, double mr, double wr){
 
   double mx  = sqrt(mx2);
   double mr2 = mr*mr;
   double g1 = wr/anywid1640(mr2);
   double wid0 = anywid1640(mx2)*g1;
 
   double denom_real = mr2-mx2;
   double denom_imag = mr*wid0;//real-i*imag;
 
   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;
 
}
 
// PiPi-S wave K-Matrix
double EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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 EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::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> EvtD0To2pip2pim::CalD0Amp(){
   return Amp(m_Pip1, m_Pim1, m_Pip2, m_Pim2); 
}
complex<double> EvtD0To2pip2pim::CalDbAmp(){
 
   vector<double> cpPip1; cpPip1.clear();
   vector<double> cpPip2; cpPip2.clear();
   vector<double> cpPim1; cpPim1.clear();
   vector<double> cpPim2; cpPim2.clear();
 
   cpPip1.push_back(-m_Pim1[0]); cpPim1.push_back(-m_Pip1[0]); cpPip2.push_back(-m_Pim2[0]); cpPim2.push_back(-m_Pip2[0]);
   cpPip1.push_back(-m_Pim1[1]); cpPim1.push_back(-m_Pip1[1]); cpPip2.push_back(-m_Pim2[1]); cpPim2.push_back(-m_Pip2[1]);
   cpPip1.push_back(-m_Pim1[2]); cpPim1.push_back(-m_Pip1[2]); cpPip2.push_back(-m_Pim2[2]); cpPim2.push_back(-m_Pip2[2]);
   cpPip1.push_back(m_Pim1[3]);  cpPim1.push_back(m_Pip1[3]);  cpPip2.push_back(m_Pim2[3]);  cpPim2.push_back(m_Pip2[3]);
 
   return Amp(cpPip1, cpPim1, cpPip2, cpPim2); 
}
 
complex<double> EvtD0To2pip2pim::Amp(vector<double> Pip1, vector<double> Pim1, vector<double> Pip2, vector<double> Pim2)
{
 
   vector<double> Pip1Pim1; Pip1Pim1.clear();
   vector<double> Pip1Pim2; Pip1Pim2.clear();
   vector<double> Pip2Pim1; Pip2Pim1.clear();
   vector<double> Pip2Pim2; Pip2Pim2.clear();
 
   Pip1Pim1 = sum_tensor(Pip1, Pim1);
   Pip1Pim2 = sum_tensor(Pip1, Pim2);
   Pip2Pim1 = sum_tensor(Pip2, Pim1);
   Pip2Pim2 = sum_tensor(Pip2, Pim2);
 
   vector<double> Pip1Pip2Pim1; Pip1Pip2Pim1.clear();
   vector<double> Pip1Pip2Pim2; Pip1Pip2Pim2.clear();
   vector<double> Pim1Pim2Pip1; Pim1Pim2Pip1.clear();
   vector<double> Pim1Pim2Pip2; Pim1Pim2Pip2.clear();
 
   Pip1Pip2Pim1 = sum_tensor(Pip1Pim1, Pip2);
   Pip1Pip2Pim2 = sum_tensor(Pip1Pim2, Pip2);
   Pim1Pim2Pip1 = sum_tensor(Pip1Pim1, Pim2);
   Pim1Pim2Pip2 = sum_tensor(Pip2Pim1, Pim2);
 
   vector<double> D0; D0.clear();
   D0 = sum_tensor(Pip1Pip2Pim1, Pim2);
 
   double M2_Pip1Pim1 = contract_11_0(Pip1Pim1, Pip1Pim1);
   double M2_Pip1Pim2 = contract_11_0(Pip1Pim2, Pip1Pim2);
   double M2_Pip2Pim1 = contract_11_0(Pip2Pim1, Pip2Pim1);
   double M2_Pip2Pim2 = contract_11_0(Pip2Pim2, Pip2Pim2);
 
   double M2_Pip1Pip2Pim1 = contract_11_0(Pip1Pip2Pim1, Pip1Pip2Pim1);
   double M2_Pip1Pip2Pim2 = contract_11_0(Pip1Pip2Pim2, Pip1Pip2Pim2);
   double M2_Pim1Pim2Pip1 = contract_11_0(Pim1Pim2Pip1, Pim1Pim2Pip1);
   double M2_Pim1Pim2Pip2 = contract_11_0(Pim1Pim2Pip2, Pim1Pim2Pip2);
   double M2_D0 = contract_11_0(D0, D0);
 
   complex<double> GS_rho770_11 = GS(M2_Pip1Pim1, m0_rho770, w0_rho770, m2_Pi, m2_Pi, rRes, 1);
   complex<double> GS_rho770_12 = GS(M2_Pip1Pim2, m0_rho770, w0_rho770, m2_Pi, m2_Pi, rRes, 1);
   complex<double> GS_rho770_21 = GS(M2_Pip2Pim1, m0_rho770, w0_rho770, m2_Pi, m2_Pi, rRes, 1);
   complex<double> GS_rho770_22 = GS(M2_Pip2Pim2, m0_rho770, w0_rho770, m2_Pi, m2_Pi, rRes, 1);
 
   complex<double> GS_rho1450_11 = GS(M2_Pip1Pim1, m0_rho1450, w0_rho1450, m2_Pi, m2_Pi, rRes, 1);
   complex<double> GS_rho1450_12 = GS(M2_Pip1Pim2, m0_rho1450, w0_rho1450, m2_Pi, m2_Pi, rRes, 1);
   complex<double> GS_rho1450_21 = GS(M2_Pip2Pim1, m0_rho1450, w0_rho1450, m2_Pi, m2_Pi, rRes, 1);
   complex<double> GS_rho1450_22 = GS(M2_Pip2Pim2, m0_rho1450, w0_rho1450, m2_Pi, m2_Pi, rRes, 1);
 
   complex<double> RBW_f21270_11 = RBW(M2_Pip1Pim1, m0_f21270, w0_f21270, m2_Pi, m2_Pi, rRes, 2);
   complex<double> RBW_f21270_12 = RBW(M2_Pip1Pim2, m0_f21270, w0_f21270, m2_Pi, m2_Pi, rRes, 2);
   complex<double> RBW_f21270_21 = RBW(M2_Pip2Pim1, m0_f21270, w0_f21270, m2_Pi, m2_Pi, rRes, 2);
   complex<double> RBW_f21270_22 = RBW(M2_Pip2Pim2, m0_f21270, w0_f21270, m2_Pi, m2_Pi, rRes, 2);
 
   complex<double> PiPiS_11_0 = Fvector(M2_Pip1Pim1, s0_prod, 0);
   complex<double> PiPiS_12_0 = Fvector(M2_Pip1Pim2, s0_prod, 0);
   complex<double> PiPiS_21_0 = Fvector(M2_Pip2Pim1, s0_prod, 0);
   complex<double> PiPiS_22_0 = Fvector(M2_Pip2Pim2, s0_prod, 0);
 
   complex<double> PiPiS_11_1 = Fvector(M2_Pip1Pim1, s0_prod, 1);
   complex<double> PiPiS_12_1 = Fvector(M2_Pip1Pim2, s0_prod, 1);
   complex<double> PiPiS_21_1 = Fvector(M2_Pip2Pim1, s0_prod, 1);
   complex<double> PiPiS_22_1 = Fvector(M2_Pip2Pim2, s0_prod, 1);
 
   complex<double> PiPiS_11_5 = Fvector(M2_Pip1Pim1, s0_prod, 5);
   complex<double> PiPiS_12_5 = Fvector(M2_Pip1Pim2, s0_prod, 5);
   complex<double> PiPiS_21_5 = Fvector(M2_Pip2Pim1, s0_prod, 5);
   complex<double> PiPiS_22_5 = Fvector(M2_Pip2Pim2, s0_prod, 5);
 
   complex<double> PiPiS_11_6 = Fvector(M2_Pip1Pim1, s0_prod, 6);
   complex<double> PiPiS_12_6 = Fvector(M2_Pip1Pim2, s0_prod, 6);
   complex<double> PiPiS_21_6 = Fvector(M2_Pip2Pim1, s0_prod, 6);
   complex<double> PiPiS_22_6 = Fvector(M2_Pip2Pim2, s0_prod, 6);
 
   complex<double> RBW_a11260p_1 = RBWa1260(M2_Pip1Pip2Pim1, m0_a11260, g1_a11260, g2_a11260);
   complex<double> RBW_a11260p_2 = RBWa1260(M2_Pip1Pip2Pim2, m0_a11260, g1_a11260, g2_a11260);
   complex<double> RBW_a11260m_1 = RBWa1260(M2_Pim1Pim2Pip1, m0_a11260, g1_a11260, g2_a11260);
   complex<double> RBW_a11260m_2 = RBWa1260(M2_Pim1Pim2Pip2, m0_a11260, g1_a11260, g2_a11260);
 
   complex<double> RBW_a21320p_1 = RBW(M2_Pip1Pip2Pim1, m0_a21320, w0_a21320, -1.0, -1.0, -1, -1);
   complex<double> RBW_a21320p_2 = RBW(M2_Pip1Pip2Pim2, m0_a21320, w0_a21320, -1.0, -1.0, -1, -1);
   complex<double> RBW_a21320m_1 = RBW(M2_Pim1Pim2Pip1, m0_a21320, w0_a21320, -1.0, -1.0, -1, -1);
   complex<double> RBW_a21320m_2 = RBW(M2_Pim1Pim2Pip2, m0_a21320, w0_a21320, -1.0, -1.0, -1, -1);
 
   complex<double> RBW_pi1300p_1 = RBWpi1300(M2_Pip1Pip2Pim1, m0_pi1300, w0_pi1300);
   complex<double> RBW_pi1300p_2 = RBWpi1300(M2_Pip1Pip2Pim2, m0_pi1300, w0_pi1300);
   complex<double> RBW_pi1300m_1 = RBWpi1300(M2_Pim1Pim2Pip1, m0_pi1300, w0_pi1300);
   complex<double> RBW_pi1300m_2 = RBWpi1300(M2_Pim1Pim2Pip2, m0_pi1300, w0_pi1300);
 
   complex<double> RBW_a11420p_1 = RBW(M2_Pip1Pip2Pim1, m0_a11420, w0_a11420, -1.0, -1.0, -1, -1);
   complex<double> RBW_a11420p_2 = RBW(M2_Pip1Pip2Pim2, m0_a11420, w0_a11420, -1.0, -1.0, -1, -1);
   complex<double> RBW_a11420m_1 = RBW(M2_Pim1Pim2Pip1, m0_a11420, w0_a11420, -1.0, -1.0, -1, -1);
   complex<double> RBW_a11420m_2 = RBW(M2_Pim1Pim2Pip2, m0_a11420, w0_a11420, -1.0, -1.0, -1, -1);
 
   // D->XP 1-rank Projection
   vector<double> Proj1_3p1; Proj1_3p1.clear();
   vector<double> Proj1_3p2; Proj1_3p2.clear();
   vector<double> Proj1_3m1; Proj1_3m1.clear();
   vector<double> Proj1_3m2; Proj1_3m2.clear();
 
   Proj1_3p1 = ProjectionTensors(Pip1Pip2Pim1,1);
   Proj1_3p2 = ProjectionTensors(Pip1Pip2Pim2,1);
   Proj1_3m1 = ProjectionTensors(Pim1Pim2Pip1,1);
   Proj1_3m2 = ProjectionTensors(Pim1Pim2Pip2,1);
 
   // D->XP 2-rank Projection   
   vector<double> Proj2_3p1; Proj2_3p1.clear();
   vector<double> Proj2_3p2; Proj2_3p2.clear();
   vector<double> Proj2_3m1; Proj2_3m1.clear();
   vector<double> Proj2_3m2; Proj2_3m2.clear();
 
   Proj2_3p1 = ProjectionTensors(Pip1Pip2Pim1,2);
   Proj2_3p2 = ProjectionTensors(Pip1Pip2Pim2,2);
   Proj2_3m1 = ProjectionTensors(Pim1Pim2Pip1,2);
   Proj2_3m2 = ProjectionTensors(Pim1Pim2Pip2,2);
 
   // X->PP Orbital
   vector<double> T1_Pip1Pim1; T1_Pip1Pim1.clear();
   vector<double> T1_Pip1Pim2; T1_Pip1Pim2.clear();
   vector<double> T1_Pip2Pim1; T1_Pip2Pim1.clear();
   vector<double> T1_Pip2Pim2; T1_Pip2Pim1.clear();
 
   T1_Pip1Pim1 = OrbitalTensors(Pip1Pim1, Pip1, Pim1, rRes, 1);
   T1_Pip1Pim2 = OrbitalTensors(Pip1Pim2, Pip1, Pim2, rRes, 1);
   T1_Pip2Pim1 = OrbitalTensors(Pip2Pim1, Pip2, Pim1, rRes, 1);
   T1_Pip2Pim2 = OrbitalTensors(Pip2Pim2, Pip2, Pim2, rRes, 1);
 
   vector<double> T1_Pim1Pip1; T1_Pim1Pip1.clear();
   vector<double> T1_Pim1Pip2; T1_Pim1Pip2.clear();
   vector<double> T1_Pim2Pip1; T1_Pim2Pip1.clear();
   vector<double> T1_Pim2Pip2; T1_Pim2Pip2.clear();
 
   T1_Pim1Pip1 = OrbitalTensors(Pip1Pim1, Pim1, Pip1, rRes, 1);
   T1_Pim1Pip2 = OrbitalTensors(Pip2Pim1, Pim1, Pip2, rRes, 1);
   T1_Pim2Pip1 = OrbitalTensors(Pip1Pim2, Pim2, Pip1, rRes, 1);
   T1_Pim2Pip2 = OrbitalTensors(Pip2Pim2, Pim2, Pip2, rRes, 1);
 
   vector<double> T2_Pip1Pim1; T2_Pip1Pim1.clear();
   vector<double> T2_Pip1Pim2; T2_Pip1Pim2.clear();
   vector<double> T2_Pip2Pim1; T2_Pip2Pim1.clear();
   vector<double> T2_Pip2Pim2; T2_Pip2Pim1.clear();
 
   T2_Pip1Pim1 = OrbitalTensors(Pip1Pim1, Pip1, Pim1, rRes, 2);
   T2_Pip1Pim2 = OrbitalTensors(Pip1Pim2, Pip1, Pim2, rRes, 2);
   T2_Pip2Pim1 = OrbitalTensors(Pip2Pim1, Pip2, Pim1, rRes, 2);
   T2_Pip2Pim2 = OrbitalTensors(Pip2Pim2, Pip2, Pim2, rRes, 2);
 
   // X->YP Orbital
   vector<double> T1_Pip1Pim1Pip2; T1_Pip1Pim1Pip2.clear();
   vector<double> T1_Pip2Pim1Pip1; T1_Pip2Pim1Pip1.clear();
   vector<double> T1_Pip1Pim2Pip2; T1_Pip1Pim2Pip2.clear();
   vector<double> T1_Pip2Pim2Pip1; T1_Pip2Pim2Pip1.clear();
   vector<double> T1_Pip1Pim1Pim2; T1_Pip1Pim1Pim2.clear();
   vector<double> T1_Pip1Pim2Pim1; T1_Pip1Pim2Pim1.clear();
   vector<double> T1_Pip2Pim1Pim2; T1_Pip2Pim1Pim2.clear();
   vector<double> T1_Pip2Pim2Pim1; T1_Pip2Pim2Pim1.clear();
 
   T1_Pip1Pim1Pip2 = OrbitalTensors(Pip1Pip2Pim1, Pip1Pim1, Pip2, rRes, 1);
   T1_Pip2Pim1Pip1 = OrbitalTensors(Pip1Pip2Pim1, Pip2Pim1, Pip1, rRes, 1);
   T1_Pip1Pim2Pip2 = OrbitalTensors(Pip1Pip2Pim2, Pip1Pim2, Pip2, rRes, 1);
   T1_Pip2Pim2Pip1 = OrbitalTensors(Pip1Pip2Pim2, Pip2Pim2, Pip1, rRes, 1);
   T1_Pip1Pim1Pim2 = OrbitalTensors(Pim1Pim2Pip1, Pip1Pim1, Pim2, rRes, 1);
   T1_Pip2Pim1Pim2 = OrbitalTensors(Pim1Pim2Pip2, Pip2Pim1, Pim2, rRes, 1);
   T1_Pip1Pim2Pim1 = OrbitalTensors(Pim1Pim2Pip1, Pip1Pim2, Pim1, rRes, 1);
   T1_Pip2Pim2Pim1 = OrbitalTensors(Pim1Pim2Pip2, Pip2Pim2, Pim1, rRes, 1);
 
   vector<double> T2_Pip1Pim1Pip2; T2_Pip1Pim1Pip2.clear();
   vector<double> T2_Pip2Pim1Pip1; T2_Pip2Pim1Pip1.clear();
   vector<double> T2_Pip1Pim2Pip2; T2_Pip1Pim2Pip2.clear();
   vector<double> T2_Pip2Pim2Pip1; T2_Pip2Pim2Pip1.clear();
   vector<double> T2_Pip1Pim1Pim2; T2_Pip1Pim1Pim2.clear();
   vector<double> T2_Pip2Pim1Pim2; T2_Pip2Pim1Pim2.clear();
   vector<double> T2_Pip1Pim2Pim1; T2_Pip1Pim2Pim1.clear();
   vector<double> T2_Pip2Pim2Pim1; T2_Pip2Pim2Pim1.clear();
 
   T2_Pip1Pim1Pip2 = OrbitalTensors(Pip1Pip2Pim1, Pip1Pim1, Pip2, rRes, 2);
   T2_Pip2Pim1Pip1 = OrbitalTensors(Pip1Pip2Pim1, Pip2Pim1, Pip1, rRes, 2);
   T2_Pip1Pim2Pip2 = OrbitalTensors(Pip1Pip2Pim2, Pip1Pim2, Pip2, rRes, 2);
   T2_Pip2Pim2Pip1 = OrbitalTensors(Pip1Pip2Pim2, Pip2Pim2, Pip1, rRes, 2);
   T2_Pip1Pim1Pim2 = OrbitalTensors(Pim1Pim2Pip1, Pip1Pim1, Pim2, rRes, 2);
   T2_Pip2Pim1Pim2 = OrbitalTensors(Pim1Pim2Pip2, Pip2Pim1, Pim2, rRes, 2);
   T2_Pip1Pim2Pim1 = OrbitalTensors(Pim1Pim2Pip1, Pip1Pim2, Pim1, rRes, 2);
   T2_Pip2Pim2Pim1 = OrbitalTensors(Pim1Pim2Pip2, Pip2Pim2, Pim1, rRes, 2);
 
   // D->XX Orbital
   vector<double> T1_2z11; T1_2z11.clear();
   vector<double> T1_2z12; T1_2z12.clear();
   vector<double> T1_2z21; T1_2z21.clear();
   vector<double> T1_2z22; T1_2z22.clear();
 
   T1_2z11 = OrbitalTensors(D0, Pip1Pim1, Pip2Pim2, rD, 1);
   T1_2z12 = OrbitalTensors(D0, Pip2Pim2, Pip1Pim1, rD, 1);
   T1_2z21 = OrbitalTensors(D0, Pip1Pim2, Pip2Pim1, rD, 1);
   T1_2z22 = OrbitalTensors(D0, Pip2Pim1, Pip1Pim2, rD, 1);
 
   vector<double> T2_2z11; T2_2z11.clear();
   vector<double> T2_2z12; T2_2z12.clear();
   vector<double> T2_2z21; T2_2z21.clear();
   vector<double> T2_2z22; T2_2z22.clear();
 
   T2_2z11 = OrbitalTensors(D0, Pip1Pim1, Pip2Pim2, rD, 2);
   T2_2z12 = OrbitalTensors(D0, Pip2Pim2, Pip1Pim1, rD, 2);
   T2_2z21 = OrbitalTensors(D0, Pip1Pim2, Pip2Pim1, rD, 2);
   T2_2z22 = OrbitalTensors(D0, Pip2Pim1, Pip1Pim2, rD, 2);
 
   // D->XP Orbital
   vector<double> T1_3p1; T1_3p1.clear();
   vector<double> T1_3p2; T1_3p2.clear();
   vector<double> T1_3m1; T1_3m1.clear();
   vector<double> T1_3m2; T1_3m2.clear();
 
   T1_3p1 = OrbitalTensors(D0, Pip1Pip2Pim1, Pim2, rD, 1);
   T1_3p2 = OrbitalTensors(D0, Pip1Pip2Pim2, Pim1, rD, 1);
   T1_3m1 = OrbitalTensors(D0, Pim1Pim2Pip1, Pip2, rD, 1);
   T1_3m2 = OrbitalTensors(D0, Pim1Pim2Pip2, Pip1, rD, 1);
 
   vector<double> T2_3p1; T2_3p1.clear();
   vector<double> T2_3p2; T2_3p2.clear();
   vector<double> T2_3m1; T2_3m1.clear();
   vector<double> T2_3m2; T2_3m2.clear();
 
   T2_3p1 = OrbitalTensors(D0, Pip1Pip2Pim1, Pim2, rD, 2);
   T2_3p2 = OrbitalTensors(D0, Pip1Pip2Pim2, Pim1, rD, 2);
   T2_3m1 = OrbitalTensors(D0, Pim1Pim2Pip1, Pip2, rD, 2);
   T2_3m2 = OrbitalTensors(D0, Pim1Pim2Pip2, Pip1, rD, 2);
 
   complex<double> amplitude(0,0);
 
   // D0 -> a1(1260)+ {rho(770)0 pi+[S]} pi-
   double SF_Ap_S_VP_1 = contract_11_0(contract_21_1(Proj1_3p1, T1_Pip2Pim1), T1_3p1);
   double SF_Ap_S_VP_2 = contract_11_0(contract_21_1(Proj1_3p1, T1_Pip1Pim1), T1_3p1);
   double SF_Ap_S_VP_3 = contract_11_0(contract_21_1(Proj1_3p2, T1_Pip2Pim2), T1_3p2);
   double SF_Ap_S_VP_4 = contract_11_0(contract_21_1(Proj1_3p2, T1_Pip1Pim2), T1_3p2);
   amplitude += fitpara[0]*(SF_Ap_S_VP_1*RBW_a11260p_1*GS_rho770_21 + SF_Ap_S_VP_2*RBW_a11260p_1*GS_rho770_11 + SF_Ap_S_VP_3*RBW_a11260p_2*GS_rho770_22 + SF_Ap_S_VP_4*RBW_a11260p_2*GS_rho770_12);
 
   // D0 -> a1(1260)+ {rho(770)0 pi+[D]} pi-
   double SF_Ap_D_VP_1 = contract_11_0(contract_21_1(T2_Pip2Pim1Pip1, T1_Pip2Pim1), T1_3p1);
   double SF_Ap_D_VP_2 = contract_11_0(contract_21_1(T2_Pip1Pim1Pip2, T1_Pip1Pim1), T1_3p1);
   double SF_Ap_D_VP_3 = contract_11_0(contract_21_1(T2_Pip2Pim2Pip1, T1_Pip2Pim2), T1_3p2);
   double SF_Ap_D_VP_4 = contract_11_0(contract_21_1(T2_Pip1Pim2Pip2, T1_Pip1Pim2), T1_3p2);
 
   amplitude += fitpara[1]*(SF_Ap_D_VP_1*RBW_a11260p_1*GS_rho770_21 + SF_Ap_D_VP_2*RBW_a11260p_1*GS_rho770_11 + SF_Ap_D_VP_3*RBW_a11260p_2*GS_rho770_22 + SF_Ap_D_VP_4*RBW_a11260p_2*GS_rho770_12);
 
   // D0 -> a1(1260)+ {f2(1270) pi+ [P]} pi-
   double SF_Ap_P_TP_1 = contract_11_0(contract_21_1(contract_42_2(Proj2_3p1, T2_Pip2Pim1), T1_Pip2Pim1Pip1), T1_3p1);
   double SF_Ap_P_TP_2 = contract_11_0(contract_21_1(contract_42_2(Proj2_3p1, T2_Pip1Pim1), T1_Pip1Pim1Pip2), T1_3p1);
   double SF_Ap_P_TP_3 = contract_11_0(contract_21_1(contract_42_2(Proj2_3p2, T2_Pip2Pim2), T1_Pip2Pim2Pip1), T1_3p2);
   double SF_Ap_P_TP_4 = contract_11_0(contract_21_1(contract_42_2(Proj2_3p2, T2_Pip1Pim2), T1_Pip1Pim2Pip2), T1_3p2);
 
   amplitude += fitpara[2]*(SF_Ap_P_TP_1*RBW_a11260p_1*RBW_f21270_21 + SF_Ap_P_TP_2*RBW_a11260p_1*RBW_f21270_11 + SF_Ap_P_TP_3*RBW_a11260p_2*RBW_f21270_22 + SF_Ap_P_TP_4*RBW_a11260p_2*RBW_f21270_12);
 
   // D0 -> a1(1260)+ {f0 pi+ [P]} pi-
   double SF_Ap_P_SP_1 = contract_11_0(T1_3p1, T1_Pip2Pim1Pip1);
   double SF_Ap_P_SP_2 = contract_11_0(T1_3p1, T1_Pip1Pim1Pip2);
   double SF_Ap_P_SP_3 = contract_11_0(T1_3p2, T1_Pip2Pim2Pip1);
   double SF_Ap_P_SP_4 = contract_11_0(T1_3p2, T1_Pip1Pim2Pip2);
 
   amplitude += fitpara[3]*(SF_Ap_P_SP_1*RBW_a11260p_1*PiPiS_21_0 + SF_Ap_P_SP_2*RBW_a11260p_1*PiPiS_11_0 + SF_Ap_P_SP_3*RBW_a11260p_2*PiPiS_22_0 + SF_Ap_P_SP_4*RBW_a11260p_2*PiPiS_12_0);
   amplitude += fitpara[4]*(SF_Ap_P_SP_1*RBW_a11260p_1*PiPiS_21_1 + SF_Ap_P_SP_2*RBW_a11260p_1*PiPiS_11_1 + SF_Ap_P_SP_3*RBW_a11260p_2*PiPiS_22_1 + SF_Ap_P_SP_4*RBW_a11260p_2*PiPiS_12_1);
   amplitude += fitpara[5]*(SF_Ap_P_SP_1*RBW_a11260p_1*PiPiS_21_5 + SF_Ap_P_SP_2*RBW_a11260p_1*PiPiS_11_5 + SF_Ap_P_SP_3*RBW_a11260p_2*PiPiS_22_5 + SF_Ap_P_SP_4*RBW_a11260p_2*PiPiS_12_5);
 
   // D0->a1(1260)- {rho(770)0 pi- [S]} pi+
   double SF_Am_S_VP_1 = contract_11_0(contract_21_1(Proj1_3m1, T1_Pim2Pip1), T1_3m1);
   double SF_Am_S_VP_2 = contract_11_0(contract_21_1(Proj1_3m1, T1_Pim1Pip1), T1_3m1);
   double SF_Am_S_VP_3 = contract_11_0(contract_21_1(Proj1_3m2, T1_Pim2Pip2), T1_3m2);
   double SF_Am_S_VP_4 = contract_11_0(contract_21_1(Proj1_3m2, T1_Pim1Pip2), T1_3m2);
 
   amplitude += fitpara[0]*fitpara[6]*(SF_Am_S_VP_1*RBW_a11260m_1*GS_rho770_12 + SF_Am_S_VP_2*RBW_a11260m_1*GS_rho770_11 + SF_Am_S_VP_3*RBW_a11260m_2*GS_rho770_22 + SF_Am_S_VP_4*RBW_a11260m_2*GS_rho770_21);
 
   // D0 -> a1(1260)- {rho(770)0 pi-[D]} pi+
   double SF_Am_D_VP_1 = contract_11_0(contract_21_1(T2_Pip1Pim2Pim1, T1_Pim2Pip1), T1_3m1);
   double SF_Am_D_VP_2 = contract_11_0(contract_21_1(T2_Pip1Pim1Pim2, T1_Pim1Pip1), T1_3m1);
   double SF_Am_D_VP_3 = contract_11_0(contract_21_1(T2_Pip2Pim2Pim1, T1_Pim2Pip2), T1_3m2);
   double SF_Am_D_VP_4 = contract_11_0(contract_21_1(T2_Pip2Pim1Pim2, T1_Pim1Pip2), T1_3m2);
 
   amplitude += fitpara[1]*fitpara[6]*(SF_Am_D_VP_1*RBW_a11260m_1*GS_rho770_12 + SF_Am_D_VP_2*RBW_a11260m_1*GS_rho770_11 + SF_Am_D_VP_3*RBW_a11260m_2*GS_rho770_22 + SF_Am_D_VP_4*RBW_a11260m_2*GS_rho770_21);
 
   // D0 -> a1(1260)- {f2(1270) pi- [P]} pi+
   double SF_Am_P_TP_1 = contract_11_0(contract_21_1(contract_42_2(Proj2_3m1, T2_Pip1Pim2), T1_Pip1Pim2Pim1), T1_3m1);
   double SF_Am_P_TP_2 = contract_11_0(contract_21_1(contract_42_2(Proj2_3m1, T2_Pip1Pim1), T1_Pip1Pim1Pim2), T1_3m1);
   double SF_Am_P_TP_3 = contract_11_0(contract_21_1(contract_42_2(Proj2_3m2, T2_Pip2Pim2), T1_Pip2Pim2Pim1), T1_3m2);
   double SF_Am_P_TP_4 = contract_11_0(contract_21_1(contract_42_2(Proj2_3m2, T2_Pip2Pim1), T1_Pip2Pim1Pim2), T1_3m2);
 
   amplitude += fitpara[2]*fitpara[6]*(SF_Am_P_TP_1*RBW_a11260m_1*RBW_f21270_12 + SF_Am_P_TP_2*RBW_a11260m_1*RBW_f21270_11 + SF_Am_P_TP_3*RBW_a11260m_2*RBW_f21270_22 + SF_Am_P_TP_4*RBW_a11260m_2*RBW_f21270_21);
 
   // D0 -> a1(1260)- {f0 pi- [P]} pi+
   double SF_Am_P_SP_1 = contract_11_0(T1_3m1, T1_Pip1Pim2Pim1);
   double SF_Am_P_SP_2 = contract_11_0(T1_3m1, T1_Pip1Pim1Pim2);
   double SF_Am_P_SP_3 = contract_11_0(T1_3m2, T1_Pip2Pim2Pim1);
   double SF_Am_P_SP_4 = contract_11_0(T1_3m2, T1_Pip2Pim1Pim2);
 
   amplitude += fitpara[3]*fitpara[6]*(SF_Am_P_SP_1*RBW_a11260m_1*PiPiS_12_0 + SF_Am_P_SP_2*RBW_a11260m_1*PiPiS_11_0 + SF_Am_P_SP_3*RBW_a11260m_2*PiPiS_22_0 + SF_Am_P_SP_4*RBW_a11260m_2*PiPiS_21_0);
   amplitude += fitpara[4]*fitpara[6]*(SF_Am_P_SP_1*RBW_a11260m_1*PiPiS_12_1 + SF_Am_P_SP_2*RBW_a11260m_1*PiPiS_11_1 + SF_Am_P_SP_3*RBW_a11260m_2*PiPiS_22_1 + SF_Am_P_SP_4*RBW_a11260m_2*PiPiS_21_1);
   amplitude += fitpara[5]*fitpara[6]*(SF_Am_P_SP_1*RBW_a11260m_1*PiPiS_12_5 + SF_Am_P_SP_2*RBW_a11260m_1*PiPiS_11_5 + SF_Am_P_SP_3*RBW_a11260m_2*PiPiS_22_5 + SF_Am_P_SP_4*RBW_a11260m_2*PiPiS_21_5);
 
   // D0 -> a1(1420)+ {rho(770)0 pi+[S]} pi-
//   amplitude += fitpara[7]*(SF_Ap_S_VP_1*RBW_a11420p_1*GS_rho770_21 + SF_Ap_S_VP_2*RBW_a11420p_1*GS_rho770_11 + SF_Ap_S_VP_3*RBW_a11420p_2*GS_rho770_22 + SF_Ap_S_VP_4*RBW_a11420p_2*GS_rho770_12);
 
   // D0 -> a1(1420)+ {f0 pi+[P]} pi-
   amplitude += fitpara[7]*(SF_Ap_P_SP_1*RBW_a11420p_1*PiPiS_21_5 + SF_Ap_P_SP_2*RBW_a11420p_1*PiPiS_11_5 + SF_Ap_P_SP_3*RBW_a11420p_2*PiPiS_22_5 + SF_Ap_P_SP_4*RBW_a11420p_2*PiPiS_12_5);
   amplitude += fitpara[8]*(SF_Ap_P_SP_1*RBW_a11420p_1*PiPiS_21_6 + SF_Ap_P_SP_2*RBW_a11420p_1*PiPiS_11_6 + SF_Ap_P_SP_3*RBW_a11420p_2*PiPiS_22_6 + SF_Ap_P_SP_4*RBW_a11420p_2*PiPiS_12_6);
 
   // D0 -> a2(1320)+ {rho(770)0 pi+[D]} pi-
   double SF_Tp_D_VP_1 = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3p1, T1_Pip2Pim1)), Pip1Pip2Pim1), contract_42_2(Proj2_3p1, T2_3p1)), T2_Pip2Pim1Pip1);
   double SF_Tp_D_VP_2 = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3p1, T1_Pip1Pim1)), Pip1Pip2Pim1), contract_42_2(Proj2_3p1, T2_3p1)), T2_Pip1Pim1Pip2);
   double SF_Tp_D_VP_3 = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3p2, T1_Pip2Pim2)), Pip1Pip2Pim2), contract_42_2(Proj2_3p2, T2_3p2)), T2_Pip2Pim2Pip1);
   double SF_Tp_D_VP_4 = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3p2, T1_Pip1Pim2)), Pip1Pip2Pim2), contract_42_2(Proj2_3p2, T2_3p2)), T2_Pip1Pim2Pip2);
 
   amplitude += fitpara[9]*(SF_Tp_D_VP_1*RBW_a21320p_1*GS_rho770_21 + SF_Tp_D_VP_2*RBW_a21320p_1*GS_rho770_11 + SF_Tp_D_VP_3*RBW_a21320p_2*GS_rho770_22 + SF_Tp_D_VP_4*RBW_a21320p_2*GS_rho770_12);
 
   // D0 -> a2(1320)- {rho(770)0 pi-[D]} pi+
   double SF_Tm_D_VP_1 = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3m1, T1_Pim2Pip1)), Pim1Pim2Pip1), contract_42_2(Proj2_3m1, T2_3m1)), T2_Pip1Pim2Pim1);
   double SF_Tm_D_VP_2 = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3m1, T1_Pim1Pip1)), Pim1Pim2Pip1), contract_42_2(Proj2_3m1, T2_3m1)), T2_Pip1Pim1Pim2);
   double SF_Tm_D_VP_3 = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3m2, T1_Pim2Pip2)), Pim1Pim2Pip2), contract_42_2(Proj2_3m2, T2_3m2)), T2_Pip2Pim2Pim1);
   double SF_Tm_D_VP_4 = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3m2, T1_Pim1Pip2)), Pim1Pim2Pip2), contract_42_2(Proj2_3m2, T2_3m2)), T2_Pip2Pim1Pim2);
 
   amplitude += fitpara[10]*(SF_Tm_D_VP_1*RBW_a21320m_1*GS_rho770_12 + SF_Tm_D_VP_2*RBW_a21320m_1*GS_rho770_11 + SF_Tm_D_VP_3*RBW_a21320m_2*GS_rho770_22 + SF_Tm_D_VP_4*RBW_a21320m_2*GS_rho770_21);
 
   // D0 -> pi(1300)- {rho(770)0 pi-} pi+
   double SF_Pm_P_VP_1 = contract_11_0(T1_Pim2Pip1,T1_Pip1Pim2Pim1);
   double SF_Pm_P_VP_2 = contract_11_0(T1_Pim1Pip1,T1_Pip1Pim1Pim2);
   double SF_Pm_P_VP_3 = contract_11_0(T1_Pim2Pip2,T1_Pip2Pim2Pim1);
   double SF_Pm_P_VP_4 = contract_11_0(T1_Pim1Pip2,T1_Pip2Pim1Pim2);
 
   amplitude += fitpara[11]*(SF_Pm_P_VP_1*GS_rho770_12*RBW_pi1300m_1 + SF_Pm_P_VP_2*GS_rho770_11*RBW_pi1300m_1 + SF_Pm_P_VP_3*GS_rho770_22*RBW_pi1300m_2 + SF_Pm_P_VP_4*GS_rho770_21*RBW_pi1300m_2);
/*
   // D0 -> pi(1300)- {f2(1270)0 pi-} pi+
   double SF_Pm_D_TP_1 = contract_22_0(T2_Pip1Pim2,T2_Pip1Pim2Pim1);
   double SF_Pm_D_TP_2 = contract_22_0(T2_Pip1Pim1,T2_Pip1Pim1Pim2);
   double SF_Pm_D_TP_3 = contract_22_0(T2_Pip2Pim2,T2_Pip2Pim2Pim1);
   double SF_Pm_D_TP_4 = contract_22_0(T2_Pip2Pim1,T2_Pip2Pim1Pim2);
 
   amplitude += fitpara[12]*fitpara[11]*(SF_Pm_D_TP_1*RBW_f21270_12*RBW_pi1300m_1 + SF_Pm_D_TP_2*RBW_f21270_11*RBW_pi1300m_1 + SF_Pm_D_TP_3*RBW_f21270_22*RBW_pi1300m_2 + SF_Pm_D_TP_4*RBW_f21270_21*RBW_pi1300m_2);
*/   
   // D0 -> pi(1300)- {f0 pi-} pi+
   amplitude += fitpara[12]*fitpara[11]*(RBW_pi1300m_1*PiPiS_12_0 + RBW_pi1300m_1*PiPiS_11_0 + RBW_pi1300m_2*PiPiS_22_0 + RBW_pi1300m_2*PiPiS_21_0);
//   amplitude += fitpara[13]*fitpara[11]*(RBW_pi1300m_1*PiPiS_12_5 + RBW_pi1300m_1*PiPiS_11_5 + RBW_pi1300m_2*PiPiS_22_5 + RBW_pi1300m_2*PiPiS_21_5);
 
   amplitude += fitpara[13]*fitpara[11]*(RBW_pi1300m_1*PiPiS_12_6 + RBW_pi1300m_1*PiPiS_11_6 + RBW_pi1300m_2*PiPiS_22_6 + RBW_pi1300m_2*PiPiS_21_6);
 
   // D0 -> pi(1300)+ {rho(770)0 pi+} pi-
   double SF_Pp_P_VP_1 = contract_11_0(T1_Pip2Pim1,T1_Pip2Pim1Pip1);
   double SF_Pp_P_VP_2 = contract_11_0(T1_Pip1Pim1,T1_Pip1Pim1Pip2);
   double SF_Pp_P_VP_3 = contract_11_0(T1_Pip2Pim2,T1_Pip2Pim2Pip1);
   double SF_Pp_P_VP_4 = contract_11_0(T1_Pip1Pim2,T1_Pip1Pim2Pip2);
 
   amplitude += fitpara[14]*(SF_Pp_P_VP_1*GS_rho770_21*RBW_pi1300p_1 + SF_Pp_P_VP_2*GS_rho770_11*RBW_pi1300p_1 + SF_Pp_P_VP_3*GS_rho770_22*RBW_pi1300p_2 + SF_Pp_P_VP_4*GS_rho770_12*RBW_pi1300p_2);
/*
   // D0 -> pi(1300)+ {f2(1270)0 pi+} pi-
   double SF_Pp_D_TP_1 = contract_22_0(T2_Pip2Pim1,T2_Pip2Pim1Pip1);
   double SF_Pp_D_TP_2 = contract_22_0(T2_Pip1Pim1,T2_Pip1Pim1Pip2);
   double SF_Pp_D_TP_3 = contract_22_0(T2_Pip2Pim2,T2_Pip2Pim2Pip1);
   double SF_Pp_D_TP_4 = contract_22_0(T2_Pip1Pim2,T2_Pip1Pim2Pip2);
 
   amplitude += fitpara[12]*fitpara[15]*(SF_Pp_D_TP_1*RBW_f21270_21*RBW_pi1300p_1 + SF_Pp_D_TP_2*RBW_f21270_11*RBW_pi1300p_1 + SF_Pp_D_TP_3*RBW_f21270_22*RBW_pi1300p_2 + SF_Pp_D_TP_4*RBW_f21270_12*RBW_pi1300p_2);
*/   
   // D0 -> pi(1300)+ {f0 pi+} pi-
   amplitude += fitpara[12]*fitpara[14]*(RBW_pi1300p_1*PiPiS_21_0 + RBW_pi1300p_1*PiPiS_11_0 + RBW_pi1300p_2*PiPiS_22_0 + RBW_pi1300p_2*PiPiS_12_0);
//   amplitude += fitpara[13]*fitpara[15]*(RBW_pi1300p_1*PiPiS_21_5 + RBW_pi1300p_1*PiPiS_11_5 + RBW_pi1300p_2*PiPiS_22_5 + RBW_pi1300p_2*PiPiS_12_5);
 
   amplitude += fitpara[13]*fitpara[14]*(RBW_pi1300p_1*PiPiS_21_6 + RBW_pi1300p_1*PiPiS_11_6 + RBW_pi1300p_2*PiPiS_22_6 + RBW_pi1300p_2*PiPiS_12_6);
 
   // D0 -> rho rho [S]
   double SF_VV_S_1 = contract_11_0(T1_Pip1Pim1, T1_Pip2Pim2);
   double SF_VV_S_3 = contract_11_0(T1_Pip1Pim2, T1_Pip2Pim1);
 
   amplitude += fitpara[15]*(SF_VV_S_1*GS_rho770_11*GS_rho770_22+SF_VV_S_3*GS_rho770_12*GS_rho770_21);
 
   // D0 -> rho rho [P]
   double SF_VV_P_1 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, T1_Pip1Pim1),T1_Pip2Pim2),T1_2z11), D0);
   double SF_VV_P_2 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, T1_Pip2Pim2),T1_Pip1Pim1),T1_2z12), D0);
   double SF_VV_P_3 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, T1_Pip1Pim2),T1_Pip2Pim1),T1_2z21), D0);
   double SF_VV_P_4 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, T1_Pip2Pim1),T1_Pip1Pim2),T1_2z22), D0);
 
   amplitude += fitpara[16]*(SF_VV_P_1*GS_rho770_11*GS_rho770_22+SF_VV_P_3*GS_rho770_12*GS_rho770_21);
 
   // D0 -> rho rho [D]
   double SF_VV_D_1 = contract_11_0(contract_21_1(T2_2z11,T1_Pip2Pim2), T1_Pip1Pim1);
   double SF_VV_D_3 = contract_11_0(contract_21_1(T2_2z21,T1_Pip2Pim1), T1_Pip1Pim2);
 
   amplitude += fitpara[17]*(SF_VV_D_1*GS_rho770_11*GS_rho770_22+SF_VV_D_3*GS_rho770_12*GS_rho770_21);
 
   // D0 -> rho rho' [P]
   amplitude += fitpara[18]*(SF_VV_P_1*GS_rho770_11*GS_rho1450_22+SF_VV_P_2*GS_rho770_22*GS_rho1450_11 + SF_VV_P_3*GS_rho770_12*GS_rho1450_21 + SF_VV_P_3*GS_rho770_21*GS_rho1450_12);
 
   // D0 ->rho f0
   double SF_VS_P_1 = contract_11_0(T1_Pip1Pim1,T1_2z11);
   double SF_VS_P_2 = contract_11_0(T1_Pip2Pim2,T1_2z12);
   double SF_VS_P_3 = contract_11_0(T1_Pip1Pim2,T1_2z21);
   double SF_VS_P_4 = contract_11_0(T1_Pip2Pim1,T1_2z22);
 
   amplitude += fitpara[19]*(SF_VS_P_1*GS_rho770_11*PiPiS_22_0 + SF_VS_P_2*GS_rho770_22*PiPiS_11_0 + SF_VS_P_3*GS_rho770_12*PiPiS_21_0 + SF_VS_P_4*GS_rho770_21*PiPiS_12_0);
   amplitude += fitpara[20]*(SF_VS_P_1*GS_rho770_11*PiPiS_22_5 + SF_VS_P_2*GS_rho770_22*PiPiS_11_5 + SF_VS_P_3*GS_rho770_12*PiPiS_21_5 + SF_VS_P_4*GS_rho770_21*PiPiS_12_5);
   amplitude += fitpara[21]*(SF_VS_P_1*GS_rho770_11*PiPiS_22_6 + SF_VS_P_2*GS_rho770_22*PiPiS_11_6 + SF_VS_P_3*GS_rho770_12*PiPiS_21_6 + SF_VS_P_4*GS_rho770_21*PiPiS_12_6);
 
   // D0 -> f0f0
   //00
   amplitude += fitpara[22]*(PiPiS_11_0*PiPiS_22_0 + PiPiS_12_0*PiPiS_21_0 + PiPiS_22_0*PiPiS_11_0 + PiPiS_21_0*PiPiS_12_0);
   amplitude += fitpara[23]*(PiPiS_11_0*PiPiS_22_1 + PiPiS_12_0*PiPiS_21_1 + PiPiS_22_0*PiPiS_11_1 + PiPiS_21_0*PiPiS_12_1);
   amplitude += fitpara[24]*(PiPiS_11_1*PiPiS_22_1 + PiPiS_12_1*PiPiS_21_1 + PiPiS_22_1*PiPiS_11_1 + PiPiS_21_1*PiPiS_12_1);
   amplitude += fitpara[25]*(PiPiS_11_1*PiPiS_22_5 + PiPiS_12_1*PiPiS_21_5 + PiPiS_22_1*PiPiS_11_5 + PiPiS_21_1*PiPiS_12_5);
   amplitude += fitpara[26]*(PiPiS_11_5*PiPiS_22_5 + PiPiS_12_5*PiPiS_21_5 + PiPiS_22_5*PiPiS_11_5 + PiPiS_21_5*PiPiS_12_5);
   amplitude += fitpara[27]*(PiPiS_11_5*PiPiS_22_6 + PiPiS_12_5*PiPiS_21_6 + PiPiS_22_5*PiPiS_11_6 + PiPiS_21_5*PiPiS_12_6);
 
   // D0 -> f2 f0
   double SF_TS_D_1 = contract_22_0(T2_Pip1Pim1,T2_2z11);
   double SF_TS_D_2 = contract_22_0(T2_Pip2Pim2,T2_2z12);
   double SF_TS_D_3 = contract_22_0(T2_Pip1Pim2,T2_2z21);
   double SF_TS_D_4 = contract_22_0(T2_Pip2Pim1,T2_2z22);
         
   amplitude += fitpara[28]*(SF_TS_D_1*RBW_f21270_11*PiPiS_22_5 + SF_TS_D_2*RBW_f21270_22*PiPiS_11_5 + SF_TS_D_3*RBW_f21270_12*PiPiS_21_5 + SF_TS_D_4*RBW_f21270_21*PiPiS_12_5);
   amplitude += fitpara[29]*(SF_TS_D_1*RBW_f21270_11*PiPiS_22_6 + SF_TS_D_2*RBW_f21270_22*PiPiS_11_6 + SF_TS_D_3*RBW_f21270_12*PiPiS_21_6 + SF_TS_D_4*RBW_f21270_21*PiPiS_12_6);
 
   return amplitude;
 
}
 
int EvtD0To2pip2pim::CalAmp()
{
   m_AmpD0 = CalD0Amp();
   m_AmpDb = CalDbAmp();
   return 0;
}
 
double EvtD0To2pip2pim::mag2(complex<double> x)
{
   double temp = x.real()*x.real() + x.imag()*x.imag();
   return temp;
}
 
double EvtD0To2pip2pim::arg(complex<double> x)
{
   double temp = atan(x.imag()/x.real());
   if(x.real()<0) temp=temp+TMath::Pi();
   return temp;
}
double EvtD0To2pip2pim::Get_strongPhase()
{
   double temp = arg(m_AmpD0) - arg(m_AmpDb);
   while (temp < -TMath::Pi()){
      temp += 2.0*TMath::Pi();
   }
   while (temp > TMath::Pi()){
      temp -= 2.0*TMath::Pi();
   }  
   return temp;
}
 
double EvtD0To2pip2pim::AmplitudeSquare(int charm, int tagmode){
 
   EvtVector4R dp1=GetDaugMomLab(0),dp2=GetDaugMomLab(1),dp3=GetDaugMomLab(2), dp4=GetDaugMomLab(3); // pi+ pi+ pi- pi-
   EvtVector4R mp = dp1 + dp2 + dp3 + dp4;
 
   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);
   EvtVector4R dp4bst = boostTo(dp4, boostp3);
 
   double p4pip1[4], p4pip2[4], p4pim1[4], p4pim2[4];
   for(int i=0; i<3; i++){
      p4pip1[i]=dp1bst.get(i+1);
      p4pip2[i]=dp2bst.get(i+1);
      p4pim1[i]=dp3bst.get(i+1);
      p4pim2[i]=dp4bst.get(i+1);
   }
   p4pip1[3]=dp1bst.get(0);
   p4pip2[3]=dp2bst.get(0);
   p4pim1[3]=dp3bst.get(0);
   p4pim2[3]=dp4bst.get(0);
 
   setInput(p4pip1, p4pim1, p4pip2, p4pim2);
   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;
}