EvtD0Topippim2pi0.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:  EvtD0Topippim2pi0.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/EvtD0Topippim2pi0.hh"
#include <stdlib.h>
#include <iostream>
#include <string>
#include <complex>
#include <vector>
#include <math.h>
#include "TMath.h"
using namespace std;  //::endl;
 
EvtD0Topippim2pi0::~EvtD0Topippim2pi0() {}
 
void EvtD0Topippim2pi0::getName(std::string& model_name){
  model_name="D0Topippim2pi0";     
}
 
EvtDecayBase* EvtD0Topippim2pi0::clone(){
  return new EvtD0Topippim2pi0;
}
 
void EvtD0Topippim2pi0::init(){
   checkNArg(2);
   checkNDaug(4);
   charm   = getArg(0);
   tagmode = getArg(1);
   //std::cout<<"Initializing EvtD0Topippim2pi0: charm="<<charm<<" tagmode= "<<tagmode<<std::endl;
 
   double mag[35], pha[35];
   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]= 0.584175; pha[7]= -2.89693;
   mag[8]= 127.158;  pha[8]= -2.47773;
   mag[9]= 339.914;  pha[9]= 2.22856;
   mag[10]=0.320888; pha[10]=-2.6194;
   mag[11]=0.366283; pha[11]=-0.26867;
   mag[12]=14.1344;  pha[12]=-0.41164;
   mag[13]=86.0865;  pha[13]=-2.49649;
   mag[14]=6.1541;   pha[14]=-1.18299;
   mag[15]=56.6067;  pha[15]=0.142977;
   mag[16]=92.3073;  pha[16]=-2.15881;
   mag[17]=80.9453;  pha[17]=0.825815;
   mag[18]=16.9555;  pha[18]=-2.98994;
   mag[19]=9.72524;  pha[19]=-1.39929;
   mag[20]=5.71448;  pha[20]=0.271902;
   mag[21]=21.4195;  pha[21]=-1.23701;
   mag[22]=56.8867;  pha[22]=-0.385837;
   mag[23]=231.626;  pha[23]=2.14842;
   mag[24]=2938.45;  pha[24]=-0.693491;
   mag[25]=7252.7;   pha[25]=2.23659;
   mag[26]=5165.87;  pha[26]=0.913557;
   mag[27]=11508.6;  pha[27]=-1.07187;
   mag[28]=2461.86;  pha[28]=1.8709;
   mag[29]=8757.75;  pha[29]=2.40756;
   mag[30]=19.7413;  pha[30]=-1.0753;
   mag[31]=66.3826;  pha[31]=2.34666;
   mag[32]=11.2904;  pha[32]=-0.822345;
   mag[33]=2.04576;  pha[33]=-0.281429;
   mag[34]=0.57927;  pha[34]=2.7182;
 
   fitpara.clear();
   for(int i=0; i<35; 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;
   m_Pi0  = 0.134977;
   m2_Pi = m_Pi*m_Pi;
   m2_Pi0 = m_Pi0*m_Pi0;
 
   m0_rho7700 = 0.77526;
   w0_rho7700 = 0.1478;
 
   m0_rho770p = 0.77511;
   w0_rho770p = 0.1491;
 
   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_pi21670 = 1.6706;
   w0_pi21670 = 0.258;
 
   m0_pi11400 = 1.354;
   w0_pi11400 = 0.330;
 
   m0_h11170  = 1.166;
   w0_h11170  = 0.375;
 
   m0_omega   = 0.78265;
   w0_omega   = 0.00849;
 
   m0_phi     = 1.019461;
   w0_phi     = 0.004249;
 
   s0_prod = -5.0;
 
}
 
void EvtD0Topippim2pi0::initProbMax(){
  setProbMax(1989484027.5);
}
 
void EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::setInput(double* pip, double* pim, double* pi01, double* pi02){
   m_Pip.clear(); m_Pim.clear(); m_Pi01.clear(); m_Pi02.clear();
   m_Pip.push_back(pip[0]); m_Pim.push_back(pim[0]); m_Pi01.push_back(pi01[0]); m_Pi02.push_back(pi02[0]);
   m_Pip.push_back(pip[1]); m_Pim.push_back(pim[1]); m_Pi01.push_back(pi01[1]); m_Pi02.push_back(pi02[1]);
   m_Pip.push_back(pip[2]); m_Pim.push_back(pim[2]); m_Pi01.push_back(pi01[2]); m_Pi02.push_back(pi02[2]);
   m_Pip.push_back(pip[3]); m_Pim.push_back(pim[3]); m_Pi01.push_back(pi01[3]); m_Pi02.push_back(pi02[3]);
}
 
vector<double> EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::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> EvtD0Topippim2pi0::CalD0Amp(){
   return Amp(m_Pip, m_Pim, m_Pi01, m_Pi02); 
}
complex<double> EvtD0Topippim2pi0::CalDbAmp(){
 
   vector<double> cpPip;  cpPip.clear();
   vector<double> cpPim;  cpPim.clear();
   vector<double> cpPi01; cpPi01.clear();
   vector<double> cpPi02; cpPi02.clear();
 
   cpPip.push_back(-m_Pim[0]); cpPim.push_back(-m_Pip[0]); cpPi01.push_back(-m_Pi01[0]); cpPi02.push_back(-m_Pi02[0]);
   cpPip.push_back(-m_Pim[1]); cpPim.push_back(-m_Pip[1]); cpPi01.push_back(-m_Pi01[1]); cpPi02.push_back(-m_Pi02[1]);
   cpPip.push_back(-m_Pim[2]); cpPim.push_back(-m_Pip[2]); cpPi01.push_back(-m_Pi01[2]); cpPi02.push_back(-m_Pi02[2]);
   cpPip.push_back( m_Pim[3]); cpPim.push_back( m_Pip[3]); cpPi01.push_back( m_Pi01[3]); cpPi02.push_back( m_Pi02[3]);
 
   return Amp(cpPip, cpPim, cpPi01, cpPi02); 
}
 
