NuContour.cxx

Go to the documentation of this file.

#include <iostream>

#include "TF1.h"
#include "TH2.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TAxis.h"
#include "TGaxis.h"
#include "TLine.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TStyle.h"
#include "MessageService/MsgService.h"


#include "NtupleUtils/NuContour.h"

using namespace std;
int NuContour::nFillCols = 10;
int NuContour::usedFills =  0;
Int_t* NuContour::colFill =  0;
CVSID("$Id: ");


//____________________________________________________________________72
NuContour::NuContour() 
  : lin(true), loglin(false), log(false), 
    fillLevel(-1), contLevel(0), neg(1),
    hdivide(9999), ymid(0), bmid(250), bmax(500)
{
  if (!colFill) {
    colFill = new Int_t[10];
    for (int i = 0; i < 10; i++) {
      colFill[i] = 0;
    }
  }
  
  hist = 0;
  disp_hist = 0;
  hist_high = 0;
  
  graph = 0;
  disp_graph = 0;
}


//____________________________________________________________________72
NuContour::NuContour(const TH2D *hin) 
  : lin(true), loglin(false), log(false),
    fillLevel(-1), contLevel(0),  neg(1),
    hdivide(9999), ymid(0), bmid(250), bmax(500)
{
  if (!colFill) {
    colFill = new Int_t[10];
    for (int i = 0; i < 10; i++) {
      colFill[i] = 0;
    }
  }
  
  hist = new TH2D(*hin);
  disp_hist = 0;
  hist_high = 0;
  
  graph = 0;
  disp_graph = 0;
}


//____________________________________________________________________72
NuContour::NuContour(const TGraph *g) 
  : lin(true), loglin(false), log(false),
    fillLevel(-1), contLevel(0),  neg(1),
    hdivide(9999), ymid(0), bmid(250), bmax(500)
{
  
  if (!colFill) {
    colFill = new Int_t[10];
    for (int i = 0; i < 10; i++) {
      colFill[i] = 0;
    }
  }
  
  hist = 0;
  disp_hist = 0;
  hist_high = 0;
  
  graph = new TGraph(*g);
  disp_graph = 0;
}


//____________________________________________________________________72
NuContour::NuContour(const NuContour &rhs) 
{
  Reset();
  if (rhs.hist) hist = new TH2D(*rhs.hist);
  if (rhs.disp_hist) disp_hist =  new TH2D(*rhs.disp_hist);
  if (rhs.hist_high) hist_high =  new TH2D(*rhs.hist_high);
  if (rhs.graph) graph = new TGraph(*rhs.graph);
  if (rhs.disp_graph) disp_graph = new TGraph(*rhs.disp_graph);

  lin = rhs.lin;
  loglin = rhs.loglin;
  log = rhs.log;
  fillLevel = rhs.fillLevel;
  contLevel = rhs.contLevel;
  neg = rhs.neg;
    hdivide = rhs.hdivide;
  ymid = rhs.ymid;
  ymax = rhs.ymax;
  bmid = rhs.bmid;
  bmax = rhs.bmax;
}


//____________________________________________________________________72
NuContour & NuContour::operator=(const NuContour &rhs)
{
  // Check for self assignment
  if (this == &rhs) return *this;

  Reset();
  if (rhs.hist) hist = new TH2D(*rhs.hist);
  if (rhs.disp_hist) disp_hist =  new TH2D(*rhs.disp_hist);
  if (rhs.hist_high) hist_high =  new TH2D(*rhs.hist_high);
  if (rhs.graph) graph = new TGraph(*rhs.graph);
  if (rhs.disp_graph) disp_graph = new TGraph(*rhs.disp_graph);
  
  lin = rhs.lin;
  loglin = rhs.loglin;
  log = rhs.log;
  fillLevel = rhs.fillLevel;
  contLevel = rhs.contLevel;
  neg = rhs.neg;
  hdivide = rhs.hdivide;
  ymid = rhs.ymid;
  ymax = rhs.ymax;
  bmid = rhs.bmid;
  bmax = rhs.bmax;
  
  return *this;
}


//____________________________________________________________________72
NuContour::~NuContour()
{
    Reset();
}


//____________________________________________________________________72
void NuContour::Reset() 
{
    /*
    if (hist) delete hist;
    if (disp_hist) delete disp_hist;
    if (hist_high) delete hist_high;
    if (graph) delete graph;
    if (disp_graph) delete disp_graph;
    */
    
    hist = 0;
    disp_hist = 0;
    hist_high = 0;
    graph = 0;
    disp_graph = 0; 
    
}

