CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
asimovutils Namespace Reference

Functions

RooAbsData * asimovDatasetNominal (RooStats::ModelConfig *mc, double poiValue=0.0, int verbose=0)
 Generate asimov dataset from nominal value of nuisance parameters. More...
 
RooAbsData * asimovDatasetWithFit (RooStats::ModelConfig *mc, RooAbsData &realdata, RooAbsCollection &snapshot, double poiValue=0.0, int verbose=0)
 Generate asimov dataset from best fit value of nuisance parameters, and fill in snapshot of corresponding global observables. More...
 

Function Documentation

RooAbsData * asimovutils::asimovDatasetNominal ( RooStats::ModelConfig *  mc,
double  poiValue = 0.0,
int  verbose = 0 
)

Generate asimov dataset from nominal value of nuisance parameters.

Definition at line 15 of file AsimovUtils.cc.

References toymcoptutils::SimPdfGenInfo::generateAsimov(), and alignCSCRings::r.

Referenced by Asymptotic::asimovDataset().

15  {
16  RooArgSet poi(*mc->GetParametersOfInterest());
17  RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());
18  r->setConstant(true); r->setVal(poiValue);
19  toymcoptutils::SimPdfGenInfo newToyMC(*mc->GetPdf(), *mc->GetObservables(), false);
20  RooRealVar *weightVar = 0;
21  RooAbsData *asimov = newToyMC.generateAsimov(weightVar);
22  delete weightVar;
23  return asimov;
24 }
toymcoptutils::SimPdfGenInfo::generateAsimov
RooAbsData * generateAsimov(RooRealVar *&weightVar)
Definition: ToyMCSamplerOpt.cc:368
RooAbsData * asimovutils::asimovDatasetWithFit ( RooStats::ModelConfig *  mc,
RooAbsData &  realdata,
RooAbsCollection &  snapshot,
double  poiValue = 0.0,
int  verbose = 0 
)

Generate asimov dataset from best fit value of nuisance parameters, and fill in snapshot of corresponding global observables.

Definition at line 26 of file AsimovUtils.cc.

References a, dtNoiseDBValidation_cfg::cerr, CascadeMinimizer::Constrained, gather_cfg::cout, funct::false, toymcoptutils::SimPdfGenInfo::generateAsimov(), utils::makeNuisancePdf(), match(), parseEventContent::prod, alignCSCRings::r, and utils::setAllConstant().

Referenced by Asymptotic::asimovDataset(), and Combine::run().

