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; 
  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]);
  }
  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,6);
 
  // 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");
  f1->SetParameter(0,slice->GetFunction("gaus")->GetParameter(0));
  f1->SetParameter(1,0.783);
  f1->SetParameter(2,0.01);
  f1->SetParameter(4,1);
  f1->SetParameter(5,-0.05);
  f1->SetParameter(5,-0.2);


  slice->GetXaxis()->SetRangeUser(0.7,1.1);
  slice->Fit("f1");
  
  TF1 *f2=new TF1("f2",my3gfunc,0,2.5,10);
  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,10); 
  f3->SetParameter(0,f1->GetParameter(0));
  f3->SetParameter(1,f1->GetParameter(1));
  f3->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.005 GeV");
  slice->SetMinimum(0.);
  slice->Draw();
  f2->SetLineColor(4);
  f2->Draw("same");
  f3->SetLineColor(3);
  f3->Draw("same"); 

  yields[1]=f3->Integral(0.7,0.9)/0.005;
  uncs[1]=sqrt(slice->Integral(141,180)+f2->Integral(0.7,0.9)/0.005);
  cout <<"Omega: " << yields[1]<< "+/-" << uncs[1] << endl; 

 }