ProtoDUNEFitUtils.cxx
Go to the documentation of this file.
1 #include "ProtoDUNEFitUtils.h"
2 
3 #include <TIterator.h>
4 #include <TLegend.h>
5 #include <TLegendEntry.h>
6 #include <TLine.h>
7 #include <TF1.h>
8 #include <TFile.h>
9 #include <TGaxis.h>
10 #include <TKey.h>
11 
12 #include <RooDataSet.h>
13 #include <RooSimultaneous.h>
14 #include <RooCategory.h>
15 #include <RooRealIntegral.h>
16 #include <RooRealVar.h>
17 #include <RooAbsReal.h>
18 #include <RooAbsArg.h>
19 #include <Roo1DTable.h>
20 #include <RooRealSumPdf.h>
21 #include <RooProduct.h>
22 #include <RooProdPdf.h>
23 #include <RooHist.h>
24 
25 #include <RooStats/HistFactory/PiecewiseInterpolation.h>
26 #include <RooStats/ModelConfig.h>
27 
28 //********************************************************************
29 double GetArgonNumberDensity(double argon_density, double argon_molecularmass){
30  //********************************************************************
31 
32  // argon_density: g/cm3, argon_molecularmass: g/mol
33  double avogadro_number = 6.0221*10e23;
34  return (avogadro_number*argon_density/argon_molecularmass);
35 
36 }
37 
38 //********************************************************************
39 std::vector<TH1*> protoana::ProtoDUNEFitUtils::GetSystHistograms(std::string name){
40  //********************************************************************
41 
42  std::vector<TH1*> vec;
43 
44  TFile* f = new TFile(name.c_str(), "READ");
45  if(!f) return vec;
46 
47  TIter next(f->GetListOfKeys());
48  TKey *key;
49 
50  while ((key = (TKey*)next())){
51  TString classname(key->GetClassName());
52  if(!classname.Contains("TH1")) continue;
53 
54  TH1* h = (TH1*)key->ReadObj();
55  h->SetDirectory(0);
56  vec.push_back(h);
57  }
58 
59  f->Close();
60 
61  return vec;
62 
63 }
64 
65 //********************************************************************
66 TH1* protoana::ProtoDUNEFitUtils::GetStatsSystHistogram(TH1* nominal){
67  //********************************************************************
68 
69  TString name = TString(nominal->GetName()) + TString("_mcshapestats");
70  TH1* histo = (TH1*)nominal->Clone(name.Data());
71 
72  for(Int_t i=1; i <= nominal->GetNbinsX(); i++){
73  Double_t mcerror = sqrt(nominal->GetBinContent(i));
74 
75  Double_t ratio = 0.0;
76  if(nominal->GetBinContent(i) > 0.)
77  ratio = mcerror/nominal->GetBinContent(i);
78 
79  histo->SetBinContent(i, ratio);
80  histo->SetBinError(i, 0.0);
81  }
82 
83  return histo;
84 
85 }
86 
87 //********************************************************************
88 TH1* protoana::ProtoDUNEFitUtils::GetSystematicHistoFromNominal(TH1* nominal, TH1* syst, TString name){
89  //********************************************************************
90 
91  TString histoname = TString(syst->GetName()) + TString("_") + name;
92  TH1D* histo = new TH1D(histoname.Data(), histoname.Data(), nominal->GetNbinsX(), nominal->GetBinLowEdge(1), nominal->GetBinLowEdge(nominal->GetNbinsX()+1));
93  if(nominal->GetNbinsX() != syst->GetNbinsX()){
94  std::cout << "WARNING::Nominal and systematic histograms does not have the same number of bins. No systematic is applied!" << std::endl;
95  return histo;
96  }
97 
98  for(Int_t i=0; i < nominal->GetNbinsX(); i++){
99  double nom = nominal->GetBinContent(i+1);
100  double sys = syst->GetBinContent(i+1);
101 
102  if(nom <= 0.){
103  std::cout << "WARNING::No events in histogram " << nominal->GetName() << " for bin " << i << ". Skip systematics propagation for systematic histogram " << syst->GetName() << std::endl;
104  histo->SetBinContent(i+1, 0.0);
105  continue;
106  }
107 
108  if(sys < 0.)
109  std::cout << "WARNING::Negative fractional error for histogram " << syst->GetName() << " and bin " << i << std::endl;
110 
111  // Avoid large fluctuation in the systematics
112  if(sys > 1.0) sys = 1.0;
113 
114  double w = 1.0;
115  if(name == "up" || name == "Up" || name == "UP" || name == "UP2" || name == "HIGH" || name == "High" || name == "high"){
116  w = nom + nom*sys;
117  }
118  else if(name == "down" || name == "Down" || name == "DOWN" || name == "DOWN2" || name == "LOW" || name == "Low" || name == "low"){
119  w = nom - nom*sys;
120  }
121 
122  if(w < 0.) w = 0.;
123 
124  histo->SetBinContent(i+1, w);
125  }
126 
127  return histo;
128 
129 }
130 
131 //********************************************************************
132 bool protoana::ProtoDUNEFitUtils::IsSingleBinHisto(TH1* histo){
133  //********************************************************************
134 
135  int nbinswithevents = 0;
136  for(int j = 0; j < histo->GetNbinsX(); j++){
137  if(histo->GetBinContent(j) > 0.0)
138  nbinswithevents++;
139  }
140 
141  if(nbinswithevents == 1) return true;
142 
143  return false;
144 
145 }
146 
147 //********************************************************************
148 TH2* protoana::ProtoDUNEFitUtils::GetFitCovariance(RooFitResult* result){
149  //********************************************************************
150 
151  RooArgList floatParsFinal = result->floatParsFinal();
152  const TMatrixDSym& covmatrix = result->covarianceMatrix();
153  Int_t n = covmatrix.GetNcols();
154  TH2D* CovarianceHisto = new TH2D("FitCovariance", "FitCovariance", n,0,n, n,0,n);
155  CovarianceHisto->GetXaxis()->SetLabelSize(0.01);
156  CovarianceHisto->GetYaxis()->SetLabelSize(0.01);
157  for(Int_t j = 0 ; j<n ; j++) {
158  for(Int_t k = 0 ; k<n; k++) {
159  CovarianceHisto->Fill(j+0.5,n-k-0.5, covmatrix(j,k));
160  }
161  CovarianceHisto->GetXaxis()->SetBinLabel(j+1, floatParsFinal.at(j)->GetName());
162  CovarianceHisto->GetYaxis()->SetBinLabel(n-j, floatParsFinal.at(j)->GetName());
163  }
164  CovarianceHisto->SetMinimum(-1);
165  CovarianceHisto->SetMaximum(+1);
166 
167  return CovarianceHisto;
168 
169 }
170 
171 //********************************************************************
172 TH2* protoana::ProtoDUNEFitUtils::GetFitCorrelationMatrix(RooFitResult* result){
173  //********************************************************************
174 
175  TH2* corrhisto = result->correlationHist();
176  corrhisto->GetXaxis()->SetLabelSize(0.01);
177  corrhisto->GetYaxis()->SetLabelSize(0.01);
178 
179  return corrhisto;
180 
181 }
182 
183 //********************************************************************
184 void protoana::ProtoDUNEFitUtils::SaveSnapshot(RooWorkspace* ws, TString snapshotname){
185  //********************************************************************
186 
187  if(!ws){
188  std::cerr << "ERROR::Workspace not found. Will not save snapshot." << std::endl;
189  return;
190  }
191 
192  RooSimultaneous* simpdf = (RooSimultaneous*)ws->pdf("simPdf");
193  if(!simpdf)
194  simpdf = (RooSimultaneous*)ws->pdf("combPdf");
195  if(!simpdf){
196  std::cerr << "ERROR::Pdf not found in workspace. Will not save snapshot." << std::endl;
197  return;
198  }
199 
200  RooAbsData* obsdata = (RooAbsData*)ws->data("obsData");
201  RooArgSet* parameters = (RooArgSet*)simpdf->getParameters(*obsdata);
202  if(!ws->loadSnapshot(snapshotname.Data())){
203  ws->saveSnapshot(snapshotname.Data(),*parameters);
204  }
205  else{
206  std::cout << "WARNING::Snapshot " << snapshotname.Data() << " found in workspace. Will not overwrite it." << std::endl;
207  }
208 
209  //gDirectory->Add(ws);
210 
211 }
212 
213 //********************************************************************
214 bool protoana::ProtoDUNEFitUtils::LoadSnapshot(RooWorkspace* ws, TString snapshotname){
215  //********************************************************************
216 
217  if(!ws){
218  std::cerr << "ERROR::Workspace not found. Can't find snapshot." << std::endl;
219  return false;
220  }
221 
222  if(!ws->loadSnapshot(snapshotname)){
223  std::cerr << "ERROR::Snapshot not found in workspace." << std::endl;
224  return false;
225  }
226 
227  ws->loadSnapshot(snapshotname);
228  std::cout << "INFO::Workspace snapshot " << snapshotname.Data() << " loaded." << std::endl;
229 
230  return true;
231 
232 }
233 
234 //********************************************************************
235 void protoana::ProtoDUNEFitUtils::SaveWorkspace(RooWorkspace* ws, TString outFileName){
236  //********************************************************************
237 
238  if(!ws){
239  std::cerr << "ERROR::Workspace not found. Won't save." << std::endl;
240  return;
241  }
242 
243  ws->writeToFile(outFileName.Data());
244  std::cout << "INFO::Workspace written to file " << outFileName.Data() << std::endl;
245 
246  return;
247 
248 }
249 
250 //********************************************************************
251 void protoana::ProtoDUNEFitUtils::SetInterpolationCode(RooWorkspace* ws, Int_t code){
252  //********************************************************************
253 
254  if(!ws){
255  std::cerr << "ERROR::NULL workspace. No interpolation code added." << std::endl;
256  return;
257  }
258 
259  RooArgSet fun = ws->allFunctions();
260  TIterator* iter = fun.createIterator();
261 
262  RooAbsArg* arg(0);
263  while( (arg=(RooAbsArg*)iter->Next()) ) {
264  if(arg->ClassName()!=TString("PiecewiseInterpolation") ) continue;
265  PiecewiseInterpolation* p = (PiecewiseInterpolation*)ws->function(arg->GetName() );
266  p->setAllInterpCodes(code);
267  }
268 
269  delete iter;
270 
271 }
272 
273 //********************************************************************
274 void protoana::ProtoDUNEFitUtils::RemoveEmptyBins(RooPlot* frame){
275  //********************************************************************
276 
277  const char* histname = 0;
278 
279  // Find histogram
280  RooHist* histo = (RooHist*) frame->findObject(histname, RooHist::Class());
281  if(!histo) return;
282 
283  for(Int_t i=0; i < histo->GetN(); i++){
284  Double_t x,y;
285  histo->GetPoint(i,x,y);
286 
287  if(fabs(y) < 0.00001 && histo->GetErrorYhigh(i) > 0.){
288  histo->RemovePoint(i);
289  if(i != histo->GetN()) --i;
290  }
291  }
292 
293  return;
294 
295 }
296 
297 //********************************************************************
298 double protoana::ProtoDUNEFitUtils::GetDataMCChi2(RooWorkspace *work, TString channelname, RooAbsData* data){
299  //********************************************************************
300 
301  if(!work){
302  std::cerr << "ERROR::NULL workspace. No chi2 is computed!" << std::endl;
303  return -1.0;
304  }
305 
306  // Silence output
307  RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);
308  RooMsgService::instance().getStream(1).removeTopic(RooFit::Plotting);
309 
310  // Get pdf from workspace
311  RooSimultaneous* pdf = (RooSimultaneous*)work->pdf("simPdf");
312  if(!pdf){
313  std::cerr << "ERROR::No pdf found in workspace. No chi2 is computed!" << std::endl;
314  return -1.0;
315  }
316 
317  // If not provided, get data from workspace
318  if(!data) data = work->data("obsData");
319 
320  // Get category components
321  RooCategory* categories = work->cat("channelCat");
322 
323  std::vector<TString> categoriesName;
324  TIterator* iter = categories->typeIterator();
325  RooCatType* catType;
326  while( (catType = (RooCatType*) iter->Next())) {
327  TString catname = catType->GetName();
328  categoriesName.push_back(catname);
329  }
330 
331  double chi2 = -999.0;
332  for(UInt_t i = 0; i < categoriesName.size(); i++){
333  TString catname = categoriesName[i];
334 
335  if(categories->setLabel(catname)){
336  std::cout << "WARNING::Category " << catname.Data() << " is not a member of channelCat. Will skip." << std::endl;
337  continue;
338  }
339 
340  if(catname == channelname){
341 
342  RooAbsPdf* subpdf = (RooAbsPdf*)pdf->getPdf(catname.Data());
343  if(!subpdf){
344  std::cout << "WARNING::Can't find sub-pdf for region " << catname.Data() << ". Will skip." << std::endl;
345  continue;
346  }
347 
348  TString subdataset_str = Form("channelCat==channelCat::%s",catname.Data());
349  RooAbsData* subdataset = (RooAbsData*) data->reduce(subdataset_str.Data());
350  if(!subdataset){
351  std::cout << "WARNING::Can't find sub-dataset for region " << catname.Data() << ". Will skip." << std::endl;
352  continue;
353  }
354 
355  RooRealVar* var =(RooRealVar*) ((RooArgSet*) subpdf->getObservables(*subdataset))->find(Form("obs_x_%s", catname.Data()));
356  RooPlot* frame = var->frame();
357  frame->SetName(Form("chi2frame_%s", catname.Data()));
358 
359  subdataset->plotOn(frame, RooFit::DataError(RooAbsData::Poisson), RooFit::MarkerColor(1), RooFit::LineColor(1));
360  subpdf->plotOn(frame,RooFit::Normalization(1,RooAbsReal::RelativeExpected),RooFit::Precision(1e-5));
361  chi2 = frame->chiSquare();
362  break;
363  }
364 
365  }
366 
367  return chi2;
368 
369 }
370 
371 std::vector<TH1 *> protoana::ProtoDUNEFitUtils::PlotXSecs(
372  RooWorkspace * work, std::string name, /*std::string error,*/
373  std::vector<TString> binnames, std::vector<double> recobins,
374  std::vector<TString> incidentBinNames,
375  RooAbsData * data,
376  RooFitResult * result) {
377 
378  std::vector<TH1 *> xsecs;
379 
380  if (!work) {
381  std::cerr << "ERROR:NULL dir. Will return empty canvas" << std::endl;
382  return xsecs;
383  }
384 
385  // Silence output
386  RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);
387  RooMsgService::instance().getStream(1).removeTopic(RooFit::Plotting);
388 
389  // Get pdf from workspace
390  RooSimultaneous* pdf = (RooSimultaneous*)work->pdf("simPdf");
391  if(!pdf){
392  std::cerr << "ERROR::No pdf found in workspace. Will return empty vector!"
393  << std::endl;
394  std::cerr << "ERROR::No pdf found in workspace. Will return empty vector!"
395  << std::endl;
396  std::cerr << "ERROR::No pdf found in workspace. Will return empty vector!"
397  << std::endl;
398  return xsecs;
399  }
400 
401 
402 
403  // Get category components
404  // i.e. Incident, Abs, Cex
405  RooCategory* categories = work->cat("channelCat");
406 
407  std::vector<TString> categoriesName;
408  TIterator* iter = categories->typeIterator();
409  RooCatType* catType;
410  while( (catType = (RooCatType*) iter->Next())) {
411  TString catname = catType->GetName();
412  categoriesName.push_back(catname);
413  std::cout << catname << std::endl;
414  }
415 
416  for (size_t i = 0; i < categoriesName.size(); ++i) {
417 
418  TString catname = categoriesName[i];
419 
420  RooAbsPdf* subpdf = (RooAbsPdf*)pdf->getPdf(catname.Data());
421  if(!subpdf){
422  std::cout << "WARNING::Can't find sub-pdf for region " << catname.Data()
423  << ". Will skip." << std::endl;
424  continue;
425  }
426 
427  TString RRSumPdfName = Form("%s_model",catname.Data());
428  RooRealSumPdf* RRSumPdf =
429  (RooRealSumPdf*)subpdf->getComponents()->find(RRSumPdfName);
430  RooArgList RRSumComponentsList = RRSumPdf->funcList();
431  RooLinkedListIter iter = RRSumComponentsList.iterator();
432  RooProduct* component;
433 
434  std::vector<TString> compNameVec;
435  while( (component = (RooProduct*) iter.Next()) ){
436  TString componentName = component->GetName();
437  std::cout << componentName << std::endl;
438  }
439  }
440 
441  /*
442  TH1 * incident_signal = 0x0;
443  std::vector<TH1 *> incident_backgrounds;
444 
445  for (const auto && key : *keys) {
446  std::string name = key->GetName();
447 
448  auto find_Incident = name.find("MC_ChannelIncident_Pions");
449  if (find_Incident != std::string::npos) {
450  incident_signal = (TH1*)key->Clone();
451  }
452  else {
453  incident_backgrounds.push_back((TH1*)key->Clone());
454  }
455  }
456  std::cout << incident_signal << " " << incident_backgrounds.size() << std::endl;
457 
458  for (const std::string & main_channel : {"ABS", "CEX"}) {
459  std::cout << main_channel << std::endl;
460 
461  for (const auto && key : *keys) {
462  std::string name = key->GetName();
463  std::cout << name << std::endl;
464  }
465  }
466  */
467 
468  return xsecs;
469 }
470 
471 //********************************************************************
472 std::vector<TCanvas*> protoana::ProtoDUNEFitUtils::PlotDatasetsAndPdfs(
473  RooWorkspace *work, TString name, TString error, TString plottodraw,
474  std::vector<TString> binnames, std::vector<double> recobins,
475  std::vector<TString> incidentBinNames,
476  std::vector<TString> sidebandBinNames, std::vector<double> sidebandBins,
477  TString measurement,
478  bool doNegativeReco, RooAbsData* data, RooFitResult* result) {
479 //********************************************************************
480 
481  std::vector<TCanvas*> rooplots;
482 
483  if(!work){
484  std::cerr << "ERROR::NULL workspace. Will return empty vector!" << std::endl;
485  return rooplots;
486  }
487 
488  // Silence output
489  RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);
490  RooMsgService::instance().getStream(1).removeTopic(RooFit::Plotting);
491 
492  // Get pdf from workspace
493  RooSimultaneous* pdf = (RooSimultaneous*)work->pdf("simPdf");
494  if(!pdf){
495  std::cerr << "ERROR::No pdf found in workspace. Will return empty vector!" << std::endl;
496  std::cerr << "ERROR::No pdf found in workspace. Will return empty vector!" << std::endl;
497  std::cerr << "ERROR::No pdf found in workspace. Will return empty vector!" << std::endl;
498  return rooplots;
499  }
500 
501  // If not provide, get data from workspace
502  if(!data) data = work->data("obsData");
503 
504  // Get category components
505  RooCategory* categories = work->cat("channelCat");
506 
507  // Print table with the number of data entries per category (channel)
508  data->table(*((RooAbsCategory*)categories))->Print("v");
509 
510  std::vector<TString> categoriesName;
511  TIterator* iter = categories->typeIterator();
512  RooCatType* catType;
513  while( (catType = (RooCatType*) iter->Next())) {
514  TString catname = catType->GetName();
515  categoriesName.push_back(catname);
516  }
517 
518  for(UInt_t i = 0; i < categoriesName.size(); i++){
519  TString catname = categoriesName[i];
520 
521  if(categories->setLabel(catname)){
522  std::cout << "WARNING::Category " << catname.Data() << " is not a member of channelCat. Will skip." << std::endl;
523  continue;
524  }
525 
526  RooAbsPdf* subpdf = (RooAbsPdf*)pdf->getPdf(catname.Data());
527  if(!subpdf){
528  std::cout << "WARNING::Can't find sub-pdf for region " << catname.Data() << ". Will skip." << std::endl;
529  continue;
530  }
531 
532  TString subdataset_str = Form("channelCat==channelCat::%s",catname.Data());
533  RooAbsData* subdataset = (RooAbsData*) data->reduce(subdataset_str.Data());
534  if(!subdataset){
535  std::cout << "WARNING::Can't find sub-dataset for region " << catname.Data() << ". Will skip." << std::endl;
536  continue;
537  }
538 
539  RooRealVar* var =(RooRealVar*) ((RooArgSet*) subpdf->getObservables(*subdataset))->find(Form("obs_x_%s", catname.Data()));
540 
541  // Define canvas first
542  TString canName = Form("canvas_%s_%s_%s", catname.Data(), name.Data(), plottodraw.Data());
543  TCanvas* c = new TCanvas(canName, canName);
544 
545  // Legend
546  TLegend* legend = new TLegend(0.10,0.70,0.90,0.90,"");
547  legend->SetFillStyle(0);
548  legend->SetFillColor(0);
549  legend->SetBorderSize(0);
550  legend->SetTextSize(0.036);
551 
552  legend->SetNColumns(3);
553 
554  TLegendEntry* entry = legend->AddEntry("","Data","pl") ;
555  entry->SetMarkerColor(kBlack);
556  entry->SetMarkerStyle(20);
557  entry=legend->AddEntry("","MC Total","lf") ;
558  entry->SetLineColor(kBlack);
559  entry->SetFillColor(kBlue);
560  entry->SetFillStyle(3444);
561 
562  RooPlot* frame = var->frame();
563  frame->SetName(Form("frame_%s_%s_%s", catname.Data(), name.Data(), plottodraw.Data()));
564 
565  // Draw data and pdf on top of error band - change to RooAbsData::SumW2 if data is weighted
566  if(error.Contains("SumW2") || error.Contains("sumw2") || error.Contains("sumW2"))
567  subdataset->plotOn(frame, RooFit::DataError(RooAbsData::SumW2), RooFit::MarkerColor(1), RooFit::LineColor(1));
568  else
569  subdataset->plotOn(frame, RooFit::DataError(RooAbsData::Poisson), RooFit::MarkerColor(1), RooFit::LineColor(1));
570 
571  TString RRSumPdfName = Form("%s_model",catname.Data());
572  RooRealSumPdf* RRSumPdf = (RooRealSumPdf*) subpdf->getComponents()->find(RRSumPdfName);
573 
574  RooArgList RRSumComponentsList = RRSumPdf->funcList();
575 
576  TString binWidthName = Form("binWidth_obs_x_%s_0",catname.Data());
577  RooRealVar* regionBinWidth = ((RooRealVar*) RRSumPdf->getVariables()->find(Form("binWidth_obs_x_%s_0",catname.Data()))) ;
578  if(regionBinWidth == NULL)
579  std::cout << "WARNING::bindWidth variable not found for region(" << catname << "). PLOTTING COMPONENTS WILL BE WRONG!" << std::endl;
580 
581  // Normalize data to expected number of events
582  double normCount = subpdf->expectedEvents(*var);
583  if(result)
584  std::cout << "INFO::MC events after fit = " << normCount << " for " << subpdf->GetName() << std::endl;
585 
586  RooLinkedListIter iter = RRSumComponentsList.iterator() ;
587  RooProduct* component;
588 
589  std::vector<TString> compNameVec;
590  std::vector <double> compFracVec;
591  std::vector<TString> compStackNameVec;
592  std::vector <double> compStackFracVec;
593  compNameVec.clear();
594  compStackNameVec.clear();
595  compFracVec.clear();
596  compStackFracVec.clear();
597 
598  while( (component = (RooProduct*) iter.Next()) ){
599  TString componentName = component->GetName();
600  TString stackComponentName = componentName;
601 
602  if(!compStackNameVec.empty())
603  stackComponentName = Form("%s,%s",compStackNameVec.back().Data() ,componentName.Data());
604  compNameVec.push_back(componentName);
605  compStackNameVec.push_back(stackComponentName);
606 
607  RooAbsReal* i_RRSumPdf = ((RooAbsPdf*)work->pdf(RRSumPdfName))->createIntegral(RooArgSet(*var));
608  Double_t Int_RRSumPdf = i_RRSumPdf->getVal();
609  RooAbsReal* i_component = ((RooProduct*)work->obj(componentName))->createIntegral(RooArgSet(*var));
610  Double_t Int_component = i_component->getVal();
611 
612  Double_t componentFrac = 0.;
613  if(Int_RRSumPdf != 0.)
614  componentFrac = Int_component * regionBinWidth->getVal() / Int_RRSumPdf;
615 
616  double stackComponentFrac = componentFrac;
617  if(!compStackFracVec.empty())
618  stackComponentFrac = compStackFracVec.back() + componentFrac;
619 
620  compFracVec.push_back(componentFrac);
621  compStackFracVec.push_back(stackComponentFrac);
622  }
623 
624  Int_t counter = 0;
625  Int_t sigcolor[13] = {2,3,4,5,6,7,8,9,kMagenta, 1, kGreen+2, kTeal, kOrange+10};
626  for(int i = (compFracVec.size()-1); i > -1; i--){
627  Int_t compPlotColor = i;
628  if(compNameVec[i].Contains("ChannelABS_CEX") || compNameVec[i].Contains("ChannelCEX_ABS")){
629  compPlotColor = 0;
630  }
631  else if(compNameVec[i].Contains("RecoBin")){
632  compPlotColor = 40;
633  }
634  else{
635  compPlotColor = sigcolor[counter];
636  counter++;
637  }
638  std::cout << "INFO::Drawing " << compNameVec[i] << " " << counter << " " << compPlotColor << std::endl;
639 
640  subpdf->plotOn(frame,RooFit::Components(compStackNameVec[i].Data()),RooFit::FillColor(compPlotColor),RooFit::FillStyle(3001),RooFit::DrawOption("F"),RooFit::Normalization(compStackFracVec[i]*normCount,RooAbsReal::NumEvent),RooFit::Precision(1e-5));
641  }
642 
643  if(result)
644  subpdf->plotOn(frame, RooFit::Normalization(1,RooAbsReal::RelativeExpected), RooFit::Precision(1e-5), RooFit::VisualizeError(*result), RooFit::FillColor(kBlue), RooFit::FillStyle(3444));
645 
646  // Plot again so that it is on top of errors
647  subpdf->plotOn(frame, RooFit::Normalization(1,RooAbsReal::RelativeExpected), RooFit::Precision(1e-5), RooFit::LineColor(1));
648  if(error.Contains("SumW2") || error.Contains("sumw2") || error.Contains("sumW2"))
649  subdataset->plotOn(frame, RooFit::DataError(RooAbsData::SumW2), RooFit::MarkerColor(1), RooFit::LineColor(1));
650  else
651  subdataset->plotOn(frame, RooFit::DataError(RooAbsData::Poisson), RooFit::MarkerColor(1), RooFit::LineColor(1));
652 
653  // Remove empty data bins
654  RemoveEmptyBins(frame);
655 
656  Int_t counter2 = 0;
657  bool found = false; bool found2 = false;
658  //for(int i = (compNameVec.size()-1) ; i > -1; i--)
659  for(unsigned int i = 0; i < compFracVec.size(); i++){
660  std::cout << "NAME: " << compNameVec[i] << std::endl;
661  Int_t compPlotColor = i;
662  if(compNameVec[i].Contains("ChannelABS_CEX") && !found){
663  compPlotColor = 0;
664  entry=legend->AddEntry("","CEX Bkg","f");
665  entry->SetLineColor(46);
666  entry->SetFillColor(compPlotColor);
667  entry->SetFillStyle(3001);
668  found = true;
669  continue;
670  }
671  else if(compNameVec[i].Contains("ChannelCEX_ABS") && !found){
672  compPlotColor = 0;
673  entry=legend->AddEntry("","ABS Bkg","f");
674  entry->SetLineColor(46);
675  entry->SetFillColor(compPlotColor);
676  entry->SetFillStyle(3001);
677  found = true;
678  continue;
679  }
680  else if(compNameVec[i].Contains("RecoBin") && !found2){
681  compPlotColor = 40;
682  entry=legend->AddEntry("","Signal","f");
683  entry->SetLineColor(40);
684  entry->SetFillColor(compPlotColor);
685  entry->SetFillStyle(3001);
686  found2 = true;
687  continue;
688  }
689  else{
690  counter--;
691  compPlotColor = sigcolor[counter];
692  }
693 
694  if(counter2 < (int)binnames.size()){
695  TString legName = binnames[counter2];
696  if (compNameVec[i].Contains("Incident"))
697  legName = incidentBinNames[counter2];
698  if (compNameVec[i].Contains("Sideband"))
699  legName = sidebandBinNames[counter2];
700 
701  entry=legend->AddEntry("",legName.Data(),"f");
702  entry->SetLineColor(compPlotColor);
703  entry->SetFillColor(compPlotColor);
704  entry->SetFillStyle(3001);
705  counter2++;
706  }
707  }
708 
709  // two pads, one for standard plot, one for data/MC ratio
710  float yMinP1=0.305;
711  float bottomMarginP1=0.015;
712  TPad *pad1 = new TPad(Form("%s_pad1",canName.Data()),Form("%s_pad1",canName.Data()),0.,yMinP1,.99,1);
713  pad1->SetBottomMargin(bottomMarginP1);
714  pad1->SetFillColor(kWhite);
715  pad1->SetTickx();
716  pad1->SetTicky();
717  TPad *pad2 = new TPad(Form("%s_pad2",canName.Data()),Form("%s_pad2",canName.Data()),0.,0.01,.99,0.295);
718  pad2->SetTopMargin(0.021);
719  pad2->SetBottomMargin(0.3);
720  pad2->SetFillColor(kWhite);
721 
722  pad1->Draw();
723  pad2->Draw();
724  frame->GetXaxis()->SetLabelSize(0.);
725 
726  pad1->cd();
727  frame->SetTitle(measurement.Data());
728  frame->Draw(); //Draw("same");
729  frame->GetXaxis()->SetTitle("Reco Bin");
730  frame->GetYaxis()->SetTitle("Events");
731 
732 
733  std::string chi2_text = "#chi^{2}: " +
734  std::to_string(GetDataMCChi2(work, catname, data));
735  legend->AddEntry((TObject*)0x0, chi2_text.c_str(), "");
736 
737  legend->Draw();
738 
739  pad2->cd();
740  RooPlot* frame_dummy = var->frame();
741 
742  subdataset->plotOn(frame_dummy,RooFit::Cut(subdataset_str),RooFit::DataError(RooAbsData::Poisson));
743 
744  // Get RooHist of the data - will be used later in the ratio plot
745  const char* curvename = 0;
746  RooHist* nominal_hist = (RooHist*) frame_dummy->findObject(curvename, RooHist::Class());
747  if(!nominal_hist){
748  std::cerr << "ERROR::Drawing the data/MC histogram. Can't find nominal histogram." << std::endl;
749  return rooplots;
750  }
751 
752  // Normalize pdf to number of expected events, not to number of events in dataset
753  subpdf->plotOn(frame_dummy,RooFit::Normalization(1,RooAbsReal::RelativeExpected),RooFit::Precision(1e-5));
754 
755  // Get RooCurve of the pdf - will be used later in the ratio plot
756  RooCurve* nominal_curve = (RooCurve*) frame_dummy->findObject(curvename, RooCurve::Class());
757  if(!nominal_curve){
758  std::cerr << "ERROR::Drawing the data/MC histogram. Can't find nominal curve." << std::endl;
759  return rooplots;
760  }
761 
762  // frame_dummy->Print("v");
763 
764  RooHist* hratio = NULL;
765  RooCurve* ratio_curve = new RooCurve;
766 
767  if(plottodraw == "pull"){
768  hratio = (RooHist*) frame_dummy->pullHist();
769  hratio->SetTitle("Pull Distribution");
770  }
771  else if(plottodraw == "residuals"){
772  hratio = (RooHist*) frame_dummy->residHist();
773  hratio->SetTitle("Residual Distribution");
774  }
775  else if(plottodraw == "ratio"){
776  if(result)
777  subpdf->plotOn(frame_dummy, RooFit::Normalization(1,RooAbsReal::RelativeExpected), RooFit::Precision(1e-5), RooFit::VisualizeError(*result, 1), RooFit::FillColor(kBlue-5), RooFit::FillStyle(3004));
778 
779  // Get error plot
780  RooCurve* error_curve = (RooCurve*)frame_dummy->findObject(curvename, RooCurve::Class());
781  if(!error_curve){
782  std::cerr << "ERROR::Drawing the data/MC histogram. Can't find error curve." << std::endl;
783  return rooplots;
784  }
785 
786  ratio_curve->SetName(Form("%s_ratiocurve",nominal_curve->GetName()));
787  ratio_curve->SetLineColor(kBlue-5);
788  ratio_curve->SetFillColor(kBlue-5);
789  ratio_curve->SetFillStyle(3004);
790  ratio_curve->SetLineWidth(1);
791 
792  // Fill error curve
793  Int_t j = 0;
794  bool bottomCurve = false;
795  for(Int_t i=1; i < error_curve->GetN()-1; i++){
796  Double_t x = 0.;
797  Double_t y = 0.;
798  error_curve->GetPoint(i,x,y) ;
799 
800  if( i >= (nominal_curve->GetN()-1) ) bottomCurve = true;
801 
802  Double_t xNom = x;
803  Double_t yNom = y;
804 
805  if( i == (nominal_curve->GetN() - 1) || i == nominal_curve->GetN() ){
806  ratio_curve->addPoint(x, 0.);
807  continue;
808  }
809 
810  if(bottomCurve){
811  nominal_curve->GetPoint(j,xNom,yNom);
812  j--;
813  }
814  else{
815  j++;
816  nominal_curve->GetPoint(j,xNom,yNom);
817  }
818 
819  if(fabs(yNom) > 0.00001){
820  ratio_curve->addPoint(x, (y / yNom));
821  }
822  else{
823  ratio_curve->addPoint(x, 0.);
824  }
825  }
826 
827  // Define ratio plot
828  hratio = new RooHist(nominal_hist->getNominalBinWidth());
829 
830  // Determine range of curve
831  Double_t xstart, xstop, y;
832  nominal_curve->GetPoint(2,xstart,y);
833  nominal_curve->GetPoint(nominal_curve->GetN()-1,xstop,y);
834 
835  for(Int_t i=0; i < nominal_hist->GetN(); i++){
836  Double_t x,ypoint;
837  nominal_hist->GetPoint(i,x,ypoint);
838 
839  if(x < xstart || x > xstop) continue;
840  if(fabs(ypoint) < 0.000001 ) continue;
841 
842  Double_t yy;
843  Double_t yerrorl = nominal_hist->GetErrorYlow(i);
844  Double_t yerrorh = nominal_hist->GetErrorYhigh(i);
845 
846  //Double_t xerrorl = nominal_curve->GetErrorXlow(i);
847  //Double_t xerrorh = nominal_curve->GetErrorXhigh(i);
848  //if(xerrorl<=0 ) xerrorl = nominal_curve->GetErrorX(i);
849  //if(xerrorh<=0 ) xerrorh = nominal_curve->GetErrorX(i);
850  //if(xerrorl<=0 ) xerrorl = 0.5*nominal_hist->getNominalBinWidth();
851  //if(yerrorh<=0 ) yerrorh = 0.5*nominal_hist->getNominalBinWidth();
852 
853  //yy = ypoint / nominal_curve->average(x-exl,x+exh);
854  //yerrorl /= nominal_curve->average(x-xerrorl,x+xerrorh);
855  //yerrorh /= nominal_curve->average(x-xerrorl,x+xerrorh);
856 
857  yy = ypoint / nominal_curve->interpolate(x);
858  yerrorl /= nominal_curve->interpolate(x);
859  yerrorh /= nominal_curve->interpolate(x);
860 
861  hratio->addBinWithError(x,yy,yerrorl,yerrorh);
862  }
863  }
864 
865  if(!hratio){
866  std::cerr << "ERROR::Drawing data and MC histogram failed. RooHist is not found." << std::endl;
867  return rooplots;
868  }
869 
870  hratio->SetMarkerColor(kRed);
871  hratio->SetLineColor(kRed);
872 
873  // Create a new frame to draw the residual/pull/ratio distribution and add the distribution to the frame
874  RooPlot* frame2 = var->frame();
875  if(plottodraw == "ratio" && result) frame2->addPlotable(ratio_curve,"F");
876  frame2->addPlotable(hratio,"P");
877 
878  if (catname.Contains("Sideband")) {
879  //frame2->GetXaxis()->SetBinLabel(1, "Generic Label");
880  frame2->GetXaxis()->SetTitle("Generic Label");
881  for(size_t i = 1; i < sidebandBins.size(); i++){
882  TString ibinstr = Form("%.1f-%.1f",sidebandBins[i-1],sidebandBins[i]);
883  frame2->GetXaxis()->SetBinLabel(i, ibinstr.Data());
884  }
885  }
886  else {
887  if (doNegativeReco) {
888  frame2->GetXaxis()->SetBinLabel(1, "< 0.");
889  for(size_t i = 1; i < recobins.size(); i++){
890  TString ibinstr = Form("%.1f-%.1f",recobins[i-1],recobins[i]);
891  frame2->GetXaxis()->SetBinLabel(i+1, ibinstr.Data());
892  }
893  }
894  else {
895  for(size_t i = 1; i < recobins.size(); i++){
896  TString ibinstr = Form("%.1f-%.1f",recobins[i-1],recobins[i]);
897  frame2->GetXaxis()->SetBinLabel(i, ibinstr.Data());
898  }
899  }
900  frame2->GetXaxis()->SetTitle("E_{reco} [MeV]");
901  }
902 
903  // Cosmetics
904  int firstbin = frame_dummy->GetXaxis()->GetFirst();
905  int lastbin = frame_dummy->GetXaxis()->GetLast();
906  double xmax = frame_dummy->GetXaxis()->GetBinUpEdge(lastbin);
907  double xmin = frame_dummy->GetXaxis()->GetBinLowEdge(firstbin);
908 
909  if(plottodraw=="pull"){
910  TLine* lp1 = new TLine(xmin,1.,xmax,1.);
911  TLine* lp2 = new TLine(xmin,2.,xmax,2.);
912  TLine* lp3 = new TLine(xmin,3.,xmax,3.);
913  TLine* lp4 = new TLine(xmin,4.,xmax,4.);
914 
915  TLine* lp6 = new TLine(xmin,-1.,xmax,-1.);
916  TLine* lp7 = new TLine(xmin,-2.,xmax,-2.);
917  TLine* lp8 = new TLine(xmin,-3.,xmax,-3.);
918  TLine* lp9 = new TLine(xmin,-4.,xmax,-4.);
919 
920  TLine* lp10 = new TLine(xmin,0,xmax,0);
921 
922  lp1->SetLineStyle(3);
923  lp2->SetLineStyle(3);
924  lp3->SetLineStyle(3);
925  lp4->SetLineStyle(3);
926  lp6->SetLineStyle(3);
927  lp7->SetLineStyle(3);
928  lp8->SetLineStyle(3);
929  lp9->SetLineStyle(3);
930  lp10->SetLineStyle(3);
931 
932  frame2->addObject(lp1);
933  frame2->addObject(lp2);
934  frame2->addObject(lp3);
935  frame2->addObject(lp4);
936  frame2->addObject(lp6);
937  frame2->addObject(lp7);
938  frame2->addObject(lp8);
939  frame2->addObject(lp9);
940  frame2->addObject(lp10);
941 
942  frame2->SetMinimum(-4.0);
943  frame2->SetMaximum(4.0);
944 
945  frame2->GetYaxis()->SetTitle("Pull");
946  }
947  else if(plottodraw == "residuals"){
948  TLine* l = new TLine(xmin,0.,xmax,0.);
949  l->SetLineWidth(1);
950  l->SetLineStyle(2);
951  frame2->addObject(l);
952  frame2->GetYaxis()->SetTitle("Residual");
953  }
954  else if(plottodraw == "ratio"){
955  TLine* lp1 = new TLine(xmin,1.,xmax,1.);
956  TLine* lp2 = new TLine(xmin,0.5,xmax,0.5);
957  TLine* lp3 = new TLine(xmin,1.5,xmax,1.5);
958  TLine* lp4 = new TLine(xmin,2.,xmax,2.);
959  TLine* lp5 = new TLine(xmin,2.5,xmax,2.5);
960  lp1->SetLineWidth(1);
961  lp1->SetLineStyle(3);
962  lp2->SetLineStyle(3);
963  lp3->SetLineStyle(3);
964  lp4->SetLineStyle(3);
965  lp5->SetLineStyle(3);
966  frame2->addObject(lp1);
967  frame2->addObject(lp2);
968  frame2->addObject(lp3);
969  frame2->addObject(lp4);
970  frame2->addObject(lp5);
971 
972  frame2->SetMinimum(.0);
973  frame2->SetMaximum(3.0);
974 
975  frame2->GetYaxis()->SetTitle("Data / MC");
976  }
977 
978  frame2->GetYaxis()->SetLabelSize(0.10);
979  frame2->GetYaxis()->SetNdivisions(504);
980  frame2->GetXaxis()->SetLabelSize(0.10);
981  frame2->GetYaxis()->SetTitleSize(0.10);
982  frame2->GetXaxis()->SetTitleSize(0.10);
983  frame2->GetYaxis()->SetTitleOffset(0.35);
984  frame2->GetXaxis()->SetTitleOffset(1.);
985  frame2->GetYaxis()->SetLabelOffset(0.01);
986  frame2->GetXaxis()->SetLabelOffset(0.03);
987  frame2->GetXaxis()->SetTickLength(0.06);
988 
989  frame2->SetTitle("");
990  frame2->GetYaxis()->CenterTitle();
991  frame2->Draw();
992 
993  rooplots.push_back(c);
994  }
995 
996  return rooplots;
997 
998 }
999 
1000 //********************************************************************
1001 std::vector<TCanvas*> protoana::ProtoDUNEFitUtils::PlotNLL(RooWorkspace *work, TString name, RooFitResult* result, bool plotPLL){
1002  //********************************************************************
1003 
1004  std::vector<TCanvas*> rooplots;
1005 
1006  if(!work){
1007  std::cerr << "ERROR::NULL workspace. Will return empty vector!" << std::endl;
1008  return rooplots;
1009  }
1010 
1011  if(!result){
1012  std::cerr << "ERROR::NULL fit result. Will return empty vector!" << std::endl;
1013  return rooplots;
1014  }
1015 
1016  RooSimultaneous* pdf = (RooSimultaneous*) work->pdf("simPdf");
1017  if(!pdf){
1018  std::cout << "ERROR::No pdf found in workspace. Will return empty vector!" << std::endl;
1019  return rooplots;
1020  }
1021 
1022  // Get data from workspace
1023  RooAbsData* data = work->data("obsData");
1024 
1025  // Get category components
1026  RooCategory* categories = work->cat("channelCat");
1027  // Print table with the number of data entries per category (channel)
1028  data->table(*((RooAbsCategory*)categories))->Print("v");
1029 
1030  RooStats::ModelConfig *combined_config = (RooStats::ModelConfig*)work->obj("ModelConfig");
1031  if(!combined_config){
1032  std::cout << "ERROR::No model config " << " ModelConfig " << " in workspace. Will return empty vector!" << std::endl;
1033  return rooplots;
1034  }
1035 
1036  const RooArgSet* obs = combined_config->GetGlobalObservables();
1037  if(!obs){
1038  std::cout << "ERROR::No observables found in ModelConfig. Will return empty vector!" << std::endl;
1039  return rooplots;
1040  }
1041 
1042  RooArgList floatParsFinal = result->floatParsFinal();
1043 
1044  // Create NLL
1045  RooAbsReal* nll = pdf->createNLL(*data, RooFit::NumCPU(2), RooFit::GlobalObservables(*obs), RooFit::Offset(true));
1046 
1047  for(Int_t i = 0; i < floatParsFinal.getSize(); i++){
1048  RooAbsArg* arg = floatParsFinal.at(i);
1049  if(!arg->InheritsFrom("RooRealVar")) continue;
1050 
1051  RooRealVar* par = (RooRealVar*) arg;
1052  TString parName = par->GetName();
1053 
1054  // set parameter range to readable range
1055  double minRange = par->getMin();
1056  double maxRange = par->getMax();
1057  if(minRange < 0.){
1058  par->setMin(-3.);
1059  par->setMax(3.);
1060  }
1061  else{
1062  par->setMin(minRange);
1063  par->setMax(2.);
1064  }
1065 
1066  RooPlot* frame = par->frame();
1067  nll->plotOn(frame, RooFit::ShiftToZero());
1068  frame->SetMinimum(0.);
1069  // To be able to see the 1/2 sigma
1070  frame->SetMaximum(2.5);
1071 
1072  RooAbsReal* pll = NULL;
1073  if(plotPLL){
1074  pll = nll->createProfile(*par) ;
1075  pll->plotOn(frame, RooFit::LineColor(kRed), RooFit::LineStyle(kDashed), RooFit::NumCPU(4));
1076  }
1077 
1078  TString canName=Form("Canvas_NLL_%s_%s", name.Data(), parName.Data());
1079  TCanvas* c = new TCanvas(canName,canName,600,600);
1080  c->cd();
1081  frame->Draw();
1082 
1083  TLegend* legend = new TLegend(0.55,0.65,0.85,0.95,"");
1084  legend->SetFillStyle(0);
1085  legend->SetFillColor(0);
1086  legend->SetBorderSize(0);
1087  legend->SetTextSize(0.036);
1088  TLegendEntry* entry=legend->AddEntry("","NLL","l") ;
1089  entry->SetLineColor(kBlue);
1090  if(plotPLL){
1091  entry=legend->AddEntry("","PLL","l") ;
1092  entry->SetLineColor(kRed);
1093  entry->SetLineStyle(kDashed);
1094  }
1095  legend->Draw();
1096 
1097  // reset parameter range to previous values
1098  par->setMin(minRange);
1099  par->setMax(maxRange);
1100 
1101  if(pll) delete pll;
1102 
1103  rooplots.push_back(c);
1104  }
1105 
1106  return rooplots;
1107 
1108 }
1109 
1110 //********************************************************************
1111 std::vector<TCanvas*> protoana::ProtoDUNEFitUtils::PlotNuisanceParametersImpact(RooWorkspace *work, RooFitResult* fitresult, TString snapshotload, std::vector<std::string> SystsToConsider){
1112  //********************************************************************
1113 
1114  std::vector<TCanvas*> plots;
1115 
1116  if(!work){
1117  std::cerr << "ERROR::NULL workspace. Will return empty vector!" << std::endl;
1118  return plots;
1119  }
1120 
1121  // Get the configuration model from workspace
1122  RooStats::ModelConfig *combined_config = (RooStats::ModelConfig*)work->obj("ModelConfig");
1123  if(!combined_config){
1124  std::cerr << "ERROR::No model config " << "ModelConfig " << " in workspace. Return empty vector!" << std::endl;
1125  return plots;
1126  }
1127 
1128  // Dataset
1129  RooAbsData *obsdata = work->data("obsData");
1130 
1131  if(!fitresult){
1132  std::cerr << "ERROR::NULL fit result. Will return empty vector!" << std::endl;
1133  return plots;
1134  }
1135 
1136  protoana::ProtoDUNEFitUtils::LoadSnapshot(work, snapshotload.Data());
1137  protoana::ProtoDUNEFitUtils::MakeNuisanceParamsConstant(work,"");
1138 
1139  RooArgList floatParList = protoana::ProtoDUNEFitUtils::GetPostfitPOIList(fitresult->floatParsFinal());
1140  //floatParList.Print("v");
1141 
1142  const int nSystToConsider = SystsToConsider.size();
1143  const int nPOI = floatParList.getSize();
1144 
1145  std::vector<TH1D*> systhistovec; std::vector<TH1D*> systhistovec_prefitup; std::vector<TH1D*> systhistovec_prefitdown; std::vector<TH1D*> systhistovec_postfitup; std::vector<TH1D*> systhistovec_postfitdown;
1146 
1147  for(int i=0; i < nPOI; i++){
1148  TH1D* systhisto = new TH1D(Form("systhisto%i",i), Form("systhisto%i",i), nSystToConsider, 0, nSystToConsider);
1149  TH1D* systhisto_prefitup = new TH1D(Form("systhisto%i_prefitup",i), Form("systhisto%i_prefitup",i), nSystToConsider, 0, nSystToConsider);
1150  TH1D* systhisto_prefitdown = new TH1D(Form("systhisto%i_prefitdown",i), Form("systhisto%i_prefitdown",i), nSystToConsider, 0, nSystToConsider);
1151  TH1D* systhisto_postfitup = new TH1D(Form("systhisto%i_postfitup",i), Form("systhisto%i_postfitup",i), nSystToConsider, 0, nSystToConsider);
1152  TH1D* systhisto_postfitdown = new TH1D(Form("systhisto%i_postfitdown",i),Form("systhisto%i_postfitdown",i), nSystToConsider, 0, nSystToConsider);
1153 
1154  systhistovec.push_back(systhisto);
1155  systhistovec_prefitup.push_back(systhisto_prefitup);
1156  systhistovec_prefitdown.push_back(systhisto_prefitdown);
1157  systhistovec_postfitup.push_back(systhisto_postfitup);
1158  systhistovec_postfitdown.push_back(systhisto_postfitdown);
1159  }
1160 
1161  for(int i=0; i < nSystToConsider; i++){
1162 
1163  TString systname = TString("alpha_") + TString(SystsToConsider[i].c_str());
1164  RooRealVar* var = work->var(systname.Data());
1165  if(!var){
1166  std::cerr << "ERROR::Variable " << systname.Data() << " not found in workspace. Skip!" << std::endl;
1167  continue;
1168  }
1169 
1170  double originalval = var->getVal();
1171  double originalerror = var->getError();
1172 
1173  for(int j=0; j < nPOI; j++){
1174  systhistovec[j]->SetBinContent(i+1, originalval);
1175  systhistovec[j]->SetBinError(i+1, originalerror);
1176  systhistovec[j]->GetXaxis()->SetBinLabel(i+1, SystsToConsider[i].c_str());
1177  systhistovec[j]->SetTitle( (floatParList.at(j))->GetName() );
1178  }
1179 
1180  // Using the pre-fit error
1181  var->setVal(originalval + 1.0);
1182  RooFitResult *tempfitresult = combined_config->GetPdf()->fitTo(*obsdata, RooFit::Minimizer("Minuit2"), RooFit::Save(), RooFit::Strategy(1), RooFit::Minos(false), RooFit::PrintLevel(1), RooFit::Constrain(*combined_config->GetPdf()->getParameters(*obsdata)), RooFit::Extended(), RooFit::Offset(true) );
1183  RooArgList paramsfit = protoana::ProtoDUNEFitUtils::GetPostfitPOIList(tempfitresult->floatParsFinal());
1184  //paramsfit.Print("v");
1185 
1186  var->setVal(originalval - 1.0);
1187  RooFitResult *tempfitresult2 = combined_config->GetPdf()->fitTo(*obsdata, RooFit::Minimizer("Minuit2"), RooFit::Save(), RooFit::Strategy(1), RooFit::Minos(false), RooFit::PrintLevel(1), RooFit::Constrain(*combined_config->GetPdf()->getParameters(*obsdata)), RooFit::Extended(), RooFit::Offset(true) );
1188  RooArgList paramsfit2 = protoana::ProtoDUNEFitUtils::GetPostfitPOIList(tempfitresult2->floatParsFinal());
1189 
1190  var->setVal(originalval + originalerror);
1191  RooFitResult *tempfitresult3 = combined_config->GetPdf()->fitTo(*obsdata, RooFit::Minimizer("Minuit2"), RooFit::Save(), RooFit::Strategy(1), RooFit::Minos(false), RooFit::PrintLevel(1), RooFit::Constrain(*combined_config->GetPdf()->getParameters(*obsdata)), RooFit::Extended(), RooFit::Offset(true) );
1192  RooArgList paramsfit3 = protoana::ProtoDUNEFitUtils::GetPostfitPOIList(tempfitresult3->floatParsFinal());
1193 
1194  var->setVal(originalval - originalerror);
1195  RooFitResult *tempfitresult4 = combined_config->GetPdf()->fitTo(*obsdata, RooFit::Minimizer("Minuit2"), RooFit::Save(), RooFit::Strategy(1), RooFit::Minos(false), RooFit::PrintLevel(1), RooFit::Constrain(*combined_config->GetPdf()->getParameters(*obsdata)), RooFit::Extended(), RooFit::Offset(true) );
1196  RooArgList paramsfit4 = protoana::ProtoDUNEFitUtils::GetPostfitPOIList(tempfitresult4->floatParsFinal());
1197 
1198  TIterator* Inititer = floatParList.createIterator();
1199  RooRealVar* Initvar(0);
1200  for(int j=0; (Initvar = (RooRealVar*)Inititer->Next()); j++){
1201  double DMu = 0;
1202  TString Initvarname(Initvar->GetName());
1203  if(Initvarname.Contains("alpha_") || Initvarname.Contains("gamma_")) continue;
1204  //Double_t Initvarerror = Initvar->getError();
1205  // Protection
1206  //if(Initvarerror == 0) Initvarerror = 1.0;
1207 
1208  TIterator* iter = paramsfit.createIterator();
1209  RooRealVar* Postvar(0);
1210  for(int k=0; (Postvar = (RooRealVar*)iter->Next()); k++){
1211  TString varname(Postvar->GetName());
1212  if(varname == Initvarname){
1213  DMu = Postvar->getVal() - Initvar->getVal();
1214  break;
1215  }
1216  }
1217  systhistovec_postfitup[j]->SetBinContent(i+1, DMu);
1218 
1219  DMu = 0;
1220  TIterator* iter2 = paramsfit2.createIterator();
1221  RooRealVar* Postvar2(0);
1222  for(int k=0; (Postvar2 = (RooRealVar*)iter2->Next()); k++){
1223  TString varname(Postvar2->GetName());
1224  if(varname == Initvarname){
1225  //DTheta = (Postvar->getVal() - Initvar->getVal())/Initvarerror;
1226  DMu = Postvar2->getVal() - Initvar->getVal();
1227  break;
1228  }
1229  }
1230  systhistovec_postfitdown[j]->SetBinContent(i+1, DMu);
1231 
1232  DMu = 0;
1233  TIterator* iter3 = paramsfit3.createIterator();
1234  RooRealVar* Postvar3(0);
1235  for(int k=0; (Postvar3 = (RooRealVar*)iter3->Next()); k++){
1236  TString varname(Postvar3->GetName());
1237  if(varname == Initvarname){
1238  DMu = Postvar3->getVal() - Initvar->getVal();
1239  break;
1240  }
1241  }
1242  systhistovec_prefitup[j]->SetBinContent(i+1, DMu);
1243 
1244  DMu = 0;
1245  TIterator* iter4 = paramsfit4.createIterator();
1246  RooRealVar* Postvar4(0);
1247  for(int k=0; (Postvar4 = (RooRealVar*)iter4->Next()); k++){
1248  TString varname(Postvar4->GetName());
1249  if(varname == Initvarname){
1250  DMu = Postvar4->getVal() - Initvar->getVal();
1251  break;
1252  }
1253  }
1254  systhistovec_prefitdown[j]->SetBinContent(i+1, DMu);
1255 
1256  delete iter;
1257  delete iter2;
1258  delete iter3;
1259  delete iter4;
1260  }
1261 
1262  delete Inititer;
1263 
1264  // Back to original values
1265  var->setVal(originalval);
1266  var->setError(originalerror);
1267 
1268  }
1269 
1270  TLine* lp1 = new TLine(0,-1.,nSystToConsider,-1.);
1271  lp1->SetLineStyle(kDashed);
1272  TLine* lp2 = new TLine(0,1.,nSystToConsider,1.);
1273  lp2->SetLineStyle(kDashed);
1274  TLine* lp3 = new TLine(0,0,nSystToConsider,0);
1275  lp3->SetLineColor(kBlack);
1276 
1277  TLegend* legend = new TLegend(0.20,0.80,0.85,0.90,"");
1278  legend->SetFillStyle(0);
1279  legend->SetFillColor(0);
1280  legend->SetBorderSize(0);
1281  legend->SetTextSize(0.036);
1282 
1283  legend->SetNColumns(3);
1284  TLegendEntry* entry = legend->AddEntry("","Pull","pl") ;
1285  entry->SetMarkerColor(kBlue);
1286  entry->SetMarkerStyle(20);
1287  entry=legend->AddEntry("","Pre-fit impact on #mu","f") ;
1288  entry->SetFillColor(kYellow);
1289  entry->SetLineColor(kYellow);
1290  entry->SetFillStyle(3001);
1291  entry=legend->AddEntry("","Post-fit impact on #mu","f") ;
1292  entry->SetFillColor(kBlue);
1293  entry->SetFillStyle(3005);
1294 
1295  TGaxis *Maxis = new TGaxis(nSystToConsider,-1.4, nSystToConsider, 1.4,-1.4,1.4,510,"+L");
1296  Maxis->SetLabelSize(0.03);
1297  Maxis->SetTextFont(72);
1298  Maxis->SetLabelOffset(0.025);
1299  Maxis->SetTitleOffset(1.1);
1300  Maxis->SetTitle("#Delta#mu");
1301 
1302  for(int i=0; i < nPOI; i++){
1303  systhistovec_prefitup[i]->SetLineColor(kBlue);
1304  systhistovec_prefitup[i]->SetFillColor(kBlue);
1305  systhistovec_prefitup[i]->SetFillStyle(3005);
1306  systhistovec_prefitdown[i]->SetLineColor(kBlue);
1307  systhistovec_prefitdown[i]->SetFillColor(kBlue);
1308  systhistovec_prefitdown[i]->SetFillStyle(3005);
1309  systhistovec_postfitup[i]->SetLineColor(kYellow);
1310  systhistovec_postfitup[i]->SetFillColor(kYellow);
1311  systhistovec_postfitdown[i]->SetLineColor(kYellow);
1312  systhistovec_postfitdown[i]->SetFillColor(kYellow);
1313  systhistovec[i]->SetMarkerColor(4);
1314  systhistovec[i]->SetMarkerStyle(20);
1315  systhistovec[i]->SetLineWidth(2);
1316  systhistovec[i]->SetStats(false);
1317  systhistovec[i]->GetYaxis()->SetRangeUser(-1.4,1.4);
1318  systhistovec[i]->GetYaxis()->SetTitle("(#theta - #theta_{0}) / #Delta#theta");
1319  systhistovec[i]->GetYaxis()->SetTitleOffset(1.15);
1320 
1321  TCanvas* can = new TCanvas(Form("can_impact%i",i),Form("can_impact%i",i));
1322  systhistovec[i]->Draw();
1323  systhistovec_postfitup[i]->Draw("same");
1324  systhistovec_postfitdown[i]->Draw("same");
1325  systhistovec_prefitup[i]->Draw("same");
1326  systhistovec_prefitdown[i]->Draw("same");
1327  systhistovec[i]->Draw("e1same");
1328  lp1->Draw("same");
1329  lp2->Draw("same");
1330  lp3->Draw("same");
1331  can->Update();
1332  gPad->RedrawAxis();
1333  Maxis->Draw("same");
1334  legend->Draw();
1335 
1336  plots.push_back(can);
1337  }
1338 
1339  return plots;
1340 
1341 }
1342 
1343 //********************************************************************
1344 TCanvas* protoana::ProtoDUNEFitUtils::PlotParametersPull(TTree* tree, RooWorkspace* ws){
1345  //********************************************************************
1346 
1347  if(!tree || !ws){
1348  std::cerr << "ERROR::No tree or workspace found. Pull plot failed!" << std::endl;
1349  return NULL;
1350  }
1351 
1352  if(tree->GetEntries() < 2){
1353  std::cerr << "ERROR::Tree has less than two entries. No pull plots! Ignore if you run on data!" << std::endl;
1354  return NULL;
1355  }
1356 
1357  RooStats::ModelConfig *combined_config = (RooStats::ModelConfig*)ws->obj("ModelConfig");
1358  if(!combined_config){
1359  std::cerr << "ERROR::No model config " << " ModelConfig " << " in workspace. No pull plots!" << std::endl;
1360  return NULL;
1361  }
1362 
1363  std::vector<Float_t> varValsPull;
1364  const RooArgSet* floatParsList = combined_config->GetParametersOfInterest();
1365  const Int_t n = floatParsList->getSize();
1366  if(n == 0){
1367  std::cerr << "ERROR::No floating parameters found. Pull plot failed!" << std::endl;
1368  return NULL;
1369  }
1370  varValsPull.resize( floatParsList->getSize(), -999. );
1371 
1372  TCanvas* cpull = new TCanvas("PullMeanSigma","PullMeanSigma");
1373 
1374  TH1F* pullmeanhisto = new TH1F("pullmeanplot", "", n+1, 0, n+1);
1375  TH1F* pullsigmahisto = new TH1F("pullsigmaplot", "", n+1, 0, n+1);
1376 
1377  Int_t status;
1378  tree->SetBranchAddress("status", &status);
1379 
1380  RooRealVar* var(0);
1381  TIterator* Itr = floatParsList->createIterator();
1382  Int_t counter = 0;
1383  for (Int_t i=0; (var = (RooRealVar*)Itr->Next()); ++i) {
1384  if(var->isConstant()) continue;
1385  TString varName = var->GetName() + TString("pull");
1386 
1387  if(varName.Contains("Lumi") || varName.Contains("binWidth") || varName.Contains("corr")) continue;
1388  tree->SetBranchAddress(varName.Data(), &varValsPull[i]);
1389 
1390  TH1F* pullhisto = new TH1F(Form("pullhisto%i",i),Form("pullhisto%i",i),100,-5,5);
1391  for(Int_t j=0; j < tree->GetEntries(); j++){
1392  tree->GetEntry(j);
1393  if(status%1000 == 0)
1394  pullhisto->Fill(varValsPull[i]);
1395  }
1396 
1397  if(pullhisto->GetEntries() == 0){
1398  //delete pullhisto;
1399  continue;
1400  }
1401 
1402  pullhisto->Fit("gaus","Q");
1403  TF1* fit = pullhisto->GetFunction("gaus");
1404  if(!fit){
1405  //delete pullhisto;
1406  continue;
1407  }
1408 
1409  pullmeanhisto->SetBinContent(counter+1, fit->GetParameter(1));
1410  pullmeanhisto->SetBinError(counter+1, fit->GetParError(1));
1411  pullmeanhisto->GetXaxis()->SetBinLabel(counter+1, varName.Data());
1412 
1413  pullsigmahisto->SetBinContent(counter+1, fit->GetParameter(2));
1414  pullsigmahisto->SetBinError(counter+1, fit->GetParError(2));
1415  pullsigmahisto->GetXaxis()->SetBinLabel(counter+1, varName.Data());
1416 
1417  counter++;
1418 
1419  //delete pullhisto;
1420  }
1421 
1422  delete Itr;
1423 
1424  // Decorate histograms
1425  pullmeanhisto->SetLineColor(2);
1426  pullmeanhisto->SetMarkerStyle(20);
1427  pullmeanhisto->SetMarkerColor(2);
1428 
1429  pullsigmahisto->SetLineColor(1);
1430  pullsigmahisto->SetMarkerStyle(21);
1431  pullsigmahisto->SetMarkerColor(1);
1432 
1433  pullmeanhisto->GetYaxis()->SetRangeUser(-3, 3);
1434  pullsigmahisto->GetYaxis()->SetRangeUser(-3, 3);
1435 
1436  pullmeanhisto->GetXaxis()->SetLabelSize(0.015);
1437  pullsigmahisto->GetXaxis()->SetLabelSize(0.015);
1438 
1439  pullmeanhisto->SetStats(0);
1440  pullsigmahisto->SetStats(0);
1441 
1442  TLegend* legend = new TLegend(0.55,0.65,0.85,0.95,"");
1443  legend->SetFillStyle(0);
1444  legend->SetFillColor(0);
1445  legend->SetBorderSize(0);
1446  legend->SetTextSize(0.036);
1447  legend->AddEntry(pullmeanhisto, "Pull Mean", "l");
1448  legend->AddEntry(pullsigmahisto,"Pull Sigma", "l");
1449 
1450  TLine *mline = new TLine(0,0.0,n,0.0);
1451  mline->SetLineColor(kBlue);
1452  TLine *sline = new TLine(0,1.0,n,1.0);
1453  sline->SetLineColor(kBlue);
1454 
1455  cpull->cd();
1456  pullsigmahisto->Draw("e");
1457  pullmeanhisto->Draw("esame");
1458  mline->Draw();
1459  sline->Draw();
1460  legend->Draw();
1461 
1462  return cpull;
1463 
1464 }
1465 
1466 //********************************************************************
1467 TCanvas* protoana::ProtoDUNEFitUtils::PlotNuisanceParametersPull(TTree* tree, RooWorkspace* ws){
1468  //********************************************************************
1469 
1470  if(!tree || !ws){
1471  std::cerr << "ERROR::No tree or workspace found. Pull plot failed!" << std::endl;
1472  return NULL;
1473  }
1474 
1475  if(tree->GetEntries() < 2){
1476  std::cerr << "ERROR::Tree has less than two entries. No pull plots! Ignore if you run on data!" << std::endl;
1477  return NULL;
1478  }
1479 
1480  RooStats::ModelConfig *combined_config = (RooStats::ModelConfig*)ws->obj("ModelConfig");
1481  if(!combined_config){
1482  std::cerr << "ERROR::No model config " << " ModelConfig " << " in workspace. No pull plots!" << std::endl;
1483  return NULL;
1484  }
1485 
1486  std::vector<Float_t> varValsPull;
1487  const RooArgSet* nuisanceParsList = combined_config->GetNuisanceParameters();
1488  if(!nuisanceParsList){
1489  std::cerr << "ERROR::No nuisance parameters found. Nuisance parameters plot failed!" << std::endl;
1490  return NULL;
1491  }
1492  const Int_t n = nuisanceParsList->getSize();
1493  varValsPull.resize( nuisanceParsList->getSize(), -999. );
1494 
1495  TCanvas* cpullnui = new TCanvas("NuisancePullMeanSigma","NuisancePullMeanSigma");
1496 
1497  TH1F* nuispullmeanhisto = new TH1F("nuispullmeanplot", "", n+1, 0, n+1);
1498  TH1F* nuispullsigmahisto = new TH1F("nuispullsigmaplot", "", n+1, 0, n+1);
1499 
1500  Int_t status;
1501  tree->SetBranchAddress("status", &status);
1502 
1503  RooRealVar* var(0);
1504  TIterator* Itr = nuisanceParsList->createIterator();
1505  Int_t counter = 0;
1506  int ngammas = 0;
1507  for (Int_t i=0; (var = (RooRealVar*)Itr->Next()); ++i) {
1508  if(var->isConstant()) continue;
1509  TString varName = var->GetName() + TString("pull");
1510 
1511  if(!varName.Contains("alpha") && !varName.Contains("gamma_stat")) continue;
1512  if(varName.Contains("Lumi") || varName.Contains("binWidth") || varName.Contains("corr")) continue;
1513  tree->SetBranchAddress(varName.Data(), &varValsPull[i]);
1514 
1515  if (varName.Contains("gamma_stat"))
1516  ++ngammas;
1517 
1518  TH1F* nuispullhisto = new TH1F(Form("nuispullhisto%i",i),Form("nuispullhisto%i",i),100,-5,5);
1519  for(Int_t j=0; j < tree->GetEntries(); j++){
1520  tree->GetEntry(j);
1521  if(status%1000 == 0)
1522  nuispullhisto->Fill(varValsPull[i]);
1523  }
1524 
1525  if(nuispullhisto->GetEntries() == 0){
1526  //delete pullhisto;
1527  continue;
1528  }
1529 
1530  nuispullhisto->Fit("gaus","Q");
1531  TF1* fit = nuispullhisto->GetFunction("gaus");
1532  if(!fit){
1533  //delete pullhisto;
1534  continue;
1535  }
1536 
1537  nuispullmeanhisto->SetBinContent(counter+1, fit->GetParameter(1));
1538  nuispullmeanhisto->SetBinError(counter+1, fit->GetParError(1));
1539  //nuispullmeanhisto->GetXaxis()->SetBinLabel(counter+1, varName.Data());
1540 
1541  nuispullsigmahisto->SetBinContent(counter+1, fit->GetParameter(2));
1542  nuispullsigmahisto->SetBinError(counter+1, fit->GetParError(2));
1543  //nuispullsigmahisto->GetXaxis()->SetBinLabel(counter+1, varName.Data());
1544 
1545  counter++;
1546 
1547  //delete pullhisto;
1548  }
1549 
1550  delete Itr;
1551 
1552  // Decorate histograms
1553  nuispullmeanhisto->SetLineColor(2);
1554  nuispullmeanhisto->SetMarkerStyle(20);
1555  nuispullmeanhisto->SetMarkerColor(2);
1556 
1557  nuispullsigmahisto->SetLineColor(1);
1558  nuispullsigmahisto->SetMarkerStyle(21);
1559  nuispullsigmahisto->SetMarkerColor(1);
1560 
1561  nuispullmeanhisto->GetYaxis()->SetRangeUser(-3, 3);
1562  nuispullsigmahisto->GetYaxis()->SetRangeUser(-3, 3);
1563 
1564  nuispullmeanhisto->SetStats(0);
1565  nuispullsigmahisto->SetStats(0);
1566 
1567  nuispullmeanhisto->GetXaxis()->SetTitle("Systematic ID");
1568  nuispullsigmahisto->GetXaxis()->SetTitle("Systematic ID");
1569 
1570  TLegend* legend = new TLegend(0.55,0.65,0.85,0.95,"");
1571  legend->SetFillStyle(0);
1572  legend->SetFillColor(0);
1573  legend->SetBorderSize(0);
1574  legend->SetTextSize(0.036);
1575  legend->AddEntry(nuispullmeanhisto, "Pull Mean", "l");
1576  legend->AddEntry(nuispullsigmahisto,"Pull Sigma", "l");
1577 
1578  TLine *mline = new TLine(0,0.0,n,0.0);
1579  mline->SetLineColor(kBlue);
1580  TLine *sline = new TLine(0,1.0,n,1.0);
1581  sline->SetLineColor(kBlue);
1582  TLine *vline = new TLine(n-ngammas, -3., n - ngammas, 3.);
1583  vline->SetLineColor(kGreen);
1584 
1585  cpullnui->cd();
1586  nuispullsigmahisto->Draw("e");
1587  nuispullmeanhisto->Draw("esame");
1588  mline->Draw();
1589  sline->Draw();
1590  legend->Draw();
1591  vline->Draw();
1592 
1593  return cpullnui;
1594 
1595 }
1596 
1597 //********************************************************************
1598 TCanvas* protoana::ProtoDUNEFitUtils::PlotNuisanceParameters(TTree* tree, RooWorkspace* ws){
1599  //********************************************************************
1600 
1601  if(!tree || !ws){
1602  std::cerr << "ERROR::No tree or workspace found. Pull plot failed!" << std::endl;
1603  return NULL;
1604  }
1605 
1606  RooStats::ModelConfig *combined_config = (RooStats::ModelConfig*)ws->obj("ModelConfig");
1607  if(!combined_config){
1608  std::cerr << "ERROR::No model config " << " ModelConfig " << " in workspace. No nuisance parameters plots!" << std::endl;
1609  return NULL;
1610  }
1611 
1612  Int_t status;
1613  tree->SetBranchAddress("status", &status);
1614 
1615  RooRealVar* var(0);
1616  std::vector<Float_t> varValsNuisance;
1617  const RooArgSet* nuisanceParsList = combined_config->GetNuisanceParameters();
1618  if(!nuisanceParsList){
1619  std::cerr << "ERROR::No nuisance parameters found. Nuisance parameters plot failed!" << std::endl;
1620  return NULL;
1621  }
1622  const Int_t n = nuisanceParsList->getSize();
1623  varValsNuisance.resize( nuisanceParsList->getSize(), -999. );
1624 
1625  TH1F* nuisancehisto = new TH1F("nuisanceplot", "", n+1, 0, n+1);
1626 
1627  TIterator* Itr = nuisanceParsList->createIterator();
1628  for (Int_t i=0; (var = (RooRealVar*)Itr->Next()); ++i) {
1629  TString varName = var->GetName();
1630  if(!varName.Contains("alpha") && !varName.Contains("gamma_stat")) continue;
1631  if(varName.Contains("Lumi") || varName.Contains("nom") || varName.Contains("binWidth") || varName.Contains("err") || varName.Contains("pull")) continue;
1632 
1633  tree->SetBranchAddress(varName.Data(), &varValsNuisance[i]);
1634 
1635  //std::cout << "INFO::Plotting nuisance parameter " << varName << std::endl;
1636  TH1F* temphisto = new TH1F("temphisto", "temphisto", 100, -5, 5);
1637 
1638  for(Int_t j=0; j < tree->GetEntries(); j++){
1639  tree->GetEntry(j);
1640  if(status%1000 == 0)
1641  temphisto->Fill(varValsNuisance[i]);
1642  //nuisancehisto->SetBinContent(i+1, varValsNuisance[i]);
1643  }
1644  nuisancehisto->SetBinContent(i+1, temphisto->GetMean());
1645  delete temphisto;
1646  }
1647 
1648  delete Itr;
1649 
1650  Int_t ngammas = 0;
1651  // Fill the error
1652  TIterator* Itr2 = nuisanceParsList->createIterator();
1653  for (Int_t i=0; (var = (RooRealVar*)Itr2->Next()); ++i) {
1654  TString varName = var->GetName() + TString("err");
1655 
1656  if(!varName.Contains("alpha") && !varName.Contains("gamma_stat")) continue;
1657  if(!varName.Contains("err")) continue;
1658  if(varName.Contains("Lumi") || varName.Contains("nom") || varName.Contains("binWidth") || varName.Contains("pull")) continue;
1659 
1660  tree->SetBranchAddress(varName.Data(), &varValsNuisance[i]);
1661 
1662  if(varName.Contains("gamma_stat"))
1663  ngammas++;
1664 
1665  TH1F* temphisto = new TH1F("temphisto", "temphisto", 100, -5, 5);
1666 
1667  for(Int_t j=0; j < tree->GetEntries(); j++){
1668  tree->GetEntry(j);
1669  if(status%1000 == 0)
1670  temphisto->Fill(varValsNuisance[i]);
1671  //nuisancehisto->SetBinError(i+1, varValsNuisance[i]);
1672  }
1673 
1674  nuisancehisto->SetBinError(i+1, temphisto->GetMean());
1675  nuisancehisto->GetXaxis()->SetBinLabel(i+1, varName);
1676  delete temphisto;
1677  }
1678 
1679  delete Itr2;
1680 
1681  nuisancehisto->SetLineColor(1);
1682  nuisancehisto->SetMarkerStyle(21);
1683  nuisancehisto->SetMarkerColor(1);
1684  nuisancehisto->GetYaxis()->SetRangeUser(-2.0,2.0);
1685  nuisancehisto->GetXaxis()->SetTitle("Systematic ID");
1686  nuisancehisto->GetYaxis()->SetTitle("Fit Result");
1687  if(tree->GetEntries() > 1)
1688  nuisancehisto->GetYaxis()->SetTitle("Mean Fit Result From Toys");
1689  nuisancehisto->SetTitle("Nuisance Parameters");
1690  nuisancehisto->SetStats(false);
1691 
1692  TLine *mline = new TLine(0,0.0,n,0.0);
1693  mline->SetLineColor(kRed);
1694  TLine *sline = new TLine(0,1.0,n,1.0);
1695  sline->SetLineColor(kRed);
1696  TLine *s1line = new TLine(0,-1.0,n,-1.0);
1697  s1line->SetLineColor(kRed);
1698  TLine *vline = new TLine(n-ngammas,-2,n-ngammas,2.0);
1699  vline->SetLineColor(kGreen);
1700 
1701  TCanvas* cNuisanceParameters = new TCanvas("cNuisanceParameters","cNuisanceParamters");
1702  nuisancehisto->Draw("e");
1703  mline->Draw("same");
1704  sline->Draw("same");
1705  s1line->Draw("same");
1706  vline->Draw("same");
1707 
1708  return cNuisanceParameters;
1709 
1710 }
1711 
1712 //********************************************************************
1713 TCanvas* protoana::ProtoDUNEFitUtils::PlotAverageResultsFromToys(TTree* tree, RooWorkspace* ws, TString channelname, TString catname, RooArgList * PreFit_POI){
1714 //********************************************************************
1715 
1716  if(!tree || !ws){
1717  std::cerr << "ERROR::No tree or workspace found. Plot failed!" << std::endl;
1718  return NULL;
1719  }
1720 
1721  RooStats::ModelConfig *combined_config = (RooStats::ModelConfig*)ws->obj("ModelConfig");
1722  if(!combined_config){
1723  std::cerr << "ERROR::No model config " << " ModelConfig " << " in workspace. No plots!" << std::endl;
1724  return NULL;
1725  }
1726 
1727  std::vector<Float_t> varVals;
1728  std::vector<Float_t> varErrVals;
1729  const RooArgSet* floatParsList = combined_config->GetParametersOfInterest();
1730  if(floatParsList->getSize() == 0){
1731  std::cerr << "ERROR::No floating parameters found. Plot failed!" << std::endl;
1732  return NULL;
1733  }
1734 
1735  varVals.resize( floatParsList->getSize(), -999. );
1736  varErrVals.resize( floatParsList->getSize(), -999. );
1737 
1738  Int_t status;
1739  tree->SetBranchAddress("status", &status);
1740 
1741  Int_t NGoodToys = 0;
1742  for(Int_t j=0; j < tree->GetEntries(); j++){
1743  tree->GetEntry(j);
1744  if(status%1000 == 0)
1745  NGoodToys++;
1746  }
1747 
1748  if(NGoodToys == 0){
1749  std::cerr << "ERROR::Number of good fit quality is " << NGoodToys << ". Plot failed!" << std::endl;
1750  return NULL;
1751  }
1752 
1753  RooRealVar* var1(0);
1754  TIterator* Itr1 = floatParsList->createIterator();
1755  std::vector<TString> namevec;
1756 
1757  for(Int_t i=0; (var1 = (RooRealVar*)Itr1->Next()); ++i) {
1758  TString varName = TString(var1->GetName());
1759  if(!varName.Contains(channelname.Data())) continue;
1760  if(!varName.Contains(catname.Data())) continue;
1761  //if(varName.Contains("Lumi") || varName.Contains("binWidth") || varName.Contains("corr") || varName.Contains("Gamma")) continue;
1762  if(var1->isConstant()) varName += TString("_constant");
1763  namevec.push_back(varName);
1764  }
1765  delete Itr1;
1766 
1767  Int_t n = namevec.size();
1768  TString histoname = channelname + "_MomHisto_" + catname;
1769  TH1D* histo = new TH1D(histoname.Data(), histoname.Data(), n/*+1*/, 0, n/*+1*/);
1770 
1771  TH1D * prefit_histo = 0x0;
1772  if (PreFit_POI) {
1773  TString prefit_name = "PreFit_" + histoname;
1774  prefit_histo = (TH1D*)histo->Clone(prefit_name.Data());
1775  }
1776 
1777  RooRealVar* var(0);
1778  TIterator* Itr = floatParsList->createIterator();
1779  Int_t counter = n;
1780 
1781  for(Int_t i=0; (var = (RooRealVar*)Itr->Next()); ++i) {
1782  TString varName = TString(var->GetName());
1783  std::cout << "Ave name: " << varName << std::endl;
1784  if(!varName.Contains(channelname.Data())) continue;
1785  if(!varName.Contains(catname.Data())) continue;
1786  //if(varName.Contains("Lumi") || varName.Contains("binWidth") || varName.Contains("corr") || varName.Contains("Gamma")) continue;
1787 
1788  tree->SetBranchAddress(varName.Data(), &varVals[i]);
1789  TString varName2 = varName + TString("err");
1790  tree->SetBranchAddress(varName2.Data(), &varErrVals[i]);
1791 
1792  Float_t av = 0.0;
1793  Float_t averr = 0.0;
1794  for(Int_t j=0; j < tree->GetEntries(); j++){
1795  tree->GetEntry(j);
1796  if(status%1000 == 0){
1797  av += varVals[i]/NGoodToys;
1798  averr += varErrVals[i]/NGoodToys;
1799  }
1800  }
1801 
1802  histo->SetBinContent(counter, av);
1803  histo->SetBinError(counter, averr);
1804 
1805  if (PreFit_POI) {
1806  RooRealVar * var2(0);
1807  TIterator * prefit_itr = PreFit_POI->createIterator();
1808  while ((var2 = (RooRealVar*)prefit_itr->Next())) {
1809  std::cout << "Prefit:" << var2->GetName() << std::endl;
1810  TString var2Name(var2->GetName());
1811  if (var2Name == varName) {
1812  std::cout << "Found match " << var2->GetName() << std::endl;
1813  prefit_histo->SetBinContent(counter, var2->getVal());
1814  }
1815  }
1816  }
1817 
1818  counter--;
1819  }
1820 
1821  delete Itr;
1822 
1823  TLine *line = new TLine(0,1.0,n,1.0);
1824  line->SetLineColor(kRed);
1825 
1826  //TLine* line1 = new TLine(1,-2,1,5);
1827  //line1->SetLineColor(kGreen);
1828  //line1->SetLineStyle(kDashed);
1829 
1830  //TLine* line2 = new TLine(5,-2,5,5);
1831  //line2->SetLineColor(kGreen);
1832  //line2->SetLineStyle(kDashed);
1833 
1834  //TLine* line3 = new TLine(4,-2,4,5);
1835  //line3->SetLineColor(kGreen);
1836  //line3->SetLineStyle(kDashed);
1837 
1838  histoname = TString("canvas_") + histoname;
1839  histo->SetTitle(channelname.Data());
1840  //histo->GetXaxis()->SetTitle("True Momentum bin");
1841  histo->GetYaxis()->SetTitle("Fit result");
1842  if(tree->GetEntries() > 1)
1843  histo->GetYaxis()->SetTitle("Average Fit result");
1844  histo->GetYaxis()->SetRangeUser(-1.5,4.0);
1845  for(Int_t i=n; i > 0; i--)
1846  histo->GetXaxis()->SetBinLabel(i, namevec[n-i]);
1847  histo->SetMarkerStyle(20);
1848  histo->SetStats(false);
1849  histo->GetXaxis()->SetLabelSize(0.02);
1850 
1851  TCanvas* c = new TCanvas(histoname.Data(), histoname.Data());
1852  histo->SetFillStyle(3144);
1853  histo->SetFillColor(kRed);
1854  histo->Draw("e2");
1855  line->Draw("same");
1856 
1857  //line1->Draw("same");
1858  //line2->Draw("same");
1859  //line3->Draw("same");
1860 
1861  TLegend * leg = new TLegend(.65, .65, .85, .85);
1862  leg->AddEntry(histo, "Postfit", "lpf");
1863 
1864  if (PreFit_POI) {
1865  prefit_histo->SetMarkerColor(kBlue);
1866  prefit_histo->SetMarkerStyle(20);
1867  prefit_histo->Draw("p same");
1868  leg->AddEntry(prefit_histo, "Prefit", "lpf");
1869  }
1870 
1871  leg->Draw("same");
1872 
1873  return c;
1874 
1875 }
1876 
1877 //********************************************************************
1878 RooArgList protoana::ProtoDUNEFitUtils::GetPostfitPOIList(RooArgList paramsfit, bool print){
1879  //********************************************************************
1880 
1881  if(print){
1882  std::cout << std::endl;
1883  std::cout << "INFO::Printing fractional error for each fir parameter:-" << std::endl;
1884  }
1885 
1886  RooArgList poilist;
1887  TIterator* iter = paramsfit.createIterator();
1888  RooRealVar* var(0);
1889  for(Int_t i=0; (var = (RooRealVar*)iter->Next()); ++i){
1890  TString varname(var->GetName());
1891  if(varname.Contains("alpha_") || varname.Contains("gamma_")) continue;
1892  if(varname.Contains("POI_")){
1893  poilist.add(*var);
1894  if(print)
1895  std::cout << varname.Data() << ": Fit result = " << var->getVal() << " +/- " << var->getError() << ". Fractional error = " << var->getError()/var->getVal() << std::endl;
1896  }
1897  }
1898 
1899  delete iter;
1900 
1901  return poilist;
1902 
1903 }
1904 
1905 //********************************************************************
1906 void protoana::ProtoDUNEFitUtils::MakeNuisanceParamsConstant(RooWorkspace* ws, TString exceptPar){
1907  //********************************************************************
1908 
1909  RooStats::ModelConfig *combined_config = (RooStats::ModelConfig*)ws->obj("ModelConfig");
1910  if(!combined_config){
1911  std::cerr << "ERROR::No model config " << " ModelConfig " << " in workspace. Exit!" << std::endl;
1912  return;
1913  }
1914 
1915  RooAbsPdf* pdf = combined_config->GetPdf();
1916  if(!pdf){
1917  std::cerr << "ERROR::No pdf found in ModelConfig. Exit!" << std::endl;
1918  return;
1919  }
1920 
1921  const RooArgSet* obs = combined_config->GetNuisanceParameters();
1922  if(!obs){
1923  std::cerr << "ERROR::No nuisance found in ModelConfig. Exit!" << std::endl;
1924  return;
1925  }
1926 
1927  RooArgList floatParList;
1928  if(obs)
1929  floatParList.add(*obs);
1930 
1931  TIterator* iter = floatParList.createIterator();
1932  RooAbsArg* arg;
1933  while((arg=(RooAbsArg*)iter->Next())){
1934  TString varname;
1935  if(arg->InheritsFrom("RooRealVar") && !arg->isConstant()){
1936  varname = TString(arg->GetName());
1937  }
1938  else{
1939  continue;
1940  }
1941 
1942  if(varname == exceptPar) continue;
1943 
1944  RooRealVar* var = ws->var(varname.Data());
1945  if(!var){
1946  std::cout << "WARNING::Can't find parameter " << varname.Data() << " in workspace." << std::endl;
1947  continue;
1948  }
1949 
1950  if(varname == ""){
1951  std::cout << "WARNING::Empty variable name. Skipping reset value." << std::endl;
1952  continue;
1953  }
1954  else if(varname.Contains("alpha")){
1955  var->setConstant(kTRUE);
1956  }
1957  else if(varname.Contains("gamma_stat")){
1958  var->setConstant(kTRUE);
1959  }
1960  else{
1961  continue;
1962  }
1963  }
1964 
1965  delete iter;
1966 
1967 }
1968 
1969 //********************************************************************
1970 void protoana::ProtoDUNEFitUtils::ResetValues(RooWorkspace* ws, const RooArgList& parList){
1971  //********************************************************************
1972 
1973  TIterator* iter = parList.createIterator();
1974  RooAbsArg* arg;
1975  while((arg=(RooAbsArg*)iter->Next())){
1976  TString varname;
1977  if(arg->InheritsFrom("RooRealVar") && !arg->isConstant()){
1978  varname = TString(arg->GetName());
1979  }
1980  else{
1981  continue;
1982  }
1983 
1984  RooRealVar* var = ws->var(varname.Data());
1985  if(!var){
1986  std::cout << "WARNING::Can't find parameter " << varname.Data() << " in workspace." << std::endl;
1987  continue;
1988  }
1989 
1990  if(varname == ""){
1991  std::cout << "WARNING::Empty variable name. Skipping reset value." << std::endl;
1992  continue;
1993  }
1994  else if(varname.Contains("alpha")){
1995  var->setVal(0.0);
1996  }
1997  else if(varname.Contains("gamma_stat")){
1998  var->setVal(1.0);
1999  }
2000  else{
2001  var->setVal(1.0);
2002  }
2003  }
2004 
2005  delete iter;
2006 
2007 }
2008 
2009 //********************************************************************
2010 void protoana::ProtoDUNEFitUtils::ResetValuesToNominal(RooWorkspace* ws, const RooArgSet& parSet){
2011  //********************************************************************
2012 
2013  TIterator* iter = parSet.createIterator();
2014  RooAbsArg* arg;
2015  while((arg=(RooAbsArg*)iter->Next())){
2016  TString varname;
2017  if(arg->InheritsFrom("RooRealVar") && !arg->isConstant()){
2018  varname = TString(arg->GetName());
2019  }
2020  else{
2021  continue;
2022  }
2023 
2024  RooRealVar* var = ws->var(varname.Data());
2025  if(!var){
2026  std::cout << "WARNING::Can't find parameter " << varname.Data() << " in workspace." << std::endl;
2027  continue;
2028  }
2029 
2030  if(varname == ""){
2031  std::cout << "WARNING::Empty variable name. Skipping reset value to nominal." << std::endl;
2032  continue;
2033  }
2034  else if(varname.Contains("nom_gamma_stat")){
2035  var->setVal(1.0);
2036  }
2037  else if(varname.Contains("nom")){
2038  var->setVal(0.0);
2039  }
2040  else{
2041  var->setVal(0.0);
2042  }
2043  }
2044 
2045  delete iter;
2046 
2047 }
2048 
2049 //********************************************************************
2050 void protoana::ProtoDUNEFitUtils::ResetError(RooWorkspace* ws, const RooArgList& parList){
2051  //********************************************************************
2052 
2053  TIterator* iter = parList.createIterator();
2054  RooAbsArg* arg;
2055  while((arg=(RooAbsArg*)iter->Next())){
2056  TString varname;
2057  if(arg->InheritsFrom("RooRealVar") && !arg->isConstant()){
2058  varname = TString(arg->GetName());
2059  }
2060  else{
2061  continue;
2062  }
2063 
2064  RooRealVar* var = ws->var(varname.Data());
2065  if(!var){
2066  std::cout << "WARNING::Can't find parameter " << varname.Data() << " in workspace." << std::endl;
2067  continue;
2068  }
2069 
2070  if(varname == ""){
2071  std::cout << "WARNING::Empty variable name. Skipping reset error." << std::endl;
2072  continue;
2073  }
2074  else if(varname.Contains("alpha")){
2075  var->setError(1.0);
2076  if(var->getMin() < var->getVal() - 6.) var->setMin(var->getVal() - 6.);
2077  if(var->getMax() > var->getVal() + 6.) var->setMax(var->getVal() + 6.);
2078  }
2079  else if(varname.Contains("gamma_stat")){
2080  // Constraint could be either Gaus or Poisson
2081  RooAbsReal* constraint = (RooAbsReal*) ws->obj(Form("%s_constraint",varname.Data()));
2082  if(!constraint){
2083  std::cout << "WARNING::Constraint for variable " << varname.Data() << " not found. Skip reset error." << std::endl;
2084  continue;
2085  }
2086 
2087  TString constraintString = TString(constraint->IsA()->GetName());
2088  if(constraintString == "") continue;
2089  else if(constraintString.Contains("RooGaussian")){
2090  RooAbsReal* ErrorVar = (RooAbsReal*)ws->obj(Form("%s_sigma",varname.Data()));
2091  if(!ErrorVar){
2092  std::cout << "WARNING::Constraint type " << constraintString.Data() << " for variable " << varname.Data() << " not found. Skip reset error." << std::endl;
2093  continue;
2094  }
2095  Double_t err = ErrorVar->getVal();
2096  var->setError(err);
2097  if(var->getMin() < var->getVal() - 6.) var->setMin(var->getVal() - 6.);
2098  if(var->getMax() > var->getVal() + 6.) var->setMax(var->getVal() + 6.);
2099  }
2100  else if(constraintString.Contains("RooPoisson")){
2101  RooAbsReal* ErrorVar = (RooAbsReal*)ws->obj(Form("nom_%s",varname.Data()));
2102  if(!ErrorVar){
2103  std::cout << "WARNING::Constraint type " << constraintString.Data() << " for variable " << varname.Data() << " not found. Skip reset error." << std::endl;
2104  continue;
2105  }
2106  Double_t err = 1/sqrt(ErrorVar->getVal());
2107  var->setError(err);
2108  if(var->getMin() < var->getVal() - 6.) var->setMin(var->getVal() - 6.);
2109  if(var->getMax() > var->getVal() + 6.) var->setMax(var->getVal() + 6.);
2110  }
2111  else{
2112  std::cout << "WARNING::Unknown constraint type " << constraintString.Data() << ". Set prefit uncertainty to 0.00001." << std::endl;
2113  }
2114  }
2115 
2116  }
2117 
2118  delete iter;
2119 
2120 }
2121 
2122 //********************************************************************
2123 void protoana::ProtoDUNEFitUtils::ResetAllValuesAndErrors(RooWorkspace* ws){
2124  //********************************************************************
2125 
2126  RooStats::ModelConfig *combined_config = (RooStats::ModelConfig*)ws->obj("ModelConfig");
2127  if(!combined_config){
2128  std::cerr << "ERROR::No model config " << " ModelConfig " << " in workspace. Will not reset values and errors." << std::endl;
2129  return;
2130  }
2131 
2132  RooAbsPdf* pdf = combined_config->GetPdf();
2133  if(!pdf){
2134  std::cerr << "ERROR::No pdf found in ModelConfig. Will not reset values and errors." << std::endl;
2135  return;
2136  }
2137 
2138  const RooArgSet* obs = combined_config->GetObservables();
2139  if(!obs){
2140  std::cerr << "ERROR::No observables found in ModelConfig. Will not reset values and errors." << std::endl;
2141  return;
2142  }
2143 
2144  const RooArgSet* pdfpars = pdf->getParameters(obs);
2145  if(!pdfpars){
2146  std::cerr << "ERROR::No pdf parameters found. Will not reset values and errors." << std::endl;
2147  return;
2148  }
2149 
2150  const RooArgSet* globalobs = combined_config->GetGlobalObservables();
2151  if(!globalobs){
2152  std::cerr << "ERROR::No global observables found in ModelConfig. Will not reset values and errors." << std::endl;
2153  return;
2154  }
2155 
2156  RooArgList floatParamsList;
2157 
2158  TIterator* iter = pdfpars->createIterator() ;
2159  RooAbsArg* absarg;
2160  while( (absarg=(RooAbsArg*)iter->Next()) ) {
2161  if(absarg->InheritsFrom("RooRealVar") && !absarg->isConstant()){
2162  floatParamsList.add(*absarg);
2163  }
2164  }
2165  delete iter;
2166 
2167  // Now reset all values and error
2168  protoana::ProtoDUNEFitUtils::ResetValues(ws, floatParamsList);
2169  protoana::ProtoDUNEFitUtils::ResetError(ws, floatParamsList);
2170  protoana::ProtoDUNEFitUtils::ResetValuesToNominal(ws, *globalobs);
2171 
2172 }
name
static QCString name
Definition: declinfo.cpp:673
protoana::ProtoDUNEFitUtils::PlotDatasetsAndPdfs
std::vector< TCanvas * > PlotDatasetsAndPdfs(RooWorkspace *work, TString name, TString error, TString plottodraw, std::vector< TString > binnames, std::vector< double > recobins, std::vector< TString > incidentBinNames, std::vector< TString > sidebandBinNames, std::vector< double > sidebandBins, TString measurement="PDFit", bool doNegativeReco=false, RooAbsData *data=NULL, RooFitResult *result=NULL)
Definition: ProtoDUNEFitUtils.cxx:472
GetArgonNumberDensity
double GetArgonNumberDensity(double argon_density, double argon_molecularmass)
Definition: ProtoDUNEFitUtils.cxx:29
protoana::ProtoDUNEFitUtils::GetPostfitPOIList
RooArgList GetPostfitPOIList(RooArgList paramsfit, bool print=false)
Definition: ProtoDUNEFitUtils.cxx:1878
protoana::ProtoDUNEFitUtils::ResetValues
void ResetValues(RooWorkspace *ws, const RooArgList &parList)
Definition: ProtoDUNEFitUtils.cxx:1970
entry
QList< Entry > entry
Definition: tclscanner.cpp:1150
protoana::ProtoDUNEFitUtils::SetInterpolationCode
void SetInterpolationCode(RooWorkspace *ws, int code)
keras_to_tensorflow.f
f
Definition: keras_to_tensorflow.py:162
varName
static QCString varName
Definition: vhdljjparser.cpp:49
result
static QCString result
Definition: fortranscanner.cpp:56614
protoana::ProtoDUNEFitUtils::SaveWorkspace
void SaveWorkspace(RooWorkspace *ws, TString outFileName)
Definition: ProtoDUNEFitUtils.cxx:235
protoana::ProtoDUNEFitUtils::GetStatsSystHistogram
TH1 * GetStatsSystHistogram(TH1 *nominal)
Definition: ProtoDUNEFitUtils.cxx:66
string
std::string string
Definition: nybbler.cc:12
instance
const std::string instance
Definition: BlipMaker_module.cc:35
error
error
Definition: include.cc:26
QList::find
int find(const type *d) const
Definition: qlist.h:88
protoana::ProtoDUNEFitUtils::LoadSnapshot
bool LoadSnapshot(RooWorkspace *ws, TString snapshotname)
Definition: ProtoDUNEFitUtils.cxx:214
l
static QStrList * l
Definition: config.cpp:1044
work
void work()
Definition: CurrentProcessingContext_t.cpp:80
protoana::ProtoDUNEFitUtils::MakeNuisanceParamsConstant
void MakeNuisanceParamsConstant(RooWorkspace *ws, TString exceptPar)
Definition: ProtoDUNEFitUtils.cxx:1906
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
protoana::ProtoDUNEFitUtils::PlotXSecs
std::vector< TH1 * > PlotXSecs(RooWorkspace *work, std::string name, std::vector< TString > binnames, std::vector< double > recobins, std::vector< TString > incidentBinNames, RooAbsData *data=0x0, RooFitResult *result=0x0)
Definition: ProtoDUNEFitUtils.cxx:371
e
const double e
Definition: gUpMuFluxGen.cxx:165
var1
int var1
Definition: 029_hideinit.c:7
wirecell.dfp.graph.key
def key(type, name=None)
Definition: graph.py:13
protoana::ProtoDUNEFitUtils::GetSystematicHistoFromNominal
TH1 * GetSystematicHistoFromNominal(TH1 *nominal, TH1 *syst, TString name)
Definition: ProtoDUNEFitUtils.cxx:88
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
protoana::ProtoDUNEFitUtils::RemoveEmptyBins
void RemoveEmptyBins(RooPlot *frame)
Definition: ProtoDUNEFitUtils.cxx:274
create_imaplist.status
int status
Definition: create_imaplist.py:5
protoana::ProtoDUNEFitUtils::PlotNuisanceParametersImpact
std::vector< TCanvas * > PlotNuisanceParametersImpact(RooWorkspace *work, RooFitResult *result, TString snapshotload, std::vector< std::string > SystsToConsider)
Definition: ProtoDUNEFitUtils.cxx:1111
protoana::ProtoDUNEFitUtils::GetFitCovariance
TH2 * GetFitCovariance(RooFitResult *result)
Definition: ProtoDUNEFitUtils.cxx:148
test.p
p
Definition: test.py:223
protoana::ProtoDUNEFitUtils::PlotParametersPull
TCanvas * PlotParametersPull(TTree *tree, RooWorkspace *ws)
Definition: ProtoDUNEFitUtils.cxx:1344
protoana::ProtoDUNEFitUtils::GetDataMCChi2
double GetDataMCChi2(RooWorkspace *work, TString channelname, RooAbsData *data=NULL)
Definition: ProtoDUNEFitUtils.cxx:298
code
CodeOutputInterface * code
Definition: tclscanner.cpp:1114
make_THn_beam_input.counter
int counter
Definition: make_THn_beam_input.py:56
err
void err(const char *fmt,...)
Definition: message.cpp:226
make_THn_beam_input.data
data
Definition: make_THn_beam_input.py:80
protoana::ProtoDUNEFitUtils::ResetValuesToNominal
void ResetValuesToNominal(RooWorkspace *ws, const RooArgSet &parSet)
Definition: ProtoDUNEFitUtils.cxx:2010
Data
Fw2dFFT::Data Data
Definition: test_Fw2dFFT.cxx:33
Offset
FixedTimeOffsetTool::Offset Offset
Definition: FixedTimeOffsetTool_tool.cc:12
protoana::ProtoDUNEFitUtils::PlotNuisanceParameters
TCanvas * PlotNuisanceParameters(TTree *tree, RooWorkspace *ws)
Definition: ProtoDUNEFitUtils.cxx:1598
var
int var
Definition: 018_def.c:9
protoana::ProtoDUNEFitUtils::PlotNLL
std::vector< TCanvas * > PlotNLL(RooWorkspace *work, TString name, RooFitResult *result, bool plotPLL=false)
Definition: ProtoDUNEFitUtils.cxx:1001
protoana::ProtoDUNEFitUtils::GetFitCorrelationMatrix
TH2 * GetFitCorrelationMatrix(RooFitResult *result)
Definition: ProtoDUNEFitUtils.cxx:172
pduneana::line
void line(double t, double *p, double &x, double &y, double &z)
Definition: PDSPAnalyzer_module.cc:4741
protoana::ProtoDUNEFitUtils::SaveSnapshot
void SaveSnapshot(RooWorkspace *ws, TString snapshotname)
Definition: ProtoDUNEFitUtils.cxx:184
train.x
list x
Definition: train.py:276
protoana::ProtoDUNEFitUtils::ResetError
void ResetError(RooWorkspace *ws, const RooArgList &parList)
Definition: ProtoDUNEFitUtils.cxx:2050
protoana::ProtoDUNEFitUtils::PlotAverageResultsFromToys
TCanvas * PlotAverageResultsFromToys(TTree *tree, RooWorkspace *ws, TString channelname, TString catname, RooArgList *PreFit_POI=0x0)
Definition: ProtoDUNEFitUtils.cxx:1713
art::to_string
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
protoana::ProtoDUNEFitUtils::PlotNuisanceParametersPull
TCanvas * PlotNuisanceParametersPull(TTree *tree, RooWorkspace *ws)
Definition: ProtoDUNEFitUtils.cxx:1467
var2
int var2
Definition: 029_hideinit.c:12
make_resolution.tree
tree
Definition: make_resolution.py:22
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
protoana::ProtoDUNEFitUtils::GetSystHistograms
std::vector< TH1 * > GetSystHistograms(std::string name)
Definition: ProtoDUNEFitUtils.cxx:39
protoana::ProtoDUNEFitUtils::ResetAllValuesAndErrors
void ResetAllValuesAndErrors(RooWorkspace *ws)
Definition: ProtoDUNEFitUtils.cxx:2123
Generated by   doxygen 1.8.11