FCSurface.cxx

Go to the documentation of this file.

#include "CAFAna/FC/FCSurface.h"

#include "CAFAna/FC/FCPoint.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/Utilities.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TH2.h"
#include "TTree.h"
#include "TVirtualPad.h"

#include <cassert>
#include <cmath>
#include <iostream>

namespace ana
{
  //----------------------------------------------------------------------
  FCSurface::FCSurface(int nbinsx, double xmin, double xmax,
                       int nbinsy, double ymin, double ymax)
  {
    DontAddDirectory guard;

    fBinning = new TH2F(UniqueName().c_str(), "",
                        nbinsx, xmin, xmax, nbinsy, ymin, ymax);

    fBins.resize((nbinsx+2)*(nbinsy+2));
    fPrintedWarning = false;
  }

  //----------------------------------------------------------------------
  FCSurface::~FCSurface()
  {
    delete fBinning;
  }

  //----------------------------------------------------------------------
  void FCSurface::Add(const FCPoint& pt, std::string plot)
  {
    if(pt.Up()<0) 
      return;

    int bin;
    if(plot == "delta")
      bin = fBinning->FindBin(pt.TrueDelta(), 0);
    else if(plot == "ssth23") 
      bin = fBinning->FindBin(pt.TrueTh23(), 0);
    else if(plot == "th13")
      bin = fBinning->FindBin(pt.TrueTh13(), 0);
    else if(plot == "dmsq")
      bin = fBinning->FindBin(pt.TrueDmsq(), 0);
    else if(plot == "deltassth23")
      bin = fBinning->FindBin(pt.TrueDelta(), pt.TrueTh23());
    else if(plot == "ssth23dmsq32") {
      bin = fBinning->FindBin(pt.TrueTh23(), pt.TrueDmsq());
    }
    else if(plot == "deltath13")
      bin = fBinning->FindBin(pt.TrueTh13(), pt.TrueDelta());
    else if(plot == "mass")
      bin = fBinning->FindBin(pt.ProfDmsq(), 0);
    else {
      std::cout << "I don't know what plot you're trying to make." << std::endl;
      std::cout << "You might have to add the option to FCSurface::Add(). " << std::endl;
      exit(1);
    }
    
    // For deltassth23, since the first and last x bins cover the same range
    // in sin(delta_CP), we want to merge these two columns of bins
    if(plot == "deltassth23"){
      // Get binwidth so we can decide whether trueDelta falls in either
      // of these bins
      double binWidth = fBinning->GetXaxis()->GetBinWidth(1);
      // If trueDelta falls into the first bin, also add it to the last
      if(pt.TrueDelta() > -1*binWidth/2 && pt.TrueDelta() < binWidth/2) {
        int mergeBin = fBinning->FindBin(pt.TrueDelta()+2, pt.TrueTh23());
        fBins[mergeBin].Add(pt);
      }
      // If trueDelta falls into the last bin, also add it to the first
      if(pt.TrueDelta() > 2-binWidth/2 && pt.TrueDelta() < 2+binWidth/2) {
        int mergeBin = fBinning->FindBin(pt.TrueDelta()-2, pt.TrueTh23());
        fBins[mergeBin].Add(pt);
      }   
    }
    
    fBins[bin].Add(pt);
  }

  //----------------------------------------------------------------------
  void FCSurface::Add(const FCCollection& coln, std::string plot)
  {
    Progress prog("Filling FCSurface from collection");
    const unsigned int N = coln.NPoints();

    for(unsigned int n = 0; n < N; ++n){
      if(n%100000 == 0) prog.SetProgress(double(n)/N);
      Add(coln.Point(n), plot);
    }
    prog.Done();
  }

  //----------------------------------------------------------------------
  TH2* FCSurface::SurfaceForSignificance(double sig)
  {
    DontAddDirectory guard;

    TH2* hRet = new TH2F(*fBinning);
    hRet->SetName(UniqueName().c_str());

    // Don't do overflow and underflow
    for(int x = 1; x < hRet->GetNbinsX()+1; ++x){
      for(int y = 1; y < hRet->GetNbinsY()+1; ++y){
        const int bin = hRet->GetBin(x, y);
        if(fBins[bin].Empty()){
          double crit2D = -2*log(1-sig);
          hRet->SetBinContent(bin,crit2D);
          if (!fPrintedWarning) std::cout << "WARNING, NO FCPOINTS IN THIS BIN, fall back on 2D Gaussian up value." << std::endl;
          fPrintedWarning = true;
          continue;
        }
        const double up = fBins[bin].UpValueForSignificance(sig);
        assert(!std::isinf(up) && !std::isnan(up));
        hRet->SetBinContent(bin, up);
      }
    }

    hRet->SetMinimum(0);

    return hRet;
  }

  //----------------------------------------------------------------------
  void FCSurface::DrawBins()
  {
    TVirtualPad* c = gPad;
    if(!c) c = new TCanvas;

    const int X = fBinning->GetNbinsX();
    const int Y = fBinning->GetNbinsY();
    c->Divide(X, Y);

    for(int x = 1; x < X+1; ++x){
      for(int y = 1; y < Y+1; ++y){
        c->cd(x+(y-1)*X);
        const int bin = fBinning->GetBin(x, y);
        fBins[bin].Draw();
      }
    }
  }

  //----------------------------------------------------------------------
  const FCBin* FCSurface::GetBin(int x, int y) const
  {
    return &fBins[fBinning->GetBin(x, y)];
  }

  //----------------------------------------------------------------------
  void FCSurface::SaveToFile(const std::string& fname) const
  {
    TFile fout(fname.c_str(), "RECREATE");
    fBinning->Write("bins");

    TTree tr("fcsurf", "fcsurf");
    FCBin::InitToTree(&tr);
    for(const FCBin& bin : fBins) bin.ToTree(&tr);
    tr.Write();
  }

  //----------------------------------------------------------------------
  std::unique_ptr<FCSurface> FCSurface::LoadFromFile(const std::string& fname)
  {
    DontAddDirectory guard;

    std::unique_ptr<FCSurface> ret(new FCSurface);

    TFile* fin = TFile::Open(pnfs2xrootd(fname).c_str());
    assert(!fin->IsZombie());

    ret->fBinning = (TH2F*)(fin->Get("bins")->Clone());

    TTree* tr = (TTree*)fin->Get("fcsurf");
    assert(tr);

    FCBin::InitFromTree(tr);

    for(int i = 0; i < tr->GetEntries(); ++i)
      ret->fBins.push_back(FCBin::FromTree(tr, i));

    delete fin;

    return ret;
  }

} // namespace