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.00849);
  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],par[9],par[10]);
  return f;




}

void calc_3pi_vs_t_yields(int myEbin,int mytbin,double yields[],double uncs[]){
  char myhist[80];
  sprintf(myhist,"PipPimPi0_vs_t_%3.1fGeV",7.3+0.6*myEbin);
  TH2F *h3pi=(TH2F*)_file0->FindObjectAny(myhist);

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

    char name[80];
    sprintf(name,"slice%d",mytbin);
    
    
    TH1D *slice=h3pi->ProjectionY(name,mytbin,mytbin+1);
    slice->GetXaxis()->SetRangeUser(0.5,0.6);
    slice->Fit("gaus","q");  
    f1->SetParameter(0,0.01);
    //    f1->SetParameter(0,slice->GetFunction("gaus")->GetParameter(2));
    TF1 *f4=new TF1("f4","[0]*TMath::Gaus(x,[1],[2])",0.5,0.6);
    if (slice->GetFunction("gaus")){
      f4->SetParameter(0,slice->GetFunction("gaus")->GetParameter(0));
      f4->SetParameter(1,slice->GetFunction("gaus")->GetParameter(1));
      f4->SetParameter(2,slice->GetFunction("gaus")->GetParameter(2));
    }
	else{
	f4->SetParameter(0,1.);
	f4->SetParameter(1,0.55);
	f4->SetParameter(2,0.01);
	}
    slice->GetXaxis()->SetRangeUser(0.75,0.83);
    slice->Fit("gaus","q");
    if (slice->GetFunction("gaus")){
      f1->SetParameter(1,slice->GetFunction("gaus")->GetParameter(0));
	}
else f1->SetParameter(1,1.);
      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));
	}
	else f1->SetParameter(3,1.);
    f1->SetParameter(4,1.02);
    
    f1->SetParameter(5,1000);
    f1->SetParameter(6,-0.05);
    f1->SetParameter(7,-0.2);
    //.f1->SetParameter(8,0.1*f1->GetParameter(1));
    //f1->SetParameter(9,0.78);
    //f1->SetParameter(10,0.02);
    
    slice->GetXaxis()->SetRangeUser(0.7,1.5);
    slice->Fit("f1");
    
    TF1 *f2=new TF1("f2",my3pifunc,0,2.5,11);
    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,11); 
    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,11); 
    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 -t=%5.3f GeV^{2}",0.02+0.02*(mytbin-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"); 
    f5->SetLineColor(7);
    f5->Draw("same");
    f4->SetLineColor(1);
    f4->Draw("same");
    
     yields[0]=0.,yields[1]=0.,yields[2]=0.;
     uncs[0]=0.,uncs[1]=0.,uncs[2]=0.;
    if (slice->Integral(101,120)>0. && f3->Integral(0.5,0.6)>0.){
    yields[0]=slice->Integral(101,120)-f3->Integral(0.5,0.6)/0.005;
    uncs[0]=sqrt(slice->Integral(101,120)+f3->Integral(0.5,0.6)/0.005);
}
    // yields[1]=f3->Integral(0.7,0.9)/0.005;
    if (slice->Integral(141,180)>0. && f2->Integral(0.7,0.9)>0.){
    yields[1]=slice->Integral(141,180)-f2->Integral(0.7,0.9)/0.005;
    uncs[1]=sqrt(slice->Integral(141,180)+f2->Integral(0.7,0.9)/0.005);
	}
    //yields[2]=f5->Integral(0.96,1.08)/0.005;
    if (slice->Integral(193,216)>0. && f2->Integral(0.96,1.08)>=0.
	&& f3->Integral(0.96,1.08)>=0.
	){
    yields[2]=slice->Integral(193,216)-f2->Integral(0.96,1.08)/0.005
      -f3->Integral(0.96,1.08)/0.005;
    uncs[2]=sqrt(slice->Integral(193,216)+f2->Integral(0.96,1.08)/0.005
		 +f3->Integral(0.96,1.08)/0.005);
 }
if (yields[0]<0.){ yields[0]=0.;
	uncs[0]=sqrt(slice->Integral(101,120)/0.005);
}
if (yields[1]<0.){ yields[1]=0.;
	uncs[1]=sqrt(slice->Integral(141,180)/0.005);
	}
if (yields[2]<0.){ yields[2]=0.;
	uncs[2]=sqrt(slice->Integral(193,216)/0.005);
	}

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