complex<double> EvtD0Topippim2pi0::Amp(vector<double> Pip, vector<double> Pim, vector<double> Pi01, vector<double> Pi02)
{
 
   vector<double> PipPim; PipPim.clear();
   vector<double> PipPi01; PipPi01.clear();
   vector<double> PipPi02; PipPi02.clear();
   vector<double> PimPi01; PimPi01.clear();
   vector<double> PimPi02; PimPi02.clear();
   vector<double> Pi01Pi02; Pi01Pi02.clear();
 
   PipPim   = sum_tensor(Pip, Pim);
   PipPi01  = sum_tensor(Pip, Pi01);
   PipPi02  = sum_tensor(Pip, Pi02);
   PimPi01  = sum_tensor(Pim, Pi01);
   PimPi02  = sum_tensor(Pim, Pi02);
   Pi01Pi02 = sum_tensor(Pi01, Pi02);
 
   vector<double> PipPimPi01; PipPimPi01.clear();
   vector<double> PipPimPi02; PipPimPi02.clear();
   vector<double> PipPi01Pi02; PipPi01Pi02.clear();
   vector<double> PimPi01Pi02; PimPi01Pi02.clear();
 
   PipPimPi01 = sum_tensor(PipPim, Pi01);
   PipPimPi02 = sum_tensor(PipPim, Pi02);
   PipPi01Pi02 = sum_tensor(PipPi01, Pi02);
   PimPi01Pi02 = sum_tensor(PimPi01, Pi02);
 
   vector<double> D0; D0.clear();
   D0 = sum_tensor(PipPimPi01, Pi02);
 
   double M2_PipPim  = contract_11_0(PipPim, PipPim);
   double M2_PipPi01 = contract_11_0(PipPi01, PipPi01);
   double M2_PipPi02  = contract_11_0(PipPi02, PipPi02);
   double M2_PimPi01  = contract_11_0(PimPi01, PimPi01);
   double M2_PimPi02  = contract_11_0(PimPi02, PimPi02);
   double M2_Pi01Pi02  = contract_11_0(Pi01Pi02, Pi01Pi02);
 
   double M2_PipPimPi01 = contract_11_0(PipPimPi01, PipPimPi01);
   double M2_PipPimPi02 = contract_11_0(PipPimPi02, PipPimPi02);
   double M2_PipPi01Pi02 = contract_11_0(PipPi01Pi02, PipPi01Pi02);
   double M2_PimPi01Pi02 = contract_11_0(PimPi01Pi02, PimPi01Pi02);
   double M2_D0 = contract_11_0(D0, D0);
 
   complex<double> GS_rho770_pm = GS(M2_PipPim,  m0_rho7700, w0_rho7700, m2_Pi, m2_Pi, rRes, 1);
   complex<double> GS_rho770_p1 = GS(M2_PipPi01, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
   complex<double> GS_rho770_p2 = GS(M2_PipPi02, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
   complex<double> GS_rho770_m1 = GS(M2_PimPi01, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
   complex<double> GS_rho770_m2 = GS(M2_PimPi02, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
 
   complex<double> RBW_f21270_pm = RBW(M2_PipPim, m0_f21270, w0_f21270, m2_Pi, m2_Pi, rRes, 2);
   complex<double> RBW_f21270_00 = RBW(M2_Pi01Pi02, m0_f21270, w0_f21270, m2_Pi0, m2_Pi0, rRes, 2);
 
   complex<double> PiPiS_pm_0 = Fvector(M2_PipPim, s0_prod, 0);
   complex<double> PiPiS_00_0 = Fvector(M2_Pi01Pi02, s0_prod, 0);
 
   complex<double> PiPiS_pm_1 = Fvector(M2_PipPim, s0_prod, 1);
   complex<double> PiPiS_00_1 = Fvector(M2_Pi01Pi02, s0_prod, 1);
 
   complex<double> PiPiS_pm_5 = Fvector(M2_PipPim, s0_prod, 5);
   complex<double> PiPiS_00_5 = Fvector(M2_Pi01Pi02, s0_prod, 5);
 
   complex<double> PiPiS_pm_6 = Fvector(M2_PipPim, s0_prod, 6);
   complex<double> PiPiS_00_6 = Fvector(M2_Pi01Pi02, s0_prod, 6);
 
   complex<double> RBW_a11260_p = RBWa1260(M2_PipPi01Pi02, m0_a11260, g1_a11260, g2_a11260);
   complex<double> RBW_a11260_m = RBWa1260(M2_PimPi01Pi02, m0_a11260, g1_a11260, g2_a11260);
   complex<double> RBW_a11260_01 = RBWa1260(M2_PipPimPi01, m0_a11260, g1_a11260, g2_a11260);
   complex<double> RBW_a11260_02 = RBWa1260(M2_PipPimPi02, m0_a11260, g1_a11260, g2_a11260);
 
   complex<double> RBW_a11420_p = RBW(M2_PipPi01Pi02, m0_a11420, w0_a11420,-1,-1,-1,-1);
   complex<double> RBW_a11420_m = RBW(M2_PimPi01Pi02, m0_a11420, w0_a11420,-1,-1,-1,-1);
   complex<double> RBW_a11420_01 = RBW(M2_PipPimPi01, m0_a11420, w0_a11420,-1,-1,-1,-1);
   complex<double> RBW_a11420_02 = RBW(M2_PipPimPi02, m0_a11420, w0_a11420,-1,-1,-1,-1);
 
   complex<double> RBW_omega_01 = RBW(M2_PipPimPi01, m0_omega, w0_omega,-1,-1,-1,-1);
   complex<double> RBW_omega_02 = RBW(M2_PipPimPi02, m0_omega, w0_omega,-1,-1,-1,-1);
 
   complex<double> RBW_phi_01 = RBW(M2_PipPimPi01, m0_phi, w0_phi,-1,-1,-1,-1);
   complex<double> RBW_phi_02 = RBW(M2_PipPimPi02, m0_phi, w0_phi,-1,-1,-1,-1);
 
   complex<double> RBW_a21320_p = RBW(M2_PipPi01Pi02, m0_a21320, w0_a21320,-1,-1,-1,-1);
   complex<double> RBW_a21320_m = RBW(M2_PimPi01Pi02, m0_a21320, w0_a21320,-1,-1,-1,-1);
 
   complex<double> RBW_pi1300_p = RBWpi1300(M2_PipPi01Pi02, m0_pi1300, w0_pi1300);
   complex<double> RBW_pi1300_m = RBWpi1300(M2_PimPi01Pi02, m0_pi1300, w0_pi1300);
   complex<double> RBW_pi1300_01 = RBWpi1300(M2_PipPimPi01, m0_pi1300, w0_pi1300);
   complex<double> RBW_pi1300_02 = RBWpi1300(M2_PipPimPi02, m0_pi1300, w0_pi1300);
 
   complex<double> RBW_h11170_01 = RBW(M2_PipPimPi01, m0_h11170, w0_h11170,-1,-1,-1,-1);
   complex<double> RBW_h11170_02 = RBW(M2_PipPimPi02, m0_h11170, w0_h11170,-1,-1,-1,-1);
 
   complex<double> RBW_pi21670_01 = RBW(M2_PipPimPi01, m0_pi21670, w0_pi21670,-1,-1,-1,-1);
   complex<double> RBW_pi21670_02 = RBW(M2_PipPimPi02, m0_pi21670, w0_pi21670,-1,-1,-1,-1);
 
   // D->XX Projection
   vector<double> Proj1_3p; Proj1_3p.clear();
   vector<double> Proj1_3m; Proj1_3m.clear();
   vector<double> Proj1_3z1; Proj1_3z1.clear();
   vector<double> Proj1_3z2; Proj1_3z2.clear();
 
   Proj1_3p  = ProjectionTensors(PipPi01Pi02,1);
   Proj1_3m  = ProjectionTensors(PimPi01Pi02,1);
   Proj1_3z1 = ProjectionTensors(PipPimPi01,1);
   Proj1_3z2 = ProjectionTensors(PipPimPi02,1);
 
   vector<double> Proj2_3p; Proj2_3p.clear();
   vector<double> Proj2_3m; Proj2_3m.clear();
   vector<double> Proj2_3z1; Proj2_3z1.clear();
   vector<double> Proj2_3z2; Proj2_3z2.clear();
 
   Proj2_3p  = ProjectionTensors(PipPi01Pi02,2);
   Proj2_3m  = ProjectionTensors(PimPi01Pi02,2);
   Proj2_3z1 = ProjectionTensors(PipPimPi01,2);
   Proj2_3z2 = ProjectionTensors(PipPimPi02,2);
 
   // X->PP Orbital
   vector<double> T1_PipPim; T1_PipPim.clear();
   vector<double> T1_PipPi01; T1_PipPi01.clear();
   vector<double> T1_PipPi02; T1_PipPi02.clear();
   vector<double> T1_PimPi01; T1_PimPi01.clear();
   vector<double> T1_PimPi02; T1_PimPi02.clear();
   vector<double> T1_Pi01Pi02; T1_Pi01Pi02.clear();
 
   T1_PipPim = OrbitalTensors(PipPim, Pip, Pim, rRes, 1);
   T1_PipPi01 = OrbitalTensors(PipPi01, Pip, Pi01, rRes, 1);
   T1_PipPi02 = OrbitalTensors(PipPi02, Pip, Pi02, rRes, 1);
   T1_PimPi01 = OrbitalTensors(PimPi01, Pim, Pi01, rRes, 1);
   T1_PimPi02 = OrbitalTensors(PimPi02, Pim, Pi02, rRes, 1);
   T1_Pi01Pi02 = OrbitalTensors(Pi01Pi02, Pi01, Pi02, rRes, 1);
 
   vector<double> T2_PipPim; T2_PipPim.clear();
   vector<double> T2_Pi01Pi02; T2_Pi01Pi02.clear();
 
   T2_PipPim = OrbitalTensors(PipPim, Pip, Pim, rRes, 2);
   T2_Pi01Pi02 = OrbitalTensors(Pi01Pi02, Pi01, Pi02, rRes, 2);
 
   // X->YP Orbital
   vector<double> T1_PipPimPi01; T1_PipPimPi01.clear();
   vector<double> T1_PipPimPi02; T1_PipPimPi02.clear();
   vector<double> T1_PipPi01Pi02; T1_PipPi01Pi02.clear();
   vector<double> T1_PipPi02Pi01; T1_PipPi02Pi01.clear();
   vector<double> T1_PimPi01Pi02; T1_PimPi01Pi02.clear();
   vector<double> T1_PimPi02Pi01; T1_PimPi02Pi01.clear();
   vector<double> T1_PipPi01Pim; T1_PipPi01Pim.clear();
   vector<double> T1_PipPi02Pim; T1_PipPi02Pim.clear();
   vector<double> T1_PimPi01Pip; T1_PimPi01Pip.clear();
   vector<double> T1_PimPi02Pip; T1_PimPi02Pip.clear();
   vector<double> T1_Pi01Pi02Pip; T1_Pi01Pi02Pip.clear();
   vector<double> T1_Pi01Pi02Pim; T1_Pi01Pi02Pim.clear();
 
   T1_PipPimPi01 = OrbitalTensors(PipPimPi01, PipPim, Pi01, rRes, 1);
   T1_PipPimPi02 = OrbitalTensors(PipPimPi02, PipPim, Pi02, rRes, 1);
   T1_PipPi01Pi02 = OrbitalTensors(PipPi01Pi02, PipPi01, Pi02, rRes, 1);
   T1_PipPi02Pi01 = OrbitalTensors(PipPi01Pi02, PipPi02, Pi01, rRes, 1);
   T1_PimPi01Pi02 = OrbitalTensors(PimPi01Pi02, PimPi01, Pi02, rRes, 1);
   T1_PimPi02Pi01 = OrbitalTensors(PimPi01Pi02, PimPi02, Pi01, rRes, 1);
   T1_PipPi01Pim = OrbitalTensors(PipPimPi01, PipPi01, Pim, rRes, 1);
   T1_PipPi02Pim = OrbitalTensors(PipPimPi02, PipPi02, Pim, rRes, 1);
   T1_PimPi01Pip = OrbitalTensors(PipPimPi01, PimPi01, Pip, rRes, 1);
   T1_PimPi02Pip = OrbitalTensors(PipPimPi02, PimPi02, Pip, rRes, 1);
   T1_Pi01Pi02Pip = OrbitalTensors(PipPi01Pi02, Pi01Pi02, Pip, rRes, 1);
   T1_Pi01Pi02Pim = OrbitalTensors(PimPi01Pi02, Pi01Pi02, Pim, rRes, 1);
 
   vector<double> T2_PipPimPi01; T2_PipPimPi01.clear();
   vector<double> T2_PipPimPi02; T2_PipPimPi02.clear();
   vector<double> T2_PipPi01Pi02; T2_PipPi01Pi02.clear();
   vector<double> T2_PipPi02Pi01; T2_PipPi02Pi01.clear();
   vector<double> T2_PimPi01Pi02; T2_PimPi01Pi02.clear();
   vector<double> T2_PimPi02Pi01; T2_PimPi02Pi01.clear();
   vector<double> T2_PipPi01Pim; T2_PipPi01Pim.clear();
   vector<double> T2_PipPi02Pim; T2_PipPi02Pim.clear();
   vector<double> T2_PimPi01Pip; T2_PimPi01Pip.clear();
   vector<double> T2_PimPi02Pip; T2_PimPi02Pip.clear();
   vector<double> T2_Pi01Pi02Pip; T2_Pi01Pi02Pip.clear();
   vector<double> T2_Pi01Pi02Pim; T2_Pi01Pi02Pim.clear();
 
   T2_PipPimPi01 = OrbitalTensors(PipPimPi01, PipPim, Pi01, rRes, 2);
   T2_PipPimPi02 = OrbitalTensors(PipPimPi02, PipPim, Pi02, rRes, 2);
   T2_PipPi01Pi02 = OrbitalTensors(PipPi01Pi02, PipPi01, Pi02, rRes, 2);
   T2_PipPi02Pi01 = OrbitalTensors(PipPi01Pi02, PipPi02, Pi01, rRes, 2);
   T2_PimPi01Pi02 = OrbitalTensors(PimPi01Pi02, PimPi01, Pi02, rRes, 2);
   T2_PimPi02Pi01 = OrbitalTensors(PimPi01Pi02, PimPi02, Pi01, rRes, 2);
   T2_PipPi01Pim = OrbitalTensors(PipPimPi01, PipPi01, Pim, rRes, 2);
   T2_PipPi02Pim = OrbitalTensors(PipPimPi02, PipPi02, Pim, rRes, 2);
   T2_PimPi01Pip = OrbitalTensors(PipPimPi01, PimPi01, Pip, rRes, 2);
   T2_PimPi02Pip = OrbitalTensors(PipPimPi02, PimPi02, Pip, rRes, 2);
   T2_Pi01Pi02Pip = OrbitalTensors(PipPi01Pi02, Pi01Pi02, Pip, rRes, 2);
   T2_Pi01Pi02Pim = OrbitalTensors(PimPi01Pi02, Pi01Pi02, Pim, rRes, 2);
 
   // D->XX Orbital
   vector<double> T1_2pm12; T1_2pm12.clear();
   vector<double> T1_2p1m2; T1_2p1m2.clear();
   vector<double> T1_2p2m1; T1_2p2m1.clear();
 
   T1_2pm12 = OrbitalTensors(D0, PipPim, Pi01Pi02, rD, 1);
   T1_2p1m2 = OrbitalTensors(D0, PipPi01, PimPi02, rD, 1);
   T1_2p2m1 = OrbitalTensors(D0, PipPi02, PimPi01, rD, 1);
 
   vector<double> T2_2pm12; T2_2pm12.clear();
   vector<double> T2_2p1m2; T2_2p1m2.clear();
   vector<double> T2_2p2m1; T2_2p2m1.clear();
 
   T2_2pm12 = OrbitalTensors(D0, PipPim, Pi01Pi02, rD, 2);
   T2_2p1m2 = OrbitalTensors(D0, PipPi01, PimPi02, rD, 2);
   T2_2p2m1 = OrbitalTensors(D0, PipPi02, PimPi01, rD, 2);
 
   // D->XP Orbital
   vector<double> T1_3pm; T1_3pm.clear();
   vector<double> T1_3mp; T1_3mp.clear();
   vector<double> T1_3z12; T1_3z12.clear();
   vector<double> T1_3z21; T1_3z21.clear();
 
   T1_3pm  = OrbitalTensors(D0, PipPi01Pi02, Pim, rD, 1);
   T1_3mp  = OrbitalTensors(D0, PimPi01Pi02, Pip, rD, 1);
   T1_3z12 = OrbitalTensors(D0, PipPimPi01, Pi02, rD, 1);
   T1_3z21 = OrbitalTensors(D0, PipPimPi02, Pi01, rD, 1);
 
   vector<double> T2_3pm;  T2_3pm.clear();
   vector<double> T2_3mp;  T2_3mp.clear();
   vector<double> T2_3z12; T2_3z12.clear();
   vector<double> T2_3z21; T2_3z21.clear();
 
   T2_3pm  = OrbitalTensors(D0, PipPi01Pi02, Pim, rD, 2);
   T2_3mp  = OrbitalTensors(D0, PimPi01Pi02, Pip, rD, 2);
   T2_3z12 = OrbitalTensors(D0, PipPimPi01, Pi02, rD, 2);
   T2_3z21 = OrbitalTensors(D0, PipPimPi02, Pi01, rD, 2);
 
   complex<double> amplitude(0,0);
 
   // D0 -> a1(1260)+ {rho(770)+ pi0[S]} pi-
   double SF_Ap_S_Vp1P = contract_11_0(contract_21_1(Proj1_3p, T1_PipPi01), T1_3pm);
   double SF_Ap_S_Vp2P = contract_11_0(contract_21_1(Proj1_3p, T1_PipPi02), T1_3pm);
 
   amplitude += fitpara[0]*(SF_Ap_S_Vp1P*RBW_a11260_p*GS_rho770_p1 + SF_Ap_S_Vp2P*RBW_a11260_p*GS_rho770_p2);
 
   // D0 -> a1(1260)+ {rho(770)+ pi0[D]} pi-
   double SF_Ap_D_Vp1P = contract_11_0(contract_21_1(T2_PipPi01Pi02, T1_PipPi01), T1_3pm);
   double SF_Ap_D_Vp2P = contract_11_0(contract_21_1(T2_PipPi02Pi01, T1_PipPi02), T1_3pm);
 
//   cout<<SF_Ap_D_VP_1<<","<<SF_Ap_D_VP_2<<","<<SF_Ap_D_VP_3<<","<<SF_Ap_D_VP_4<<","<<endl;
//   cout<<"-------"<<endl;
   amplitude += fitpara[1]*(SF_Ap_D_Vp1P*RBW_a11260_p*GS_rho770_p1 + SF_Ap_D_Vp2P*RBW_a11260_p*GS_rho770_p2);
 
   // D0 -> a1(1260)+ {f2(1270) pi+ [P]} pi-
   double SF_Ap_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3p, T2_Pi01Pi02), T1_Pi01Pi02Pip), T1_3pm);
 
   amplitude += fitpara[2]*(SF_Ap_P_TP*RBW_a11260_p*RBW_f21270_00);
 
   // D0 -> a1(1260)+ {f0 pi+ [P]} pi-
   double SF_Ap_P_SP = contract_11_0(T1_3pm, T1_Pi01Pi02Pip);
 
   amplitude += fitpara[3]*(SF_Ap_P_SP*RBW_a11260_p*PiPiS_00_0);
   amplitude += fitpara[4]*(SF_Ap_P_SP*RBW_a11260_p*PiPiS_00_1);
   amplitude += fitpara[5]*(SF_Ap_P_SP*RBW_a11260_p*PiPiS_00_5);
 
   // D0 -> a1(1260)- { rho(770)- pi0 [S]} pi+
   double SF_Am_S_Vm1P = contract_11_0(contract_21_1(Proj1_3m, T1_PimPi01), T1_3mp);
   double SF_Am_S_Vm2P = contract_11_0(contract_21_1(Proj1_3m, T1_PimPi02), T1_3mp);
 
   amplitude += fitpara[6]*fitpara[0]*(SF_Am_S_Vm1P*RBW_a11260_m*GS_rho770_m1 + SF_Am_S_Vm2P*RBW_a11260_m*GS_rho770_m2);
 
   // D0 -> a1(1260)- {rho(770)- pi0[D]} pi+
   double SF_Am_D_Vm1P = contract_11_0(contract_21_1(T2_PimPi01Pi02, T1_PimPi01), T1_3mp);
   double SF_Am_D_Vm2P = contract_11_0(contract_21_1(T2_PimPi02Pi01, T1_PimPi02), T1_3mp);
 
   amplitude += fitpara[6]*fitpara[1]*(SF_Am_D_Vm1P*RBW_a11260_m*GS_rho770_m1 + SF_Am_D_Vm2P*RBW_a11260_m*GS_rho770_m2);
 
   // D0 -> a1(1260)- {f2(1270) pi- [P]} pi+
   double SF_Am_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3m, T2_Pi01Pi02), T1_Pi01Pi02Pim), T1_3mp);
 
   amplitude += fitpara[6]*fitpara[2]*(SF_Am_P_TP*RBW_a11260_m*RBW_f21270_00);
 
   // D0 -> a1(1260)- {f0 pi- [P]} pi+
   double SF_Am_P_SP = contract_11_0(T1_3mp, T1_Pi01Pi02Pim);
 
   amplitude += fitpara[6]*fitpara[3]*(SF_Am_P_SP*RBW_a11260_m*PiPiS_00_0);
   amplitude += fitpara[6]*fitpara[4]*(SF_Am_P_SP*RBW_a11260_m*PiPiS_00_1);
   amplitude += fitpara[6]*fitpara[5]*(SF_Am_P_SP*RBW_a11260_m*PiPiS_00_5);
   
   // D -> a1(1260)0 pi0
   double SF_A01_S_Vp1P = contract_11_0(contract_21_1(Proj1_3z1, T1_PipPi01), T1_3z12);
   double SF_A02_S_Vp2P = contract_11_0(contract_21_1(Proj1_3z2, T1_PipPi02), T1_3z21);
   double SF_A01_S_Vm1P = contract_11_0(contract_21_1(Proj1_3z1, T1_PimPi01), T1_3z12);
   double SF_A02_S_Vm2P = contract_11_0(contract_21_1(Proj1_3z2, T1_PimPi02), T1_3z21);
   double SF_A01_S_VzP  = contract_11_0(contract_21_1(Proj1_3z1, T1_PipPim), T1_3z12);
   double SF_A02_S_VzP  = contract_11_0(contract_21_1(Proj1_3z2, T1_PipPim), T1_3z21);
 
   amplitude += fitpara[7]*fitpara[0]*(SF_A01_S_Vp1P*RBW_a11260_01*GS_rho770_p1 + SF_A02_S_Vp2P*RBW_a11260_02*GS_rho770_p2 + SF_A01_S_Vm1P*RBW_a11260_01*GS_rho770_m1 + SF_A02_S_Vm2P*RBW_a11260_02*GS_rho770_m2);
 
   double SF_A01_D_Vp1P = contract_11_0(contract_21_1(T2_PipPi01Pim, T1_PipPi01), T1_3z12);
   double SF_A02_D_Vp2P = contract_11_0(contract_21_1(T2_PipPi02Pim, T1_PipPi02), T1_3z21);
   double SF_A01_D_Vm1P = contract_11_0(contract_21_1(T2_PimPi01Pip, T1_PimPi01), T1_3z12);
   double SF_A02_D_Vm2P = contract_11_0(contract_21_1(T2_PimPi02Pip, T1_PimPi02), T1_3z21);
 
   amplitude += fitpara[7]*fitpara[1]*(SF_A01_D_Vp1P*RBW_a11260_01*GS_rho770_p1 + SF_A02_D_Vp2P*RBW_a11260_02*GS_rho770_p2 + SF_A01_D_Vm1P*RBW_a11260_01*GS_rho770_m1 + SF_A02_D_Vm2P*RBW_a11260_02*GS_rho770_m2);
 
   double SF_A01_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3z1,T2_PipPim), T1_PipPimPi01), T1_3z12);
   double SF_A02_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3z2,T2_PipPim), T1_PipPimPi02), T1_3z21);
 
   amplitude += fitpara[7]*fitpara[2]*(-1.0)*(SF_A01_P_TP*RBW_a11260_01*RBW_f21270_pm + SF_A02_P_TP*RBW_a11260_02*RBW_f21270_pm);
 
   double SF_A01_P_SP = contract_11_0(T1_3z12, T1_PipPimPi01);
   double SF_A02_P_SP = contract_11_0(T1_3z21, T1_PipPimPi02);
 
   amplitude += fitpara[7]*fitpara[3]*(-1.0)*(SF_A01_P_SP*RBW_a11260_01*PiPiS_pm_0 + SF_A02_P_SP*RBW_a11260_02*PiPiS_pm_0);
   amplitude += fitpara[7]*fitpara[4]*(-1.0)*(SF_A01_P_SP*RBW_a11260_01*PiPiS_pm_1 + SF_A02_P_SP*RBW_a11260_02*PiPiS_pm_1);
   amplitude += fitpara[7]*fitpara[5]*(-1.0)*(SF_A01_P_SP*RBW_a11260_01*PiPiS_pm_5 + SF_A02_P_SP*RBW_a11260_02*PiPiS_pm_5);
 
   // D0 -> a1(1420)+ {rho(770)+ pi0[S]} pi-
