double BTagEffService::GetEfficiency(BTagger const &bTagger, double pt, double eta,
  unsigned flavour) const
{
    // Find the appropriate efficiency histogram. Load it if needed
    auto histGroupIt = effHists.find(bTagger);
    TH2 *hist = nullptr;
    
    if (histGroupIt == effHists.end())
    {
        // Try to load histograms for the given b-tagger
        const_cast<BTagEffService *>(this)->LoadEfficiencies(bTagger);
        hist = effHists.at(bTagger).at(flavour).get();
    }
    else
        hist = histGroupIt->second.at(flavour).get();
    
    
    // Make sure the histogram exists
    if (not hist)
    {
        std::ostringstream message;
        message << "BTagEffService[\"" << GetName() << "\"]::GetEfficiency: " <<
          "Failed to find an efficiency histogram for b tagger " << bTagger.GetTextCode() <<
          ", efficiency label \"" << curEffLabel << "\", jet flavour " << flavour << ".";
        throw std::runtime_error(message.str());
    }

    
    // Return the efficiency
    return hist->GetBinContent(hist->FindFixBin(pt, eta));
}

void HcalElectronicsSelector::fill(const MyElectronicsId& id, double value) {
 
  TH2* h = m_hist;
  int htrChan=(id.fiber-1)*3+1+id.fiberChan;
  double fpga= (id.tb==1)?(id.Slot+0.5):(id.Slot);

  h->SetBinContent(h->FindBin(fpga*1.0,htrChan*1.0,0.0),value);
}

void QAvertex(const Char_t *fdata, const Char_t *fmc)
{
  style();

  TFile *fdtin = TFile::Open(fdata);
  TList *ldtin = (TList *)fdtin->Get("clist");
  TH2 *hdtin = (TH2 *)ldtin->FindObject("zv");
  TH1 *hdt = (TH1 *)ldtin->FindObject("zvNoSel");
  SetHistoStyle(hdt, 20, kRed+1);
  hdt->Scale(1. / hdt->Integral());

  TH1 *hdt0010 = hdtin->ProjectionX("hdt0010", 1, 4);
  SetHistoStyle(hdt0010, 20, kRed+1);
  hdt0010->Scale(1. / hdt0010->Integral());

  TH1 *hdt7080 = hdtin->ProjectionX("hdt7080", 11, 11);
  SetHistoStyle(hdt7080, 25, kAzure-3);
  hdt7080->Scale(1. / hdt7080->Integral());
  
  TFile *fmcin = TFile::Open(fmc);
  TList *lmcin = (TList *)fmcin->Get("clist");
  TH1 *hmc = (TH1 *)lmcin->FindObject("zvNoSel");
  SetHistoStyle(hmc, 25, kAzure-3);
  hmc->Scale(1. / hmc->Integral());

  TCanvas *c = new TCanvas("cVertex", "cVertex", 800, 800);
  TH1 * hfr = c->DrawFrame(-20., 0., 20., 0.1);
  hfr->SetTitle(";#it{z}_{vtx};");
  hdt0010->Draw("same");
  hdt7080->Draw("same");
  TLegend *legend = new TLegend(0.20, 0.18+0.60, 0.50, 0.30+0.60);
  legend->SetFillColor(0);
  legend->SetBorderSize(0);
  legend->SetTextFont(42);
  legend->SetTextSize(0.04);
  legend->AddEntry(hdt0010, "0-10%", "p");
  legend->AddEntry(hdt7080, "70-80%", "p");
  legend->Draw("same");
  c->SaveAs(canvasPrefix+"vertex.pdf");
  
  TCanvas *c1 = new TCanvas("cVertexDataMC", "cVertexDataMC", 800, 800);
  hfr = c1->DrawFrame(-20., 0., 20., 0.1);
  hfr->SetTitle(";#it{z}_{vtx};");
  hdt->Draw("same");
  hmc->Draw("same");
  legend = new TLegend(0.20, 0.18+0.60, 0.50, 0.30+0.60);
  legend->SetFillColor(0);
  legend->SetBorderSize(0);
  legend->SetTextFont(42);
  legend->SetTextSize(0.04);
  legend->AddEntry(hdt, "data", "p");
  legend->AddEntry(hmc, "Monte Carlo", "p");
  legend->Draw("same");
  c1->SaveAs(canvasPrefix+"vertexDataMC.pdf");
  
  //return 0;  
}

