NuMu2020_MakePIDPlots.C

Go to the documentation of this file.

/*
  Karl W: Written for the 2020 PID NuMu Optimisation.
  To run:
     root
     .x NuMu2020_MakePIDPlots.C( <bool for FD> , 
                                 <bool for Horn Current>, 
                                 <enum for Quantile to use>,   { 0 = AllQuants, 123 = Quants 1,2,3 and 4 = Just Quant 4 }
                                 <bool for ReMId/ProngCVN>,    Defaulted to true
                                 <bool for LoosePTP/OldPresel> Defaulted to true )
  Will run over either the ND or the FD.
  Will run over either FHC or RHC data. 
  Will fetch either ReMId or Prong CVN Scores
  Will fetch either LoosePTP or OldPresel trained CVN Scores.
*/

#include "TFile.h"
#include "TH1.h"
#include "TH3.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"
#include "TSystem.h"

#include <iostream>
#include <iomanip> 
#include <stdio.h>
#include <string> 

#include <TStyle.h>
#include <TCanvas.h>

//#include "Utilities/rootlogon.C"

using namespace std;

// ---------------------------------------------------------------------------------------------------------------------
// --- What are the scales to use?
// The official POTs and Livetimes are taken from Analysis/Exposures.h, and are the Ana19 values.
// The production POTs and Livetimes are taken from a Spectra ran over the entire dataset in NuMu2020_BasicPlots_*.C. 
const double FHC_POT_Scale = 9.47964e20 / 1.7980696e24;
const double FHC_Liv_Scale = 438.1584   / 57385.145;

const double RHC_POT_Scale = 12.3319e20 / 4.5758019e24;
const double RHC_Liv_Scale = 317.0466   / 43780.253;

double MyPOTScale = 0;
double MyLivScale = 0;
// ---------------------------------------------------------------------------------------------------------------------
const double MinTrPID  = 0.3, MaxTrPID  = 0.95; const int nbinTrPID  = 65;
const double MinCVN    = 0.4, MaxCVN    = 0.95; const int nbinCVN    = 55;
const double MinCosRej = 0.4, MaxCosRej = 0.65; const int nbinCosRej = 25;
// ---------------------------------------------------------------------------------------------------------------------
const bool bPrintInfo = false; // Dirty, but easy to toggle how much info I want...

// ---------------------------------------------------------------------------------------------------------------------
//                                A struct to store key FOM values.
struct MyFOMParams {
  // Numbers for the full dist.
  double FOM_Full;
  double TotalSig_Full;
  double TotalBkg_Full;
  // Numbers for the dip region.
  double FOM_Dip;
  double TotalSig_Dip;
  double TotalBkg_Dip;
};
// ---------------------------------------------------------------------------------------------------------------------

// ---------------------------------------------------------------------------------------------------------------------
//                                Declare any functions that I'm going to use later
// ---------------------------------------------------------------------------------------------------------------------
void PrintInfo( int Entry, int TotEnt, double IsNuMuCC, double TrueNuE, double RecoNuE, double EventWgt, double CosVeto,
                double CVNMuonLooPTP, double ReMId, double CVNProngMuon, double CosRej2020 );
// ---------------------------------------------------------------------------------------------------------------------
void CalcFOM( TH3D* Sig, TH3D* Bkg, TH3D* Data, int bintrpid, int bincvn, int bincos, bool IsFD, 
              double &FOM, double &TotSig, double &TotBkg );
// ---------------------------------------------------------------------------------------------------------------------
void PrintMyMap( std::map< std::string, MyFOMParams > MainMap, std::string Key );
// ---------------------------------------------------------------------------------------------------------------------
bool PassCutPoint( double PIDScore, double PIDCut, double CVNScore, double CVNCut, 
                   int Quant, int CosVeto, std::vector<int> cQuantVec, double RecoE, bool isDip );
// ---------------------------------------------------------------------------------------------------------------------
void FillMyPlots ( double CVN, double ReMId, double CosRej, int Quant, TH1D* Hists[9], double RecoE, double Wgt );
// ---------------------------------------------------------------------------------------------------------------------
// --- Make the axis for my TH1Ds if making Reco Nu E plots.
const int HistBins = 19;
double RecoEBins[HistBins+1] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 };