//   amplitude += fitpara[8]*(SF_Ap_S_Vp1P*RBW_a11420_p*GS_rho770_p1 + SF_Ap_S_Vp2P*RBW_a11420_p*GS_rho770_p2);
//   amplitude += fitpara[9]*(SF_A01_S_Vp1P*RBW_a11420_01*GS_rho770_p1 + SF_A02_S_Vp2P*RBW_a11420_02*GS_rho770_p2 + SF_A01_S_Vm1P*RBW_a11420_01*GS_rho770_m1 + SF_A02_S_Vm2P*RBW_a11420_02*GS_rho770_m2);
   
   // D0 -> a1(1420)+ {f0 pi0+[S]} pi-
   amplitude += fitpara[8]*(SF_Ap_P_SP*RBW_a11420_p*PiPiS_00_5);
   amplitude += fitpara[9]*(SF_Ap_P_SP*RBW_a11420_p*PiPiS_00_6);
 
   // D0 -> a2(1320)+ {rho+ pi0} pi-
   double SF_Tp_D_Vp1P = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3p, T1_PipPi01)), PipPi01Pi02), contract_42_2(Proj2_3p, T2_3pm)), T2_PipPi01Pi02);
   double SF_Tp_D_Vp2P = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3p, T1_PipPi02)), PipPi01Pi02), contract_42_2(Proj2_3p, T2_3pm)), T2_PipPi02Pi01);
 
   amplitude += fitpara[10]*(SF_Tp_D_Vp1P*GS_rho770_p1*RBW_a21320_p + SF_Tp_D_Vp2P*GS_rho770_p2*RBW_a21320_p);
 
   // D0 -> a2(1320)- {rho- pi0} pi+
   double SF_Tm_D_Vm1P = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3m, T1_PimPi01)), PimPi01Pi02), contract_42_2(Proj2_3m, T2_3mp)), T2_PimPi01Pi02);
   double SF_Tm_D_Vm2P = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3m, T1_PimPi02)), PimPi01Pi02), contract_42_2(Proj2_3m, T2_3mp)), T2_PimPi02Pi01);
   amplitude += fitpara[11]*(SF_Tm_D_Vm1P*GS_rho770_m1*RBW_a21320_m + SF_Tm_D_Vm2P*GS_rho770_m2*RBW_a21320_m);
 
   // D0 -> h1(1170)0 {rho  pi0} pi0
   amplitude += fitpara[12]*(SF_A01_S_Vp1P*RBW_h11170_01*GS_rho770_p1 + SF_A02_S_Vp2P*RBW_h11170_02*GS_rho770_p2 - SF_A01_S_Vm1P*RBW_h11170_01*GS_rho770_m1 - SF_A02_S_Vm2P*RBW_h11170_02*GS_rho770_m2 - SF_A01_S_VzP*RBW_h11170_01*GS_rho770_pm - SF_A02_S_VzP*RBW_h11170_02*GS_rho770_pm);
 
   // D0 -> pi(1300)- {rho(770)- pi0} pi+
   double SF_Pm_P_Vm1P = contract_11_0(T1_PimPi01,T1_PimPi01Pi02);
   double SF_Pm_P_Vm2P = contract_11_0(T1_PimPi02,T1_PimPi02Pi01);
 
   amplitude += fitpara[13]*(SF_Pm_P_Vm1P*GS_rho770_m1*RBW_pi1300_m + SF_Pm_P_Vm2P*GS_rho770_m2*RBW_pi1300_m);
 
   // D0 -> pi(1300)- {f2 pi-} pi+