void plot_corr(TH2* phin, float cmin, float cmax) {
  ostringstream ssnam;
  ssnam << phin->GetName() << "_" << cmin << "_" << cmax;
  new TCanvas;
  TH2* ph = phin->Clone(ssnam.str().c_str());
  ph->GetXaxis()->SetRangeUser(cmin, cmax);
  ph->GetYaxis()->SetRangeUser(cmin, cmax);
  ph->Draw("colz");
  addaxis();
}

void QAoccupancy(const Char_t *fdata, const Char_t *fmc)
{
  style();

  TFile *fdtin = TFile::Open(fdata);
  TList *ldtin = (TList *)fdtin->Get("clist");
  TH2 *hdtin = (TH2 *)ldtin->FindObject("NClustersSPD2");
  TProfile *pdtin = hdtin->ProfileY("pdtin_clusters");
  pdtin->SetMarkerStyle(20);
  pdtin->SetMarkerSize(2);
  pdtin->SetMarkerColor(kAzure-3);
  
  TFile *fmcin = TFile::Open(fmc);
  TList *lmcin = (TList *)fmcin->Get("clist");
  TH2 *hmcin = (TH2 *)lmcin->FindObject("NClustersSPD2");
  TProfile *pmcin = hmcin->ProfileY("pmcin_clusters");
  pmcin->SetMarkerStyle(25);
  pmcin->SetMarkerSize(2);
  pmcin->SetMarkerColor(kRed+1);

  TCanvas *c = new TCanvas("cOccupancy", "cOccupancy", 800, 800);
  c->SetLogy();
  TH1 * hfr = c->DrawFrame(-0.5, 2., 10.5, 500.);
  DrawBinLabelsX(hfr, kTRUE);
  hfr->SetTitle(";;#LT#it{N}_{clusters,SPD-1}#GT");
  pdtin->DrawCopy("same");
  pmcin->DrawCopy("same");
  TLegend *legend = new TLegend(0.20, 0.18, 0.50, 0.30);
  legend->SetFillColor(0);
  legend->SetBorderSize(0);
  legend->SetTextFont(42);
  legend->SetTextSize(0.04);
  legend->AddEntry(pdtin, "data", "pl");
  legend->AddEntry(pmcin, "Monte Carlo", "pl");
  legend->Draw("same");
  c->SaveAs(canvasPrefix+"occupancy.pdf");
  return;
  TCanvas *cr = new TCanvas("cOccupancyr", "cOccupancyr", 800, 800);
  // hfr = cr->DrawFrame(-0.5, 0.75, 10.5, 1.25);
  // DrawBinLabelsX(hfr, kTRUE);
  // hfr->SetTitle(";;#LT#it{N}_{clusters,SPD-1}#GT ratio");
  pdtin->SetLineColor(kAzure-3);
  pdtin->SetLineWidth(3);
  pdtin->Divide(pmcin);
  pdtin->Draw("same,histo"); 
  legend = new TLegend(0.505025, 0.760673, 0.805276, 0.930142);
  legend->SetFillColor(0);
  legend->SetBorderSize(0);
  legend->SetTextFont(42);
  legend->SetTextSize(0.04);
  legend->AddEntry(pdtin, "data / Monte Carlo", "l");
  legend->Draw("same");
  cr->SaveAs(canvasPrefix+"occupancyr.pdf");
  
}