// ---------------------------------------------------------------------------------------------------------------------
//                                        The main class
// ---------------------------------------------------------------------------------------------------------------------
void NuMu2020_MakePIDPlots( bool IsFD, bool IsFHC, int QuantEnum, bool IsReMId=true, bool IsLoosePTP=true ) {
  // --- Quant Enum
  // 0   - all quants together.
  // 123 - quantiles 1, 2 and 3 together.
  // 4   - quantile 4 by itself.

  // Remove stat boxes.
  gStyle->SetOptStat(0);

  // FD or ND? This affects whether the "data" is real data, or cosmic data.
  std::string  sIsFD = "FD", sData = "Cosm", sMC = "MC";
  if (!IsFD) { sIsFD = "ND", sData = "Data"; }
  // FHC or RHC? 
  std::string sIsFHC = "FHC";
  if (!IsFHC) sIsFHC = "RHC";
  // Am I using ReMId or Prong CVN?
  std::string   sTrPID = "ReMId";
  if (!IsReMId) sTrPID = "ProngCVN";
  // Am I using the LoosePTP CVN or the OldPresel CVN?
  std::string      sCVN = "LoosePTP";
  if (!IsLoosePTP){sCVN = "OldPresel"; std::cout << "\nHave you reenabled OldPresel CVN?\n" << std::endl; return; }
  // What is my quantile enum?
  std::string sQuant = "";
  std::vector<int> cQuantVec;
  if      (QuantEnum == 0  ) { sQuant = "AllQuant"; cQuantVec = {1,2,3,4}; }
  else if (QuantEnum == 123) { sQuant = "Quant123"; cQuantVec = {1,2,3  }; }
  else if (QuantEnum == 4  ) { sQuant = "Quant4"  ; cQuantVec = {      4}; }
  else { std::cout << "\n\nQuantEnum must be 0, 123 or 4. You gave me " << QuantEnum << "\n\n" << std::endl; return; }
  // Put it all together to make a string
  std::string PIDOutString = sIsFD+"-"+sIsFHC+"-"+sTrPID+"-"+sCVN+"-"+sQuant;

  // Where are the files that I'm loading in?
  std::string MCFile   = "/nova/ana/users/karlwarb/Analysis2020/PIDOpt/TrimCAF/200204/TrimCAF-2020-"+sIsFD+"-"+sIsFHC+"-"+sMC  +".root";
  std::string DataFile = "/nova/ana/users/karlwarb/Analysis2020/PIDOpt/TrimCAF/200204/TrimCAF-2020-"+sIsFD+"-"+sIsFHC+"-"+sData+".root";

  // --> Set my POT and Livetime accordingly to the above..
  if (IsFHC) { MyPOTScale = FHC_POT_Scale; MyLivScale = FHC_Liv_Scale; }
  else       { MyPOTScale = RHC_POT_Scale; MyLivScale = RHC_Liv_Scale; }

  // --- Print out some key info to check things are correct.
  std::cout << "\n Is this FD? " << IsFD << "\t Is this FHC? " << IsFHC << "\t QuantEnum? " << QuantEnum
            << "\n\t --> POT Scale is " << MyPOTScale << ", Livetime Scale is " << MyLivScale
            << "\n\t MCFile   is " << MCFile
            << "\n\t DataFile is " << DataFile
            << "\n\t Am I using ReMId as my Tracking PID?  " << IsReMId    << " --> sTrPID is " << sTrPID
            << "\n\t Am I using the Loose PTP trained CVN? " << IsLoosePTP << " --> sCVN is   " << sCVN
            << "\n Therefore PIDOutString is " << PIDOutString
            << "\n" << std::endl;

  // --- Define my OutDir and make sure it exists.
  std::string OutDir = "Plots-"+PIDOutString+"/";
  gSystem -> MakeDirectory( OutDir.c_str() );

  // ---- Define my TFiles. 
  TFile *MyMCSamp   = new TFile( MCFile  .c_str() );
  TFile *MyDataSamp = new TFile( DataFile.c_str() );

  // --- Load the MC TTreeReader, and set its branches
  TTreeReader MC_Reader("TrimTree", MyMCSamp);
  const int MC_TotalTree = MC_Reader.GetEntries(true);
  TTreeReaderValue<double> MC_Run          (MC_Reader, "Run"          );
  TTreeReaderValue<double> MC_Evt          (MC_Reader, "Evt"          );
  TTreeReaderValue<double> MC_PPFXWgt      (MC_Reader, "PPFXWgt"      );
  TTreeReaderValue<double> MC_XSecWgt      (MC_Reader, "XSecWgt"      );
  TTreeReaderValue<double> MC_TrueNuE      (MC_Reader, "TrueNuE"      );
  TTreeReaderValue<double> MC_IsNuMuCC     (MC_Reader, "IsNuMuCC"     );
  TTreeReaderValue<double> MC_IntMode      (MC_Reader, "IntMode"      );
  TTreeReaderValue<double> MC_DumbOsc      (MC_Reader, "DumbOsc"      );
  TTreeReaderValue<double> MC_BestFit2019  (MC_Reader, "BestFit2019"  );
  TTreeReaderValue<double> MC_RecoNuE      (MC_Reader, "RecoNuE"      );
  //TTreeReaderValue<double> MC_OldPreQuant  (MC_Reader, "OldPreQuant"  );
  //TTreeReaderValue<double> MC_CVNMuonOldPre(MC_Reader, "CVNMuonOldPre");
  //TTreeReaderValue<double> MC_CVNCosmOldPre(MC_Reader, "CVNCosmOldPre");
  TTreeReaderValue<double> MC_LooPTPQuant  (MC_Reader, "LooPTPQuant"  );
  TTreeReaderValue<double> MC_CosVeto      (MC_Reader, "CosVeto"      );
  TTreeReaderValue<double> MC_CVNMuonLooPTP(MC_Reader, "CVNMuonLooPTP");
  TTreeReaderValue<double> MC_CVNCosmLooPTP(MC_Reader, "CVNCosmLooPTP");  
  TTreeReaderValue<double> MC_ReMId        (MC_Reader, "ReMId"        );
  TTreeReaderValue<double> MC_CVNProngMuon (MC_Reader, "CVNProngMuon" );
  TTreeReaderValue<double> MC_CosRej2020   (MC_Reader, "CosRej2020"   );
  //TTreeReaderValue<double> MC_CosRej2018   (MC_Reader, "CosRej2018"   );

  // --- Load the Data TTreeReader, and set its branches
  TTreeReader Data_Reader("TrimTree", MyDataSamp);
  const int Data_TotalTree = Data_Reader.GetEntries(true);
  TTreeReaderValue<double> Data_Run          (Data_Reader, "Run"          );
  TTreeReaderValue<double> Data_Evt          (Data_Reader, "Evt"          );
  TTreeReaderValue<double> Data_PPFXWgt      (Data_Reader, "PPFXWgt"      );
  TTreeReaderValue<double> Data_XSecWgt      (Data_Reader, "XSecWgt"      );
  TTreeReaderValue<double> Data_TrueNuE      (Data_Reader, "TrueNuE"      );
  TTreeReaderValue<double> Data_IsNuMuCC     (Data_Reader, "IsNuMuCC"     );
  TTreeReaderValue<double> Data_IntMode      (Data_Reader, "IntMode"      );
  TTreeReaderValue<double> Data_DumbOsc      (Data_Reader, "DumbOsc"      );
  TTreeReaderValue<double> Data_BestFit2019  (Data_Reader, "BestFit2019"  );
  TTreeReaderValue<double> Data_RecoNuE      (Data_Reader, "RecoNuE"      );
  //TTreeReaderValue<double> Data_OldPreQuant  (Data_Reader, "OldPreQuant"  );
  //TTreeReaderValue<double> Data_CVNMuonOldPre(Data_Reader, "CVNMuonOldPre");
  //TTreeReaderValue<double> Data_CVNCosmOldPre(Data_Reader, "CVNCosmOldPre");
  TTreeReaderValue<double> Data_LooPTPQuant  (Data_Reader, "LooPTPQuant"  );
  TTreeReaderValue<double> Data_CosVeto      (Data_Reader, "CosVeto"      );
  TTreeReaderValue<double> Data_CVNMuonLooPTP(Data_Reader, "CVNMuonLooPTP");
  TTreeReaderValue<double> Data_CVNCosmLooPTP(Data_Reader, "CVNCosmLooPTP");  
  TTreeReaderValue<double> Data_ReMId        (Data_Reader, "ReMId"        );
  TTreeReaderValue<double> Data_CVNProngMuon (Data_Reader, "CVNProngMuon" );
  TTreeReaderValue<double> Data_CosRej2020   (Data_Reader, "CosRej2020"   );
  //TTreeReaderValue<double> Data_CosRej2018   (Data_Reader, "CosRej2018"   );

  // --- Make my 3D histograms.
  // First the full 0 to 5 GeV histograms.
  TH3D* hMC_Sig_FullE_Hist = new TH3D( "MCSig-FullE-Hist", "Monte Carlo Signal; TrPID Score; CVN Loose PTP Score; CosRej",
                                      nbinTrPID, MinTrPID, MaxTrPID, nbinCVN, MinCVN, MaxCVN, nbinCosRej, MinCosRej, MaxCosRej );
  TH3D* hMC_Bkg_FullE_Hist = new TH3D( "MCBkg-FullE-Hist", "Monte Carlo Background; TrPID Score; CVN Loose PTP Score; CosRej", 
                                      nbinTrPID, MinTrPID, MaxTrPID, nbinCVN, MinCVN, MaxCVN, nbinCosRej, MinCosRej, MaxCosRej );
  TH3D* hData_FullE_Hist   = new TH3D( "Data-FullE-Hist" , "Cosmic Events; TrPID Score; CVN Loose PTP Score; CosRej", 
                                      nbinTrPID, MinTrPID, MaxTrPID, nbinCVN, MinCVN, MaxCVN, nbinCosRej, MinCosRej, MaxCosRej );
  // Next the dip region hists - 1 to 2 GeV. 
  TH3D* hMC_Sig_Dip_Hist   = new TH3D( "MCSig-Dip-Hist", "Monte Carlo Signal; TrPID Score; CVN Loose PTP Score; CosRej",
                                      nbinTrPID, MinTrPID, MaxTrPID, nbinCVN, MinCVN, MaxCVN, nbinCosRej, MinCosRej, MaxCosRej );
  TH3D* hMC_Bkg_Dip_Hist   = new TH3D( "MCBkg-Dip-Hist", "Monte Carlo Background; TrPID Score; CVN Loose PTP Score; CosRej", 
                                      nbinTrPID, MinTrPID, MaxTrPID, nbinCVN, MinCVN, MaxCVN, nbinCosRej, MinCosRej, MaxCosRej );
  TH3D* hData_Dip_Hist     = new TH3D( "Data-Dip-Hist" , "Cosmic Events; TrPID Score; CVN Loose PTP Score; CosRej", 
                                      nbinTrPID, MinTrPID, MaxTrPID, nbinCVN, MinCVN, MaxCVN, nbinCosRej, MinCosRej, MaxCosRej );


  TH1D* hRecoNuE_Sig[9];
  TH1D* hRecoNuE_Bkg[9];
  TH1D* hRecoNuE_Cos[9];
  for (int hh=0; hh<9; ++hh) {
    std::string app = "";
    int Col = 0;
    if      (hh==0) { app = "FOMFull-AllQ"; Col = 1; }
    else if (hh==1) { app = "FOMFull-Q123"; Col = 1; }
    else if (hh==2) { app = "FOMFull-Q4"  ; Col = 1; }
    else if (hh==3) { app = "FOMDip-AllQ" ; Col = 2; }
    else if (hh==4) { app = "FOMDip-Q123" ; Col = 2; }
    else if (hh==5) { app = "FOMDip-Q4"   ; Col = 2; }
    else if (hh==6) { app = "BadCut-AllQ" ; Col = 4; }
    else if (hh==7) { app = "BadCut-Q123" ; Col = 4; }
    else if (hh==8) { app = "BadCut-Q4"   ; Col = 4; }
    std::string axis = ";Reconstructed Neutrino Enegry (GeV);Events / 0.1 GeV";
    hRecoNuE_Sig[hh] = new TH1D( TString("hSig-")+TString(app), axis.c_str(), HistBins, RecoEBins );
    hRecoNuE_Bkg[hh] = new TH1D( TString("hBkg-")+TString(app), axis.c_str(), HistBins, RecoEBins );
    hRecoNuE_Cos[hh] = new TH1D( TString("hCos-")+TString(app), axis.c_str(), HistBins, RecoEBins );
    hRecoNuE_Sig[hh] -> SetLineColor( 1 ); 
    hRecoNuE_Bkg[hh] -> SetLineColor( 2 ); 
    hRecoNuE_Cos[hh] -> SetLineColor( 4 );
  }

  // Declare two maps to store the FOMs.
  std::map< std::string, MyFOMParams > FOMMap;

  int NFills = 0;

  // Just access the data as if myPx and myPy were iterators (note the '*' in front of them):
  // ---- Monte Carlo Reader Loop
  while (MC_Reader.Next()) {
    // What entry am I looking at?
    int Entry = MC_Reader.GetCurrentEntry();
    // What is my Event Weight?
    //double EventWgt = *MC_PPFXWgt * *MC_XSecWgt * *MC_DumbOsc;
    double EventWgt = *MC_PPFXWgt * *MC_XSecWgt * *MC_BestFit2019;
    //double EventWgt = *MC_PPFXWgt * *MC_XSecWgt;
    // Do I want to write out some useful information?
    if ( bPrintInfo && Entry % 50000 == 0 ) {
      PrintInfo( Entry, MC_TotalTree, *MC_IsNuMuCC, *MC_TrueNuE, *MC_RecoNuE, EventWgt, *MC_CosVeto, 
                 *MC_CVNMuonLooPTP, *MC_ReMId, *MC_CVNProngMuon, *MC_CosRej2020 );
    }

    // --- Fill my Reco plots...
    if (*MC_IsNuMuCC) FillMyPlots( *MC_CVNMuonLooPTP, *MC_ReMId, *MC_CosRej2020, (int)*MC_LooPTPQuant, hRecoNuE_Sig, *MC_RecoNuE, EventWgt );
    else              FillMyPlots( *MC_CVNMuonLooPTP, *MC_ReMId, *MC_CosRej2020, (int)*MC_LooPTPQuant, hRecoNuE_Bkg, *MC_RecoNuE, EventWgt );

    // Set the Track PID score that I want to look at...
    double           MyTrPIDScore = *MC_ReMId;
    if (!IsReMId)    MyTrPIDScore = *MC_CVNProngMuon;
    // Set the CVN score that I want to look at...
    double           MyCVNScore = *MC_CVNMuonLooPTP;
    //if (!IsLoosePTP) MyCVNScore = *MC_CVNMuonOldPre;
    // ---- Does this event pass the cut points? If so, fill the Full E plots.
    if ( PassCutPoint( MyTrPIDScore, MinTrPID, MyCVNScore, MinCVN, *MC_LooPTPQuant, *MC_CosVeto, cQuantVec, *MC_RecoNuE, false ) ) {
      ++NFills; // increment the NFills Var.
      if (*MC_IsNuMuCC) hMC_Sig_FullE_Hist->Fill( MyTrPIDScore, MyCVNScore, *MC_CosRej2020, EventWgt );
      else              hMC_Bkg_FullE_Hist->Fill( MyTrPIDScore, MyCVNScore, *MC_CosRej2020, EventWgt );
    }
    // ---- Does this event pass the cut points? If so, fill the Dip E plots. 
    if ( PassCutPoint( MyTrPIDScore, MinTrPID, MyCVNScore, MinCVN, *MC_LooPTPQuant, *MC_CosVeto, cQuantVec, *MC_RecoNuE, true ) ) {
      if (*MC_IsNuMuCC) hMC_Sig_Dip_Hist->Fill( MyTrPIDScore, MyCVNScore, *MC_CosRej2020, EventWgt );
      else              hMC_Bkg_Dip_Hist->Fill( MyTrPIDScore, MyCVNScore, *MC_CosRej2020, EventWgt );
    } // Only fill if passes a loose presel.
  } // MC Reader

  // ---- Data Reader Loop
  while (Data_Reader.Next()) {
    // What entry am I looking at?
    int Entry = Data_Reader.GetCurrentEntry();
    // What is my Event Weight? --> Always 1 for Data.
    double EventWgt = 1;
    // Do I want to write out some useful information?
    if ( bPrintInfo && Entry % 10000 == 0 ) {
      PrintInfo( Entry, Data_TotalTree, *Data_IsNuMuCC, *Data_TrueNuE, *Data_RecoNuE, EventWgt, *Data_CosVeto,
                 *Data_CVNMuonLooPTP, *Data_ReMId, *Data_CVNProngMuon, *Data_CosRej2020 );
    }

    // --- Fill my Reco plots...
    FillMyPlots( *Data_CVNMuonLooPTP, *Data_ReMId, *Data_CosRej2020, (int)*Data_LooPTPQuant, hRecoNuE_Cos, *Data_RecoNuE, EventWgt );

    // Set the Track PID score that I want to look at...
    double           MyTrPIDScore = *Data_ReMId;
    if (!IsReMId)    MyTrPIDScore = *Data_CVNProngMuon;
    // Set the CVN score that I want to look at...
    double           MyCVNScore = *Data_CVNMuonLooPTP;
    //if (!IsLoosePTP) MyCVNScore = *Data_CVNMuonOldPre;
    // ---- Does this event pass the cut points? If so, fill the Full E plots.
    if ( PassCutPoint( MyTrPIDScore, MinTrPID, MyCVNScore, MinCVN, *Data_LooPTPQuant, *Data_CosVeto, cQuantVec, *Data_RecoNuE, false ) ) {
      hData_FullE_Hist->Fill( MyTrPIDScore, MyCVNScore, *Data_CosRej2020, EventWgt );
    }
    // ---- Does this event pass the cut points? If so, fill the Dip E plots. 
    if ( PassCutPoint( MyTrPIDScore, MinTrPID, MyCVNScore, MinCVN, *Data_LooPTPQuant, *Data_CosVeto, cQuantVec, *Data_RecoNuE, true ) ) {
      hData_Dip_Hist->Fill  ( MyTrPIDScore, MyCVNScore, *Data_CosRej2020, EventWgt );
    }
  } // Data Reader

  // ----- Now do my dirty plotting things....

  for (int hh=0; hh<9; ++hh) {
    hRecoNuE_Sig[hh] -> Scale( MyPOTScale );
    hRecoNuE_Bkg[hh] -> Scale( MyPOTScale );
    hRecoNuE_Cos[hh] -> Scale( MyLivScale );
  }
  
  std::cout << setprecision(3) << "\n After scaling, the different cut hists are as follows;"
            << "\n\n\t FOM Full --> ReMId 0.3, CVN 0.85, CosRej 0.45."
            << "\n\t\t AllQuant: Sig " << hRecoNuE_Sig[0]->Integral() << ", Bkg " << hRecoNuE_Bkg[0]->Integral() << ", Cos " << hRecoNuE_Cos[0]->Integral()
            << "\n\t\t Quant123: Sig " << hRecoNuE_Sig[1]->Integral() << ", Bkg " << hRecoNuE_Bkg[1]->Integral() << ", Cos " << hRecoNuE_Cos[1]->Integral()
            << "\n\t\t Quant4  : Sig " << hRecoNuE_Sig[2]->Integral() << ", Bkg " << hRecoNuE_Bkg[2]->Integral() << ", Cos " << hRecoNuE_Cos[2]->Integral()
            << "\n\t\t --> Q4Frac is " << hRecoNuE_Sig[2]->Integral() / hRecoNuE_Sig[0]->Integral()

            << "\n\n\t FOM Dip --> ReMId 0.32, CVN 0.91, CosRej 0.45."
            << "\n\t\t AllQuant: Sig " << hRecoNuE_Sig[3]->Integral() << ", Bkg " << hRecoNuE_Bkg[3]->Integral() << ", Cos " << hRecoNuE_Cos[3]->Integral()
            << "\n\t\t Quant123: Sig " << hRecoNuE_Sig[4]->Integral() << ", Bkg " << hRecoNuE_Bkg[4]->Integral() << ", Cos " << hRecoNuE_Cos[4]->Integral()
            << "\n\t\t Quant4  : Sig " << hRecoNuE_Sig[5]->Integral() << ", Bkg " << hRecoNuE_Bkg[5]->Integral() << ", Cos " << hRecoNuE_Cos[5]->Integral()
            << "\n\t\t --> Q4Frac is " << hRecoNuE_Sig[5]->Integral() / hRecoNuE_Sig[3]->Integral()

            << "\n\n\t Finally, current cut --> ReMId 0.7, CVN 0.82, CosRej 0.47."
            << "\n\t\t AllQuant: Sig " << hRecoNuE_Sig[6]->Integral() << ", Bkg " << hRecoNuE_Bkg[6]->Integral() << ", Cos " << hRecoNuE_Cos[6]->Integral()
            << "\n\t\t Quant123: Sig " << hRecoNuE_Sig[7]->Integral() << ", Bkg " << hRecoNuE_Bkg[7]->Integral() << ", Cos " << hRecoNuE_Cos[7]->Integral()
            << "\n\t\t Quant4  : Sig " << hRecoNuE_Sig[8]->Integral() << ", Bkg " << hRecoNuE_Bkg[8]->Integral() << ", Cos " << hRecoNuE_Cos[8]->Integral()
            << "\n\t\t --> Q4Frac is " << hRecoNuE_Sig[8]->Integral() / hRecoNuE_Sig[6]->Integral()
            << "\n\n" << std::endl;

  TCanvas *cFOMFull = new TCanvas(TString(PIDOutString)+TString("-FOMFull_QuickRecoE"), "Reco E for the cut optimised for FOM over All E");
  hRecoNuE_Sig[0]->Scale(0.1,"width"); hRecoNuE_Sig[0] -> Draw();
  hRecoNuE_Bkg[0]->Scale(0.1,"width"); hRecoNuE_Bkg[0] -> Draw("same");
  hRecoNuE_Cos[0]->Scale(0.1,"width"); hRecoNuE_Cos[0] -> Draw("same");
  cFOMFull    -> SaveAs( TString(OutDir)+TString(cFOMFull   ->GetName())+TString(".pdf") ); 

  TCanvas *cFOMDip = new TCanvas(TString(PIDOutString)+TString("-FOMDip_QuickRecoE"), "Reco E for the cut optimised for FOM in Dip");
  hRecoNuE_Sig[3]->Scale(0.1,"width"); hRecoNuE_Sig[3] -> Draw();
  hRecoNuE_Bkg[3]->Scale(0.1,"width"); hRecoNuE_Bkg[3] -> Draw("same");
  hRecoNuE_Cos[3]->Scale(0.1,"width"); hRecoNuE_Cos[3] -> Draw("same");
  cFOMDip     -> SaveAs( TString(OutDir)+TString(cFOMDip    ->GetName())+TString(".pdf") ); 

  TCanvas *cFOMCurrCut = new TCanvas(TString(PIDOutString)+TString("-FOMCurrCut_QuickRecoE"), "Reco E for the current cut");
  hRecoNuE_Sig[6]->Scale(0.1,"width"); hRecoNuE_Sig[6] -> Draw();
  hRecoNuE_Bkg[6]->Scale(0.1,"width"); hRecoNuE_Bkg[6] -> Draw("same");
  hRecoNuE_Cos[6]->Scale(0.1,"width"); hRecoNuE_Cos[6] -> Draw("same");
  cFOMCurrCut -> SaveAs( TString(OutDir)+TString(cFOMCurrCut->GetName())+TString(".pdf") ); 

  TCanvas *cFOMFull_Q123 = new TCanvas(TString(PIDOutString)+TString("-FOMFull_Q123"), "Reco E for the cut optimised for FOM over All E - Q123");
  hRecoNuE_Sig[1]->Scale(0.1,"width"); hRecoNuE_Sig[1] -> Draw();
  hRecoNuE_Bkg[1]->Scale(0.1,"width"); hRecoNuE_Bkg[1] -> Draw("same");
  hRecoNuE_Cos[1]->Scale(0.1,"width"); hRecoNuE_Cos[1] -> Draw("same");
  cFOMFull_Q123-> SaveAs( TString(OutDir)+TString(cFOMFull_Q123->GetName())+TString(".pdf") ); 

  TCanvas *cFOMFull_Q4 = new TCanvas(TString(PIDOutString)+TString("-FOMFull_Q4"), "Reco E for the cut optimised for FOM over All E - Q4");
  hRecoNuE_Sig[2]->Scale(0.1,"width"); hRecoNuE_Sig[2] -> Draw();
  hRecoNuE_Bkg[2]->Scale(0.1,"width"); hRecoNuE_Bkg[2] -> Draw("same");
  hRecoNuE_Cos[2]->Scale(0.1,"width"); hRecoNuE_Cos[2] -> Draw("same");
  cFOMFull_Q4-> SaveAs( TString(OutDir)+TString(cFOMFull_Q4->GetName())+TString(".pdf") ); 

  //return;
  
  // ---- AAAnnnd scale the plots...
  hMC_Sig_FullE_Hist->Scale(MyPOTScale); hMC_Sig_Dip_Hist->Scale(MyPOTScale);
  hMC_Bkg_FullE_Hist->Scale(MyPOTScale); hMC_Bkg_Dip_Hist->Scale(MyPOTScale);
  hData_FullE_Hist  ->Scale(MyLivScale); hData_Dip_Hist  ->Scale(MyLivScale);

  // --- Declare some 2D FOM vs background plots.
  double FOMMin_f, FOMMax_f, FOMMin_d, FOMMax_d;
  int    nbinsFOM_f, nbinsFOM_d;
  if      ( IsFHC && QuantEnum == 0   ) { FOMMin_f = 6.; FOMMax_f = 12.; nbinsFOM_f = 60; FOMMin_d = 3.; FOMMax_d = 5.5; nbinsFOM_d = 70; }
  else if ( IsFHC && QuantEnum == 123 ) { FOMMin_f = 6.; FOMMax_f = 10.; nbinsFOM_f = 40; FOMMin_d = 3.; FOMMax_d = 4.2; nbinsFOM_d = 24; }
  else if ( IsFHC && QuantEnum == 4   ) { FOMMin_f = 3.; FOMMax_f = 6.5; nbinsFOM_f = 35; FOMMin_d = 1.; FOMMax_d = 3.5; nbinsFOM_d = 50; }
  else if (!IsFHC && QuantEnum == 0   ) { FOMMin_f = 6.; FOMMax_f = 10.; nbinsFOM_f = 40; FOMMin_d = 2.; FOMMax_d = 4.5; nbinsFOM_d = 50; }
  else if (!IsFHC && QuantEnum == 123 ) { FOMMin_f = 6.; FOMMax_f = 8.5; nbinsFOM_f = 25; FOMMin_d = 2.; FOMMax_d = 4  ; nbinsFOM_d = 40; }
  else if (!IsFHC && QuantEnum == 4   ) { FOMMin_f = 3.; FOMMax_f = 5.5; nbinsFOM_f = 25; FOMMin_d = 0.; FOMMax_d = 2.5; nbinsFOM_d = 50; }
  TH2 *hFOMFull_Bkg = new TH2D("FOMFull_Bkg_Hist", "Full FOM vs Number of Bkg Events; Total Background; FOM", 150, 0, 30, nbinsFOM_f, FOMMin_f, FOMMax_f);
  TH2 *hFOMDip_Bkg  = new TH2D("FOMDip_Bkg_Hist" , "Dip FOM vs Number of Bkg Events; Total Background; FOM" , 150, 0, 15, nbinsFOM_d, FOMMin_d, FOMMax_d);

  // --- Declare some VarX vs Var Y plots to give a feel for the density of good cut points.
  // First the full E Dist.
  TH2 *hGoodFOM_TrPID_CVN    = new TH2D("GoodFOM_TrPID_CVN"   , "Density of TrackPID vs CVN Cuts above FOM Thresholds; Track PID Cut; CVN Cut", 
                                        nbinTrPID, MinTrPID, MaxTrPID, nbinCVN   , MinCVN   , MaxCVN );
  TH2 *hGoodFOM_TrPID_CosRej = new TH2D("GoodFOM_TrPID_CosRej", "Density of TrackPID vs CosRej Cuts above FOM Thresholds; Track PID Cut; CosRej Cut", 
                                        nbinTrPID, MinTrPID, MaxTrPID, nbinCosRej, MinCosRej, MaxCosRej );
  TH2 *hGoodFOM_CVN_CosRej   = new TH2D("GoodFOM_CVN_CosRej"  , "Density of CosRej vs CVN Cuts above FOM Thresholds; CVN Cut; CosRej Cut", 
                                        nbinCVN  , MinCVN  , MaxCVN  , nbinCosRej, MinCosRej, MaxCosRej );

  // ---- Figure out my FOMs.
  double MaxFOM_Full = 0, MaxFOM_Dip = 0;
  // ---- What are my FOM thresholds for an "interesting event"
  double FOMInt_Full = 0, FOMInt_Dip = 0;
  // FHC
  if      ( IsFHC && QuantEnum == 0   ) { FOMInt_Full = 10.8; FOMInt_Dip = 4.8; }
  else if ( IsFHC && QuantEnum == 123 ) { FOMInt_Full = 9.1 ; FOMInt_Dip = 3.9; }
  else if ( IsFHC && QuantEnum == 4   ) { FOMInt_Full = 5.9 ; FOMInt_Dip = 2.8; }
  // RHC
  else if (!IsFHC && QuantEnum == 0   ) { FOMInt_Full = 9.5 ; FOMInt_Dip = 3.9; }
  else if (!IsFHC && QuantEnum == 123 ) { FOMInt_Full = 7.9 ; FOMInt_Dip = 3.4; }
  else if (!IsFHC && QuantEnum == 4   ) { FOMInt_Full = 5.2 ; FOMInt_Dip = 2.0; }
  // And loop through my combinations.
  for (    double trpid = MinTrPID ; trpid < MaxTrPID ; trpid=trpid+0.01) {
    std::cout << "\tNow looking at combinations with trpid " << trpid << ". FOMMap currently has size " << FOMMap.size() << std::endl;
    for (  double cvn   = MinCVN   ; cvn   < MaxCVN   ; cvn  =cvn  +0.01) {
      for (double cos   = MinCosRej; cos   < MaxCosRej; cos  =cos  +0.01) {
        // Calculate the bins to integrate over.
        int bin_trpid = hData_FullE_Hist->GetXaxis()->FindBin( trpid );
        int bin_cvn   = hData_FullE_Hist->GetYaxis()->FindBin( cvn   );
        int bin_cos   = hData_FullE_Hist->GetZaxis()->FindBin( cos   );
        // Set the variables to send to CalcFOM.
        MyFOMParams TempFOM;
        // Calculate my FOMs. 
        CalcFOM( hMC_Sig_FullE_Hist, hMC_Bkg_FullE_Hist, hData_FullE_Hist, bin_trpid, bin_cvn, bin_cos, false, TempFOM.FOM_Full, TempFOM.TotalSig_Full, TempFOM.TotalBkg_Full );
        
        CalcFOM( hMC_Sig_Dip_Hist  , hMC_Bkg_Dip_Hist  , hData_Dip_Hist  , bin_trpid, bin_cvn, bin_cos, true , TempFOM.FOM_Dip , TempFOM.TotalSig_Dip , TempFOM.TotalBkg_Dip  );
        // Is this the largest FOM?
        if (TempFOM.FOM_Full > MaxFOM_Full) MaxFOM_Full = TempFOM.FOM_Full;
        if (TempFOM.FOM_Dip  > MaxFOM_Dip ) MaxFOM_Dip  = TempFOM.FOM_Dip;
        // Fill my 2D plots.
        hFOMFull_Bkg -> Fill( TempFOM.TotalBkg_Full, TempFOM.FOM_Full );
        hFOMDip_Bkg  -> Fill( TempFOM.TotalBkg_Dip , TempFOM.FOM_Dip  );
        
        // --- If the FOM is good, and only if it is good, add it to my Map for closer inspection later...
        if ( TempFOM.FOM_Full > FOMInt_Full || TempFOM.FOM_Dip > FOMInt_Dip ) {
          // What is my cut string for this combination?
          std::ostringstream CutString;
          CutString << std::setprecision(2) << std::fixed << "TrPID-" << trpid << "-CVN-" << cvn << "-CosRej-" << cos;
          // Write out the cut poisiton info.
          if (bPrintInfo 
              || ( bin_trpid==1  && bin_cvn==46 && bin_cos==6 ) // 0.3, 0.85, 0.45
              || ( bin_trpid==1  && bin_cvn==52 && bin_cos==6 ) // 0.3, 0.91, 0.45
              || ( bin_trpid==41 && bin_cvn==43 && bin_cos==8 ) // 0.7, 0.82, 0.47
              ) {
            std::cout << setprecision(4)
                      << "Looking at My TrPID " << trpid     << ", My CVN     " << cvn     << ", CosRej2020 " << cos << " --> " << CutString.str()
                      << "\n\t bin_rem is     " << bin_trpid << ", bin_cvn is " << bin_cvn << ", bin_cos is " << bin_cos
                      << "\n\t TotSig_Full is " << TempFOM.TotalSig_Full << ", TotBkg_Full is " << TempFOM.TotalBkg_Full << " ---> FOM_Full is " << TempFOM.FOM_Full
                      << "\n\t TotSig_Dip is  " << TempFOM.TotalSig_Dip  << ", TotBkg_Dip is  " << TempFOM.TotalBkg_Dip  << " ---> FOM_Dip is  " << TempFOM.FOM_Dip
                      << std::endl;
          }
          // Fill my maps.
          FOMMap[ CutString.str() ] = TempFOM;
          // Fill the 2D VarX vs VarY plots.
          hGoodFOM_TrPID_CVN    -> Fill( trpid, cvn );
          hGoodFOM_TrPID_CosRej -> Fill( trpid, cos );
          hGoodFOM_CVN_CosRej   -> Fill( cvn  , cos );
        }

      } // loop over cosrej
    }   // loop over cvn
  }     // loop over remid 

  // --- Canvas for the FOM vs Bkg plots.
  TCanvas* cFOMFull_Bkg = new TCanvas(TString(PIDOutString)+TString("-FOMFull_Bkg"), "FOM over All E vs Total Bkg");
  hFOMFull_Bkg->Draw("colz");
  TCanvas* cFOMDip_Bkg  = new TCanvas(TString(PIDOutString)+TString("-FOMDip_Bkg") , "FOM in Dip Region vs Total Bkg");
  hFOMDip_Bkg ->Draw("colz");

  // --- Canvas for the VarX vs VarY Good FOM density plots.
  TCanvas *cGoodFOM_TrPID_CVN    = new TCanvas(TString(PIDOutString)+TString("-GoodFOM_TrPID_CVN")   , "Density of TrackPID vs CVN Cuts above FOM Thresholds"    );
  hGoodFOM_TrPID_CVN   -> Draw("colz");
  TCanvas *cGoodFOM_TrPID_CosRej = new TCanvas(TString(PIDOutString)+TString("-GoodFOM_TrPID_CosRej"), "Density of TrackPID vs CosRej Cuts above FOM Thresholds" );
  hGoodFOM_TrPID_CosRej -> Draw("colz");
  TCanvas *cGoodFOM_CVN_CosRej   = new TCanvas(TString(PIDOutString)+TString("-GoodFOM_CVN_CosRej")  , "Density of CVN vs CosRej Cuts above FOM Thresholds"      );
  hGoodFOM_CVN_CosRej   -> Draw("colz");  

  // --- Save the canvases.
  cFOMFull_Bkg         -> SaveAs( TString(OutDir)+TString(cFOMFull_Bkg         ->GetName())+TString(".pdf") );
  cFOMDip_Bkg          -> SaveAs( TString(OutDir)+TString(cFOMDip_Bkg          ->GetName())+TString(".pdf") );
  cGoodFOM_TrPID_CVN   -> SaveAs( TString(OutDir)+TString(cGoodFOM_TrPID_CVN   ->GetName())+TString(".pdf") );
  cGoodFOM_TrPID_CosRej-> SaveAs( TString(OutDir)+TString(cGoodFOM_TrPID_CosRej->GetName())+TString(".pdf") );
  cGoodFOM_CVN_CosRej  -> SaveAs( TString(OutDir)+TString(cGoodFOM_CVN_CosRej  ->GetName())+TString(".pdf") );


  // ---- After all of that how many combinations were considered "good"?
  std::cout << "\n After all of that. The MaxFOM_Full is " << MaxFOM_Full << ", the MaxFOM_Dip is " << MaxFOM_Dip
            << "\n I considered " << FOMMap.size() << " of " << nbinTrPID*nbinCVN*nbinCosRej << " cut combinations as good."
            << " Lets look at what cuts are best for different criteria."
            << std::endl;
  const int PrintLim = 10; int Printed;

  // ------------------------------------------------
  // --- Invert the map to be <FOM_Full, CutString>.
  std::cout << "\nFirst, Lets sort the Map by FOM_Full..." << std::endl;
  Printed = 0;
  std::multimap< double, std::string, std::greater<double> > FOMMap_SortFullFOM;
  for (std::pair< std::string, MyFOMParams > FOMMapIt : FOMMap) {
    FOMMap_SortFullFOM.insert( std::make_pair( FOMMapIt.second.FOM_Full, FOMMapIt.first ) );
  }
  for (std::pair< double, std::string > SortFOMIt : FOMMap_SortFullFOM) {
    if (Printed > PrintLim) break;
    PrintMyMap( FOMMap, SortFOMIt.second );
    ++Printed;
  }

  // ------------------------------------------------
  // --- Invert the map to be <FOM_Dip, CutString>.
  std::cout << "\nNext, Lets sort the Map by FOM_Dip..." << std::endl;
  Printed = 0;
  std::multimap< double, std::string, std::greater<double> > FOMMap_SortDipFOM;
  for (std::pair< std::string, MyFOMParams > FOMMapIt : FOMMap) {
    FOMMap_SortDipFOM.insert( std::make_pair( FOMMapIt.second.FOM_Dip, FOMMapIt.first ) );
  }
  for (std::pair< double, std::string > SortFOMIt : FOMMap_SortDipFOM) {
    if (Printed > PrintLim) break;
    PrintMyMap( FOMMap, SortFOMIt.second );
    ++Printed;
  }

  // ------------------------------------------------
  // --- Invert the map to be <TotalSig_Dip, CutString>.
  std::cout << "\nNext, Lets sort the Map by TotalSig_Dip..." << std::endl;
  Printed = 0;
  std::multimap< double, std::string, std::greater<double> > FOMMap_SortSig_Dip;
  for (std::pair< std::string, MyFOMParams > FOMMapIt : FOMMap) {
    FOMMap_SortSig_Dip.insert( std::make_pair( FOMMapIt.second.TotalSig_Dip, FOMMapIt.first ) );
  }
  for (std::pair< double, std::string > SortFOMIt : FOMMap_SortSig_Dip) {
    if (Printed > PrintLim) break;
    PrintMyMap( FOMMap, SortFOMIt.second );
    ++Printed;
  }

  // ------------------------------------------------
  // --- Invert the map to be <TotalBkg_Dip, CutString>.
  std::cout << "\nNext, Lets sort the Map by TotalBkg_Dip..." << std::endl;
  Printed = 0;
  std::multimap< double, std::string > FOMMap_SortBkg_Dip;
  for (std::pair< std::string, MyFOMParams > FOMMapIt : FOMMap) {
    FOMMap_SortBkg_Dip.insert( std::make_pair( FOMMapIt.second.TotalBkg_Dip, FOMMapIt.first ) );
  }
  for (std::pair< double, std::string > SortFOMIt : FOMMap_SortBkg_Dip) {
    if (Printed > PrintLim) break;
    PrintMyMap( FOMMap, SortFOMIt.second );
    ++Printed;
  }

  return;
} // Main Class.

