MRDiFStudy_RHC_Step2.C

Go to the documentation of this file.

#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;

//need global histos that'll turn into the angle weights
TH1D* hPreselGENIE_forCRY=0;
TH1D* hPreselGENIE_forData=0;
//TH1D* hPreselCRY=0;
//TH1D* hPreselData=0;



void MRDiFStudy_RHC_Step2()
{

  //datasets to run over
  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";
  
  //input step 1 file - has the angle distributions used to make the weights
  //TFile *f1 = new TFile("3FlavorAna/Ana2020/MRstudies/MRDiF_FD_RHC_Step1.root", "READ");
  TFile *f1 = TFile::Open(pnfs2xrootd("/pnfs/nova/scratch/users/amhall/DiF/RHC/MRDiF_FD_RHC_Step1.root").c_str());
  f1->cd();

  //pull in those angle spectra
  auto sShwAnglePreselCosmicsData=*Spectrum::LoadFrom(f1, "sShwAngleCorePreselCosmicsData");
  auto sShwAnglePreselCRYMC=*Spectrum::LoadFrom(f1, "sShwAngleCorePreselCRYMC");
  auto sShwAnglePreselGENIEFluxswap=*Spectrum::LoadFrom(f1, "sShwAngleCorePreselGENIEFluxswap");

  //get applicable exposures for normalizing
  double mylivetime=sShwAnglePreselCosmicsData.Livetime();
  double myPOT=sShwAnglePreselGENIEFluxswap.POT();

  hPreselGENIE_forCRY=sShwAnglePreselGENIEFluxswap.ToTH1(myPOT);
  hPreselGENIE_forData=sShwAnglePreselGENIEFluxswap.ToTH1(myPOT);
  
  TH1D* hPreselCRY=sShwAnglePreselCRYMC.ToTH1(mylivetime,kLivetime);
  TH1D* hPreselData=sShwAnglePreselCosmicsData.ToTH1(mylivetime,kLivetime);
  
  hPreselGENIE_forCRY->Rebin();
  hPreselGENIE_forData->Rebin();
  hPreselCRY->Rebin();
  hPreselData->Rebin();
  
  hPreselGENIE_forCRY->Scale(hPreselCRY->Integral()/hPreselGENIE_forCRY->Integral());
  hPreselGENIE_forCRY->Divide(hPreselCRY);
  
  hPreselGENIE_forData->Scale(hPreselData->Integral()/hPreselGENIE_forData->Integral());
  hPreselGENIE_forData->Divide(hPreselData);
  
  //now we have to define some variables before we can define the weights properly
  
  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;
                     
                   });

  //now we can define the weights as Vars
  //CRY  weights
  const Var kCRYMC1Dweights( [kShwAngle, kShwE](const caf::SRProxy* sr)
                                          {
                                            double myweight=1.0;
                                            double myWgtVar;
                                            myWgtVar = kShwAngle(sr);
                                            int bin=hPreselGENIE_forCRY->FindBin(myWgtVar);
                                            myweight=hPreselGENIE_forCRY->GetBinContent(bin);
                                            if(myweight>=0)return myweight;
                                            else return 1.0;
                                          });


  //Data weights
  const Var kData1Dweights( [kShwAngle, kShwE](const caf::SRProxy* sr)
                                    {
                                      double myweight=1.0;
                                      double myWgtVar;
                                      myWgtVar = kShwAngle(sr);
                                      int bin=hPreselGENIE_forData->FindBin(myWgtVar);
                                      myweight=hPreselGENIE_forData->GetBinContent(bin);
                                      if(myweight>=0)return myweight;
                                      else return 1.0;
                                    });


  //-----------------------------------------------------------------------------------

  //now the rest of the desired variables we need to define

