double my3pifunc(double *x,double *par){
  double f=0.;
  double scale1=TMath::Voigt(0.,par[0],0.00849);
  f+=par[1]*TMath::Voigt(x[0]-par[2],par[0],0.00849)/scale1; 
  double scale2=TMath::Voigt(0.,par[0],0.00426);
  f+=par[3]*TMath::Voigt(x[0]-par[4],par[0],0.00426)/scale2; 
  if (x[0]>0.7){
    f+=par[5]*pow(x[0]-0.7,3.)*exp(1.+par[6]*x[0]+par[7]*x[0]*x[0]);
  }
  //  f+=par[8]*TMath::Gaus(x[0],0.8,par[10]);
  return f;




}

void calc_3pi_yields(int mybin,double yields[],double uncs[]){
  TH2F *h3pi=(TH2F*)_file0->FindObjectAny("PipPimPi0Mass_vs_Ebeam");

  TF1 *f1=new TF1("f1",my3pifunc,0,2.5,8);
 
  TCanvas *c1=new TCanvas("c1","c1");

  char name[80];
  sprintf(name,"slice%d",mybin);
    
  
  TH1D *slice=h3pi->ProjectionY(name,mybin,mybin);
  slice->GetXaxis()->SetRangeUser(0.5,0.6);
  slice->Fit("gaus","q");
  TF1 *f4=new TF1("f4","[0]*TMath::Gaus(x,[1],[2])",0.5,0.6);
  if (slice->GetFunction("gaus")){  
    f1->SetParameter(0,slice->GetFunction("gaus")->GetParameter(2));
    f4->SetParameter(0,slice->GetFunction("gaus")->GetParameter(0));
    f4->SetParameter(1,slice->GetFunction("gaus")->GetParameter(1));
    f4->SetParameter(2,slice->GetFunction("gaus")->GetParameter(2));
  }
  slice->GetXaxis()->SetRangeUser(0.75,0.83);
  slice->Fit("gaus","q");
  if (slice->GetFunction("gaus")){
    f1->SetParameter(1,slice->GetFunction("gaus")->GetParameter(0));
    f1->SetParameter(2,0.783);
  }
  slice->GetXaxis()->SetRangeUser(0.95,1.05);
  slice->Fit("gaus","q");
  if (slice->GetFunction("gaus")){
    f1->SetParameter(3,slice->GetFunction("gaus")->GetParameter(0));
    f1->SetParameter(4,1.02);
  }
  f1->SetParameter(5,1000);
  f1->SetParameter(6,-0.05);
  f1->SetParameter(7,-0.2);
  /*  f1->SetParameter(8,0.2*f1->GetParameter(1));
  f1->SetParameter(9,0.8);
  f1->SetParameter(10,0.02);
  */
  slice->GetXaxis()->SetRangeUser(0.7,1.5);
  slice->Fit("f1");
  
  TF1 *f2=new TF1("f2",my3pifunc,0,2.5,8);
  f2->SetParameter(5,f1->GetParameter(5));
  f2->SetParameter(6,f1->GetParameter(6));
  f2->SetParameter(7,f1->GetParameter(7));
  //f2->SetParameter(8,f1->GetParameter(8));
  //f2->SetParameter(9,f1->GetParameter(9));
  //f2->SetParameter(10,f1->GetParameter(10));
  TF1 *f3=new TF1("f3",my3pifunc,0,2.5,8); 
  f3->SetParameter(0,f1->GetParameter(0));
  f3->SetParameter(1,f1->GetParameter(1));
  f3->SetParameter(2,f1->GetParameter(2)); 
 
  TF1 *f5=new TF1("f5",my3pifunc,0,2.5,8); 
  f5->SetParameter(0,f1->GetParameter(0));
  f5->SetParameter(3,f1->GetParameter(3));
  f5->SetParameter(4,f1->GetParameter(4));
		  
  slice->GetXaxis()->SetRangeUser(0.5,1.5);
  
  char myname[80];
  sprintf(myname,"#pi^{+}#pi^{-}#pi^{0} yield for E_{#gamma}=%3.1f GeV",2.9+0.2*(mybin-1));
  slice->SetTitle(myname);
  slice->SetYTitle("counts / 0.0025 GeV");
  slice->SetMinimum(0.);
  slice->Draw();
  f2->SetLineColor(4);
  f2->Draw("same");
  f3->SetLineColor(3);
  f3->Draw("same"); 
  f5->SetLineColor(7);
  f5->Draw("same");
  f4->SetLineColor(1);
  f4->Draw("same");

  //  yields[0]=f4->Integral(0.5,0.6)/0.0025;
  double yield_unc=0.;
  yields[0]=slice->IntegralAndError(201,240,yield_unc)-f3->Integral(0.5,0.6)/0.0025;
  uncs[0]=sqrt(yield_unc*yield_unc+f3->Integral(0.5,0.6)/0.0025);
  // yields[1]=f3->Integral(0.7,0.9)/0.0025;
  yields[1]=slice->IntegralAndError(281,360,yield_unc)-f2->Integral(0.7,0.9)/0.0025;
  uncs[1]=sqrt(yield_unc*yield_unc+f2->Integral(0.7,0.9)/0.0025);
  //yields[2]=f5->Integral(0.96,1.08)/0.0025;
  yields[2]=slice->IntegralAndError(385,432,yield_unc)-f2->Integral(0.96,1.08)/0.0025
    -f3->Integral(0.96,1.08)/0.0025;
  uncs[2]=sqrt(yield_unc*yield_unc+fabs(f2->Integral(0.96,1.08))/0.0025
	       +fabs(f3->Integral(0.96,1.08))/0.0025);


  cout <<"Eta:   " << yields[0] << "+/-" << uncs[0] << endl;
  cout <<"Omega: " << yields[1]<< "+/-" << uncs[1] << endl; 
  cout <<"Phi:   " << yields[2]<< "+/-" << uncs[2] << endl; 
 }