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 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 = nullptr;
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 = nullptr;
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] arg5,arg6,arg7,arg8 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) {
234 return fSData;
235 } else
236 return nullptr;
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 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 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
300double SPlot::GetSumOfEventSWeight(Int_t numEvent) const
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 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 SPlot::GetYieldFromSWeight(const char* sVariable) const
332{
333
334 double 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 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(true), RooFit::SumW2Error(true), 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] arg5,arg6,arg7,arg8 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, true, true);
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)
453 pdf->fitTo(*fSData, RooFit::Extended(true), RooFit::SumW2Error(true), RooFit::PrintLevel(-1), RooFit::PrintEvalErrors(-1), arg5, arg6, arg7, arg8);
454
455 // The list of variables to normalize over when calculating PDF values.
456 RooArgSet vars(*fSData->get() );
457 vars.remove(projDeps, true, true);
458
459 // Hold the value of the fitted yields
460 std::vector<double> yieldsHolder;
461
462 for(Int_t i = 0; i < yieldsTmp.getSize(); i++)
463 yieldsHolder.push_back(static_cast<RooAbsReal*>(yieldsTmp.at(i))->getVal(&vars));
464
465 const Int_t nspec = yieldsTmp.getSize();
466 RooArgList yields = *(RooArgList*)yieldsTmp.snapshot(false);
467
469 coutI(InputArguments) << "Printing Yields" << endl;
470 yields.Print();
471 }
472
473
474 // first calculate the pdf values for all species and all events
475 std::vector<RooAbsRealLValue*> yieldvars ;
476 RooArgSet pdfServers;
477 pdf->treeNodeServerList(&pdfServers);
478
479 std::vector<double> yieldvalues ;
480 for (Int_t k = 0; k < nspec; ++k) {
481 auto thisyield = static_cast<const RooAbsReal*>(yields.at(k)) ;
482 auto yieldinpdf = static_cast<RooAbsRealLValue*>(pdfServers.find(thisyield->GetName()));
483 assert(pdf->dependsOn(*yieldinpdf));
484
485 if (yieldinpdf) {
486 coutI(InputArguments)<< "yield in pdf: " << yieldinpdf->GetName() << " " << thisyield->getVal(&vars) << endl;
487
488 yieldvars.push_back(yieldinpdf) ;
489 yieldvalues.push_back(thisyield->getVal(&vars)) ;
490 }
491 }
492
493 Int_t numevents = fSData->numEntries() ;
494
495
496
497
498 // set all yield to zero
499 for(Int_t m=0; m<nspec; ++m) {
500 auto theVar = static_cast<RooAbsRealLValue*>(yieldvars[m]);
501 theVar->setVal(0) ;
502
503 //Check that range of yields is at least (0,1), and fix otherwise
504 if (theVar->getMin() > 0) {
505 coutE(InputArguments) << "Yield variables need to have a range that includes at least [0, 1]. Minimum for "
506 << theVar->GetName() << " is " << theVar->getMin() << std::endl;
507 if (RooRealVar* realVar = dynamic_cast<RooRealVar*>(theVar)) {
508 coutE(InputArguments) << "Setting min range to 0" << std::endl;
509 realVar->setMin(0);
510 } else {
511 throw std::invalid_argument(std::string("Yield variable ") + theVar->GetName() + " must have a range that includes 0.");
512 }
513 }
514
515 if (theVar->getMax() < 1) {
516 coutW(InputArguments) << "Yield variables need to have a range that includes at least [0, 1]. Maximum for "
517 << theVar->GetName() << " is " << theVar->getMax() << std::endl;
518 if (RooRealVar* realVar = dynamic_cast<RooRealVar*>(theVar)) {
519 coutE(InputArguments) << "Setting max range to 1" << std::endl;
520 realVar->setMax(1);
521 } else {
522 throw std::invalid_argument(std::string("Yield variable ") + theVar->GetName() + " must have a range that includes 1.");
523 }
524 }
525 }
526
527
528 // For every event and for every species,
529 // calculate the value of the component pdf for that specie
530 // by setting the yield of that specie to 1
531 // and all others to 0. Evaluate the pdf for each event
532 // and store the values.
533
534 RooArgSet * pdfvars = pdf->getVariables();
535 std::vector<std::vector<double> > pdfvalues(numevents,std::vector<double>(nspec,0)) ;
536
537 for (Int_t ievt = 0; ievt <numevents; ievt++)
538 {
539 //WVE: FIX THIS PART, EVALUATION PROGRESS!!
540
541 pdfvars->assign(*fSData->get(ievt));
542
543 for(Int_t k = 0; k < nspec; ++k) {
544 auto theVar = static_cast<RooAbsRealLValue*>(yieldvars[k]);
545
546 // set this yield to 1
547 theVar->setVal( 1 ) ;
548 // evaluate the pdf
549 double f_k = pdf->getVal(&vars) ;
550 pdfvalues[ievt][k] = f_k ;
551 if( !(f_k>1 || f_k<1) )
552 coutW(InputArguments) << "Strange pdf value: " << ievt << " " << k << " " << f_k << std::endl ;
553 theVar->setVal( 0 ) ;
554 }
555 }
556 delete pdfvars;
557
558 // check that the likelihood normalization is fine
559 std::vector<double> norm(nspec,0) ;
560 for (Int_t ievt = 0; ievt <numevents ; ievt++)
561 {
562 double dnorm(0) ;
563 for(Int_t k=0; k<nspec; ++k) dnorm += yieldvalues[k] * pdfvalues[ievt][k] ;
564 for(Int_t j=0; j<nspec; ++j) norm[j] += pdfvalues[ievt][j]/dnorm ;
565 }
566
567 coutI(Contents) << "likelihood norms: " ;
568
569 for(Int_t k=0; k<nspec; ++k) coutI(Contents) << norm[k] << " " ;
570 coutI(Contents) << std::endl ;
571
572 // Make a TMatrixD to hold the covariance matrix.
573 TMatrixD covInv(nspec, nspec);
574 for (Int_t i = 0; i < nspec; i++) for (Int_t j = 0; j < nspec; j++) covInv(i,j) = 0;
575
576 coutI(Contents) << "Calculating covariance matrix";
577
578
579 // Calculate the inverse covariance matrix, using weights
580 for (Int_t ievt = 0; ievt < numevents; ++ievt)
581 {
582
583 fSData->get(ievt) ;
584
585 // Calculate contribution to the inverse of the covariance
586 // matrix. See BAD 509 V2 eqn. 15
587
588 // Sum for the denominator
589 double dsum(0);
590 for(Int_t k = 0; k < nspec; ++k)
591 dsum += pdfvalues[ievt][k] * yieldvalues[k] ;
592
593 for(Int_t n=0; n<nspec; ++n)
594 for(Int_t j=0; j<nspec; ++j)
595 {
596 if(includeWeights)
597 covInv(n,j) += fSData->weight()*pdfvalues[ievt][n]*pdfvalues[ievt][j]/(dsum*dsum) ;
598 else
599 covInv(n,j) += pdfvalues[ievt][n]*pdfvalues[ievt][j]/(dsum*dsum) ;
600 }
601
602 //ADDED WEIGHT ABOVE
603
604 }
605
606 // Covariance inverse should now be computed!
607
608 // Invert to get the covariance matrix
609 if (covInv.Determinant() <=0)
610 {
611 coutE(Eval) << "SPlot Error: covariance matrix is singular; I can't invert it!" << std::endl;
612 covInv.Print();
613 return;
614 }
615
616 TMatrixD covMatrix(TMatrixD::kInverted,covInv);
617
618 //check cov normalization
619 if (RooMsgService::instance().isActive(this, RooFit::Eval, RooFit::DEBUG)) {
620 coutI(Eval) << "Checking Likelihood normalization: " << std::endl;
621 coutI(Eval) << "Yield of specie Sum of Row in Matrix Norm" << std::endl;
622 for(Int_t k=0; k<nspec; ++k)
623 {
624 double covnorm(0) ;
625 for(Int_t m=0; m<nspec; ++m) covnorm += covInv[k][m]*yieldvalues[m] ;
626 double sumrow(0) ;
627 for(Int_t m = 0; m < nspec; ++m) sumrow += covMatrix[k][m] ;
628 coutI(Eval) << yieldvalues[k] << " " << sumrow << " " << covnorm << endl ;
629 }
630 }
631
632 // calculate for each event the sWeight (BAD 509 V2 eq. 21)
633 coutI(Eval) << "Calculating sWeight" << std::endl;
634 std::vector<RooRealVar*> sweightvec ;
635 std::vector<RooRealVar*> pdfvec ;
636 RooArgSet sweightset ;
637
638 // Create and label the variables
639 // used to store the SWeights
640
642
643 for(Int_t k=0; k<nspec; ++k)
644 {
645 std::string wname = std::string(yieldvars[k]->GetName()) + "_sw";
646 RooRealVar* var = new RooRealVar(wname.c_str(),wname.c_str(),0) ;
647 sweightvec.push_back( var) ;
648 sweightset.add(*var) ;
649 fSWeightVars.add(*var);
650
651 wname = "L_" + std::string(yieldvars[k]->GetName());
652 var = new RooRealVar(wname.c_str(),wname.c_str(),0) ;
653 pdfvec.push_back( var) ;
654 sweightset.add(*var) ;
655 }
656
657 // Create and fill a RooDataSet
658 // with the SWeights
659
660 RooDataSet* sWeightData = new RooDataSet("dataset", "dataset with sWeights", sweightset);
661
662 for(Int_t ievt = 0; ievt < numevents; ++ievt)
663 {
664
665 fSData->get(ievt) ;
666
667 // sum for denominator
668 double dsum(0);
669 for(Int_t k = 0; k < nspec; ++k) dsum += pdfvalues[ievt][k] * yieldvalues[k] ;
670 // covariance weighted pdf for each specief
671 for(Int_t n=0; n<nspec; ++n)
672 {
673 double nsum(0) ;
674 for(Int_t j=0; j<nspec; ++j) nsum += covMatrix(n,j) * pdfvalues[ievt][j] ;
675
676
677 //Add the sWeights here!!
678 //Include weights,
679 //ie events weights are absorbed into sWeight
680
681
682 if(includeWeights) sweightvec[n]->setVal(fSData->weight() * nsum/dsum) ;
683 else sweightvec[n]->setVal( nsum/dsum) ;
684
685 pdfvec[n]->setVal( pdfvalues[ievt][n] ) ;
686
687 if( !(std::abs(nsum/dsum)>=0 ) )
688 {
689 coutE(Contents) << "error: " << nsum/dsum << endl ;
690 return;
691 }
692 }
693
694 sWeightData->add(sweightset) ;
695 }
696
697
698 // Add the SWeights to the original data set
699
700 fSData->merge(sWeightData);
701
702 delete sWeightData;
703
704 //Restore yield values
705
706 for(Int_t i = 0; i < yieldsTmp.getSize(); i++)
707 static_cast<RooAbsRealLValue*>(yieldsTmp.at(i))->setVal(yieldsHolder.at(i));
708
709 //Make any variables that were forced to constant no longer constant
710
711 for(Int_t i=0; i < (Int_t) constVarHolder.size(); i++)
712 constVarHolder.at(i)->setConstant(false);
713
714 return;
715
716}
#define coutI(a)
Definition: RooMsgService.h:34
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
int Int_t
Definition: RtypesCore.h:45
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
char name[80]
Definition: TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2468
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:72
bool dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=nullptr, bool valueOnly=false) const
Test whether we depend on (ie, are served by) any object in the specified collection.
Definition: RooAbsArg.cxx:827
bool isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:365
RooArgSet * getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:2083
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:563
void treeNodeServerList(RooAbsCollection *list, const RooAbsArg *arg=nullptr, bool doBranch=true, bool doLeaf=true, bool valueOnly=false, bool recurseNonDerived=false) const
Fill supplied list with nodes of the arg tree, following all server links, starting with ourself as t...
Definition: RooAbsArg.cxx:521
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
RooAbsCollection * snapshot(bool deepCopy=true) const
Take a snap shot of current collection contents.
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
const_iterator end() const
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
virtual RooAbsArg * addClone(const RooAbsArg &var, bool silent=false)
Add a clone of the specified argument to list.
const_iterator begin() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:374
virtual RooFitResult * fitTo(RooAbsData &data, const RooLinkedList &cmdList={})
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1611
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
void setConstant(bool value=true)
virtual void setVal(double value)=0
Set the current value of the object. Needs to be overridden by implementations.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:104
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:56
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:26
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
void add(const RooArgSet &row, double weight=1.0, double weightError=0.0) override
Add one ore more rows of data.
bool merge(RooDataSet *data1, RooDataSet *data2=nullptr, RooDataSet *data3=nullptr, RooDataSet *data4=nullptr, RooDataSet *data5=nullptr, RooDataSet *data6=nullptr)
double weight() const override
Return event weight of current event.
Definition: RooDataSet.cxx:953
static RooMsgService & instance()
Return reference to singleton instance.
bool 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:40
A class to calculate "sWeights" used to create an "sPlot".
Definition: SPlot.h:32
double GetSWeight(Int_t numEvent, const char *sVariable) const
Retrieve an s weight.
Definition: SPlot.cxx:251
RooDataSet * fSData
Definition: SPlot.h:82
void AddSWeight(RooAbsPdf *pdf, const RooArgList &yieldsTmp, const RooArgSet &projDeps=RooArgSet(), bool includeWeights=true, 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
RooArgList GetSWeightVars() const
Return a RooArgList containing all paramters that have s weights.
Definition: SPlot.cxx:372
double GetSumOfEventSWeight(Int_t numEvent) const
Sum the SWeights for a particular event.
Definition: SPlot.cxx:300
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
double GetYieldFromSWeight(const char *sVariable) const
Sum the SWeights for a particular species over all events.
Definition: SPlot.cxx:331
RooArgList fSWeightVars
Definition: SPlot.h:78
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
Double_t Determinant() const override
Return the matrix determinant.
Definition: TMatrixT.cxx:1362
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
virtual void Clear(Option_t *="")
Definition: TObject.h:119
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:201
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:774
RooCmdArg SumW2Error(bool flag)
RooCmdArg PrintEvalErrors(Int_t numErrors)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg Extended(bool flag=true)
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1778
const Int_t n
Definition: legend1.C:16
@ InputArguments
Definition: RooGlobalFunc.h:63
Namespace for the RooStats classes.
Definition: Asimov.h:19
TMarker m
Definition: textangle.C:8