PositionComparison.C
Go to the documentation of this file.
1 //..................................
3 //#include "CAFAna/Analysis/Plots.h"
5 
6 #include "CAFAna/Cuts/Cuts.h"
10 
11 #include "CAFAna/Core/Binning.h"
12 #include "CAFAna/Core/Loaders.h"
14 #include "CAFAna/Core/Var.h"
15 #include "CAFAna/Core/Spectrum.h"
17 
19 #include "CAFAna/Vars/HistAxes.h"
22 
23 #include "Utilities/rootlogon.C"
24 
25 #include "TAttLine.h"
26 #include "TAttMarker.h"
27 #include "TCanvas.h"
28 #include "TColor.h"
29 #include "TFile.h"
30 #include "TGraph.h"
31 #include "TH1.h"
32 #include "TH2.h"
33 #include "TLatex.h"
34 #include "TLegend.h"
35 #include "TLine.h"
36 #include "TStyle.h"
37 
38 #include <cmath>
39 #include <iostream>
40 #include <fstream>
41 #include <vector>
42 #include <string>
43 #include <TROOT.h>
44 #include <TStyle.h>
45 //..................................
46 using namespace ana;
47 
48 std::string PlotDir = "plots_position/";
49 
50 //..................................
51 void MakeTextFile( std::string Plane, std::string Overlay ) {
52  // Determine output name, and caption.
53  std::string FNa = PlotDir+Plane+"_"+Overlay+".txt";
54  std::string Cap = "Plot showing the distribution of NuMi data events passing the NuMu FD selection in the "+Plane+" plane, with "+Overlay;
55  // Write to file.
56  std::ofstream TxtOut ( FNa.c_str(), std::ofstream::out );
57  TxtOut << Cap;
58  TxtOut.close();
59  // Done.
60  return;
61 }
62 
63 //..................................
64 void PositionComparison( bool Test=false )
65 {
66  // ---- First off, lets set the styles...
67  gROOT->SetStyle("novaStyle");
68  //gStyle->SetPalette(87); // https://root.cern.ch/doc/master/classTColor.html
69 
70  // --- Define my definitions
71  std::string AllData = "/pnfs/nova/persistent/users/bzamoran/CheckEvents/fd_data.root";
72  std::string AllMC = "prod_sumdecaf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1_numu2017";
73  if (Test) AllMC = "prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1";
74 
75  // --- Define my binnnings.
76  const Binning XBins = Binning::Simple(20, -10, 10 );
77  const Binning YBins = Binning::Simple(20, -10, 10 );
78  const Binning ZBins = Binning::Simple(31, -2 , 60);
79 
80  // --- Define my loaders
81  SpectrumLoader allswap_MC( AllMC );
82 
83  // --- Define my Spill cuts
84  allswap_MC.SetSpillCut(kStandardSpillCuts);
85 
86  // --- Define my spectra.
87  // X vs Y
88  Spectrum* XvsY_MC = new Spectrum("Positions in X vs Y (m)", allswap_MC, XBins, kTrkStartX, YBins, kTrkStartY, kNumuCutFD2017, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
89  // Y vs Z
90  Spectrum* ZvsY_MC = new Spectrum("Positions in Z vs Y (m)", allswap_MC, ZBins, kTrkStartZ, YBins, kTrkStartY, kNumuCutFD2017, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
91  // X Z vs X
92  Spectrum* ZvsX_MC = new Spectrum("Positions in Z vs X (m)", allswap_MC, ZBins, kTrkStartZ, XBins, kTrkStartX, kNumuCutFD2017, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
93 
94  allswap_MC.Go();
95 
96  // .............................. Y vs X ..............................
97  TH2* hXvsY_MC = XvsY_MC -> ToTH2(kAna2017POT);
98  hXvsY_MC->GetXaxis()->SetTitle("Vertex X position (m)");
99  hXvsY_MC->GetYaxis()->SetTitle("Vertex Y position (m)");
100  hXvsY_MC->GetYaxis()->SetTitleOffset(0.7);
101  hXvsY_MC->GetZaxis()->SetTitle("Arbitrary units");
102  CenterTitles(hXvsY_MC);
103 
105  gXvsY->SetMarkerStyle(20);
106  gXvsY->SetMarkerSize(2.5);
107  gXvsY->GetXaxis()->SetTitle("Vertex X position (m)");
108  gXvsY->GetXaxis()->CenterTitle();
109  gXvsY->GetYaxis()->SetTitle("Vertex Y position (m)");
110  gXvsY->GetYaxis()->SetTitleOffset(0.7);
111  gXvsY->GetYaxis()->CenterTitle();
112 
113  // First draw canvas with data only.
114  TCanvas* cXvsY = new TCanvas("cXvsY","X pos vs Y pos no MC", 900, 800);
115  gXvsY -> Draw("AP");
116  Preliminary();
117  // Define my lines.
118  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");
119  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");
120  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");
121  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");
122  gXvsY -> GetXaxis()->SetRangeUser(-10, 10);
123  gXvsY -> GetYaxis()->SetRangeUser(-10, 10);
124  cXvsY -> SaveAs( TString(PlotDir)+TString("XvsY_noMCOverlay.pdf") );
125  MakeTextFile( "XvsY", "noMCOverlay" );
126 
127  // Scale histogram
128  double max1 = hXvsY_MC->GetBinContent(hXvsY_MC->GetMaximumBin());
129  hXvsY_MC->Scale(1./max1);
130  // Set range
131  hXvsY_MC->GetXaxis()->SetRangeUser(-8, 8);
132  hXvsY_MC->GetYaxis()->SetRangeUser(-8, 7);
133  // Now draw the canvas with MC background, and zoom in on detector.
134  TCanvas* cXvsY_MC = new TCanvas("cXvsY_MC","X pos vs Y pos MC Overlay", 900, 800);
135  cXvsY_MC -> SetRightMargin (.15);
136  hXvsY_MC -> Draw("COLZ");
137  gXvsY -> Draw("P");
138  Preliminary();
139  cXvsY_MC -> SaveAs( TString(PlotDir)+TString("XvsY_MCOverlay.pdf") );
140  MakeTextFile( "XvsY", "MCOverlay" );
141 
142  // .............................. Z vs Y ..............................
143  TH2* hZvsY_MC = ZvsY_MC->ToTH2(kAna2017POT);
144  hZvsY_MC->GetYaxis()->SetTitle("Vertex Y position (m)");
145  hZvsY_MC->GetYaxis()->SetTitleOffset(0.5);
146  hZvsY_MC->GetXaxis()->SetTitle("Vertex Z position (m)");
147  hZvsY_MC->GetZaxis()->SetTitle("Arbitrary units");
148  CenterTitles(hZvsY_MC);
149 
151  gZvsY->SetMarkerStyle(20);
152  gZvsY->SetMarkerSize(2.5);
153  gZvsY->GetXaxis()->SetTitle("Vertex Z position (m)");
154  gZvsY->GetXaxis()->CenterTitle();
155  gZvsY->GetYaxis()->SetTitle("Vertex Y position (m)");
156  gZvsY->GetYaxis()->SetTitleOffset(0.5);
157  gZvsY->GetYaxis()->CenterTitle();
158 
159  // First draw canvas with data only.
160  TCanvas* cZvsY = new TCanvas("cZvsY","Z pos vs Y pos no MC", 2500, 800);
161  gZvsY -> Draw("AP");
162  Preliminary();
163  // Define my lines.
164  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");
165  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");
166  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");
167  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");
168  gZvsY -> GetXaxis()->SetRangeUser( -5 , 60 ); // Z Pos
169  gZvsY -> GetYaxis()->SetRangeUser( -10, 10 ); // Y Pos
170  cZvsY -> SaveAs( TString(PlotDir)+TString("ZvsY_noMCOverlay.pdf") );
171  MakeTextFile( "ZvsY", "noMCOverlay" );
172 
173  // Scale the histogram
174  double max2 = hZvsY_MC->GetBinContent(hZvsY_MC->GetMaximumBin());
175  hZvsY_MC->Scale(1./max2);
176  // Set the range
177  hZvsY_MC->GetXaxis()->SetRangeUser( 0 , 58 ); // Z Pos
178  hZvsY_MC->GetYaxis()->SetRangeUser( -8, 7 ); // Y Pos
179  // Now draw the canvas with MC background, and zoom in on detector.
180  TCanvas* cZvsY_MC = new TCanvas("cZvsY_MC","Z pos vs Y pos MC Overlay", 2500, 800);
181  cZvsY_MC -> SetRightMargin (.15);
182  hZvsY_MC -> Draw("COLZ");
183  gZvsY -> Draw("P");
184  Preliminary();
185  cZvsY_MC -> SaveAs( TString(PlotDir)+TString("ZvsY_MCOverlay.pdf") );
186  MakeTextFile( "ZvsY", "MCOverlay" );
187 
188  // .............................. Z vs X ..............................
189  TH2* hZvsX_MC = ZvsX_MC->ToTH2(kAna2017POT);
190  hZvsX_MC->GetXaxis()->SetTitle("Vertex Z position (m)");
191  hZvsX_MC->GetYaxis()->SetTitle("Vertex X position (m)");
192  hZvsX_MC->GetYaxis()->SetTitleOffset(0.5);
193  hZvsX_MC->GetZaxis()->SetTitle("Arbitrary units");
194  CenterTitles(hZvsX_MC);
195 
197  gZvsX->SetMarkerStyle(20);
198  gZvsX->SetMarkerSize(2.5);
199  gZvsX->GetXaxis()->SetTitle("Vertex Z position (m)");
200  gZvsX->GetXaxis()->CenterTitle();
201  gZvsX->GetYaxis()->SetTitle("Vertex X position (m)");
202  gZvsX->GetYaxis()->SetTitleOffset(0.5);
203 
204  // First draw canvas with data only.
205  TCanvas* cZvsX = new TCanvas("cZvsX","Z pos vs X pos no MC", 2500, 800);
206  gZvsX -> Draw("AP");
207  Preliminary();
208  // Define my lines.
209  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");
210  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");
211  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");
212  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");
213  gZvsX -> GetXaxis()->SetRangeUser( -5 , 60 ); // Z Pos
214  gZvsX -> GetYaxis()->SetRangeUser( -10, 10 ); // Y Pos
215  cZvsX -> SaveAs( TString(PlotDir)+TString("ZvsX_noMCOverlay.pdf") );
216  MakeTextFile( "ZvsX", "noMCOverlay" );
217 
218  // Scale the histogram
219  double max3 = hZvsX_MC->GetBinContent(hZvsX_MC->GetMaximumBin());
220  hZvsX_MC->Scale(1./max3);
221  // Set the range
222  hZvsX_MC->GetXaxis()->SetRangeUser( 0 , 58 ); // Z Pos
223  hZvsX_MC->GetYaxis()->SetRangeUser( -8, 7 ); // Y Pos
224  // Now draw the canvas with MC background, and zoom in on detector.
225  TCanvas* cZvsX_MC = new TCanvas("cZvsX_MC","Z pos vs X pos MC Overlay", 2500, 800);
226  cZvsX_MC -> SetRightMargin (.15);
227  hZvsX_MC -> Draw("COLZ");
228  gZvsX -> Draw("P");
229  Preliminary();
230  cZvsX_MC -> SaveAs( TString(PlotDir)+TString("ZvsX_MCOverlay.pdf") );
231  MakeTextFile( "ZvsX", "MCOverlay" );
232 
233  // .............................. Define an outfile ..............................
234  std::string OutName = PlotDir+"PosComparisons.root";
235  TFile *OutFile = new TFile(OutName.c_str(), "RECREATE");
236  OutFile->cd();
237  // Y vs X
238  gXvsY -> Write("gXvsY" );
239  hXvsY_MC -> Write("hXvsY_MC");
240  cXvsY_MC -> Write("cXvsY_MC");
241  cXvsY -> Write("cXvsY" );
242  // Y vs Z
243  gZvsY -> Write("gZvsY" );
244  hZvsY_MC -> Write("hZvsY_MC");
245  cZvsY_MC -> Write("cZvsY_MC");
246  cZvsY -> Write("cZvsY" );
247  // X vs Z
248  gZvsX -> Write("gZvsX" );
249  hZvsX_MC -> Write("hZvsX_MC");
250  cZvsX_MC -> Write("cZvsX_MC");
251  cZvsX -> Write("cZvsX" );
252 
253  return;
254 }
tree Draw("slc.nhit")
enum BeamMode kRed
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
correl_xv GetYaxis() -> SetDecimals()
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kNumuCutFD2017
Definition: NumuCuts2017.h:39
void MakeTextFile(std::string Plane, std::string Overlay)
const Var kTrkStartY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Y()/100;})
Definition: NumuVars.h:52
const double kAna2017POT
Definition: Exposures.h:174
void SetSpillCut(const SpillCut &cut)
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1483
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void PositionComparison(bool Test=false)
correl_xv GetXaxis() -> SetDecimals()
const Var kTrkStartZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Z()/100;})
Definition: NumuVars.h:53
std::string PlotDir
virtual void Go() override
Load all the registered spectra.
std::vector< float > Spectrum
Definition: Constants.h:759
void Preliminary()
const Var kTrkStartX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.X()/100;})
Definition: NumuVars.h:51
const SystShifts kNoShift
Definition: SystShifts.cxx:22
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
Definition: TimingCuts.h:8
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
TFile * OutFile
const Var kXSecCVWgt2017
Definition: XsecTunes.h:36
cosmicTree SaveAs("cosmicTree.root")
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
TGraph * MakeGraph(const std::string &wildcard, const Cut &cut, const Var &VarX, const Var &VarY, const SpillCut *spillCut)
Get a graph with two variables when a histogram is not the best choice.
Definition: GraphDrawer.cxx:13
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: UtilsExt.cxx:115
gm Write()
enum BeamMode string