const Var kShwStartX(
                              [](const caf::SRProxy *sr)
                              {

                                double maxe = -99.0;
                                int maxi = 999;
                                double startx = -99999;

                                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;
                                    startx = sr->vtx.elastic.fuzzyk.png[i].shwlid.start.x;
                                  }
                                }

                                return startx;

                              });



 const Var kShwStartY(
                     [](const caf::SRProxy *sr)
                     {

                       double maxe = -99.0;
                       int maxi = 999;
                       double starty = -99999;

                       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;
                           starty = sr->vtx.elastic.fuzzyk.png[i].shwlid.start.y;
                         }
                       }

                       return starty;

                     });



 const Var kShwStartZ(
                              [](const caf::SRProxy *sr)
                              {

                                double maxe = -99.0;
                                int maxi = 999;
                                double startz = -99999;

                                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;
                                    startz = sr->vtx.elastic.fuzzyk.png[i].shwlid.start.z;
                                  }
                                }

                                return startz;

                              });


  const Var kShwWidth(
                  [](const caf::SRProxy* sr)
                  {

                    double maxe = -99.0;
                    int maxi = 999;
                    double width = -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;
                        width = sr->vtx.elastic.fuzzyk.png[i].shwlid.width;
                      }
                    }

                    return width;

                  });




  const Var kShwLength(
                  [](const caf::SRProxy* sr)
                  {

                    double maxe = -99.0;
                    int maxi = 999;
                    double length = -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;
                        length = sr->vtx.elastic.fuzzyk.png[i].shwlid.len;
                      }
                    }

                    return length;

                  });


  const Var kEPerHit = kCaloE / kNHit;

  const Var kShwNHit(
                  [](const caf::SRProxy* sr)
                  {

                    double maxe = -99.0;
                    int maxi = 999;
                    double nhit = -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;
                        nhit = sr->vtx.elastic.fuzzyk.png[i].shwlid.nhit;
                      }
                    }

                    return nhit;

                  }); 