void DynamicExec()
{
   // Example of function called when a mouse event occurs in a pad.
   // When moving the mouse in the canvas, a second canvas shows the
   // projection along X of the bin corresponding to the Y position
   // of the mouse. The resulting histogram is fitted with a gaussian.
   // A "dynamic" line shows the current bin position in Y.
   // This more elaborated example can be used as a starting point
   // to develop more powerful interactive applications exploiting CINT
   // as a development engine.
   //
   // Author:  Rene Brun
   
   TObject *select = gPad->GetSelected();
   if(!select) return;
   if (!select->InheritsFrom("TH2")) {gPad->SetUniqueID(0); return;}
   TH2 *h = (TH2*)select;
   gPad->GetCanvas()->FeedbackMode(kTRUE);

   //erase old position and draw a line at current position
   int pyold = gPad->GetUniqueID();
   int px = gPad->GetEventX();
   int py = gPad->GetEventY();
   float uxmin = gPad->GetUxmin();
   float uxmax = gPad->GetUxmax();
   int pxmin = gPad->XtoAbsPixel(uxmin);
   int pxmax = gPad->XtoAbsPixel(uxmax);
   if(pyold) gVirtualX->DrawLine(pxmin,pyold,pxmax,pyold);
   gVirtualX->DrawLine(pxmin,py,pxmax,py);
   gPad->SetUniqueID(py);
   Float_t upy = gPad->AbsPixeltoY(py);
   Float_t y = gPad->PadtoY(upy);

   //create or set the new canvas c2
   TVirtualPad *padsav = gPad;
   TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
   if(c2) delete c2->GetPrimitive("Projection");
   else   c2 = new TCanvas("c2","Projection Canvas",710,10,700,500);
   c2->SetGrid();
   c2->cd();

   //draw slice corresponding to mouse position
   Int_t biny = h->GetYaxis()->FindBin(y);
   TH1D *hp = h->ProjectionX("",biny,biny);
   hp->SetFillColor(38);
   char title[80];
   sprintf(title,"Projection of biny=%d",biny);
   hp->SetName("Projection");
   hp->SetTitle(title);
   hp->Fit("gaus","ql");
   hp->GetFunction("gaus")->SetLineColor(kRed);
   hp->GetFunction("gaus")->SetLineWidth(6);
   c2->Update();
   padsav->cd();
}

TH2* RootWriter::CreateTH2(Histogram2D* h)
{
    const Axis& xax = h->GetAxisX(), yax = h->GetAxisY();
    const int xchannels = xax.GetBinCount();
    const int ychannels = yax.GetBinCount();
    TH2* mat = new TH2F( h->GetName().c_str(), h->GetTitle().c_str(),
                         xchannels, xax.GetLeft(), xax.GetRight(),
                         ychannels, yax.GetLeft(), yax.GetRight() );
    mat->SetOption( "colz" );
    mat->SetContour( 64 );

    TAxis* rxax = mat->GetXaxis();
    rxax->SetTitle(xax.GetTitle().c_str());
    rxax->SetTitleSize(0.03);
    rxax->SetLabelSize(0.03);

    TAxis* ryax = mat->GetYaxis();
    ryax->SetTitle(yax.GetTitle().c_str());
    ryax->SetTitleSize(0.03);
    ryax->SetLabelSize(0.03);
    ryax->SetTitleOffset(1.3);

    TAxis* zax = mat->GetZaxis();
    zax->SetLabelSize(0.025);

    for(int iy=0; iy<ychannels+2; ++iy)
        for(int ix=0; ix<xchannels+2; ++ix)
            mat->SetBinContent(ix, iy, h->GetBinContent(ix, iy));
    mat->SetEntries( h->GetEntries() );

    return mat;
}

