double mykpkmfunc(double *x,double *par){
  double f=0.;
  double mK=0.494;
  if (x[0]>2.*mK){
    // if (false){
      //f+=par[0]*TMath::Gaus(x[0],par[1],par[2]);
    f+=par[0]*TMath::Voigt(x[0]-par[1],par[2],0.00425);
      /*
    }
    else{
       f+=par[0]*TMath::Gaus(x[0],par[1],par[5]);
    }
      */
    //    f+=par[6]*TMath::Gaus(x[0],par[7],par[8]);
    if (par[3]>0)
      f+=par[3]*pow(x[0]-2.*mK,par[4]);
  }
  //f+=par[8]*TMath::Gaus(x[0],par[9],par[10]);
  return f;




}

void calc_kpkm_yield(int mybin,double &yield,double &unc){
  TH2F *hkpkm=(TH2F*)_file0->FindObjectAny("KpKmMass_vs_Ebeam");

  TF1 *f1=new TF1("f1",mykpkmfunc,0.99,1.1,5);
  f1->SetParameters(1.,1.02,0.005,50,0.5);
  
  
  TCanvas *c1=new TCanvas("c1","c1");

  char name[80];
  sprintf(name,"slice%d",mybin);
 
  
  TH1D *slice=hkpkm->ProjectionY(name,mybin,mybin); 
  f1->SetParameter(0,slice->GetMaximum()/100.);

  slice->GetXaxis()->SetRangeUser(0.98,1.1);  
  f1->SetParameter(1,slice->GetBinCenter(slice->GetMaximumBin()));
  f1->SetParameter(0,0.1*slice->GetMaximum());
  f1->SetParameter(2,0.001);
  f1->SetParameter(3,1.);
  f1->SetParameter(4,0.5);
  slice->Fit("f1","r");
  
  char myname[80];
  sprintf(myname,"K^{+}K^{-} yield for E_{#gamma}=%3.1f GeV",2.9+0.2*(mybin-1));
  slice->SetTitle(myname);
  slice->SetYTitle("counts / 0.001 GeV");
  slice->SetMinimum(0.);
  slice->Draw();

  TF1 *f2=new TF1("f2",mykpkmfunc,0.99,1.1,6);

  f2->SetParameter(4,f1->GetParameter(4));
  f2->SetParameter(5,f1->GetParameter(5));

  f2->SetLineStyle(2);
  f2->Draw("same");
  
  double yield_unc=0.;
  yield=slice->IntegralAndError(41,110,yield_unc)-f2->Integral(0.99,1.06)/0.001;
  unc=sqrt(yield_unc*yield_unc+f2->Integral(0.99,1.06)/0.001);

  cout <<"Phi:   " << yield << "+/-" << unc << endl; 
  
 }