//____________________________________________________________________72
void NuContour::Reset(TH2D *hin)
{
    Reset();
    hist = new TH2D(*hin);
    if (IsLogLin()) DoLogLin();
    if (IsLog())    DoLog();
}

//____________________________________________________________________72
void NuContour::Reset(TGraph *g)
{
    Reset();
    graph = new TGraph(*g);
    if (IsLogLin()) DoLogLin();
    if (IsLog())    DoLog();
}

//____________________________________________________________________72
void NuContour::AddExtraHist(TH2D *hin, double crossover)
{
    hist_high = new TH2D(*hin);
    hdivide = crossover;

    Double_t diff = MinBin(*hist) - MinBin(*hist_high);
    if (diff < 0) ShiftHist(hist_high, diff);
    else          ShiftHist(hist,     -diff);

    if (IsLogLin()) DoLogLin();
    if (IsLog())    DoLog();
}


//____________________________________________________________________72
void NuContour::Combine(TH2D* correction)
{
    double min = 0;

    if (hist) {
        min = MinBin(*hist);
        ShiftHist(hist, -min);
        Add(hist, correction);
    }
    if (hist_high) {
        ShiftHist(hist_high, -min);
        Add(hist_high, correction);
    }
    if (disp_hist) {
        ShiftHist(disp_hist, -min);
        Add(disp_hist, correction, 1, true);
    }
    
}


//____________________________________________________________________72
void NuContour::Combine(TH2D* correction_low, TH2D* correction_high)
{
    double min = 0;
    
    if (hist) {
        min = MinBin(*hist);
        ShiftHist(hist, -min);
        Add(hist, correction_low);
    }
    if (hist_high) {
        ShiftHist(hist_high, -min);
        Add(hist_high, correction_high);
        AddExtraHist(hist_high, hdivide);
    }
}


//____________________________________________________________________72
void NuContour::SetGauss(double level)
{
    double min = 0;
    contLevel = 0;
    
    if (hist) {
        min = MinBin(*hist);
        ShiftHist(hist, -min-level);
        ChooseContour(hist);
    }
    if (hist_high) {
        ShiftHist(hist_high, -min-level);
        ChooseContour(hist_high);
    }
    if (disp_hist) {
        ShiftHist(disp_hist, -min-level);
        ChooseContour(disp_hist);
    }
}




//____________________________________________________________________72
void NuContour::SetFC(TH2D* correction)
{
   double min = 0;
   contLevel = 0;

   if (hist) {
        min = MinBin(*hist);
        ShiftHist(hist, -min);
        Add(hist, correction, -1);
        ChooseContour(hist);
    }
    if (hist_high) {
        ShiftHist(hist_high, -min);
        Add(hist_high, correction, -1);
        ChooseContour(hist_high);
    }
    if (disp_hist) {
        ShiftHist(disp_hist, -min);
        Add(disp_hist, correction, -1, true);
        ChooseContour(disp_hist);
    }
    
}


//____________________________________________________________________72

void NuContour::Add(TH2D *surface, TH2D* correction, double scale, bool convert)
{
  double xval, yval, val, corr;

  for (int xbin = 1; xbin <= surface->GetNbinsX(); xbin++) {
    for (int ybin = 1; ybin <= surface->GetNbinsY(); ybin++) {
      val = surface->GetBinContent(xbin, ybin);

      xval = surface->GetXaxis()->GetBinCenter(xbin);
      yval = surface->GetYaxis()->GetBinCenter(ybin);
      if (convert) yval = Can2Cont(yval);
      
      corr = interpolate(correction, xval, yval, 0);

      surface->SetBinContent(xbin, ybin, val + scale*corr);
    }
  }
}


//____________________________________________________________________72
void NuContour::ContourStyle(int color, int width, int style)
{
    if (graph) {
        graph->SetLineWidth(width);
        graph->SetLineColor(color);
        graph->SetLineStyle(style);
    }
    if (disp_graph) {
        disp_graph->SetLineWidth(width);
        disp_graph->SetLineColor(color);
        disp_graph->SetLineStyle(style);
    }
    
    if (hist) {
        if (style != 1 && !graph && !disp_graph) {
            MSG("NuContour",Msg::kWarning) 
            << "Line styles other than solid do not come out well on histogram." << endl
            << "For better results call ConvertToGraph() on this contour." << endl;
        }
        
        hist->SetLineWidth(width);
        hist->SetLineColor(color);
        hist->SetLineStyle(style);
        ChooseContour(hist);
    }
    if (hist_high) {
        hist_high->SetLineWidth(width);
        hist_high->SetLineColor(color);
        hist_high->SetLineStyle(style);
        ChooseContour(hist_high);
    }
    if (disp_hist) {
        disp_hist->SetLineWidth(width);
        disp_hist->SetLineColor(color);
        disp_hist->SetLineStyle(style);
        ChooseContour(disp_hist);
    }
}