//________________________________________________________________
void KVCanvas::DynamicZoom(Int_t Sign, Int_t px, Int_t py)
{
   // Zoom in or out of histogram with mouse wheel

   //    Info("DynamicZoom","px=%d py=%d",px,py);

   if (!fSelected) return;
   TH2* TheHisto = (TH2*) FindHisto();//fSelected;

   Double_t percent = 0.15 - Sign * 0.05;

   Int_t dX = 0;
   Int_t dY = 0;

   Double_t ppx = AbsPixeltoX(px);
   Double_t ppy = AbsPixeltoY(py);

   TAxis* ax = TheHisto->GetXaxis();
   Int_t NbinsXtmp = ax->GetNbins();
   Int_t X0tmp = ax->GetFirst();
   Int_t X1tmp = ax->GetLast();
   Int_t step = TMath::Min(TMath::Max(1, (Int_t)(percent * (X1tmp - X0tmp))), NbinsXtmp / 2);
   step *= Sign;
   X0tmp = TMath::Min(TMath::Max(X0tmp + step, 1), X1tmp - step);
   X1tmp = TMath::Max(TMath::Min(X1tmp - step, NbinsXtmp), X0tmp);
   if (X0tmp >= X1tmp) X0tmp = X1tmp - 1;
   if (Sign > 0) dX = (Int_t)(X0tmp + (X1tmp - X0tmp) * 0.5 - ax->FindBin(ppx));
   if ((X0tmp - dX) < 0) ax->SetRange(0, X1tmp - X0tmp);
   else if ((X1tmp - dX) > ax->GetNbins()) ax->SetRange(ax->GetNbins() - (X1tmp - X0tmp), ax->GetNbins());
   else ax->SetRange(X0tmp - dX, X1tmp - dX);

   ax = TheHisto->GetYaxis();
   Int_t NbinsYtmp = ax->GetNbins();
   Int_t Y0tmp = ax->GetFirst();
   Int_t Y1tmp = ax->GetLast();
   step = TMath::Min(TMath::Max(1, (Int_t)(percent * (Y1tmp - Y0tmp))), NbinsYtmp / 2);
   step *= Sign;
   Y0tmp = TMath::Min(TMath::Max(Y0tmp + step, 1), Y1tmp - step);
   Y1tmp = TMath::Max(TMath::Min(Y1tmp - step, NbinsYtmp), Y0tmp);
   if (Y0tmp >= Y1tmp) Y0tmp = Y1tmp - 1;
   if (Sign > 0) dY = (Int_t)(Y0tmp + (Y1tmp - Y0tmp) * 0.5 - ax->FindBin(ppy));
   if ((Y0tmp - dY) < 0) ax->SetRange(0, Y1tmp - Y0tmp);
   else if ((Y1tmp - dY) > ax->GetNbins()) ax->SetRange(ax->GetNbins() - (Y1tmp - Y0tmp), ax->GetNbins());
   else ax->SetRange(Y0tmp - dY, Y1tmp - dY);

   Modified();
   Update();
   return;
}

void
HistoCompare::PVCompute(TH2 * oldHisto , TH2 * newHisto , TText * te )
{

  myoldHisto2 = oldHisto;
  mynewHisto2 = newHisto;
  myte = te;

  Double_t *res ;
  Double_t mypv = myoldHisto2->Chi2Test(mynewHisto2,"WW",res);
  TString title = myoldHisto2->GetName();
  printRes(title, mypv, myte);
  return;

}

void exec2()
{

   if (!gPad) {
      Error("exec2", "gPad is null, you are not supposed to run this macro");
      return;
   }


   TObject *select = gPad->GetSelected();
   if(!select) return;
   if (!select->InheritsFrom(TH2::Class())) {gPad->SetUniqueID(0); return;}
   gPad->GetCanvas()->FeedbackMode(kTRUE);

   //erase old position and draw a line at current position
   int pyold = gPad->GetUniqueID();
   int px = gPad->GetEventX();
   int py = gPad->GetEventY();
   float uxmin = gPad->GetUxmin();
   float uxmax = gPad->GetUxmax();
   int pxmin = gPad->XtoAbsPixel(uxmin);
   int pxmax = gPad->XtoAbsPixel(uxmax);
   if(pyold) gVirtualX->DrawLine(pxmin,pyold,pxmax,pyold);
   gVirtualX->DrawLine(pxmin,py,pxmax,py);
   gPad->SetUniqueID(py);
   Float_t upy = gPad->AbsPixeltoY(py);
   Float_t y = gPad->PadtoY(upy);

   //create or set the new canvas c2
   TVirtualPad *padsav = gPad;
   TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
   if(c2) delete c2->GetPrimitive("Projection");
   else   c2 = new TCanvas("c2");
   c2->cd();

   //draw slice corresponding to mouse position
   TH2 *h = (TH2*)select;
   Int_t biny = h->GetYaxis()->FindBin(y);
   TH1D *hp = h->ProjectionX("",biny,biny);
   char title[80];
   sprintf(title,"Projection of biny=%d",biny);
   hp->SetName("Projection");
   hp->SetTitle(title);
   hp->Fit("gaus","ql");
   c2->Update();
   padsav->cd();
}