// ---------------------------------------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------------------------------------
void PrintInfo( int Entry, int TotEnt, double IsNuMuCC, double TrueNuE, double RecoNuE, double EventWgt, double CosVeto, 
                double CVNMuonLooPTP, double ReMId, double CVNProngMuon, double CosRej2020 ) {
  
  std::cout << setprecision(4) 
            << "\n Looking at entry " << Entry << " of " << TotEnt << " ---> " << 100*(double)Entry/(double)TotEnt << "% complete... \n\t"
            << "  IsNuMuCC " << IsNuMuCC << ", TrueNuE " << TrueNuE << ", RecoNuE " << RecoNuE << ", EventWgt " << EventWgt
            << ", Pass CosVeto? " << CosVeto
            << ", CVNMuonLooPTP " << CVNMuonLooPTP
    //<< ", CVNMuonOldPre " << CVNMuonOldPre
            << ", ReMId "         << ReMId
            << ", CVNProngMuon "  << CVNProngMuon
            << ", CosRej2020 "    << CosRej2020
            << std::endl;
} // PrintInfo
// ---------------------------------------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------------------------------------
void CalcFOM( TH3D* Sig, TH3D* Bkg, TH3D* Data, int bintrpid, int bincvn, int bincos, bool IsFD, double &FOM, 
              double &TotSig, double &TotBkg ) {
  // Calc the integrals for each histogram.
  double SigInt  = Sig  -> Integral( bintrpid, nbinTrPID+1, bincvn, nbinCVN+1, bincos, nbinCosRej+1 );
  double BkgInt  = Bkg  -> Integral( bintrpid, nbinTrPID+1, bincvn, nbinCVN+1, bincos, nbinCosRej+1 );
  double DataInt = Data -> Integral( bintrpid, nbinTrPID+1, bincvn, nbinCVN+1, bincos, nbinCosRej+1 );
  
  // Set the inputs.
  TotSig = SigInt;
  TotBkg = BkgInt + DataInt;

  // Calc the FOM
  FOM = TotSig / std::pow( TotSig + TotBkg, 0.5 );

} // CalcFOM
// ---------------------------------------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------------------------------------
void PrintMyMap( std::map< std::string, MyFOMParams > MainMap, std::string Key ) {
  std::cout << "\t"<<Key << setprecision(4)
            << ", FOM_Full "       << MainMap[ Key ].FOM_Full
            << ", TotalSig_Full: " << MainMap[ Key ].TotalSig_Full
            << ", TotalBkg_Full: " << MainMap[ Key ].TotalBkg_Full
            << ", FOM_Dip "        << MainMap[ Key ].FOM_Dip
            << ", TotalSig_Dip: "  << MainMap[ Key ].TotalSig_Dip
            << ", TotalBkg_Dip: "  << MainMap[ Key ].TotalBkg_Dip
            << std::endl;
} // PrintMyMap
// ---------------------------------------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------------------------------------
bool PassCutPoint( double PIDScore, double PIDCut, double CVNScore, double CVNCut, 
                   int Quant, int CosVeto, std::vector<int> cQuantVec, double RecoE, bool isDip ) {

  // Easiest check is on the TrPID and CVN score.
  if ( PIDScore >= PIDCut && CVNScore >= CVNCut
       && CosVeto == 1
       ) {
    // Loop through cQuantVec
    for (size_t cq=0; cq < cQuantVec.size(); ++cq) {
      // Check that this quant is the one I'm interested in.
      if (Quant == cQuantVec[cq]) {
        // Finally, if in the dip check the energy range, if not, return true.
        if ( (RecoE > 1 && RecoE < 2) || !isDip ) return true;
      }
    }
  }
  return false;
} // PassCutPoint
// ---------------------------------------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------------------------------------
void FillMyPlots ( double CVN, double ReMId, double CosRej, int Quant, TH1D* Hists[9], double RecoE, double Wgt ) {
  // -- First check FOM Full
  if ( ReMId > 0.3 && CVN > 0.85 && CosRej > 0.45 ) {
    Hists[0] -> Fill( RecoE, Wgt );
    if (Quant != 4) Hists[1] -> Fill( RecoE, Wgt ); 
    else            Hists[2] -> Fill( RecoE, Wgt ); 
  }
  // -- Then check FOM Dip
  if ( ReMId > 0.32 && CVN > 0.91 && CosRej > 0.45 ) {
    Hists[3] -> Fill( RecoE, Wgt );
    if (Quant != 4) Hists[4] -> Fill( RecoE, Wgt ); 
    else            Hists[5] -> Fill( RecoE, Wgt ); 
  }
  // -- Finally, check the currently implemented cut
  if ( ReMId > 0.7 && CVN > 0.82 && CosRej > 0.47 ) {
    Hists[6] -> Fill( RecoE, Wgt );
    if (Quant != 4) Hists[7] -> Fill( RecoE, Wgt ); 
    else            Hists[8] -> Fill( RecoE, Wgt ); 
  }
  return;
}
// ---------------------------------------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------------------------------------