savePrediction_complete.C

Go to the documentation of this file.

#ifdef __CINT__
void savePrediction_complete(std::string& fEpoch){
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#include "../WorkshopIncludes.h"

void savePrediction_complete(std::string& fEpoch)
{
  std::cout << "Drawing plots..." << std::endl;

  // Production Concats
  std::string directoryD   = "/pnfs/nova/persistent/production/concat/R16-03-03-prod2reco.d/";
  std::string directoryF   = "/pnfs/nova/persistent/production/concat/R16-03-03-prod2reco.f/";
  // Data
  std::string nd_data      = directoryD+"prod_decaf_R16-03-03-prod2reco.d_nd_numi_fhc_full_numu_contain_v1_goodruns_prod2-snapshot.root";

  // Baseline MC                                                                                    
  std::string nd_nonswap  = directoryD+"prod_decaf_R16-03-03-prod2reco.d_nd_genie_nonswap_genierw_f\
hc_nova_v08_epoch1-3c_numu_contain_v1_prod2-snapshot/*.root";
  std::string fd_nonswap  = directoryF+"/"+fEpoch+"/prod_decaf_R16-03-03-prod2reco.f_fd_genie_nonsw\
ap_fhc_nova_v08_"+fEpoch+"_numu_contain_v1_prod2-snapshot.root";
  std::string fd_fluxswap = directoryF+"/"+fEpoch+"/prod_decaf_R16-03-03-prod2reco.f_fd_genie_fluxs\
wap_fhc_nova_v08_"+fEpoch+"_numu_contain_v1_prod2-snapshot.root";
  std::string fd_tau      = directoryF+"/"+fEpoch+"/prod_decaf_R16-03-03-prod2reco.f_fd_genie_tau_f\
hc_nova_v08_"+fEpoch+"_numu_contain_v1_prod2-snapshot.root";

  Loaders loaders;
  loaders.SetLoaderPath( nd_data,     caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap );
  loaders.SetLoaderPath( nd_nonswap,  caf::kNEARDET, Loaders::kMC,   ana::kBeam, Loaders::kNonSwap );
  loaders.SetLoaderPath( fd_nonswap,  caf::kFARDET,  Loaders::kMC,   ana::kBeam, Loaders::kNonSwap );
  loaders.SetLoaderPath( fd_fluxswap, caf::kFARDET,  Loaders::kMC,   ana::kBeam, Loaders::kFluxSwap);
  loaders.SetLoaderPath( fd_tau,      caf::kFARDET,  Loaders::kMC,   ana::kBeam, Loaders::kTauSwap );

  loaders.SetSpillCut(kStandardSpillCuts);

  osc::OscCalculatorPMNSOpt calc;
  ana::ResetOscCalcToDefault(&calc); // This seeds to default (in particular max-mixing)

  calc.SetTh23(TMath::ASin(sqrt(0.403))); // Second analysis best fit                             
  calc.SetDmsq32(2.67e-3); // Same

  std::vector<const ISyst*> systs = getAllNumuSysts();

  TFile* f = new TFile(("../Predictions/Prediction_"+fEpoch+"_systs_kCCE.root").c_str(),"RECREATE");

  const Cut kUglyRemID(
    {"sel.remid.pid"},
    [](const caf::StandardRecord* sr)
      { return (sr->sel.remid.pid > 0.25 && sr->sel.remid.pid < 0.75 );
      }
  );

  const Cut kBadRemID(
    {"sel.remid.pid"},
    [](const caf::StandardRecord* sr)
      { return (sr->sel.remid.pid >= 0.75 && sr->sel.remid.pid < 1.0 );
      }
  );

  const Cut kGoodRemID(
    {"sel.remid.pid"},
    [](const caf::StandardRecord* sr)
      { return (sr->sel.remid.pid == 1.0 );
      }
  );

  const Cut kUglyFD = kNumuQuality && kNumuContainFD && kUglyRemID && kNumuCosmicRej;
  const Cut kBadFD  = kNumuQuality && kNumuContainFD && kBadRemID  && kNumuCosmicRej;
  const Cut kGoodFD = kNumuQuality && kNumuContainFD && kGoodRemID && kNumuCosmicRej;

  const Cut kUglyND = kNumuQuality && kNumuContainND && kUglyRemID;
  const Cut kBadND  = kNumuQuality && kNumuContainND && kBadRemID ;
  const Cut kGoodND = kNumuQuality && kNumuContainND && kGoodRemID;

  HistAxis axis( ";Reconstructed neutrino energy (GeV);Events / (0.25 GeV)", kNumuEnergyBinning, kCCE );

  NumuExtrapGenerator extrap_ugly( axis, kUglyFD, kUglyND, kNoShift, kTuftsWeightCC );
  NumuExtrapGenerator extrap_bad ( axis, kBadFD,  kBadND,  kNoShift, kTuftsWeightCC );
  NumuExtrapGenerator extrap_good( axis, kGoodFD, kGoodND, kNoShift, kTuftsWeightCC );

  NumuExtrapGenerator extrap_std ( axis, kNumuFD, kNumuND, kNoShift, kTuftsWeightCC );

  PredictionInterp* prediction_ugly = new PredictionInterp(systs, &calc, extrap_ugly, loaders);
  PredictionInterp* prediction_bad  = new PredictionInterp(systs, &calc, extrap_bad,  loaders);
  PredictionInterp* prediction_good = new PredictionInterp(systs, &calc, extrap_good, loaders);

  PredictionInterp* prediction_std  = new PredictionInterp(systs, &calc, extrap_std,  loaders);

  // GO!
  loaders.Go();

  prediction_ugly->SaveTo(f->mkdir( "pred_ugly" ));
  prediction_bad->SaveTo(f->mkdir( "pred_bad" ));
  prediction_good->SaveTo(f->mkdir( "pred_good" ));

  prediction_std->SaveTo(f->mkdir( "pred_std" ));

  f->Close();

  std::cout << "Predictions saved!" << std::endl;

} // End of function

#endif