//   double SF_Pm_D_TP = contract_22_0(T2_Pi01Pi02,T2_Pi01Pi02Pim);
//   amplitude += fitpara[14]*fitpara[13]*(SF_Pm_D_TP*RBW_f21270_00*RBW_pi1300_m);
 
   // D0 -> pi(1300)- {f0 pi-} pi+
   amplitude += fitpara[14]*fitpara[13]*(RBW_pi1300_m*PiPiS_00_0);
//   amplitude += fitpara[15]*fitpara[13]*(RBW_pi1300_m*PiPiS_00_5);
   amplitude += fitpara[15]*fitpara[13]*(RBW_pi1300_m*PiPiS_00_6);
 
   // D0 -> pi(1300)+ {rho(770)+ pi0} pi-
   double SF_Pp_P_Vp1P = contract_11_0(T1_PipPi01,T1_PipPi01Pi02);
   double SF_Pp_P_Vp2P = contract_11_0(T1_PipPi02,T1_PipPi02Pi01);
 
   amplitude += fitpara[16]*(SF_Pp_P_Vp1P*GS_rho770_p1*RBW_pi1300_p + SF_Pp_P_Vp2P*GS_rho770_p2*RBW_pi1300_p);
 
   // D0 -> pi(1300)+ {f2 pi-} pi+