void
polar()
{
   TCanvas *c1 = new TCanvas("c1","c1",500,500);
   TGraphPolar * grP1 = new TGraphPolar();
   grP1->SetTitle("TGraphPolar example");

   for (Int_t i = 0; i < 2; i++) {
     grP1->SetPoint(i, (i+1) * TMath::Pi() / 4, (i+1) * 0.05);
     grP1->SetPointError(i, 9*TMath::Pi()/180, 0.0);
     Double_t rr = grP1->GetY()[i];
     Double_t tt = grP1->GetX()[i];
     Double_t x  = rr * TMath::Cos(tt);
     Double_t y  = rr * TMath::Sin(tt);
     Printf("(x,y)=(%f,%f)", x, y);
   }
   Double_t r = 1;
   TH2* frame = new TH2F("frame", "Frame", 100, -r,r, 100, -r, r);
   frame->Draw();

   grP1->SetMarkerStyle(1);
   grP1->SetMarkerSize(1.);
   grP1->SetMarkerColor(4);
   grP1->SetLineColor(4);
   grP1->SetLineWidth(3);
   grP1->SetFillColor(kRed+1);
   grP1->SetFillStyle(3001);
   grP1->Draw("PNEF same");
   // grP1->Draw("APNEF");

   // Update, otherwise GetPolargram returns 0
   c1->Update();
   TGraphPolargram* gram = grP1->GetPolargram();
   gram->SetLineWidth(0);
   // gram->SetLineColor(kWhite);
   gram->SetNdivPolar(20);
   gram->SetNdivRadial(10);   
   gram->SetTextSize(0);
   gram->SetRadialLabelSize(0);
   gram->SetPolarLabelSize(0);
   gram->SetAxisAngle(9*TMath::Pi()/180);
   gram->SetTwoPi();
   gram->SetToRadian();
   c1->SetGridx();
   c1->SetGridy();
   c1->Update();
}

//______________________________________________________________________________
TH2* CreateHisto(const Char_t* name, Int_t nPar, const Double_t* par)
{
    TH2* h;
    if (gIsTAPS)
        h = new TH2F(name, name, 4000, -100, 100, nPar, 0, nPar);
    else
        h = new TH2F(name, name, 2000, -100, 100, nPar, 0, nPar);

    // fill histo
    for (Long64_t i = 0; i < gEntries; i++)
    {
        const Double_t time_1 = gGain[gElem_1[i]] * (gRaw_Time_1[i] - par[gElem_1[i]]);
        const Double_t time_2 = gGain[gElem_2[i]] * (gRaw_Time_2[i] - par[gElem_2[i]]);
        h->Fill(time_1 - time_2, gElem_1[i]);
        h->Fill(time_2 - time_1, gElem_2[i]);
    }

    return h;
}

