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




}

void calc_3g_yields(int mybin,double yields[],double uncs[]){
  TH2F *h3g=(TH2F*)_file0->FindObjectAny("ThreeGammaMass_vs_Ebeam");

  TF1 *f1=new TF1("f1",my3gfunc,0,2.5,11);
 
  // TCanvas *c1=new TCanvas("c1","c1");

  char name[80];
  sprintf(name,"slice%d",mybin);
    
  
  TH1D *slice=h3g->ProjectionY(name,mybin,mybin);
  slice->GetXaxis()->SetRangeUser(0.75,0.83);
  slice->Fit("gaus","q");
  if (slice->GetFunction("gaus")){
    f1->SetParameter(0,slice->GetFunction("gaus")->GetParameter(0));
  }
  f1->SetParameter(1,0.783);
  f1->SetParameter(2,0.01);
  f1->SetParameter(3,1);
  f1->SetParameter(4,-0.05);
  f1->SetParameter(5,-0.2); 
  slice->GetXaxis()->SetRangeUser(0.96,1.08);
  slice->Fit("gaus","q"); 
  if (slice->GetFunction("gaus")){
    f1->SetParameter(6,slice->GetFunction("gaus")->GetParameter(0)); 
    f1->SetParameter(7,1.02);
    f1->SetParameter(8,0.1*f1->GetParameter(0));
    f1->SetParameter(9,0.783);
    f1->SetParameter(10,0.02);
  }

  slice->GetXaxis()->SetRangeUser(0.7,1.1);
  slice->Fit("f1");
  
  TF1 *f2=new TF1("f2",my3gfunc,0,2.5,11);
  f2->SetParameter(3,f1->GetParameter(3));
  f2->SetParameter(4,f1->GetParameter(4));
  f2->SetParameter(5,f1->GetParameter(5)); 
  TF1 *f3=new TF1("f3",my3gfunc,0,2.5,11); 
  f3->SetParameter(0,f1->GetParameter(0));
  f3->SetParameter(1,f1->GetParameter(1));
  f3->SetParameter(2,f1->GetParameter(2)); 
  TF1 *f4=new TF1("f4",my3gfunc,0,2.5,11); 
  f4->SetParameter(6,f1->GetParameter(6));
  f4->SetParameter(7,f1->GetParameter(7));
  f4->SetParameter(2,f1->GetParameter(2));
		  
  slice->GetXaxis()->SetRangeUser(0.5,1.5);
  
  char myname[80];
  sprintf(myname,"3#gamma 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"); 
  f4->SetLineColor(5);
  f4->Draw("same"); 


  //yields[1]=f3->Integral(0.7,0.9)/0.0025; 
  double yield_unc=0.;
  yields[1]=slice->IntegralAndError(281,360,yield_unc)-f2->Integral(0.7,0.9)/0.0025; 
  uncs[1]=sqrt(yield_unc*yield_unc+fabs(f2->Integral(0.7,0.9))/0.0025); 
  // yields[2]=f4->Integral(0.96,1.08)/0.0025; 
  yields[2]=slice->IntegralAndError(385,432,yield_unc)-f2->Integral(0.96,1.08)/0.0025;
  uncs[2]=sqrt(yield_unc*yield_unc+fabs(f2->Integral(0.96,1.08))/0.0025);
  cout <<"Omega: " << yields[1]<< "+/-" << uncs[1] << endl; 
  cout <<"Phi:   " << yields[2]<< "+/-" << uncs[2] << endl; 

 }