26  {
27  RooArgSet poi(*mc->GetParametersOfInterest());
28  RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());
29  r->setConstant(true); r->setVal(poiValue);
30  {
31  CloseCoutSentry sentry(verbose < 3);
32  bool needsFit = false;
33  if (mc->GetNuisanceParameters()) {
34  needsFit = true;
35  } else {
36  // Do we have free parameters anyway that need fitting?
37  std::auto_ptr<RooArgSet> params(mc->GetPdf()->getParameters(realdata));
38  std::auto_ptr<TIterator> iter(params->createIterator());
39  for (RooAbsArg *a = (RooAbsArg *) iter->Next(); a != 0; a = (RooAbsArg *) iter->Next()) {
40  RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);
41  if ( rrv != 0 && rrv->isConstant() == false ) { needsFit = true; break; }
42  }
43  }
44  if (needsFit) {
45  //mc->GetPdf()->fitTo(realdata, RooFit::Minimizer("Minuit2","minimize"), RooFit::Strategy(1), RooFit::Constrain(*mc->GetNuisanceParameters()));
46  const RooCmdArg &constrain = (mc->GetNuisanceParameters() ? RooFit::Constrain(*mc->GetNuisanceParameters()) : RooCmdArg());
47  std::auto_ptr<RooAbsReal> nll(mc->GetPdf()->createNLL(realdata, constrain, RooFit::Extended(mc->GetPdf()->canBeExtended())));
48  CascadeMinimizer minim(*nll, CascadeMinimizer::Constrained);
49  minim.setStrategy(1);
50  minim.minimize(verbose-1);
51  }
52  }
53  if (mc->GetNuisanceParameters() && verbose > 1) {
54  std::cout << "Nuisance parameters after fit for asimov dataset: " << std::endl;
55  mc->GetNuisanceParameters()->Print("V");
56  }
57  toymcoptutils::SimPdfGenInfo newToyMC(*mc->GetPdf(), *mc->GetObservables(), false);
58  RooRealVar *weightVar = 0;
59  RooAbsData *asimov = newToyMC.generateAsimov(weightVar);
60  delete weightVar;
61 
62  // NOW SNAPSHOT THE GLOBAL OBSERVABLES
63  if (mc->GetGlobalObservables() && mc->GetGlobalObservables()->getSize() > 0) {
64  RooArgSet gobs(*mc->GetGlobalObservables());
65 
66  // snapshot data global observables
67  RooArgSet snapGlobalObsData;
68  utils::setAllConstant(gobs, true);
69  gobs.snapshot(snapGlobalObsData);
70 
71  RooArgSet nuis(*mc->GetNuisanceParameters());
72  std::auto_ptr<RooAbsPdf> nuispdf(utils::makeNuisancePdf(*mc));
73  RooProdPdf *prod = dynamic_cast<RooProdPdf *>(nuispdf.get());
74  if (prod == 0) throw std::runtime_error("AsimovUtils: the nuisance pdf is not a RooProdPdf!");
75  std::auto_ptr<TIterator> iter(prod->pdfList().createIterator());
76  for (RooAbsArg *a = (RooAbsArg *) iter->Next(); a != 0; a = (RooAbsArg *) iter->Next()) {
77  RooAbsPdf *cterm = dynamic_cast<RooAbsPdf *>(a);
78  if (!cterm) throw std::logic_error("AsimovUtils: a factor of the nuisance pdf is not a Pdf!");
79  if (!cterm->dependsOn(nuis)) continue; // dummy constraints
80  if (typeid(*cterm) == typeid(RooUniform)) continue;
81  std::auto_ptr<RooArgSet> cpars(cterm->getParameters(&gobs));
82  std::auto_ptr<RooArgSet> cgobs(cterm->getObservables(&gobs));
83  if (cgobs->getSize() != 1) {
84  throw std::runtime_error(Form("AsimovUtils: constraint term %s has multiple global observables", cterm->GetName()));
85  }
86  RooRealVar &rrv = dynamic_cast<RooRealVar &>(*cgobs->first());
87 
88  RooAbsReal *match = 0;
89  if (cpars->getSize() == 1) {
90  match = dynamic_cast<RooAbsReal *>(cpars->first());
91  } else {
92  std::auto_ptr<TIterator> iter2(cpars->createIterator());
93  for (RooAbsArg *a2 = (RooAbsArg *) iter2->Next(); a2 != 0; a2 = (RooAbsArg *) iter2->Next()) {
94  RooRealVar *rrv2 = dynamic_cast<RooRealVar *>(a2);
95  if (rrv2 != 0 && !rrv2->isConstant()) {
96  if (match != 0) throw std::runtime_error(Form("AsimovUtils: constraint term %s has multiple floating params", cterm->GetName()));
97  match = rrv2;
98  }
99  }
100  }
101  if (match == 0) {
102  std::cerr << "ERROR: AsimovUtils: can't find nuisance for constraint term " << cterm->GetName() << std::endl;
103  std::cerr << "Parameters: " << std::endl;
104  cpars->Print("V");
105  std::cerr << "Observables: " << std::endl;
106  cgobs->Print("V");
107  throw std::runtime_error(Form("AsimovUtils: can't find nuisance for constraint term %s", cterm->GetName()));
108  }
109  std::string pdfName(cterm->ClassName());
110  if (pdfName == "RooGaussian" || pdfName == "RooPoisson") {
111  // this is easy
112  rrv.setVal(match->getVal());
113  } else if (pdfName == "RooGamma") {
114  // notation as in http://en.wikipedia.org/wiki/Gamma_distribution
115  // nuisance = x
116  // global obs = kappa ( = observed sideband events + 1)
117  // scaling = theta ( = extrapolation from sideband to signal)
118  // we want to set the global obs to a value for which the current value
119  // of the nuisance is the best fit one.
120  // best fit x = (k-1)*theta ----> k = x/theta + 1
121  RooArgList leaves; cterm->leafNodeServerList(&leaves);
122  std::auto_ptr<TIterator> iter2(leaves.createIterator());
123  RooAbsReal *match2 = 0;
124  for (RooAbsArg *a2 = (RooAbsArg *) iter2->Next(); a2 != 0; a2 = (RooAbsArg *) iter2->Next()) {
125  RooAbsReal *rar = dynamic_cast<RooAbsReal *>(a2);
126  if (rar == 0 || rar == match || rar == &rrv) continue;
127  if (!rar->isConstant()) throw std::runtime_error(Form("AsimovUtils: extra floating parameter %s of RooGamma %s.", rar->GetName(), cterm->GetName()));
128  if (rar->getVal() == 0) continue; // this could be mu
129  if (match2 != 0) throw std::runtime_error(Form("AsimovUtils: extra constant non-zero parameter %s of RooGamma %s.", rar->GetName(), cterm->GetName()));
130  match2 = rar;
131  }
132  if (match2 == 0) throw std::runtime_error(Form("AsimovUtils: could not find the scaling term for RooGamma %s.", cterm->GetName()));
133  //std::cout << " nuisance " << match->GetName() << " = x = " << match->getVal() << std::endl;
134  //std::cout << " scaling param " << match2->GetName() << " = theta = " << match2->getVal() << std::endl;
135  //std::cout << " global obs " << rrv.GetName() << " = kappa = " << rrv.getVal() << std::endl;
136  //std::cout << " new value for global obs " << rrv.GetName() << " = x/theta+1 = " << (match->getVal()/match2->getVal() + 1) << std::endl;
137  rrv.setVal(match->getVal()/match2->getVal() + 1.);
138  } else {
139  throw std::runtime_error(Form("AsimovUtils: can't handle constraint term %s of type %s", cterm->GetName(), cterm->ClassName()));
140  }
141  }
142 
143  // snapshot
144  snapshot.removeAll();
145  utils::setAllConstant(gobs, true);
146  gobs.snapshot(snapshot);
147 
148  // revert things to normal
149  gobs = snapGlobalObsData;
150 
151  if (verbose > 1) {
152  std::cout << "Global observables for data: " << std::endl;
153  snapGlobalObsData.Print("V");
154  std::cout << "Global observables for asimov: " << std::endl;
155  snapshot.Print("V");
156  }
157  }
158 
159  return asimov;
160 }
utils::setAllConstant
bool setAllConstant(const RooAbsCollection &coll, bool constant=true)
set all RooRealVars to constants. return true if at least one changed status
Definition: utils.cc:248
toymcoptutils::SimPdfGenInfo::generateAsimov
RooAbsData * generateAsimov(RooRealVar *&weightVar)
Definition: ToyMCSamplerOpt.cc:368
utils::makeNuisancePdf
RooAbsPdf * makeNuisancePdf(RooStats::ModelConfig &model, const char *name="nuisancePdf")
Definition: utils.cc:200
dtNoiseDBValidation_cfg.cerr
tuple cerr
Definition: dtNoiseDBValidation_cfg.py:22
a
double a
Definition: hdecay.h:121
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
match
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:6
Generated for CMSSW Reference Manual by doxygen 1.8.5