#include <math.h>

// x-axis range for resolution plots
const Double_t xmin = -0.01;
const Double_t xmax =  0.01;
const Int_t nsigma  = 3;

using namespace RooFit;

void batchAnalysis_rf4 ( Int_t runNumber, Int_t detType, bool log = 0, bool fix = 0, Double_t SigmaDVal = 1.3e-3 ) {

  if ( detType != 0 && detType !=1 ) {
    cout << "Side can be tele ( <=> 0 ) or striplet ( <=> 1 )\n";
    cout << "detType : " << detType << "\n";
    exit( 1 );
  }

  if ( !TClass::GetDict("DoubleGauss") ) {
    gROOT->ProcessLine(".L DoubleGauss.C++");
  }



  RooRealVar u("u","Residuals on p side",xmin,xmax,"cm");
  RooRealVar v("v","Residuals on n side",xmin,xmax,"cm");

  RooRealVar uMean("uMean","U mean",0.,-2.e-4,2.e-4,"cm");
  RooRealVar uSigma("uSigma","U sigma",1.3e-3,+1.e-4,+2.5e-3,"cm");

  RooRealVar vMean("vMean","V mean",0.,-2.e-4,2.e-4,"cm");
  RooRealVar vSigma("vSigma","V sigma",1.3e-3,+1.e-4,+2.5e-3,"cm");

  RooGaussian *uPdf = new RooGaussian("uPdf","U pdf",u,uMean,uSigma);
  RooGaussian *vPdf = new RooGaussian("vPdf","V pdf",v,vMean,vSigma);

  RooRealVar uMeanD("uMeanD","U double mean",-2.5e-3,-3.e-3,-1.e-3,"cm");
  RooRealVar uSigmaD("uSigmaD","U double sigma",1.3e-3,1.e-3,2.e-3,"cm");

  RooRealVar vMeanD("vMeanD","V double mean",-2.5e-3,-5.e-3,-1.e-3,"cm");
  RooRealVar vSigmaD("vSigmaD","V double sigma",1.3e-3,.5e-3,2.5e-3,"cm");

  RooRealVar nUsig("nUsig","# of signal events U",1000.,0.,1000000.);
  RooRealVar nUd("nUd","# of backgd events U",200.,0.,100000.);

  RooRealVar nVsig("nVsig","# of signal events V",1000.,0.,100000.);
  RooRealVar nVd("nVd","# of backgd events V",200.,0.,10000.);

  DoubleGauss* doubleGaussU = new DoubleGauss("doubleGaussU","U double gauss pdf",u,uMeanD,uSigmaD);
  DoubleGauss* doubleGaussV = new DoubleGauss("doubleGaussV","V double gauss pdf",v,vMeanD,vSigmaD);

  
  RooAddPdf uFcn("uFcn","total residual U pdf",RooArgList(*uPdf,*doubleGaussU),RooArgList(nUsig,nUd));
  RooAddPdf vFcn("vFcn","total residual V pdf",RooArgList(*vPdf,*doubleGaussV),RooArgList(nVsig,nVd));

  if ( fix ) {
    uMeanD.setConstant();
    uSigmaD.setVal(SigmaDVal);
    uSigmaD.setConstant();


    vMeanD.setConstant();
    vSigmaD.setVal(SigmaDVal);
    vSigmaD.setConstant();
  }

  TString fileName ="";
  fileName.Form("rootuple/Sbt_Run%d_out.root",runNumber);
  
  TFile file(fileName.Data(),"read");
  file.cd();


  std::cout << "File open: " << fileName.Data() << "\n";



  TString saveNamepng;

  // Sytle and Canvas (to plot and save)
  gStyle->SetOptStat();
  gStyle->SetOptFit(1111);


  const Int_t nSlices = 20;

  TH1D* ResU_Vs_PH[nSlices];
  TH1D* ResV_Vs_PH[nSlices];
  
  TH1D* ResU_Vs_PH_CS1[nSlices];
  TH1D* ResV_Vs_PH_CS1[nSlices];
  
  TH1D* ResU_Vs_PH_CS2[nSlices];
  TH1D* ResV_Vs_PH_CS2[nSlices];
  
  TPaveLabel* t1[nSlices];
  
  TCanvas *c1 = new TCanvas("c1","",1024,768);
  c1->cd();
  TString hName;

  std::ofstream outlog;
  outlog.open("log/cazzivari.txt");


  saveNamepng.Form("plots/ResU_Vs_PH_CS1_Run%d.ps[",runNumber);
  c1->SaveAs(saveNamepng.Data());

  for ( int i = 0; i < 20; i++ ) {

    hName.Form("_resU_Vs_PH_%i_CS1_Det4",i);
    ResU_Vs_PH_CS1[i] = (TH1D*)gROOT->FindObject(hName.Data());
    if ( ResU_Vs_PH_CS1[i]->GetEntries() < 1000. ) {
      if ( i == 19 ) continue;
      hName.Form("_resU_Vs_PH_%i_CS1_Det4",++i);
      ResU_Vs_PH_CS1[--i]->Add((TH1D*)gROOT->FindObject(hName.Data()) );
      ResU_Vs_PH_CS1[i+1] = ResU_Vs_PH_CS1[i];
      i++;
      if ( ResU_Vs_PH_CS1[i]->GetEntries() < 1000. ) continue;
    }
    ResU_Vs_PH_CS1[i]->GetXaxis()->SetRangeUser(xmin,xmax);
    ResU_Vs_PH_CS1[i]->Draw("e");

    hName.Form("Residual U for CS = 1, bin %i",i);
    RooDataHist uData_CS1("uData_CS1",hName.Data(),u,ResU_Vs_PH_CS1[i]);
    RooPlot* uFrame = u.frame();
    uData_CS1.plotOn(uFrame);

    RooFitResult* r_ml_wgt = uFcn.fitTo(uData_CS1,Extended());
    uFcn.plotOn(uFrame);
    uFcn.paramOn(uFrame,Format(("NEU"),AutoPrecision(1)),Layout(.6,.99,.99));
    c1->SetLogy(log);

    hName.Form("_resU_Vs_PH_%i_CS1_Det4",i);
    Double_t signalchi = uFrame->chiSquare("uPdf",hName.Data(),6);
    outlog << "signal chi2 for bin " << i << " = " << signalchi << "\n";
    t1[i] = new TPaveLabel(0.1,0.1,0.3,0.28, Form("#chi^{2} = %f", signalchi),"brNDC");
    t1[i]->Draw();
    uFrame->addObject(t1[i]); 


    uFrame->Draw();
    saveNamepng.Form("plots/ResU_Vs_PH_CS1_Run%d.ps",runNumber);
    c1->SaveAs(saveNamepng.Data());

  }
  saveNamepng.Form("plots/ResU_Vs_PH_CS1_Run%d.ps]",runNumber);
  c1->SaveAs(saveNamepng.Data());
   
  // closing file
  outlog.close();
}