Logo ROOT  
Reference Guide
Helper.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id: cranmer $
2 // Author: Kyle Cranmer, Akira Shibata
3 /*************************************************************************
4  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5  * All rights reserved. *
6  * *
7  * For the licensing terms see $ROOTSYS/LICENSE. *
8  * For the list of contributors see $ROOTSYS/README/CREDITS. *
9  *************************************************************************/
10 
11 ////////////////////////////////////////////////////////////////////////////////
12 
13 /** \class RooStats::HistFactory::Helper
14  * \ingroup HistFactory
15  */
16 
17 
18 
19 //#define DEBUG
20 
21 #include <sstream>
22 #include "Helper.h"
23 #include "RooStats/ModelConfig.h"
24 
25 #include "RooArgSet.h"
26 #include "RooRealVar.h"
27 #include "RooMinimizer.h"
28 #include "RooSimultaneous.h"
29 #include "RooCategory.h"
30 #include "RooFitResult.h"
31 
32 using namespace std;
33 
34 namespace RooStats {
35  namespace HistFactory {
36 
37 
38  vector< pair<std::string, std::string> > get_comb(vector<std::string> names){
39  vector< pair<std::string, std::string> > list;
40  for(vector<std::string>::iterator itr=names.begin(); itr!=names.end(); ++itr){
41  vector<std::string>::iterator itr2=itr;
42  for(++itr2; itr2!=names.end();++itr2){
43  list.push_back(pair<std::string, std::string>(*itr, *itr2));
44  }
45  }
46  return list;
47  }
48 
49  vector<EstimateSummary>* loadSavedInputs(TFile* outFile, std::string channel ){
50  outFile->ShowStreamerInfo();
51 
52  vector<EstimateSummary>* summaries = new vector<EstimateSummary>;
53  outFile->cd(channel.c_str());
54 
55  // loop through estimate summaries
56  TIter next(gDirectory->GetListOfKeys());
57  EstimateSummary* summary;
58  while ((summary=(EstimateSummary*) next())) {
59  // if(summary){
60  summary->Print();
61  cout << "was able to read summary with name " << summary->name << std::endl;
62  cout << " nominal hist = " << summary->nominal << std::endl;
63  if(summary->nominal)
64  cout << " hist name = " << summary->nominal->GetName() <<endl;
65  cout << "still ok" << std::endl;
66 
67  summaries->push_back(*summary);
68 
69  //L.M. This code cannot be reached- remove it
70  // }
71  // else{
72  // cout << "was not able to read summary" << std::endl;
73  // }
74  }
75  return summaries;
76  }
77 
78 
79  void saveInputs(TFile* outFile, std::string channel, vector<EstimateSummary> summaries){
80  vector<EstimateSummary>::iterator it = summaries.begin();
81  vector<EstimateSummary>::iterator end = summaries.end();
82  vector<TH1*>::iterator histIt;
83  vector<TH1*>::iterator histEnd;
84  outFile->mkdir(channel.c_str());
85 
86  for(; it!=end; ++it){
87  if(channel != it->channel){
88  cout << "channel mismatch" << std::endl;
89  return;
90  }
91  outFile->cd(channel.c_str());
92 
93  // write the EstimateSummary object
94  it->Write();
95 
96  gDirectory->mkdir(it->name.c_str());
97  gDirectory->cd(it->name.c_str());
98 
99  it->nominal->Write();
100 
101  histIt = it->lowHists.begin();
102  histEnd = it->lowHists.end();
103  for(; histIt!=histEnd; ++histIt)
104  (*histIt)->Write();
105 
106  histIt = it->highHists.begin();
107  histEnd = it->highHists.end();
108  for(; histIt!=histEnd; ++histIt)
109  (*histIt)->Write();
110 
111  }
112  }
113 
114 
115  TH1 * GetHisto( TFile * inFile, const std::string name ){
116 
117  if(!inFile || name.empty()){
118  cerr << "Not all necessary info are set to access the input file. Check your config" << std::endl;
119  cerr << "fileptr: " << inFile
120  << "path/obj: " << name << std::endl;
121  return 0;
122  }
123 #ifdef DEBUG
124  cout << "Retrieving " << name ;
125 #endif
126  TH1 * ptr = (TH1 *) (inFile->Get( name.c_str() )->Clone());
127 #ifdef DEBUG
128  cout << " found at " << ptr << " with integral " << ptr->Integral() << " and mean " << ptr->GetMean() << std::endl;
129 #endif
130  if (ptr) ptr->SetDirectory(0); // for the current histogram h
131  //TH1::AddDirectory(kFALSE);
132  return ptr;
133 
134  }
135 
136  TH1 * GetHisto( const std::string file, const std::string path, const std::string obj){
137 
138 #ifdef DEBUG
139  cout << "Retrieving " << file << ":" << path << obj ;
140 #endif
141 
142  TFile inFile(file.c_str());
143  TH1 * ptr = (TH1 *) (inFile.Get( (path+obj).c_str() )->Clone());
144 
145 #ifdef DEBUG
146  cout << " found at " << ptr << " with integral " << ptr->Integral()
147  << " and mean " << ptr->GetMean() << std::endl;
148 #endif
149  // if(file.empty() || path.empty() || obj.empty()){
150 
151  if(!ptr){
152  cerr << "Not all necessary info are set to access the input file. Check your config"
153  << std::endl;
154  cerr << "filename: " << file
155  << "path: " << path
156  << "obj: " << obj << std::endl;
157  }
158  else
159  ptr->SetDirectory(0); // for the current histogram h
160 
161  return ptr;
162 
163  }
164 
165  void AddSubStrings( vector<std::string> & vs, std::string s){
166  const std::string delims("\\ ");
167  std::string::size_type begIdx, endIdx;
168  begIdx=s.find_first_not_of(delims);
169  while(begIdx!=string::npos){
170  endIdx=s.find_first_of(delims, begIdx);
171  if(endIdx==string::npos) endIdx=s.length();
172  vs.push_back(s.substr(begIdx,endIdx-begIdx));
173  begIdx=s.find_first_not_of(delims, endIdx);
174  }
175  }
176 
177  // Turn a std::string of "children" (space separated items)
178  // into a vector of std::strings
179  std::vector<std::string> GetChildrenFromString( std::string str ) {
180 
181  std::vector<std::string> child_vec;
182 
183  const std::string delims("\\ ");
184  std::string::size_type begIdx, endIdx;
185  begIdx=str.find_first_not_of(delims);
186  while(begIdx!=string::npos){
187  endIdx=str.find_first_of(delims, begIdx);
188  if(endIdx==string::npos) endIdx=str.length();
189  std::string child_name = str.substr(begIdx,endIdx-begIdx);
190  child_vec.push_back(child_name);
191  begIdx=str.find_first_not_of(delims, endIdx);
192  }
193 
194  return child_vec;
195  }
196 
197  // Turn a std::string of "children" (space separated items)
198  // into a vector of std::strings
199  void AddParamsToAsimov( RooStats::HistFactory::Asimov& asimov, std::string str ) {
200 
201  // First, split the string into a list
202  // each describing a parameter
203  std::vector<std::string> string_list = GetChildrenFromString( str );
204 
205  // Next, go through each one and split based
206  // on the '=' to separate the name from the val
207  // and fill the map
208  std::map<std::string, double> param_map;
209 
210  for( unsigned int i=0; i < string_list.size(); ++i) {
211 
212  std::string param = string_list.at(i);
213  // Split the string
214  size_t eql_location = param.find("=");
215 
216  // If there is no '=' deliminator, we only
217  // set the variable constant
218  if( eql_location==string::npos ) {
219  asimov.SetFixedParam(param);
220  }
221  else {
222 
223  std::string param_name = param.substr(0,eql_location);
224  double param_val = atof( param.substr(eql_location+1, param.size()).c_str() );
225 
226  std::cout << "ASIMOV - Param Name: " << param_name
227  << " Param Val: " << param_val << std::endl;
228  // Give the params a value AND set them constant
229  asimov.SetParamValue(param_name, param_val);
230  asimov.SetFixedParam(param_name);
231  }
232 
233  }
234 
235  return;
236 
237  }
238 
239 
240 
241  /*
242  bool AddSummaries( vector<EstimateSummary> & channel, vector<vector<EstimateSummary> > &master){
243  std::string channel_str=channel[0].channel;
244  for( unsigned int proc=1; proc < channel.size(); proc++){
245  if(channel[proc].channel != channel_str){
246  std::cout << "Illegal channel description, should be " << channel_str << " but found " << channel.at(proc).channel << std::endl;
247  std::cout << "name " << channel.at(proc).name << std::endl;
248  exit(1);
249  }
250  master.push_back(channel);
251  }
252  return true;
253  }*/
254 
255 
256  // Looking to deprecate this function and convert entirely to Measurements
257  std::vector<EstimateSummary> GetChannelEstimateSummaries(Measurement& measurement,
258  Channel& channel) {
259 
260  // Convert a "Channel" into a list of "Estimate Summaries"
261  // This should only be a temporary function, as the
262  // EstimateSummary class should be deprecated
263 
264  std::vector<EstimateSummary> channel_estimateSummary;
265 
266  std::cout << "Processing data: " << std::endl;
267 
268  // Add the data
269  EstimateSummary data_es;
270  data_es.name = "Data";
271  data_es.channel = channel.GetName();
272  TH1* data_hist = (TH1*) channel.GetData().GetHisto();
273  if( data_hist != NULL ) {
274  //data_es.nominal = (TH1*) channel.GetData().GetHisto()->Clone();
275  data_es.nominal = data_hist;
276  channel_estimateSummary.push_back( data_es );
277  }
278 
279  // Add the samples
280  for( unsigned int sampleItr = 0; sampleItr < channel.GetSamples().size(); ++sampleItr ) {
281 
282  EstimateSummary sample_es;
283  RooStats::HistFactory::Sample& sample = channel.GetSamples().at( sampleItr );
284 
285  std::cout << "Processing sample: " << sample.GetName() << std::endl;
286 
287  // Define the mapping
288  sample_es.name = sample.GetName();
289  sample_es.channel = sample.GetChannelName();
290  sample_es.nominal = (TH1*) sample.GetHisto()->Clone();
291 
292  std::cout << "Checking NormalizeByTheory" << std::endl;
293 
294  if( sample.GetNormalizeByTheory() ) {
295  sample_es.normName = "" ; // Really bad, confusion convention
296  }
297  else {
298  TString lumiStr;
299  lumiStr += measurement.GetLumi();
300  lumiStr.ReplaceAll(' ', TString());
301  sample_es.normName = lumiStr ;
302  }
303 
304  std::cout << "Setting the Histo Systs" << std::endl;
305 
306  // Set the Histo Systs:
307  for( unsigned int histoItr = 0; histoItr < sample.GetHistoSysList().size(); ++histoItr ) {
308 
309  RooStats::HistFactory::HistoSys& histoSys = sample.GetHistoSysList().at( histoItr );
310 
311  sample_es.systSourceForHist.push_back( histoSys.GetName() );
312  sample_es.lowHists.push_back( (TH1*) histoSys.GetHistoLow()->Clone() );
313  sample_es.highHists.push_back( (TH1*) histoSys.GetHistoHigh()->Clone() );
314 
315  }
316 
317  std::cout << "Setting the NormFactors" << std::endl;
318 
319  for( unsigned int normItr = 0; normItr < sample.GetNormFactorList().size(); ++normItr ) {
320 
321  RooStats::HistFactory::NormFactor& normFactor = sample.GetNormFactorList().at( normItr );
322 
323  EstimateSummary::NormFactor normFactor_es;
324  normFactor_es.name = normFactor.GetName();
325  normFactor_es.val = normFactor.GetVal();
326  normFactor_es.high = normFactor.GetHigh();
327  normFactor_es.low = normFactor.GetLow();
328  normFactor_es.constant = normFactor.GetConst();
329 
330  sample_es.normFactor.push_back( normFactor_es );
331 
332  }
333 
334  std::cout << "Setting the OverallSysList" << std::endl;
335 
336  for( unsigned int sysItr = 0; sysItr < sample.GetOverallSysList().size(); ++sysItr ) {
337 
338  RooStats::HistFactory::OverallSys& overallSys = sample.GetOverallSysList().at( sysItr );
339 
340  std::pair<double, double> DownUpPair( overallSys.GetLow(), overallSys.GetHigh() );
341  sample_es.overallSyst[ overallSys.GetName() ] = DownUpPair; //
342 
343  }
344 
345  std::cout << "Checking Stat Errors" << std::endl;
346 
347  // Do Stat Error
348  sample_es.IncludeStatError = sample.GetStatError().GetActivate();
349 
350  // Set the error and error threshold
352  if( sample.GetStatError().GetErrorHist() ) {
353  sample_es.relStatError = (TH1*) sample.GetStatError().GetErrorHist()->Clone();
354  }
355  else {
356  sample_es.relStatError = NULL;
357  }
358 
359 
360  // Set the constraint type;
362 
363  // Set the default
365 
366  if( type == Constraint::Gaussian) {
367  std::cout << "Using Gaussian StatErrors" << std::endl;
369  }
370  if( type == Constraint::Poisson ) {
371  std::cout << "Using Poisson StatErrors" << std::endl;
373  }
374 
375 
376  std::cout << "Getting the shape Factor" << std::endl;
377 
378  // Get the shape factor
379  if( sample.GetShapeFactorList().size() > 0 ) {
380  sample_es.shapeFactorName = sample.GetShapeFactorList().at(0).GetName();
381  }
382  if( sample.GetShapeFactorList().size() > 1 ) {
383  std::cout << "Error: Only One Shape Factor currently supported" << std::endl;
384  throw hf_exc();
385  }
386 
387 
388  std::cout << "Setting the ShapeSysts" << std::endl;
389 
390  // Get the shape systs:
391  for( unsigned int shapeItr=0; shapeItr < sample.GetShapeSysList().size(); ++shapeItr ) {
392 
393  RooStats::HistFactory::ShapeSys& shapeSys = sample.GetShapeSysList().at( shapeItr );
394 
395  EstimateSummary::ShapeSys shapeSys_es;
396  shapeSys_es.name = shapeSys.GetName();
397  shapeSys_es.hist = shapeSys.GetErrorHist();
398 
399  // Set the constraint type;
400  Constraint::Type systype = shapeSys.GetConstraintType();
401 
402  // Set the default
404 
405  if( systype == Constraint::Gaussian) {
407  }
408  if( systype == Constraint::Poisson ) {
409  shapeSys_es.constraint = EstimateSummary::Poisson;
410  }
411 
412  sample_es.shapeSysts.push_back( shapeSys_es );
413 
414  }
415 
416  std::cout << "Adding this sample" << std::endl;
417 
418  // Push back
419  channel_estimateSummary.push_back( sample_es );
420 
421  }
422 
423  return channel_estimateSummary;
424 
425  }
426 
427 
428 
429  /*
430  void unfoldConstraints(RooArgSet& initial, RooArgSet& final, RooArgSet& obs, RooArgSet& nuis, int& counter) {
431  //
432  // Loop through an initial set of constraints
433  // and create a final list of constraints
434  // This is done recursively
435 
436 
437  if (counter > 50) {
438  cout << "ERROR::Couldn't unfold constraints!" << std::endl;
439  cout << "Initial: " << std::endl;
440  initial.Print("v");
441  cout << std::endl;
442  cout << "Final: " << std::endl;
443  final.Print("v");
444  return;
445  }
446 
447  TIterator* itr = initial.createIterator();
448  RooAbsPdf* pdf;
449  while ((pdf = (RooAbsPdf*)itr->Next())) {
450  RooArgSet nuis_tmp = nuis;
451  RooArgSet constraint_set(*pdf->getAllConstraints(obs, nuis_tmp, false));
452  //if (constraint_set.getSize() > 1)
453  //{
454  string className(pdf->ClassName());
455  if (className != "RooGaussian" && className != "RooLognormal" && className != "RooGamma" && className != "RooPoisson" && className != "RooBifurGauss") {
456  counter++;
457  unfoldConstraints(constraint_set, final, obs, nuis, counter);
458  }
459  else {
460  final.add(*pdf);
461  }
462  }
463  delete itr;
464  }
465  */
466  /*
467  RooAbsData* makeAsimovData(ModelConfig* mcInWs, bool doConditional, RooWorkspace* combWS,
468  RooAbsPdf* combPdf, RooDataSet* combData, bool b_only, double doMuHat,
469  double muVal, bool signalInjection, bool doNuisPro) {
470  ////////////////////
471  //make asimov data//
472  ////////////////////
473 
474  // mcInWs -> The ModelCofig for this likelihood
475  // doConditional -> Minimize parameters for asimov quantities
476  // b_only -> Make the 'background only' asimov dataset, ie mu=0 (set muVal = 0)
477  // doNuisPro -> Set all nuisance parameters to '0' and to constant
478  // before minimizing. This should be done with *care*!!
479  // i.e. It should probably be removed as an option.
480  // signalInjection -> If true, then do the following:
481  // Perform the fit with m=0
482  // After the fit, set the value to mu=muVal
483  // so that the asimov is created with that value of mu fixed
484  // doMuHat -> Set 'mu' to be float during the fit (in the range -10 to 100)
485  // Even if it floats in the fit, it is still set to
486  // 'muVal' before the dataset is made (so the only difference
487  // comes from the other parameters that can float simultaneously with mu
488 
489  // Defaults:
490  // double doMuHat = false
491  // double muVal = -999,
492  // bool signalInjection = false
493  // bool doNuisPro = true
494 
495  if( b_only ) muVal = 0.0;
496 
497  int _printLevel = 0;
498 
499  // If using signal injection or a non-zero mu value,
500  // add a suffix showing that value
501  std::stringstream muStr;
502  if(signalInjection || !b_only) {
503  muStr << "_" << muVal;
504  }
505 
506  // Create the name of the resulting dataset
507  std::string dataSetName;
508  if(signalInjection) dataSetName = "signalInjection" + muStr.str();
509  else dataSetName = "asimovData" + muStr.str();
510 
511  // Set the parameter of interest
512  // to the 'background' value
513  RooRealVar* mu = (RooRealVar*) mcInWs->GetParametersOfInterest()->first();
514  if(signalInjection){
515  std::cout << "Asimov: Setting " << mu->GetName() << " value to 0 for fit" << std::endl;
516  mu->setVal(0);
517  }
518  else {
519  std::cout << "Asimov: Setting " << mu->GetName() << " value to " << muVal << " for fit" << std::endl;
520  mu->setVal(muVal);
521  }
522 
523  // Get necessary info from the ModelConfig
524  RooArgSet mc_obs = *mcInWs->GetObservables();
525  RooArgSet mc_globs = *mcInWs->GetGlobalObservables();
526  RooArgSet mc_nuis = *mcInWs->GetNuisanceParameters();
527 
528  // Create a set of constraint terms, which
529  // is stored in 'constraint_set'
530  // Make some temporary variables and use the
531  // unfoldConstrants function to do this.
532  RooArgSet constraint_set;
533  int counter_tmp = 0;
534  RooArgSet mc_nuis_tmp = mc_nuis;
535  RooArgSet constraint_set_tmp(*combPdf->getAllConstraints(mc_obs, mc_nuis_tmp, false));
536  unfoldConstraints(constraint_set_tmp, constraint_set, mc_obs, mc_nuis_tmp, counter_tmp);
537 
538  // Now that we have the constraint terms, we
539  // can create the full lists of nuisance parameters
540  // and global variables
541  RooArgList nui_list("ordered_nuis");
542  RooArgList glob_list("ordered_globs");
543 
544  TIterator* cIter = constraint_set.createIterator();
545  RooAbsArg* arg;
546  while ((arg = (RooAbsArg*)cIter->Next())) {
547  RooAbsPdf* pdf = (RooAbsPdf*) arg;
548  if (!pdf) continue;
549 
550  TIterator* nIter = mc_nuis.createIterator();
551  RooRealVar* thisNui = NULL;
552  RooAbsArg* nui_arg;
553  while((nui_arg = (RooAbsArg*)nIter->Next())) {
554  if(pdf->dependsOn(*nui_arg)) {
555  thisNui = (RooRealVar*) nui_arg;
556  break;
557  }
558  }
559  delete nIter;
560 
561  // need this in case the observable isn't fundamental.
562  // in this case, see which variable is dependent on the nuisance parameter and use that.
563  RooArgSet* components = pdf->getComponents();
564  components->remove(*pdf);
565  if(components->getSize()) {
566  TIterator* itr1 = components->createIterator();
567  RooAbsArg* arg1;
568  while ((arg1 = (RooAbsArg*)itr1->Next())) {
569  TIterator* itr2 = components->createIterator();
570  RooAbsArg* arg2;
571  while ((arg2 = (RooAbsArg*)itr2->Next())) {
572  if(arg1 == arg2) continue;
573  if(arg2->dependsOn(*arg1)) {
574  components->remove(*arg1);
575  }
576  }
577  delete itr2;
578  }
579  delete itr1;
580  }
581  if (components->getSize() > 1) {
582  std::cout << "ERROR::Couldn't isolate proper nuisance parameter" << std::endl;
583  return NULL;
584  }
585  else if (components->getSize() == 1) {
586  thisNui = (RooRealVar*)components->first();
587  }
588 
589  TIterator* gIter = mc_globs.createIterator();
590  RooRealVar* thisGlob = NULL;
591  RooAbsArg* glob_arg;
592  while ((glob_arg = (RooAbsArg*)gIter->Next()))
593  {
594  if (pdf->dependsOn(*glob_arg))
595  {
596  thisGlob = (RooRealVar*)glob_arg;
597  break;
598  }
599  }
600  delete gIter;
601 
602  if (!thisNui || !thisGlob)
603  {
604  std::cout << "WARNING::Couldn't find nui or glob for constraint: " << pdf->GetName() << std::endl;
605  continue;
606  }
607 
608  if (_printLevel >= 1) std::cout << "Pairing nui: " << thisNui->GetName() << ", with glob: " << thisGlob->GetName() << ", from constraint: " << pdf->GetName() << std::endl;
609 
610  nui_list.add(*thisNui);
611  glob_list.add(*thisGlob);
612 
613  } // End Loop over Constraint Terms
614  delete cIter;
615 
616  //save the snapshots of nominal parameters
617  combWS->saveSnapshot("nominalGlobs",glob_list);
618  combWS->saveSnapshot("nominalNuis", nui_list);
619 
620  RooArgSet nuiSet_tmp(nui_list);
621 
622  // Interesting line here:
623  if(!doMuHat) {
624  std::cout << "Asimov: Setting mu constant in fit" << std::endl;
625  mu->setConstant(true);
626  }
627  else {
628  std::cout << "Asimov: Letting mu float in fit (muHat)" << std::endl;
629  mu->setRange(-10,100);
630  }
631 
632  // Conditional: "Minimize the parameters"
633  if(doConditional) {
634 
635  std::cout << "Starting minimization.." << std::endl;
636 
637  // Consider removing this option:
638  if(!doNuisPro) {
639  std::cout << "Asimov: Setting nuisance parameters constant in the fit (ARE YOU SURE YOU WANT THIS)" << std::endl;
640  TIterator* nIter = nuiSet_tmp.createIterator();
641  RooRealVar* thisNui = NULL;
642  while((thisNui = (RooRealVar*) nIter->Next())) {
643  thisNui->setVal(0);
644  thisNui->setConstant();
645  }
646  delete nIter;
647  // This should be checked, we don't want to
648  if(combWS->var("Lumi")) {
649  combWS->var("Lumi")->setVal(1);
650  combWS->var("Lumi")->setConstant();
651  }
652  }
653 
654  // Create the nll and its minimizer
655  RooAbsReal* nll = combPdf->createNLL(*combData, RooFit::Constrain(nuiSet_tmp));
656  RooMinimizer minim(*nll);
657  minim.setStrategy(2);
658  minim.setPrintLevel(999);
659 
660  // Do the minimization
661  std::cout << "Minimizing to make Asimov dataset:" << std::endl;
662  int status = minim.minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), "migrad");
663  if (status == 0) {
664  // status==0 means fit was successful
665  std::cout << "Successfully minimized to make Asimov dataset:" << std::endl;
666  RooFitResult* fit_result = minim.lastMinuitFit();
667  std::cout << "Asimov: Final Fitted Parameters" << std::endl;
668  fit_result->Print("V");
669  } else{
670  std::cout << "Fit failed for mu = " << mu->getVal() << " with status " << status << std::endl;
671  std::cout << "Trying minuit..." << std::endl;
672  status = minim.minimize("Minuit", "migrad");
673  if (status != 0) {
674  cout << "Fit failed for mu = " << mu->getVal() << " with status " << status << std::endl;
675  throw hf_exc();
676  }
677  }
678 
679  // Undo the 'doNuisPro' part
680  // Again, may want to remove this
681  if (!doNuisPro) {
682  TIterator* nIter = nuiSet_tmp.createIterator();
683  RooRealVar* thisNui = NULL;
684  while ((thisNui = (RooRealVar*)nIter->Next())) {
685  thisNui->setConstant(false);
686  }
687  delete nIter;
688  if (combWS->var("Lumi")) {
689  combWS->var("Lumi")->setConstant(false);
690  }
691  }
692 
693  std::cout << "Done" << std::endl;
694  } // END: DoConditional
695 
696  mu->setConstant(false);
697 
698  //loop over the nui/glob list, grab the corresponding variable from the tmp ws,
699  // and set the glob to the value of the nui
700  int nrNuis = nui_list.getSize();
701  if (nrNuis != glob_list.getSize()) {
702  std::cout << "ERROR::nui_list.getSize() != glob_list.getSize()!" << std::endl;
703  return NULL;
704  }
705 
706  for(int i=0; i<nrNuis; i++) {
707  RooRealVar* nui = (RooRealVar*) nui_list.at(i);
708  RooRealVar* glob = (RooRealVar*) glob_list.at(i);
709 
710  if (_printLevel >= 1) std::cout << "nui: " << nui << ", glob: " << glob << std::endl;
711  if (_printLevel >= 1) std::cout << "Setting glob: " << glob->GetName() << ", which had previous val: " << glob->getVal() << ", to conditional val: " << nui->getVal() << std::endl;
712 
713  glob->setVal(nui->getVal());
714  }
715 
716  //save the snapshots of conditional parameters
717  //cout << "Saving conditional snapshots" << std::endl;
718  combWS->saveSnapshot(("conditionalGlobs"+muStr.str()).c_str(),glob_list);
719  combWS->saveSnapshot(("conditionalNuis" +muStr.str()).c_str(), nui_list);
720 
721  if(!doConditional) {
722  combWS->loadSnapshot("nominalGlobs");
723  combWS->loadSnapshot("nominalNuis");
724  }
725 
726  //cout << "Making asimov" << std::endl;
727  //make the asimov data (snipped from Kyle)
728 
729  // Restore the value of mu to the target value
730  mu->setVal(muVal);
731 
732  ModelConfig* mc = mcInWs;
733 
734  int iFrame=0;
735 
736  const char* weightName = "weightVar";
737  RooArgSet obsAndWeight;
738  obsAndWeight.add(*mc->GetObservables());
739 
740  // Get the weightVar, or create one if necessary
741  RooRealVar* weightVar = combWS->var(weightName); // NULL;
742  // if (!(weightVar = combWS->var(weightName)))
743  if( weightVar==NULL ) {
744  combWS->import(*(new RooRealVar(weightName, weightName, 1,0,1000)));
745  weightVar = combWS->var(weightName);
746  }
747  if (_printLevel >= 1) std::cout << "weightVar: " << weightVar << std::endl;
748  obsAndWeight.add(*combWS->var(weightName));
749 
750  //////////////////////////////////////////////////////
751  //////////////////////////////////////////////////////
752  //////////////////////////////////////////////////////
753  //////////////////////////////////////////////////////
754  //////////////////////////////////////////////////////
755  // MAKE ASIMOV DATA FOR OBSERVABLES
756 
757  // dummy var can just have one bin since it's a dummy
758  if(combWS->var("ATLAS_dummyX")) combWS->var("ATLAS_dummyX")->setBins(1);
759 
760  if (_printLevel >= 1) std::cout << " check expectedData by category" << std::endl;
761  RooSimultaneous* simPdf = dynamic_cast<RooSimultaneous*>(mc->GetPdf());
762 
763  // Create the pointer to be returned
764  RooDataSet* asimovData=NULL;
765 
766  // If the pdf isn't simultaneous:
767  if(!simPdf) {
768 
769  // Get pdf associated with state from simpdf
770  RooAbsPdf* pdftmp = mc->GetPdf();//simPdf->getPdf(channelCat->getLabel()) ;
771 
772  // Generate observables defined by the pdf associated with this state
773  RooArgSet* obstmp = pdftmp->getObservables(*mc->GetObservables()) ;
774 
775  if (_printLevel >= 1) {
776  obstmp->Print();
777  }
778 
779 
780  asimovData = new RooDataSet(dataSetName.c_str(), dataSetName.c_str(),
781  RooArgSet(obsAndWeight), RooFit::WeightVar(*weightVar));
782 
783  RooRealVar* thisObs = ((RooRealVar*)obstmp->first());
784  double expectedEvents = pdftmp->expectedEvents(*obstmp);
785  double thisNorm = 0;
786  for(int jj=0; jj<thisObs->numBins(); ++jj){
787  thisObs->setBin(jj);
788 
789  thisNorm=pdftmp->getVal(obstmp)*thisObs->getBinWidth(jj);
790  if (thisNorm*expectedEvents <= 0)
791  {
792  cout << "WARNING::Detected bin with zero expected events (" << thisNorm*expectedEvents << ") ! Please check your inputs. Obs = " << thisObs->GetName() << ", bin = " << jj << std::endl;
793  }
794  if (thisNorm*expectedEvents > 0 && thisNorm*expectedEvents < pow(10.0, 18)) asimovData->add(*mc->GetObservables(), thisNorm*expectedEvents);
795  }
796 
797  if (_printLevel >= 1)
798  {
799  asimovData->Print();
800  std::cout <<"sum entries "<<asimovData->sumEntries()<<endl;
801  }
802  if(asimovData->sumEntries()!=asimovData->sumEntries()){
803  std::cout << "sum entries is nan"<<endl;
804  exit(1);
805  }
806 
807  // don't import, return (of course)
808  //combWS->import(*asimovData);
809  }
810  else
811  {
812 
813  // If it IS a simultaneous pdf
814 
815  std::cout << "found a simPdf: " << simPdf << std::endl;
816  map<std::string, RooDataSet*> asimovDataMap;
817 
818  RooCategory* channelCat = (RooCategory*)&simPdf->indexCat();
819  TIterator* iter = channelCat->typeIterator() ;
820  RooCatType* tt = NULL;
821  int nrIndices = 0;
822  while((tt=(RooCatType*) iter->Next())) {
823  nrIndices++;
824  }
825 
826  for (int i=0;i<nrIndices;i++){
827 
828  channelCat->setIndex(i);
829 
830  std::cout << "Checking channel: " << channelCat->getLabel() << std::endl;
831  iFrame++;
832  // Get pdf associated with state from simpdf
833  RooAbsPdf* pdftmp = simPdf->getPdf(channelCat->getLabel()) ;
834 
835  // Generate observables defined by the pdf associated with this state
836  RooArgSet* obstmp = pdftmp->getObservables(*mc->GetObservables()) ;
837 
838  if (_printLevel >= 1) {
839  obstmp->Print();
840  cout << "on type " << channelCat->getLabel() << " " << iFrame << std::endl;
841  }
842 
843  RooDataSet* obsDataUnbinned = new RooDataSet(Form("combAsimovData%d",iFrame),Form("combAsimovData%d",iFrame),
844  RooArgSet(obsAndWeight,*channelCat), RooFit::WeightVar(*weightVar));
845  RooRealVar* thisObs = ((RooRealVar*)obstmp->first());
846  double expectedEvents = pdftmp->expectedEvents(*obstmp);
847  double thisNorm = 0;
848  TString pdftmp_name = pdftmp->GetName();
849 
850  if (!expectedEvents) {
851  std::cout << "Not expected events" << std::endl;
852  if (pdftmp_name == "model_E")
853  ((RooRealVar*)obstmp->first())->setVal(combWS->function("p_e")->getVal());
854 
855  else if (pdftmp_name == "model_MU")
856  ((RooRealVar*)obstmp->first())->setVal(combWS->function("p_mu")->getVal());
857 
858  else if ((pdftmp_name == "model_ratio_ELMU") || (pdftmp_name == "model_comb")) {
859  //((RooRealVar*)obstmp->first())->setVal(combWS->function("p_comb")->getVal());
860  double p_asimov_val = combWS->var("p_asimov")->getVal();
861  std::cout << "p_asimov val: " << p_asimov_val << std::endl;
862  ((RooRealVar*)obstmp->first())->setVal(combWS->var("p_asimov")->getVal());
863  }
864 
865  else {
866  std::cout << "Failed to set asimov data for non-extended pdf" << std::endl;
867  exit(1);
868  }
869  obsDataUnbinned->add(*mc->GetObservables());
870 
871  }
872  else {
873  std::cout << "expected events" << std::endl;
874  for(int jj=0; jj<thisObs->numBins(); ++jj){
875  thisObs->setBin(jj);
876 
877  thisNorm=pdftmp->getVal(obstmp)*thisObs->getBinWidth(jj);
878  if (thisNorm*expectedEvents <= 0)
879  {
880  std::cout << "WARNING::Detected bin with zero expected events (" << thisNorm*expectedEvents << ") ! Please check your inputs. Obs = " << thisObs->GetName() << ", bin = " << jj << std::endl;
881  }
882  if (thisNorm*expectedEvents > pow(10.0, -9) && thisNorm*expectedEvents < pow(10.0, 9)) obsDataUnbinned->add(*mc->GetObservables(), thisNorm*expectedEvents);
883  }
884  }
885 
886  if (_printLevel >= 1)
887  {
888  obsDataUnbinned->Print();
889  std::cout <<"sum entries "<<obsDataUnbinned->sumEntries()<<endl;
890  }
891  if(obsDataUnbinned->sumEntries()!=obsDataUnbinned->sumEntries()){
892  cout << "sum entries is nan"<<endl;
893  exit(1);
894  }
895 
896  asimovDataMap[string(channelCat->getLabel())] = obsDataUnbinned;//tempData;
897 
898  if (_printLevel >= 1)
899  {
900  std::cout << "channel: " << channelCat->getLabel() << ", data: ";
901  obsDataUnbinned->Print();
902  std::cout << std::endl;
903  }
904  }
905 
906  channelCat->setIndex(0);
907 
908  asimovData = new RooDataSet(dataSetName.c_str(),dataSetName.c_str(),
909  RooArgSet(obsAndWeight,*channelCat), RooFit::Index(*channelCat),
910  RooFit::Import(asimovDataMap), RooFit::WeightVar(*weightVar));
911 
912  // Don't import, return (of course)
913  //combWS->import(*asimovData);
914  } // End if over simultaneous pdf
915 
916  combWS->loadSnapshot("nominalNuis");
917  combWS->loadSnapshot("nominalGlobs");
918 
919  return asimovData;
920 
921  }
922  */
923  }
924 }
RooStats::HistFactory::HistogramUncertaintyBase::GetHistoLow
const TH1 * GetHistoLow() const
Definition: Systematics.h:159
RooStats::HistFactory::EstimateSummary::shapeSysts
std::vector< ShapeSys > shapeSysts
Definition: EstimateSummary.h:81
RooStats::HistFactory::Sample
Definition: Sample.h:27
RooStats::HistFactory::Sample::GetShapeSysList
std::vector< RooStats::HistFactory::ShapeSys > & GetShapeSysList()
Definition: Sample.h:113
RooStats::HistFactory::EstimateSummary::channel
std::string channel
Definition: EstimateSummary.h:62
RooStats::HistFactory::EstimateSummary::NormFactor::constant
bool constant
Definition: EstimateSummary.h:33
RooStats::HistFactory::NormFactor::GetLow
double GetLow() const
Definition: Systematics.h:94
RooStats::HistFactory::StatErrorConfig::GetRelErrorThreshold
double GetRelErrorThreshold() const
Definition: Systematics.h:385
RooStats::HistFactory::EstimateSummary::name
std::string name
Definition: EstimateSummary.h:61
RooSimultaneous.h
RooStats::HistFactory::Sample::GetHisto
const TH1 * GetHisto() const
Definition: Sample.cxx:99
RooStats::HistFactory::Channel::GetSamples
std::vector< RooStats::HistFactory::Sample > & GetSamples()
get vector of samples for this channel
Definition: Channel.h:76
RooStats::HistFactory::EstimateSummary::shapeFactorName
std::string shapeFactorName
Definition: EstimateSummary.h:80
RooArgSet.h
RooStats::HistFactory::Asimov::SetFixedParam
void SetFixedParam(const std::string &param, bool constant=true)
Definition: Asimov.h:34
RooStats::HistFactory::Sample::GetOverallSysList
std::vector< RooStats::HistFactory::OverallSys > & GetOverallSysList()
Definition: Sample.h:109
RooStats::HistFactory::StatError::GetActivate
bool GetActivate() const
Definition: Systematics.h:341
RooStats::HistFactory::OverallSys::GetHigh
double GetHigh() const
Definition: Systematics.h:61
RooStats::HistFactory::EstimateSummary::ShapeSys::constraint
ConstraintType constraint
Definition: EstimateSummary.h:42
RooStats::HistFactory::Constraint::Poisson
@ Poisson
Definition: Systematics.h:25
RooStats::HistFactory::Sample::GetHistoSysList
std::vector< RooStats::HistFactory::HistoSys > & GetHistoSysList()
Definition: Sample.h:111
RooStats::HistFactory::EstimateSummary::relStatError
TH1 * relStatError
Definition: EstimateSummary.h:77
RooStats::HistFactory::NormFactor
Configuration for an un- constrained overall systematic to scale sample normalisations.
Definition: Systematics.h:77
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
RooStats::HistFactory::HistoSys
Configuration for a constrained, coherent shape variation of affected samples.
Definition: Systematics.h:205
TDirectoryFile::cd
Bool_t cd(const char *path=nullptr) override
Change current directory to "this" directory.
Definition: TDirectoryFile.cxx:340
RooStats::HistFactory::EstimateSummary::nominal
TH1 * nominal
Definition: EstimateSummary.h:64
RooStats::HistFactory::EstimateSummary::IncludeStatError
bool IncludeStatError
Definition: EstimateSummary.h:74
RooStats::HistFactory::hf_exc
Definition: HistFactoryException.h:23
RooStats::HistFactory::EstimateSummary::highHists
std::vector< TH1 * > highHists
Definition: EstimateSummary.h:67
RooStats::HistFactory::EstimateSummary::NormFactor::val
double val
Definition: EstimateSummary.h:32
RooStats::HistFactory::Sample::GetChannelName
std::string GetChannelName() const
get name of associated channel
Definition: Sample.h:103
RooStats::HistFactory::OverallSys::GetLow
double GetLow() const
Definition: Systematics.h:60
RooStats::HistFactory::GetChildrenFromString
std::vector< std::string > GetChildrenFromString(std::string str)
Definition: Helper.cxx:179
RooStats::HistFactory::NormFactor::GetHigh
double GetHigh() const
Definition: Systematics.h:95
RooStats::HistFactory::ShapeSys
Constrained bin-by-bin variation of affected histogram.
Definition: Systematics.h:225
RooStats::HistFactory::EstimateSummary::ShapeSys::name
std::string name
Definition: EstimateSummary.h:40
RooStats::HistFactory::loadSavedInputs
vector< EstimateSummary > * loadSavedInputs(TFile *outFile, std::string channel)
Definition: Helper.cxx:49
RooStats::HistFactory::AddParamsToAsimov
void AddParamsToAsimov(RooStats::HistFactory::Asimov &asimov, std::string str)
Definition: Helper.cxx:199
TString
Basic string class.
Definition: TString.h:136
RooStats::HistFactory::Sample::GetStatError
RooStats::HistFactory::StatError & GetStatError()
Definition: Sample.h:125
RooFit::HistFactory
@ HistFactory
Definition: RooGlobalFunc.h:69
TString::ReplaceAll
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:692
RooStats::HistFactory::get_comb
vector< pair< std::string, std::string > > get_comb(vector< std::string > names)
Definition: Helper.cxx:38
RooStats::HistFactory::EstimateSummary::NormFactor
Definition: EstimateSummary.h:30
RooStats::HistFactory::Sample::GetNormalizeByTheory
bool GetNormalizeByTheory() const
does the normalization scale with luminosity
Definition: Sample.h:79
RooStats::HistFactory::Channel
This class encapsulates all information for the statistical interpretation of one experiment.
Definition: Channel.h:30
TH1::Clone
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2740
RooStats::HistFactory::Asimov
TODO Here, we are missing some documentation.
Definition: Asimov.h:22
TH1::GetMean
virtual Double_t GetMean(Int_t axis=1) const
For axis = 1,2 or 3 returns the mean value of the histogram along X,Y or Z axis.
Definition: TH1.cxx:7428
RooStats::HistFactory::HistogramUncertaintyBase::GetName
const std::string & GetName() const
Definition: Systematics.h:163
RooStats::HistFactory::EstimateSummary::NormFactor::high
double high
Definition: EstimateSummary.h:32
RooStats::HistFactory::Measurement
The RooStats::HistFactory::Measurement class can be used to construct a model by combining multiple R...
Definition: Measurement.h:30
TDirectoryFile::Get
TObject * Get(const char *namecycle) override
Return pointer to object identified by namecycle.
Definition: TDirectoryFile.cxx:909
RooStats::HistFactory::EstimateSummary::RelErrorThreshold
Double_t RelErrorThreshold
Definition: EstimateSummary.h:76
RooStats::HistFactory::OverallSys::GetName
const std::string & GetName() const
Definition: Systematics.h:56
RooStats::HistFactory::EstimateSummary::normFactor
std::vector< NormFactor > normFactor
Definition: EstimateSummary.h:70
RooStats::HistFactory::EstimateSummary::Print
void Print(const char *opt=0) const
This method must be overridden when a class wants to print itself.
Definition: EstimateSummary.cxx:39
ModelConfig.h
size
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
RooStats::HistFactory::Constraint::Gaussian
@ Gaussian
Definition: Systematics.h:25
RooStats::HistFactory::EstimateSummary::systSourceForHist
std::vector< std::string > systSourceForHist
Definition: EstimateSummary.h:65
TFile::ShowStreamerInfo
virtual void ShowStreamerInfo()
Show the StreamerInfo of all classes written to this file.
Definition: TFile.cxx:3664
RooStats::HistFactory::Measurement::GetLumi
double GetLumi()
retrieve integrated luminosity
Definition: Measurement.h:88
RooStats::HistFactory::Constraint::Type
Type
Definition: Systematics.h:25
TObject::Clone
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TObject.cxx:146
gDirectory
#define gDirectory
Definition: TDirectory.h:327
RooStats::HistFactory::ShapeSys::GetErrorHist
const TH1 * GetErrorHist() const
Definition: Systematics.h:256
RooStats::HistFactory::EstimateSummary::overallSyst
std::map< std::string, std::pair< double, double > > overallSyst
Definition: EstimateSummary.h:68
RooStats::HistFactory::OverallSys
Configuration for a constrained overall systematic to scale sample normalisations.
Definition: Systematics.h:49
RooStats::HistFactory::HistogramUncertaintyBase::GetHistoHigh
const TH1 * GetHistoHigh() const
Definition: Systematics.h:160
TH1::SetDirectory
virtual void SetDirectory(TDirectory *dir)
By default when an histogram is created, it is added to the list of histogram objects in the current ...
Definition: TH1.cxx:8777
RooStats::HistFactory::Asimov::SetParamValue
void SetParamValue(const std::string &param, double value)
Definition: Asimov.h:35
RooCategory.h
RooStats::HistFactory::Channel::GetName
std::string GetName() const
get name of channel
Definition: Channel.h:43
RooStats::HistFactory::NormFactor::GetVal
double GetVal() const
Definition: Systematics.h:87
RooRealVar.h
RooStats::HistFactory::GetHisto
TH1 * GetHisto(const std::string file, const std::string path, const std::string obj)
Definition: Helper.cxx:136
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
RooStats::HistFactory::Sample::GetShapeFactorList
std::vector< RooStats::HistFactory::ShapeFactor > & GetShapeFactorList()
Definition: Sample.h:114
RooFitResult.h
RooStats::HistFactory::Sample::GetName
std::string GetName() const
get name of sample
Definition: Sample.h:83
RooStats::HistFactory::NormFactor::GetName
std::string GetName() const
Definition: Systematics.h:84
RooStats::HistFactory::Data::GetHisto
TH1 * GetHisto()
Definition: Data.cxx:28
RooStats::HistFactory::AddSubStrings
void AddSubStrings(vector< std::string > &vs, std::string s)
Definition: Helper.cxx:165
RooStats::HistFactory::saveInputs
void saveInputs(TFile *outFile, std::string channel, vector< EstimateSummary > summaries)
Definition: Helper.cxx:79
RooMinimizer.h
RooStats::HistFactory::Sample::GetNormFactorList
std::vector< RooStats::HistFactory::NormFactor > & GetNormFactorList()
Definition: Sample.h:110
RooStats::HistFactory::NormFactor::GetConst
bool GetConst() const
Definition: Systematics.h:90
RooStats
Namespace for the RooStats classes.
Definition: Asimov.h:19
file
Definition: file.py:1
TH1
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
name
char name[80]
Definition: TGX11.cxx:110
Helper.h
TIter
Definition: TCollection.h:233
RooStats::HistFactory::EstimateSummary::normName
std::string normName
Definition: EstimateSummary.h:63
RooStats::HistFactory::EstimateSummary::ShapeSys::hist
const TH1 * hist
Definition: EstimateSummary.h:41
RooStats::HistFactory::EstimateSummary
Definition: EstimateSummary.h:26
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
type
int type
Definition: TGX11.cxx:121
RooStats::HistFactory::StatErrorConfig::GetConstraintType
Constraint::Type GetConstraintType() const
Definition: Systematics.h:388
RooStats::HistFactory::EstimateSummary::StatConstraintType
ConstraintType StatConstraintType
Definition: EstimateSummary.h:75
TH1::Integral
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition: TH1.cxx:7834
RooStats::HistFactory::EstimateSummary::NormFactor::low
double low
Definition: EstimateSummary.h:32
RooStats::HistFactory::EstimateSummary::ShapeSys
Definition: EstimateSummary.h:37
RooStats::HistFactory::EstimateSummary::NormFactor::name
std::string name
Definition: EstimateSummary.h:31
RooStats::HistFactory::Channel::GetData
RooStats::HistFactory::Data & GetData()
get data object
Definition: Channel.h:59
RooStats::HistFactory::EstimateSummary::lowHists
std::vector< TH1 * > lowHists
Definition: EstimateSummary.h:66
RooStats::HistFactory::ShapeSys::GetConstraintType
Constraint::Type GetConstraintType() const
Definition: Systematics.h:264
RooStats::HistFactory::StatError::GetErrorHist
const TH1 * GetErrorHist() const
Definition: Systematics.h:356
RooStats::HistFactory::Channel::GetStatErrorConfig
HistFactory::StatErrorConfig & GetStatErrorConfig()
get information about threshold for statistical uncertainties and constraint term
Definition: Channel.h:71
TDirectoryFile::mkdir
TDirectory * mkdir(const char *name, const char *title="", Bool_t returnExistingDirectory=kFALSE) override
Create a sub-directory "a" or a hierarchy of sub-directories "a/b/c/...".
Definition: TDirectoryFile.cxx:1207
RooStats::HistFactory::GetChannelEstimateSummaries
std::vector< RooStats::HistFactory::EstimateSummary > GetChannelEstimateSummaries(RooStats::HistFactory::Measurement &measurement, RooStats::HistFactory::Channel &channel)
Definition: Helper.cxx:257