//calculate Plab for given E in CM [GeV]
const double eff_sig = 0.70;
const double eff_bkg = 1e-7;
const double br_sig = 5e-3;
const double sig_sig = 1e-7;
const double sig_bkg = 5e-2;
const double Lumi = 1e8; //1e32*1e-24;
const double mp = 0.938272;//proton mass, GeV
const double mp2 = TMath::Power(mp,2);
const double mh = 3.52593;
const double G = 7e-4;
const double G2 = TMath::Power(G,2);
Double_t time_sig(Double_t *x, Double_t *par){
  double N = eff_sig*br_sig*sig_sig*Lumi*x[0];
  return N;
}

Double_t time_bkg(Double_t *x, Double_t *par){
  double N = eff_bkg*sig_bkg*Lumi*x[0];
  return N;
}

Double_t time_5Sigma(Double_t *x, Double_t *par){
  double sig_contr = eff_sig*br_sig*sig_sig;
  double bkg_contr = eff_bkg*sig_bkg;
  double N = 25.*(sig_contr+bkg_contr)/(sig_contr*sig_contr*Lumi);
  return N;
}

Double_t time_SIGevents(Double_t *x, Double_t *par){
  double N = par[0]/(eff_sig*br_sig*sig_sig*Lumi);
  return N;
}




Double_t myfunction(Double_t *x, Double_t *par){
  //norm factor to have 1e4 events in the peak: 1.64652418744031340e+05 Br=5e-3
  //norm factor to have 1e4 events in the peak: 2.05815523430039204e+05  Br=4e-3
  //norm facror to have 69469.4 events in the peak: 1.42978809233707632e+06 (Br=4e-3)
  double Ecm1 = x[0];
  double Ecm12 = TMath::Power(Ecm1,2)/(2*mp);
  double first_coef = 3*4*TMath::Pi()/(Ecm12-4*mp2);
  double Er = Ecm1 - mh;
  double Er2 = TMath::Power(Er,2);
  double reson_depend = br_sig*G2/(4*Er2+G2);
  double sigma_BW = first_coef*reson_depend;
  //sigma_BW *=1.64652418744031340e+05;
  //  sigma_BW *= 2.05815523430039204e+05;
  // sigma_BW *=50./0.060734;
  sigma_BW *=1.42978809233707632e+06;
  return sigma_BW;
  //  sigma_BW *=50./0.060734;
}

void N_BW(double Ecm=3.52593){
 
  double Ecm2 = TMath::Power(Ecm,2)/(2*mp);
  //  cout<<" Ecm2 = "<<Ecm2<<endl;
  double term = Ecm2-mp;
  double plab2 = TMath::Power(term,2)-mp2;
  double plab = TMath::Power(plab2,0.5);
  cout<<"For Ecm = "<<Ecm<<" Plab = "<<plab<<endl;
  double xn[2],parn[2];

  double t5s = time_5Sigma(xn,parn);
  cout<<"t for 5 sigma significance = "<<t5s<<" = "<<t5s/(60.*60.*24.*7)<<endl;
  xn[0] = t5s;
  cout<<"Number of signal events = "<<time_sig(xn,parn)<<" Number of bkg events = "<<time_bkg(xn,parn)<<endl;
  parn[0] = 1e4*eff_sig;
  double tSIG = time_SIGevents(xn,parn);
  cout<<"t for "<<parn[0]<<" signal events = "<<tSIG<<endl;
  xn[0] = tSIG;
  cout<<"Number of signal events = "<<time_sig(xn,parn)<<" Number of bkg events = "<<time_bkg(xn,parn)<<endl;
  TF1 *f0 = new TF1("timeSig",time_sig,0,1e2,0);
  f0->GetYaxis()->SetTitle("N^{REC}_{ev}");
  f0->GetXaxis()->SetTitle("time, s");
  f0->SetLineColor(2);
  TF1 *f01 = new TF1("timeBkg",time_bkg,0,1e2,0);
  // f01->GetYaxis()->SetTitle("N^{REC}_{ev}");
  // f01->GetXaxis()->SetTitle("time, s");
  f01->SetLineColor(3);
  // f01->Draw();
  // f0->Draw("same");

  TF1 *f1 = new TF1("myfunc",myfunction,3.51,3.54,0);
  f1->GetXaxis()->SetTitle("E_{CM}, GeV");
  f1->GetYaxis()->SetTitle("N_{ev}");
  //  f1->Draw();
  double x[2];
  x[0] = Ecm;
  double par[2];
  double sigma_BW = myfunction(x,par);
  cout<<"sigma_BW = "<<sigma_BW<<endl;
  // cout<<"Number of events to have 1e4 events in the peak: "<<int(sigma_BW*1.64652418744031340e+05)<<endl;
 // double mh = 3.52593;
 // double G = 7e-4;
  double Ecmpoins[20], nBW[20];
  double errEcmpoins[20];
  for(int i=0;i<20;i++){
    Ecmpoins[i] = 3.52593-i*0.5*G+5*G;
    double x1[2];
    x1[0] =  Ecmpoins[i];
    nBW[i] = int(myfunction(x1,par));
    //  nBW[i] += 1.92e+08;
    Ecm2 = TMath::Power(Ecmpoins[i],2)/(2*mp);
    term = Ecm2-mp;
    plab2 = TMath::Power(term,2)-mp2;
   plab = TMath::Power(plab2,0.5);
   cout<<" Ecm = "<< Ecmpoins[i]<<" plab = "<<plab<<" nBW = "<<nBW[i]<<endl;
   errEcmpoins[i] = 1e-5;
  }
  TGraph *gr = new TGraph(20,Ecmpoins,nBW);
 gr->SetMarkerColor(2);
 gr->SetMarkerStyle(20);
 gr->Draw("AP");
 double Nexp[20]={72, 89, 76, 141, 207, 272, 404, 692, 1367, 3456, 6764, 3513, 1389, 635, 409, 218, 222, 121, 113, 82};
 //double Nexp[20]={33,60,63,85,106,164,229,427,879,2130,4318,2190,855,408,257,176,110,84,72,52};
 for(int i=0;i<20;i++){
   Nexp[i] /=5.64965996451803620e-02;
   Nexp[i] += 5459.18;
 }
 TGraphErrors *gr2 = new TGraphErrors(20,Ecmpoins,Nexp);
 gr2->SetMarkerColor(4);
 gr2->SetMarkerStyle(21);
 //gr2->Draw("AP");
 // TMultiGraph *mrg = new TMultiGraph();
 // mrg->Add(gr);
 // mrg->Add(gr2);
 // mrg->Draw("AP");
 // double eff_sig = 0.45;
 // double br_sig = 4e-3;
 // double sig_sg = 100.0*1e-9;
 // double Lumi = 1e32*1e-24;
 // double time_5sig = 1.51200000000000000e+07;
 // double Nsig = eff_sig*Lumi*sig_sg*br_sig*time_5sig;
 // cout<<" Nsig = "<<Nsig<<endl;
 // double time_Nsig = 1e4/(eff_sig*Lumi*sig_sg*br_sig);
 // cout<<" time_Nsig (1e4) = "<<time_Nsig<<endl;
 // double eff_bkg = 1e-7;
 // double sig_bkg = 50*1e-3;
 // double Nbkg = eff_bkg*sig_bkg*Lumi*time_5sig;
 // cout<<"Nbkg = "<<Nbkg<<endl;
}