//____________________________________________________________________72
void NuContour::ContourFillStyle(int color, int style)
{
 
    if (graph || disp_graph) {
        fillLevel = 999;
        if (graph) {
            graph->SetFillColor(color);
            graph->SetFillStyle(style);
        }
        if (disp_graph) {
            disp_graph->SetFillColor(color);
            disp_graph->SetFillStyle(style);
        }
    }
    else {
        neg = -1;
        fillLevel = usedFills;
        usedFills++;
        
        MSG("NuContour",Msg::kInfo) << "Setting countour level " << fillLevel <<  " to color " << color << endl;
        MSG("NuContour",Msg::kInfo) << "Current palette: " << endl;
        colFill[nFillCols-fillLevel-1] = color;
        gStyle->SetPalette(nFillCols,colFill);
        for (int i = nFillCols - fillLevel - 1; i < nFillCols; i++) {
            MSG("NuContour",Msg::kInfo) << "  colFill[" << i << "] = " << colFill[i] << endl;
        }
        
        
        if (hist) {
            ChooseContour(hist);
            hist->SetFillColor(color);
            hist->SetFillStyle(style);
        }
        if (hist_high) {
            ChooseContour(hist_high);
            hist_high->SetFillColor(color);
            hist_high->SetFillStyle(style);
        }
        if (disp_hist) {
            ChooseContour(disp_hist);
            disp_hist->SetFillColor(color);
            disp_hist->SetFillStyle(style);
        }
    }
    
  //ContourStyle(lc, width, style);
}

//____________________________________________________________________72
void NuContour::DoLin()
{
  lin = true;
  loglin = false;
  log = false;
  bmid = bmax;

  hist->GetXaxis()->CenterTitle();
  TAxis *a1 = hist->GetYaxis();
  a1->SetLabelSize(0);
  a1->SetTickLength(0);
  a1->SetTitleOffset(1.2);
  a1->CenterTitle();
  a1->SetTitleOffset(1.2);
  hist->SetTitle(";sin^{2}(2#bar{#theta});|#Delta#bar{m}^{2}| (10^{-3} eV^{2})");
}


//____________________________________________________________________72
void NuContour::DoLogLin()
{
  lin = false;
  loglin = true;
  log = false;
  bmid = bmax / 2;
  GenerateLogLin();
}


//____________________________________________________________________72
void NuContour::DoLog()
{
  lin = false;
  loglin = false;
  log = true;
  bmid = 0;
  GenerateLog();
}


//____________________________________________________________________72
TH2D * NuContour::GetHist()
{
    if (!hist) {
      MSG("NuContour",Msg::kWarning)
        << "Histogram not defined for this contour." << endl;
        return 0;
    }
    if (IsLin())         {
        //cout << "GetHist returning hist" << endl;        
        return hist;
    }
    else if (IsLogLin()) {
        //cout << "GetHist returning disp_hist (loglin)" << endl;       
        return disp_hist;
    }
    else if (IsLog()) {
        //cout << "GetHist returning disp_hist (log)" << endl;       
        return disp_hist;
    }
    else {
        MSG("NuContour",Msg::kError)
        << "Contour is in undefined mode.  This shouldn't happen." << endl;
        return 0;
    }
}


//____________________________________________________________________72
TGraph * NuContour::GetGraph()
{
    if (!graph && !disp_graph) {
        MSG("NuContour",Msg::kWarning)
        << "Graph not defined for this contour." << endl;
        return 0;
    }
    if (IsLin()) {
        //cout << "GetGraph returning graph" << endl;
        return graph;
    }
    else if (IsLogLin() || IsLog()) {
        //cout << "GetGraph returning disp_graph" << endl;
        return disp_graph;
    }
    else {
        MSG("NuContour",Msg::kError)
        << "Contour is in undefined mode.  This shouldn't happen." << endl;
        return 0;
    }
}



