double myfunc(double x[],double par[]){
  double f=0;
  double thresh=0.135+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.00849); 
  // par[11]=1.019;
  f+=par[9]*TMath::Voigt(x[0]-par[10],par[6],0.00424);
  f+=par[11]*TMath::Voigt(x[0]-par[12],par[6],0.110);
  //  par[7]=14000;
  //  par[8]=1.02;
  //par[9]=0.015;
  
  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_pippimpi0(){
  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,1200);
  f1->SetParameter(14,1.6);
  f1->SetParameter(15,0.25);
  /*  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("PipPimPi0Mass_vs_Ebeam");
  //TH1D *hproj=(TH1D *)hist->ProjectionY("hproj",17,-1);
  //hproj->Draw();
  // hproj->Fit(f1);
  PipPimPi0Mass_kf->Fit("f1");
  f1->SetParameter(4,6000);
  f1->SetParameter(5,0.55);
  f1->SetParameter(6,0.002);
  f1->SetParameter(7,50000);
  f1->SetParameter(8,0.783);
  f1->SetParameter(9,200000);
  f1->SetParameter(10,1.02);
  f1->SetParameter(11,10000);
  f1->SetParameter(12,1.32);
  PipPimPi0Mass_kf->Fit("f1");

  f1->SetParameter(4,0);
  f1->SetParameter(7,0);
  f1->SetParameter(9,0);
  f1->SetLineStyle(2);
  f1->Draw("same");

  double eta_all=PipPimPi0Mass_kf->Integral(105,116);
  double eta_back=f1->Integral(0.52,0.58)/0.005;
  double omega_all=PipPimPi0Mass_kf->Integral(149,166);
  double omega_back=f1->Integral(0.74,0.83)/0.005;
  double phi_all=PipPimPi0Mass_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;


}