reduce_the_predfile_nue_numu.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/Vars.h"

#include "TFile.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TGaxis.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TObjectTable.h"
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <cstddef> 

#include "syst_variations.h"

using namespace ana;

void reduce_the_predfile_nue_numu(const bool ReduceNue = false, const TString suffix = "v2", const int thisquantile=-1)
{
  TString outfilename;
  std::vector<TString>  filenames;
  std::vector<TString>  xptype;
  std::vector<SNameDef> systnames;
  if (ReduceNue) {
    // --- If reducing for nue
    outfilename = "pred_nue_reduced_" + suffix + ".root";
    TString folder =  "/nova/ana/nu_e_ana/Ana2017/Predictions/provisional_singles/latest/";
    filenames = { 
     folder+"pred_nxp_EnergyCVN2D_full_nueconcat_other.root",
     folder+"pred_nxp_EnergyCVN2D_full_nueconcat_XSecSyst_part1.root",
     folder+"pred_nxp_EnergyCVN2D_full_nueconcat_XSecSyst_part2.root",
     folder+"pred_xp_prop_EnergyCVN2D_full_nueconcat_other.root",
     folder+"pred_xp_prop_EnergyCVN2D_full_nueconcat_XSecSyst_part1.root",
     folder+"pred_xp_prop_EnergyCVN2D_full_nueconcat_XSecSyst_part2.root",
     folder+"pred_xp_combo_EnergyCVN2D_full_nueconcat_other.root",
     folder+"pred_xp_combo_EnergyCVN2D_full_nueconcat_XSecSyst_part1.root",
     folder+"pred_xp_combo_EnergyCVN2D_full_nueconcat_XSecSyst_part2.root"
    };
    xptype = {"pred_nxp"/*, "pred_xp_prop"*/, "pred_xp_combo"};
    systnames = GetShiftNames2017("allNueReduced");
  } else {
    // --- If reducing for numu.
    outfilename = "pred_numu_reduced_" + suffix + ".root";
    filenames   = { 
      "/nova/ana/nu_mu_ana/Ana2017/Predictions/provisional/singles/fromShiqi/TO_NUMU_XP.root",
      "/nova/ana/nu_mu_ana/Ana2017/Predictions/provisional/pred_nxp_Energy_AllQuantiles_full_numuconcat_numuconcat_Everything_v2.root"
    };
    //filenames   = { "/nova/ana/users/karlwarb/Systematic_Studies/pred_nxp_xp_numu_Energy_AllQuants_full_numudecaf_CalShifts_2sigmaNoTau.root" };
    xptype      = {"pred_nxp","pred_xp_numu"};
    systnames   = GetShiftNames2017("allNumuReduced");
  }
  
  std::cout << "\n\t ReduceNue is " << ReduceNue << " ==> outfilename is " << outfilename << "\n" << std::endl;

  // --- Now open the output file.
  auto outfile = new TFile(outfilename,"recreate");
  // --- Push back the nominal sets to this vector.
  systnames.push_back({"Nominal"   ,{"noShift"}}); // Real nominal
  
  for(auto const & filename:filenames){
    std::cout << "\n Open file " << filename << std::endl;
    auto infile = new TFile (filename,"read");
    for ( auto const & S:systnames ){ 
      
      std::vector<TString> Sigmas = S.sigma_names;

      //HACK HACK
      //if(!ReduceNue && ! S.syst_name.Contains("SumSmallXSecJoint2017")) continue;
      // --- Make sure that Lightlevel has the noshift sample!!
      if (S.syst_name == "Lightlevel") Sigmas = {"plusOne","noShift","minusOne"};
      
      for ( auto const & xp:xptype ){ 
        auto dir = infile->GetDirectory(xp + "_" + S.syst_name);
        if(!dir) continue;
        if(!ReduceNue && S.syst_name.Contains("SumSmallXSecJoint2017") && xp == "pred_xp_numu" && !filename.Contains("Shiqi")) {
          std::cerr << "\nSkipping " << S.syst_name << " from " << filename << std::endl;
          continue;
        }

        auto outdir = outfile->mkdir(xp + "_" + S.syst_name);
        std::cout << "\nCreate " << xp + "_" + S.syst_name << std::endl;

        for ( uint s = 0; s<Sigmas.size(); ++s){
          // === This section is a direct + lazy copy for ReduceNue and ReduceNumu. If you change something in the schema of
          // ===   Nue/Numu please change the same thing in the other one.
          // --- If Reduce Nue
          if (ReduceNue) {
            TString VarNa = "nue_pred_EnergyCVN2D";
            auto dirSigma = dir->GetDirectory(Sigmas[s] + "/" + VarNa );
            if(!dirSigma) {
              std::cerr << "\n\nI can't find " << Sigmas[s] + "/" + VarNa << ". Check the name" << std::endl;
              continue;
              //std::cerr << "Problem " << dirSigma->GetPath() << std::endl;
              //abort();              
            }
            auto temp = outdir->mkdir(Sigmas[s]);
            auto pred = ana::LoadFrom<IPrediction> (dirSigma);
            pred->SaveTo( temp,  VarNa ) ;
            std::cout << "wrote " <<  temp->GetPath() << std::endl;
          } else {
            // --- If Reduce Numu
            std::cerr << "ReduceNue is false so doing NuMu path..." << std::endl;
            // Now loop through quantiles.
            for (uint qu=1; qu<5; ++qu) {

              if(thisquantile > 0 && qu !=thisquantile) continue;

              TString QuantNa = "numu_pred_Energy_Quant" + std::to_string(qu);
              TString DirNa   = Sigmas[s] + "_Quant" + std::to_string(qu);
              
              auto dirSigma = dir->GetDirectory( Sigmas[s] + "/" + QuantNa );
              //HACK HACK
              if(filename.Contains("Shiqi"))  dirSigma = dir->GetDirectory( Sigmas[s] + "_Quant" + std::to_string(qu) + "/" + QuantNa );
              
              if(!dirSigma) {
                std::cerr << "\n\nI can't find " << Sigmas[s] + "/" + QuantNa << ". Check the name" << std::endl;
                continue;
                //abort();
              }

              auto temp = outdir->mkdir( DirNa );
              auto pred = ana::LoadFrom<IPrediction> (dirSigma);
              pred->SaveTo( temp,  QuantNa ) ;
              
              std::cout << "wrote " <<  temp->GetPath() << std::endl;
              delete dirSigma; 
            }
          }
        }
      }//end xp
    }//end systs;
    std::cout << "Done with " << filename << std::endl;
    delete infile;
  }
  std::cout << "Done. Now closing file " << outfilename << std::endl;
  delete outfile;

  std::cout << "Saved to "  << outfilename << std::endl;
}