tex_tables_util.h
Go to the documentation of this file.
1 #include "TStyle.h"
2 #include "TFile.h"
3 #include "TH1.h"
4 #include "TCanvas.h"
5 #include "TGraph.h"
6 #include "TGraphAsymmErrors.h"
7 #include "TGaxis.h"
8 #include "TLatex.h"
9 #include "TLegend.h"
10 #include "TLine.h"
11 #include "TSystem.h"
12 
13 #include <iostream>
14 #include <fstream>
15 #include <iomanip>
16 
18 
19 namespace ana
20 {
21 
22  double CalcChi2(TH1* hmc, TH1* hdata)
23  {
24  double chi2 = 0;
25  for (int i = 1; i < hmc->GetNbinsX() + 1; ++i) {
26  double ei = hmc->GetBinContent(i);
27  double oi = hdata->GetBinContent(i);
28  double sigma = sqrt(ei);
29  if(ei>0) chi2 += util::sqr(ei - oi) / util::sqr(sigma);
30  }
31  return chi2;
32  }
33 
34 TString MakeLatexCommandName(TString mystring){
35  //There must be a better way
36  std::vector<TString> in = {"0","1","2","3","4","5",
37  "6","7","8","9","_","#"," ",".","-","(",")"};
38  std::vector<TString> out = {"zero","one","two","three","four","five",
39  "six","seven","eight","nine","","","","","XX","",""};
40  for(unsigned int i=0;i<in.size();i++)
41  mystring.ReplaceAll(in[i],out[i]);
42  return mystring;
43 
44 }
45 
46 TString FixPlotName(TString mystring){
47  //There must be a better way
48  std::vector<TString> in = {"#"," ",".","_"};
49  std::vector<TString> out = {"","","",""};
50  for(unsigned int i=0;i<in.size();i++)
51  mystring.ReplaceAll(in[i],out[i]);
52  return mystring;
53 
54 }
55 
56 
57 void PrintLatexFigure(TString name, TString plotdir="plots/"){
58  output << " \\includegraphics[width=0.49\\textwidth]{"<<plotdir<<name<<".png}\n";
59 }
60 
61 
62 
63 void PrintTableHeader(TString ltxcommand="",
64  TString labelcol_style= "l",
65  TString bincol_style = "| r r r",
66  int nbins = 1 )
67 {
68  if(ltxcommand == "") { std::cerr << "WARNING Problem with PrintTableHeader ltxcommand\n"; return;}
69 
70  output << "\n % " << std::string(60,'-');
71  output << "\n \\vspace{30 pt} \n";
72  output << "\\setlength{\\tabcolsep}{1pt} \\renewcommand{\\arraystretch}{1.2} \n";
73  output << "\\" << MakeLatexCommandName(ltxcommand) << "\n\n";
74 
75  tables << "\\newcommand{\\" << MakeLatexCommandName(ltxcommand) << "}\n";
76  tables << "{\\centerline{ \n \\begin{tabular}{ "
77  << labelcol_style << " ";
78  for (int n =0; n<nbins; ++n) tables << bincol_style << (n +1 < nbins ? " | " : "}\n");
79 }
80 
81 void PrintTableFooter(){
82  tables << "\\end{tabular} } \n } \n"; //end centerline end command
83 }
84 
85 struct SystResults{
86  TString syst;
87  //plus minus, sig bkg
88  // First index --> plus/minus
89  // Second index
90  // 0 = bin1
91  // 1 = bin2
92  // 2 = bin3
93  // 3 = bin4 (Numu ONLY --> 4 quantiles, Nue 2 CVN + 1 peripheral bins)
94  // 4 = Total
95  double diff_sig[2][5];
96  double diff_bkg[2][5];
97  double diff_tot[2][5];
98  double chisq_sig[2];
99  double chisq_bkg[2];
100  double chisq_tot[2];
101 };
102 
103 bool sort_chisq_sig (SystResults a, SystResults b) {
104  double aa = max (a.chisq_sig[0], a.chisq_sig[1]);
105  double bb = max (b.chisq_sig[0], b.chisq_sig[1]);
106  return aa > bb;
107 }
108 
109 bool sort_chisq_bkg (SystResults a, SystResults b) {
110  double aa = max (a.chisq_bkg[0], a.chisq_bkg[1]);
111  double bb = max (b.chisq_bkg[0], b.chisq_bkg[1]);
112  return aa > bb;
113 }
114 
116  double aa = max (a.chisq_tot[0], a.chisq_tot[1]);
117  double bb = max (b.chisq_tot[0], b.chisq_tot[1]);
118  return aa > bb;
119 }
120 
121 // std::vector< std::pair<TString,TH1D*> > allQhist_xp_wtv;
122 std::vector <SystResults> systresults_xp_prop;
123 std::vector <SystResults> systresults_xp_combo;
124 std::vector <SystResults> systresults_nxp;
125 
126 TH1D* systsum_plu_xp_prop[6] = {0,0,0,0,0,0};
127 TH1D* systsum_plu_xp_combo[6] = {0,0,0,0,0,0};
128 TH1D* systsum_plu_nxp[6] = {0,0,0,0,0,0};
129 TH1D* systsum_min_xp_prop[6] = {0,0,0,0,0,0};
130 TH1D* systsum_min_xp_combo[6] = {0,0,0,0,0,0};
131 TH1D* systsum_min_nxp[6] = {0,0,0,0,0,0};
132 
133 /*
134 void AddToEnvelop(TString ltxcommand, TString systName,
135  const std::vector<TH1*> & mcnom, const std::vector<TH1*> & mcplu, const std::vector<TH1*> & mcmin){
136 
137  if ( ltxcommand.Contains("ND") ) return;
138  if ( !(ltxcommand.Contains("EnergyCVN2D")) ) return;
139  if ( ltxcommand.Contains("3bin") ) return;
140 
141  int Nbins = mcnom[0]->GetNbinsX();
142  double binW = mcnom[0]->GetBinWidth(1);
143 
144  std::vector<TString> labels = {"$\\nu_e$ signal","Tot beam bkg",
145  "Beam $\\nu_{e}$CC","NC", "$\\nu_{\\mu}$CC",
146  "$\\nu_{\\tau}$CC"};
147 
148  std::vector<TString> tags = {"nueSignal","TotBeamBkg",
149  "BeamnueCC","NC","numuCC","nutauCC"};
150 
151  // Initialize histos if first time
152  for ( int s = 0; s < 6; ++s ){
153  if ( systsum_plu_xp_prop[0] == 0 ){
154 
155  systsum_plu_xp_prop[s] = new TH1D( "Envelope_xpp_plu_"+tags[s], "Envelope_xpp_plu"+labels[s],
156  Nbins, 0, (Nbins*binW) );
157  systsum_plu_xp_combo[s] = new TH1D( "Envelope_xpc_plu_"+tags[s], "Envelope_xpc_plu"+labels[s],
158  Nbins, 0, (Nbins*binW) );
159  systsum_plu_nxp[s] = new TH1D( "Envelope_nxp_plu_"+tags[s], "Envelope_nxp_plu"+labels[s],
160  Nbins, 0, (Nbins*binW) );
161  systsum_min_xp_prop[s] = new TH1D( "Envelope_xpp_min_"+tags[s], "Envelope_xpp_min"+labels[s],
162  Nbins, 0, (Nbins*binW) );
163  systsum_min_xp_combo[s] = new TH1D( "Envelope_xpc_min_"+tags[s], "Envelope_xpc_min"+labels[s],
164  Nbins, 0, (Nbins*binW) );
165  systsum_min_nxp[s] = new TH1D( "Envelope_nxp_min_"+tags[s], "Envelope_nxp_min"+labels[s],
166  Nbins, 0, (Nbins*binW) );
167  }
168 
169  TH1D* absPlu = mcplu[s]->Clone();
170  TH1D* absMin = mcmin[s]->Clone();
171 
172  // Make each bin content in + larger than -
173  FormatErrorBand(absPlu, absMin, false, true);
174 
175  // Add % diff in quadrature for + and -
176  for ( int b = 1; b <= Nbins; ++b ){
177 
178  double binNom = mcnom[s]->GetBinContent(b);
179  double binPlu = absPlu ->GetBinContent(b);
180  double binMin = absMin ->GetBinContent(b);
181 
182  double shiftPlu = (binPlu - binNom) / binNom;
183  double shiftMin = (binMin - binNom) / binNom;
184 
185  // +_ shift bands on one side of the nominal are covered
186  // by the largest one. Truncate to +diff >= 0 -diff<= 0
187  shiftPlu = std::max( shiftPlu, 0 );
188  shiftMin = std::min( shiftMin, 0 );
189 
190  // sum bin in quadrature
191  if ( ltxcommand.Contains("nxp") ){
192 
193  shiftPlu = systsum_plu_nxp[s]->GetBinContent(b) + (shiftPlu*shiftPlu);
194  shiftMin = systsum_min_nxp[s]->GetBinContent(b) + (shiftMin*shiftMin);
195  systsum_plu_nxp[s]->SetBinContent(b, shiftPlu );
196  systsum_min_nxp[s]->SetBinContent(b, shiftMin );
197  }
198  if ( ltxcommand.Contains("xp_prop") ){
199 
200  shiftPlu = systsum_plu_nxp[s]->GetBinContent(b) + (shiftPlu*shiftPlu);
201  shiftMin = systsum_min_nxp[s]->GetBinContent(b) + (shiftMin*shiftMin);
202  systsum_plu_nxp[s]->SetBinContent(b, shiftPlu );
203  systsum_min_nxp[s]->SetBinContent(b, shiftMin );
204  }
205  if ( ltxcommand.Contains("xp_combo") ){
206 
207  shiftPlu = systsum_plu_nxp[s]->GetBinContent(b) + (shiftPlu*shiftPlu);
208  shiftMin = systsum_min_nxp[s]->GetBinContent(b) + (shiftMin*shiftMin);
209  systsum_plu_nxp[s]->SetBinContent(b, shiftPlu );
210  systsum_min_nxp[s]->SetBinContent(b, shiftMin );
211  }
212 
213  }// bins
214  }// components
215 
216 }
217 */
218 
219 //////////////////////////////////////////////////////////////////////
220  void SetVecValues_Diff( SystResults &MyVec, int Ind, double Sig0, double Sig1, double Bkg0, double Bkg1, double Tot0, double Tot1 ) {
221  MyVec.diff_sig[0][Ind] = Sig0;
222  MyVec.diff_sig[1][Ind] = Sig1;
223  MyVec.diff_bkg[0][Ind] = Bkg0;
224  MyVec.diff_bkg[1][Ind] = Bkg1;
225  MyVec.diff_tot[0][Ind] = Tot0;
226  MyVec.diff_tot[1][Ind] = Tot1;
227 
228  if (isNumuAna && Ind == 3) {
229  /*
230  std::cout << "\n\n\tLooking at a numu analysis in the 4th quantile...Now I can fill out the total bin!!\n\n" << std::endl;
231  std::cout << "\tFor sig0 I have; "
232  << MyVec.diff_sig[0][0] << " && " << MyVec.diff_sig[0][1] << " && "
233  << MyVec.diff_sig[0][2] << " && " << MyVec.diff_sig[0][3] << "\n"
234  << "\tFor Sig1 I have; "
235  << MyVec.diff_sig[1][0] << " && " << MyVec.diff_sig[1][1] << " && "
236  << MyVec.diff_sig[1][2] << " && " << MyVec.diff_sig[1][3] << "\n"
237  << std::endl;
238  */
239  MyVec.diff_sig[0][4] = (MyVec.diff_sig[0][0] + MyVec.diff_sig[0][1] + MyVec.diff_sig[0][2] + MyVec.diff_sig[0][3]) / 4;
240  MyVec.diff_sig[1][4] = (MyVec.diff_sig[1][0] + MyVec.diff_sig[1][1] + MyVec.diff_sig[1][2] + MyVec.diff_sig[1][3]) / 4;
241  MyVec.diff_bkg[0][4] = (MyVec.diff_bkg[0][0] + MyVec.diff_bkg[0][1] + MyVec.diff_bkg[0][2] + MyVec.diff_bkg[0][3]) / 4;
242  MyVec.diff_bkg[1][4] = (MyVec.diff_bkg[1][0] + MyVec.diff_bkg[1][1] + MyVec.diff_bkg[1][2] + MyVec.diff_bkg[1][3]) / 4;
243  MyVec.diff_tot[0][4] = (MyVec.diff_tot[0][0] + MyVec.diff_tot[0][1] + MyVec.diff_tot[0][2] + MyVec.diff_tot[0][3]) / 4;
244  MyVec.diff_tot[1][4] = (MyVec.diff_tot[1][0] + MyVec.diff_tot[1][1] + MyVec.diff_tot[1][2] + MyVec.diff_tot[1][3]) / 4;
245  //std::cout << "After that I have;\t sig0 " << MyVec.diff_sig[0][4] << " && sig1 " << MyVec.diff_sig[1][4] << std::endl;
246  }
247  }
248 
249  void SetVecValues_ChiSq( SystResults &MyVec, double Sig0, double Sig1, double Bkg0, double Bkg1, double Tot0, double Tot1 ) {
250  /*
251  std::cout << "\nFor ChiSq I had ;\n\t"
252  << MyVec.chisq_sig[0] << " && " << MyVec.chisq_sig[1] << " && "
253  << MyVec.chisq_bkg[0] << " && " << MyVec.chisq_bkg[1] << " && "
254  << MyVec.chisq_tot[0] << " && " << MyVec.chisq_tot[1] << std::endl;
255  */
256  MyVec.chisq_sig[0] += Sig0;
257  MyVec.chisq_sig[1] += Sig1;
258  MyVec.chisq_bkg[0] += Bkg0;
259  MyVec.chisq_bkg[1] += Bkg1;
260  MyVec.chisq_tot[0] += Tot0;
261  MyVec.chisq_tot[1] += Tot1;
262  /*
263  std::cout << "Now I have ;\n\t"
264  << MyVec.chisq_sig[0] << " && " << MyVec.chisq_sig[1] << " && "
265  << MyVec.chisq_bkg[0] << " && " << MyVec.chisq_bkg[1] << " && "
266  << MyVec.chisq_tot[0] << " && " << MyVec.chisq_tot[1] << std::endl;
267  */
268  }
269 //////////////////////////////////////////////////////////////////////
270 
271 void AddToSummaryV(TString ltxcommand, TString systName,
272  const std::vector<TH1*> & mcnom,
273  const std::vector<TH1*> & mcplu,
274  const std::vector<TH1*> & mcmin)
275 {
276  if( !(isNumuAna) && !(ltxcommand.Contains("Nue2020Axis"))) return;
277  //if( !(isNumuAna) && ltxcommand.Contains("3bin") ) return;
278  //if (systName.Contains("RelativeCalib") ) return;
279 
281  int VecIndex_xp = -1;
282  int VecIndex_nxp = -1;
283  int VecIndex_com = -1;
284  TString FillSyst = systName;
285  // If Numu I may want to use something one which is already in my vector...
286  if (isNumuAna)
287  {
288  if (systName.Contains("NumuEnergy")) {std::cout << "\n\n I don't want NumuEnergy in the table.\n\n" << std::endl; return; }// For now just don't include the full spectrum in the table...
289  int QLoc = systName.Index("Ehad");
290  /*std::max(
291  systName.Index("Ehad"), std::max(
292  systName.Index("Q2"), std::max(
293  systName.Index("Q3"), systName.Index("Q4") ) ) );*/
294  if (QLoc > 0) FillSyst = systName(0,QLoc);
295  if (ltxcommand.Contains("nxp") )
296  {
297  for (size_t resLoop=0; resLoop<systresults_nxp.size(); ++resLoop)
298  {
299  if (systresults_nxp[resLoop].syst == FillSyst)
300  {
301  VecIndex_nxp = resLoop;
302  break;
303  }
304  }
305  }
306  else if (ltxcommand.Contains("_xp_prop") ||
307  ltxcommand.Contains("_xp_numu") )
308  {
309  for (size_t resLoop=0; resLoop<systresults_xp_prop.size(); ++resLoop)
310  {
311  if (systresults_xp_prop[resLoop].syst == FillSyst)
312  {
313  VecIndex_xp = resLoop;
314  break;
315  }
316  }
317  }
318  else if (ltxcommand.Contains("_xp_combo") )
319  {
320  for (size_t resLoop=0; resLoop<systresults_xp_combo.size(); ++resLoop)
321  {
322  if (systresults_xp_combo[resLoop].syst == FillSyst)
323  {
324  VecIndex_com = resLoop;
325  break;
326  }
327  }
328  } // Which Vector Index do I fill?
329  }
330  ret.syst = FillSyst;
331 
332  // Number of bins from prediction histogram
333  int Nbins=mcnom[0]->GetNbinsX();
334  // TODO: can we define this somewhere in analysis header files?
335  // Would be better not to have it hard coded
336  int TotBins = 4; // for nue, 2 CVN + 1 peripheral + total
337  //int binNbins = Nbins/TotBins;
338  int binLow[3] = {1, 10, 19};
339  int binHigh[3] = {9, 18, mcnom[0]->GetNbinsX()};
340 
341  if (isNumuAna)
342  { // If numu
343  TotBins = 1;
344  //binNbins = Nbins;
345  }
346 
347  bool hasmin = mcmin.size()>0;
348 
349  // ### Fill the main bins.
350  for(int n=0; n < TotBins; ++n)
351  {
352 
353  int ArIndex = n; // for nue.
354  if (isNumuAna && systName.Contains("Q1") ) ArIndex = 0;
355  if (isNumuAna && systName.Contains("Q2") ) ArIndex = 1;
356  if (isNumuAna && systName.Contains("Q3") ) ArIndex = 2;
357  if (isNumuAna && systName.Contains("Q4") ) ArIndex = 3;
358 
359  //int bin1 = 1 + (n*binNbins);
360  //int bin2 = (n+1) * binNbins;
361  int bin1 = binLow[n];
362  int bin2 = binHigh[n];
363  if(n==TotBins-1) {bin1 = 1; bin2 = Nbins;}
364 
365  double Sig0 = -111;
366  double Sig1 = -111;
367 
368  if(mcnom[0] -> Integral(bin1,bin2) >0 )
369  {
370  if(mcplu[0])
371  Sig0 = mcplu[0]->Integral(bin1,bin2)/
372  mcnom[0]->Integral(bin1,bin2) - 1;
373  if(hasmin)
374  Sig1 = mcmin[0]->Integral(bin1,bin2)/
375  mcnom[0] -> Integral(bin1,bin2) - 1;
376  }
377 
378  double Bkg0 = -111;
379  double Bkg1 = -111;
380 
381  if(mcnom[1] -> Integral(bin1,bin2) >0 )
382  {
383  if(mcplu[1])
384  Bkg0 = mcplu[1]->Integral(bin1,bin2)/
385  mcnom[1] -> Integral(bin1,bin2) - 1;
386  if(hasmin)
387  Bkg1 = mcmin[1]->Integral(bin1,bin2)/
388  mcnom[1] -> Integral(bin1,bin2) - 1;
389  }
390 
391  double Tot0 = -111;
392  double Tot1 = -111;
393 
394  if( (mcnom[0] -> Integral(bin1,bin2) + mcnom[1]->Integral(bin1,bin2)) > 0 )
395  {
396  if(mcplu[0] && mcplu[1])
397  Tot0 =((mcplu[0]->Integral(bin1,bin2) + mcplu[1]->Integral(bin1,bin2))/
398  (mcnom[0]->Integral(bin1,bin2) + mcnom[1]->Integral(bin1,bin2)) -1 );
399  if(hasmin)
400  Tot1 =((mcmin[0]->Integral(bin1,bin2) + mcmin[1]->Integral(bin1,bin2))/
401  (mcnom[0]->Integral(bin1,bin2) + mcnom[1]->Integral(bin1,bin2)) -1 );
402  }
403 
404  // ### If numu how do I want to fill the bins?
405  if (VecIndex_xp == -1 && VecIndex_nxp == -1)
406  {
407  SetVecValues_Diff(ret, ArIndex, Sig0, Sig1, Bkg0, Bkg1, Tot0, Tot1 );
408  }
409  else if (ltxcommand.Contains("nxp"))
410  {
411  SetVecValues_Diff(systresults_nxp[VecIndex_nxp], ArIndex, Sig0,
412  Sig1, Bkg0, Bkg1, Tot0, Tot1 );
413  }
414  else if (ltxcommand.Contains("_xp_prop") ||
415  ltxcommand.Contains("_xp_numu"))
416  {
417  SetVecValues_Diff(systresults_xp_prop[VecIndex_xp], ArIndex, Sig0,
418  Sig1, Bkg0, Bkg1, Tot0, Tot1 );
419  }
420  else if (ltxcommand.Contains("_xp_combo"))
421  {
422  SetVecValues_Diff(systresults_xp_combo[VecIndex_com], ArIndex, Sig0,
423  Sig1, Bkg0, Bkg1, Tot0, Tot1 );
424  }
425  }//end bins
426 
427  // #### Chi-Square calculation section
428  double chisq_Sig0 = -111, chisq_Sig1 = -111;
429  double chisq_Bkg0 = -111, chisq_Bkg1 = -111;
430  double chisq_Tot0 = -111, chisq_Tot1 = -111;
431 
432  auto temp_nom = (TH1*)mcnom[0]->Clone(UniqueName().c_str());
433  temp_nom->Add (mcnom[1]);
434  auto temp_plu = (TH1*)mcplu[0]->Clone(UniqueName().c_str());
435  temp_plu->Add (mcplu[1]);
436 
437  TH1* temp_min;
438  if(hasmin)
439  {
440  temp_min = (TH1*)mcmin[0]->Clone(UniqueName().c_str());
441  temp_min->Add (mcmin[1]);
442  }
443 
444  chisq_Sig0 = CalcChi2(mcnom[0],mcplu[0]);
445  chisq_Bkg0 = CalcChi2(mcnom[1],mcplu[1]);
446  chisq_Tot0 = CalcChi2(temp_nom,temp_plu);
447  if(hasmin){
448  chisq_Sig1 = CalcChi2(mcnom[0],mcmin[0]);
449  chisq_Bkg1 = CalcChi2(mcnom[1],mcmin[1]);
450  chisq_Tot1 = CalcChi2(temp_nom,temp_min);
451  }
452 
453  if (VecIndex_xp == -1 && VecIndex_nxp == -1)
454  {
455  SetVecValues_ChiSq(ret, chisq_Sig0, chisq_Sig1, chisq_Bkg0,
456  chisq_Bkg1, chisq_Tot0, chisq_Tot1 );
457  }
458  else if (ltxcommand.Contains("nxp"))
459  {
460  SetVecValues_ChiSq(systresults_nxp[VecIndex_nxp], chisq_Sig0, chisq_Sig1,
461  chisq_Bkg0, chisq_Bkg1, chisq_Tot0, chisq_Tot1 );
462  }
463  else if (ltxcommand.Contains("_xp_prop") || ltxcommand.Contains("_xp_numu"))
464  {
466  systresults_xp_prop[VecIndex_xp], chisq_Sig0,
467  chisq_Sig1, chisq_Bkg0, chisq_Bkg1, chisq_Tot0, chisq_Tot1 );
468  }
469  else if (ltxcommand.Contains("_xp_combo"))
470  {
472  systresults_xp_combo[VecIndex_com], chisq_Sig0,
473  chisq_Sig1, chisq_Bkg0, chisq_Bkg1, chisq_Tot0, chisq_Tot1 );
474  }
475 
476  // ### Junk output
477  /*
478  SystResults OutRet = ret;
479  if (VecIndex_nxp!=-1) {
480  OutRet = systresults_nxp[VecIndex_nxp];
481  std::cout << "I am looking at a NO extrap result..." << std::endl;
482  } else if (VecIndex_xp!=-1) {
483  OutRet = systresults_xp_prop [VecIndex_xp ];
484  std::cout << "I am looking at an extrap result..." << std::endl;
485  } else {
486  std::cout << "This is the first time I've seen this one..." << std::endl;
487  }
488  for (uint zz=0; zz<4; ++zz) {
489  std::cout << "\nLooking at ret["<<zz<<"];"
490  << "\n\t It has sig values " << OutRet.diff_sig[0][zz] << " && " << OutRet.diff_sig[1][zz]
491  << "\n\t It has bkg values " << OutRet.diff_bkg[0][zz] << " && " << OutRet.diff_bkg[1][zz]
492  << "\n\t It has tot values " << OutRet.diff_tot[0][zz] << " && " << OutRet.diff_tot[1][zz]
493  << std::endl;
494  }
495  */
496  // ### Prepare to push back if neccessary.
497  delete temp_nom,temp_plu,temp_min;
498 
499  if(ltxcommand.Contains("ND")) return;
500  if(VecIndex_nxp==-1 && ltxcommand.Contains("nxp"))
501  {
502  systresults_nxp.push_back(ret);
503  return;
504  }
505  if(VecIndex_xp ==-1 && (ltxcommand.Contains("_xp_prop") ||
506  ltxcommand.Contains("_xp_numu") ) )
507  {
508  systresults_xp_prop.push_back(ret);
509  return;
510  }
511  if(VecIndex_com == -1 && ltxcommand.Contains("xp_combo"))
512  {
513  systresults_xp_combo.push_back(ret);
514  return;
515  }
516 }
517 
518 //////////////////////////////////////////////////////////////////////
519 void ComparisonTableOne(const std::vector<TH1*> & mcnom, const std::vector<TH1*> & mcshift,
520  const std::vector<TString> & labels, const TString ltxcommand =""){
521  unsigned int nentries = mcnom.size();
522  if(mcshift.size()!=nentries || labels.size()!=nentries){std::cerr << "%incompatible sizes\n"; return;}
523 
524  // Extract Sytematic name from title
525  TString title = mcnom[0]->GetTitle();
526  TString systName = FixPlotName(title(0, title.Index(" ")));
527 
528  PrintTableHeader(ltxcommand, "| l |"," | r | r | r || r r | ", 1);
529  tables << TString::Format("%20s & %10s & %10s & %10s & %10s & %10s \\\\ \n",
530  " ", "Nominal",
531  "Sh(+)",
532  "Sh(-)",
533  "\\%(+)",
534  "\\%(-)"
535  );
536  tables << "\\hline \n";
537  TString sum0, sum1;
538  for(unsigned int i=0;i<nentries;i++){
539 
540  TString base = (mcnom[i]->Integral()>100 ?
541  "%20s & %10.0f & %10.0f & %10s & %+10.0f & %10s \\\\ \n":
542  "%20s & %10.3f & %10.3f & %10s & %+10.1f & %10s \\\\ \n");
543  auto temp = TString::Format(base.Data(),
544  labels[i].Data(),
545  mcnom[i]->Integral(),
546  mcshift[i]->Integral(), "--",
547  ((mcnom[i]->Integral()>0) ? ((mcshift[i]->Integral()/mcnom[i]->Integral()) -1.) : 0) *100.,
548  "--");
549  tables << temp;
550  if(i == 0) sum0 = temp;
551  if(i == 1) sum1 = temp;
552  }
554  AddToSummaryV (ltxcommand, systName, mcnom, mcshift,{});
555 }
556 
557 
558 
559 //////////////////////////////////////////////////////////////////////
560 
561 void ComparisonTableOneNbins(const std::vector<TH1*> & mcnom, const std::vector<TH1*> & mcshift,
562  const std::vector<TString> & labels, TString ltxcommand, AxisType pidaxis){
563  //if(N<2) return;
564  int binLow[3] = {1, 10, 19};
565  int binHigh[3] = {9, 18, mcnom[0]->GetNbinsX()};
566  if(pidaxis == kNueNDbin) binLow[2] = 0;
567  tables << "\n \\vspace{30 pt} \n";
568 
569  //int Nbins=mcnom[0]->GetNbinsX();
570  //int binNbins=Nbins/N;
571 
572  if(ltxcommand=="") {std::cerr<< "WARNING ComparisonTableOneNbins ltxcommand\n"; return;}
573  output << "\\setlength{\\tabcolsep}{1pt} \\renewcommand{\\arraystretch}{1.2} \n";
574  output << "\n{\\small \\vspace{30pt} \n";
575  output << "\\" << MakeLatexCommandName(ltxcommand+"bins") << "\n\n";
576 
577  tables << "\\newcommand{\\" << MakeLatexCommandName(ltxcommand+"bins") << "}\n";
578  tables << "{\\centerline{\n";
579 
580 
581  for ( int n = 0; n<3; ++n){
582 
583  int bin1 = binLow[n];
584  int bin2 = binHigh[n];
585  if(bin1 == 0) continue;
586 
587  unsigned int nentries = mcnom.size();
588  if(mcshift.size()!=nentries || labels.size()!=nentries){std::cerr << "%incompatible sizes\n"; return;}
589  tables << "\\begin{tabular}{" << (n==0? "l":"") <<"| r r | r |} \n"
590  << std::setw(20) << (n==0?"&":"") << "\\multicolumn{3}{|c|}{Bin " << n+1 << "}\\\\ \\hline "
591  << (n==0?" & ":" ")
592  << std::setw(10) << "Nom." << " & "
593  << std::setw(10) << "Shift" << " & "
594  << std::setw(12) << "\\% Diff."
595  <<"\t \\\\ \\hline \n";
596  for(unsigned int i=0;i<nentries;i++){
597  tables.unsetf(std::ios_base::floatfield);
598 
599  if(mcnom[i] ->Integral(bin1,bin2) > 100) tables << std::fixed << std::setprecision(0);
600  else tables << std::setprecision(2);
601 
602  tables << std::setw(20); if( n==0 ) tables << labels[i] << " & ";
603  tables << std::setw(10) << mcnom[i] ->Integral(bin1,bin2) << " & "
604  << std::setw(10) << mcshift[i]->Integral(bin1,bin2) << " & "
605  << std::setw(10) << std::setprecision(1) << std::fixed << std::showpos
606  << ( mcnom[i]->Integral()>0 ? (((mcshift[i]->Integral(bin1,bin2)/mcnom[i]->Integral(bin1,bin2)) -1.)*100.):0.)
607  << std::noshowpos
608  <</* "\\%*/" \t \\\\ \n";
609  }
610  tables << "\\end{tabular} \n";
611 
612  }
613  tables << "} } \n";
614  output << "}\n";
615 
616 }
617 
618 void ComparisonTable(const std::vector<TH1*> & mcnom, const std::vector<TH1*> & mcplus, const std::vector<TH1*> & mcminus,
619  std::vector<TString> labels, TString ltxcommand = ""){
620  // Extract Sytematic name from title
621  output << " \n % " << std::string(60,'-') << std::endl;
622  TString title = mcnom[0]->GetTitle();
623  TString systName = FixPlotName(title(0, title.Index(" ")));
624 
625  unsigned int nentries = mcnom.size();
626  if(mcplus.size()!=nentries || mcminus.size()!=nentries || labels.size()!=nentries){std::cerr << "incompatible sizes\n"; return;}
627  if(ltxcommand == "") {std::cerr<< "WARNING Problem with comparisonTable ltx\n"; return;}
628 
629  output << "\\setlength{\\tabcolsep}{1pt} \\renewcommand{\\arraystretch}{1.2} \n";
630  output << "\\"<< MakeLatexCommandName(ltxcommand) << "\n\n";
631 
632  tables << "\\newcommand{\\" << MakeLatexCommandName(ltxcommand) << "}\n";
633  tables << "{\\centerline { \n \\begin{tabular}{ l | r r r | r r } \n"
634  << std::setw(25) << "Integral shift " << " & "
635  << std::setw(10) << "Nominal" << " & "
636  << std::setw(10) << "Shift (+)" << " & "
637  << std::setw(10) << "Shift (-)" << " & "
638  << std::setw(12) << "\\% Diff. (+)" << " & "
639  << std::setw(12) << "\\% Diff. (-)"
640  <<"\t \\\\ \\hline \n";
641  for(unsigned int i=0;i<nentries;i++){
642  tables.unsetf(std::ios_base::floatfield);
643  if(mcnom[i] ->Integral() > 100) tables << std::fixed << std::setprecision(0);
644  else tables << std::setprecision(2);
645  tables << std::setw(25) << labels[i] << " & "
646  << std::setw(10) << mcnom[i] ->Integral() << " & "
647  << std::setw(10) << mcplus[i] ->Integral() << " & "
648  << std::setw(10) << mcminus[i]->Integral() << " & "
649  << std::setw(10) << std::setprecision(1) << std::fixed << std::showpos
650  << ( mcnom[i]->Integral()>0 ? (((mcplus [i]->Integral()/mcnom[i]->Integral()) -1.)*100.):0.)
651  << /*"\\%"*/" & "
652  << std::setw(10)
653  << ( mcnom[i]->Integral()>0 ? (((mcminus[i]->Integral()/mcnom[i]->Integral()) -1.)*100.):0.)
654  << std::noshowpos
655  << /*"\\%"*/" \t \\\\ \n";
656  }
657  tables << "\\end{tabular} } \n } \n";
658 
659  // Can I add a table about the mean here???
660  // NOOO
661  // TODO fix latex command
662  if(isNumuAna){
663 
664  output << "\\setlength{\\tabcolsep}{1pt} \\renewcommand{\\arraystretch}{1.2} \n";
665  output << "\\vspace{25pt} \\" << MakeLatexCommandName(ltxcommand+"mean") << "\n\n";
666 
667  tables << "\\newcommand{\\" << MakeLatexCommandName(ltxcommand + "mean") << "}\n";
668  tables << "{\\centerline { \n \\begin{tabular}{ l | r r r | r r } \n"
669  << std::setw(25) << "Mean Value shift " << " & "
670  << std::setw(10) << "Nominal" << " & "
671  << std::setw(10) << "Shift (+)" << " & "
672  << std::setw(10) << "Shift (-)" << " & "
673  << std::setw(12) << "\\% Diff. (+)" << " & "
674  << std::setw(12) << "\\% Diff. (-)"
675  <<"\t \\\\ \\hline \n";
676  for(unsigned int i=0;i<nentries;i++){
677  tables.unsetf(std::ios_base::floatfield);
678  tables << std::setprecision(3);
679  tables << std::setw(25) << labels[i] << " & "
680  << std::setw(10) << mcnom[i] ->GetMean() << " & "
681  << std::setw(10) << mcplus[i] ->GetMean() << " & "
682  << std::setw(10) << mcminus[i]->GetMean() << " & "
683  << std::setw(10) << std::setprecision(1) << std::fixed << std::showpos
684  << ( mcnom[i]->Integral()>0 ? (((mcplus [i]->GetMean()/mcnom[i]->GetMean()) -1.)*100.):0.)
685  << /*"\\%"*/" & "
686  << std::setw(10)
687  << ( mcnom[i]->Integral()>0 ? (((mcminus[i]->GetMean()/mcnom[i]->GetMean()) -1.)*100.):0.)
688  << std::noshowpos
689  << /*"\\%"*/" \t \\\\ \n";
690  }
691  tables << "\\end{tabular} } \n } \n";
692 
693  // Can I add a table about the RMS here???
694  output << "\\setlength{\\tabcolsep}{1pt} \\renewcommand{\\arraystretch}{1.2} \n";
695  output << "\\vspace{25pt} \\" << MakeLatexCommandName(ltxcommand+"rms") << "\n\n";
696 
697  tables << "\\newcommand{\\" << MakeLatexCommandName(ltxcommand + "rms") << "}\n";
698  tables << "{\\centerline { \n \\begin{tabular}{ l | r r r | r r } \n"
699  << std::setw(25) << "RMS Value shift " << " & "
700  << std::setw(10) << "Nominal" << " & "
701  << std::setw(10) << "Shift (+)" << " & "
702  << std::setw(10) << "Shift (-)" << " & "
703  << std::setw(12) << "\\% Diff. (+)" << " & "
704  << std::setw(12) << "\\% Diff. (-)"
705  <<"\t \\\\ \\hline \n";
706  for(unsigned int i=0;i<nentries;i++){
707  tables.unsetf(std::ios_base::floatfield);
708  tables << std::setprecision(3);
709  tables << std::setw(25) << labels[i] << " & "
710  << std::setw(10) << mcnom[i] ->GetRMS() << " & "
711  << std::setw(10) << mcplus[i] ->GetRMS() << " & "
712  << std::setw(10) << mcminus[i]->GetRMS() << " & "
713  << std::setw(10) << std::setprecision(1) << std::fixed << std::showpos
714  << ( mcnom[i]->Integral()>0 ? (((mcplus [i]->GetRMS()/mcnom[i]->GetRMS()) -1.)*100.):0.)
715  << /*"\\%"*/" & "
716  << std::setw(10)
717  << ( mcnom[i]->Integral()>0 ? (((mcminus[i]->GetRMS()/mcnom[i]->GetRMS()) -1.)*100.):0.)
718  << std::noshowpos
719  << /*"\\%"*/" \t \\\\ \n";
720  }
721  tables << "\\end{tabular} } \n } \n";
722  }
723 
724  AddToSummaryV (ltxcommand, systName, mcnom, mcplus, mcminus);
725  //AddToEnvelop (ltxcommand, systName, mcnom, mcplus, mcminus);
726 
727 }
728 
729 TString FixMeSummary (double number, bool isdiff = true){
730  if(number <-100) return "--";
731  if(!isdiff) return TString::Format("%.3f",number);
732  TString pre="";
733  if(fabs(number*100) >= 7) pre = "\\cellcolor{\\myhighlightC}" ;
734  else if(fabs(number*100) >= 3) pre = "\\cellcolor{\\myhighlightB}" ;
735  else if(fabs(number*100) >= 1) pre = "\\cellcolor{\\myhighlightA}" ;
736  return pre + TString::Format("%+.2f",number*100);
737 
738 }
739 
740 void MakeSummaryV(TString sample = "tot", TString xptype = "nxp", bool sorted = false){
741 
742  auto results = systresults_nxp;
743  if(xptype.Contains ("prop")) results = systresults_xp_prop;
744  if(xptype.Contains ("combo")) results = systresults_xp_combo;
745  if(xptype.Contains ("pred_xp_numu")) results = systresults_xp_prop;
746  bool extrap = (xptype != "nxp");
747 
748  if(sorted && (sample == "sig"))std::sort(results.begin(), results.end(), sort_chisq_sig);
749  if(sorted && (sample == "bkg"))std::sort(results.begin(), results.end(), sort_chisq_bkg);
750  if(sorted && (sample == "tot"))std::sort(results.begin(), results.end(), sort_chisq_tot);
751 
752  TString summarylatex = MakeLatexCommandName("summaryT"+ xptype + sample + (sorted?"sorted":""));
753 
754  output << "\\"<< summarylatex <<"\n\n";
755 
756  summary << "\\newcommand{\\" << summarylatex << "}{\n";
757  summary << "\\hspace*{-1em} \n";
758  if (!isNumuAna)
759  summary << "\\begin{longtable}{p{0.3\\textwidth}||cc|cc|cc||cc||cc||}\n";
760  else
761  summary << "\\begin{longtable}{p{0.3\\textwidth}||cc|cc|cc||cc||cc||}\n";
762  summary << (extrap?"Extrap. ":"FDMC ") << sample << "\n";
763  TString BinQuant = isNumuAna? "Quant" : "Bin";
764  if (isNumuAna)
765  {
766  // --- Table for Numu
767  summary << " & \\multicolumn{2}{c|}{"<<BinQuant<<"1} & "
768  << " \\multicolumn{2}{c|}{"<<BinQuant<<"2} & "
769  << " \\multicolumn{2}{c|}{"<<BinQuant<<"3} & "
770  << " \\multicolumn{2}{c||}{"<<BinQuant<<"4} &"
771  << " \\multicolumn{2}{c||}{Chi2Test} \\\\ \n";
772  summary << "Syst & +(\\%) & -(\\%) & +(\\%) & -(\\%) & +(\\%) & -(\\%) & "
773  << " +(\\%) & - (\\%) & + & - \\\\ \n \\hline \n";
774  }
775  else
776  {
777  // --- Table for nue
778  summary << " & \\multicolumn{2}{c|}{"<<BinQuant<<"1} & "
779  << " \\multicolumn{2}{c|}{"<<BinQuant<<"2} &"
780  << " \\multicolumn{2}{c||}{"<<BinQuant<<"3} & "
781  << " \\multicolumn{2}{c||}{Tot} & "
782  << " \\multicolumn{2}{c||}{Chi2Test} \\\\ \n";
783  summary << "Syst & +(\\%) & -(\\%) & +(\\%) & -(\\%) & +(\\%) & -(\\%) & "
784  << " +(\\%) & -(\\%) & + & - \\\\ \n \\hline \n";
785  }
786  for (auto const & res:results){
787  summary.unsetf(std::ios_base::floatfield);
788  TString section = "\\hyperref[sec:"+MakeLatexCommandName(res.syst)+"]{"+res.syst+"} ";
789  if(sample == "sig")
790  {
791  summary << std::setprecision(3) << std::showpos
792  << section << " & "
793  << FixMeSummary(res.diff_sig[0][0]) << " & \t"
794  << FixMeSummary(res.diff_sig[1][0]) << " & \t"
795  << FixMeSummary(res.diff_sig[0][1]) << " & \t"
796  << FixMeSummary(res.diff_sig[1][1]) << " & \t"
797  << FixMeSummary(res.diff_sig[0][2]) << " & \t"
798  << FixMeSummary(res.diff_sig[1][2]) << " & \t";
799  if (!isNumuAna) // print Nue total column
800  summary << FixMeSummary(res.diff_sig[0][3]) << " & \t"
801  << FixMeSummary(res.diff_sig[1][3]) << " & \t";
802  else // print Numu quant 4
803  summary << FixMeSummary(res.diff_sig[0][3]) << " & \t"
804  << FixMeSummary(res.diff_sig[1][3]) << " & \t";
805  summary << std::noshowpos // for both Nue and Numu
806  << FixMeSummary(res.chisq_sig[0],false) << "& \t"
807  << FixMeSummary(res.chisq_sig[1],false) << " \t \\\\ \n";
808  }
809  else if(sample == "bkg")
810  {
811  summary << std::setprecision(3) << std::showpos
812  << section << " & "
813  << FixMeSummary(res.diff_bkg[0][0]) << " & \t"
814  << FixMeSummary(res.diff_bkg[1][0]) << " & \t"
815  << FixMeSummary(res.diff_bkg[0][1]) << " & \t"
816  << FixMeSummary(res.diff_bkg[1][1]) << " & \t"
817  << FixMeSummary(res.diff_bkg[0][2]) << " & \t"
818  << FixMeSummary(res.diff_bkg[1][2]) << " & \t";
819  if (!isNumuAna) // print Nue total
820  summary << FixMeSummary(res.diff_bkg[0][3]) << " & \t"
821  << FixMeSummary(res.diff_bkg[1][3]) << " & \t";
822  else // print Numu quant 4
823  summary << FixMeSummary(res.diff_bkg[0][3]) << " & \t"
824  << FixMeSummary(res.diff_bkg[1][3]) << " & \t";
825  summary << std::noshowpos
826  << FixMeSummary(res.chisq_bkg[0],false) << "& \t"
827  << FixMeSummary(res.chisq_bkg[1],false) << " \t \\\\ \n";
828  }
829  else
830  {
831  summary << std::setprecision(3) << std::showpos
832  << section << " & "
833  << FixMeSummary(res.diff_tot[0][0]) << " & \t"
834  << FixMeSummary(res.diff_tot[1][0]) << " & \t"
835  << FixMeSummary(res.diff_tot[0][1]) << " & \t"
836  << FixMeSummary(res.diff_tot[1][1]) << " & \t"
837  << FixMeSummary(res.diff_tot[0][2]) << " & \t"
838  << FixMeSummary(res.diff_tot[1][2]) << " & \t";
839  if (!isNumuAna) // print Nue total
840  summary << FixMeSummary(res.diff_tot[0][3]) << " & \t"
841  << FixMeSummary(res.diff_tot[1][3]) << " & \t";
842  else // print Numu quant 4
843  summary << FixMeSummary(res.diff_tot[0][3]) << " & \t"
844  << FixMeSummary(res.diff_tot[1][3]) << " & \t";
845  summary << std::noshowpos
846  << FixMeSummary(res.chisq_tot[0],false) << "& \t"
847  << FixMeSummary(res.chisq_tot[1],false) << " \t \\\\ \n";
848  }
849  }
850  summary << "\\end{longtable} \n } \n\n";
851 }
852 
853 void ComparisonTableNbins(const std::vector<TH1*> & mcnom, const std::vector<TH1*> & mcplus, const std::vector<TH1*> & mcminus,
854  std::vector<TString> labels,TString ltxcommand = "", AxisType pidaxis = kNue3bin){
855 
856  //if(N<3) return;
857  if(pidaxis != kNue3bin) return;
858  int N = 3;
859 
860  int binLow[3] = {1, 10, 19};
861  int binHigh[3] = {9, 18, mcnom[0]->GetNbinsX()};
862  if(pidaxis == kNueNDbin) binLow[2] = 0;
863  output << "\n % " << std::string(60,'-') << std::endl;
864 
865  //int Nbins=mcnom[0]->GetNbinsX();
866  //int binNbins=Nbins/N;
867 
868  unsigned int nentries = mcnom.size();
869  if(mcplus.size()!=nentries || mcminus.size()!=nentries || labels.size()!=nentries){std::cerr << "WARNING incompatible sizes\n"; return;}
870  if(ltxcommand == "") {std::cerr << "WARNING Problem with ComparisonTableNbins ltxcommand\n"; return;}
871 
872  output << "\n\n{\\footnotesize \\vspace{30pt} \n";
873  output << "\\setlength{\\tabcolsep}{1.5pt} \\renewcommand{\\arraystretch}{1.2} \n";
874  output << "\\" << MakeLatexCommandName(ltxcommand+"bins") << "\n\n";
875 
876  tables << "\\newcommand{\\" << MakeLatexCommandName(ltxcommand+"bins") << "}\n";
877  tables << "{\\centerline{\n";
878  tables << "\\begin{tabular}{ l ";
879  for (int n = 1; n <= N; ++n ) tables << "| r r r | r r ";
880  tables << "} \n & ";
881  for (int n = 1; n <= N; ++n ){tables << "\\multicolumn{5}{c|}{Bin " << n << "} " <<( n < N ? " & ":" \\\\ \\hline \n");}
882  tables << std::setw(20) << " & ";
883  for (int n = 1; n <= N; ++n ){
884  tables << std::setw(3) << "Nom."<< "&"
885  << std::setw(3) << "(+)" << "&"
886  << std::setw(3) << "(-)" << "&"
887  << std::setw(5) << "\\%(+)" << "&"
888  << std::setw(5) << "\\%(-)"
889  <<( n < N ? " & ":" \\\\ \\hline \n");
890  }
891  for(unsigned int i=0;i<nentries;i++){
892  tables << std::setw(20) << labels[i] << " & ";
893 
894  for ( int n = 0; n<N; ++n){
895 
896  //int bin1 = 1 + (n*binNbins);
897  //int bin2 = (n+1)*binNbins;
898  int bin1 = binLow[n];
899  int bin2 = binHigh[n];
900  if(bin1 == 0) continue;
901 
902  tables.unsetf(std::ios_base::floatfield);
903  if(mcnom[i] ->Integral(bin1,bin2) > 100) tables << std::fixed << std::setprecision(0);
904  else tables << std::setprecision(2);
905  tables << std::setw(5) << mcnom[i] ->Integral(bin1,bin2) << "&"
906  << std::setw(5) << mcplus[i] ->Integral(bin1,bin2) << "&"
907  << std::setw(5) << mcminus[i]->Integral(bin1,bin2) << "&"
908  << std::setw(5) << std::setprecision(1) << std::fixed << std::showpos
909  << ( mcnom[i]->Integral()>0 ? (((mcplus [i]->Integral(bin1,bin2)/mcnom[i]->Integral(bin1,bin2)) -1.)*100.):0.)
910  << "&"
911  << std::setw(5)
912  << ( mcnom[i]->Integral()>0 ? (((mcminus[i]->Integral(bin1,bin2)/mcnom[i]->Integral(bin1,bin2)) -1.)*100.):0.)
913  << std::noshowpos
914  << ( n + 1 < N ? " & ":" \\\\ \n");
915  }
916  }
917  tables << "\\end{tabular} \n";
918  tables << "}}\n";
919  output << "}\n"; //end footnotesize
920 
921 }
922 //////////////////////////////////////////////////////////////////////
923 void PrintRawEventCounts(TDirectory * dpred, TString title){
924 
925  std::vector<int> ndcounts;
926  std::vector<int> fdcounts;
927  std::vector<TString> labels;
928  //Data
929  labels.push_back("Data - $\\nu_\\mu $");
930  ndcounts.push_back(((TH1*)dpred->Get("extrap/MEextrap/Decomp/data_comp/hist"))->GetEntries());
931  fdcounts.push_back(0);
932  labels.push_back("Data - $\\nu_e $");
933  ndcounts.push_back(((TH1*)dpred->Get("extrap/EEextrap/Decomp/data_comp/hist"))->GetEntries());
934  fdcounts.push_back(0);
935 
936  //Signal
937  labels.push_back("Signal - $\\nu_\\mu \\to \\nu_e $");
938  labels.push_back("Signal - $\\bar{\\nu}_\\mu \\to \\bar{\\nu}_e $");
939  for(TString xp:{"ME","MEAnti"}){
940  ndcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoToTrueND/hist"))->GetEntries());
941  fdcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/TrueToRecoFD/hist"))->GetEntries());
942  }
943 
944  //Major bkg
945  labels.push_back("Bkg. - $\\nu_e \\to \\nu_e $");
946  labels.push_back("Bkg. - $\\nu_\\mu \\to \\nu_\\mu $");
947  labels.push_back("Bkg. - $NC \\to NC $");
948  for(TString xp:{"EE","MM","NC"}){
949  ndcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoND/hist"))->GetEntries());
950  fdcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/TrueToRecoFD/hist"))->GetEntries());
951  }
952 
953  //Minor bkg
954  labels.push_back("Bkg. - $\\nu_e \\to \\nu_\\tau $");
955  labels.push_back("Bkg. - $\\nu_\\mu \\to \\nu_\\tau $");
956  for(TString xp: {"ET","MT"}){
957  ndcounts.push_back(0);
958  fdcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoFD/hist"))->GetEntries());
959  }
960  labels.push_back("Bkg. - Other");
961  double fdcount_minor=0;
962  for(TString xp: {"EM","EEAnti","EMAnti","ETAnti","MMAnti","MTAnti"}){
963  fdcount_minor+=(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoFD/hist"))->GetEntries());
964  }
965  ndcounts.push_back(0);
966  fdcounts.push_back(fdcount_minor);
967 
968  //Print
969  output << "% " << std::string(60,'-') << std::endl;
970  output << "\n " << title << "\n";
971  output << "\\setlength{\\tabcolsep}{1.5pt} \\renewcommand{\\arraystretch}{1.2} \n";
972  output << "\\vspace{20pt} \n \\centering{\\begin{tabular}{ l | r | r } \n"
973  << std::setw(45) << " " << " & "
974  << std::setw(10) << "ND" << " & "
975  << std::setw(10) << "FD"
976  << " \\\\ \\hline \n";
977  for(unsigned int i=0;i<labels.size();i++){
978  output << std::setw(45) << labels[i] << " & "
979  << std::setw(10)<< ndcounts[i] << " & "
980  << std::setw(10)<< fdcounts[i] << " \\\\ ";
981  if(i==1 || i==3 || i==6) output << "\\hline \n";
982  else output << "\n";
983  }
984  output << "\\end{tabular} \n} \n";
985  output << "% " << std::string(60,'-') << std::endl;
986 }
987 
988 void PrintOscilationParameters (osc::IOscCalc * calc, bool usingdumb)
989 {
990  if(usingdumb){
991  output << "\n\\bigskip\n"
992  << "Using dumb oscillations - no solar term, mass hierarchy nor delta. \n"
993  << "$L= 810 km,$ "
994  << "$\\Delta m^{2} = 2.40e-3 eV^{-2},$ "
995  << "$\\sin^{2}2\\theta_{13} = 0.1,$ "
996  << "$\\sin^{2}\\theta_{23} = 0.5$ "
997  << ".\\\\\n\\medskip\n\n";
998  }
999  else{
1000  auto calc2 = dynamic_cast<osc::IOscCalcAdjustable *> (calc);
1001  output << "\n\\bigskip\n"
1002  << "Oscillation parameters:\n"
1003  << TString::Format("$L= %.0f km,$ ",calc2->GetL())
1004  << TString::Format("$\\rho= %.2f g/cm^3,$ ",calc2->GetRho())
1005  << TString::Format("$\\Delta m^{2} = %.2e eV^{-2},$ ",calc2->GetDmsq32())
1006  << TString::Format("$\\Delta m^{2} = %.2e eV^{-2},$ ",calc2->GetDmsq21())
1007  << TString::Format("$\\sin^{2}2\\theta_{12} = %.3f,$ ",util::sqr(sin(2*calc2->GetTh12())))
1008  << TString::Format("$\\sin^{2}2\\theta_{13} = %.3f,$ ",util::sqr(sin(2*calc2->GetTh13())))
1009  << TString::Format("$\\sin^{2}\\theta_{23} = %.3f,$ ",util::sqr(sin(calc2->GetTh23())))
1010  << TString::Format("$\\delta_{CP} = %.2f\\pi.$ ",calc2->GetdCP()/M_PI)
1011  << "\\\\\n\\medskip\n\n";
1012  }
1013 
1014 }
1015 
1016 }
TH1D * systsum_plu_xp_combo[6]
const XML_Char * name
Definition: expat.h:151
ofstream output
void PrintTableHeader(TString ltxcommand="", TString labelcol_style="l", TString bincol_style="| r r r", int nbins=1)
double diff_bkg[2][5]
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ComparisonTableOneNbins(const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcshift, const std::vector< TString > &labels, TString ltxcommand, int N=3)
TString FixPlotName(TString mystring)
double CalcChi2(TH1 *hmc, TH1 *hdata)
TH1D * systsum_min_xp_combo[6]
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double Integral(const Spectrum &s, const double pot, int cvnregion)
void SetVecValues_Diff(SystResults &MyVec, int Ind, double Sig0, double Sig1, double Bkg0, double Bkg1, double Tot0, double Tot1)
T sqrt(T number)
Definition: d0nt_math.hpp:156
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
bool sort_chisq_sig(SystResults a, SystResults b)
TH1D * systsum_plu_nxp[6]
OStream cerr
Definition: OStream.cxx:7
bool sort_chisq_bkg(SystResults a, SystResults b)
bool sort_chisq_tot(SystResults a, SystResults b)
void ComparisonTableOne(const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcshift, const std::vector< TString > &labels, const TString ltxcommand="")
osc::OscCalcDumb calc
double chi2()
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
cout<< t1-> GetEntries()
Definition: plottest35.C:29
double diff_tot[2][5]
std::vector< SystResults > systresults_nxp
#define M_PI
Definition: SbMath.h:34
TH1D * systsum_plu_xp_prop[6]
const int nbins
Definition: cellShifts.C:15
bool isNumuAna
std::vector< SystResults > systresults_xp_prop
TH1D * systsum_min_nxp[6]
void ComparisonTable(const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcplus, const std::vector< TH1 * > &mcminus, std::vector< TString > labels, TString ltxcommand="")
Long64_t nentries
const double a
void AddToSummaryV(TString ltxcommand, TString systName, const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcplu, const std::vector< TH1 * > &mcmin)
TH1D * systsum_min_xp_prop[6]
double sigma(TH1F *hist, double percentile)
TString FixMeSummary(double number, bool isdiff=true)
void SetVecValues_ChiSq(SystResults &MyVec, double Sig0, double Sig1, double Bkg0, double Bkg1, double Tot0, double Tot1)
void PrintLatexFigure(TString name, TString plotdir="plots/")
OStream cout
Definition: OStream.cxx:6
std::vector< SystResults > systresults_xp_combo
void MakeSummaryV(bool signal=true, TString xptype="nxp", bool sorted=false)
ofstream summary
TString MakeLatexCommandName(TString mystring)
ifstream in
Definition: comparison.C:7
T sin(T number)
Definition: d0nt_math.hpp:132
AxisType
void PrintTableFooter()
void PrintOscilationParameters(osc::IOscCalc *calc, bool usingdumb)
void ComparisonTableNbins(const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcplus, const std::vector< TH1 * > &mcminus, std::vector< TString > labels, TString ltxcommand="", int N=3)
ofstream tables
const hit & b
Definition: hits.cxx:21
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
double diff_sig[2][5]
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void PrintRawEventCounts(TDirectory *dpred, TString title)
enum BeamMode string