double myfunc(double *x,double *par){
  double m=x[0];  
  double m_pi=0.13957;  
  double f=0;
  double thresh=m_pi+m_pi;
  if (m>thresh){
    f+=par[0]*pow(m-thresh,par[1])*exp(1.+par[2]*m+par[3]*m*m);
  
    double m_rho=par[5];
    double a=m*m-m_rho*m_rho;
    double q2=0.25*(m*m-4.*m_pi*m_pi);
    double q2_rho=0.25*(m_rho*m_rho-4.*m_pi*m_pi);
    double R2=par[6];
    double Gamma_rho=par[7]*pow(q2/q2_rho,1.5)*m_rho/m*(1.+R2*q2_rho)/(1.+R2*q2);
    //Gamma_rho*=0.8;
    double b=m_rho*Gamma_rho;
    double b2=b*b;
    double a2=a*a;
    f+=par[4]*b2/q2/(a2+b2);
  
    double m_omega=par[9];
  //par[4]=2.;
    double R2o=par[10]; 
    double ao=m*m-m_omega*m_omega;
    double q2_omega=0.25*(m_omega*m_omega-4.*m_pi*m_pi);
    double Gamma_omega=par[11]*pow(q2/q2_omega,1.5)*m_omega/m*(1.+R2o*q2_omega)/(1.+R2o*q2);
    //  Gamma_omega*=4.;
    double bo=m_omega*Gamma_omega;
    double b2o=bo*bo;
    double a2o=ao*ao;
    f+=par[8]*b2o/q2/(a2o+b2o);
    f+=2.*sqrt(par[4]*par[8])*b*bo/q2/(a2+b2)/(a2o+b2o)*(cos(par[12])*(a*ao-b*bo)+sin(par[12])*(a*bo-ao*b)); 
    
    
    double m_f2=par[14];
    double af=m*m-m_f2*m_f2;
    double q2_f2=0.25*(m_f2*m_f2-4.*m_pi*m_pi);
    double z=q2*par[15];
    double z0=q2_f2*par[15];
    double Gamma_f2=par[16]*pow(q2/q2_f2,2.5)*m_f2/m
      *((z0-3)*(z0-3)+9*z0)/((z-3)*(z-3)+9*z);
    double bf=m_f2*Gamma_f2;
    double b2f=bf*bf;
    double a2f=af*af;
    f+=par[13]*b2f/q2/(a2f+b2f);

    double m_rp=par[18];
    double R2rp=par[19]; 
    double arp=m*m-m_rp*m_rp;
    double q2_rp=0.25*(m_rp*m_rp-4.*m_pi*m_pi);
    double Gamma_rp=par[20]*pow(q2/q2_rp,1.5)*m_omega/m*(1.+R2rp*q2_rp)/(1.+R2rp*q2);
    //  Gamma_omega*=4.;
    double brp=m_rp*Gamma_rp;
    double b2rp=brp*brp;
    double a2rp=arp*arp;
    f+=par[17]*b2rp/q2/(a2rp+b2rp);
  }

  return f;
}

void fit_rho(){
  TF1 *f1=new TF1("f1",myfunc,0.28,2.,21);
  f1->SetParameters(5000,5.,1.8,-2.);
  f1->SetParameter(4,10000.);
  f1->SetParameter(5,0.76);
  f1->SetParameter(6,4.);
  f1->SetParameter(7,0.15);
  f1->SetParameter(8,1000.);
  f1->SetParameter(9,0.78);
  f1->SetParameter(10,2.);
  f1->SetParameter(11,0.008);
  f1->SetParameter(12,1.);
  f1->SetParameter(13,150);
  f1->SetParameter(14,1.27);
  f1->SetParameter(15,1.4);
  f1->SetParameter(16,0.2);
  f1->SetParameter(17,30);
  f1->SetParameter(18,1.65);
  f1->SetParameter(19,4.);
  f1->SetParameter(20,0.3);
  //f1->Draw("same");
  f1->SetRange(0.38,1.8);
  PipPimMass_kf->Fit("f1","r");
  //// f1->SetRange(0.28,2.);
  //PipPimMass_kf->Fit("f1","r");
  TF1 *f2=f1;
  TF1 *f3=f1;
  f2->SetParameter(0,0);
  f2->SetParameter(8,0);
  f2->SetParameter(13,0);
  f2->SetParameter(17,0);
  f2->SetLineColor(4);
  f2->Draw("same"); 
  
  cout << f2->Integral(0.28,4.)/0.0025 << endl;

  /*
  f3->SetParameter(0,0);
  f3->SetParameter(4,0);
  f2->SetParameter(13,0);
  f2->SetParameter(17,0);
  f3->SetLineColor(6);
  f3->Draw("same");
  */

}