make_nue_ana2019_rhc_ub_pot.C

Go to the documentation of this file.

#include <bitset>
#include <iostream>
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Cuts/NueCutsFirstAna.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "CAFAna/Vars/HistAxes.h"
#include "CAFAna/Core/EventList.h"
#include "StandardRecord/Proxy/SRProxy.h"
using namespace ana;

#include "TFile.h"
#include "TH1D.h"

//fiducial mass accounting from docdb-13556

const SpillVar kSpillTime([](const caf::SRSpillProxy* spill)
                         {
                           return spill->spilltimesec;
                         });
const SpillVar kRun(     [](const caf::SRSpillProxy* spill)
                         {
                           return spill->run;
                         });
const SpillVar kSpillPOT([](const caf::SRSpillProxy* spill)
                         {
                           return spill->spillpot;
                         });


const SpillVar kLiveMass(
                       [](const caf::SRSpillProxy* spill)
                       {
                          std::string binary = std::bitset<14>(spill->dibmask).to_string();
                          double live = 12e-6;
                          if (spill->run < 17946) live = 24e-6;
                          double sum = 0;
                          int dbcount = 1;
                          int mydb = 0;
                          for (int i=13; i>-1; --i){
                            int test = binary[i] - '0';
                            if (dbcount >= spill->dibfirst && 
                                dbcount <= spill->diblast && 
                                test == 1) mydb++;
                            if (dbcount >= spill->dibfirst && 
                                dbcount <= spill->diblast && 
                                test == 0){
                              if (mydb >= 4) sum+= live*kAna2017FiducialMass[mydb];
                              mydb = 0;
                            }
                            dbcount++;
                         }
                         if (mydb >= 4) sum += live*kAna2017FiducialMass[mydb];
                         return sum;
                       });
const SpillVar kLive(  [](const caf::SRSpillProxy* spill)
                       {
                         if (spill->run < 17946)
                           return 24e-6;
                         else
                           return 12e-6;
                       });
const SpillVar kPOTMass(
                       [](const caf::SRSpillProxy* spill)
                       {
                          std::string binary = std::bitset<14>(spill->dibmask).to_string();
                          double pot = 0;
                          int dbcount = 1;
                          int mydb = 0;
                          for (int i=13; i>-1; --i){
                            int test = binary[i] - '0';
                            if (dbcount >= spill->dibfirst && 
                                dbcount <= spill->diblast && 
                                test == 1) mydb++;
                            if (dbcount >= spill->dibfirst && 
                                dbcount <= spill->diblast && 
                                test == 0){
                              if (mydb >= 4) pot+= spill->spillpot*kAna2017FiducialMass[mydb];
                              mydb = 0;
                            }
                            dbcount++;
                         }
                         if (mydb >= 4) pot += spill->spillpot*kAna2017FiducialMass[mydb];
                         return pot;
                       });

const Var kSliceTime(
                    [](const caf::SRProxy* sr)
                    {
                      return sr->hdr.unixtime;
                    });


void make_nue_ana2019_rhc_ub_pot()
{
  TH1::SetDefaultSumw2();

  //binned in unix time encompassing first analysis period
  TH1D* hLiveMass = new TH1D("hLiveMass","livetime exposure;run;exposure (s #times kT)",19458, 12939.5, 32397.5);
  TH1D* hLive = new TH1D("hLive","livetime;run;livetime (s)",19458, 12939.5, 32397.5);
  TH1D* hPOTMass = new TH1D("hPOTMass","exposure;run;exposure (POT #times kT)",19458, 12939.5, 32397.5);
  TH1D* hPOT = new TH1D("hPOT","POT;run;POT",19458, 12939.5, 32397.5);
  TH1D* hLiveMassT = new TH1D("hLiveMassT","livetime exposure;run;exposure (s #times kT)",16900,1391000000,1560000000);
  TH1D* hLiveT = new TH1D("hLiveT","livetime;run;livetime (s)",16900,1391000000,1560000000);
  TH1D* hPOTMassT = new TH1D("hPOTMassT","exposure;run;exposure (POT #times kT)",16900,1391000000,1560000000);
  TH1D* hPOTT = new TH1D("hPOTT","POT;run;POT",16900,1391000000,1560000000);
  

  //load data files
  SpectrumLoader loaderFDData("prod4_sumrestricteddecaf_fd_numi_rhc_full_goodruns_nue2019", kCosmic);

  loaderFDData.SetSpillCut(kStandardSpillCuts);
  loaderFDData.AddSpillHistogram(hLiveMass,kRun,kStandardSpillCuts, kLiveMass);
  loaderFDData.AddSpillHistogram(hLive,kRun,kStandardSpillCuts, kLive);
  loaderFDData.AddSpillHistogram(hPOTMass,kRun,kStandardSpillCuts, kPOTMass);
  loaderFDData.AddSpillHistogram(hPOT,kRun,kStandardSpillCuts, kSpillPOT);
  loaderFDData.AddSpillHistogram(hLiveMassT,kSpillTime,kStandardSpillCuts, kLiveMass);
  loaderFDData.AddSpillHistogram(hLiveT,kSpillTime,kStandardSpillCuts, kLive);
  loaderFDData.AddSpillHistogram(hPOTMassT,kSpillTime,kStandardSpillCuts, kPOTMass);
  loaderFDData.AddSpillHistogram(hPOTT,kSpillTime,kStandardSpillCuts, kSpillPOT);
  loaderFDData.Go();

  TFile* outFile = new TFile("nue_ana2019_rhc_pot.root", "RECREATE");
  hLiveMass->Write();
  hLive->Write();
  hPOTMass->Write();
  hPOT->Write();
  hLiveMassT->Write();
  hLiveT->Write();
  hPOTMassT->Write();
  hPOTT->Write();
  outFile->Close();

}