make_nus17_pot.C

Go to the documentation of this file.

// Calculate analysis exposure, selected events as a function of run

#include <bitset>
#include <cstdlib>
#include <iostream>
#include "NuXAna/macros/Nus17/Nus17Exposures.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 "NuXAna/macros/Nus17/NusLoadProd3.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 fPOT(     [](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*kNus17FiducialMass[mydb];
                              mydb = 0;
                            }
                            dbcount++;
                         }
                         if (mydb >= 4) sum += live*kNus17FiducialMass[mydb];
                         return sum;
                       });
const SpillVar kLive(  [](const caf::SRSpillProxy* spill)
                       {
                         if (spill->run < 17946)
                           return 24e-6;
                         else
                           return 12e-6;
                       });
const SpillVar fPOTMass(
                       [](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*kNus17FiducialMass[mydb];
                              mydb = 0;
                            }
                            dbcount++;
                         }
                         if (mydb >= 4) pot += spill->spillpot*kNus17FiducialMass[mydb];
                         return pot;
                       });

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


void make_nus17_pot(int period)
{
  TH1::SetDefaultSumw2();

  //binned in unix time encompassing first analysis period
  TH1D* hLiveMass = new TH1D("hLiveMass","livetime exposure;run;exposure (s #times kT)",12473, 12939.5, 25412.5);
  TH1D* hLive = new TH1D("hLive","livetime;run;livetime (s)",12473, 12939.5, 25412.5);
  TH1D* hPOTMass = new TH1D("hPOTMass","exposure;run;exposure (POT #times kT)",12473, 12939.5, 25412.5);
  TH1D* hPOT = new TH1D("hPOT","POT;run;POT",12473, 12939.5, 25412.5);
  TH1D* hLiveMassT = new TH1D("hLiveMassT","livetime exposure;run;exposure (s #times kT)",9900,1391000000,1490000000);
  TH1D* hLiveT = new TH1D("hLiveT","livetime;run;livetime (s)",9900,1391000000,1490000000);
  TH1D* hPOTMassT = new TH1D("hPOTMassT","exposure;run;exposure (POT #times kT)",9900,1391000000,1490000000);
  TH1D* hPOTT = new TH1D("hPOTT","POT;run;POT",9900,1391000000,1490000000);
  
  const std::vector<std::string> fNumiUnblind(MakeNumiUnblindList());

  std::vector<std::string> fPeriod;
  if(period == 0) fPeriod = fNumiUnblind;
  else if(period == 1) fPeriod = fnumi_period1;
  else if(period == 2) fPeriod = fnumi_period2;
  else if(period == 3) fPeriod = fnumi_period3;
  else if(period == 5) fPeriod = fnumi_period5;
  else {
    std::cout << "Invalid period. Select a number from 0-5, excl. 4" << std::endl;
    std::cout << "0 selects all periods." << std::endl;
    abort();
  }
  //load data files
  SpectrumLoader loaderFDData(fPeriod, kCosmic);

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


  std::cout << "********* Period " << std::to_string(period)
            << " *********" << std::endl;
  
  std::cout << "Live Integral:      " << hLive->Integral() << std::endl;
  std::cout << "LiveT Integral:     " << hLiveT->Integral() << std::endl;
  std::cout << "LiveMass Integral:  " << hLiveMass->Integral() << std::endl;
  std::cout << "LiveMassT Integral: " << hLiveMassT->Integral() << std::endl;
  std::cout << "POT Integral:      " << hPOT->Integral() << std::endl;
  std::cout << "POTT Integral:     " << hPOTT->Integral() << std::endl;
  std::cout << "POTMass Integral:  " << hPOTMass->Integral() << std::endl;
  std::cout << "POTMassT Integral: " << hPOTMassT->Integral() << std::endl;


  std::string outname = "nus17_pot_period"+std::to_string(period)+".root";
  
  TFile* outFile = new TFile(outname.c_str(), "RECREATE");
  hLiveMass->Write();
  hLive->Write();
  hPOTMass->Write();
  hPOT->Write();
  hLiveMassT->Write();
  hLiveT->Write();
  hPOTMassT->Write();
  hPOTT->Write();
  outFile->Close();

}