//____________________________________________________________________72
void NuContour::ConvertToGraph()
{
    if (graph && !hist) {
        MSG("NuContour",Msg::kWarning)
        << "Already a graph contour -- cannot be converted" << endl;
        return;
    }
    else if (!hist && !disp_hist) {
        MSG("NuContour",Msg::kWarning)
        << "No histogram -- cannot be converted" << endl;
        return;
    } 
    else if (IsLogLin()) {
        MSG("NuContour",Msg::kWarning)
        << "Convert to graph not supported for LogLin histograms" << endl;
        return;
    }
    else {
        MSG("NuContour",Msg::kInfo)
        << "Converting histogram to graph" << endl;
    }
    
    TVirtualPad *curr = gPad;
    TCanvas *ctemp = new TCanvas("ctemp","ctemp");
    
//    if (hist_high)  {
//        DoLog();
//        disp_hist->Draw("cont list");
//    }
//    else
    if (!hist_high) {
        hist->Draw("cont list");
        
        ctemp->Update();
        ctemp->Paint();
        
        TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->At(0);
        TList *lcontour1 = (TList*)contours->At(0);
        graph = new TGraph();
        graph->Merge(lcontour1);
//        TGraph *gc1 = (TGraph*)lcontour1->First();
//        graph = new TGraph(*gc1);
        delete ctemp; ctemp = 0;
        
        if (curr) curr->cd();
     
        if (IsLog()) DoLog();
        if (IsLogLin()) DoLogLin();
    }
    else {
        if (IsLog()) DoLog();
        if (IsLogLin()) DoLogLin();
        ChooseContour(disp_hist);
        disp_hist->Draw("cont list");
        
        ctemp->Update();
        ctemp->Paint();
        
        TObjArray *contours = (TObjArray*)gROOT->GetListOfSpecials()->At(0);
        TList *lcontour1 = (TList*)contours->At(0);
        disp_graph = new TGraph();
        graph = new TGraph();
        disp_graph->Merge(lcontour1);
        double px, py, y;
        for (int i = 0; i < disp_graph->GetN(); i++) {
            disp_graph->GetPoint(i, px, py);
            y = Can2Cont(py);
            graph->SetPoint(i, px, y);
        }
        
//        TString name = hist->GetName();
//        name += "_graph";
//        TGraph *gc = new TGraph();
//        gc->SetName(name);
//        TGraph *gci;
//        double px, py, y;
//        int n = 0;
//        
//        for (int i = 0; i < 10; i++) {
//            gci = (TGraph*)lcontour1->At(i);
//            if (!gci) {
//                MSG("NuContour",Msg::kInfo) << "  Merged " << i << " graphs" << endl;
//                break;
//            }
//            for (int j = 0; j < gci->GetN(); j++) {
//                gci->GetPoint(j, px, py);
//                y = Can2Cont(py);
//                gc->SetPoint(n, px, y);
//                n++;
//            }
//        }
//        if (addCorner) gc->SetPoint(n++, 1, ymax);
//        if (n != gc->GetN()) {
//            MSG("NuContour",Msg::kError) << "gc has the wrong number of points.  Should have "
//            << n << " but actually has " << gc->GetN() << endl;
//            assert(false);
//        }
//        graph = new TGraph(*gc);
//
//        if (IsLog()) DoLog();
//        if (IsLogLin()) DoLogLin();

        delete ctemp; ctemp = 0;
        
        if (curr) curr->cd();
        
    }
}

//____________________________________________________________________72
void NuContour::AddCorner(bool upper, bool right) 
{
    double x, y, yd;
    if (upper) {
        y = ymax;
        yd = bmax;
    }
    else {
        if (IsLog())  y = ymid;
        else          y = 0;
        yd = 0;
    }
    if (right) x = 1;
    else       x = 0;
       
    disp_graph->SetPoint(disp_graph->GetN(), x, yd);
    graph->SetPoint(graph->GetN(), x, y);
}


//____________________________________________________________________72
void NuContour::CompleteGraph() 
{
    double x, y;
    
    disp_graph->GetPoint(0, x, y);
    disp_graph->SetPoint(disp_graph->GetN(), x, y);

    graph->GetPoint(0, x, y);
    graph->SetPoint(graph->GetN(), x, y);
    

}


//____________________________________________________________________72
void NuContour::Draw(TString opt) 
{
    if (graph || disp_graph) {
        opt += "l";
        if (fillLevel >= 0) opt+="cf";
        MSG("NuContour",Msg::kInfo) << "Drawing graph with opt=" << opt <<  endl;
        GetGraph()->Draw(opt);
    }
    else {
        opt += "cont3";
        if (fillLevel >= 0) { 
            opt+="col";
        }
        MSG("NuContour",Msg::kInfo) << "Drawing histogram with opt=" << opt <<  endl;
        TH2D *htemp = new TH2D(*GetHist());
        htemp->Scale(neg);
        htemp->Draw(opt);
    }
}


//____________________________________________________________________72
void NuContour::GenerateLogLin()
{
    if (hist) {
        GenerateLogLinHist();
    }
    if (graph) {
        GenerateLogLinGraph();
    }
}


//____________________________________________________________________72
void NuContour::GenerateLog()
{
    GenerateLogLin();
//    if (hist) {
//        GenerateLogHist();
//        
//    }
}
    

