NumuCC2p2hAnalysis.cxx

Go to the documentation of this file.

#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Unfold/UnfoldIterative.h"
#include "CAFAna/Analysis/Plots.h"

#include "NDAna/numucc_2p2h/NumuCC2p2hAnalysis.h"

#include "TH3.h"                                                                   
#include "TFile.h"
#include "TCanvas.h"
#include "TVectorD.h" 

#include <cassert>
#include <vector>
#include <string>       // std::string
#include <iostream>     // std::cout
#include <sstream>      // std::stringstream
#include <iomanip>      // std::setprecision

namespace ana
{
  TH2* NumuCC2p2hAnalysis::Result2D()
  {
    Spectrum bkgd = fBkgdEst->Background();
    Spectrum sig  = SignalEst(bkgd);
    Spectrum sig_unfolded = UnfoldedSignal(sig);
    Ratio eff = Efficiency();
    double ntargs = NTargets();
    //double flux = fFlux->Integral(fData.POT());                             
    TH1* hflux = fFlux->ToTH1(fData.POT());
    hflux->Scale(1e-4); // Convert nu/m^2 to nu/cm^2
    double flux = hflux->Integral();
    TH2* xsec = (sig_unfolded/eff).ToTH2(fData.POT());
    for(int iBinX = 1; iBinX<=xsec->GetXaxis()->GetNbins(); iBinX++)
      {
        for(int iBinY = 1; iBinY<=xsec->GetYaxis()->GetNbins(); iBinY++)
          {
            double binwidth = xsec->GetXaxis()->GetBinWidth(iBinX)*xsec->GetYaxis()->GetBinWidth(iBinY);
            double newcontent = xsec->GetBinContent(iBinX,iBinY)/binwidth;
            xsec->SetBinContent(iBinX,iBinY,newcontent);
          }
      }
    //xsec->Divide(hflux);  
    xsec->Scale(1./ntargs);
    xsec->Scale(1./flux);
    //xsec->Scale(1./(nnuc*flux));        

    return xsec;
  }

  TH1* NumuCC2p2hAnalysis::Result1DFluxInt()
  {
    Spectrum bkgd = fBkgdEst->Background();
    Spectrum sig  = SignalEst(bkgd);
    Spectrum sig_unfolded = UnfoldedSignal(sig);
    Ratio eff = Efficiency();
    double ntargs = NTargets();
    //double flux = fFlux->Integral(fData.POT());    
    TH1* hflux = fFlux->ToTH1(fData.POT());
    hflux->Scale(1e-4); // Convert nu/m^2 to nu/cm^2      
    double flux = hflux->Integral();
    TH1* xsec = (sig_unfolded/eff).ToTH2(fData.POT());
    for(int iBinX = 1; iBinX<=xsec->GetXaxis()->GetNbins(); iBinX++)
      {
        double binwidth = xsec->GetXaxis()->GetBinWidth(iBinX);
        double newcontent = xsec->GetBinContent(iBinX)/binwidth;
        xsec->SetBinContent(iBinX,newcontent);
      }
    
    //xsec->Divide(hflux);         
    xsec->Scale(1./ntargs);
    xsec->Scale(1./flux);
    //xsec->Scale(1./(nnuc*flux));   

    return xsec;
  }


  TH2* NumuCC2p2hAnalysis::PlotEfficiency2D()
  {
    return Efficiency().ToTH2();
  }

  //----------------------------------------------------------------------                                
  TH2* NumuCC2p2hAnalysis::PlotSignalEstimate2D()
  {
    return SignalEst(fBkgdEst->Background()).ToTH2(fData.POT());
  }

  //----------------------------------------------------------------------                  

  TH2* NumuCC2p2hAnalysis::PlotUnfoldedSignal2D()
  {
    Spectrum bkgd = fBkgdEst->Background();
    return UnfoldedSignal(SignalEst(bkgd)).ToTH2(fData.POT());
  }

  //----------------------------------------------------------------------                        

