uncertainty1png.C

Go to the documentation of this file.

 

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"
#include "TFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "TH1.h"
#include "TH2.h"
#include "3FlavorAna/Vars/NumuVars.h"

#include "NDAna/ncpi0_SemiInclusive/ncpi0_1prong/Headers/NCPi0SingleProngVars.h"
#include "CAFAna/Vars/NDNCPi0Xsec.h"

#include "CAFAna/Analysis/Plots.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TColor.h"
#include <stdio.h>
#include <string.h>

using namespace ana;

std::string fname = "BDTGOutput_optimizationTest.root";

TFile fout(fname.c_str(),"RECREATE");

/******************************************************************/

//Declare vectors of all Spectra that we want to generate
std::vector<std::string> parse_as_vector(std::string s)
{
  std::vector<std::string> ret;
  char *pch;
  char *cstr = &s[0u];
  // strtok splits at delimiter, but input has to be char *
  pch = strtok (cstr,",");

  while (pch != NULL)
    {
      //   cout << pch << endl; 
      ret.push_back(pch);
      pch = strtok (NULL, ",");
    }
  return ret;
}

/****************************************************************/

void uncertainty1png()
{
  //Define the datasets and indices to be used.
  const std::string& datasetString = 
    "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1,"
    "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1,"
    "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1,"
    "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1,"
    "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-neg-offset_v1,"
    "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-pos-offset_v1,"
    "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.j_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1";

  const int dataIdx = 0;
  const int nominalIdx = 1;

  /*****************************************************************/

  //Label and load in the spectrum objects for each dataset.
  std::vector<std::string> dataset_labels = {"nominal", "lightdown", "lightup", "ckv", "calibneg", "calibpos", "calibshape"};
  std::vector<ana::SpectrumLoaderBase*> loaders;
  std::vector<std::string> datasets = parse_as_vector(datasetString);

  // Print all the inputs for easy checking.
  std::cout<<"\nDataset names:\n";
  for(auto &id: datasets)
    std::cout<<id<<"\n";

  /*****************************************************************/

  /*Define vectors of spectra of different variables for different loaders.*/
  
  //Total truth information.
  std::vector<std::vector<ana::Spectrum*> > spec_totmc;

  //True signal with fiducial cuts only.
  std::vector<std::vector<ana::Spectrum*> > spec_Nsig_fid_only;

  //True signal that passes preselection cuts.
  std::vector<std::vector<ana::Spectrum*> > spec_NSel_Cuts;

  //Total background.
  std::vector<std::vector<ana::Spectrum*> > spec_Nbkg_Cuts;

  /******************************************************************/

  /*Define vectors of spectra of different variables for
    nominal loader but different flux universes.*/

  //Total flux truth information.
  std::vector<std::vector<ana::Spectrum*> > spec_totmc_flux;

  //True signal with fiducial cuts only. 
  std::vector<std::vector<ana::Spectrum*> > spec_Nsig_fid_only_flux;
  
  //True signal that passes preselection cuts.
  std::vector<std::vector<ana::Spectrum*> > spec_NSel_Cuts_flux; 

  //True signal that passes all cuts (Total Background)
  std::vector<std::vector<ana::Spectrum*> > spec_Nbkg_Cuts_flux;

  /*******************************************************************/

  /*Define vectors of spectra of different variables
    for nominal loader, but different xsec universes.*/

  //Total MC Preselection.
  std::vector<std::vector<ana::Spectrum*> > spec_totmc_xsec; 

  //True signal with fiducial cuts only.
  std::vector<std::vector<ana::Spectrum*> > spec_Nsig_fid_only_xsec;
  
  //Signal that passes preselection cuts.
  std::vector<std::vector<ana::Spectrum*> > spec_NSel_Cuts_xsec; 

  //True signal that passes all cuts (Total Background)
  std::vector<std::vector<ana::Spectrum*> > spec_Nbkg_Cuts_xsec;    

  /********************************************************************/

  //Set spill cuts on all spectra by looping over all loaders.
  for(unsigned int idataset = 0; idataset < datasets.size(); ++idataset)
    {
      SpectrumLoaderBase*  loader;
      loader = new SpectrumLoader(datasets[idataset]);
      loader->SetSpillCut(kStandardDQCuts && kTightBeamQualityCuts);
      loaders.push_back(loader);
    }

  const int nvar = 1;

  const Binning BDTBins = Binning::Simple(100,-1,1);

  HistAxis taxes[nvar] = 
    {
      HistAxis ("Single Prong BDTG Output", BDTBins, kNCPi0BDTID)
    };

  /************************************************************************/

  //Begin filling spectra with multiverse parameters.
  std::vector<std::string> taxes_labels = {"BDTG Output"};
  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
  GenieMultiverseParameters verse(100,  "/nova/app/users/rbowles/tag_releaseS18-04-23/CAFAna/XSec/Utilities/knob_config.txt");
  std::vector<std::map<const ISyst*, double> > genie_shifts = verse.ShiftTables();
  
  //Loop over loaders and variables to fill spectra.
  for(unsigned int iload = 0; iload < loaders.size(); ++iload)
    {
      Var weight = kXSecCVWgt2017 * kPPFXFluxCVWgt;
      // SpillTruthVar weightST = kXSecCVWgt2017ST * kPPFXFluxCVWgtST;
      // one vector per loader for all the vars

      //Total Preselection.
      std::vector<ana::Spectrum*> vspec_totmc;

      //Denominator: True signal fiducial.
      std::vector<ana::Spectrum*>  vspec_Nsig_fid_only;

      //Signal numerator: Signal cuts.
      std::vector<ana::Spectrum*>  vspec_NSel_Cuts; 
 
      //Preselection + !Signal
      std::vector<ana::Spectrum*>  vspec_Nbkg_Cuts;
    
      //Loop over axes.
      for(int ivar = 0; ivar < nvar; ivar++)
        {
          //Total Preselection including RemID.
          vspec_totmc.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection,kNoShift,weight));

          //Total signal.
          vspec_NSel_Cuts.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Signal,kNoShift,weight));

          //True Signal with only true fiducial (Fiducial, pi0, NCurrent)
          vspec_Nsig_fid_only.push_back(new Spectrum(*loaders[iload],taxes[ivar], kTrueFiducial && Signal,kNoShift,weight));

          //Total Background
          vspec_Nbkg_Cuts.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Bkg,kNoShift,weight));


          //Only build universes for the nominal loader.
          if(iload == 0)
            {
              /*Define flux spectra.*/

              //Total flux.
              std::vector<ana::Spectrum*>  vspec_totmc_flux;

              //Total signal flux.
              std::vector<ana::Spectrum*>  vspec_Nsig_fid_only_flux; 

              //Signal that passes preselection.
              std::vector<ana::Spectrum*>  vspec_NSel_Cuts_flux;

              //Total background.
              std::vector<ana::Spectrum*>  vspec_Nbkg_Cuts_flux;

              /*****************************************************************/
              /*Define cross section spectra.*/
              
              //Total cross section.
              std::vector<ana::Spectrum*>  vspec_totmc_xsec;

              //Total signal cross section.
              std::vector<ana::Spectrum*>  vspec_Nsig_fid_only_xsec; 

              //Signal that passes preselection.
              std::vector<ana::Spectrum*>  vspec_NSel_Cuts_xsec;

              //Total background.
              std::vector<ana::Spectrum*>  vspec_Nbkg_Cuts_xsec;

              //Loop over all 100 universes.
              for(int iuniv = 0; iuniv < 100; iuniv++)
                {
                  /**********Fill spectra from ppfx universes.***********/
                  vspec_totmc_flux.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection,kNoShift,ppfx_univ[iuniv]));

                  vspec_NSel_Cuts_flux.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Signal,kNoShift,ppfx_univ[iuniv]));
                  
                  vspec_Nsig_fid_only_flux.push_back(new Spectrum(*loaders[iload],taxes[ivar],kTrueFiducial && Signal,kNoShift,ppfx_univ[iuniv]));        
                  
                  vspec_Nbkg_Cuts_flux.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Bkg,kNoShift,ppfx_univ[iuniv]));

                  /**********Fill spectra from genie universes.**********/
                  vspec_totmc_xsec.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection,genie_shifts[iuniv],weight));
                  
                  vspec_NSel_Cuts_xsec.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Signal,genie_shifts[iuniv],weight) );
                  
                  vspec_Nsig_fid_only_xsec.push_back(new Spectrum(*loaders[iload],taxes[ivar],kTrueFiducial && Signal,genie_shifts[iuniv],weight));
                  
                  vspec_Nbkg_Cuts_xsec.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Bkg,genie_shifts[iuniv],weight)); 
                }
              spec_totmc_xsec.push_back(vspec_totmc_xsec);
              spec_Nsig_fid_only_xsec.push_back(vspec_Nsig_fid_only_xsec);
              spec_NSel_Cuts_xsec.push_back(vspec_NSel_Cuts_xsec);
              spec_Nbkg_Cuts_xsec.push_back(vspec_Nbkg_Cuts_xsec);
              
              spec_totmc_flux.push_back(vspec_totmc_flux);
              spec_Nsig_fid_only_flux.push_back(vspec_Nsig_fid_only_flux);
              spec_NSel_Cuts_flux.push_back(vspec_NSel_Cuts_flux);
              spec_Nbkg_Cuts_flux.push_back(vspec_Nbkg_Cuts_flux);
            }
        }
      spec_totmc.push_back(vspec_totmc);
      spec_Nsig_fid_only.push_back(vspec_Nsig_fid_only);
      spec_NSel_Cuts.push_back(vspec_NSel_Cuts);
      spec_Nbkg_Cuts.push_back(vspec_Nbkg_Cuts); 
    }
  /************************************************************************/
   
  for(unsigned int iload = 0; iload < loaders.size(); ++iload)
    {
      loaders[iload]->Go();
    }

  //Loop over loaders
  for(unsigned int iload = 0; iload < loaders.size(); ++iload)
    {
      TDirectory* d = fout.mkdir(dataset_labels[iload].c_str());      
      
      //Loop over variables.
      for(int ivar = 0; ivar < nvar; ++ivar)
        {
          std::cout<<"loader "<<dataset_labels[iload]<<" var "<<taxes_labels[ivar]<<"\n";
          spec_totmc[iload][ivar]->SaveTo(d->mkdir( ("totmc_"+taxes_labels[ivar]).c_str()));
          spec_Nsig_fid_only[iload][ivar]->SaveTo(d->mkdir( ("Nsig_fid_only_"+taxes_labels[ivar]).c_str()));
          spec_NSel_Cuts[iload][ivar]->SaveTo(d->mkdir( ("NSel_Cuts_"+taxes_labels[ivar]).c_str()));
          spec_Nbkg_Cuts[iload][ivar]->SaveTo(d->mkdir( ("Nbkg_Cuts_"+taxes_labels[ivar]).c_str()));
      
          if(iload == 0)
            {
              TDirectory* dftmc = fout.mkdir(("flux_totmc_"+taxes_labels[ivar]).c_str()); 
              TDirectory* dft = fout.mkdir(("flux_Nsig_fid_only_"+taxes_labels[ivar]).c_str());     
              TDirectory* dfp = fout.mkdir(("flux_NSel_Cuts_"+taxes_labels[ivar]).c_str());      
              TDirectory* dfb = fout.mkdir(("flux_Nbkg_Cuts_"+taxes_labels[ivar]).c_str());      
              
              TDirectory* dxtmc = fout.mkdir(("xsec_totmc_"+taxes_labels[ivar]).c_str()); 
              TDirectory* dxt = fout.mkdir(("xsec_Nsig_fid_only_"+taxes_labels[ivar]).c_str());      
              TDirectory* dxp = fout.mkdir(("xsec_NSel_Cuts_"+taxes_labels[ivar]).c_str());      
              TDirectory* dxb = fout.mkdir(("xsec_Nbkg_Cuts_"+taxes_labels[ivar]).c_str());      
              
              for(int iuniv = 0; iuniv < 100; iuniv++)
                {
                  spec_totmc_flux[ivar][iuniv]->SaveTo(dftmc->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
                  spec_Nsig_fid_only_flux[ivar][iuniv]->SaveTo(dft->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
                  spec_NSel_Cuts_flux[ivar][iuniv]->SaveTo(dfp->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
                  spec_Nbkg_Cuts_flux[ivar][iuniv]->SaveTo(dfb->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
                  spec_totmc_xsec[ivar][iuniv]->SaveTo(dxtmc->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
                  spec_Nsig_fid_only_xsec[ivar][iuniv]->SaveTo(dxt->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
                  spec_NSel_Cuts_xsec[ivar][iuniv]->SaveTo(dxp->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
                  spec_Nbkg_Cuts_xsec[ivar][iuniv]->SaveTo(dxb->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
                }
            }     
        }
    }
  std::cout << "\nOutput saved to BDTGOutput_optimization.root" <<std::endl;
}