//____________________________________________________________________72
void NuContour::GenerateLogLinHist()
{
    static int n = -1;
    n++;
    TString num(n);
    disp_hist = new TH2D("hCombo"+num,";sin^{2}(2#bar{#theta});|#Delta#bar{m}^{2}| (10^{-3} eV^{2})",
                      500, 0, 1,
                      bmax, 0, bmax);
    disp_hist->GetXaxis()->CenterTitle();
    disp_hist->GetYaxis()->CenterTitle();

    
    TAxis *a1 = disp_hist->GetYaxis();
    a1->SetLabelSize(0);
    a1->SetTickLength(0);
    a1->SetTitleOffset(1.2);
    
    //int old_bx, old_by;
    double val;
    double x, y;
    TH2D * hOld = 0; 
    
    for (int by = 1; by <= disp_hist->GetNbinsY(); by++) {
        y = Can2Cont(by);
        
        if (hist_high) {
            if (y > hdivide) hOld = hist_high;
            else             hOld = hist;
        }
        else {
            hOld = hist;
        }
        
        for (int bx = 1; bx <= disp_hist->GetNbinsX(); bx++) {
            x = disp_hist->GetBinCenter(bx);
            
            //old_bx = hOld->GetXaxis()->FindFixBin(x);
            //old_by = hOld->GetYaxis()->FindFixBin(y);
            //val = hOld->GetBinContent(old_bx, old_by);
            val = interpolate(hOld, x, y);
            
            disp_hist->SetBinContent(bx, by, val);
        }
    }
}


//____________________________________________________________________72
void NuContour::GenerateLogHist()
{
    static int n = -1;
    n++;
    TString num(n);
    disp_hist = new TH2D("hLog"+num,";sin^{2}(2#bar{#theta});|#Delta#bar{m}^{2}| (10^{-3} eV^{2})",
                         100, 0, 1,
                         500, ymid, ymax);
    disp_hist->GetXaxis()->CenterTitle();
    disp_hist->GetYaxis()->CenterTitle();
    
    
    TAxis *a1 = disp_hist->GetYaxis();
    //a1->SetLabelSize(0);
    //a1->SetTickLength(0);
    a1->SetTitleOffset(1.2);
    
    //int old_bx, old_by;
    double val;
    double x, y;
    TH2D * hOld = 0; 
    
    for (int by = 1; by <= disp_hist->GetNbinsY(); by++) {
        y = Can2Cont(by);
        
        if (hist_high) {
            if (y > hdivide) hOld = hist_high;
            else             hOld = hist;
        }
        else {
            hOld = hist;
        }
        
        for (int bx = 1; bx <= disp_hist->GetNbinsX(); bx++) {
            x = disp_hist->GetBinCenter(bx);
            
            //old_bx = hOld->GetXaxis()->FindFixBin(x);
            //old_by = hOld->GetYaxis()->FindFixBin(y);
            //val = hOld->GetBinContent(old_bx, old_by);
            val = interpolate(hOld, x, y);
            
            disp_hist->SetBinContent(bx, by, val);
        }
    }
}


//____________________________________________________________________72
void NuContour::GenerateLogLinGraph()
{
    if (graph == disp_graph) {
        MSG("NuContour",Msg::kInfo) << "Converting disp_graph has precednece -- no conversion" << endl;        
    }
    //cout << "Converting graph to disp_graph" << endl;
    disp_graph = new TGraph(*graph);
    TString name = graph->GetName();
    disp_graph->SetName(name+"LL");
    double px, py, y;
    
    for (int i = 0; i < disp_graph->GetN(); i++) {
        disp_graph->GetPoint(i, px, py);
        y = Cont2Can(py);

        disp_graph->SetPoint(i, px, y);
    }        
}


//____________________________________________________________________72
void NuContour::FixCanvas(TCanvas *c1)
{
    c1->SetCanvasSize(700,1000);
    c1->SetWindowSize(700,1000);    
    c1->SetLeftMargin(0.2);
    c1->SetRightMargin(0.05);
}


