// Kshort
demo->cd();

hm4rec2OSpipitt05acc->Draw();
// do not use range before fitting
//hm4rec2OSpipitt05acc->GetXaxis()->SetRangeUser(0.,2.0);
hm4rec2OSpipitt05acc->GetXaxis()->SetTitle("M (GeV/c^{2})");
hm4rec2OSpipitt05acc->GetXaxis()->CenterTitle(true);
hm4rec2OSpipitt05acc->GetYaxis()->SetTitle("#frac{d#sigma}{dM} (pb/0.01GeV/c^{2})");
hm4rec2OSpipitt05acc->GetYaxis()->CenterTitle(true);

//TH1D* suma = new TH1D( "suma", "suma", 2500, 0, 5 );
//auto suma = *hm4rec2OSpipitt05acc + *hm4rec2OSpiptt05acc;
//suma.Draw();

gStyle->SetOptStat(1111111);


//Breit-Wigner function
//                 .....it must be here ) {
Double_t mybw(Double_t* x, Double_t* par) {
  Double_t arg1 = 14.0/22.0;                     // 2 over pi
  Double_t arg2 = par[1]*par[1]*par[2]*par[2];   // Gamma=par[1]  M=par[2]  Factor=par[0]
  Double_t arg3 = ((x[0]*x[0]) - (par[2]*par[2]))*((x[0]*x[0]) - (par[2]*par[2]));
  Double_t arg4 = x[0]*x[0]*x[0]*x[0]*((par[1]*par[1])/(par[2]*par[2]));
  return par[0]*arg1*arg2/(arg3 + arg4);
}

//non-resonant background (4pi or lampipi) of the form ...see Robert's draft
// A · (m - B)^C · exp^(D·m)  *** 4 parameters
  Double_t nonresonantbackground(Double_t *x, Double_t *par) {
  return par[0]*pow((x[0]-par[1]),par[2])*TMath::Exp(par[3]*x[0]);
}

// Sum of non-resonant background and breit-wigner peak functions - Kshort channel
  Double_t mylamfit(Double_t *x, Double_t *par) {
  return  mybw(x,par) + nonresonantbackground(x,&par[3]);
}
   
   TCanvas *c1 = new TCanvas("c1","#Kshort #arrow p#pi",10,10,1300,1300);
   c1->SetFillColor(33);
   c1->SetFrameFillColor(10); // white
   //c1->SetGrid();
   
//gStyle->SetOptFit(1011);

//gEnv->Print();

//...not obsolete
//gEnv->GetValue(“Canvas.SavePrecision”, -1);
//gEnv->SetValue(“Canvas.SavePrecision”, 16);

   hm4rec2OSpipitt05acc->SetMarkerStyle(21);
   hm4rec2OSpipitt05acc->SetMarkerSize(0.8);
// hm4rec2OSpipitt05acc->SetStats(0);

   Double_t BIN_SIZE = hm4rec2OSpipitt05acc->GetBinWidth(0);

   Double_t xmin = 0.280;
   //Double_t xmin = 0.20;
   Double_t xmax = 4.000;

   // create a TF1 with 10 parameters
   TF1 *fitlam = new TF1("fitlam",mylamfit,xmin,xmax,7);
   fitlam->SetNpx(2500);
   fitlam->SetLineWidth(2);
   fitlam->SetLineColor(kMagenta );

   // first try without starting values for the parameters
   // This defaults to 1 for each param.
   // this results in an ok fit for the polynomial function
   // however the non-linear part (lorenzian) does not
   // respond well.
   fitlam->SetParameters(1,1,1,1,1,1,1);
   hm4rec2OSpipitt05acc->Fit("fitlam","0");

  //...only 10 names allowed
 // fitlam->SetParNames(
//  "signal 1 p0","signal 1 p1","signal 1 p2",
//  "signal 2 p3","f_{0}(1500) #Gamma","f_{0}(1500) M",
//  "signal 3 p6","f_{0}(1710) #Gamma","f_{0}(1710) M",
//  "backgd p9");

// Robert's background
// a =  0.180  factor
// b =  1.476  cortex
// c =  9.369  peek
// d = -3.845  shape

   // second try: set start values for some parameter
   // signal
   fitlam->SetParameter(0,3477.);    // bw 1
   fitlam->SetParLimits(0,1.04,60000.);
   fitlam->SetParameter(1,0.0304733);   // gamma 1
   fitlam->SetParLimits(1,0.020,0.100);   // gamma 1
   fitlam->SetParameter(2,0.498907);   // mean 1 
   fitlam->SetParLimits(2,0.489,0.503);   // mean 1
   //
   fitlam->SetParameter(3,21460.);   // a factor
   //fitlam->SetParLimits(3,1.,15000.);
   fitlam->SetParameter(4,0.284);       // b cortex
   fitlam->SetParLimits(4,0.280,0.301); // b
   fitlam->SetParameter(5,0.300);        // c peek
   fitlam->SetParLimits(5,0.290,0.415);  // c
   fitlam->SetParameter(6,-3.0);    // d shape
   fitlam->SetParLimits(6,-4.0,-1.0);    // d
   //
   hm4rec2OSpipitt05acc->Fit("fitlam","V+","ep",xmin,xmax);
//   hm4rec2OSpipitt05acc->Fit("fitlam","R","ep",xmin,xmax);
// QWMRL

   // improve the picture:
   TF1 *signal1 = new TF1("signal1",mybw,xmin,xmax,3);
   signal1->SetLineWidth(2);
   signal1->SetLineColor(kBlue);
   signal1->SetNpx(25000);
   TF1 *backg = new TF1("backg",nonresonantbackground,xmin,xmax,4);
   backg->SetLineWidth(2);
   backg->SetLineColor(kRed);
   backg->SetNpx(25000);
   Double_t par[7];

// writes the fit results into the par array
   fitlam->GetParameters(par);

   signal1->SetParameters(par);
   signal1->Draw("same");
   backg->SetParameters(&par[3]);
   backg->Draw("same");

   // draw the legend
   TLegend *legend=new TLegend(0.485,0.283,0.779,0.450);
   legend->SetTextFont(42);
   legend->SetTextSize(0.02);
   legend->AddEntry(hm4rec2OSpipitt05acc,"data","lp");
   legend->AddEntry(signal1,"BW signal 1: 9488.7 pb","l");
   legend->AddEntry(fitlam,"global fit: 30093.8 pb","l");
   legend->AddEntry(backg,"background: 20605.1 pb","l");
   legend->Draw();

//(TObject*)0, "Some text", "");

//double sg1 = signal1->Integral(xmin,xmax);
double sg1 = signal1->Integral(xmin,xmax)/BIN_SIZE;
std::cout << " sg1 = " << sg1 << std::endl;
//double bkg = backg->Integral(xmin,xmax);
double bkg = backg->Integral(xmin,xmax)/BIN_SIZE;
std::cout << " bkg = " << bkg << std::endl;
double global = fitlam->Integral(xmin,xmax)/BIN_SIZE;
std::cout << " global = " << global << std::endl;









//double rrrr = fitlam->Integral(1.0,3.0);
//std::cout << " rrrr = " << rrrr << std::endl

//main->Fit("func","QWEMR"); // or with likelihood: "LQWEMR"

//Double_t factor = 1.;
//hm4rec2OS14lam->Scale(factor/hm4rec2OS14lam->Integral());

//BIN_SIZE = my_hist->GetBinWidth(0);
//fhist->Integral(LOW_LIMIT, fit_limit)/BIN_SIZE;