void root_hist_test()
{
  // Project3D uses overflows?
  TH3* h3D = new TH3D("h3D", "h3D", 2,0,2, 2,0,2, 2,0,2);
  for(int i=0; i<=h3D->GetNbinsX()+1; i++) {
    for(int j=0; j<=h3D->GetNbinsY()+1; j++) {
      for(int k=0; k<=h3D->GetNbinsZ()+1; k++) {
        h3D->SetBinContent(i,j,k, 1);
      }
    }
  }

  printf("h3D integral = %d\tentries = %d\n", h3D->Integral(), h3D->GetEntries());

  h3D->GetZaxis()->SetRange(3,3);
  TH2* hYX = (TH2*) h3D->Project3D("yx");

  printf("hYX integral = %d\tentries = %d\n", hYX->Integral(), hYX->GetEntries());
  hYX->Draw("colz");
}

  void Draw() 
  {
    float min = 100000;
    float max = 0;
    std::vector<TGraph*> gs;
    float_array bg(fTKine.size());
    for (size_t i = 0; i < fTKine.size(); i++) bg[i] = fTKine[i]/fPDG->Mass();
    for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
      if (!(*j)->fStatus) continue;
      for (size_t i = 0; i < fTKine.size(); i++) {
	if ((*j)->fValues[i] == 0) continue;
	min = std::min(min, (*j)->fValues[i]);
	max = std::max(max, (*j)->fValues[i]);
      }
    }

    TCanvas* c = new TCanvas("c", "C", 700, 700);
    c->SetFillColor(0);
    c->SetLogy();
    c->SetLogx();
    c->SetGridy();
    c->SetGridx();
    float_array y(101);
    float ymin = log10(min);
    float dy   = (log10(max)+.5 - log10(min)) / y.size();
    for (size_t i = 1; i < y.size(); i++)  y[i] = pow(10, ymin + i * dy);
    TH2* f = new TH2F("x", "X-sec",bg.size()-1,&(bg[0]),y.size()-1,&(y[0]));
    f->SetXTitle("#beta#gamma");
    f->SetDirectory(0);
    f->SetStats(kFALSE);
    f->Draw();
    TLegend* l = new TLegend(0.45, 0.125, 0.90, 0.45);
    l->SetFillColor(0);
    // l->SetFillStyle(0);
    for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
      if (!(*j)->fStatus) continue;
      TGraph& g = (*j)->Draw();
      l->AddEntry(&g, g.GetTitle(), "l");
    }
    l->Draw("same");
  }

void smoothHisto(TH2* hist, int Niter) {

  for(int iter=0; iter<Niter; iter++) {
    if(iter>2) {
      TH2* clone = (TH2*) hist->Clone("clone");
      for(int ix=2; ix<hist->GetNbinsX()-1; ix++) {
	for(int iy=2; iy<hist->GetNbinsY()-1; iy++) {
	  double nn1 = clone->GetBinContent(ix-1,iy);
	  double nn2 = clone->GetBinContent(ix,iy-1);
	  double nn3 = clone->GetBinContent(ix+1,iy);
	  double nn4 = clone->GetBinContent(ix,iy+1);
	  hist->SetBinContent(ix,iy,0.25*(nn1+nn2+nn3+nn4));
	}
      }
      delete clone;
    } else {
      TH2* clone = (TH2*) hist->Clone("clone");
      for(int ix=2; ix<hist->GetNbinsX()-1; ix++) {
	for(int iy=2; iy<hist->GetNbinsY()-1; iy++) {
	  double nn1 = clone->GetBinContent(ix-1,iy-1);
	  double nn2 = clone->GetBinContent(ix-1,iy+1);
	  double nn3 = clone->GetBinContent(ix+1,iy-1);
	  double nn4 = clone->GetBinContent(ix+1,iy+1);
	  if(nn1+nn2+nn3+nn4>2 && hist->GetBinContent(ix,iy))
	    hist->SetBinContent(ix,iy,1.);
	  else
	    hist->SetBinContent(ix,iy,0.25*(nn1+nn2+nn3+nn4));
	}
      }
      delete clone;      
    }
  }
}

