CorrectedSpectrum.cxx
Go to the documentation of this file.
2 
3 #include "TMath.h"
4 
6 
7 namespace fnex{
8 
9  //----------------------------------------------------------------------------
11  : fIsConfigured(false)
12  {
13  fDoOscillate[novadaq::cnv::kNEARDET] = pset.get< bool >("OscillateND", false);
14  fDoOscillate[novadaq::cnv::kFARDET] = pset.get< bool >("OscillateFD", true );
15 
16  fDisplayRangeMin = pset.get< double >("DisplayRangeMin", 0);
17  fDisplayRangeMax = pset.get< double >("DisplayRangeMax", 30);
18  fFitRangeMin = pset.get< double >("FitRangeMin", 0);
19  fFitRangeMax = pset.get< double >("FitRangeMax", 5);
20 
21  // Grab the FitEval to determine how to bin energy
22  fFitEvalType = pset.get< std::string >("FitEvaluationType", "ChiSq");
23 
24  for(size_t i = 0; i < fnex::kNumOscParams; ++i)
26 
27  }
28 
29  //----------------------------------------------------------------------------
31  {
33 
34  if(fFitEvalType == "UBL" && (fVarName == "Nu_RecoE" ||
35  fVarName == "Lep_RecoE" ||
36  fVarName == "Had_RecoE" ))
37  binning = fnex::VarBinning::Instance()->BinInfo(fVarName + "_UBL");
38  else if(md.IsNuMuQuantiles()){
39  LOG_DEBUG("CorrSpec_SimpleExtra")
40  << "Getting the binning for "
41  << fVarName
42  << " and numu quantiles";
43 
44  binning = fnex::VarBinning::Instance()->BinInfo(fVarName + "_Quantiles_NuMu");
45  }
46  //Call special single bin binning for NuE Peripheral cosmics
48  binning = fnex::VarBinning::Instance()->BinInfo(fVarName + "_NuE_Peripheral");
49 
50 
51  LOG_DEBUG("CorrectedSpectrum::SpectrumBinning")
52  << "Inside CorrectedSpectrum::SpectrumBinning -- "
53  << " fVarName = "
54  << fVarName
55  << " Metadata: "
56  << md.ToString();
57 
58  return binning;
59  }
60 
61  //----------------------------------------------------------------------------
63  return fIsConfigured;
64  }
65 
66  //----------------------------------------------------------------------------
69  }
70 
71  //----------------------------------------------------------------------------
72  // Range information
73  void CorrectedSpectrum::SetDisplayRange (double pMin, double pMax){
74  fDisplayRangeMin = pMin;
75  fDisplayRangeMax = pMax;
77  }
78 
79  //----------------------------------------------------------------------------
80  void CorrectedSpectrum::SetFitRange(double pMin, double pMax){
81  fFitRangeMin = pMin;
82  fFitRangeMax = pMax;
84  }
85 
86  //----------------------------------------------------------------------------
88  return true; //(fFitRangeMin > fDisplayRangeMin && fFitRangeMax < fDisplayRangeMax);
89  }
90 
91  //----------------------------------------------------------------------------
93  return true; //(fFitRangeMax < fDisplayRangeMax && fFitRangeMin > fDisplayRangeMin);
94  }
95 
96  // Spill fUncorrectedHistMap contents
97  //----------------------------------------------------------------------------
99  HistMap & pHistMap,
100  DetHistMap & pHist_All_Data,
101  DetHistMap & pHist_All_MC,
102  DetHistMap & pHist_All_Data_Over_MC,
103  std::string pSuffix){
104 
105  LOG_DEBUG("CorrectedSpectrum")
106  << "Drawing stacks."
107  << "Inside CorrectedSpectrum::DrawStacks(....) ----- Drawing Stacks";
108 
109  // \TODO Remember to convert into general "input HistMap" format
110  // when this is more put together
111 
112  // Unique id to prevent ROOT from sabotaging us
114 
117  int detEpochKey;
119  std::string namePrefix;
121  std::string chisqtitle;
122 
123  // need to make a new histogram for each interaction type
124  // The interaction type is usually a good indication of Data/MC because
125  // all data has the kUnknownInteraction type while MC is specified
126  std::map<novadaq::cnv::DetId, std::map<fnex::OscType_t, TH1D*> > intHists;
127 
128  // we also want the total data and MC spectrum for each epoch
129  std::map<int, TH1D*> detEpochHists;
130 
131  // get the binning for this corrected spectrum
132  fnex::Binning const* binning = nullptr;
133 
134  for( auto & tHistMap : pHistMap){
135  md = tHistMap.first;
136  det = md.detector;
137  namePrefix = md.DetectorString();
138  detEpochKey = (md.EpochKey() + md.DetectorKey()) * md.MCKey();
139  binning = this->SpectrumBinning(md);
140  LOG_DEBUG("CorrectedSpectrum::DrawStacks")
141  << "Inside CorrectedSpectrum::DrawStacks(....) "
142  << "----- Number of Bins from SpectrumBinning = "
143  << binning->nBins();
144 
145  if(md.isMC){
146 
147  namePrefix += "MC";
148 
149  //std::cout << " Looking at oscType == " << md.OscType() << " with string " << fnex::cOscType_Strings[md.OscType()] << std::endl;
150 
151  if(intHists[det].count(md.OscType()) < 1){
152  name = namePrefix + fnex::cOscType_Strings[md.OscType()] + pSuffix;
153 
154  LOG_DEBUG("CorrectedSpectrum")
155  << "Making histogram for "
156  << name
157  << " Metadata: "
158  << md.ToString();
159 
160  title = ";";
161  title += tHistMap.second.GetXaxis()->GetTitle();
162  title += ";Events";
163  intHists[det].insert(std::make_pair(md.OscType(),
164  tfd->make<TH1D>(name.c_str(),
165  title.c_str(),
166  binning->nBins(),
167  &(binning->BinLowEdges()[0]))
168  ));
169  }
170  LOG_DEBUG("CorrectedSpectrum::DrawStacks")
171  << "Inside CorrectedSpectrum::DrawStacks(....) Adding Histos --- Metadata: = "
172  << md.ToString()
173  << " LHS bins = "
174  << intHists[det][md.OscType()]->GetNbinsX()
175  << " RHS bins = "
176  << (&(tHistMap.second))->GetNbinsX();
177 
178  intHists[det][md.OscType()]->Add( &(tHistMap.second) );
179  intHists[det][md.OscType()]->SetOption("HIST");
180  intHists[det][md.OscType()]->SetFillStyle(cOscType_FillStyles[md.OscType()]);
181  intHists[det][md.OscType()]->SetFillColor(cOscType_Colors[md.OscType()]);
182 
183  } // end if MC
184  else{
185  namePrefix += "Data";
186  } // end if data
187 
188  if(detEpochHists.count(detEpochKey) < 1){
189  name = namePrefix + "Epoch" + md.EpochString() + pSuffix;
190  title = ";";
191  title += tHistMap.second.GetXaxis()->GetTitle();
192  title += ";Events";
193  detEpochHists.insert(std::make_pair(detEpochKey,
194  tfd->make<TH1D>(name.c_str(),
195  title.c_str(),
196  binning->nBins(),
197  &(binning->BinLowEdges()[0]))
198  ));
199  } // end if we need to make the epoch histogram
200  detEpochHists[detEpochKey]->Add( &(tHistMap.second) );
201  LOG_DEBUG("CorrectedSpectrum::DrawStacks")
202  << "Inside CorrectedSpectrum::DrawStacks(....) Adding Histos for epoch histogram --- Metadata: = "
203  << md.ToString()
204  << " LHS bins = "
205  << detEpochHists[detEpochKey]->GetNbinsX()
206  << " RHS bins = "
207  << ( &(tHistMap.second) )->GetNbinsX();
208 
209  } // end loop over input histograms
210 
212 
213  for(int i=0; i<2; ++i){
214  det = ids[i];
215 
216  // Make a canvas for each detector
217  TCanvas *can = nullptr;
218  THStack *stack = nullptr;
219  TH1D *data = nullptr;
220  TH1D *mc = nullptr;
221  TH1D *binchisq = nullptr;
222  TH1D *binchisqsum = nullptr;
223  double chisq = 0.;
224  double bin_data = 0.;
225  double bin_mc = 0.;
227 
228  // loop for the data for this detector
229  if(pHist_All_Data.count(det) < 1)
230  throw cet::exception("CorrectedSpectrum")
231  << "no data for detector " << det;
232 
233  name = (det == novadaq::cnv::kFARDET) ? "FD" : "ND";
234  title = ";"; title += pHist_All_Data.find(det)->second.GetXaxis()->GetTitle(); title += ";Events;";
235  chisqtitle = ";"; title += pHist_All_Data.find(det)->second.GetXaxis()->GetTitle(); title += ";Bin #chi^2";
236 
237  data = tfd->make<TH1D>((name + "Data" + pSuffix).c_str(),
238  title.c_str(),
239  binning->nBins(),
240  &(binning->BinLowEdges()[0]));
241  mc = tfd->make<TH1D>((name + "MC" + pSuffix).c_str(),
242  title.c_str(),
243  binning->nBins(),
244  &(binning->BinLowEdges()[0]));
245 
246  binchisq = tfd->make<TH1D>((name + "ChiSq" + pSuffix).c_str(),
247  chisqtitle.c_str(),
248  binning->nBins(),
249  &(binning->BinLowEdges()[0]));
250 
251  binchisqsum = tfd->make<TH1D>((name + "ChiSqSum" + pSuffix).c_str(),
252  (chisqtitle + " Sum").c_str(),
253  binning->nBins(),
254  &(binning->BinLowEdges()[0]));
255 
256  LOG_DEBUG("CorrectedSpectrum::DrawStacks")
257  << "Inside CorrectedSpectrum::DrawStacks(....) Adding Data Histos --- LHS bins = "
258  << data->GetNbinsX()
259  << " RHS bins = "
260  << ( &(pHist_All_Data.find(det)->second) )->GetNbinsX();
261 
262  data->Add( &(pHist_All_Data.find(det)->second) );
263 
264  LOG_DEBUG("CorrectedSpectrum::DrawStacks")
265  << "Inside CorrectedSpectrum::DrawStacks(....) Adding MC Histos --- LHS bins = "
266  << mc ->GetNbinsX()
267  << " RHS bins = "
268  << ( &(pHist_All_MC .find(det)->second) )->GetNbinsX();
269 
270  mc->Add( &(pHist_All_MC.find(det)->second) );
271 
272  if(data->GetNbinsX() < 1)
273  throw cet::exception("CorrectedSpectrum")
274  << "hData must have > 0 bins!";
275 
276  if(data->GetNbinsX() != mc->GetNbinsX())
277  throw cet::exception("CorrectedSpectrum")
278  << "hData and hMC must have same number of bins!";
279 
280 
281  // Find draw ranges that incorporate all bins used in the fit
282  int fitRangeBinMin = pHist_All_Data[novadaq::cnv::kNEARDET].GetXaxis()->FindBin(fFitRangeMin);
283  int fitRangeBinMax = pHist_All_Data[novadaq::cnv::kNEARDET].GetXaxis()->FindBin(fFitRangeMax)-1;
284 
285  int dispRangeBinMin = data->GetXaxis()->FindBin(fDisplayRangeMin);
286  int dispRangeBinMax = data->GetXaxis()->FindBin(fDisplayRangeMax)-1;
287 
288  double rangeUserMin = data->GetBinLowEdge(dispRangeBinMin);
289  double rangeUserMax = (data->GetBinLowEdge(dispRangeBinMax) + data->GetBinWidth(dispRangeBinMax));
290 
291  // Find fit range that incorporate all bins used in the fit
292  double fitRangeUserMin = data->GetBinLowEdge(data->GetXaxis()->FindBin(fFitRangeMin));
293  double fitRangeUserMax = (data->GetBinLowEdge(data->GetXaxis()->FindBin(fFitRangeMax) - 1)
294  + data->GetBinWidth(data->GetXaxis()->FindBin(fFitRangeMax) - 1));
295 
296  // Enforce these draw ranges on all histograms used in the draw
297  mc->GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
298  data->GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
299  if(fFitEvalType != "UBL") {
300  double minBinWidthData = 99999.;
301  for(int b = 0; b < data->GetXaxis()->GetNbins(); ++b){
302  if( minBinWidthData > data->GetBinWidth(b) )
303  minBinWidthData = data->GetBinWidth(b);
304  }
305  for(int b = 0; b < data->GetXaxis()->GetNbins(); ++b){
306  LOG_DEBUG("CorrectedSpectrum")
307  << "Normalising data bin "
308  << b
309  << " with width "
310  << data->GetBinWidth(b)
311  << " and contents "
312  << data->GetBinContent(b)
313  << " -> "
314  << data->GetBinContent(b) / data->GetBinWidth(b);
315 
316  data->SetBinContent(b, data->GetBinContent(b) * (minBinWidthData /
317  data->GetBinWidth(b)));
318  }
319  }
320 
321  for( auto & tHistMap : pHistMap){
322  md = tHistMap.first;
323  if(!md.isMC) continue;
324  intHists[md.detector][md.OscType()]->GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
325  }
326 
327 
328  // Make the stack
329  stack = tfd->make<THStack>((name + "Stacks" + pSuffix).c_str(), title.c_str());
330  title = fVarName + " (" + pSuffix + ") (Stacks)" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD");
331  name = fVarName + "_" + pSuffix + "_Stacks_" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD") + id->getUniqID();
332  LOG_DEBUG("CorrectedSpectrum") << "Making stack " << name;
333 
334  // Make the canvas
335  can = tfd->make<TCanvas>(name.c_str(), title.c_str(), 1600, 800);
336  can->Divide(2,1);
337  can->cd(1);
338 
339  // Set up the upper pad
340  name = "Upper";
341  TPad padUpper(name.c_str(), name.c_str(), 0, 0.3, 1, 1.0);
342  padUpper.SetBottomMargin(0.016);
343  //padLower.SetGridx(); // Jenny hates this, it seems. Boo Jenny.
344  padUpper.Draw();
345  padUpper.cd();
346 
347 
348  // Make an index to organize the histograms by increasing integral
349  std::vector<unsigned int> hist_index (intHists[det].size(), 0);
350  std::vector<fnex::osc_type> hist_osctype(intHists[det].size(), kOscTypeNoneSet);
351  std::vector<double> hist_size (intHists[det].size(), 0);
352  unsigned int curIndex = 0;
353  for( auto & intHist : intHists[det]){
354  hist_index [curIndex] = curIndex;
355  hist_osctype[curIndex] = intHist.first;
356  hist_size [curIndex] = intHist.second->Integral();
357  curIndex++;
358  }
359  sort(hist_index.begin(), hist_index.end(),
360  [&](const int& a, const int& b) {
361  return (hist_size[a] < hist_size[b]);
362  }
363  );
364 
365  // What are the signal osc types? We want to plot them last
366  // no matter what
367  std::vector<fnex::osc_type> signal_osctype;
368 
369  if(det == novadaq::cnv::kFARDET){
370  if(md.IsNuMuSelected()){
371  signal_osctype.push_back(kMCmm);
372  signal_osctype.push_back(kMCmmbar);
373  }
374  else if(md.IsNuESelected()){
375  signal_osctype.push_back(kMCme);
376  signal_osctype.push_back(kMCmebar);
377  }
378  }
379 
380  // Add histogram to the stack in order of increasing integral, ignoring the signal
381  // for now -- we'll add it last
382  LOG_DEBUG("CorrectedSpectrum")
383  << "Adding "
384  << hist_index.size()
385  << " to stack";
386 
387 
388  TH1D* hist = nullptr;
389  for(auto ihist : hist_index){
390  hist = intHists[det][hist_osctype[ihist]];
391 
392  LOG_DEBUG("CorrectedSpectrum::DrawStacks")
393  << "Inside CorrectedSpectrum::DrawStacks(....) Adding Histos to stacks --- Bins = "
394  << hist->GetNbinsX()
395  << " "
396  << hist->Integral()
397  << " "
398  << *(hist->GetName());
399 
400  max += hist->GetMaximum();
401 
402  if(fFitEvalType != "UBL") {
403  double minBinWidthMC = 99999.;
404  for(int b = 0; b < hist->GetXaxis()->GetNbins(); ++b){
405  if( minBinWidthMC > hist->GetBinWidth(b) )
406  minBinWidthMC = hist->GetBinWidth(b);
407  }
408  for(int b = 0; b < hist->GetXaxis()->GetNbins(); ++b){
409  LOG_DEBUG("CorrectedSpectrum")
410  << "Normalising hist bin "
411  << b
412  << " with width "
413  << hist->GetBinWidth(b)
414  << " and contents "
415  << hist->GetBinContent(b)
416  << " -> "
417  << hist->GetBinContent(b) / hist->GetBinWidth(b);
418 
419  hist->SetBinContent(b, hist->GetBinContent(b) * ( minBinWidthMC /
420  hist->GetBinWidth(b)));
421  }
422  }
423 
424  if(std::find(signal_osctype.begin(),
425  signal_osctype.end(),
426  hist_osctype[ihist]) != signal_osctype.end()) continue;
427 
428  stack->Add( hist );
429 
430  LOG_DEBUG("CorrectedSpectrum")
431  << "Added hist "
432  << *(hist->GetName());
433  }
434 
435  // Add the signal last
436  for(auto isig : signal_osctype){
437  LOG_DEBUG("CorrectedSpectrum")
438  << "Adding signal "
439  << *(intHists[det][isig]->GetName());
440 
441  stack->Add( intHists[det][isig] );
442  }
443 
444  // Make legend and add intHists entries to it
445  TLegend leg(0.55, 0.3, 0.9, 0.9);
446  leg.SetBorderSize(0);
447  leg.SetFillStyle(0);
448  leg.AddEntry(data, "Data", "lep");
449  for( auto & intHist : intHists[det]){
450  leg.AddEntry(intHist.second, fnex::cOscType_LatexStrings[intHist.first].c_str(), "lfp" );
451  }
452 
453  // In UBL fit we also want the spline drawn
454  TSpline3 spline;
455  if(fFitEvalType == "UBL" && det == novadaq::cnv::kFARDET) {
456  spline = this->DrawSpline(mc, fitRangeUserMin, fitRangeUserMax);
457  TH1F *dummyHist = new TH1F("dummyHist", "dummyHist", 10, 0, 10);
458  dummyHist->SetLineWidth(3);
459  dummyHist->SetLineColor(kOrange+7);
460  spline.SetLineWidth(3);
461  spline.SetLineColor(kOrange+7);
462  leg.AddEntry(dummyHist, "Spline Fit", "l");
463  }
464 
465 
466  // Draw stack so that we can start manipulating it
467  stack->Draw("HIST");
468  stack->GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
469  stack->Draw("HIST");
470 
471  // Get/set the max range
472  max = stack->GetMaximum();
473  if(max < data->GetMaximum()) max = data->GetMaximum();
474  data ->SetMaximum(max * 1.25);
475  stack->SetMaximum(max * 1.25);
476 
477  // Remove labels from stack (want these on ratio plot)
478  stack->GetXaxis()->SetLabelSize(0);
479 
480  // Make lines to show fit range
481  double lineYMin = 0;
482  double lineYMax = max;
483  TLine fitMinLine(fitRangeUserMin, lineYMin, fitRangeUserMin, lineYMax);
484  TLine fitMaxLine(fitRangeUserMax, lineYMin, fitRangeUserMax, lineYMax);
485  fitMinLine.SetLineColor(kRed);
486  //fitMinLine.SetLineStyle(6);
487  fitMinLine.SetLineWidth(2);
488  fitMaxLine.SetLineColor(kRed);
489  //fitMaxLine.SetLineStyle(6);
490  fitMaxLine.SetLineWidth(2);
491  leg.AddEntry(&fitMinLine, "FitRegion", "l");
492 
493  // Draw everything for final presentation
494  stack->Draw("HIST");
495  data ->SetLineWidth(2);
496  data ->SetLineColor(kBlack);
497  data ->SetMarkerColor(kBlack);
498  data ->SetMarkerSize(1);
499  data ->SetMarkerStyle(20);
500  data ->Draw("E0 P0 SAME");
501  if(this->DrawLowerFitLine()) fitMinLine.Draw("SAME");
502  if(this->DrawUpperFitLine()) fitMaxLine.Draw("SAME");
503  leg.Draw();
504 
505  // If we're running the UBL fit we want to add the spline
506  if(fFitEvalType == "UBL" && det == novadaq::cnv::kFARDET)
507  spline.Draw("C same");
508 
509  // Update the canvas to reflect these changes
510  can->Update();
511 
512  // Let's go back to the parent canvas so that we can move to the lower pad
513  can->cd(1);
514 
515  // Set up the lower pad
516  name = "Lower";
517  TPad padLower(name.c_str(), name.c_str(), 0, 0.05, 1, 0.3);
518  padLower.SetTopMargin(0);
519  padLower.SetBottomMargin(0.25);
520  //padLower.SetGridx(); // Jenny hates this, it seems. Boo Jenny.
521  padLower.Draw();
522  padLower.cd();
523 
524  // Draw the ratio
525  TH1D pDMCRat(pHist_All_Data_Over_MC[det]);
526  pDMCRat.GetYaxis()->SetTitle("Data/MC");
527  if(fFitEvalType == "UBL" && det == novadaq::cnv::kFARDET) {
528  pDMCRat.Reset();
529  pDMCRat.GetYaxis()->SetTitle("Data/Spline");
530  unsigned int nbins = data->GetNbinsX();
531  for(unsigned int ibin = 1; ibin < nbins + 1; ++ibin) {
532  double dataVal = data->GetBinContent(ibin);
533  double splineVal = spline.Eval(data->GetBinCenter(ibin));
534  double dataErr = data->GetBinError(ibin);
535  double splineErr = mc->GetBinError(ibin);
536 
537  double ratDS;
538  if(splineVal == 0) ratDS = 0;
539  else ratDS = dataVal / splineVal;
540 
541  double ratErrDS;
542  if(splineVal == 0 || dataVal == 0) ratErrDS = 0;
543  else ratErrDS = ratDS * std::sqrt(std::pow(dataErr / dataVal, 2) +
544  std::pow(splineErr / splineVal, 2));
545  pDMCRat.SetBinContent(ibin, ratDS);
546  pDMCRat.SetBinError (ibin, ratErrDS);
547  }
548  }
549  pDMCRat.SetLineColor(kBlack);
550  pDMCRat.SetLineWidth(2);
551  pDMCRat.SetMarkerColor(kBlack);
552  pDMCRat.SetTitle("");
553 
554  pDMCRat.GetYaxis()->CenterTitle();
555  pDMCRat.GetYaxis()->SetNdivisions(505);
556  pDMCRat.GetYaxis()->SetTitleSize(20);
557  pDMCRat.GetYaxis()->SetTitleFont(43);
558  pDMCRat.GetYaxis()->SetTitleOffset(1.55);
559  pDMCRat.GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
560  pDMCRat.GetYaxis()->SetLabelSize(15);
561 
562  pDMCRat.GetXaxis()->SetTitle(data->GetXaxis()->GetTitle());
563  pDMCRat.GetXaxis()->SetTitleSize(20);
564  pDMCRat.GetXaxis()->SetTitleFont(43);
565  pDMCRat.GetXaxis()->SetTitleOffset(4.);
566  pDMCRat.GetXaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
567  pDMCRat.GetXaxis()->SetLabelSize(15);
568 
569  pDMCRat.SetStats(false);
570  pDMCRat.GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
571  pDMCRat.Draw("E1");
572 
573  // Set up the Equal1 line
574  TLine yEquals1;
575  yEquals1.SetLineWidth(1);
576  yEquals1.SetLineStyle(2);
577  yEquals1.SetLineColor(kBlack);
578  yEquals1.DrawLine(rangeUserMin, 1.0,
579  rangeUserMax, 1.0 );
580 
581 
582  // Set up the fit range lines
583  // KM June 10 2016
584  double smallLineYMin = pDMCRat.GetMinimum();
585  double smallLineYMax = pDMCRat.GetMaximum();
586  TLine fitMinLineSmall(fitRangeUserMin, smallLineYMin, fitRangeUserMin, smallLineYMax);
587  TLine fitMaxLineSmall(fitRangeUserMax, smallLineYMin, fitRangeUserMax, smallLineYMax);
588  fitMinLineSmall.SetLineColor(kRed);
589  //fitMinLineSmall.SetLineStyle(6);
590  fitMinLineSmall.SetLineWidth(2);
591  fitMaxLineSmall.SetLineColor(kRed);
592  //fitMaxLineSmall.SetLineStyle(6);
593  fitMaxLineSmall.SetLineWidth(2);
594  if(this->DrawLowerFitLine()) fitMinLineSmall.Draw("SAME");
595  if(this->DrawUpperFitLine()) fitMaxLineSmall.Draw("SAME");
596 
597  // Update the canvas with the ratio pad
598  can->Update();
599 
600  // Make and draw point legend
601  can -> cd(2);
602 
603  // Get fixed parameters and systematics
604  auto const& fixedParameters = fRefInputPoint.FixedParameters(det);
605  //auto const& fixedSystematics = fRefInputPoint.FixedSystematics(det);
606 
607  // Get chisq and number of degrees of freedom (ndof)
608  int nFreeOscPars = fRefInputPoint.ParameterSpaceLocation(det).size() - fixedParameters.size();
609  //int nFreeSystPars = fRefInputPoint.SystematicPulls(det).size() - fixedSystematics.size();
610  int nFreePars = nFreeOscPars;// + nFreeSystPars;//Free systematic should always be paired with nuisance parameters,
611  //i.e. they don't contribute to chi sq
612 
613  //std::cout << "FIXEDPARSIZE: " << nFreeOscPars << std::endl;
614  //std::cout << "FIXEDSYSTSIZE: " << nFreeSystPars << std::endl;
615  int nNonZeroBins = 0;
616  chisq = 0.0;
617  double thischisq = 0.0;
618  for(int ibin = fitRangeBinMin; ibin <= fitRangeBinMax; ++ibin){
619  bin_data = data->GetBinContent(ibin);
620  bin_mc = pHist_All_MC[det].GetBinContent(ibin);
621  if(bin_mc <= 0) continue;
622  ++nNonZeroBins;
623  thischisq = ( bin_mc - bin_data );
624  if(bin_data != 0) thischisq += bin_data * log( bin_data / bin_mc );
625  binchisq->SetBinContent(ibin, 2.0 * thischisq);
626  chisq += 2.0 * thischisq;
627  binchisqsum->SetBinContent(ibin, chisq);
628  //std::cout << " Bin: " << ibin << " ThisChisq: " << 2.0*thischisq << " chisq: " << chisq << std::endl;
629  }
630  // LOG_VERBATIM("CorrectedSpectrum")
631  // << "NNONZEROBINS: " << nNonZeroBins
632  // << " NFREEPARS: " << nFreePars;
633  int ndof = nNonZeroBins - nFreePars;
634  //std::cout << "NDOF: " << ndof << std::endl;
635 
636  // Make parameter legend
637  bool displayChisq = true;
638  if(fFitEvalType == "UBL" && det == novadaq::cnv::kFARDET) displayChisq = false;
639  auto oscLeg = this->MakeOscillationLegend(chisq, ndof, det, displayChisq);
640 
641  // Make a pad for this one
642  TPad oscPad("oscPad", "oscPad", 0., 0.5, 1., 1.0);
643  oscPad.SetBottomMargin(0);
644  oscPad.Draw();
645  oscPad.cd();
646  oscLeg.Draw();
647 
648  can->cd(2);
649 
650  // Describe systematics space in the same way
651  auto sysLeg = this->MakeSystematicsLegend(chisq, ndof, det, false);
652 
653  // Make a pad for this one
654  TPad sysPad("sysPad", "sysPad", 0., 0.0, 1., 0.5);
655  sysPad.SetTopMargin(0);
656  sysPad.Draw();
657  sysPad.cd();
658  sysLeg.Draw();
659 
660  // We're done here. Draw the canvas describing par / syst space
661  // and update the canvas
662  can->Update();
663 
664  // Save the canvas to file
665  can->Write();
666 
667  }
668 
669  // We're done here
670  return;
671 
672  }
673 
674  // Spill fUncorrectedHistMap contents
675  //----------------------------------------------------------------------------
677  DetHistMap & pHist_All_Data,
678  DetHistMap & pHist_All_MC,
679  DetHistMap & pHist_All_Data_Over_MC,
680  std::string pSuffix)
681  {
682  return;
683  }
684 
685  // Spill fCorrectedHistMap contents
686  //----------------------------------------------------------------------------
688 
689  // Call DrawStacks on Uncorrected and Corrected results,
690  // each in their own folder
691  art::TFileDirectory uncorrDir = tfd->mkdir( "UncorrectedStacks" );
692  art::TFileDirectory corrDir = tfd->mkdir( "CorrectedStacks" );
693 
694  LOG_DEBUG("CorrectedSpectrum::DrawStacks")
695  << "Inside: CorrectedSpectrum::DrawStacks(tfd) "
696  << "--- Calling DrawStacks for Uncorrected Stack";
697 
698  DrawStacks(&uncorrDir,
703  "Uncorrected");
704  DrawStacks(&corrDir,
709  "Corrected" );
710 
711  return;
712  }
713 
714  // Spill fCorrectedHist_All contents
715  //----------------------------------------------------------------------------
717  DetHistMap & pHist_All_Data,
718  DetHistMap & pHist_All_MC,
719  DetHistMap & pHist_All_Data_Over_MC,
720  std::string pSuffix ){
721 
722 
723  LOG_DEBUG("DrawDataVsMC")
724  << " Running drawdatavsmc for"
725  << pSuffix;
726 
727  // Unique id to prevent ROOT from sabotaging us
729 
733 
734  // Make a canvas for each detector
735  TCanvas *can = nullptr;
736  TPad *pad1 = nullptr;
737  TPad *pad2 = nullptr;
738  TH1D *data = nullptr;
739  TH1D *mc = nullptr;
740  TH1D *rat = nullptr;
741  TH1D *ghost = nullptr;
742 
744 
746  for(int i = 0; i < 2; ++i){
747 
748  det = ids[i];
749 
750  // loo for the data for this detector
751  if(pHist_All_Data.count(det) < 1){
752  throw cet::exception("CorrectedSpectrum")
753  << "no data for detector " << det;
754  }
755  name = (det == novadaq::cnv::kFARDET) ? "FD" : "ND";
756  title = ";"; title += pHist_All_Data.find(det)->second.GetXaxis()->GetTitle(); title += ";Events";
757 
758  data = tfd->make<TH1D>((name + "Data" + pSuffix).c_str(),
759  title.c_str(),
760  pHist_All_Data.find(det)->second.GetXaxis()->GetNbins(),
761  pHist_All_Data.find(det)->second.GetXaxis()->GetXmin(),
762  pHist_All_Data.find(det)->second.GetXaxis()->GetXmax());
763  mc = tfd->make<TH1D>((name + "MC" + pSuffix).c_str(),
764  title.c_str(),
765  pHist_All_MC.find(det)->second.GetXaxis()->GetNbins(),
766  pHist_All_MC.find(det)->second.GetXaxis()->GetXmin(),
767  pHist_All_MC.find(det)->second.GetXaxis()->GetXmax());
768 
769 
770  data->Add( &(pHist_All_Data.find(det)->second) );
771  mc ->Add( &(pHist_All_MC.find(det)->second) );
772 
773  // Find draw ranges that incorporate all bins used in the fit
774  double rangeUserMin = data->GetBinLowEdge(data->GetXaxis()->FindBin(fDisplayRangeMin));
775  double rangeUserMax = (data->GetBinLowEdge(data->GetXaxis()->FindBin(fDisplayRangeMax)-1)
776  + data->GetBinWidth(data->GetXaxis()->FindBin(fDisplayRangeMax)-1));
777 
778  // Find fit range that incorporate all bins used in the fit
779  double fitRangeUserMin = data->GetBinLowEdge(data->GetXaxis()->FindBin(fFitRangeMin));
780  double fitRangeUserMax = (data->GetBinLowEdge(data->GetXaxis()->FindBin(fFitRangeMax)-1)
781  + data->GetBinWidth(data->GetXaxis()->FindBin(fFitRangeMax)-1));
782 
783  // Enforce these draw ranges on all histograms used in the draw
784  data -> GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
785  mc -> GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
786 
787 
788  if(data->GetNbinsX() < 1)
789  throw cet::exception("CorrectedSpectrum")
790  << "hData must have > 0 bins!";
791 
792  if(data->GetNbinsX() != mc->GetNbinsX())
793  throw cet::exception("CorrectedSpectrum")
794  << "hData and hMC must have same number of bins!";
795 
796  ghost=(TH1D *)data->Clone("ghost");
797  ghost->Reset(0);
798  ghost->SetLineColor(10);
799  ghost->SetMarkerColor(10);
800 
801 
802  title = fVarName + " (" + pSuffix + ") (Data Vs MC)" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD");
803  name = fVarName + "_" + pSuffix + "_DataVsMC_" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD") + id->getUniqID();
804 
805 
806  can = tfd->make<TCanvas>(name.c_str(), title.c_str(), 700, 800);
807  pad1 = new TPad(name.c_str(),name.c_str(),0.0,0.37,1.0,1.0,0);
808  pad2 = new TPad(name.c_str(),name.c_str(),0.0,0.0,1.0,0.37,0.);
809 
810  pad1->SetBottomMargin(0.005);
811  pad1->Draw();
812  pad2->SetTopMargin(0.0);
813  pad2->SetBottomMargin(0.25);
814  pad2->Draw();
815 
816  pad1->cd();
817 
818  if(i==1){
819  mc->Scale(0.001);
820  data->Scale(0.001);
821  mc->GetYaxis()->SetTitle("Events #times 10^{3}");
822  }
823 
824  // Make and draw legend
825  TLegend leg(0.6, 0.6, 0.85, 0.85);
826  if(i==0) leg.AddEntry(ghost, "NO#nuA FD", "p");
827  if(i==1) leg.AddEntry(ghost, "NO#nuA ND", "p");
828  // leg.SetBorderWidth(0);
829  leg.AddEntry(data, "Data", "pe");
830  leg.AddEntry(mc , "MC", "l");
831 
832  max = data->GetMaximum();
833  if(max > mc->GetMaximum()) mc->SetMaximum(max);
834 
835  mc->GetYaxis()->SetRangeUser(-0.02*max,1.2*max);
836 
837  mc->GetXaxis()->SetTitle("Reconstructed Energy (GeV)");
838  mc->GetYaxis()->CenterTitle();
839  mc->GetYaxis()->SetTitleOffset(1.45);
840  mc->GetYaxis()->SetTitleFont(43); // Absolute font size in pixel (precision 3)
841  mc->GetYaxis()->SetTitleSize(24);
842  mc->GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
843  mc->GetYaxis()->SetLabelSize(20);
844 
845 
846  mc->Draw("HIST");
847  data->Draw("E1 SAME");
848 
849  mc->SetStats(kFALSE);
850 
851  // Make lines to show fit range
852  double lineYMin = 0;
853  double lineYMax = max;
854  TLine fitMinLine(fitRangeUserMin, lineYMin, fitRangeUserMin, lineYMax);
855  TLine fitMaxLine(fitRangeUserMax, lineYMin, fitRangeUserMax, lineYMax);
856  fitMinLine.SetLineColor(kRed);
857  //fitMinLine.SetLineStyle(6);
858  fitMinLine.SetLineWidth(2);
859  fitMaxLine.SetLineColor(kRed);
860  //fitMaxLine.SetLineStyle(6);
861  fitMaxLine.SetLineWidth(2);
862  if(this->DrawLowerFitLine()) fitMinLine.Draw();
863  if(this->DrawUpperFitLine()) fitMaxLine.Draw();
864  leg.AddEntry(&fitMinLine, "FitRegion", "l");
865  leg.Draw();
866 
867  mc ->SetLineColor(kMagenta+1);
868  data->SetLineColor(kBlack);
869  mc ->SetMarkerColor(kMagenta+1);
870  data->SetMarkerColor(kBlack);
871  mc ->SetMarkerStyle(20);
872  data->SetMarkerStyle(20);
873  mc ->SetLineWidth(2);
874  data->SetLineWidth(4);
875  data->SetMarkerSize(4);
876 
877  pad2->cd();
878 
879  // Draw the ratio
880  TH1D pDMCRat(pHist_All_Data_Over_MC[det]);
881  pDMCRat.SetLineColor(kBlack);
882  pDMCRat.SetLineWidth(2);
883  pDMCRat.SetMarkerColor(kBlack);
884  pDMCRat.SetTitle("");
885 
886 
887  pDMCRat.GetYaxis()->SetTitle("Data/MC");
888  pDMCRat.GetYaxis()->CenterTitle();
889  pDMCRat.GetYaxis()->SetNdivisions(505);
890  pDMCRat.GetYaxis()->SetTitleSize(20);
891  pDMCRat.GetYaxis()->SetTitleFont(43);
892  pDMCRat.GetYaxis()->SetTitleOffset(1.55);
893  pDMCRat.GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
894  pDMCRat.GetYaxis()->SetLabelSize(15);
895 
896  pDMCRat.GetXaxis()->SetTitle(data->GetXaxis()->GetTitle());
897  pDMCRat.GetXaxis()->SetTitleSize(20);
898  pDMCRat.GetXaxis()->SetTitleFont(43);
899  pDMCRat.GetXaxis()->SetTitleOffset(4.);
900  pDMCRat.GetXaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
901  pDMCRat.GetXaxis()->SetLabelSize(15);
902 
903  pDMCRat.SetStats(false);
904  pDMCRat.GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
905  pDMCRat.GetYaxis()->SetDecimals();
906  pDMCRat.SetStats(kFALSE);
907  pDMCRat.Draw("E1");
908 
909  //crosscheck of ratio - so recalculating it to make sure the one returned in the code is actually right
910  rat=(TH1D *)data->Clone("rat");
911  rat->Reset();
912  rat->Divide(data,mc);
913 
914  rat->GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
915  rat->Draw("e1 same");
916  rat->SetLineColor(kMagenta+1);
917  rat->SetMarkerStyle(20);
918  rat->SetMarkerColor(kMagenta+1);
919  rat->SetMarkerSize(0.5);
920  rat->SetLineWidth(1);
921  rat->SetStats(kFALSE);
922 
923  // Set up the Equal1 line
924  TLine yEquals1;
925  yEquals1.SetLineWidth(1);
926  yEquals1.SetLineStyle(2);
927  yEquals1.SetLineColor(kBlack);
928  yEquals1.DrawLine(rangeUserMin, 1.0,
929  rangeUserMax, 1.0 );
930 
931  // Set up the fit range lines
932  double smallLineYMin = rat->GetMinimum();
933  double smallLineYMax = rat->GetMaximum();
934  TLine fitMinLineSmall(fitRangeUserMin, smallLineYMin, fitRangeUserMin, smallLineYMax);
935  TLine fitMaxLineSmall(fitRangeUserMax, smallLineYMin, fitRangeUserMax, smallLineYMax);
936  fitMinLineSmall.SetLineColor(kRed);
937  //fitMinLineSmall.SetLineStyle(6);
938  fitMinLineSmall.SetLineWidth(2);
939  fitMaxLineSmall.SetLineColor(kRed);
940  //fitMaxLineSmall.SetLineStyle(6);
941  fitMaxLineSmall.SetLineWidth(2);
942  if(this->DrawLowerFitLine()) fitMinLineSmall.Draw("SAME");
943  if(this->DrawUpperFitLine()) fitMaxLineSmall.Draw("SAME");
944 
945  can->Update();
946 
947  // Save the canvas to file
948  can->Write();
949 
950  }
951 
952  // We're done here
953  return;
954  }
955 
956  // Spill fUncorrectedHist_All contents
957  //----------------------------------------------------------------------------
959 
960  // Get ART file service
961  //art::ServiceHandle<art::TFileService> tfs;
962 
963  // Call DrawStacks on Uncorrected and Corrected results,
964  // each in their own folder
965  art::TFileDirectory uncorrDir = tfd->mkdir( "UncorrectedDataVsMC" );
966  art::TFileDirectory corrDir = tfd->mkdir( "CorrectedDataVsMC" );
967  DrawDataVsMC(&uncorrDir,
971  "Uncorrected");
972  DrawDataVsMC(&corrDir,
976  "Corrected");
977 
978  return;
979  }
980 
981  // Spill fUncorrectedHist_All contents
982  //----------------------------------------------------------------------------
984  {
985 
986  // Get ART file service
987  //art::ServiceHandle<art::TFileService> tfs;
988  DrawCorrVsUncorr(tfd,
992  "CorrVsUncorr");
993  return;
994  }
995 
996 
997 
998  // Record histogram counts
999  //----------------------------------------------------------------------------
1001  HistMap & pHistMap,
1002  DetHistMap & pHist_All_Data,
1003  DetHistMap & pHist_All_MC,
1004  DetHistMap & pHist_All_Data_Over_MC,
1005  std::string pSuffix){
1006 
1007 
1008  // \TODO Remember to convert into general "input HistMap" format
1009  // when this is more put together
1010 
1011  // Find fit range that incorporate all bins used in the fit
1012  int fitRangeBinMin = pHist_All_Data[novadaq::cnv::kNEARDET].GetXaxis()->FindBin(fFitRangeMin);
1013  int fitRangeBinMax = pHist_All_Data[novadaq::cnv::kNEARDET].GetXaxis()->FindBin(fFitRangeMax) - 1;
1014  double fitRangeUserMin = pHist_All_Data[novadaq::cnv::kNEARDET].GetBinLowEdge(fitRangeBinMin);
1015  double fitRangeUserMax = (pHist_All_Data[novadaq::cnv::kNEARDET].GetBinLowEdge(fitRangeBinMax) +
1016  pHist_All_Data[novadaq::cnv::kNEARDET].GetBinWidth(fitRangeBinMax) );
1017 
1018  int underFlowBinMin = 0;
1019  int underFlowBinMax = fitRangeBinMin - 1;
1020  if(underFlowBinMax < underFlowBinMin) underFlowBinMax = underFlowBinMin;
1021 
1022  int overFlowBinMax = pHist_All_Data[novadaq::cnv::kNEARDET].GetNbinsX();
1023  int overFlowBinMin = fitRangeBinMax + 1;
1024  if(overFlowBinMin > overFlowBinMax) underFlowBinMin = overFlowBinMax;
1025 
1026  // Unique id to prevent ROOT from sabotaging us
1028 
1031  std::string name;
1033 
1034  // need to make a new histogram for each interaction type
1035  // The interaction type is usually a good indication of Data/MC because
1036  // all data has the kUnknownInteraction type while MC is specified
1037  std::map<novadaq::cnv::DetId, std::map<fnex::OscType_t, TH1D*> > intHists;
1038 
1039  for( auto & tHistMap : pHistMap){
1040  md = tHistMap.first;
1041 
1042  if(md.isMC){
1043  det = md.detector;
1044 
1045  name = (det == novadaq::cnv::kNEARDET) ? "ND" : "FD";
1046 
1047  if(intHists[det].count(md.OscType()) < 1){
1048  name += fnex::cOscType_Strings[md.OscType()] + pSuffix;
1049  title = ";";
1050  title += tHistMap.second.GetXaxis()->GetTitle();
1051  title += ";Events";
1052  intHists[det].insert(std::make_pair(md.OscType(),
1053  tfd->make<TH1D>(name.c_str(),
1054  title.c_str(),
1055  1,
1056  fitRangeUserMin,
1057  fitRangeUserMax)
1058  ));
1059  }
1060 
1061  intHists[det][md.OscType()]->AddBinContent( 1, tHistMap.second.Integral(fitRangeBinMin, fitRangeBinMax ) );
1062  intHists[det][md.OscType()]->AddBinContent( 0, tHistMap.second.Integral(underFlowBinMin, underFlowBinMax));
1063  intHists[det][md.OscType()]->AddBinContent( 2, tHistMap.second.Integral(overFlowBinMin, overFlowBinMax ));
1064  intHists[det][md.OscType()]->SetOption("HIST");
1065  intHists[det][md.OscType()]->SetFillStyle(cOscType_FillStyles[md.OscType()]);
1066  intHists[det][md.OscType()]->SetFillColor(cOscType_Colors[md.OscType()]);
1067 
1068  } // end if MC
1069  } // end loop over input histograms
1070 
1071 
1072 
1073  for(int i=0; i<2; i++){
1074 
1075  if(i==0) det = novadaq::cnv::kNEARDET;
1076  else det = novadaq::cnv::kFARDET;
1077 
1078 
1079  // Make a canvas for each detector
1080  TCanvas *can = nullptr;
1081  TH1D *data = nullptr;
1082  TH1D *mc = nullptr;
1083  THStack *stack = nullptr;
1084 
1085  // loo for the data for this detector
1086  if(pHist_All_Data.count(det) < 1)
1087  throw cet::exception("CorrectedSpectrum")
1088  << "no data for detector " << det;
1089 
1090  name = (det == novadaq::cnv::kFARDET) ? "FD" : "ND";
1091  title = ";"; title += pHist_All_Data.find(det)->second.GetXaxis()->GetTitle(); title += ";Counts";
1092 
1093  data = tfd->make<TH1D>((name + "Data" + pSuffix).c_str(),
1094  title.c_str(),
1095  1,
1096  fitRangeUserMin,
1097  fitRangeUserMax);
1098  mc = tfd->make<TH1D>((name + "MC" + pSuffix).c_str(),
1099  title.c_str(),
1100  1,
1101  fitRangeUserMin,
1102  fitRangeUserMax);
1103 
1104  data->AddBinContent(1, pHist_All_Data.find(det)->second.Integral(fitRangeBinMin, fitRangeBinMax ));
1105  data->AddBinContent(0, pHist_All_Data.find(det)->second.Integral(underFlowBinMin, underFlowBinMax ));
1106  data->AddBinContent(2, pHist_All_Data.find(det)->second.Integral(overFlowBinMin, overFlowBinMax ));
1107 
1108  mc ->AddBinContent(1, pHist_All_MC.find(det)->second.Integral(fitRangeBinMin, fitRangeBinMax ));
1109  mc ->AddBinContent(0, pHist_All_MC.find(det)->second.Integral(underFlowBinMin, underFlowBinMax ));
1110  mc ->AddBinContent(2, pHist_All_MC.find(det)->second.Integral(overFlowBinMin, overFlowBinMax ));
1111 
1112  for( auto & tHistMap : pHistMap){
1113  md = tHistMap.first;
1114  if(!md.isMC) continue;
1115  intHists[md.detector][md.OscType()]->GetXaxis()->SetRangeUser(fitRangeUserMin, fitRangeUserMax);
1116  }
1117 
1118 
1119  // Make the stack
1120  stack = tfd->make<THStack>((name + "Counts" + pSuffix).c_str(), title.c_str());
1121  title = fVarName + " (" + pSuffix + ") (Counts)" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD");
1122  name = fVarName + "_" + pSuffix + "_Counts_" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD") + id->getUniqID();
1123 
1124  can = tfd->make<TCanvas>(name.c_str(), title.c_str(), 1600, 800);
1125  can->Divide(2,1);
1126  can->cd(1);
1127 
1128  // Make legend and fill stacks
1129 
1130  TH1D blank;
1131  blank.SetLineColor(kWhite);
1132  blank.SetMarkerColor(kWhite);
1133  blank.SetFillColor(kWhite);
1134 
1135 
1136  double max = 0.;
1137 
1138  /*
1139  // Make an index to organize the histograms by increasing integral
1140  std::vector<unsigned int> hist_index(intHists[det].size(), 0);
1141  std::vector<fnex::osc_type> hist_osctype(intHists[det].size(), kOscTypeNoneSet);
1142  std::vector<double> hist_size(intHists[det].size(), 0);
1143  unsigned int curIndex=0;
1144  for( auto & intHist : intHists[det]){
1145  hist_index[curIndex] = curIndex;
1146  hist_osctype[curIndex] = intHist.first;
1147  hist_size[curIndex] = intHist.second->Integral();
1148  curIndex++;
1149  }
1150  for (unsigned int i = 0 ; i != intHists[det].size() ; i++) hist_index[i] = i;
1151  sort(hist_index.begin(), hist_index.end(),
1152  [&](const int& a, const int& b) {
1153  return (hist_size[a] < hist_size[b]);
1154  }
1155  );
1156 
1157  // What are the signal osc types? We want to plot them last
1158  // no matter what
1159  std::vector<fnex::osc_type> signal_osctype;
1160  if(md.IsNuMuSelected()){
1161  signal_osctype.push_back(kMCmm);
1162  signal_osctype.push_back(kMCmmbar);
1163  }else if(md.IsNuESelected()){
1164  signal_osctype.push_back(kMCme);
1165  signal_osctype.push_back(kMCmebar);
1166  }
1167 
1168  // Add histogram to the stack in order of increasing integral, ignoring the signal
1169  // for now -- we'll add it last
1170  for(unsigned int ihist = 0; ihist < hist_index.size(); ihist++){
1171  intHists[det][hist_osctype[hist_index[ihist]]];
1172  max += intHists[det][hist_osctype[hist_index[ihist]]]->GetMaximum();
1173  if(std::find(signal_osctype.begin(), signal_osctype.end(), hist_osctype[hist_index[ihist]]) != signal_osctype.end()) continue;
1174  stack->Add( intHists[det][hist_osctype[hist_index[ihist]]] );
1175  }
1176 
1177  // Add the signal last
1178  for(unsigned int isig=0; isig<signal_osctype.size(); isig++) stack->Add( intHists[det][signal_osctype[isig]] );
1179  */
1180 
1181  TLegend leg(0.1, 0.1, 0.9, 0.9);
1182  leg.AddEntry(&blank, "Event Count ( + Outflow)", "p");
1183 
1184 
1185  std::stringstream ss_entry;
1186  ss_entry.str("");
1187  ss_entry << "Data" << " : "
1188  << std::fixed << std::setprecision(2) << data->GetBinContent(1)
1189  << " ( + " << std::fixed << std::setprecision(2)
1190  << (data->GetBinContent(0)+data->GetBinContent(2))
1191  << " )";
1192  leg.AddEntry(data, ss_entry.str().c_str(), "lep");
1193 
1194  ss_entry.str("");
1195  ss_entry << "MC TOT" << " : "
1196  << std::fixed << std::setprecision(2) << mc->GetBinContent(1)
1197  << " ( + " << std::fixed << std::setprecision(2)
1198  << (mc->GetBinContent(0)+mc->GetBinContent(2))
1199  << " )";
1200  leg.AddEntry(mc, ss_entry.str().c_str(), "");
1201 
1202 
1203  for( auto & intHist : intHists[det]){
1204  max += intHist.second->GetMaximum();
1205  stack->Add( intHist.second );
1206  ss_entry.str("");
1207  ss_entry << fnex::cOscType_LatexStrings[intHist.first].c_str()
1208  << " : " << std::fixed << std::setprecision(2) << intHist.second->GetBinContent(1)
1209  << " ( + " << std::fixed << std::setprecision(2)
1210  << (intHist.second->GetBinContent(0) + intHist.second->GetBinContent(2))
1211  << " )";
1212  leg.AddEntry(intHist.second, ss_entry.str().c_str(), "lfp");
1213  }
1214 
1215  // Draw stack so that we can start manipulating it
1216  stack->Draw("HIST");
1217  stack->GetXaxis()->SetRangeUser(fitRangeUserMin, fitRangeUserMax);
1218  stack->Draw("HIST");
1219 
1220 
1221  max = stack->GetMaximum();
1222  if(max < data->GetMaximum()) max = data->GetMaximum();
1223  //std::cout << " STACK MAX: " << stack->GetMaximum() << std::endl;
1224  //std::cout << " DATA MAX: " << data->GetMaximum() << std::endl;
1225  //std::cout << " TOT MAX: " << max << std::endl;
1226  data ->SetMaximum(max);
1227  stack->SetMaximum(max);
1228  stack->Draw("HIST");
1229  data ->Draw("E1 SAME");
1230 
1231  can->cd(2);
1232  leg.Draw();
1233  can->Update();
1234 
1235  // Save the canvas to file
1236  can->Write();
1237  }
1238 
1239  // We're done here
1240  return;
1241 
1242  }
1243 
1244  // Spill fCorrectedHistMap contents
1245  //----------------------------------------------------------------------------
1247 
1248  // Call DrawStacks on Uncorrected and Corrected results,
1249  // each in their own folder
1250  art::TFileDirectory uncorrDir = tfd->mkdir( "UncorrectedCounts" );
1251  art::TFileDirectory corrDir = tfd->mkdir( "CorrectedCounts" );
1252  DrawCounts(&uncorrDir,
1257  "Uncorrected");
1258  DrawCounts(&corrDir,
1263  "Corrected" );
1264 
1265  return;
1266  }
1267 
1268  //----------------------------------------------------------------------------
1269  TSpline3 CorrectedSpectrum::DrawSpline(TH1D *mc, double fitRangeUserMin, double fitRangeUserMax){
1270  std::cout << "running DrawSpline" << std::endl;
1271 
1272  // We want to generate a spline fit for the plot
1273  // Acquire knots from histograms
1274  const int firstSplineBin = mc->GetXaxis()->FindBin(fitRangeUserMin);
1275  const int lastSplineBin = mc->GetXaxis()->FindBin(fitRangeUserMax);
1276  const int nbins = lastSplineBin - firstSplineBin + 1;
1277 
1278  // Assign knots to vector
1279  double binCenter, binContent;
1280  std::vector<std::pair<double, double>> knotVector;
1281  for( int ibin = 0; ibin < nbins; ++ibin) {
1282  int bin = mc->GetBin(firstSplineBin + ibin);
1283 
1284  if(ibin == 0) binCenter = fitRangeUserMin;
1285  else if(ibin == nbins - 1) binCenter = fitRangeUserMax;
1286  else binCenter = mc->GetBinCenter(bin);
1287 
1288  binContent = mc->GetBinContent(bin);
1289  knotVector.push_back(std::make_pair(binCenter, binContent));
1290  }
1291 
1292  // Spline constructor requires array, so convert quickly
1293  const int nKnots = knotVector.size();
1294  double knot_x[nKnots], knot_y[nKnots];
1295  for(int iKnot = 0; iKnot < nKnots; ++iKnot) {
1296  knot_x[iKnot] = knotVector.at(iKnot).first;
1297  knot_y[iKnot] = knotVector.at(iKnot).second;
1298  }
1299 
1300  // Gradient for start and end points
1301  double gradStart = (knot_y[1] - knot_y[0]) / (knot_x[1] - knot_x[0]);
1302  double gradStop = (knot_y[nKnots-1] - knot_y[nKnots-2]) /
1303  (knot_x[nKnots-1] - knot_y[nKnots-2]);
1304 
1305  // Create spline
1306  TSpline3 spline("spline", knot_x, knot_y, nKnots, "b1e1", gradStart, gradStop);
1307  return spline;
1308 }
1309 
1310  // Spill fUncorrectedHistMap contents
1311  //----------------------------------------------------------------------------
1312  std::unique_ptr<TCanvas> CorrectedSpectrum::GetStacksCanvasCopy(HistMap & pHistMap,
1313  DetHistMap & pHist_All_Data,
1314  DetHistMap & pHist_All_MC,
1315  DetHistMap & pHist_All_Data_Over_MC,
1316  novadaq::cnv::DetId & pDet,
1317  std::string pSuffix){
1318 
1319  // Unique id to prevent ROOT from sabotaging us
1321  fnex::MetaData md = pHistMap.begin()->first;
1322  std::string name;
1324  std::string chisqtitle;
1325 
1326  // What are the signal osc types? We want to plot them last
1327  // no matter what
1328  std::vector<fnex::osc_type> signal_osctype;
1329  if(pDet == novadaq::cnv::kFARDET){
1330  if(md.IsNuMuSelected()){
1331  signal_osctype.push_back(kMCmm);
1332  signal_osctype.push_back(kMCmmbar);
1333  }
1334  else if(md.IsNuESelected()){
1335  signal_osctype.push_back(kMCme);
1336  signal_osctype.push_back(kMCmebar);
1337  }
1338  }
1339  // need to make a new histogram for each interaction type
1340  // The interaction type is usually a good indication of Data/MC because
1341  // all data has the kUnknownInteraction type while MC is specified
1342  std::map<fnex::OscType_t, TH1D*> intHists;
1343 
1344  fnex::Binning const* binning = nullptr;
1345  double minWidth = 1.e6;
1346 
1347  for( auto & tHistMap : pHistMap){
1348  md = tHistMap.first;
1349 
1350  if(md.detector != pDet) continue;
1351 
1352  binning = this->SpectrumBinning(md);
1353 
1354  auto const& lowEdges = binning->BinLowEdges();
1355  for(size_t b = 0; b < lowEdges.size() - 1; ++b)
1356  minWidth = std::min(lowEdges[b+1] - lowEdges[b], minWidth);
1357 
1358  if(md.isMC){
1359 
1360  name = (pDet == novadaq::cnv::kNEARDET) ? "ND" : "FD";
1361 
1362  if(intHists.count(md.OscType()) < 1){
1363  name += fnex::cOscType_Strings[md.OscType()] + pSuffix;
1364  title = ";";
1365  title += tHistMap.second.GetXaxis()->GetTitle();
1366  title += ";Events";
1367  intHists.insert(std::make_pair(md.OscType(),
1368  new TH1D(name.c_str(),
1369  title.c_str(),
1370  binning->nBins(),
1371  &(binning->BinLowEdges()[0]))
1372  ));
1373  }
1374 
1375 
1376  LOG_DEBUG("CorrectedSpectrum")
1377  << "Adding osc type "
1379  << " to interaction histograms";
1380 
1381  intHists[md.OscType()]->Add( &(tHistMap.second) );
1382  intHists[md.OscType()]->SetOption("HIST");
1383  intHists[md.OscType()]->SetFillStyle(cOscType_FillStyles[md.OscType()]);
1384  intHists[md.OscType()]->SetFillColor(cOscType_Colors[md.OscType()]);
1385 
1386  } // end if MC
1387  } // end loop over input histograms
1388 
1389  // Make a canvas
1390  std::unique_ptr<TCanvas> can;
1391  THStack *stack = nullptr;
1393 
1394  // loop for the data for this detector
1395  if(pHist_All_Data.count(pDet) < 1)
1396  throw cet::exception("CorrectedSpectrum")
1397  << "no data for detector " << pDet;
1398 
1399  name = (pDet == novadaq::cnv::kFARDET) ? "FD" : "ND";
1400  title = ";"; title += pHist_All_Data.find(pDet)->second.GetXaxis()->GetTitle(); title += ";Events";
1401  chisqtitle = ";"; title += pHist_All_Data.find(pDet)->second.GetXaxis()->GetTitle(); title += ";Bin #chi^2";
1402 
1403  TH1D* data = new TH1D((name + "Data" + pSuffix).c_str(),
1404  title.c_str(),
1405  binning->nBins(),
1406  &(binning->BinLowEdges()[0]));
1407  TH1D* mc = new TH1D((name + "MC" + pSuffix).c_str(),
1408  title.c_str(),
1409  binning->nBins(),
1410  &(binning->BinLowEdges()[0]));
1411 
1412  /*TH1D* binchisq = new TH1D((name + "ChiSq" + pSuffix).c_str(),
1413  chisqtitle.c_str(),
1414  binning->nBins(),
1415  &(binning->BinLowEdges()[0]));
1416  */
1417 
1418  data->Add( &(pHist_All_Data.find(pDet)->second) );
1419  mc ->Add( &(pHist_All_MC.find(pDet)->second) );
1420 
1421  if(data->GetNbinsX() < 1)
1422  throw cet::exception("CorrectedSpectrum")
1423  << "hData must have > 0 bins!";
1424 
1425  if(data->GetNbinsX() != mc->GetNbinsX())
1426  throw cet::exception("CorrectedSpectrum")
1427  << "hData and hMC must have same number of bins!";
1428 
1429 
1430  // Find draw ranges that incorporate all bins used in the fit
1431 // int fitRangeBinMin = pHist_All_Data[novadaq::cnv::kNEARDET].GetXaxis()->FindBin(fFitRangeMin);
1432 // int fitRangeBinMax = pHist_All_Data[novadaq::cnv::kNEARDET].GetXaxis()->FindBin(fFitRangeMax)-1;
1433 
1434  int dispRangeBinMin = pHist_All_Data[novadaq::cnv::kNEARDET].GetXaxis()->FindBin(fDisplayRangeMin);
1435  int dispRangeBinMax = pHist_All_Data[novadaq::cnv::kNEARDET].GetXaxis()->FindBin(fDisplayRangeMax);
1436 
1437  double rangeUserMin = data->GetBinLowEdge(dispRangeBinMin);
1438  double rangeUserMax = data->GetBinLowEdge(dispRangeBinMax)+ data->GetBinWidth(dispRangeBinMax);
1439 
1440  // Find fit range that incorporate all bins used in the fit
1441  double fitRangeUserMin = data->GetBinLowEdge(data->GetXaxis()->FindBin(fFitRangeMin));
1442  double fitRangeUserMax = data->GetBinLowEdge(data->GetXaxis()->FindBin(fFitRangeMax)-1)
1443  + data->GetBinWidth(data->GetXaxis()->FindBin(fFitRangeMax)-1);
1444 
1445  // Enforce these draw ranges on all histograms used in the draw
1446  data -> GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
1447  mc -> GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
1448 
1449  for( auto & tHistMap : pHistMap){
1450  md = tHistMap.first;
1451  if(!md.isMC || md.detector != pDet) continue;
1452  intHists[md.OscType()]->GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
1453  }
1454 
1455  // Make the stack
1456  stack = new THStack((name + "Stacks" + pSuffix).c_str(), title.c_str());
1457  title = fVarName + " (" + pSuffix + ") (Stacks)" + ((pDet == novadaq::cnv::kNEARDET) ? "ND" : "FD");
1458  name = fVarName + "_" + pSuffix + "_Stacks_" + ((pDet == novadaq::cnv::kNEARDET) ? "ND" : "FD") + id->getUniqID();
1459 
1460  // Make the canvas
1461  can = std::make_unique<TCanvas>(name.c_str(), title.c_str(), 800, 800);
1462  can->cd();
1463 
1464  // Set up the upper pad
1465  name = "Upper";
1466  TPad * padUpper = new TPad(name.c_str(), name.c_str(), 0, 0.3, 1, 1.0);
1467  padUpper->SetBottomMargin(0);
1468 
1469  //padLower.SetGridx(); // Jenny hates this, it seems. Boo Jenny.
1470  padUpper->Draw();
1471  padUpper->cd();
1472 
1473  // Make an index to organize the histograms by increasing integral
1474  std::vector<unsigned int> hist_index (intHists.size(), 0);
1475  std::vector<fnex::osc_type> hist_osctype(intHists.size(), kOscTypeNoneSet);
1476  std::vector<double> hist_size (intHists.size(), 0);
1477  unsigned int curIndex = 0;
1478  for( auto & intHist : intHists){
1479  hist_index [curIndex] = curIndex;
1480  hist_osctype[curIndex] = intHist.first;
1481  hist_size [curIndex] = intHist.second->Integral();
1482  ++curIndex;
1483  }
1484 
1485  sort(hist_index.begin(), hist_index.end(),
1486  [&](const int& a, const int& b) {
1487  return (hist_size[a] < hist_size[b]);
1488  }
1489  );
1490 
1491  // Add histogram to the stack in order of increasing integral, ignoring the signal
1492  // for now -- we'll add it last
1493  TH1D* hist = nullptr;
1494 
1495  for(auto histIdx : hist_index){
1496  hist = intHists[hist_osctype[histIdx]];
1497  max += hist->GetMaximum();
1498 
1499  for(int b = 0; b < hist->GetXaxis()->GetNbins(); ++b){
1500  hist->SetBinContent(b, hist->GetBinContent(b) * minWidth / hist->GetBinWidth(b));
1501  }
1502 
1503  if(std::find(signal_osctype.begin(),
1504  signal_osctype.end(),
1505  hist_osctype[histIdx]) != signal_osctype.end()) continue;
1506 
1507  stack->Add( hist );
1508  }
1509 
1510  // Add the signal last - the signal_osctype is only defined for the FD
1511  // and should have no entries for the ND, so this for loop is skipped in
1512  // that case
1513  for(auto sigT : signal_osctype){
1514  hist = intHists[sigT];
1515  max += hist->GetMaximum();
1516  stack->Add( hist );
1517  }
1518 
1519  for(int b = 0; b < data->GetXaxis()->GetNbins(); ++b){
1520  data->SetBinContent(b, data->GetBinContent(b) * minWidth / data->GetBinWidth(b));
1521  }
1522 
1523 
1524  // Make legend and add intHists entries to it
1525  TLegend * leg = new TLegend(0.525, 0.1, 0.875, 0.7);
1526  leg->AddEntry(data, "Data", "lep");
1527  for( auto & itr : intHists){
1528  leg->AddEntry(itr.second, fnex::cOscType_LatexStrings[itr.first].c_str(), "lfp" );
1529  }
1530 
1531  TLegend * leg_expname = new TLegend(0.525, 0.725, 0.875, 0.875);
1532  leg_expname->SetBorderSize(1);
1533  leg_expname->SetFillColor(0);
1534  std::string newheader = pSuffix;
1535  boost::trim(newheader);
1536  leg_expname->SetHeader(newheader.data());
1537 
1538  // Draw stack so that we can start manipulating it
1539  stack->Draw("HIST");
1540  stack->GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
1541  stack->GetYaxis()->SetTitle("Events");
1542  stack->GetYaxis()->CenterTitle();
1543  stack->GetYaxis()->SetTitleSize(20);
1544  stack->GetYaxis()->SetTitleFont(43);
1545  stack->GetYaxis()->SetTitleOffset(1.55);
1546  stack->GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
1547  stack->GetYaxis()->SetLabelSize(15);
1548 
1549  stack->Draw("HIST");
1550 
1551  // Get/set the max range
1552  max = stack->GetMaximum();
1553  if(max < data->GetMaximum()) max = data->GetMaximum();
1554  data ->SetMaximum(max);
1555  stack->SetMaximum(max);
1556 
1557  // Make lines to show fit range
1558  double lineYMin = 0;
1559  double lineYMax = max;
1560  TLine * fitMinLine = new TLine(fitRangeUserMin, lineYMin, fitRangeUserMin, lineYMax);
1561  TLine * fitMaxLine = new TLine(fitRangeUserMax, lineYMin, fitRangeUserMax, lineYMax);
1562  fitMinLine->SetLineColor(kRed);
1563  //fitMinLine.SetLineStyle(6);
1564  fitMinLine->SetLineWidth(2);
1565  fitMaxLine->SetLineColor(kRed);
1566  //fitMaxLine.SetLineStyle(6);
1567  fitMaxLine->SetLineWidth(2);
1568  leg->AddEntry(fitMinLine, "FitRegion", "l");
1569 
1570  // Draw everything for final presentation
1571  stack->Draw("HIST");
1572  data ->Draw("E1 SAME");
1573  if(this->DrawLowerFitLine()) fitMinLine->Draw("SAME");
1574  if(this->DrawUpperFitLine()) fitMaxLine->Draw("SAME");
1575  leg->Draw();
1576  leg_expname->Draw();
1577 
1578  // Update the canvas to reflect these changes
1579  can->Update();
1580 
1581  // Let's go back to the parent canvas so that we can move to the lower pad
1582  can->cd();
1583 
1584  // Set up the lower pad
1585  name = "Lower";
1586  TPad * padLower = new TPad(name.c_str(), name.c_str(), 0, 0.05, 1, 0.3);
1587  padLower->SetTopMargin(0);
1588  padLower->SetBottomMargin(0.25);
1589  //padLower.SetGridx(); // Jenny hates this, it seems. Boo Jenny.
1590  padLower->Draw();
1591  padLower->cd();
1592 
1593  // Draw the ratio
1594  TH1D * pDMCRat = new TH1D(pHist_All_Data_Over_MC[pDet]);
1595  pDMCRat->SetLineColor(kBlack);
1596  pDMCRat->SetLineWidth(2);
1597  pDMCRat->SetMarkerColor(kBlack);
1598  pDMCRat->SetTitle("");
1599 
1600  pDMCRat->GetYaxis()->SetTitle("Data/MC");
1601  pDMCRat->GetYaxis()->CenterTitle();
1602  pDMCRat->GetYaxis()->SetNdivisions(505);
1603  pDMCRat->GetYaxis()->SetTitleSize(20);
1604  pDMCRat->GetYaxis()->SetTitleFont(43);
1605  pDMCRat->GetYaxis()->SetTitleOffset(1.55);
1606  pDMCRat->GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
1607  pDMCRat->GetYaxis()->SetLabelSize(15);
1608 
1609  pDMCRat->GetXaxis()->SetTitle(data->GetXaxis()->GetTitle());
1610  pDMCRat->GetXaxis()->SetTitleSize(20);
1611  pDMCRat->GetXaxis()->SetTitleFont(43);
1612  pDMCRat->GetXaxis()->SetTitleOffset(4.);
1613  pDMCRat->GetXaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
1614  pDMCRat->GetXaxis()->SetLabelSize(15);
1615 
1616  pDMCRat->SetStats(false);
1617  pDMCRat->GetXaxis()->SetRangeUser(rangeUserMin, rangeUserMax);
1618  pDMCRat->Draw("E1");
1619 
1620  // Set up the Equal1 line
1621  TLine * yEquals1 = new TLine;
1622  yEquals1->SetLineWidth(1);
1623  yEquals1->SetLineStyle(2);
1624  yEquals1->SetLineColor(kBlack);
1625  yEquals1->DrawLine(rangeUserMin, 1.0,
1626  rangeUserMax, 1.0 );
1627 
1628 
1629  // Set up the fit range lines
1630  // KM June 10 2016
1631  double smallLineYMin = pDMCRat->GetMinimum();
1632  double smallLineYMax = pDMCRat->GetMaximum();
1633  TLine * fitMinLineSmall = new TLine(fitRangeUserMin, smallLineYMin, fitRangeUserMin, smallLineYMax);
1634  TLine * fitMaxLineSmall = new TLine(fitRangeUserMax, smallLineYMin, fitRangeUserMax, smallLineYMax);
1635  fitMinLineSmall->SetLineColor(kRed);
1636  //fitMinLineSmall.SetLineStyle(6);
1637  fitMinLineSmall->SetLineWidth(2);
1638  fitMaxLineSmall->SetLineColor(kRed);
1639  //fitMaxLineSmall.SetLineStyle(6);
1640  fitMaxLineSmall->SetLineWidth(2);
1641  if(this->DrawLowerFitLine()) fitMinLineSmall->Draw("SAME");
1642  if(this->DrawUpperFitLine()) fitMaxLineSmall->Draw("SAME");
1643 
1644  // Update the canvas with the ratio pad
1645  can->Update();
1646 
1647  //can->Write();
1648 
1649  // We're done here
1650  return can;
1651 
1652  }//end of method
1653 
1654 
1655  //----------------------------------------------------------------------------
1657  std::string pSuffix){
1658 
1659 
1660  // Make the stack
1662  TString title = fVarName + " (" + pSuffix + ") (Stacks)" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD");
1663  TString name = fVarName + "_" + pSuffix + "_Stacks_" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD") + id->getUniqID();
1664 
1665  // Make the canvas
1666  std::unique_ptr<TCanvas> can = std::make_unique<TCanvas>(name.Data(), title.Data(), 800, 800);
1667  can->cd();
1668 
1669  // Make parameter legend - we are not adding chi^2 information to it
1670  // (last parameter is false) so first two parameters don't matter
1671  TPaveText *oscLeg = new TPaveText(this->MakeOscillationLegend(0., 0, det, false));
1672 
1673  // We're done here. Draw the canvas describing par / syst space
1674  // and update the canvas
1675  oscLeg->Draw();
1676  can->Update();
1677 
1678  return can;
1679 
1680 
1681  }
1682 
1683  //----------------------------------------------------------------------------
1685  int ndof,
1686  novadaq::cnv::DetId const& det,
1687  bool addChiSqr)
1688  {
1689  TPaveText pointLegend(0.1, 0.1, 0.9, 0.9);
1690  pointLegend.SetBorderSize(1);
1691  pointLegend.SetFillColor(0);
1692 
1693  if(addChiSqr){
1694  std::stringstream ss_chi;
1695  ss_chi << "#chi^{2} / NDOF = " << chisq << " / " << ndof;
1696  pointLegend.AddText(ss_chi.str().data());
1697  pointLegend.AddText(" ");
1698  }
1699 
1700  pointLegend.AddText("Oscillation Parameters");
1701 
1702  // Describe par space, including scale factors
1703  // and whether they were free in the fit
1704  std::string scale = "";
1705  std::string parName = "";
1706  std::string freeFix = "";
1707  double scaleFac = 1.;
1708  double result = 0.;
1709 
1710  auto const& fixedParameters = fRefInputPoint.FixedParameters(det);
1711 
1712  for(const auto & oscParam : fRefInputPoint.ParameterSpaceLocation(det)){
1713 
1714  scale = "";
1715  freeFix = ( !fRefInputPoint.IsParameterFixed(oscParam.first, det) ) ? " (FREE)" : "";
1716  scaleFac = 1.;
1717  result = oscParam.second;
1718  parName = fOscParmNameToLatex[oscParam.first];
1719 
1720  LOG_DEBUG("CorrectedSpectrum")
1721  << parName
1722  << " is fixed: "
1723  << freeFix;
1724 
1725 
1726  if(oscParam.first.compare("Dmsq32") == 0){
1727  scale = " #times 10^{-3}";
1728  scaleFac = 1.e3;
1729  }
1730  else if(oscParam.first.compare("Dmsq21") == 0){
1731  scale = " #times 10^{-5}";
1732  scaleFac = 1.e5;
1733  }
1734  std::stringstream ss;
1735 
1736  // Do not print the NSI parameters we are not sensitive to
1737  if(oscParam.first.compare("Eps_ee") == 0 ||
1738  oscParam.first.compare("Eps_mumu") == 0 ||
1739  oscParam.first.compare("Eps_tautau") == 0 ||
1740  oscParam.first.compare("Eps_mutau") == 0 ||
1741  oscParam.first.compare("Delta_mutau") == 0) continue;
1742 
1743  if(oscParam.first.compare("dCP") == 0 ||
1744  oscParam.first.compare("Delta_emu") == 0 ||
1745  oscParam.first.compare("Delta_etau") == 0 ||
1746  oscParam.first.compare("Delta_mutau") == 0){
1747  ss << std::setw(5)
1748  << parName
1749  << " = "
1750  << std::setprecision(3)
1751  << std::fixed
1752  << std::setw(5)
1753  << result * scaleFac
1754  << scale
1755  << " = "
1756  << std::setprecision(3)
1757  << result/TMath::Pi()
1758  << "#pi"
1759  << freeFix;
1760  }
1761  else if(oscParam.first.compare("dCP") != 0 ||
1762  oscParam.first.compare("Delta_emu") != 0 ||
1763  oscParam.first.compare("Delta_etau") != 0 ||
1764  oscParam.first.compare("Delta_mutau") != 0){
1765  ss << std::setw(5)
1766  << parName
1767  << " = "
1768  << std::setprecision(3)
1769  << std::fixed
1770  << std::setw(5)
1771  << result * scaleFac
1772  << scale
1773  << freeFix;
1774  }
1775 
1776  pointLegend.AddText(ss.str().data());
1777 
1778  if(oscParam.first.compare("Th23") == 0){
1779  std::stringstream extra;
1780  extra
1781  << std::setw(5)
1782  << "sin^{2}(#theta_{23}) = "
1783  << std::setprecision(3)
1784  << std::fixed
1785  << std::sin(result) * std::sin(result)
1786  << freeFix;
1787  pointLegend.AddText(extra.str().data());
1788  }
1789 
1790  }
1791 
1792  // Describe systematics space in the same way
1793  if(fixedParameters.size() == 0) pointLegend.AddText("None");
1794 
1795  return pointLegend;
1796  }
1797 
1798  //----------------------------------------------------------------------------
1800  int ndof,
1801  novadaq::cnv::DetId const& det,
1802  bool addChiSqr)
1803  {
1804  TPaveText pointLegend(0.1, 0.1, 0.9, 0.9);
1805  pointLegend.SetBorderSize(1);
1806  pointLegend.SetFillColor(0);
1807 
1808  if(addChiSqr){
1809  std::stringstream ss_chi;
1810  ss_chi << "#chi^{2} / NDOF = " << chisq << " / " << ndof;
1811  pointLegend.AddText(ss_chi.str().data());
1812  pointLegend.AddText(" ");
1813  }
1814 
1815  pointLegend.AddText(" ");
1816  pointLegend.AddText("Systematic Shifts");
1817 
1818  auto const& fixedSystematics = fRefInputPoint.FixedSystematics(det);
1819 
1820  std::string freeFix = "";
1821 
1822  for (const auto & sysPull : fRefInputPoint.SystematicPulls(det)){
1823 
1824  freeFix = ( !fRefInputPoint.IsParameterFixed(sysPull.first, det) ) ? " (FREE)" : "";
1825 
1826  std::stringstream ss;
1827  ss
1828  << std::setw(15)
1829  << sysPull.first
1830  << " = "
1831  << std::setprecision(3)
1832  << std::fixed
1833  << std::setw(10)
1834  << sysPull.second
1835  << freeFix;
1836 
1837  pointLegend.AddText(ss.str().data());
1838  }
1839 
1840  if(fixedSystematics.size() == 0) pointLegend.AddText("None");
1841 
1842  return pointLegend;
1843  }
1844 
1845  //----------------------------------------------------------------------------
1847  std::string pSuffix){
1848 
1849  // Make the stack
1851  TString title = fVarName + " (" + pSuffix + ") (Stacks)" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD");
1852  TString name = fVarName + "_" + pSuffix + "_Stacks_" + ((det == novadaq::cnv::kNEARDET) ? "ND" : "FD") + id->getUniqID();
1853 
1854  // Make the canvas
1855  std::unique_ptr<TCanvas> can = std::make_unique<TCanvas>(name.Data(), title.Data(), 800, 800);
1856  can->cd();
1857 
1858  // Make parameter legend - we are not adding chi^2 information to it
1859  // (last parameter is false) so first two parameters don't matter
1860  TPaveText *sysLeg = new TPaveText(this->MakeSystematicsLegend(0., 0, det, false));
1861 
1862  // We're done here. Draw the canvas describing par / syst space
1863  // and update the canvas
1864  sysLeg->Draw();
1865  can->Update();
1866 
1867  return can;
1868 
1869  }
1870 
1871  // Spill fCorrectedHistMap contents
1872  //----------------------------------------------------------------------------
1874  std::string pSuffix){
1875  std::string newSuffix = pSuffix;
1876  newSuffix.append("_Corr");
1881  det,
1882  newSuffix );
1883  }
1884 
1885  //----------------------------------------------------------------------------
1887  std::string pSuffix){
1888 
1889  std::string newSuffix = pSuffix;
1890  newSuffix.append("_Uncorr");
1895  det,
1896  newSuffix );
1897 
1898  }
1899 
1900 
1901 }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
const XML_Char * name
Definition: expat.h:151
const std::string cOscType_LatexStrings[20]
Definition: Constants.h:292
bool IsNuESelected() const
Definition: Structs.cxx:265
DetHistMap fUncorrectedHist_All_Data_Over_MC
bool IsNuMuSelected() const
Definition: Structs.cxx:256
TSpline3 DrawSpline(TH1D *mc, double fitRangeMinUser, double fitRangeMaxUser)
std::unique_ptr< TCanvas > GetStacksCanvasCopy(HistMap &pHistMap, DetHistMap &pHist_All_Data, DetHistMap &pHist_All_MC, DetHistMap &pHist_All_Data_Over_MC, novadaq::cnv::DetId &det, std::string pSuffix)
std::unique_ptr< TCanvas > GetOscParCanvasCopy(novadaq::cnv::DetId &det, std::string pSuffix)
T sqrt(T number)
Definition: d0nt_math.hpp:156
const std::string cOscType_Strings[20]
Definition: Constants.h:282
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
constexpr T pow(T x)
Definition: pow.h:75
Float_t ss
Definition: plot.C:24
string trim(string in)
Definition: rootgINukeVal.C:65
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void SetDisplayRange(double pMin, double pMax)
long const EpochKey() const
Definition: Structs.h:58
const int cOscType_Colors[20]
Definition: Constants.h:313
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::string const EpochString() const
Definition: Structs.cxx:101
Create a list of fnex::Events to be used in fits.
std::unique_ptr< TCanvas > GetSystCanvasCopy(novadaq::cnv::DetId &det, std::string pSuffix)
void DrawDataVsMC(art::TFileDirectory *tfd, DetHistMap &pCorrectedHist_All_Data, DetHistMap &pCorrectedHist_All_MC, DetHistMap &pCorrectedHist_All_Data_Over_MC, std::string suffix)
std::vector< std::string > const & FixedParameters(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:396
std::map< std::string, std::string > fOscParmNameToLatex
const int cOscType_FillStyles[20]
Definition: Constants.h:327
const XML_Char const XML_Char * data
Definition: expat.h:268
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
bool IsNuMuQuantiles() const
Definition: Structs.cxx:247
std::vector< std::string > const & FixedSystematics(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:404
const int nbins
Definition: cellShifts.C:15
Double_t scale
Definition: plot.C:25
std::string const ToString() const
Definition: Structs.cxx:114
Far Detector at Ash River, MN.
virtual std::size_t nBins() const =0
Number of bins always means something, but data is stored differently in derived classes.
TPaveText MakeOscillationLegend(double chisq, int ndof, novadaq::cnv::DetId const &det, bool addChiSqr=true)
ParameterMap const ParameterSpaceLocation(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:330
const double a
correl_xv GetXaxis() -> SetDecimals()
T get(std::string const &key) const
Definition: ParameterSet.h:231
DetHistMap fCorrectedHist_All_Data_Over_MC
long const DetectorKey() const
Definition: Structs.h:57
std::unordered_map< novadaq::cnv::DetId, TH1D, fnex::DetectorHasher > DetHistMap
Near Detector in the NuMI cavern.
void SetFitRange(double pMin, double pMax)
static FNEXUniqID * getInstance()
Definition: FNEXUniqID.cxx:14
novadaq::cnv::DetId detector
Definition: Structs.h:50
static VarBinning * Instance()
Definition: VarBinning.cxx:15
float bin[41]
Definition: plottest35.C:14
std::unique_ptr< TCanvas > GetCorrectedStacksCanvasCopy(novadaq::cnv::DetId &det, std::string pSuffix)
void DrawCounts(art::TFileDirectory *tfd, HistMap &pHistMap, DetHistMap &pHist_All_Data, DetHistMap &pHist_All_MC, DetHistMap &pHist_All_Data_Over_MC, std::string pSuffix)
std::string const DetectorString() const
Definition: Structs.h:93
fnex::OscType_t const OscType() const
Definition: Structs.cxx:190
long const MCKey() const
Definition: Structs.h:56
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
fnex::SelectionType_t selectionType
Definition: Structs.h:52
bool IsParameterFixed(ParameterDet const &pardet) const
Definition: InputPoint.cxx:468
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
T sin(T number)
Definition: d0nt_math.hpp:132
CorrectedSpectrum(fhicl::ParameterSet const &pset)
void DrawStacks(art::TFileDirectory *tfd, HistMap &pHistMap, DetHistMap &pHist_All_Data, DetHistMap &pHist_All_MC, DetHistMap &pHist_All_Data_Over_MC, std::string pSuffix)
std::unordered_map< fnex::MetaData, TH1D, fnex::MetaDataHasher > HistMap
const std::string cOscParams_Strings[kNumOscParams]
Definition: Constants.h:196
fnex::Binning const * SpectrumBinning(fnex::MetaData const &md)
const hit & b
Definition: hits.cxx:21
const std::string cOscParams_Strings_Latex[kNumOscParams]
Definition: Constants.h:205
TPad * pad2
Definition: analysis.C:13
T min(const caf::Proxy< T > &a, T b)
std::unique_ptr< TCanvas > GetUncorrectedStacksCanvasCopy(novadaq::cnv::DetId &det, std::string pSuffix)
TPaveText MakeSystematicsLegend(double chisq, int ndof, novadaq::cnv::DetId const &det, bool addChiSqr=true)
Base class for &#39;binning&#39; objects. Pure virtual – you want either a FixedBinning or a VariableBinning...
Definition: Binning.h:20
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
virtual std::vector< double > const & BinLowEdges() const =0
fnex::Binning const * BinInfo(std::string const &varName)
Definition: VarBinning.cxx:65
c cd(1)
void DrawCorrVsUncorr(art::TFileDirectory *tfd, DetHistMap &pHist_All_Data, DetHistMap &pHist_All_MC, DetHistMap &pHist_All_Data_Over_MC, std::string pSuffix)
TPad * pad1
Definition: analysis.C:13
ParameterMap const SystematicPulls(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:364