//____________________________________________________________________72
void NuContour::DrawAxes(int size)
{
//    if (IsLin()) {
//        MSG("NuContour",Msg::kInfo) << "Not drawing axes for linear plot" << endl;
//        return;
//    }
//    if (IsLog()) {
//        MSG("NuContour",Msg::kInfo) << "Not drawing axes for log plot" << endl;
//        return;
//    }        

    double gxmin = gPad->GetUxmin();
    double gymin = gPad->GetUymin();
    double gxmax = gPad->GetUxmax();
    double gymax = gPad->GetUymax();
    MSG("NuContour",Msg::kInfo) << "Drawing axes.  Found canvas range x: " << gxmin << " - " << gxmax << ", y: " << gymin << " - " << gymax << endl;
    int div_lin = 505;
    int div_log = 505;
    
    
    double x1 = gxmin;
    double x2 = gxmax;
    double y0 = gymin;
    double y1 = gymax*bmid/bmax;
    double y2 = gymax;
    
    //cout << "Using y0 = " << y0 << ", y1 = " << y1 << ", y2 = " << y2 << endl;

    if (IsLin()) {
        ymid = y2;
    }
    
    TF1 *fLin = new TF1("fLin","x",0,ymid);
    TF1 *fLog = new TF1("fLog","TMath::Log10(x)",ymid,ymax);
    // Prevent unused variable warnings
    fLin->GetName();
    fLog->GetName();
    
    
    if (IsLog() || IsLogLin()) {
        //cout << "Drawing log axes" << endl;
        // Draw Log Axes from bmid to bmax
        TGaxis *aLeft_log = new TGaxis(x1, y1, x1, y2, "fLog", div_log,"G");
        TGaxis *aRight_log = new TGaxis(x2, y1, x2, y2, "fLog", div_log, "+G");
        
        aLeft_log->SetLabelSize(0);
        aRight_log->SetLabelSize(0);
        aRight_log->SetTickSize(0.05);
        aLeft_log->SetTickSize(0.05);
        
        aLeft_log->SetGridLength(.75);
        
        aLeft_log->Draw();
        aRight_log->Draw();
    }
    
    if (IsLin() || IsLogLin()) {
        //cout << "Drawing lin axes" << endl;
        // Draw Linear Axes from 0 to bmid
        TGaxis *aLeft_lin = new TGaxis(x1, y0, x1, y1, "fLin", div_lin,"");
        TGaxis *aRight_lin = new TGaxis(x2, y0, x2, y1, "fLin", div_lin, "+");
        
        aLeft_lin->SetLabelSize(0);    
        aRight_lin->SetLabelSize(0);    
        aLeft_lin->SetTickSize(0.05);    
        aRight_lin->SetTickSize(0.05);    
        
        aLeft_lin->SetGridLength(.75);    
        
        aLeft_lin->Draw();
        aRight_lin->Draw();    
    }
    
    
    double yl, ln, yl_p;
    int nlab = 0;
    double log_num[99];
    
    
    if (IsLin() || IsLogLin()) {
        //cout << "Drawing lin labels" << endl;
        for (ln = 0; ln/1e3 <= ymid; ln += 2) {
            //cout << "Choosing label #" << nlab << " at " << ln << endl;
            log_num[nlab++] = ln;
        }
    }
    
    if (IsLog() || IsLogLin()) {
        for (ln = 1e-5 ; ln/1e3 <= ymax; ln *= 10) {
            if (ln/1e3 > ymid)
                log_num[nlab++] = ln;
            if (ln*2/1e3 <= ymax && ln/1e3*2 > ymid) 
                log_num[nlab++] = ln*2;
            if (ln*4/1e3 <= ymax && ln/1e3*4 > ymid) 
                log_num[nlab++] = ln*4;
        }
    }
    
    TLatex *lab_tex[99];
    yl_p = -999;
    for (int i = 0; i < nlab; i++) {
        yl = Cont2Can(log_num[i]/1e3);
        if (!IsLin()) {
            if ( yl - yl_p < 0.04*bmax ) {
                if (log_num[i] < ymid) log_num[i] += 2;
                else                   log_num[i] *= 2;
                yl = Cont2Can(log_num[i]/1e3);
            }
            yl_p = yl;
        }
        //cout << "Drawing label #" << i << " at " << yl << endl;
        lab_tex[i] = new TLatex(x1-0.01, yl, TString::Format("%3g",log_num[i]));
        if (size == -1) {
            lab_tex[i]->SetTextFont(42);
            lab_tex[i]->SetTextSize(0.05);
        }
        else {
            lab_tex[i]->SetTextFont(43);
            lab_tex[i]->SetTextSize(size);
        }
        lab_tex[i]->SetTextAlign(32);
        lab_tex[i]->Draw();
    }
    

    // Put in log-lin dividing line
    if (IsLogLin()) {
      TLine *l1 = new TLine(x1, y1, x2, y1);
      l1->SetLineStyle(2);
      l1->Draw();

      TLatex *t[99];
      int nT = -1;
      
      t[++nT] = new TLatex(x1+0.08, y1+0.01*y2, "Log");
      t[nT]->SetTextAlign(11);
      double defsize = t[nT]->GetTextSize();
      
      t[++nT] = new TLatex(x1+0.08, y1-0.01*y2, "Linear");
      t[nT]->SetTextAlign(13);
      
      
      for (int i = 0; i <= nT; i++) {
        if (size == -1) {
          if (t[i]->GetTextSize() == defsize) {
            t[i]->SetTextSize(0.045);
          }
          t[i]->SetTextFont(42);
        }
        else {
          t[i]->SetTextFont(43);
          t[i]->SetTextSize(size);
        }
        t[i]->Draw();
      }
    }
}


