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