#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TVector3.h>
#include <TLorentzVector.h>
#include <TRotation.h>
#include <TLorentzRotation.h>

#include <cmath>

struct MyIntegFunc {
   MyIntegFunc(TF1 * f): fFunc(f) {}
   double Integral (double *x, double * ) const {
      double a = fFunc->GetXmin();
      return fFunc->Integral(a, *x);
   }
  TF1 * fFunc;
};

Double_t gauland(Double_t *x,Double_t *par) {

     Double_t risq2pi =  0.3989422804014;
     Int_t nstep = 10000;
     Double_t sc = 5.0;    
     Double_t xlow, xupp,step;
     Double_t sum;
     Double_t sigma,km,pga,norm;
 
     sum = 0; sigma=par[0]; km=par[1]; pga=par[2]; norm=par[3]; 
 
     xlow = x[0] - sc*sigma;
     xupp = x[0] + sc*sigma;
     step = (xupp-xlow)/(nstep*2.);

     for (Int_t ij=1;ij<=nstep;++ij){

      Double_t xxl = xlow + (ij-0.5)*step;
      Double_t xxc = xupp - (ij-0.5)*step;

      Double_t smear = exp(-(x[0]-xxl)*(x[0]-xxl)*0.5/(sigma*sigma))*risq2pi/sigma;  
      Double_t landa = pga*risq2pi*risq2pi/((xxl-km)*(xxl-km)+0.25*pga*pga);   
      sum += smear*landa;

      smear = exp(-(x[0]-xxc)*(x[0]-xxc)*0.5/(sigma*sigma))*risq2pi/sigma; 
      landa = pga*risq2pi*risq2pi/((xxc-km)*(xxc-km)+0.25*pga*pga);    
      sum += smear*landa;

     }

     Double_t fitval =  norm*step*sum;
     return fitval;
}


void piaga_ROOT() {
  
  TChain chain("h1");
  chain.Add("SelppM_new_0_9.root");
  chain.Add("SelppM_new_10_29.root");
  chain.Add("SelppM_new_30_50.root");

  TH1D *hma2 = new TH1D("hma2","M_invariant for 2-body decay", 200, 475., 512.);

  TCut cutcos="sqrt(pkvt[0][0]*pkvt[0][0]+pkvt[0][1]*pkvt[0][1])/sqrt(pkvt[0][0]*pkvt[0][0]+pkvt[0][1]*pkvt[0][1]+pkvt[0][2]*pkvt[0][2])>=0.9&&sqrt(ppvt[0][0]*ppvt[0][0]+ppvt[0][1]*ppvt[0][1])/sqrt(ppvt[0][0]*ppvt[0][0]+ppvt[0][1]*ppvt[0][1]+ppvt[0][2]*ppvt[0][2])>=0.9";    
  TCut spion="pstarpvt[0]<220.";
  TCut fkin = "sqrt(rkvt[0][0]*rkvt[0][0]+rkvt[0][1]*rkvt[0][1])>30"; 
  TCut cutsign = "chk[0]*cht[0]*flip[0]==1";
  TCut add = "sqrt(ppvt[0][0]*ppvt[0][0]+ppvt[0][1]*ppvt[0][1])>= 140.&&sqrt(ppvt[0][0]*ppvt[0][0]+ppvt[0][1]*ppvt[0][1])<=240.&&sqrt(pkvt[0][0]*pkvt[0][0]+pkvt[0][1]*pkvt[0][1])>80.&&Minvariant[0]<=520.&&Minvariant[0]>=470.";


  chain.Draw("Minvariant[0]>>hma2", cutcos&&spion&&fkin&&cutsign&&add&&"nt==1&&nkg==1");

 
  TF1 *fgauland = new TF1("gauland", gauland, 475., 512.,4);
  fgauland->SetParameters(2., 493.677, 1., 1000.);
  fgauland->SetParNames ("Sigma","Mean","Pga","Norm");
 
  
  TCanvas *c3 = new TCanvas("c3", "c3");
  
  hma2->Fit("gauland","RVE","", 483., 503.);
  hma2->Draw();
  gStyle->SetOptFit(111111);
  hma2->Draw("samehisterror");
  Double_t binwidth= 37./200.; 
  Double_t ierr;
  Double_t fcha = (483. - 475.)/binwidth + 1.;
  Double_t lcha = (503. - 475.)/binwidth + 1.; 
  Int_t ndof; 
  Double_t histen = hma2->IntegralAndError(fcha,lcha,ierr);
  
  Double_t check=fgauland->Integral(483., 503.); check=check/binwidth;
  cout<<"check function integral  "<<check<<"  "<<histen<<"  "<<ierr<<"\n";
  
  
}