  TH1* NumuCC2p2hAnalysis::PlotFluxEstimate2D()
  {
    TH1 *h = fFlux->ToTH1(1);
    //h->Scale(1e-4);             
    h->GetYaxis()->SetTitle("Flux / POT / m^{2}");
    return h;
  }

  //----------------------------------------------------------------------

  TH2* NumuCC2p2hAnalysis::PlotBackgroundEstimate2D()
  {
    Spectrum bkgd = fBkgdEst->Background();
    return bkgd.ToTH2(fData.POT());
  }

  //----------------------------------------------------------------------                               

  TH2* NumuCC2p2hAnalysis::PlotData2D()
  {
    return fData.ToTH2(fData.POT());
  }

  //----------------------------------------------------------------------                

  TH2D* NumuCC2p2hAnalysis::PlotRecoToTrueMatrix2D()
  {
    return fRecoTrue.ToTH2(fData.POT());
  }
  //----------------------------------------------------------------------
  std::unique_ptr<NumuCC2p2hAnalysis> NumuCC2p2hAnalysis::LoadFrom(TDirectory* dir){
    return LoadFrom(dir, kIterative);
  }

  //----------------------------------------------------------------------    
  std::unique_ptr<NumuCC2p2hAnalysis> NumuCC2p2hAnalysis::LoadFrom(TDirectory* dir,
                                                                   UnfoldMethod_t unfoldmethod)
  {
    Spectrum* data = ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fData")).release();
    Spectrum* mc   = ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fMC")).release();
    Spectrum* mcsig= ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fMCSig")).release();
    Spectrum* mcsignutree = ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fMCSigNuTree")).release();
    ReweightableSpectrum* rt = ana::LoadFrom<ReweightableSpectrum>
      (dir->GetDirectory("fRecoTrue")).release();

    IBkgdEstimator* bkgdest = ana::LoadFrom<IBkgdEstimator>
      (dir->GetDirectory("fBkgdEst")).release();

    TVectorD* unfoldreg = (TVectorD*)dir->Get("fUnfoldReg");
    TVector3* min = (TVector3*)dir->Get("fFidMin");
    TVector3* max = (TVector3*)dir->Get("fFidMax");

    Spectrum* flux = ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fFlux")).release();
    
    return std::unique_ptr<NumuCC2p2hAnalysis>(new NumuCC2p2hAnalysis(*data, *mc, mcsig, mcsignutree, *rt, bkgdest,
                                                                      (*unfoldreg)[0], unfoldmethod, *min, *max, flux));
  }
  //-----------------------------------------------------------------------
  std::unique_ptr<NumuCC2p2hAnalysis> NumuCC2p2hAnalysis::LoadFrom(TDirectory* dir,
                                                                   UnfoldMethod_t unfoldmethod,
                                                                   double nTargets)
  {
    Spectrum* data = ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fData")).release();
    Spectrum* mc   = ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fMC")).release();
    Spectrum* mcsig= ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fMCSig")).release();
    Spectrum* mcsignutree = ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fMCSigNuTree")).release();
    ReweightableSpectrum* rt = ana::LoadFrom<ReweightableSpectrum>
      (dir->GetDirectory("fRecoTrue")).release();

    IBkgdEstimator* bkgdest = ana::LoadFrom<IBkgdEstimator>
      (dir->GetDirectory("fBkgdEst")).release();

    TVectorD* unfoldreg = (TVectorD*)dir->Get("fUnfoldReg");
    TVector3* min = (TVector3*)dir->Get("fFidMin");
    TVector3* max = (TVector3*)dir->Get("fFidMax");

    Spectrum* flux = ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fFlux")).release();

    return std::unique_ptr<NumuCC2p2hAnalysis>(new NumuCC2p2hAnalysis(*data, *mc, mcsig, mcsignutree, *rt, bkgdest, nTargets,
                                                                      (*unfoldreg)[0], unfoldmethod, *min, *max, flux));
  }
}