syst_table_from_preds.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Prediction/PredictionInterpTemplates.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Systs/NueSystsSecondAna.h"
using namespace ana;

#include "OscLib/IOscCalc.h"
#include "OscLib/OscCalcDumb.h"

#include "TCanvas.h"

#include <iomanip>
#include <map>




// A function to make "key" string in the map unique
const std::string makeKeyUnique (const std::string& syst, const std::string& shift)
{
  std::string key;
  key = syst+"_"+shift;
  std::cout << "makeUniqueKey : " << key << std::endl;
  return key;
}

struct ColumnDef
{
  std::string name;
  Flavors::Flavors_t flav;
  Current::Current_t curr;
  Sign::Sign_t sign;
};

const std::vector<ColumnDef> cols =
  {{"Signal",     Flavors::kNuMuToNuE,  Current::kCC,   Sign::kBoth},
   {"beam NueCC", Flavors::kNuEToNuE,   Current::kCC,   Sign::kBoth},
   {"NC",         Flavors::kAll,        Current::kNC,   Sign::kBoth},
   {"Numu CC",    Flavors::kNuMuToNuMu, Current::kCC,   Sign::kBoth}};


// A function to calculate integral from a spectrum for the specific CVN region you choose
double Integral(const Spectrum& s, const double pot, int cvnregion)
{
  std::unique_ptr<TH1D> h(s.ToTH1(pot));
  
  if(cvnregion == 1)  return h->Integral(1, 9);
  if(cvnregion == 2)  return h->Integral(10, 18);
  if(cvnregion == 3)  return h->Integral(19, 27);
  if(cvnregion == 4)  return h->Integral(28, 36);
  if(cvnregion == 5)  return s.Integral(pot);
}


// A function to calculate the total bkg. for nue appearance
double TotBkg(std::vector <double> integrals)
{
  double totbkg = 0;
  for(auto it: integrals){
    totbkg += it;
  }
  return totbkg;
}


void syst_table_from_preds(int cvnregion, std::string fname = "")
{
  // A function to switch between different PID region
  switch (cvnregion)
    {
    case 1 :
      std::cout << "We are in Region 0.75 < CVN < 0.87" << std::endl;
      break;
    case 2 :
      std::cout << "We are in Region 0.87 < CVN < 0.95" << std::endl;
      break;
    case 3 :
      std::cout << "We are in Region 0.95 < CVN <= 1.0" << std::endl;
      break;
    case 4 :
      std::cout<< "We are in Peripheral Region" <<std::endl;
      break;
    case 5 :
      std::cout<< "We are in CVN Region > 0.75" <<std::endl;
      break;
    default: std::cout << "Run me like <cafe 'syst_table_from_preds.C+(1/2/3/4/5)'> for CVN 1,2,3/peripheral/CVN >.75 region you want to study. \n"
                       << "1 is for 0.75 < CVN1 < 0.87; 2 is for 0.87 < CVN < 0.95; 3 is for CVN > 1.0; 4 is for peripheral region; 5 is for CVN > 0.75. \n"
                       << "Exiting(0)" << std::endl;
      return;
    }
  
  
  osc::OscCalcDumb dumbosc;
  
  const double kPOT= 9.0000e20;   
  
  const int kNumPreds = 2;
  const std::string predNames[kNumPreds] = {"pred_propDecomp", "pred_comboDecomp"};
  
  const int kNumSysts = 5;
  const std::string systNames[kNumSysts] = {"Nominal", "Calibration", "Lightlevel", "Cherenkov", "CalibShape"};
  
  const int kNumShifts = 3;
  const std::string shiftNames[kNumShifts] = {"noShift","plusOne", "minusOne"};
  
  const std::string varName = {"EnergyCVN2D"};
  
  std::unique_ptr <PredictionExtrap> predsPropNom, predsPropShifted;
  
  // map to hold predictions from all systs and then play with it
  std::map <const std::string, std::unique_ptr <PredictionExtrap> >  map_to_preds;
  std::map <const std::string, std::unique_ptr <PredictionExtrap> >::iterator map_it;
  
  
  // I think the preds of 2017 ana will only consider propDecomp but you can loop over "kNumPreds" in future. 
  for(int predIdx = 0; predIdx < 1; ++predIdx)
    {
      const char* predName = predNames[predIdx].c_str();
      std::cout << "predName :" << predName << std::endl;
      
      predsPropNom = LoadFromFile<PredictionExtrap>(fname, TString::Format("%s/%s/%s/nue_pred_%s", predName, systNames[0].c_str(), shiftNames[0].c_str(), varName.c_str()).Data());
      
      for(int systIdx = 1; systIdx < kNumSysts; ++systIdx)
        {
          const char* systName = systNames[systIdx].c_str();
          std::cout << "systName :" << systName << std::endl;
          
          for(int shiftIdx = 1; shiftIdx < kNumShifts; ++shiftIdx)
            {
              const char* shiftName = shiftNames[shiftIdx].c_str();
              std::cout << "shifttName :" << shiftName << std::endl;
              
              if ((systIdx == 3 || systIdx == 4) && shiftIdx == 2) continue;
              
              predsPropShifted = LoadFromFile<PredictionExtrap>(fname, TString::Format("%s/%s/%s/nue_pred_%s", predName, systName, shiftName, varName.c_str()).Data());
              std ::cout << TString::Format("%s/%s/%s/nue_pred_%s", predName, systName, shiftName, varName.c_str()) << std::endl;       
              
              // Insert predictions for all syst. samples in the map.
              map_to_preds.insert(make_pair(makeKeyUnique(systNames[systIdx],shiftNames[shiftIdx]), std::move(predsPropShifted)));
              
            } // end loop of shiftIdx
        } // end loop of systIdx 
    } // end loop of predIdx
  
  
  
  std::cout << "***********************  For Nominal  ************************ " << std::endl;
  
  std::vector <double> nomintegrals;  
  for(ColumnDef col: cols)
    {
      auto spec = predsPropNom->PredictComponent(&dumbosc, col.flav, col.curr, col.sign);
      
      std::cout << col.name << " : " <<  Integral (spec,kPOT,cvnregion) << std::endl;
      
      if ((col.name == "beam NueCC") || (col.name == "NC") || (col.name == "Numu CC"))
        nomintegrals.push_back (Integral(spec,kPOT,cvnregion));

    }
  std::cout << "Total bkg." << " : " << TotBkg(nomintegrals) << std::endl;

  
  
  std::cout << "***********************  For Systematic Sample  ************************ " << std::endl;
  
  for (map_it = map_to_preds.begin(); map_it != map_to_preds.end(); ++map_it)
    {
      auto key  = map_it->first;
      auto pred = map_it->second.get();
      
      std::vector <double> sysintegrals;
      std::cout << "Key/'Name of the Syst.' = " << key << "  :: Value/'Check the valid address' = " << pred <<std::endl;
      
      for(ColumnDef col: cols)
        {
          auto spec = pred->PredictComponent(&dumbosc, col.flav, col.curr, col.sign);
          
          std::cout << col.name << " : "  << Integral (spec,kPOT,cvnregion) << std::endl;
          
          if ((col.name == "beam NueCC") || (col.name == "NC") || (col.name == "Numu CC")) 
            sysintegrals.push_back (Integral(spec,kPOT,cvnregion)); 
        }
      
      std::cout << "Total bkg." << " : " << TotBkg(sysintegrals) << std::endl;
      std::cout << " --------------------------- moving  to  next  systematics -----------------" << std::endl;
      std::cout << " ----------------------------------------------------------------------------" << std::endl;
      
    } // loop for map_it
  
  
}