Logo ROOT   master
Reference Guide
SPlot.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer 28/07/2008
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /*****************************************************************************
13  * Project: RooStats
14  * Package: RooFit/RooStats
15  *
16  * Authors:
17  * Original code from M. Pivk as part of MLFit package from BaBar.
18  * Modifications:
19  * Giacinto Piacquadio, Maurizio Pierini: modifications for new RooFit version
20  * George H. Lewis, Kyle Cranmer: generalized for weighted events
21  *
22  * Porting to RooStats (with permission) by Kyle Cranmer, July 2008
23  *
24  *****************************************************************************/
25 
26 
27 /** \class RooStats::SPlot
28  \ingroup Roostats
29 
30  A class to calculate "sWeights" used to create an "sPlot".
31  An sPlot can reweight a dataset to show different components (e.g. signal / background),
32  but it doesn't use cuts, and therefore doesn't have to sort events into signal/background (or other) categories.
33  Instead of *assigning* a category to each event in the dataset, all events are *weighted*.
34  To compute the weights, a PDF with different components is analysed, and the weights are added
35  to the dataset. When plotting the dataset with the weights of the signal or
36  background components, the data looks like "signal", but all events in the dataset are used.
37 
38  The result is similar to a likelihood projection plot, but without cuts.
39 
40  \note SPlot needs to fit the pdf to the data once, so make sure that all relevant fit arguments such as
41  the fit range are passed in the constructor.
42 
43  The code is based on
44  ``SPlot: A statistical tool to unfold data distributions,''
45  Nucl. Instrum. Meth. A 555, 356 (2005) [arXiv:physics/0402083].
46 
47  ### Creating an SPlot
48  To use this class, you first must have a pdf that includes
49  yield parameters for (possibly several) different species, for example a signal and background
50  yield. Those yields must be of type RooRealVar / RooLinearVar (or anything that derives from
51  RooAbsRealLValue). This is necessary because
52  RooStats needs to be able to set the yields to 0 and 1 to probe the PDF. After
53  constructing the s weights, the yields will be restored to their original values.
54 
55  To create an instance of the SPlot, supply a data set, the pdf to analyse,
56  and a list which parameters of the pdf are yields. The SPlot will calculate SWeights, and
57  include these as columns in the RooDataSet. The dataset will have two additional columns
58  for every yield with name "`<varname>`":
59  - `L_<varname>` is the the likelihood for each event, *i.e.*, the pdf evaluated for the given value of the variable "varname".
60  - `<varname>_sw` is the value of the sWeight for the variable "varname" for each event.
61 
62  In SPlot::SPlot(), one can choose whether columns should be added to an existing dataset or whether a copy of the dataset
63  should be created.
64 
65  ### Plotting s-weighted data
66  After computing the s weights, create a new dataset that uses the s weights of the variable of interest for weighting.
67  If the yield parameter for signal was e.g. "signalYield", the dataset can be constructed as follows:
68  ~~~{.cpp}
69  RooDataSet data_signal("<name>", "<title>", <dataWithSWeights>, <variables>, 0, "signalYield_sw");
70  ~~~
71 
72  A complete tutorial with an extensive model is rs301_splot.C
73 
74  #### Using ratios as yield parameters
75  As mentioned, RooStats needs to be able to modify the yield parameters. That means that they have to be a RooRealVar
76  of a RooLinearVar. This allows using ratio parameters as in the following example:
77  ~~~{.cpp}
78  RooRealVar x("x", "observable", 0, 0, 20);
79  RooRealVar m("m", "mean", 5., -10, 10);
80  RooRealVar s("s", "sigma", 2., 0, 10);
81  RooGaussian gaus("gaus", "gaus", x, m, s);
82 
83  RooRealVar a("a", "exp", -0.2, -10., 0.);
84  RooExponential ex("ex", "ex", x, a);
85 
86  RooRealVar common("common", "common scale", 3., 0, 10);
87  RooRealVar r1("r1", "ratio of signal events", 0.3, 0, 10);
88  RooRealVar r2("r2", "ratio of background events", 0.5, 0, 10);
89  RooLinearVar c1("c1", "c1", r1, common, RooFit::RooConst(0.));
90  RooLinearVar c2("c2", "c2", r2, common, RooFit::RooConst(0.));
91 
92  RooAddPdf sum("sum", "sum", RooArgSet(gaus, ex), RooArgSet(c1, c2));
93  auto data = sum.generate(x, 1000);
94 
95  RooStats::SPlot splot("splot", "splot", *data, &sum, RooArgSet(c1, c2));
96  ~~~
97 */
98 
99 #include <vector>
100 #include <map>
101 
102 #include "RooStats/SPlot.h"
103 #include "RooAbsPdf.h"
104 #include "RooDataSet.h"
105 #include "RooRealVar.h"
106 #include "RooGlobalFunc.h"
107 #include "RooStats/RooStatsUtils.h"
108 
109 
110 #include "TMatrixD.h"
111 
112 
114 
115 using namespace RooStats;
116 using namespace std;
117 
118 ////////////////////////////////////////////////////////////////////////////////
119 
121 {
122  if(TestBit(kOwnData) && fSData)
123  delete fSData;
124 
125 }
126 
127 ////////////////////////////////////////////////////////////////////////////////
128 /// Default constructor
129 
131  TNamed()
132 {
133  RooArgList Args;
134 
135  fSWeightVars = Args;
136 
137  fSData = NULL;
138 
139 }
140 
141 ////////////////////////////////////////////////////////////////////////////////
142 
143 SPlot::SPlot(const char* name, const char* title):
144  TNamed(name, title)
145 {
146  RooArgList Args;
147 
148  fSWeightVars = Args;
149 
150  fSData = NULL;
151 
152 }
153 
154 ////////////////////////////////////////////////////////////////////////////////
155 ///Constructor from a RooDataSet
156 ///No sWeighted variables are present
157 
158 SPlot::SPlot(const char* name, const char* title, const RooDataSet &data):
159  TNamed(name, title)
160 {
161  RooArgList Args;
162 
163  fSWeightVars = Args;
164 
165  fSData = (RooDataSet*) &data;
166 }
167 
168 ////////////////////////////////////////////////////////////////////////////////
169 /// Copy Constructor from another SPlot
170 
171 SPlot::SPlot(const SPlot &other):
172  TNamed(other)
173 {
174  RooArgList Args = (RooArgList) other.GetSWeightVars();
175 
176  fSWeightVars.addClone(Args);
177 
178  fSData = (RooDataSet*) other.GetSDataSet();
179 
180 }
181 
182 ////////////////////////////////////////////////////////////////////////////////
183 ///Construct a new SPlot instance, calculate sWeights, and include them
184 ///in the RooDataSet held by this instance.
185 ///
186 /// The constructor automatically calls AddSWeight() to add s weights to the dataset.
187 /// These can be retrieved later using GetSWeight() or GetSDataSet().
188 ///\param[in] name Name of the instance.
189 ///\param[in] title Title of the instance.
190 ///\param[in] data Dataset to fit to.
191 ///\param[in] pdf PDF to compute s weights for.
192 ///\param[in] yieldsList List of parameters in `pdf` that are yields. These must be RooRealVar or RooLinearVar, since RooStats will need to modify their values.
193 ///\param[in] projDeps Don't normalise over these parameters when calculating the sWeights. Will be passed on to AddSWeight().
194 ///\param[in] useWeights Include weights of the input data in calculation of s weights.
195 ///\param[in] cloneData Make a clone of the incoming data before adding weights.
196 ///\param[in] newName New name for the data.
197 ///\param[in] argX Additional arguments for the fitting step in AddSWeight().
198 SPlot::SPlot(const char* name, const char* title, RooDataSet& data, RooAbsPdf* pdf,
199  const RooArgList &yieldsList, const RooArgSet &projDeps,
200  bool useWeights, bool cloneData, const char* newName,
201  const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8):
202  TNamed(name, title)
203 {
204  if(cloneData == 1) {
205  fSData = (RooDataSet*) data.Clone(newName);
206  SetBit(kOwnData);
207  }
208  else
209  fSData = (RooDataSet*) &data;
210 
211  // Add check that yieldsList contains all RooRealVar / RooAbsRealLValue
212  for (const auto arg : yieldsList) {
213  if (!dynamic_cast<const RooAbsRealLValue*>(arg)) {
214  coutE(InputArguments) << "SPlot::SPlot(" << GetName() << ") input argument "
215  << arg->GetName() << " is not of type RooRealVar (or RooLinearVar)."
216  << "\nRooStats must be able to set it to 0 and to 1 to probe the PDF." << endl ;
217  throw std::invalid_argument(Form("SPlot::SPlot(%s) input argument %s is not of type RooRealVar/RooLinearVar",GetName(),arg->GetName())) ;
218  }
219  }
220 
221  //Construct a new SPlot class,
222  //calculate sWeights, and include them
223  //in the RooDataSet of this class.
224 
225  this->AddSWeight(pdf, yieldsList, projDeps, useWeights, arg5, arg6, arg7, arg8);
226 }
227 
228 ////////////////////////////////////////////////////////////////////////////////
229 /// Set dataset (if not passed in constructor).
231 {
232  if(data) {
233  fSData = (RooDataSet*) data;
234  return fSData;
235  } else
236  return NULL;
237 }
238 
239 ////////////////////////////////////////////////////////////////////////////////
240 /// Retrieve s-weighted data.
241 /// It does **not** automatically call AddSWeight(). This needs to be done manually.
243 {
244  return fSData;
245 }
246 
247 ////////////////////////////////////////////////////////////////////////////////
248 /// Retrieve an s weight.
249 /// \param[in] numEvent Event number to retrieve s weight for.
250 /// \param[in] sVariable The yield parameter to retrieve the s weight for.
251 Double_t SPlot::GetSWeight(Int_t numEvent, const char* sVariable) const
252 {
253  if(numEvent > fSData->numEntries() )
254  {
255  coutE(InputArguments) << "Invalid Entry Number" << endl;
256  return -1;
257  }
258 
259  if(numEvent < 0)
260  {
261  coutE(InputArguments) << "Invalid Entry Number" << endl;
262  return -1;
263  }
264 
265  Double_t totalYield = 0;
266 
267  std::string varname(sVariable);
268  varname += "_sw";
269 
270 
271  if(fSWeightVars.find(sVariable) )
272  {
273  RooArgSet Row(*fSData->get(numEvent));
274  totalYield += Row.getRealValue(sVariable);
275 
276  return totalYield;
277  }
278 
279  if( fSWeightVars.find(varname.c_str()) )
280  {
281 
282  RooArgSet Row(*fSData->get(numEvent));
283  totalYield += Row.getRealValue(varname.c_str() );
284 
285  return totalYield;
286  }
287 
288  else
289  coutE(InputArguments) << "InputVariable not in list of sWeighted variables" << endl;
290 
291  return -1;
292 }
293 
294 
295 ////////////////////////////////////////////////////////////////////////////////
296 /// Sum the SWeights for a particular event.
297 /// This sum should equal the total weight of that event.
298 /// This method is intended to be used as a check.
299 
301 {
302  if(numEvent > fSData->numEntries() )
303  {
304  coutE(InputArguments) << "Invalid Entry Number" << endl;
305  return -1;
306  }
307 
308  if(numEvent < 0)
309  {
310  coutE(InputArguments) << "Invalid Entry Number" << endl;
311  return -1;
312  }
313 
314  Int_t numSWeightVars = this->GetNumSWeightVars();
315 
316  Double_t eventSWeight = 0;
317 
318  RooArgSet Row(*fSData->get(numEvent));
319 
320  for (Int_t i = 0; i < numSWeightVars; i++)
321  eventSWeight += Row.getRealValue(fSWeightVars.at(i)->GetName() );
322 
323  return eventSWeight;
324 }
325 
326 ////////////////////////////////////////////////////////////////////////////////
327 /// Sum the SWeights for a particular species over all events.
328 /// This should equal the total (weighted) yield of that species.
329 /// This method is intended as a check.
330 
331 Double_t SPlot::GetYieldFromSWeight(const char* sVariable) const
332 {
333 
334  Double_t totalYield = 0;
335 
336  std::string varname(sVariable);
337  varname += "_sw";
338 
339 
340  if(fSWeightVars.find(sVariable) )
341  {
342  for(Int_t i=0; i < fSData->numEntries(); i++)
343  {
344  RooArgSet Row(*fSData->get(i));
345  totalYield += Row.getRealValue(sVariable);
346  }
347 
348  return totalYield;
349  }
350 
351  if( fSWeightVars.find(varname.c_str()) )
352  {
353  for(Int_t i=0; i < fSData->numEntries(); i++)
354  {
355  RooArgSet Row(*fSData->get(i));
356  totalYield += Row.getRealValue(varname.c_str() );
357  }
358 
359  return totalYield;
360  }
361 
362  else
363  coutE(InputArguments) << "InputVariable not in list of sWeighted variables" << endl;
364 
365  return -1;
366 }
367 
368 
369 ////////////////////////////////////////////////////////////////////////////////
370 /// Return a RooArgList containing all paramters that have s weights.
371 
373 {
374 
375  RooArgList Args = fSWeightVars;
376 
377  return Args;
378 
379 }
380 
381 ////////////////////////////////////////////////////////////////////////////////
382 /// Return the number of SWeights
383 /// In other words, return the number of
384 /// species that we are trying to extract.
385 
387 {
388  RooArgList Args = fSWeightVars;
389 
390  return Args.getSize();
391 }
392 
393 ////////////////////////////////////////////////////////////////////////////////
394 /// Method which adds the sWeights to the dataset.
395 ///
396 /// The SPlot will contain two new variables for each yield parameter:
397 /// - `L_<varname>` is the the likelihood for each event, *i.e.*, the pdf evaluated for the a given value of the variable "varname".
398 /// - `<varname>_sw` is the value of the sWeight for the variable "varname" for each event.
399 ///
400 /// Find Parameters in the PDF to be considered fixed when calculating the SWeights
401 /// and be sure to NOT include the yields in that list.
402 ///
403 /// After fixing non-yield parameters, this function will start a fit by calling
404 /// ```
405 /// pdf->fitTo(*fSData, RooFit::Extended(kTRUE), RooFit::SumW2Error(kTRUE), RooFit::PrintLevel(-1), RooFit::PrintEvalErrors(-1)).
406 /// ```
407 /// One can pass additional arguments to `fitTo`, such as `RooFit::Range("fitrange")`, as `arg5`, `arg6`, `arg7`, `arg8`.
408 ///
409 /// \note A `RooFit::Range` may be necessary to get expected results if you initially fit in a range
410 /// and/or called `pdf->fixCoefRange("fitrange")` on `pdf`.
411 /// Pass `arg5`, `arg6`, `arg7`, `arg8` AT YOUR OWN RISK.
412 ///
413 /// \param[in] pdf PDF to fit to data to compute s weights.
414 /// \param[in] yieldsTmp Yields to use to compute s weights.
415 /// \param[in] projDeps These will not be normalized over when calculating the sWeights,
416 /// and will be considered parameters, not observables.
417 /// \param[in] includeWeights Include weights of the input data in calculation of s weights.
418 /// \param[in] argX Optional additional arguments for the fitting step.
419 void SPlot::AddSWeight( RooAbsPdf* pdf, const RooArgList &yieldsTmp,
420  const RooArgSet &projDeps, bool includeWeights,
421  const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
422 {
423 
424  // Find Parameters in the PDF to be considered fixed when calculating the SWeights
425  // and be sure to NOT include the yields in that list
426  RooArgList* constParameters = (RooArgList*)pdf->getParameters(fSData) ;
427  for (unsigned int i=0; i < constParameters->size(); ++i) {
428  // Need a counting loop since collection is being modified
429  auto& par = (*constParameters)[i];
430  if (std::any_of(yieldsTmp.begin(), yieldsTmp.end(), [&](const RooAbsArg* yield){ return yield->dependsOn(par); })) {
431  constParameters->remove(par, kTRUE, kTRUE);
432  --i;
433  }
434  }
435 
436 
437  // Set these parameters constant and store them so they can later
438  // be set to not constant
439  std::vector<RooAbsRealLValue*> constVarHolder;
440 
441  for(Int_t i = 0; i < constParameters->getSize(); i++)
442  {
443  RooAbsRealLValue* varTemp = static_cast<RooAbsRealLValue*>( constParameters->at(i) );
444  if(varTemp && varTemp->isConstant() == 0 )
445  {
446  varTemp->setConstant();
447  constVarHolder.push_back(varTemp);
448  }
449  }
450 
451  // Fit yields to the data with all other variables held constant
452  // This is necessary because SPlot assumes the yields minimise -Log(likelihood)
454 
455  // Hold the value of the fitted yields
456  std::vector<double> yieldsHolder;
457 
458  for(Int_t i = 0; i < yieldsTmp.getSize(); i++)
459  yieldsHolder.push_back(static_cast<RooAbsReal*>(yieldsTmp.at(i))->getVal());
460 
461  const Int_t nspec = yieldsTmp.getSize();
462  RooArgList yields = *(RooArgList*)yieldsTmp.snapshot(kFALSE);
463 
465  coutI(InputArguments) << "Printing Yields" << endl;
466  yields.Print();
467  }
468 
469  // The list of variables to normalize over when calculating PDF values.
470 
471  RooArgSet vars(*fSData->get() );
472  vars.remove(projDeps, kTRUE, kTRUE);
473 
474  // Attach data set
475 
476  // const_cast<RooAbsPdf*>(pdf)->attachDataSet(*fSData);
477 
478  pdf->attachDataSet(*fSData);
479 
480  // first calculate the pdf values for all species and all events
481  std::vector<RooAbsRealLValue*> yieldvars ;
482  RooArgSet pdfServers;
483  pdf->treeNodeServerList(&pdfServers);
484 
485  std::vector<Double_t> yieldvalues ;
486  for (Int_t k = 0; k < nspec; ++k) {
487  auto thisyield = static_cast<const RooAbsReal*>(yields.at(k)) ;
488  auto yieldinpdf = static_cast<RooAbsRealLValue*>(pdfServers.find(thisyield->GetName()));
489  assert(pdf->dependsOn(*yieldinpdf));
490 
491  if (yieldinpdf) {
492  coutI(InputArguments)<< "yield in pdf: " << yieldinpdf->GetName() << " " << thisyield->getVal() << endl;
493 
494  yieldvars.push_back(yieldinpdf) ;
495  yieldvalues.push_back(thisyield->getVal()) ;
496  }
497  }
498 
499  Int_t numevents = fSData->numEntries() ;
500 
501 
502 
503 
504  // set all yield to zero
505  for(Int_t m=0; m<nspec; ++m) {
506  auto theVar = static_cast<RooAbsRealLValue*>(yieldvars[m]);
507  theVar->setVal(0) ;
508 
509  //Check that range of yields is at least (0,1), and fix otherwise
510  if (theVar->getMin() > 0) {
511  coutE(InputArguments) << "Yield variables need to have a range that includes at least [0, 1]. Minimum for "
512  << theVar->GetName() << " is " << theVar->getMin() << std::endl;
513  if (RooRealVar* realVar = dynamic_cast<RooRealVar*>(theVar)) {
514  coutE(InputArguments) << "Setting min range to 0" << std::endl;
515  realVar->setMin(0);
516  } else {
517  throw std::invalid_argument(std::string("Yield variable ") + theVar->GetName() + " must have a range that includes 0.");
518  }
519  }
520 
521  if (theVar->getMax() < 1) {
522  coutW(InputArguments) << "Yield variables need to have a range that includes at least [0, 1]. Maximum for "
523  << theVar->GetName() << " is " << theVar->getMax() << std::endl;
524  if (RooRealVar* realVar = dynamic_cast<RooRealVar*>(theVar)) {
525  coutE(InputArguments) << "Setting max range to 1" << std::endl;
526  realVar->setMax(1);
527  } else {
528  throw std::invalid_argument(std::string("Yield variable ") + theVar->GetName() + " must have a range that includes 1.");
529  }
530  }
531  }
532 
533 
534  // For every event and for every species,
535  // calculate the value of the component pdf for that specie
536  // by setting the yield of that specie to 1
537  // and all others to 0. Evaluate the pdf for each event
538  // and store the values.
539 
540  RooArgSet * pdfvars = pdf->getVariables();
541  std::vector<std::vector<Double_t> > pdfvalues(numevents,std::vector<Double_t>(nspec,0)) ;
542 
543  for (Int_t ievt = 0; ievt <numevents; ievt++)
544  {
545  //WVE: FIX THIS PART, EVALUATION PROGRESS!!
546 
547  RooStats::SetParameters(fSData->get(ievt), pdfvars);
548 
549  for(Int_t k = 0; k < nspec; ++k) {
550  auto theVar = static_cast<RooAbsRealLValue*>(yieldvars[k]);
551 
552  // set this yield to 1
553  theVar->setVal( 1 ) ;
554  // evaluate the pdf
555  Double_t f_k = pdf->getVal(&vars) ;
556  pdfvalues[ievt][k] = f_k ;
557  if( !(f_k>1 || f_k<1) )
558  coutW(InputArguments) << "Strange pdf value: " << ievt << " " << k << " " << f_k << std::endl ;
559  theVar->setVal( 0 ) ;
560  }
561  }
562  delete pdfvars;
563 
564  // check that the likelihood normalization is fine
565  std::vector<Double_t> norm(nspec,0) ;
566  for (Int_t ievt = 0; ievt <numevents ; ievt++)
567  {
568  Double_t dnorm(0) ;
569  for(Int_t k=0; k<nspec; ++k) dnorm += yieldvalues[k] * pdfvalues[ievt][k] ;
570  for(Int_t j=0; j<nspec; ++j) norm[j] += pdfvalues[ievt][j]/dnorm ;
571  }
572 
573  coutI(Contents) << "likelihood norms: " ;
574 
575  for(Int_t k=0; k<nspec; ++k) coutI(Contents) << norm[k] << " " ;
576  coutI(Contents) << std::endl ;
577 
578  // Make a TMatrixD to hold the covariance matrix.
579  TMatrixD covInv(nspec, nspec);
580  for (Int_t i = 0; i < nspec; i++) for (Int_t j = 0; j < nspec; j++) covInv(i,j) = 0;
581 
582  coutI(Contents) << "Calculating covariance matrix";
583 
584 
585  // Calculate the inverse covariance matrix, using weights
586  for (Int_t ievt = 0; ievt < numevents; ++ievt)
587  {
588 
589  fSData->get(ievt) ;
590 
591  // Calculate contribution to the inverse of the covariance
592  // matrix. See BAD 509 V2 eqn. 15
593 
594  // Sum for the denominator
595  Double_t dsum(0);
596  for(Int_t k = 0; k < nspec; ++k)
597  dsum += pdfvalues[ievt][k] * yieldvalues[k] ;
598 
599  for(Int_t n=0; n<nspec; ++n)
600  for(Int_t j=0; j<nspec; ++j)
601  {
602  if(includeWeights)
603  covInv(n,j) += fSData->weight()*pdfvalues[ievt][n]*pdfvalues[ievt][j]/(dsum*dsum) ;
604  else
605  covInv(n,j) += pdfvalues[ievt][n]*pdfvalues[ievt][j]/(dsum*dsum) ;
606  }
607 
608  //ADDED WEIGHT ABOVE
609 
610  }
611 
612  // Covariance inverse should now be computed!
613 
614  // Invert to get the covariance matrix
615  if (covInv.Determinant() <=0)
616  {
617  coutE(Eval) << "SPlot Error: covariance matrix is singular; I can't invert it!" << std::endl;
618  covInv.Print();
619  return;
620  }
621 
622  TMatrixD covMatrix(TMatrixD::kInverted,covInv);
623 
624  //check cov normalization
625  if (RooMsgService::instance().isActive(this, RooFit::Eval, RooFit::DEBUG)) {
626  coutI(Eval) << "Checking Likelihood normalization: " << std::endl;
627  coutI(Eval) << "Yield of specie Sum of Row in Matrix Norm" << std::endl;
628  for(Int_t k=0; k<nspec; ++k)
629  {
630  Double_t covnorm(0) ;
631  for(Int_t m=0; m<nspec; ++m) covnorm += covInv[k][m]*yieldvalues[m] ;
632  Double_t sumrow(0) ;
633  for(Int_t m = 0; m < nspec; ++m) sumrow += covMatrix[k][m] ;
634  coutI(Eval) << yieldvalues[k] << " " << sumrow << " " << covnorm << endl ;
635  }
636  }
637 
638  // calculate for each event the sWeight (BAD 509 V2 eq. 21)
639  coutI(Eval) << "Calculating sWeight" << std::endl;
640  std::vector<RooRealVar*> sweightvec ;
641  std::vector<RooRealVar*> pdfvec ;
642  RooArgSet sweightset ;
643 
644  // Create and label the variables
645  // used to store the SWeights
646 
648 
649  for(Int_t k=0; k<nspec; ++k)
650  {
651  std::string wname = std::string(yieldvars[k]->GetName()) + "_sw";
652  RooRealVar* var = new RooRealVar(wname.c_str(),wname.c_str(),0) ;
653  sweightvec.push_back( var) ;
654  sweightset.add(*var) ;
655  fSWeightVars.add(*var);
656 
657  wname = "L_" + std::string(yieldvars[k]->GetName());
658  var = new RooRealVar(wname.c_str(),wname.c_str(),0) ;
659  pdfvec.push_back( var) ;
660  sweightset.add(*var) ;
661  }
662 
663  // Create and fill a RooDataSet
664  // with the SWeights
665 
666  RooDataSet* sWeightData = new RooDataSet("dataset", "dataset with sWeights", sweightset);
667 
668  for(Int_t ievt = 0; ievt < numevents; ++ievt)
669  {
670 
671  fSData->get(ievt) ;
672 
673  // sum for denominator
674  Double_t dsum(0);
675  for(Int_t k = 0; k < nspec; ++k) dsum += pdfvalues[ievt][k] * yieldvalues[k] ;
676  // covariance weighted pdf for each specief
677  for(Int_t n=0; n<nspec; ++n)
678  {
679  Double_t nsum(0) ;
680  for(Int_t j=0; j<nspec; ++j) nsum += covMatrix(n,j) * pdfvalues[ievt][j] ;
681 
682 
683  //Add the sWeights here!!
684  //Include weights,
685  //ie events weights are absorbed into sWeight
686 
687 
688  if(includeWeights) sweightvec[n]->setVal(fSData->weight() * nsum/dsum) ;
689  else sweightvec[n]->setVal( nsum/dsum) ;
690 
691  pdfvec[n]->setVal( pdfvalues[ievt][n] ) ;
692 
693  if( !(fabs(nsum/dsum)>=0 ) )
694  {
695  coutE(Contents) << "error: " << nsum/dsum << endl ;
696  return;
697  }
698  }
699 
700  sWeightData->add(sweightset) ;
701  }
702 
703 
704  // Add the SWeights to the original data set
705 
706  fSData->merge(sWeightData);
707 
708  delete sWeightData;
709 
710  //Restore yield values
711 
712  for(Int_t i = 0; i < yieldsTmp.getSize(); i++)
713  static_cast<RooAbsRealLValue*>(yieldsTmp.at(i))->setVal(yieldsHolder.at(i));
714 
715  //Make any variables that were forced to constant no longer constant
716 
717  for(Int_t i=0; i < (Int_t) constVarHolder.size(); i++)
718  constVarHolder.at(i)->setConstant(kFALSE);
719 
720  return;
721 
722 }
RooCmdArg SumW2Error(Bool_t flag)
virtual void Clear(Option_t *="")
Definition: TObject.h:100
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:1910
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event &#39;index&#39;.
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
#define coutE(a)
Definition: RooMsgService.h:33
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
void treeNodeServerList(RooAbsCollection *list, const RooAbsArg *arg=0, Bool_t doBranch=kTRUE, Bool_t doLeaf=kTRUE, Bool_t valueOnly=kFALSE, Bool_t recurseNonDerived=kFALSE) const
Fill supplied list with nodes of the arg tree, following all server links, starting with ourself as t...
Definition: RooAbsArg.cxx:503
Storage_t::size_type size() const
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
auto * m
Definition: textangle.C:8
Bool_t dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=0, Bool_t valueOnly=kFALSE) const
Test whether we depend on (ie, are served by) any object in the specified collection.
Definition: RooAbsArg.cxx:731
virtual RooAbsArg * addClone(const RooAbsArg &var, Bool_t silent=kFALSE)
Add a clone of the specified argument to list.
#define coutI(a)
Definition: RooMsgService.h:30
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
int Int_t
Definition: CPyCppyy.h:43
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:65
RooDataSet * GetSDataSet() const
Retrieve s-weighted data.
Definition: SPlot.cxx:242
const_iterator begin() const
const_iterator end() const
static RooMsgService & instance()
Return reference to singleton instance.
STL namespace.
#define coutW(a)
Definition: RooMsgService.h:32
virtual Double_t Determinant() const
Return the matrix determinant.
Definition: TMatrixT.cxx:1364
void attachDataSet(const RooAbsData &set)
Replace server nodes with names matching the dataset variable names with those data set variables...
Definition: RooAbsArg.cxx:1468
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:695
Double_t GetYieldFromSWeight(const char *sVariable) const
Sum the SWeights for a particular species over all events.
Definition: SPlot.cxx:331
RooCmdArg PrintLevel(Int_t code)
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
RooCmdArg PrintEvalErrors(Int_t numErrors)
Double_t GetSWeight(Int_t numEvent, const char *sVariable) const
Retrieve an s weight.
Definition: SPlot.cxx:251
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
RooDataSet * fSData
Definition: SPlot.h:82
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual TObject * Clone(const char *newname="") const override
Make a clone of an object using the Streamer facility.
Definition: RooDataSet.h:68
Int_t getSize() const
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents.
virtual Double_t weight() const override
Return event weight of current event.
Definition: RooDataSet.cxx:971
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
void setConstant(Bool_t value=kTRUE)
void AddSWeight(RooAbsPdf *pdf, const RooArgList &yieldsTmp, const RooArgSet &projDeps=RooArgSet(), bool includeWeights=kTRUE, const RooCmdArg &fitToarg5=RooCmdArg::none(), const RooCmdArg &fitToarg6=RooCmdArg::none(), const RooCmdArg &fitToarg7=RooCmdArg::none(), const RooCmdArg &fitToarg8=RooCmdArg::none())
Method which adds the sWeights to the dataset.
Definition: SPlot.cxx:419
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
char * Form(const char *fmt,...)
RooArgList GetSWeightVars() const
Return a RooArgList containing all paramters that have s weights.
Definition: SPlot.cxx:372
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
Int_t GetNumSWeightVars() const
Return the number of SWeights In other words, return the number of species that we are trying to extr...
Definition: SPlot.cxx:386
virtual void setVal(Double_t value)=0
const Bool_t kFALSE
Definition: RtypesCore.h:90
Namespace for the RooStats classes.
Definition: Asimov.h:19
#define ClassImp(name)
Definition: Rtypes.h:361
RooAbsArg * find(const char *name) const
Find object with given name in list.
double Double_t
Definition: RtypesCore.h:57
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don&#39;t match any of...
Definition: RooAbsArg.cxx:545
RooCmdArg Extended(Bool_t flag=kTRUE)
void Print(Option_t *name="") const
Print the matrix as a table of elements.
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
Bool_t merge(RooDataSet *data1, RooDataSet *data2=0, RooDataSet *data3=0, RooDataSet *data4=0, RooDataSet *data5=0, RooDataSet *data6=0)
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Definition: RooArgSet.cxx:472
RooDataSet * SetSData(RooDataSet *data)
Set dataset (if not passed in constructor).
Definition: SPlot.cxx:230
SPlot()
Default constructor.
Definition: SPlot.cxx:130
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1259
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
RooArgList fSWeightVars
Definition: SPlot.h:78
Double_t GetSumOfEventSWeight(Int_t numEvent) const
Sum the SWeights for a particular event.
Definition: SPlot.cxx:300
const Bool_t kTRUE
Definition: RtypesCore.h:89
A class to calculate "sWeights" used to create an "sPlot".
Definition: SPlot.h:32
const Int_t n
Definition: legend1.C:16
char name[80]
Definition: TGX11.cxx:109
Bool_t isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:360
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:28