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"
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
32using namespace std;
33
34namespace 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 ) {
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}
#define gDirectory
Definition: TDirectory.h:223
char name[80]
Definition: TGX11.cxx:109
int type
Definition: TGX11.cxx:120
TODO Here, we are missing some documentation.
Definition: Asimov.h:23
void SetParamValue(const std::string &param, double value)
Definition: Asimov.h:36
void SetFixedParam(const std::string &param, bool constant=true)
Definition: Asimov.h:35
This class encapsulates all information for the statistical interpretation of one experiment.
Definition: Channel.h:26
HistFactory::StatErrorConfig & GetStatErrorConfig()
get information about threshold for statistical uncertainties and constraint term
Definition: Channel.h:67
RooStats::HistFactory::Data & GetData()
get data object
Definition: Channel.h:55
std::string GetName()
get name of channel
Definition: Channel.h:39
std::vector< RooStats::HistFactory::Sample > & GetSamples()
get vector of samples for this channel
Definition: Channel.h:71
Configuration for a constrained, coherent shape variation of affected samples.
Definition: Systematics.h:205
The RooStats::HistFactory::Measurement class can be used to construct a model by combining multiple R...
Definition: Measurement.h:30
double GetLumi()
retrieve integrated luminosity
Definition: Measurement.h:88
Configuration for an un- constrained overall systematic to scale sample normalisations.
Definition: Systematics.h:77
std::string GetName() const
Definition: Systematics.h:84
Configuration for a constrained overall systematic to scale sample normalisations.
Definition: Systematics.h:49
std::string GetName() const
Definition: Systematics.h:56
std::vector< RooStats::HistFactory::OverallSys > & GetOverallSysList()
Definition: Sample.h:109
std::string GetName() const
get name of sample
Definition: Sample.h:83
const TH1 * GetHisto() const
Definition: Sample.cxx:99
RooStats::HistFactory::StatError & GetStatError()
Definition: Sample.h:118
std::vector< RooStats::HistFactory::ShapeFactor > & GetShapeFactorList()
Definition: Sample.h:116
std::string GetChannelName() const
get name of associated channel
Definition: Sample.h:103
std::vector< RooStats::HistFactory::NormFactor > & GetNormFactorList()
Definition: Sample.h:110
std::vector< RooStats::HistFactory::HistoSys > & GetHistoSysList()
Definition: Sample.h:112
bool GetNormalizeByTheory() const
does the normalization scale with luminosity
Definition: Sample.h:79
std::vector< RooStats::HistFactory::ShapeSys > & GetShapeSysList()
Definition: Sample.h:115
Constrained bin-by-bin variation of affected histogram.
Definition: Systematics.h:225
Constraint::Type GetConstraintType() const
Definition: Systematics.h:257
const TH1 * GetErrorHist() const
Definition: Systematics.h:249
Constraint::Type GetConstraintType() const
Definition: Systematics.h:379
const TH1 * GetErrorHist() const
Definition: Systematics.h:349
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/...".
Bool_t cd(const char *path=nullptr) override
Change current directory to "this" directory.
TObject * Get(const char *namecycle) override
Return pointer to object identified by namecycle.
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:48
virtual void ShowStreamerInfo()
Show the StreamerInfo of all classes written to this file.
Definition: TFile.cxx:3589
The TH1 histogram class.
Definition: TH1.h:56
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:8381
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:7074
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2665
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition: TH1.cxx:7438
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TObject.cxx:144
Basic string class.
Definition: TString.h:131
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
vector< EstimateSummary > * loadSavedInputs(TFile *outFile, std::string channel)
Definition: Helper.cxx:49
std::vector< RooStats::HistFactory::EstimateSummary > GetChannelEstimateSummaries(RooStats::HistFactory::Measurement &measurement, RooStats::HistFactory::Channel &channel)
Definition: Helper.cxx:257
vector< pair< std::string, std::string > > get_comb(vector< std::string > names)
Definition: Helper.cxx:38
void AddSubStrings(vector< std::string > &vs, std::string s)
Definition: Helper.cxx:165
TH1 * GetHisto(TFile *inFile, const std::string name)
Definition: Helper.cxx:115
void saveInputs(TFile *outFile, std::string channel, vector< EstimateSummary > summaries)
Definition: Helper.cxx:79
std::vector< std::string > GetChildrenFromString(std::string str)
Definition: Helper.cxx:179
void AddParamsToAsimov(RooStats::HistFactory::Asimov &asimov, std::string str)
Definition: Helper.cxx:199
Namespace for the RooStats classes.
Definition: Asimov.h:20
static constexpr double s
Definition: file.py:1
void Print(const char *opt=0) const
This method must be overridden when a class wants to print itself.
std::vector< NormFactor > normFactor
std::vector< std::string > systSourceForHist
std::map< std::string, std::pair< double, double > > overallSyst