const Var kVtxX([](const caf::SRProxy* sr) 
                 { 
                   if (!sr->vtx.elastic.IsValid || !(sr->vtx.elastic.fuzzyk.png.size() >0)) return -99999.0;
                   double vtx = sr->vtx.elastic.vtx.x;
                   return vtx;
                 });

 const Var kVtxY([](const caf::SRProxy* sr) 
                 { 
                   if (!sr->vtx.elastic.IsValid || !(sr->vtx.elastic.fuzzyk.png.size() >0)) return -99999.0;
                   double vtx =  sr->vtx.elastic.vtx.y;
                   return vtx;
                 });

 const Var kVtxZ([](const caf::SRProxy* sr) 
                 { 
                   if (!sr->vtx.elastic.IsValid || !(sr->vtx.elastic.fuzzyk.png.size() >0)) return -99999.0;
                   double vtx =  sr->vtx.elastic.vtx.z;
                   return vtx;
                 });

 //cut to ignore angle bins where we had data points, but not MC. 
  const Cut kShwAngleCut = (kShwAngle >= 0.5 && kShwAngle < 0.9667);

   //Binning
  const Binning ShwEbin= Binning::Simple(120, 0.0, 5.0);
  const Binning ShwWidthbin= Binning::Simple(100, 0.0, 10.0);
  const Binning ShwLengthbin= Binning::Simple(100, 0.0, 1000.0);
  const Binning ShwAnglebin= Binning::Simple(30, 0.0, 1.0);
  const Binning ShwStartXbin= Binning::Simple(800, -800, 800);
  const Binning ShwStartYbin= Binning::Simple(800, -800, 800);
  const Binning ShwStartZbin= Binning::Simple(3000,0,6000);
  const Binning ShwCVNebin= Binning::Simple(120, 0.0, 1.000);
  const Binning NHitbin= Binning::Simple(120,0.0,150);
  const Binning EPerHitbin= Binning::Simple(200, 0, 0.5);
  const Binning distVerticalbin = Binning::Simple(3260, 0, 1600);
  const Binning distHorizontalbin= Binning::Simple(320, 0, 1600);
  const Binning distLengthbin= Binning::Simple(3000, 0, 6000);


  //some structure to make making the sprectra easier
  struct Plot{const char* name; Var var; const Binning binning;};

  const std::vector<Plot> plots = {
    {"Shower Angle (cos_theta_z)",  kShwAngle,      ShwAnglebin},
    {"Number of Hits in Slice",     kNHit,          NHitbin },
    {"Number of Hits in Shower",    kShwNHit,       NHitbin},
    {"Energy in Shower (GeV)",      kShwE,          ShwEbin },
    {"Energy per Hit (CaloE/NHit)", kEPerHit,       EPerHitbin},
    {"CVNe - loosptpt",             kCVNe_looseptp, ShwCVNebin},
    {"Shower Width",                kShwWidth,      ShwWidthbin},
    {"Shower Length",               kShwLength,     ShwLengthbin},
    {"Shower Start X",              kShwStartX,     ShwStartXbin},
    {"Shower Start Y",              kShwStartY,     ShwStartYbin},
    {"Shower Start Z",              kShwStartZ,     ShwStartZbin},
    {"Vtx X",                       kVtxX,          ShwStartXbin},
    {"Vtx Y",                       kVtxY,          ShwStartYbin},
    {"Vtx Z",                       kVtxZ,          ShwStartZbin},
    {"Distance to Top",             kdistTop,       distVerticalbin},
    {"Distance to Bottom",          kdistBottom,    distVerticalbin},
    {"Distance to East",            kdistEast,      distHorizontalbin},
    {"Distance to West",            kdistWest,      distHorizontalbin},
    {"Distance to Front",           kdistFront,     distLengthbin},
    {"Distance to Back",            kdistBack,      distLengthbin}
  };

  //and define our loaders
  SpectrumLoader loaderCosmicsData(filenameCosmicsData, kCosmic);
  SpectrumLoader loaderCRYMC(filenameCRYMC, kCosmic);
  SpectrumLoader loaderGENIE(filenameGENIEFluxswap);
  loaderGENIE.SetSpillCut(kStandardSpillCuts);

  std::vector<Spectrum*> sDataPresel, sCRYMCPresel, sGENIEPresel, sData, sCRYMC, sGENIE;
  std::vector<Spectrum*> sDataPreselWtd, sCRYMCPreselWtd, sDataWtd, sCRYMCWtd;

  for(const Plot& p: plots){

    sDataPresel.push_back(    new Spectrum(p.name, p.binning, loaderCosmicsData, p.var, kDiFPreselectionCut && kShwAngleCut && kInCosmicTimingWindow_FD_MR )); 
    sData.push_back(          new Spectrum(p.name, p.binning, loaderCosmicsData, p.var, kDiFSelectionCut    && kShwAngleCut && kInCosmicTimingWindow_FD_MR )); 
    sDataPreselWtd.push_back( new Spectrum(p.name, p.binning, loaderCosmicsData, p.var, kDiFPreselectionCut && kShwAngleCut && kInCosmicTimingWindow_FD_MR, kNoShift, kData1Dweights)); 
    sDataWtd.push_back(       new Spectrum(p.name, p.binning, loaderCosmicsData, p.var, kDiFSelectionCut    && kShwAngleCut && kInCosmicTimingWindow_FD_MR, kNoShift, kData1Dweights)); 

    sCRYMCPresel.push_back(   new Spectrum(p.name, p.binning, loaderCRYMC,       p.var, kDiFPreselectionCut && kShwAngleCut && kInCosmicTimingWindow_FD_MR )); 
    sCRYMC.push_back(         new Spectrum(p.name, p.binning, loaderCRYMC,       p.var, kDiFSelectionCut    && kShwAngleCut && kInCosmicTimingWindow_FD_MR )); 
    sCRYMCPreselWtd.push_back(new Spectrum(p.name, p.binning, loaderCRYMC,       p.var, kDiFPreselectionCut && kShwAngleCut && kInCosmicTimingWindow_FD_MR, kNoShift, kCRYMC1Dweights)); 
    sCRYMCWtd.push_back(      new Spectrum(p.name, p.binning, loaderCRYMC,       p.var, kDiFSelectionCut    && kShwAngleCut && kInCosmicTimingWindow_FD_MR, kNoShift, kCRYMC1Dweights)); 

    sGENIEPresel.push_back(   new Spectrum(p.name, p.binning, loaderGENIE,       p.var, kDiFPreselectionCut && kShwAngleCut )); 
    sGENIE.push_back(         new Spectrum(p.name, p.binning, loaderGENIE,       p.var, kDiFSelectionCut    && kShwAngleCut )); 

  }

  loaderCosmicsData.Go();
  loaderCRYMC.Go();
  loaderGENIE.Go();


  //and save the spectra to a file

  TFile *outFile = new TFile("MRDiF_FD_RHC_Step2_v2_clean.root", "recreate");
  char cname[128];

  for (unsigned int i=0; i<sDataPresel.size(); i++){
    sprintf(cname, "sDataPresel_V%d", i);
    sDataPresel[i]->SaveTo(outFile, cname);
    sprintf(cname, "sData_V%d", i);
    sData[i]->SaveTo(outFile, cname);
    sprintf(cname, "sDataPreselWtd_V%d", i);
    sDataPreselWtd[i]->SaveTo(outFile, cname);
    sprintf(cname, "sDataWtd_V%d", i);
    sDataWtd[i]->SaveTo(outFile, cname);

    sprintf(cname, "sCRYMCPresel_V%d", i);
    sCRYMCPresel[i]->SaveTo(outFile, cname);
    sprintf(cname, "sCRYMC_V%d", i);
    sCRYMC[i]->SaveTo(outFile, cname);
    sprintf(cname, "sCRYMCPreselWtd_V%d", i);
    sCRYMCPreselWtd[i]->SaveTo(outFile, cname);
    sprintf(cname, "sCRYMCWtd_V%d", i);
    sCRYMCWtd[i]->SaveTo(outFile, cname);

    sprintf(cname, "sGENIEPresel_V%d", i);
    sGENIEPresel[i]->SaveTo(outFile, cname);
    sprintf(cname, "sGENIE_V%d", i);
    sGENIE[i]->SaveTo(outFile, cname);
  }

}