//____________________________________________________________________72
void NuContour::DispBins(int n) 
{
    bmax = n;
    
    if (IsLogLin()) DoLogLin();
    if (IsLog())    DoLog();
}

//____________________________________________________________________72
void NuContour::LogMax(double y) 
{
    ymax = y;
    if (IsLogLin()) DoLogLin();
    if (IsLog())    DoLog();
}


//____________________________________________________________________72
void NuContour::LogMin(double y) 
{
    ymid = y;
    if (IsLogLin()) DoLogLin();
    if (IsLog())    DoLog();
}


//____________________________________________________________________72
double NuContour::interpolate(TH2D *h, double x, double y, double def)
{
   // Given a point P(x,y), Interpolate approximates the value via bilinear 
   // interpolation based on the four nearest bin centers 
   // see Wikipedia, Bilinear Interpolation 
   // Andy Mastbaum 10/8/2008 
   // vaguely based on R.Raja 6-Sep-2008 

   Double_t f=0; 
   Double_t x1=0,x2=0,y1=0,y2=0; 
   Double_t dx,dy; 
   
   TAxis *fXaxis = h->GetXaxis();
   TAxis *fYaxis = h->GetYaxis();

   Int_t bin_x = fXaxis->FindBin(x); 
   Int_t bin_y = fYaxis->FindBin(y); 
   if(bin_x<1 || bin_x>h->GetNbinsX() || bin_y<1 || bin_y>h->GetNbinsY()) {
     MAXMSG("NuContour",Msg::kWarning,5) 
       << "Cannot interpolate outside histogram domain, using default value " 
       << def << endl;
      return 999.;     
   }

   Int_t quadrant = 0; // CCW from UR 1,2,3,4 
   // which quadrant of the bin (bin_P) are we in? 
   dx = fXaxis->GetBinUpEdge(bin_x)-x; 
   dy = fYaxis->GetBinUpEdge(bin_y)-y; 
   if (dx<=fXaxis->GetBinWidth(bin_x)/2 && dy<=fYaxis->GetBinWidth(bin_y)/2) 
     quadrant = 1; // upper right 
   if (dx>fXaxis->GetBinWidth(bin_x)/2 && dy<=fYaxis->GetBinWidth(bin_y)/2) 
     quadrant = 2; // upper left 
   if (dx>fXaxis->GetBinWidth(bin_x)/2 && dy>fYaxis->GetBinWidth(bin_y)/2) 
     quadrant = 3; // lower left 
   if (dx<=fXaxis->GetBinWidth(bin_x)/2 && dy>fYaxis->GetBinWidth(bin_y)/2) 
     quadrant = 4; // lower right 

   switch(quadrant) { 
   case 1: 
      x1 = fXaxis->GetBinCenter(bin_x); 
      y1 = fYaxis->GetBinCenter(bin_y); 
      x2 = fXaxis->GetBinCenter(bin_x+1); 
      y2 = fYaxis->GetBinCenter(bin_y+1); 
      break; 
   case 2: 
      x1 = fXaxis->GetBinCenter(bin_x-1); 
      y1 = fYaxis->GetBinCenter(bin_y); 
      x2 = fXaxis->GetBinCenter(bin_x); 
      y2 = fYaxis->GetBinCenter(bin_y+1); 
      break; 
   case 3: 
      x1 = fXaxis->GetBinCenter(bin_x-1); 
      y1 = fYaxis->GetBinCenter(bin_y-1); 
      x2 = fXaxis->GetBinCenter(bin_x); 
      y2 = fYaxis->GetBinCenter(bin_y); 
      break; 
   case 4: 
      x1 = fXaxis->GetBinCenter(bin_x); 
      y1 = fYaxis->GetBinCenter(bin_y-1); 
      x2 = fXaxis->GetBinCenter(bin_x+1); 
      y2 = fYaxis->GetBinCenter(bin_y); 
      break; 
   } 

   Int_t bin_x1 = fXaxis->FindBin(x1);
   if(bin_x1<1) bin_x1=1;
   Int_t bin_x2 = fXaxis->FindBin(x2);
   if(bin_x2>h->GetNbinsX()) bin_x2=h->GetNbinsX();
   Int_t bin_y1 = fYaxis->FindBin(y1);
   if(bin_y1<1) bin_y1=1;
   Int_t bin_y2 = fYaxis->FindBin(y2);
   if(bin_y2>h->GetNbinsY()) bin_y2=h->GetNbinsY();

   Int_t bin_q22 = h->GetBin(bin_x2,bin_y2);
   Int_t bin_q12 = h->GetBin(bin_x1,bin_y2);
   Int_t bin_q11 = h->GetBin(bin_x1,bin_y1);
   Int_t bin_q21 = h->GetBin(bin_x2,bin_y1);
   Double_t q11 = h->GetBinContent(bin_q11); 
   Double_t q12 = h->GetBinContent(bin_q12); 
   Double_t q21 = h->GetBinContent(bin_q21); 
   Double_t q22 = h->GetBinContent(bin_q22); 
   Double_t d = 1.0*(x2-x1)*(y2-y1);

   f = 1.0*q11/d*(x2-x)*(y2-y)+1.0*q21/d*(x-x1)*(y2-y)+1.0*q12/d*(x2-x)*(y-y1)+1.0*q22/d*(x-x1)*(y-y1);
   return f; 

}



