pass1_plots.C
Go to the documentation of this file.
1 const char* gPrefix = 0; // Plot files will carry this prefix
2 int gMode = 1; // We are comparing Trident MC to GENIE MC
3 // int gMode = 2; // ................ ND Data to GENIE MC (?)
4 int gPrintPlots = 1; // Should we save plots to .pdf?
5 int gPlotIt = 1; // Turns plots on/off during development
6 
7 //......................................................................
8 
9 TFile* fpSIG = 0;
10 TFile* fpBKG = 0;
11 
12 TTree* ntSIG = 0;
13 TTree* ntBKG = 0;
14 
15 void open_files()
16 {
17  if (gMode==1) {
18  std::cout << "Comparing trident MC to GENIE MC" << std::endl;
19  gPrefix = "p1_T2G";
20  fpSIG = TFile::Open("dimuon_signal_hist.root");
21  fpBKG = TFile::Open("dimuon_background_hist.root");
22 
23  ntSIG = (TTree*)fpSIG->FindObjectAny("p1");
24  ntBKG = (TTree*)fpBKG->FindObjectAny("p1");
25  }
26  else {
27  std::cout << "Ugh oh - this mode isn't supported..." << std::endl;
28  }
29 }
30 
31 //......................................................................
32 //
33 // A lot of the plotting is repetitive. Make some functions to do the
34 // routine stuff and to enforce a good style.
35 //
36 TCanvas* canvas(const char* nm, const char* ti, const char* a)
37 {
38  if (a == TString("square")) return new TCanvas(nm, nm, 600,600);
39  if (a == TString("portrait")) return new TCanvas(nm, nm, 500,700);
40  if (a == TString("landscape")) return new TCanvas(nm, nm, 700,500);
41  return 0;
42 }
43 void area_normalize(TH1F* h1, TH1F* h2)
44 {
45  h2->Scale(h1->GetEntries()/h2->GetEntries());
46 }
47 float plot_max(TH1F* h1, TH1F* h2, float r, int islog=0)
48 {
49  float m1 = h1->GetMaximum();
50  float m2 = h2->GetMaximum();
51  if (islog) {
52  m1 = log10(m1);
53  m2 = log10(m2);
54  }
55 
56  float mx = 0;
57  if (m1>m2) mx = (std::floor(m1/r)+1)*r;
58  else mx = (std::floor(m2/r)+1)*r;
59  if (islog) mx = pow(10.0, mx);
60 
61  h1->SetMaximum(mx);
62  h2->SetMaximum(mx);
63 
64  if (islog) {
65  m1 = h1->GetMinimum();
66  m2 = h2->GetMinimum();
67  if (m1<=1 || m2<=1) h1->SetMinimum(0.1);
68  else h1->SetMinimum(1.0);
69  }
70  else {
71  h1->SetMinimum(0);
72  }
73 
74  return mx;
75 }
76 void show_x_cut(float x, float dx, float y1, float y2, float fy=0.82)
77 {
78  TLine* l = new TLine(x, y1, x, y2);
79  l->SetLineColor(kBlack);
80  l->SetLineStyle(3);
81  l->Draw();
82 
83  double y = y1 + fy*(y2-y1);
84  TArrow* a = new TArrow(x, y, x+dx, y, 0.02);
85 
86  a->SetLineColor(kBlack);
87  a->Draw();
88 }
89 void sig_style(TH1F* h)
90 {
91  h->SetFillColor(kAzure-9);
92  h->GetXaxis()->SetNdivisions(505);
93  h->GetYaxis()->SetNdivisions(505);
94 }
95 void bkg_style(TH1F* h)
96 {
97  h->SetLineColor(kRed);
98  h->SetLineWidth(1);
99 }
100 
101 //......................................................................
102 //
103 // Again, the projecting and plotting of variables is repetitive. Code
104 // up the common features.
105 //
106 void proj_1D(TH1F** h1,
107  TH1F** h2,
108  const char* base_name,
109  const char* title,
110  int n,
111  float x1,
112  float x2,
113  const char* var,
114  TCut cut)
115 {
116  char sig_name[256];
117  char bkg_name[256];
118  sprintf(sig_name, "sig%s", base_name);
119  sprintf(bkg_name, "bkg%s", base_name);
120 
121  TH1F* sig = new TH1F(sig_name, title, n, x1, x2);
122  TH1F* bkg = new TH1F(bkg_name, title, n, x1, x2);
123  (*h1) = sig;
124  (*h2) = bkg;
125 
126  ntSIG->Project(sig_name, var, cut);
127  ntBKG->Project(bkg_name, var, cut);
128 
129  std::cout << "|--* "
130  << base_name << " & "
131  << sig->GetEntries() << " & "
132  << bkg->GetEntries() << " & "
133  << cut << " \\\\"
134  << std::endl;
135 
136  area_normalize(sig, bkg);
137 }
138 
139 //......................................................................
140 
141 void plot_1D(TH1F* h1,
142  TH1F* h2,
143  const char* title,
144  const char* aspect,
145  int islogy,
146  float r,
147  float cut_left=-99,
148  float cut_right=-99)
149 {
150  char ti[256];
151  sprintf(ti, "%s_%s", gPrefix, title);
152  TCanvas* c = canvas(ti, 0, aspect);
153  if (islogy==1) c->SetLogy();
154 
155  float mx = plot_max(h1, h2, r, islogy);
156 
157  sig_style(h1);
158  bkg_style(h2);
159 
160  float xlo = h1->GetXaxis()->GetXmin();
161  float xhi = h1->GetXaxis()->GetXmax();
162  if (fabs((xlo+xhi)/(xhi-xlo))<0.05) {
163  h1->GetXaxis()->CenterTitle();
164  }
165  h1->Draw();
166  h2->Draw("hist,same");
167 
168  float dx = 0.05*(xhi-xlo);
169  if (cut_left !=-99) { show_x_cut(cut_left, dx, 0, mx); }
170  if (cut_right!=-99) { show_x_cut(cut_right, -dx, 0, mx); }
171 
172  if (gPrintPlots) {
173  char fname[256];
174  sprintf(fname,"%s.pdf",c->GetName());
175  c->SaveAs(fname);
176  }
177 }
178 
179 //......................................................................
180 TCut gCutNHIT;
181 TCut gCutVETO;
182 TCut gCutVTX;
185 
186 //......................................................................
187 //
188 // Here's where I walk through the distributions and make some cut
189 // decisions.
190 //
191 // The cuts proceed in roughly two big passes. First, are some basic
192 // quality cuts which ensure that we have events with enough
193 // information to have a good chance at a correct reconstruction and
194 // containment. We really don't have much choice about these. The
195 // second pass is to try to do some signal/background separation based
196 // on the event topology. Here I've chosen to try to be ~99%
197 // efficient (as a fraction of the selections above) while rejecting
198 // as much background as possible.
199 //
200 // Note: If you turn the plotting off, this will set all the global
201 // cuts.
202 //
203 void make_plots()
204 {
205  int liny = 0; // Put the plot on linear y scale
206  int logy = 1; // Put the plot on log y scale
207 
208  /*====================================================================
209  First level of cuts is to make sure we have enough information to
210  work with in the event.
211 
212  We are searching for two muon tracks. Backgrounds are likely to be
213  long muons (rock and cosmics) and events with long pion tracks.
214 
215  The interaction length in the detector is lam_I = ~40 cm so pions
216  are very unlikely to travel more than 3 lam_I ~= 120 cm or about
217  20 planes in the detecor. Rock muons and cosmics may traverse the
218  entire detector length, ~220 planes.
219 
220  If I have two muon tracks traveling across 20 planes I should
221  expect ~40 hits in the event. If the events are oriented along the
222  beam direction, there should be roughly equal number of hits in
223  each view.
224  */
225  TCut cutNONE = TCut("1");
226  if (gPlotIt) {
227  TH1F* sigNpl = 0;
228  TH1F* bkgNpl = 0;
229  proj_1D(&sigNpl, &bkgNpl, "Nplane", ";number of hit planes;slice counts/5 planes", 44, 0, 220, "np", cutNONE);
230  plot_1D( sigNpl, bkgNpl, "01_Nplane", "portrait", logy, 0.5, 20, 200);
231  }
232  TCut cutNp("np>=20&&np<=200");
233 
234  if (gPlotIt) {
235  TH1F* sigNhit = 0;
236  TH1F* bkgNhit = 0;
237  proj_1D(&sigNhit, &bkgNhit, "Nhit", ";number of hit cells;slice counts/5 hits", 80, 0, 400, "(nhitx+nhity)", cutNONE);
238  plot_1D( sigNhit, bkgNhit, "01_Nhit", "portrait", logy, 0.5, 30, -99);
239  }
240  TCut cutNhit("(nhitx+nhity)>=30");
241 
242  if (gPlotIt) {
243  TH1F* sigAnhit = 0;
244  TH1F* bkgAnhit = 0;
245  proj_1D(&sigAnhit, &bkgAnhit, "Anhit", ";(N_{x}-N_{y})/(N_{x}+N_{y});slice counts/0.05", 20, -1, 1, "(nhitx-nhity)/(nhitx+nhity)", cutNp&&cutNhit);
246  plot_1D( sigAnhit, bkgAnhit, "01_Anhit", "portrait", logy, 0.5, -0.5, 0.5);
247  }
248  TCut cutAnhit("fabs(nhitx-nhity)/(nhitx+nhity)<=0.5");
249 
250  if (gPlotIt) {
251  TH1F* sigN2p = 0;
252  TH1F* bkgN2p = 0;
253  proj_1D(&sigN2p, &bkgN2p, "N2p", ";number of planes with 2 hits;slice counts", 100, 0, 100, "n2p", cutNp&&cutNhit&&cutAnhit);
254  plot_1D( sigN2p, bkgN2p, "01_N2p", "portrait", liny, 20, 8, -99);
255  }
256  TCut cutN2p("n2p>=8");
257 
258  if (gPlotIt) {
259  TH1F* sigN3p = 0;
260  TH1F* bkgN3p = 0;
261  proj_1D(&sigN3p, &bkgN3p, "N3p", ";number of planes with 3+ hits;slice counts", 100, 0, 100, "n3p", cutNp&&cutNhit&&cutAnhit);
262  plot_1D( sigN3p, bkgN3p, "01_N3p", "portrait", logy, 0.5, -99, 40);
263  }
264  TCut cutN3p("n3p<=40");
265 
266  //
267  // Put this level of selection into a single cut to be used later.
268  //
269  gCutNHIT = TCut(cutNp&&cutNhit&&cutAnhit&&cutN2p&&cutN3p);
270  //====================================================================
271 
272  /*====================================================================
273  Next look at containment. This done by drawing a buffer zone of 5
274  planes on the front face, and 5 cells on the top, bottom, left,
275  and right and counting how many hits occur in those zones. Then we
276  look where we guessed the event vertex to be and require it to be
277  inside a fiducial volume.
278  */
279  if (gPlotIt) {
280  TH1F* sigNp5 = 0;
281  TH1F* bkgNp5 = 0;
282  proj_1D(&sigNp5, &bkgNp5, "Np5", ";number of hits on front face;slice counts/hit", 51, -0.5, 50.5, "np5", gCutNHIT);
283  plot_1D( sigNp5, bkgNp5, "02_Np5", "portrait", logy, 0.5, -99, 1.5);
284  }
285  TCut cutNp5("np5<=1");
286 
287  if (gPlotIt) {
288  TH1F* sigNwall = 0;
289  TH1F* bkgNwall = 0;
290  proj_1D(&sigNwall, &bkgNwall, "Nwall", ";number of side walls hit;slice counts", 5, -0.5, 4.5, "(nc5t>1)+(nc5b>1)+(nc5l>1)+(nc5r>1)", gCutNHIT);
291  plot_1D( sigNwall, bkgNwall, "02_Nwall", "portrait", logy, 1.0, -99, 2.5);
292  }
293  TCut cutNwall("((nc5t>1)+(nc5b>1)+(nc5l>1)+(nc5r>1))<=2");
294 
295  //
296  // Put this level of selection into a single cut to be used later.
297  //
298  gCutVETO = TCut(cutNp5&&cutNwall);
299  //====================================================================
300 
301  /*====================================================================
302  Select verticies inside a fiducial region.
303  */
304  if (gPlotIt) {
305  TH1F* sigVx = 0;
306  TH1F* bkgVx = 0;
307  proj_1D(&sigVx, &bkgVx, "Vx", ";vertex x (cm);slice counts/ 10 cm", 40, -200, 200, "x", gCutNHIT&&gCutVETO);
308  plot_1D( sigVx, bkgVx, "03_Vx", "portrait", liny, 20, -160, 160);
309 
310  TH1F* sigVy = 0;
311  TH1F* bkgVy = 0;
312  proj_1D(&sigVy, &bkgVy, "Vy", ";vertex y (cm);slice counts/ 10 cm", 40, -200, 200, "y", gCutNHIT&&gCutVETO);
313  plot_1D( sigVy, bkgVy, "03_Vy", "portrait", liny, 20, -160, 160);
314 
315  TH1F* sigVz = 0;
316  TH1F* bkgVz = 0;
317  proj_1D(&sigVz, &bkgVz, "Vz", ";vertex z (cm);slice counts/ 50 cm", 28, 0, 1400, "z", gCutNHIT&&gCutVETO);
318  plot_1D( sigVz, bkgVz, "03_Vz", "landscape", liny, 20, 50, 1000);
319  }
320  //
321  // Summary of vertex cut
322  //
323  gCutVTX = TCut("abs(x)<=160&&abs(y)<=160&&z>=50&&z<=1000");
324  //====================================================================
325 
326  /*====================================================================
327  Make some assessment of the overall shape of the event. An
328  event with two muon tracks and no hadronic activity should have ~2
329  hits per plane consistent along its length. Other neutrino events
330  will be "front loaded" with hadronic activity.
331  */
332  if (gPlotIt) {
333  TH1F* sigHperP = 0;
334  TH1F* bkgHperP = 0;
335  proj_1D(&sigHperP, &bkgHperP, "HperP", ";Nhit/Nplane;slice counts/0.2", 25, 0, 5, "(nhitx+nhity)/np", gCutNHIT&&gCutVETO&&gCutVTX);
336  plot_1D( sigHperP, bkgHperP, "04_HperP", "portrait", liny, 450, -99, 2.5);
337  }
338  TCut cutHperP("(nhitx+nhity)/np<=2.5");
339 
340  if (gPlotIt) {
341  TH1F* sigNc = 0;
342  TH1F* bkgNc = 0;
343  proj_1D(&sigNc, &bkgNc, "NcNp", ";(N^{2}_{cx}+N^{2}_{cy})^{1/2}/N_{planes};slice counts/0.1", 25, 0, 2.5, "sqrt(ncx**2+ncy**2)/np", gCutNHIT&&gCutVETO&&gCutVTX);
344  plot_1D( sigNc, bkgNc, "04_NcNp", "portrait", logy, 0.5, -99, 1.2);
345  }
346  TCut cutNcNp("sqrt(ncx**2+ncy**2)/np<=1.2");
347 
348  TH1F* sigN10 = 0;
349  TH1F* bkgN10 = 0;
350  if (gPlotIt) {
351  proj_1D(&sigN10, &bkgN10, "N10", ";number of hits in first ten planes;slice counts", 61, -0.5, 60.5, "n10", gCutNHIT&&gCutVETO&&gCutVTX);
352  plot_1D( sigN10, bkgN10, "04_N10", "portrait", logy, 0.5, -99, 25.5);
353  }
354  TCut cutN10("n10<=25");
355 
356  if (gPlotIt) {
357  TH1F* sigPks = 0;
358  TH1F* bkgPks = 0;
359  proj_1D(&sigPks, &bkgPks, "Pks", ";log_{10}(K-S probability);slice counts/0.2", 40, -8, 0, "log10(pks)", gCutNHIT&&gCutVETO&&gCutVTX);
360  plot_1D( sigPks, bkgPks, "04_Pks", "portrait", logy, 0.5, -4.5, -99);
361  }
362  TCut cutPks("log10(pks)>=-4.5");
363  //
364  // Put all the shape cuts together...
365  //
366  gCutSHAPE = TCut(cutHperP&&cutNcNp&&cutN3p&&cutN2p&&cutN10&&cutPks);
367  //====================================================================
368 
369  //
370  // Make a plot of the N plane distribution we started with now with
371  // all cuts just to close out the counts
372  //
374  if (gPlotIt) {
375  TH1F* sigNpFinal = 0;
376  TH1F* bkgNpFinal = 0;
377  proj_1D(&sigNpFinal, &bkgNpFinal, "NpFinal", ";number of hit planes;slice counts/5 planes", 44, 0, 220, "np", gCutPASS1);
378  plot_1D( sigNpFinal, bkgNpFinal, "05_NpFinal", "portrait", logy, 1.0, 20, 200);
379  }
380 }
381 
382 //......................................................................
383 //
384 // For each level of cuts, print out some events numbers and relevant
385 // information for signal and background events that pass and
386 // fail. Use these to check that things are working as we expect by
387 // looking at the events in the event display.
388 //
389 void scan(const char* sample, // Signal or background?
390  const char* cut_name, // What cut are we inspecting?
391  TTree* nt, // Which ntuple?
392  int n0, // Where to start the scan
393  int n1, // How many entries to scan for "passed"
394  int n2, // How many entries to scan for "failed"
395  const char* var, // Variables to print
396  TCut c, // Base set of cuts
397  TCut s) // The cut to toggle on and off
398 {
399  //
400  // Show all rows, no prompt to continue, during Scan()
401  //
402  nt->SetScanField(0);
403 
404  std::cout << "|== " << sample << " passing " << cut_name << std::endl;
405  nt->Scan(var, c&&s, "", n1, n0);
406 
407  std::cout << "|== " << sample << " failing " << cut_name << std::endl;
408  nt->Scan(var, c&&!s, "", n2, n0);
409 }
410 void nt_scan()
411 {
412  TCut noCut("1");
413 
414  scan("SIGNAL", "NHIT", ntSIG, 0, 50, 50, "run:event:slice:np:nhitx:nhity:n2p:n3p", noCut, gCutNHIT);
415  scan("BACKGROUND", "NHIT", ntBKG, 0, 50, 50, "run:event:slice:np:nhitx:nhity:n2p:n3p", noCut, gCutNHIT);
416 
417  scan("SIGNAL", "VETO", ntSIG, 50, 30, 30, "run:event:slice:np5:nc5t:nc5b:nc5l:nc5r", gCutNHIT, gCutVETO);
418  scan("BACKGROUND", "VETO", ntBKG, 50, 400, 30, "run:event:slice:np5:nc5t:nc5b:nc5l:nc5r", gCutNHIT, gCutVETO);
419 
420  scan("SIGNAL", "VTX", ntSIG, 500, 40, 400, "run:event:slice:x:y:z", gCutNHIT&&gCutVETO, gCutVTX);
421  scan("BACKGROUND", "VTX", ntBKG, 500, 400, 1000, "run:event:slice:x:y:z", gCutNHIT&&gCutVETO, gCutVTX);
422 
423  scan("SIGNAL", "SHAPE", ntSIG, 1000, 40, 5000, "run:event:slice:nhitx:nhity:np:ncx:ncy:log10(pks)", gCutNHIT&&gCutVETO&&gCutVTX, gCutSHAPE);
424  scan("BACKGROUND", "SHAPE", ntBKG, 1000, 1000, 2000, "run:event:slice:nhitx:nhity:np:ncx:ncy:log10(pks)", gCutNHIT&&gCutVETO&&gCutVTX, gCutSHAPE);
425 }
426 
427 //......................................................................
428 
430 {
431  gStyle->SetOptStat(0);
432  gErrorIgnoreLevel = kWarning;
433 
434  open_files();
435  make_plots();
436  nt_scan();
437 }
438 
439 ////////////////////////////////////////////////////////////////////////
TCut gCutVETO
Definition: pass1_plots.C:181
enum BeamMode kRed
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
TCut gCutSHAPE
Definition: pass1_plots.C:183
TCut gCutPASS1
Definition: pass1_plots.C:184
Float_t x1[n_points_granero]
Definition: compare.C:5
constexpr T pow(T x)
Definition: pow.h:75
static constexpr Double_t nm
Definition: Munits.h:133
TTree * ntSIG
Definition: pass1_plots.C:12
TCut gCutNHIT
Definition: pass1_plots.C:180
TCanvas * canvas(const char *nm, const char *ti, const char *a)
Definition: pass1_plots.C:36
bool logy
void make_plots()
Definition: pass1_plots.C:203
const XML_Char * s
Definition: expat.h:262
int gPrintPlots
Definition: pass1_plots.C:4
const char * gPrefix
Definition: pass1_plots.C:1
void bkg_style(TH1F *h)
Definition: pass1_plots.C:95
TCut gCutVTX
Definition: pass1_plots.C:182
double dx[NP][NC]
void proj_1D(TH1F **h1, TH1F **h2, const char *base_name, const char *title, int n, float x1, float x2, const char *var, TCut cut)
Definition: pass1_plots.C:106
float plot_max(TH1F *h1, TH1F *h2, float r, int islog=0)
Definition: pass1_plots.C:47
const double a
TFile * fpBKG
Definition: pass1_plots.C:10
int gPlotIt
Definition: pass1_plots.C:5
void open_files()
Definition: pass1_plots.C:15
void sig_style(TH1F *h)
Definition: pass1_plots.C:89
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
void area_normalize(TH1F *h1, TH1F *h2)
Definition: pass1_plots.C:43
void pass1_plots()
Definition: pass1_plots.C:429
TH1F * h1
static constexpr Double_t m2
Definition: Munits.h:145
const Cut cut
Definition: exporter_fd.C:30
T log10(T number)
Definition: d0nt_math.hpp:120
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
TRandom3 r(0)
int gMode
Definition: pass1_plots.C:2
void nt_scan()
Definition: pass1_plots.C:410
TTree * ntBKG
Definition: pass1_plots.C:13
void scan(const char *sample, const char *cut_name, TTree *nt, int n0, int n1, int n2, const char *var, TCut c, TCut s)
Definition: pass1_plots.C:389
void plot_1D(TH1F *h1, TH1F *h2, const char *title, const char *aspect, int islogy, float r, float cut_left=-99, float cut_right=-99)
Definition: pass1_plots.C:141
TNtuple * nt
Definition: drawXsec.C:2
void show_x_cut(float x, float dx, float y1, float y2, float fy=0.82)
Definition: pass1_plots.C:76
TFile * fpSIG
Definition: pass1_plots.C:9