FD_Data_PosComp.C
Go to the documentation of this file.
1 //..................................
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 FD_Data_PosComp( bool IsFHC = true )
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 = "prod4_sumrestricteddecaf_fd_numi_fhc_full_goodruns_numu2018";
72  std::string AllMC = "prod_sumdecaf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1_numu2018_str5";
73  std::string sFHCRHC = "FHC_";
74  double MyPOT = kAna2018FHCPOT;
75  if (!IsFHC) {
76  AllData = "prod4_sumrestricteddecaf_fd_numi_rhc_full_goodruns_numu2018";
77  AllMC = "prod_sumdecaf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1_numu2018_str5";
78  sFHCRHC = "RHC_";
79  MyPOT = kAna2018RHCPOT;
80  }
81 
82  PlotDir += sFHCRHC;
83 
84  // --- Define my binnnings.
85  const Binning XBins = Binning::Simple(20, -10, 10 );
86  const Binning YBins = Binning::Simple(20, -10, 10 );
87  const Binning ZBins = Binning::Simple(31, -2 , 60);
88 
89  // --- Define my loaders
90  SpectrumLoader allswap_MC( AllMC );
91 
92  // --- Define my Spill cuts
93  allswap_MC.SetSpillCut(kStandardSpillCuts);
94 
95  // --- Define my spectra.
96  // X vs Y
97  Spectrum* XvsY_MC = new Spectrum("Positions in X vs Y (m)", allswap_MC, XBins, kTrkStartX, YBins, kTrkStartY, kNumuCutFD2018, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
98  // Y vs Z
99  Spectrum* ZvsY_MC = new Spectrum("Positions in Z vs Y (m)", allswap_MC, ZBins, kTrkStartZ, YBins, kTrkStartY, kNumuCutFD2018, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
100  // X Z vs X
101  Spectrum* ZvsX_MC = new Spectrum("Positions in Z vs X (m)", allswap_MC, ZBins, kTrkStartZ, XBins, kTrkStartX, kNumuCutFD2018, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018);
102 
103  allswap_MC.Go();
104 
105  // .............................. Y vs X ..............................
106  TH2* hXvsY_MC = XvsY_MC -> ToTH2(MyPOT);
107  hXvsY_MC->GetXaxis()->SetTitle("Vertex X position (m)");
108  hXvsY_MC->GetYaxis()->SetTitle("Vertex Y position (m)");
109  hXvsY_MC->GetYaxis()->SetTitleOffset(0.7);
110  hXvsY_MC->GetZaxis()->SetTitle("Arbitrary units");
111  CenterTitles(hXvsY_MC);
112 
114  gXvsY->SetMarkerStyle(20);
115  gXvsY->SetMarkerSize(2.5);
116  gXvsY->GetXaxis()->SetTitle("Vertex X position (m)");
117  gXvsY->GetXaxis()->CenterTitle();
118  gXvsY->GetYaxis()->SetTitle("Vertex Y position (m)");
119  gXvsY->GetYaxis()->SetTitleOffset(0.7);
120  gXvsY->GetYaxis()->CenterTitle();
121 
122  // First draw canvas with data only.
123  TCanvas* cXvsY = new TCanvas("cXvsY","X pos vs Y pos no MC", 900, 800);
124  gXvsY -> Draw("AP");
125  Preliminary();
126  // Define my lines.
127  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");
128  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");
129  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");
130  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");
131  gXvsY -> GetXaxis()->SetRangeUser(-10, 10);
132  gXvsY -> GetYaxis()->SetRangeUser(-10, 10);
133  cXvsY -> SaveAs( TString(PlotDir)+TString("XvsY_noMCOverlay.pdf") );
134  MakeTextFile( "XvsY", "noMCOverlay" );
135 
136  // Scale histogram
137  double max1 = hXvsY_MC->GetBinContent(hXvsY_MC->GetMaximumBin());
138  hXvsY_MC->Scale(1./max1);
139  // Set range
140  hXvsY_MC->GetXaxis()->SetRangeUser(-8, 8);
141  hXvsY_MC->GetYaxis()->SetRangeUser(-8, 7);
142  // Now draw the canvas with MC background, and zoom in on detector.
143  TCanvas* cXvsY_MC = new TCanvas("cXvsY_MC","X pos vs Y pos MC Overlay", 900, 800);
144  cXvsY_MC -> SetRightMargin (.15);
145  hXvsY_MC -> Draw("COLZ");
146  gXvsY -> Draw("P");
147  Preliminary();
148  cXvsY_MC -> SaveAs( TString(PlotDir)+TString("XvsY_MCOverlay.pdf") );
149  MakeTextFile( "XvsY", "MCOverlay" );
150 
151  // .............................. Z vs Y ..............................
152  TH2* hZvsY_MC = ZvsY_MC->ToTH2(MyPOT);
153  hZvsY_MC->GetYaxis()->SetTitle("Vertex Y position (m)");
154  hZvsY_MC->GetYaxis()->SetTitleOffset(0.5);
155  hZvsY_MC->GetXaxis()->SetTitle("Vertex Z position (m)");
156  hZvsY_MC->GetZaxis()->SetTitle("Arbitrary units");
157  CenterTitles(hZvsY_MC);
158 
160  gZvsY->SetMarkerStyle(20);
161  gZvsY->SetMarkerSize(2.5);
162  gZvsY->GetXaxis()->SetTitle("Vertex Z position (m)");
163  gZvsY->GetXaxis()->CenterTitle();
164  gZvsY->GetYaxis()->SetTitle("Vertex Y position (m)");
165  gZvsY->GetYaxis()->SetTitleOffset(0.5);
166  gZvsY->GetYaxis()->CenterTitle();
167 
168  // First draw canvas with data only.
169  TCanvas* cZvsY = new TCanvas("cZvsY","Z pos vs Y pos no MC", 2500, 800);
170  gZvsY -> Draw("AP");
171  Preliminary();
172  // Define my lines.
173  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");
174  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");
175  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");
176  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");
177  gZvsY -> GetXaxis()->SetRangeUser( -5 , 60 ); // Z Pos
178  gZvsY -> GetYaxis()->SetRangeUser( -10, 10 ); // Y Pos
179  cZvsY -> SaveAs( TString(PlotDir)+TString("ZvsY_noMCOverlay.pdf") );
180  MakeTextFile( "ZvsY", "noMCOverlay" );
181 
182  // Scale the histogram
183  double max2 = hZvsY_MC->GetBinContent(hZvsY_MC->GetMaximumBin());
184  hZvsY_MC->Scale(1./max2);
185  // Set the range
186  hZvsY_MC->GetXaxis()->SetRangeUser( 0 , 58 ); // Z Pos
187  hZvsY_MC->GetYaxis()->SetRangeUser( -8, 7 ); // Y Pos
188  // Now draw the canvas with MC background, and zoom in on detector.
189  TCanvas* cZvsY_MC = new TCanvas("cZvsY_MC","Z pos vs Y pos MC Overlay", 2500, 800);
190  cZvsY_MC -> SetRightMargin (.15);
191  hZvsY_MC -> Draw("COLZ");
192  gZvsY -> Draw("P");
193  Preliminary();
194  cZvsY_MC -> SaveAs( TString(PlotDir)+TString("ZvsY_MCOverlay.pdf") );
195  MakeTextFile( "ZvsY", "MCOverlay" );
196 
197  // .............................. Z vs X ..............................
198  TH2* hZvsX_MC = ZvsX_MC->ToTH2(MyPOT);
199  hZvsX_MC->GetXaxis()->SetTitle("Vertex Z position (m)");
200  hZvsX_MC->GetYaxis()->SetTitle("Vertex X position (m)");
201  hZvsX_MC->GetYaxis()->SetTitleOffset(0.5);
202  hZvsX_MC->GetZaxis()->SetTitle("Arbitrary units");
203  CenterTitles(hZvsX_MC);
204 
206  gZvsX->SetMarkerStyle(20);
207  gZvsX->SetMarkerSize(2.5);
208  gZvsX->GetXaxis()->SetTitle("Vertex Z position (m)");
209  gZvsX->GetXaxis()->CenterTitle();
210  gZvsX->GetYaxis()->SetTitle("Vertex X position (m)");
211  gZvsX->GetYaxis()->SetTitleOffset(0.5);
212 
213  // First draw canvas with data only.
214  TCanvas* cZvsX = new TCanvas("cZvsX","Z pos vs X pos no MC", 2500, 800);
215  gZvsX -> Draw("AP");
216  Preliminary();
217  // Define my lines.
218  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");
219  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");
220  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");
221  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");
222  gZvsX -> GetXaxis()->SetRangeUser( -5 , 60 ); // Z Pos
223  gZvsX -> GetYaxis()->SetRangeUser( -10, 10 ); // Y Pos
224  cZvsX -> SaveAs( TString(PlotDir)+TString("ZvsX_noMCOverlay.pdf") );
225  MakeTextFile( "ZvsX", "noMCOverlay" );
226 
227  // Scale the histogram
228  double max3 = hZvsX_MC->GetBinContent(hZvsX_MC->GetMaximumBin());
229  hZvsX_MC->Scale(1./max3);
230  // Set the range
231  hZvsX_MC->GetXaxis()->SetRangeUser( 0 , 58 ); // Z Pos
232  hZvsX_MC->GetYaxis()->SetRangeUser( -8, 7 ); // Y Pos
233  // Now draw the canvas with MC background, and zoom in on detector.
234  TCanvas* cZvsX_MC = new TCanvas("cZvsX_MC","Z pos vs X pos MC Overlay", 2500, 800);
235  cZvsX_MC -> SetRightMargin (.15);
236  hZvsX_MC -> Draw("COLZ");
237  gZvsX -> Draw("P");
238  Preliminary();
239  cZvsX_MC -> SaveAs( TString(PlotDir)+TString("ZvsX_MCOverlay.pdf") );
240  MakeTextFile( "ZvsX", "MCOverlay" );
241 
242  // .............................. Define an outfile ..............................
243  std::string OutName = PlotDir+"PosComparisons.root";
244  TFile *OutFile = new TFile(OutName.c_str(), "RECREATE");
245  OutFile->cd();
246  // Y vs X
247  gXvsY -> Write("gXvsY" );
248  hXvsY_MC -> Write("hXvsY_MC");
249  cXvsY_MC -> Write("cXvsY_MC");
250  cXvsY -> Write("cXvsY" );
251  // Y vs Z
252  gZvsY -> Write("gZvsY" );
253  hZvsY_MC -> Write("hZvsY_MC");
254  cZvsY_MC -> Write("cZvsY_MC");
255  cZvsY -> Write("cZvsY" );
256  // X vs Z
257  gZvsX -> Write("gZvsX" );
258  hZvsX_MC -> Write("hZvsX_MC");
259  cZvsX_MC -> Write("cZvsX_MC");
260  cZvsX -> Write("cZvsX" );
261 
262  return;
263 }
void FD_Data_PosComp(bool IsFHC=true)
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
bool IsFHC
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
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
const double kAna2018RHCPOT
Definition: Exposures.h:208
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
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 Var kXSecCVWgt2018
Definition: XsecTunes.h:48
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 Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
TFile * OutFile
const double kAna2018FHCPOT
Definition: Exposures.h:207
cosmicTree SaveAs("cosmicTree.root")
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
std::string PlotDir
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