make_DataMCComp_numu.C

Go to the documentation of this file.

//------------------------------------------------------------------------------
// This makes the noextrap predictions of a bunch of interesting numu variables.
// The map of variables are defined at the bottom of the header.
// Set the period arg to 0 to use the full datasets
// Plots are made using the associated plotting macro.
//
// Andrew Sutton - ats2mk@virginia.edu
//------------------------------------------------------------------------------

#pragma once

#include "DataMCComp_header.h"

void make_DataMCComp_numu(std::string horn, uint period)
{
  // Input options
  Loaders::FluxType flux;
  if      (horn == "fhc") flux = Loaders::kFHC;
  else if (horn == "rhc") flux = Loaders::kRHC;
  else {
    std::cout << "\n\n********** No valid horn **********\n" << std::endl;
    std::abort();
  }

  std::string sPeriod = "period" + std::to_string(period);
  std::string sRange  = "period" + std::to_string(period);
  Cut kPeriodCut = kNoCut;
  if (period == 0) {
    sPeriod = "full";
    sRange  = "full";
  }
  if (period > 9) {
    std::cout << "\n\n********** No valid period **********\n" << std::endl;
    std::abort();
  }

  // Initial setup
  Loaders *mcLoader = new Prod5NomLoaders(ECAFType::kNumuConcat, flux, sPeriod);
  SpectrumLoader csLoader("prod_sumdecaf_R19-11-18-prod5reco.i_fd_cosmic_" + horn + "_" + sPeriod + "_v1_goodruns_numu2020");
  SpectrumLoader dtLoader("prod_sumrestricteddecaf_R19-11-18-prod5reco.combi_opt2_fd_numi_" + horn + "_" + sRange +"_v1_numu2020_fddatav2");
  mcLoader->SetSpillCut(kStandardSpillCuts);
  csLoader.SetSpillCut(kStandardSpillCuts);
  dtLoader.SetSpillCut(kStandardSpillCuts_FD_Prod4MotivatedDQ);
  
  Cut myCut = kNumu2020FD && kPeriodCut;
  Var mcWei = kPPFXFluxCVWgt * kXSecCVWgt2020;
    
  // Get Quants
  std::string infile_quant = "/cvmfs/nova.opensciencegrid.org/externals/numudata/v00.05/NULL/lib/ana2020/Quantiles/quantiles_" + horn + "_full_numu2020.root";
  TFile* infile = new TFile(infile_quant.c_str());
  TH2* FDSpec2D = (TH2*)infile->FindObjectAny("FDSpec2D");
  std::vector<Cut> quantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4, false);
  quantCuts.push_back(kNoCut);
  infile->Close();
  delete infile;
  const std::map< std::string, Cut> quantiles =
    {
      {"Quant1"   , quantCuts[0]},
      {"Quant2"   , quantCuts[1]},
      {"Quant3"   , quantCuts[2]},
      {"Quant4"   , quantCuts[3]},
      {"AllQuants", quantCuts[4]}
    };

  // Make the preds
  std::map< std::string, PredictionNoExtrap* > mcPreds;
  std::map< std::string, Spectrum*           > csSpecs;
  std::map< std::string, Spectrum*           > dtSpecs;
  for (auto var : vars) {

    for (auto quant : quantiles) {
        mcPreds["mcPred_" + var.first.first + "_" + quant.first] =
          new PredictionNoExtrap(*mcLoader, var.first.second, var.second.bin,
                                 var.second.var, myCut && quant.second,
                                 kNoShift, mcWei);

        csSpecs["csSpec_" + var.first.first + "_" + quant.first] =
          new Spectrum(var.first.second, var.second.bin, csLoader,
                       var.second.var, myCut && quant.second && kInCosmicTimingWindow,
                       kNoShift, kUnweighted);
                                                         
        dtSpecs["dtSpec_" + var.first.first + "_" + quant.first] =
          new Spectrum(var.first.second, var.second.bin, dtLoader,
                       var.second.var, myCut && quant.second && kInBeamSpill,
                       kNoShift, kUnweighted);

      }// end quant loop

  }

  // Liftoff
  mcLoader->Go();
  csLoader.Go();
  dtLoader.Go();

  TFile *outfile = new TFile(("numuDataMC_" + horn + "_" + sRange + ".root").c_str(), "RECREATE");


  for (auto pred : mcPreds)
    pred.second->SaveTo(outfile, pred.first.c_str());
  
  for (auto spec : csSpecs)
    spec.second->SaveTo(outfile, spec.first.c_str());

  for (auto spec : dtSpecs)
    spec.second->SaveTo(outfile, spec.first.c_str());
  
}