//____________________________________________________________________72
double NuContour::Can2Cont(double bin)
{
//    if (IsLog()) {
//        return disp_hist->GetYaxis()->GetBinCenter(bin);
//    }
    if (IsLin()) return bin;
    
    double c = TMath::Log10(ymax/ymid)/(bmax-bmid);
    double y;
    if (bin < bmid) {
        y = ymid / bmid * bin;
    } 
    else {
        y = ymid*pow(10, c * (bin - bmid));
    }
    return y;
}


//____________________________________________________________________72
double NuContour::Cont2Can(double y)
{
    if (IsLin()) return y;
    
    double c = TMath::Log10(ymax/ymid)/(bmax-bmid);
    double bin;
    if (y <= ymid) {
        bin = y * bmid/ymid;
    }
    else {
        bin = 1/c * TMath::Log10(y/ymid)+bmid;
    }
    return bin;
}




//____________________________________________________________________72
Double_t NuContour::MinBin(const TH2D& hist)
{
    Double_t minimum = -1.;
    for(Int_t x = 1; x <= hist.GetNbinsX(); x++) {
        for(Int_t y = 1; y <= hist.GetNbinsY(); y++) {
            
            Double_t val = hist.GetBinContent(x, y);
            if (x == 1 && y == 1) {
                minimum = val;
            } else {
                minimum = TMath::Min(minimum, val);
            }
        }
    }
    
    return minimum;
}


//____________________________________________________________________72
void NuContour::ShiftHist(TH2D *hist, Double_t shift)
{
    for(Int_t x = 1; x <= hist->GetNbinsX(); x++) {
        for(Int_t y = 1; y <= hist->GetNbinsY(); y++) {
            Double_t val = hist->GetBinContent(x, y);
            hist->SetBinContent(x, y, val + shift);
        }
    }
}


//____________________________________________________________________72
void NuContour::ChooseContour(TH2 *h)
{   
    h->SetContour(fillLevel+1);
    
    if (fillLevel < 0) {
        h->SetContour(1);
        h->SetContourLevel(0, neg*contLevel);
    }
    else {
        h->SetContour(nFillCols);
        //h->SetContourLevel(0, -999);
        for (int i = 0; i < nFillCols; i++) {
            if (i >= nFillCols-fillLevel-1)  h->SetContourLevel(i, 0);
            else                             h->SetContourLevel(i, 999);
        }
    }
                    
//    
//    else if (fillLevel == 0) {
//        h->SetContour(3);
//        h->SetContourLevel(0, neg*contLevel);
//        h->SetContourLevel(1, neg*contLevel);
//        h->SetContourLevel(2, neg*contLevel);
//        //h->SetContourLevel(3, neg*999.);
//    }
//    else if (fillLevel == 1) {
//        h->SetContour(3);
//        h->SetContourLevel(0, neg*contLevel);
//        h->SetContourLevel(1, neg*contLevel);
//        h->SetContourLevel(2, neg*999.);
//        //h->SetContourLevel(3, neg*999.);
//    }
//    else if (fillLevel == 2) {
//        h->SetContour(3);
//        h->SetContourLevel(0, neg*contLevel);
//        h->SetContourLevel(1, neg*999.);
//        h->SetContourLevel(2, neg*999.);
//        //h->SetContourLevel(3, neg*999.);
//    }
//    
    
    h->SetDirectory(0);
}

---