Logo ROOT  
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"
108
109
110#include "TMatrixD.h"
111
112
114
115using namespace RooStats;
116using namespace std;
117
118////////////////////////////////////////////////////////////////////////////////
119
120SPlot::~SPlot()
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.assign(Args);
136
137 fSData = NULL;
138
139}
140
141////////////////////////////////////////////////////////////////////////////////
142
143SPlot::SPlot(const char* name, const char* title):
144 TNamed(name, title)
145{
146 RooArgList Args;
147
148 fSWeightVars.assign(Args);
149
150 fSData = NULL;
151
152}
153
154////////////////////////////////////////////////////////////////////////////////
155///Constructor from a RooDataSet
156///No sWeighted variables are present
157
158SPlot::SPlot(const char* name, const char* title, const RooDataSet &data):
159 TNamed(name, title)
160{
161 RooArgList Args;
162
163 fSWeightVars.assign(Args);
164
165 fSData = (RooDataSet*) &data;
166}
167
168////////////////////////////////////////////////////////////////////////////////
169/// Copy Constructor from another SPlot
170
171SPlot::SPlot(const SPlot &other):
172 TNamed(other)
173{
174 RooArgList Args = (RooArgList) other.GetSWeightVars();
175
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().
198SPlot::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) {
205 fSData = (RooDataSet*) data.Clone(newName);
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.
251Double_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
331Double_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
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{
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.
419void 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}
#define coutI(a)
Definition: RooMsgService.h:30
#define coutW(a)
Definition: RooMsgService.h:32
#define coutE(a)
Definition: RooMsgService.h:33
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
double Double_t
Definition: RtypesCore.h:59
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:364
char name[80]
Definition: TGX11.cxx:110
char * Form(const char *fmt,...)
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:72
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:799
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:536
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:2010
RooArgSet * getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:578
Bool_t isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:363
void attachDataSet(const RooAbsData &set)
Replace server nodes with names matching the dataset variable names with those data set variables,...
Definition: RooAbsArg.cxx:1573
Int_t getSize() const
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents.
virtual RooAbsArg * addClone(const RooAbsArg &var, Bool_t silent=kFALSE)
Add a clone of the specified argument to list.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
const_iterator end() const
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.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
Storage_t::size_type size() const
const_iterator begin() const
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:375
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:1519
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
virtual void setVal(Double_t value)=0
Set the current value of the object. Needs to be overridden by implementations.
void setConstant(Bool_t value=kTRUE)
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:61
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:35
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:27
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:36
virtual TObject * Clone(const char *newname="") const override
Make a clone of an object using the Streamer facility.
Definition: RooDataSet.h:72
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
Bool_t merge(RooDataSet *data1, RooDataSet *data2=0, RooDataSet *data3=0, RooDataSet *data4=0, RooDataSet *data5=0, RooDataSet *data6=0)
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 'data' argset, to the data set.
virtual Double_t weight() const override
Return event weight of current event.
static RooMsgService & instance()
Return reference to singleton instance.
Bool_t isActive(const RooAbsArg *self, RooFit::MsgTopic facility, RooFit::MsgLevel level)
Check if logging is active for given object/topic/RooFit::MsgLevel combination.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
A class to calculate "sWeights" used to create an "sPlot".
Definition: SPlot.h:32
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 GetYieldFromSWeight(const char *sVariable) const
Sum the SWeights for a particular species over all events.
Definition: SPlot.cxx:331
Double_t GetSWeight(Int_t numEvent, const char *sVariable) const
Retrieve an s weight.
Definition: SPlot.cxx:251
RooDataSet * fSData
Definition: SPlot.h:82
RooArgList GetSWeightVars() const
Return a RooArgList containing all paramters that have s weights.
Definition: SPlot.cxx:372
SPlot()
Default constructor.
Definition: SPlot.cxx:130
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
RooDataSet * SetSData(RooDataSet *data)
Set dataset (if not passed in constructor).
Definition: SPlot.cxx:230
RooDataSet * GetSDataSet() const
Retrieve s-weighted data.
Definition: SPlot.cxx:242
RooArgList fSWeightVars
Definition: SPlot.h:78
Double_t GetSumOfEventSWeight(Int_t numEvent) const
Sum the SWeights for a particular event.
Definition: SPlot.cxx:300
void Print(Option_t *name="") const
Print the matrix as a table of elements.
virtual Double_t Determinant() const
Return the matrix determinant.
Definition: TMatrixT.cxx:1362
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual void Clear(Option_t *="")
Definition: TObject.h:115
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:187
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:696
RooCmdArg SumW2Error(Bool_t flag)
RooCmdArg Extended(Bool_t flag=kTRUE)
RooCmdArg PrintEvalErrors(Int_t numErrors)
RooCmdArg PrintLevel(Int_t code)
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
@ InputArguments
Definition: RooGlobalFunc.h:61
Namespace for the RooStats classes.
Definition: Asimov.h:19
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:65
auto * m
Definition: textangle.C:8