get_numu_data_histogram.C

Go to the documentation of this file.

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Vars/HistAxes.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"

using namespace ana;

#include "TFile.h"
#include <iostream>

void get_numu_data_histogram(bool isFHC = true)
{
    /// File for making data histograms for fit
    /// Use it for extracting OLD data hitstograms from decafs made in 2019 release

    bool getfiles=0;
    std::string horn_tag = "fhc";
    if (!isFHC) horn_tag = "rhc";

    if (getfiles){

    std::string fnameData = "defname: prod4_sumrestricteddecaf_fd_numi_fhc_full_goodruns_numu2019";
    //if (!isFHC) fnameData = "defname: prod4_sumrestricteddecaf_fd_numi_rhc_Ana2018_goodruns_numu2019";
    if (!isFHC) fnameData = "defname: prod4_sumrestricteddecaf_fd_numi_rhc_full_goodruns_numu2019";

    SpectrumLoader loaderData(fnameData);
    loaderData.SetSpillCut(kStandardSpillCuts);

    std::string infile_quant = "/cvmfs/nova.opensciencegrid.org/externals/numudata/v00.00/NULL/ana2018/Quantiles/quantiles__"+horn_tag+"_full__numu2018.root";
    TFile* infile = new TFile(infile_quant.c_str());
    TH2* FDSpec2D = (TH2*)infile->FindObjectAny("FDSpec2D");
    std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
    infile->Close();
    delete infile;
    
    std::vector<Spectrum*> data;
    data.push_back(new Spectrum(loaderData, kNumuCCOptimisedAxis, kNumuCutFD2018&&kInBeamSpill,  kNoShift));
    for(auto thisCut : HadEFracQuantCuts){
        data.push_back(new Spectrum(loaderData, kNumuCCOptimisedAxis, kNumuCutFD2018 && kInBeamSpill && thisCut));
    }

    loaderData.Go();
    
    TFile* fout = new TFile(("numu_data_2019_"+horn_tag+".root").c_str(), "RECREATE");
    for(int quant=0; quant<5; quant++){
        data[quant]-> SaveTo(&fout->, Form("data_q%d", quant));
    }
    delete fout;
    }
// in order to get spectra from *.root file do:
   if (!getfiles){

    TFile * file = new TFile (("numu_data_2019_"+horn_tag+"_test.root").c_str(), "READ");
    std::vector<Spectrum*> data;
    for(int quant=0; quant<5; quant++){
         auto data = Spectrum::LoadFrom(file, Form("data_q%d", quant));
         TH1* hist = data->ToTH1(data->POT());        
         cout<<" data:  "<<hist->Integral()<<" POT "<<data->POT()<<endl;
    }
   }

}