//   double SF_Pp_D_TP = contract_22_0(T2_Pi01Pi02,T2_Pi01Pi02Pip);
//   amplitude += fitpara[14]*fitpara[17]*(SF_Pp_D_TP*RBW_f21270_00*RBW_pi1300_p);
 
   // D0 -> pi(1300)+ {f0 pi+} pi-
   amplitude += fitpara[14]*fitpara[16]*(RBW_pi1300_p*PiPiS_00_0);
//   amplitude += fitpara[15]*fitpara[17]*(RBW_pi1300_p*PiPiS_00_5);
   amplitude += fitpara[15]*fitpara[16]*(RBW_pi1300_p*PiPiS_00_6);
 
   // D0 -> pi(1300)0 {rho pi} pi0
   double SF_P01_P_Vp1P = contract_11_0(T1_PipPi01,T1_PipPi01Pim);
   double SF_P02_P_Vp2P = contract_11_0(T1_PipPi02,T1_PipPi02Pim);
   double SF_P01_P_Vm1P = contract_11_0(T1_PimPi01,T1_PimPi01Pip);
   double SF_P02_P_Vm2P = contract_11_0(T1_PimPi02,T1_PimPi02Pip);
 
   amplitude += fitpara[17]*(SF_P01_P_Vp1P*RBW_pi1300_01*GS_rho770_p1 + SF_P02_P_Vp2P*RBW_pi1300_02*GS_rho770_p2 + SF_P01_P_Vm1P*RBW_pi1300_01*GS_rho770_m1 + SF_P02_P_Vm2P*RBW_pi1300_02*GS_rho770_m2);
 
 