void DrawHistos(TString infile = "histos.root"){
    set_plot_style();
    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(0);
    
    TFile *_file0 = TFile::Open(infile);
    TH1* hHardPartonspT = (TH1*)gROOT->FindObject("hHardPartonspT");
    TH1* hHardPartonsName = (TH1*)gROOT->FindObject("hHardPartonsName");
    TH1* hParticlepT = (TH1*)gROOT->FindObject("hParticlepT");
    TH2* hParticleEtaPhi = (TH2*)gROOT->FindObject("hParticleEtaPhi");

    TCanvas* can1 = new TCanvas("can1","",600,600);
    gPad->SetLeftMargin(0.15);
    gPad->SetLogy();
    SetHist(hParticlepT,"p_{T} (GeV/c)","Counts",2);
    SetHist(hHardPartonspT,"p_{T} (GeV/c)","Counts",1);
    hParticlepT->Draw();
    hHardPartonspT->Draw("same");
    TLegend* leg1 = new TLegend(0.6,0.7,0.85,0.85);
    SetLeg(leg1);
    leg1->AddEntry(hHardPartonspT,"Hard Parton","L");
    leg1->AddEntry(hParticlepT,"Final Particles","L");
    leg1->Draw();
    can1->SaveAs("pTgraph.png");
    
    TCanvas* can2 = new TCanvas("can2","",600,600);
    gPad->SetLeftMargin(0.15);
    SetHist(hHardPartonsName,"parton ID","Counts",1);
    hHardPartonsName->Draw();
    
    TCanvas* can3 = new TCanvas("can3","",600,600);
    gPad->SetLeftMargin(0.15);
    SetHist(hParticleEtaPhi,"p_{T} (GeV/c)","Counts",1);
    hParticleEtaPhi->Draw("COLZ");
   
    
  }

TH2* smoothHoles(const TH2* originalHist)
{
  TH2* hist = (TH2*)originalHist->Clone("_smoothed");
  int xMax = hist->GetNbinsX();
  int yMax = hist->GetNbinsY();

  int xMin = 0;
  
  for(int xBin = 1; xBin <= xMax; xBin++)
    {
      for(int yBin = 1; yBin <= yMax; yBin++)
	{
	  if(hist->GetBinContent(xBin,yBin)>0)
	    {
	      xMin = xBin;
	      yBin = yMax+1;
	      xBin = xMax+1;
	    }
	}
    } 
  //  for(unsigned int i = 0; i< 1000; i++) smoothHistAcross(hist,xMin);
  for(unsigned int i = 0; i< 1000; i++) interpolateHistAcross(hist,xMin);
  return hist;
}

TList* getContours(const TH2* hist,double contourLevel,const TString filename)
{
  TH2* h = (TH2*)hist->Clone("_clone");
  double limitValue[1] = {contourLevel};
  h->SetContour(1,limitValue);
  TCanvas* c = new TCanvas("contour_canvas","Contour List",0,0,600,600);
  h->Draw("CONT LIST");
  c->Update();
  TList* contours = (TList*)((TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours"))->At(0);
  TGraph* contour = (TGraph*)contours->First();
  if(filename!="") 
    {
      for(unsigned int j = 0; j < contours->GetSize(); j++)
	{
	  TString newFilename = filename+"_";
	  newFilename+=j;
	  contour->SaveAs(newFilename+".C");
	  contour = (TGraph*)contours->After(contour); // Get Next graph
	}
    }
  delete h;
  delete c;
  return contours;
}

void
HistoCompare::PVCompute(TH2 * oldHisto , TH2 * newHisto , TText * te )
{

  myoldHisto2 = oldHisto;
  mynewHisto2 = newHisto;
  myte = te;

  Double_t *res ;
  Double_t mypvchi = myoldHisto2->Chi2Test(mynewHisto2,"WW",res);
  char str [128];
  sprintf(str,"chi^2 P Value: %f",mypvchi);
  TString title = str;
  printRes(title, mypvchi, myte);

  return;

}

void SetBorders( TH2 &hist, Double_t val=0 )
{
    int numx = hist.GetNbinsX();
    int numy = hist.GetNbinsY();

    for(int i=0; i <= numx+1 ; i++){
        hist.SetBinContent(i,0,val);
        hist.SetBinContent(i,numy+1,val);
    }
    for(int i=0; i <= numy+1 ; i++) {
        hist.SetBinContent(0,i,val);
        hist.SetBinContent(numx+1,i,val);
    }
}