double myfunc(double x[],double par[]){
  double f=0;
  double thresh=0.548+0.14+0.14;
  if (x[0]>thresh){
    f+=par[0]*pow(x[0]-thresh,par[1])*exp(1.+par[2]*x[0]+par[3]*x[0]*x[0]);
  }
  
  //par[5]=0.55;
  //par[2]=0.01;
  f+=par[4]*TMath::Gaus(x[0],par[5],par[6]);
  f+=par[7]*TMath::Voigt(x[0]-par[8],par[6],0.024); 
  // par[11]=1.019;
  //  par[7]=14000;
  //  par[8]=1.02;
  //par[9]=0.015;
  f+=par[10]*TMath::Voigt(x[0]-par[11],par[6],par[12]);
  f+=par[13]*TMath::Gaus(x[0],par[14],par[15]); 
  //f+=par[16]*TMath::Gaus(x[0],par[17],par[18]);
  
  return f;
}

void fit_pippimeta(){
  TF1 *f1=new TF1("f1",myfunc,0,3.5,16);
  f1->SetParameter(0,5);
  f1->SetParameter(1,1);
  f1->SetParameter(2,-1);
  f1->SetParameter(3,-1);
  f1->SetParameter(13,120);
  f1->SetParameter(14,1.65);
  f1->SetParameter(15,0.5);
  /*  f1->SetParameter(16,1000);
  f1->SetParameter(17,1.65);
  f1->SetParameter(18,0.25); 
  f1->SetParameter(19,1000);
  f1->SetParameter(20,1.25);
  f1->SetParameter(21,0.25);
  */
  // TH2D *hist=(TH2D *)_file0->FindObjectAny("PipPimEtaMass_vs_Ebeam");
  //TH1D *hproj=(TH1D *)hist->ProjectionY("hproj",17,-1);
  //hproj->Draw();
  // hproj->Fit(f1);
  // PipPimEtaMass_kf->Fit("f1");
  f1->SetParameter(4,800);
  f1->SetParameter(5,0.958);
  f1->SetParameter(6,0.02);
  f1->SetParameter(7,500);
  f1->SetParameter(8,1.285);
  f1->SetParameter(9,0.024);
  f1->SetParameter(10,100);
  f1->SetParameter(11,1.65);
  f1->SetParameter(12,0.15);
  PipPimEtaMass_kf->Fit("f1");
  PipPimEtaMass_kf->Fit("f1");
  /*f1->SetParameter(4,0);
  f1->SetParameter(7,0);
  //f1->SetParameter(10,0);
  f1->SetLineStyle(2);
  f1->Draw("same");
*/
  double eta_all=PipPimEtaMass_kf->Integral(105,116);
  double eta_back=f1->Integral(0.52,0.58)/0.005;
  double omega_all=PipPimEtaMass_kf->Integral(149,166);
  double omega_back=f1->Integral(0.74,0.83)/0.005;
  double phi_all=PipPimEtaMass_kf->Integral(201,208);
  double phi_back=f1->Integral(1.0,1.04)/0.005;
 

  cout << eta_all-eta_back << "+/-"<<sqrt(eta_all+eta_back)<< " "
       << omega_all-omega_back << "+/-"<<sqrt(omega_all+omega_back) << " "
       << phi_all-phi_back << "+/-"<<sqrt(phi_all+phi_back)
       <<endl;


}