//   double SF_P01_D_TP = contract_22_0(T2_PipPim,T2_PipPimPi01);
//   double SF_P02_D_TP = contract_22_0(T2_PipPim,T2_PipPimPi02);
//   amplitude += fitpara[14]*fitpara[18]*(-1.0)*(SF_P01_D_TP*RBW_f21270_pm*RBW_pi1300_01 + SF_P02_D_TP*RBW_f21270_pm*RBW_pi1300_02);
 
   amplitude += fitpara[14]*fitpara[17]*(-1.0)*(RBW_pi1300_01*PiPiS_pm_0 + RBW_pi1300_02*PiPiS_pm_0);
//   amplitude += fitpara[15]*fitpara[18]*(-1.0)*(RBW_pi1300_01*PiPiS_pm_5 + RBW_pi1300_02*PiPiS_pm_5);
   amplitude += fitpara[15]*fitpara[17]*(-1.0)*(RBW_pi1300_01*PiPiS_pm_6 + RBW_pi1300_02*PiPiS_pm_6);
   
   // D0 -> rho+ rho- [S]
   double SF_Vp1Vm2_S = contract_11_0(T1_PipPi01, T1_PimPi02);
   double SF_Vp2Vm1_S = contract_11_0(T1_PipPi02, T1_PimPi01);
 
   amplitude += fitpara[18]*(SF_Vp1Vm2_S*GS_rho770_p1*GS_rho770_m2 + SF_Vp2Vm1_S*GS_rho770_p2*GS_rho770_m1);
 
   // D0 -> rho+ rho- [P]
   double SF_Vp1Vm2_P = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, T1_PipPi01),T1_PimPi02),T1_2p1m2), D0);
   double SF_Vp2Vm1_P = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, T1_PipPi02),T1_PimPi01),T1_2p2m1), D0);
 
   amplitude += fitpara[19]*(SF_Vp1Vm2_P*GS_rho770_p1*GS_rho770_m2 + SF_Vp2Vm1_P*GS_rho770_p2*GS_rho770_m1);
 
   // D0 -> rho+ rho- [D]
   double SF_Vp1Vm2_D = contract_11_0(contract_21_1(T2_2p1m2,T1_PipPi01), T1_PimPi02);
   double SF_Vp2Vm1_D = contract_11_0(contract_21_1(T2_2p2m1,T1_PipPi02), T1_PimPi01);
   amplitude += fitpara[20]*(SF_Vp1Vm2_D*GS_rho770_p1*GS_rho770_m2 + SF_Vp2Vm1_D*GS_rho770_p2*GS_rho770_m1);
 
   // D0 -> rho0 (Pi0 Pi0)S
   double SF_VpmS12_P = contract_11_0(T1_PipPim,T1_2pm12);
 
   amplitude += fitpara[21]*(SF_VpmS12_P*GS_rho770_pm*PiPiS_00_0);
   amplitude += fitpara[22]*(SF_VpmS12_P*GS_rho770_pm*PiPiS_00_5);
   amplitude += fitpara[23]*(SF_VpmS12_P*GS_rho770_pm*PiPiS_00_6);
 
   // D0 -> f0f0
   //00
   amplitude += fitpara[24]*(PiPiS_pm_0*PiPiS_00_0 + PiPiS_00_0*PiPiS_pm_0);
   amplitude += fitpara[25]*(PiPiS_pm_0*PiPiS_00_1 + PiPiS_00_0*PiPiS_pm_1);
   amplitude += fitpara[26]*(PiPiS_pm_1*PiPiS_00_1 + PiPiS_00_1*PiPiS_pm_1);
   amplitude += fitpara[27]*(PiPiS_pm_1*PiPiS_00_5 + PiPiS_00_1*PiPiS_pm_5);
   amplitude += fitpara[28]*(PiPiS_pm_5*PiPiS_00_5 + PiPiS_00_5*PiPiS_pm_5);
   amplitude += fitpara[29]*(PiPiS_pm_5*PiPiS_00_6 + PiPiS_00_5*PiPiS_pm_6);
 
   // D0 -> f2(1270) f0
   double SF_TpmS00_D = contract_22_0(T2_PipPim, T2_2pm12);
   double SF_T00Spm_D = contract_22_0(T2_Pi01Pi02, T2_2pm12);
 
   amplitude += fitpara[30]*(SF_TpmS00_D*RBW_f21270_pm*PiPiS_00_5 + SF_T00Spm_D*RBW_f21270_00*PiPiS_pm_5);
   amplitude += fitpara[31]*(SF_TpmS00_D*RBW_f21270_pm*PiPiS_00_6 + SF_T00Spm_D*RBW_f21270_00*PiPiS_pm_6);
 
   // D0 -> pi2(1670)0[f2(1270) pi0] pi0
   double SF_PT01_S_TP = contract_22_0(contract_42_2(Proj2_3z1, T2_PipPim), T2_3z12);
   double SF_PT02_S_TP = contract_22_0(contract_42_2(Proj2_3z2, T2_PipPim), T2_3z21);
 
   amplitude += fitpara[32]*(-1.0)*(SF_PT01_S_TP*RBW_f21270_pm*RBW_pi21670_01 + SF_PT02_S_TP*RBW_f21270_pm*RBW_pi21670_02);
   
   // D -> 1--(rho pi) pi0
   double SF_V1_Vz  = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi01), T1_PipPimPi01), T1_PipPim),  contract_21_1(Proj1_3z1, T1_3z12));
   double SF_V1_Vp1 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi01), T1_PipPi01Pim), T1_PipPi01), contract_21_1(Proj1_3z1, T1_3z12));
   double SF_V1_Vm1 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi01), T1_PimPi01Pip), T1_PimPi01), contract_21_1(Proj1_3z1, T1_3z12));
 
   double SF_V2_Vz  = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi02), T1_PipPimPi02), T1_PipPim),  contract_21_1(Proj1_3z2, T1_3z21));
   double SF_V2_Vp2 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi02), T1_PipPi02Pim), T1_PipPi02), contract_21_1(Proj1_3z2, T1_3z21));
   double SF_V1_Vm2 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi02), T1_PimPi02Pip), T1_PimPi02), contract_21_1(Proj1_3z2, T1_3z21));
 
 
   // D0 -> omega(rho pi) pi0
   amplitude += (-1.0) * fitpara[33]*(SF_V1_Vp1*RBW_omega_01*GS_rho770_p1 - SF_V1_Vz*RBW_omega_01*GS_rho770_pm - SF_V1_Vm1*RBW_omega_01*GS_rho770_m1 + SF_V2_Vp2*RBW_omega_02*GS_rho770_p2 - SF_V2_Vz*RBW_omega_02*GS_rho770_pm - SF_V1_Vm2*RBW_omega_02*GS_rho770_m2);
 
 
   // D0 -> phi(rho pi) pi0
   amplitude += (-1.0) * fitpara[34]*(SF_V1_Vp1*RBW_phi_01*GS_rho770_p1 - SF_V1_Vz*RBW_phi_01*GS_rho770_pm - SF_V1_Vm1*RBW_phi_01*GS_rho770_m1 + SF_V2_Vp2*RBW_phi_02*GS_rho770_p2 - SF_V2_Vz*RBW_phi_02*GS_rho770_pm - SF_V1_Vm2*RBW_phi_02*GS_rho770_m2);
 
   return amplitude;
 
}
 
