FD_Data_PosComp.C

Go to the documentation of this file.

//..................................
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Analysis/Prod4Loaders.h"

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "CAFAna/Cuts/SpillCuts.h"

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

#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "Utilities/rootlogon.C"

#include "TAttLine.h"
#include "TAttMarker.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TStyle.h"

#include <cmath>
#include <iostream>
#include <fstream> 
#include <vector>
#include <string>
#include <TROOT.h>
#include <TStyle.h>
//..................................
using namespace ana;

std::string PlotDir = "plots_position/";

//..................................
void MakeTextFile( std::string Plane, std::string Overlay ) {
  // Determine output name, and caption.
  std::string FNa = PlotDir+Plane+"_"+Overlay+".txt";
  std::string Cap = "Plot showing the distribution of NuMi data events passing the NuMu FD selection in the "+Plane+" plane, with "+Overlay;
  // Write to file.
  std::ofstream TxtOut ( FNa.c_str(), std::ofstream::out );
  TxtOut << Cap;
  TxtOut.close();
  // Done.
  return;
}

//..................................
void FD_Data_PosComp( bool IsFHC = true )
{
  // ---- First off, lets set the styles...
  gROOT->SetStyle("novaStyle");
  //gStyle->SetPalette(87); // https://root.cern.ch/doc/master/classTColor.html

  // --- Define my definitions
  std::string AllData = "prod4_sumrestricteddecaf_fd_numi_fhc_full_goodruns_numu2018";
  std::string AllMC   = "prod_sumdecaf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1_numu2018_str5";
  std::string sFHCRHC = "FHC_";
  double      MyPOT   = kAna2018FHCPOT;
  if (!IsFHC) {
    AllData = "prod4_sumrestricteddecaf_fd_numi_rhc_full_goodruns_numu2018";
    AllMC   = "prod_sumdecaf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1_numu2018_str5";
    sFHCRHC = "RHC_";
    MyPOT   = kAna2018RHCPOT;
  }
  
  PlotDir += sFHCRHC;
  
  // --- Define my binnnings.
  const Binning XBins = Binning::Simple(20, -10, 10 );
  const Binning YBins = Binning::Simple(20, -10, 10 );
  const Binning ZBins = Binning::Simple(31, -2 , 60);

  // --- Define my loaders
  SpectrumLoader allswap_MC( AllMC );

  // --- Define my Spill cuts
  allswap_MC.SetSpillCut(kStandardSpillCuts);

  // --- Define my spectra.
  // X vs Y
  Spectrum* XvsY_MC = new Spectrum("Positions in X vs Y (m)", allswap_MC, XBins, kTrkStartX, YBins, kTrkStartY, kNumuCutFD2018, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
  // Y vs Z
  Spectrum* ZvsY_MC = new Spectrum("Positions in Z vs Y (m)", allswap_MC, ZBins, kTrkStartZ, YBins, kTrkStartY, kNumuCutFD2018, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
  // X Z vs X
  Spectrum* ZvsX_MC = new Spectrum("Positions in Z vs X (m)", allswap_MC, ZBins, kTrkStartZ, XBins, kTrkStartX, kNumuCutFD2018, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);

  allswap_MC.Go();

  // .............................. Y vs X ..............................
  TH2* hXvsY_MC = XvsY_MC -> ToTH2(MyPOT);
  hXvsY_MC->GetXaxis()->SetTitle("Vertex X position (m)");
  hXvsY_MC->GetYaxis()->SetTitle("Vertex Y position (m)");
  hXvsY_MC->GetYaxis()->SetTitleOffset(0.7);
  hXvsY_MC->GetZaxis()->SetTitle("Arbitrary units");
  CenterTitles(hXvsY_MC);

  TGraph* gXvsY = MakeGraph(AllData, kNumuCutFD2018 && kInBeamSpill, kTrkStartX, kTrkStartY, &kStandardSpillCuts);
  gXvsY->SetMarkerStyle(20);
  gXvsY->SetMarkerSize(2.5);
  gXvsY->GetXaxis()->SetTitle("Vertex X position (m)");
  gXvsY->GetXaxis()->CenterTitle();
  gXvsY->GetYaxis()->SetTitle("Vertex Y position (m)");
  gXvsY->GetYaxis()->SetTitleOffset(0.7);
  gXvsY->GetYaxis()->CenterTitle();

  // First draw canvas with data only.
  TCanvas* cXvsY    = new TCanvas("cXvsY","X pos vs Y pos no MC", 900, 800);
  gXvsY -> Draw("AP");
  Preliminary();
  // Define my lines.
  TLine *XvsY_L1 = new TLine( -8, -8, -8,  7 ); XvsY_L1->SetLineColor(kRed); XvsY_L1->SetLineWidth(4); XvsY_L1->SetLineStyle(2); XvsY_L1->Draw("same");
  TLine *XvsY_L2 = new TLine(  8, -8,  8,  7 ); XvsY_L2->SetLineColor(kRed); XvsY_L2->SetLineWidth(4); XvsY_L2->SetLineStyle(2); XvsY_L2->Draw("same");
  TLine *XvsY_L3 = new TLine( -8, -8,  8, -8 ); XvsY_L3->SetLineColor(kRed); XvsY_L3->SetLineWidth(4); XvsY_L3->SetLineStyle(2); XvsY_L3->Draw("same");
  TLine *XvsY_L4 = new TLine( -8,  7,  8,  7 ); XvsY_L4->SetLineColor(kRed); XvsY_L4->SetLineWidth(4); XvsY_L4->SetLineStyle(2); XvsY_L4->Draw("same");
  gXvsY -> GetXaxis()->SetRangeUser(-10, 10); 
  gXvsY -> GetYaxis()->SetRangeUser(-10, 10);
  cXvsY -> SaveAs( TString(PlotDir)+TString("XvsY_noMCOverlay.pdf") );
  MakeTextFile( "XvsY", "noMCOverlay" );

  // Scale histogram
  double max1 = hXvsY_MC->GetBinContent(hXvsY_MC->GetMaximumBin());
  hXvsY_MC->Scale(1./max1);
  // Set range
  hXvsY_MC->GetXaxis()->SetRangeUser(-8, 8);
  hXvsY_MC->GetYaxis()->SetRangeUser(-8, 7);
  // Now draw the canvas with MC background, and zoom in on detector.
  TCanvas* cXvsY_MC = new TCanvas("cXvsY_MC","X pos vs Y pos MC Overlay", 900, 800);
  cXvsY_MC -> SetRightMargin (.15);
  hXvsY_MC -> Draw("COLZ");
  gXvsY    -> Draw("P");
  Preliminary();
  cXvsY_MC -> SaveAs( TString(PlotDir)+TString("XvsY_MCOverlay.pdf") );
  MakeTextFile( "XvsY", "MCOverlay" );

  // .............................. Z vs Y ..............................
  TH2* hZvsY_MC = ZvsY_MC->ToTH2(MyPOT);
  hZvsY_MC->GetYaxis()->SetTitle("Vertex Y position (m)");
  hZvsY_MC->GetYaxis()->SetTitleOffset(0.5);
  hZvsY_MC->GetXaxis()->SetTitle("Vertex Z position (m)");
  hZvsY_MC->GetZaxis()->SetTitle("Arbitrary units");
  CenterTitles(hZvsY_MC);

  TGraph* gZvsY = MakeGraph(AllData, kNumuCutFD2018 && kInBeamSpill, kTrkStartZ, kTrkStartY, &kStandardSpillCuts);
  gZvsY->SetMarkerStyle(20);
  gZvsY->SetMarkerSize(2.5);
  gZvsY->GetXaxis()->SetTitle("Vertex Z position (m)");
  gZvsY->GetXaxis()->CenterTitle();
  gZvsY->GetYaxis()->SetTitle("Vertex Y position (m)");
  gZvsY->GetYaxis()->SetTitleOffset(0.5);
  gZvsY->GetYaxis()->CenterTitle();

  // First draw canvas with data only.
  TCanvas* cZvsY    = new TCanvas("cZvsY","Z pos vs Y pos no MC", 2500, 800);
  gZvsY -> Draw("AP");
  Preliminary();
  // Define my lines.
  TLine *ZvsY_L1 = new TLine(  0, -8, 0 ,  7 ); ZvsY_L1->SetLineColor(kRed); ZvsY_L1->SetLineWidth(4); ZvsY_L1->SetLineStyle(2); ZvsY_L1->Draw("same");
  TLine *ZvsY_L2 = new TLine( 58, -8, 58,  7 ); ZvsY_L2->SetLineColor(kRed); ZvsY_L2->SetLineWidth(4); ZvsY_L2->SetLineStyle(2); ZvsY_L2->Draw("same");
  TLine *ZvsY_L3 = new TLine(  0, -8, 58, -8 ); ZvsY_L3->SetLineColor(kRed); ZvsY_L3->SetLineWidth(4); ZvsY_L3->SetLineStyle(2); ZvsY_L3->Draw("same");
  TLine *ZvsY_L4 = new TLine(  0,  7, 58,  7 ); ZvsY_L4->SetLineColor(kRed); ZvsY_L4->SetLineWidth(4); ZvsY_L4->SetLineStyle(2); ZvsY_L4->Draw("same");
  gZvsY -> GetXaxis()->SetRangeUser( -5 , 60 ); // Z Pos
  gZvsY -> GetYaxis()->SetRangeUser( -10,  10 ); // Y Pos
  cZvsY -> SaveAs( TString(PlotDir)+TString("ZvsY_noMCOverlay.pdf") );
  MakeTextFile( "ZvsY", "noMCOverlay" );

  // Scale the histogram
  double max2 = hZvsY_MC->GetBinContent(hZvsY_MC->GetMaximumBin());
  hZvsY_MC->Scale(1./max2);
  // Set the range
  hZvsY_MC->GetXaxis()->SetRangeUser( 0 , 58 ); // Z Pos
  hZvsY_MC->GetYaxis()->SetRangeUser( -8, 7  ); // Y Pos
  // Now draw the canvas with MC background, and zoom in on detector.
  TCanvas* cZvsY_MC = new TCanvas("cZvsY_MC","Z pos vs Y pos MC Overlay", 2500, 800);
  cZvsY_MC -> SetRightMargin (.15);
  hZvsY_MC -> Draw("COLZ");
  gZvsY    -> Draw("P");
  Preliminary();
  cZvsY_MC -> SaveAs( TString(PlotDir)+TString("ZvsY_MCOverlay.pdf") );
  MakeTextFile( "ZvsY", "MCOverlay" );

  // .............................. Z vs X ..............................
  TH2* hZvsX_MC = ZvsX_MC->ToTH2(MyPOT);
  hZvsX_MC->GetXaxis()->SetTitle("Vertex Z position (m)");
  hZvsX_MC->GetYaxis()->SetTitle("Vertex X position (m)");
  hZvsX_MC->GetYaxis()->SetTitleOffset(0.5);
  hZvsX_MC->GetZaxis()->SetTitle("Arbitrary units");
  CenterTitles(hZvsX_MC);

  TGraph* gZvsX = MakeGraph(AllData, kNumuCutFD2018 && kInBeamSpill, kTrkStartZ, kTrkStartX, &kStandardSpillCuts);
  gZvsX->SetMarkerStyle(20);
  gZvsX->SetMarkerSize(2.5);
  gZvsX->GetXaxis()->SetTitle("Vertex Z position (m)");
  gZvsX->GetXaxis()->CenterTitle();
  gZvsX->GetYaxis()->SetTitle("Vertex X position (m)");
  gZvsX->GetYaxis()->SetTitleOffset(0.5);

  // First draw canvas with data only.
  TCanvas* cZvsX    = new TCanvas("cZvsX","Z pos vs X pos no MC", 2500, 800);
  gZvsX -> Draw("AP");
  Preliminary();
  // Define my lines.
  TLine *ZvsX_L1 = new TLine(  0, -8, 0 ,  8 ); ZvsX_L1->SetLineColor(kRed); ZvsX_L1->SetLineWidth(4); ZvsX_L1->SetLineStyle(2); ZvsX_L1->Draw("same");
  TLine *ZvsX_L2 = new TLine( 58, -8, 58,  8 ); ZvsX_L2->SetLineColor(kRed); ZvsX_L2->SetLineWidth(4); ZvsX_L2->SetLineStyle(2); ZvsX_L2->Draw("same");
  TLine *ZvsX_L3 = new TLine(  0, -8, 58, -8 ); ZvsX_L3->SetLineColor(kRed); ZvsX_L3->SetLineWidth(4); ZvsX_L3->SetLineStyle(2); ZvsX_L3->Draw("same");
  TLine *ZvsX_L4 = new TLine(  0,  8, 58,  8 ); ZvsX_L4->SetLineColor(kRed); ZvsX_L4->SetLineWidth(4); ZvsX_L4->SetLineStyle(2); ZvsX_L4->Draw("same");
  gZvsX -> GetXaxis()->SetRangeUser( -5 , 60 ); // Z Pos
  gZvsX -> GetYaxis()->SetRangeUser( -10, 10 ); // Y Pos
  cZvsX    -> SaveAs( TString(PlotDir)+TString("ZvsX_noMCOverlay.pdf") );
  MakeTextFile( "ZvsX", "noMCOverlay" );

  // Scale the histogram
  double max3 = hZvsX_MC->GetBinContent(hZvsX_MC->GetMaximumBin());
  hZvsX_MC->Scale(1./max3);
  // Set the range
  hZvsX_MC->GetXaxis()->SetRangeUser( 0 , 58 ); // Z Pos
  hZvsX_MC->GetYaxis()->SetRangeUser( -8, 7  ); // Y Pos
  // Now draw the canvas with MC background, and zoom in on detector.
  TCanvas* cZvsX_MC = new TCanvas("cZvsX_MC","Z pos vs X pos MC Overlay", 2500, 800);
  cZvsX_MC -> SetRightMargin (.15);
  hZvsX_MC -> Draw("COLZ");
  gZvsX    -> Draw("P");
  Preliminary();
  cZvsX_MC -> SaveAs( TString(PlotDir)+TString("ZvsX_MCOverlay.pdf") );
  MakeTextFile( "ZvsX", "MCOverlay" );

  // .............................. Define an outfile ..............................
  std::string OutName = PlotDir+"PosComparisons.root";
  TFile *OutFile = new TFile(OutName.c_str(), "RECREATE");
  OutFile->cd();
  // Y vs X
  gXvsY    -> Write("gXvsY"   );
  hXvsY_MC -> Write("hXvsY_MC");
  cXvsY_MC -> Write("cXvsY_MC");
  cXvsY    -> Write("cXvsY"   );
  // Y vs Z
  gZvsY    -> Write("gZvsY"   );
  hZvsY_MC -> Write("hZvsY_MC");
  cZvsY_MC -> Write("cZvsY_MC");
  cZvsY    -> Write("cZvsY"   );
  // X vs Z
  gZvsX    -> Write("gZvsX"   );
  hZvsX_MC -> Write("hZvsX_MC");
  cZvsX_MC -> Write("cZvsX_MC");
  cZvsX    -> Write("cZvsX"   );

  return;
}