MRDiFStudy_RHC_Step1.C

Go to the documentation of this file.

//macro to get spectra used to get weights for the reweighting
//process in Step 2 

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/MultiVar.h"

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"

#include "3FlavorAna/Ana2020/MRstudies/MRDiFCutsAndVars.h"
//#include "MRDiFCutsAndVars.h"

#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVars.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include <cassert>

#include "TFile.h"
#include "TCanvas.h"
#include "TH2.h"
#include "TH1.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TArrayD.h"
#include "TTree.h"
#include "TLeaf.h"
#include "stdio.h"
#include <sstream>
#include <string>
#include <iostream>
#include <fstream>
#include <string>


using namespace ana;
using namespace std;


void MRDiFStudy_RHC_Step1()
{

  const std::string filenameCosmicsData = "prod_mrdifcaf_R19-11-18-prod5reco.i_fd_cosmic_rhc_full_v1_goodruns";
  const std::string filenameCRYMC = "prod_mrdifcaf_R19-11-18-prod5reco.m_fd_cry_rhc_full_v3";
  const std::string filenameGENIEFluxswap = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_rhc_nova_v08_full_v1";

  const Var kShwAngle([](const caf::SRProxy* sr)
                      {
                        double maxe = -99.0;
                        int maxi = 999.0;
                        double cosz = -1.0;

                        if(!sr->vtx.elastic.IsValid) return -99999.0;
                        if(sr->vtx.elastic.fuzzyk.png.size() < 1) return -99999.0;

                        for(unsigned int i=0;i<sr->vtx.elastic.fuzzyk.png.size();i++){
                          if(maxe < sr->vtx.elastic.fuzzyk.png[i].shwlid.shwE){
                            maxe = sr->vtx.elastic.fuzzyk.png[i].shwlid.shwE;
                            maxi = i;
                            cosz = sr->vtx.elastic.fuzzyk.png[i].shwlid.dir.z;
                          }
                        }
                        return cosz;
                      });

  const Var kShwE( [](const caf::SRProxy* sr)
                   {
                     double maxe = -99.0;
                     int maxi = 999;
                     
                     if(!sr->vtx.elastic.IsValid) return -99999.0;
                     if(sr->vtx.elastic.fuzzyk.png.size() < 1) return -99999.0;
                     
                     for(unsigned int i=0;i<sr->vtx.elastic.fuzzyk.png.size();i++){
                       if(maxe < sr->vtx.elastic.fuzzyk.png[i].shwlid.shwE){
                         maxe = sr->vtx.elastic.fuzzyk.png[i].shwlid.shwE;
                         maxi = i;
                       }
                     }
                     return maxe;
                     
                   });


  SpectrumLoader loaderCosmicsData(filenameCosmicsData, kCosmic);
  SpectrumLoader loaderCRYMC(filenameCRYMC, kCosmic);
  SpectrumLoader loaderGENIEFluxswap(filenameGENIEFluxswap);
  loaderGENIEFluxswap.SetSpillCut(kStandardSpillCuts);


  const Binning ShwAnglebin= Binning::Simple(60, 0.0, 1.0);
  const Binning ShwEbin = Binning::Simple(120, 0, 5.0);

 //1D spectrums CorePresel  Cuts
  Spectrum sShwAngleCorePreselCosmicsData("Cos(#theta_{Z})", ShwAnglebin,loaderCosmicsData,kShwAngle,kDiFPreselectionCut&&kInCosmicTimingWindow_FD_MR);
  Spectrum sShwAngleCorePreselCRYMC("Cos(#theta_{Z})", ShwAnglebin,loaderCRYMC,kShwAngle,kDiFPreselectionCut&&kInCosmicTimingWindow_FD_MR);
  Spectrum sShwAngleCorePreselGENIEFluxswap("Cos(#theta_{Z})", ShwAnglebin,loaderGENIEFluxswap,kShwAngle,kDiFPreselectionCut);
  /*
  Spectrum sShwECorePreselCosmicsData("Cos(#theta_{Z})", ShwEbin,loaderCosmicsData,kShwE,kDiFPreselectionCut&&kInCosmicTimingWindow_FD_MR);
  Spectrum sShwECorePreselCRYMC("Cos(#theta_{Z})", ShwEbin,loaderCRYMC,kShwE,kDiFPreselectionCut&&kInCosmicTimingWindow_FD_MR);
  Spectrum sShwECorePreselGENIEFluxswap("Cos(#theta_{Z})", ShwEbin,loaderGENIEFluxswap,kShwE,kDiFPreselectionCut);
  */

  //1D spectrums Core or Peripheral Cuts
  Spectrum sShwAngleCoreCosmicsData("Cos(#theta_{Z})", ShwAnglebin,loaderCosmicsData,kShwAngle,kDiFSelectionCut&&kInCosmicTimingWindow_FD_MR);
  Spectrum sShwAngleCoreCRYMC("Cos(#theta_{Z})", ShwAnglebin,loaderCRYMC,kShwAngle,kDiFSelectionCut&&kInCosmicTimingWindow_FD_MR);
  Spectrum sShwAngleCoreGENIEFluxswap("Cos(#theta_{Z})", ShwAnglebin,loaderGENIEFluxswap,kShwAngle,kDiFSelectionCut);
  /*
  Spectrum sShwECoreCosmicsData("Cos(#theta_{Z})", ShwEbin,loaderCosmicsData,kShwE,kDiFSelectionCut&&kInCosmicTimingWindow_FD_MR);
  Spectrum sShwECoreCRYMC("Cos(#theta_{Z})", ShwEbin,loaderCRYMC,kShwE,kDiFSelectionCut&&kInCosmicTimingWindow_FD_MR);
  Spectrum sShwECoreGENIEFluxswap("Cos(#theta_{Z})", ShwEbin,loaderGENIEFluxswap,kShwE,kDiFSelectionCut);
  */


  loaderCosmicsData.Go();
  loaderCRYMC.Go();
  loaderGENIEFluxswap.Go();

  TFile *file = new TFile("MRDiF_FD_RHC_Step1.root", "recreate");

  sShwAngleCorePreselCosmicsData.SaveTo(file, "sShwAngleCorePreselCosmicsData");
  sShwAngleCorePreselCRYMC.SaveTo(file, "sShwAngleCorePreselCRYMC");
  sShwAngleCorePreselGENIEFluxswap.SaveTo(file, "sShwAngleCorePreselGENIEFluxswap");
  /*
  sShwECorePreselCosmicsData.SaveTo(file, "sShwECorePreselCosmicsData");
  sShwECorePreselCRYMC.SaveTo(file, "sShwECorePreselCRYMC");
  sShwECorePreselGENIEFluxswap.SaveTo(file, "sShwECorePreselGENIEFluxswap");
  */
  //1D spectrums Core Cuts
  sShwAngleCoreCosmicsData.SaveTo(file, "sShwAngleCoreCosmicsData");
  sShwAngleCoreCRYMC.SaveTo(file, "sShwAngleCoreCRYMC");
  sShwAngleCoreGENIEFluxswap.SaveTo(file, "sShwAngleCoreGENIEFluxswap");
  /*
  sShwECoreCosmicsData.SaveTo(file, "sShwECoreCosmicsData");
  sShwECoreCRYMC.SaveTo(file, "sShwECoreCRYMC");
  sShwECoreGENIEFluxswap.SaveTo(file, "sShwECoreGENIEFluxswap");
  */
  file->Close();
  return;

}