int EvtD0Topippim2pi0::CalAmp(){
 
   m_AmpD0 = CalD0Amp();
   m_AmpDb = CalDbAmp();
 
   return 0;
}
 
double EvtD0Topippim2pi0::mag2(complex<double> x)
{
   double temp = x.real()*x.real() + x.imag()*x.imag();
   return temp;
}
 
double EvtD0Topippim2pi0::arg(complex<double> x)
{
   double temp = atan(x.imag()/x.real());
   if(x.real()<0) temp=temp+TMath::Pi();
   return temp;
}
 
double EvtD0Topippim2pi0::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 EvtD0Topippim2pi0::AmplitudeSquare(int charm, int tagmode){
 
   EvtVector4R dp1=GetDaugMomLab(0),dp2=GetDaugMomLab(1),dp3=GetDaugMomLab(2), dp4=GetDaugMomLab(3); // pi+ pi- pi0 pi0
   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 p4pip[4], p4pim[4], p4pi01[4], p4pi02[4];
   for(int i=0; i<3; i++){
      p4pip[i]=dp1bst.get(i+1);
      p4pim[i]=dp2bst.get(i+1);
      p4pi01[i]=dp3bst.get(i+1);
      p4pi02[i]=dp4bst.get(i+1);
   }
   p4pip[3]=dp1bst.get(0);
   p4pim[3]=dp2bst.get(0);
   p4pi01[3]=dp3bst.get(0);
   p4pi02[3]=dp4bst.get(0);
 
   setInput(p4pip, p4pim, p4pi01, p4pi02);
   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;
      double 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;
}