Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 gauss("gauss", "gauss", 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(gauss, 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 std::endl;
117
118////////////////////////////////////////////////////////////////////////////////
119
121{
122 if(TestBit(kOwnData) && fSData)
123 delete fSData;
124
125}
126
127////////////////////////////////////////////////////////////////////////////////
128/// Default constructor
129
131{
132 RooArgList Args;
133
134 fSWeightVars.assign(Args);
135}
136
137////////////////////////////////////////////////////////////////////////////////
138
139SPlot::SPlot(const char *name, const char *title) : TNamed(name, title)
140{
141 RooArgList Args;
142 fSWeightVars.assign(Args);
143}
144
145////////////////////////////////////////////////////////////////////////////////
146///Constructor from a RooDataSet
147///No sWeighted variables are present
148
149SPlot::SPlot(const char *name, const char *title, const RooDataSet &data)
150 : TNamed(name, title), fSData(const_cast<RooDataSet *>(&data))
151{
152 RooArgList Args;
153
154 fSWeightVars.assign(Args);
155}
156
157////////////////////////////////////////////////////////////////////////////////
158/// Copy Constructor from another SPlot
159
160SPlot::SPlot(const SPlot &other) : TNamed(other), fSData((RooDataSet *)other.GetSDataSet())
161{
162 RooArgList Args = (RooArgList) other.GetSWeightVars();
163
165}
166
167////////////////////////////////////////////////////////////////////////////////
168///Construct a new SPlot instance, calculate sWeights, and include them
169///in the RooDataSet held by this instance.
170///
171/// The constructor automatically calls AddSWeight() to add s weights to the dataset.
172/// These can be retrieved later using GetSWeight() or GetSDataSet().
173///\param[in] name Name of the instance.
174///\param[in] title Title of the instance.
175///\param[in] data Dataset to fit to.
176///\param[in] pdf PDF to compute s weights for.
177///\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.
178///\param[in] projDeps Don't normalise over these parameters when calculating the sWeights. Will be passed on to AddSWeight().
179///\param[in] useWeights Include weights of the input data in calculation of s weights.
180///\param[in] cloneData Make a clone of the incoming data before adding weights.
181///\param[in] newName New name for the data.
182///\param[in] arg5,arg6,arg7,arg8 Additional arguments for the fitting step in AddSWeight().
183SPlot::SPlot(const char* name, const char* title, RooDataSet& data, RooAbsPdf* pdf,
184 const RooArgList &yieldsList, const RooArgSet &projDeps,
185 bool useWeights, bool cloneData, const char* newName,
186 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8):
187 TNamed(name, title)
188{
189 if(cloneData) {
190 fSData = static_cast<RooDataSet*>(data.Clone(newName));
192 }
193 else
194 fSData = (RooDataSet*) &data;
195
196 // Add check that yieldsList contains all RooRealVar / RooAbsRealLValue
197 for (const auto arg : yieldsList) {
198 if (!dynamic_cast<const RooAbsRealLValue*>(arg)) {
199 coutE(InputArguments) << "SPlot::SPlot(" << GetName() << ") input argument "
200 << arg->GetName() << " is not of type RooRealVar (or RooLinearVar)."
201 << "\nRooStats must be able to set it to 0 and to 1 to probe the PDF." << endl ;
202 throw std::invalid_argument(Form("SPlot::SPlot(%s) input argument %s is not of type RooRealVar/RooLinearVar",GetName(),arg->GetName())) ;
203 }
204 }
205
206 //Construct a new SPlot class,
207 //calculate sWeights, and include them
208 //in the RooDataSet of this class.
209
210 this->AddSWeight(pdf, yieldsList, projDeps, useWeights, arg5, arg6, arg7, arg8);
211}
212
213////////////////////////////////////////////////////////////////////////////////
214/// Set dataset (if not passed in constructor).
216{
217 if(data) {
219 return fSData;
220 } else
221 return nullptr;
222}
223
224////////////////////////////////////////////////////////////////////////////////
225/// Retrieve s-weighted data.
226/// It does **not** automatically call AddSWeight(). This needs to be done manually.
228{
229 return fSData;
230}
231
232////////////////////////////////////////////////////////////////////////////////
233/// Retrieve an s weight.
234/// \param[in] numEvent Event number to retrieve s weight for.
235/// \param[in] sVariable The yield parameter to retrieve the s weight for.
236double SPlot::GetSWeight(Int_t numEvent, const char* sVariable) const
237{
238 if(numEvent > fSData->numEntries() )
239 {
240 coutE(InputArguments) << "Invalid Entry Number" << endl;
241 return -1;
242 }
243
244 if(numEvent < 0)
245 {
246 coutE(InputArguments) << "Invalid Entry Number" << endl;
247 return -1;
248 }
249
250 double totalYield = 0;
251
252 std::string varname(sVariable);
253 varname += "_sw";
254
255
256 if(fSWeightVars.find(sVariable) )
257 {
258 RooArgSet Row(*fSData->get(numEvent));
259 totalYield += Row.getRealValue(sVariable);
260
261 return totalYield;
262 }
263
264 if( fSWeightVars.find(varname.c_str()) )
265 {
266
267 RooArgSet Row(*fSData->get(numEvent));
268 totalYield += Row.getRealValue(varname.c_str() );
269
270 return totalYield;
271 }
272
273 else
274 coutE(InputArguments) << "InputVariable not in list of sWeighted variables" << endl;
275
276 return -1;
277}
278
279
280////////////////////////////////////////////////////////////////////////////////
281/// Sum the SWeights for a particular event.
282/// This sum should equal the total weight of that event.
283/// This method is intended to be used as a check.
284
285double SPlot::GetSumOfEventSWeight(Int_t numEvent) const
286{
287 if(numEvent > fSData->numEntries() )
288 {
289 coutE(InputArguments) << "Invalid Entry Number" << endl;
290 return -1;
291 }
292
293 if(numEvent < 0)
294 {
295 coutE(InputArguments) << "Invalid Entry Number" << endl;
296 return -1;
297 }
298
299 Int_t numSWeightVars = this->GetNumSWeightVars();
300
301 double eventSWeight = 0;
302
303 RooArgSet Row(*fSData->get(numEvent));
304
305 for (Int_t i = 0; i < numSWeightVars; i++)
306 eventSWeight += Row.getRealValue(fSWeightVars.at(i)->GetName() );
307
308 return eventSWeight;
309}
310
311////////////////////////////////////////////////////////////////////////////////
312/// Sum the SWeights for a particular species over all events.
313/// This should equal the total (weighted) yield of that species.
314/// This method is intended as a check.
315
316double SPlot::GetYieldFromSWeight(const char* sVariable) const
317{
318
319 double totalYield = 0;
320
321 std::string varname(sVariable);
322 varname += "_sw";
323
324
325 if(fSWeightVars.find(sVariable) )
326 {
327 for(Int_t i=0; i < fSData->numEntries(); i++)
328 {
329 RooArgSet Row(*fSData->get(i));
330 totalYield += Row.getRealValue(sVariable);
331 }
332
333 return totalYield;
334 }
335
336 if( fSWeightVars.find(varname.c_str()) )
337 {
338 for(Int_t i=0; i < fSData->numEntries(); i++)
339 {
340 RooArgSet Row(*fSData->get(i));
341 totalYield += Row.getRealValue(varname.c_str() );
342 }
343
344 return totalYield;
345 }
346
347 else
348 coutE(InputArguments) << "InputVariable not in list of sWeighted variables" << endl;
349
350 return -1;
351}
352
353
354////////////////////////////////////////////////////////////////////////////////
355/// Return a RooArgList containing all parameters that have s weights.
356
358{
359
361
362 return Args;
363
364}
365
366////////////////////////////////////////////////////////////////////////////////
367/// Return the number of SWeights
368/// In other words, return the number of
369/// species that we are trying to extract.
370
372{
374
375 return Args.size();
376}
377
378////////////////////////////////////////////////////////////////////////////////
379/// Method which adds the sWeights to the dataset.
380///
381/// The SPlot will contain two new variables for each yield parameter:
382/// - `L_<varname>` is the likelihood for each event, i.e., the pdf evaluated for the a given value of the variable "varname".
383/// - `<varname>_sw` is the value of the sWeight for the variable "varname" for each event.
384///
385/// Find Parameters in the PDF to be considered fixed when calculating the SWeights
386/// and be sure to NOT include the yields in that list.
387///
388/// After fixing non-yield parameters, this function will start a fit by calling
389/// ```
390/// pdf->fitTo(*fSData, RooFit::Extended(true), RooFit::SumW2Error(true), RooFit::PrintLevel(-1), RooFit::PrintEvalErrors(-1)).
391/// ```
392/// One can pass additional arguments to `fitTo`, such as `RooFit::Range("fitrange")`, as `arg5`, `arg6`, `arg7`, `arg8`.
393///
394/// \note A `RooFit::Range` may be necessary to get expected results if you initially fit in a range
395/// and/or called `pdf->fixCoefRange("fitrange")` on `pdf`.
396/// Pass `arg5`, `arg6`, `arg7`, `arg8` AT YOUR OWN RISK.
397///
398/// \param[in] pdf PDF to fit to data to compute s weights.
399/// \param[in] yieldsTmp Yields to use to compute s weights.
400/// \param[in] projDeps These will not be normalized over when calculating the sWeights,
401/// and will be considered parameters, not observables.
402/// \param[in] includeWeights Include weights of the input data in calculation of s weights.
403/// \param[in] arg5,arg6,arg7,arg8 Optional additional arguments for the fitting step.
404void SPlot::AddSWeight( RooAbsPdf* pdf, const RooArgList &yieldsTmp,
405 const RooArgSet &projDeps, bool includeWeights,
406 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
407{
408
409 // Find Parameters in the PDF to be considered fixed when calculating the SWeights
410 // and be sure to NOT include the yields in that list
411 std::unique_ptr<RooArgSet> constParameters{pdf->getParameters(fSData)};
412 for (unsigned int i=0; i < constParameters->size(); ++i) {
413 // Need a counting loop since collection is being modified
414 auto& par = *(*constParameters)[i];
415 if (std::any_of(yieldsTmp.begin(), yieldsTmp.end(), [&](const RooAbsArg* yield){ return yield->dependsOn(par); })) {
416 constParameters->remove(par, true, true);
417 --i;
418 }
419 }
420
421
422 // Set these parameters constant and store them so they can later
423 // be set to not constant
424 std::vector<RooAbsRealLValue*> constVarHolder;
425
426 for(std::size_t i = 0; i < constParameters->size(); i++)
427 {
428 RooAbsRealLValue* varTemp = static_cast<RooAbsRealLValue*>((*constParameters)[i] );
429 if(varTemp && varTemp->isConstant() == 0 )
430 {
431 varTemp->setConstant();
432 constVarHolder.push_back(varTemp);
433 }
434 }
435
436 // Fit yields to the data with all other variables held constant
437 // This is necessary because SPlot assumes the yields minimise -Log(likelihood)
438 pdf->fitTo(*fSData, RooFit::Extended(true), RooFit::SumW2Error(true), RooFit::PrintLevel(-1), RooFit::PrintEvalErrors(-1), arg5, arg6, arg7, arg8);
439
440 // The list of variables to normalize over when calculating PDF values.
441 RooArgSet vars(*fSData->get() );
442 vars.remove(projDeps, true, true);
443
444 // Hold the value of the fitted yields
445 std::vector<double> yieldsHolder;
446
447 yieldsHolder.reserve(yieldsTmp.size());
448 for(std::size_t i = 0; i < yieldsTmp.size(); i++) {
449 yieldsHolder.push_back(static_cast<RooAbsReal*>(yieldsTmp.at(i))->getVal(&vars));
450 }
451
452 const Int_t nspec = yieldsTmp.size();
453 RooArgList yields = *static_cast<RooArgList*>(yieldsTmp.snapshot(false));
454
456 coutI(InputArguments) << "Printing Yields" << endl;
457 yields.Print();
458 }
459
460
461 // first calculate the pdf values for all species and all events
462 std::vector<RooAbsRealLValue*> yieldvars ;
463 RooArgSet pdfServers;
464 pdf->treeNodeServerList(&pdfServers);
465
466 std::vector<double> yieldvalues ;
467 for (Int_t k = 0; k < nspec; ++k) {
468 auto thisyield = static_cast<const RooAbsReal*>(yields.at(k)) ;
469 auto yieldinpdf = static_cast<RooAbsRealLValue*>(pdfServers.find(thisyield->GetName()));
470 assert(pdf->dependsOn(*yieldinpdf));
471
472 if (yieldinpdf) {
473 coutI(InputArguments)<< "yield in pdf: " << yieldinpdf->GetName() << " " << thisyield->getVal(&vars) << endl;
474
475 yieldvars.push_back(yieldinpdf) ;
476 yieldvalues.push_back(thisyield->getVal(&vars)) ;
477 }
478 }
479
480 Int_t numevents = fSData->numEntries() ;
481
482
483
484
485 // set all yield to zero
486 for(Int_t m=0; m<nspec; ++m) {
487 auto theVar = static_cast<RooAbsRealLValue*>(yieldvars[m]);
488 theVar->setVal(0) ;
489
490 //Check that range of yields is at least (0,1), and fix otherwise
491 if (theVar->getMin() > 0) {
492 coutE(InputArguments) << "Yield variables need to have a range that includes at least [0, 1]. Minimum for "
493 << theVar->GetName() << " is " << theVar->getMin() << std::endl;
494 if (RooRealVar* realVar = dynamic_cast<RooRealVar*>(theVar)) {
495 coutE(InputArguments) << "Setting min range to 0" << std::endl;
496 realVar->setMin(0);
497 } else {
498 throw std::invalid_argument(std::string("Yield variable ") + theVar->GetName() + " must have a range that includes 0.");
499 }
500 }
501
502 if (theVar->getMax() < 1) {
503 coutW(InputArguments) << "Yield variables need to have a range that includes at least [0, 1]. Maximum for "
504 << theVar->GetName() << " is " << theVar->getMax() << std::endl;
505 if (RooRealVar* realVar = dynamic_cast<RooRealVar*>(theVar)) {
506 coutE(InputArguments) << "Setting max range to 1" << std::endl;
507 realVar->setMax(1);
508 } else {
509 throw std::invalid_argument(std::string("Yield variable ") + theVar->GetName() + " must have a range that includes 1.");
510 }
511 }
512 }
513
514
515 // For every event and for every species,
516 // calculate the value of the component pdf for that specie
517 // by setting the yield of that specie to 1
518 // and all others to 0. Evaluate the pdf for each event
519 // and store the values.
520
521 std::unique_ptr<RooArgSet> pdfvars{pdf->getVariables()};
522 std::vector<std::vector<double> > pdfvalues(numevents,std::vector<double>(nspec,0)) ;
523
524 for (Int_t ievt = 0; ievt <numevents; ievt++)
525 {
526 //WVE: FIX THIS PART, EVALUATION PROGRESS!!
527
528 pdfvars->assign(*fSData->get(ievt));
529
530 for(Int_t k = 0; k < nspec; ++k) {
531 auto theVar = static_cast<RooAbsRealLValue*>(yieldvars[k]);
532
533 // set this yield to 1
534 theVar->setVal( 1 ) ;
535 // evaluate the pdf
536 double f_k = pdf->getVal(&vars) ;
537 pdfvalues[ievt][k] = f_k ;
538 if( !(f_k>1 || f_k<1) )
539 coutW(InputArguments) << "Strange pdf value: " << ievt << " " << k << " " << f_k << std::endl ;
540 theVar->setVal( 0 ) ;
541 }
542 }
543
544 // check that the likelihood normalization is fine
545 std::vector<double> norm(nspec,0) ;
546 for (Int_t ievt = 0; ievt <numevents ; ievt++)
547 {
548 double dnorm(0) ;
549 for(Int_t k=0; k<nspec; ++k) dnorm += yieldvalues[k] * pdfvalues[ievt][k] ;
550 for(Int_t j=0; j<nspec; ++j) norm[j] += pdfvalues[ievt][j]/dnorm ;
551 }
552
553 coutI(Contents) << "likelihood norms: " ;
554
555 for(Int_t k=0; k<nspec; ++k) coutI(Contents) << norm[k] << " " ;
556 coutI(Contents) << std::endl ;
557
558 // Make a TMatrixD to hold the covariance matrix.
559 TMatrixD covInv(nspec, nspec);
560 for (Int_t i = 0; i < nspec; i++) for (Int_t j = 0; j < nspec; j++) covInv(i,j) = 0;
561
562 coutI(Contents) << "Calculating covariance matrix";
563
564
565 // Calculate the inverse covariance matrix, using weights
566 for (Int_t ievt = 0; ievt < numevents; ++ievt)
567 {
568
569 fSData->get(ievt) ;
570
571 // Calculate contribution to the inverse of the covariance
572 // matrix. See BAD 509 V2 eqn. 15
573
574 // Sum for the denominator
575 double dsum(0);
576 for(Int_t k = 0; k < nspec; ++k)
577 dsum += pdfvalues[ievt][k] * yieldvalues[k] ;
578
579 for (Int_t n = 0; n < nspec; ++n) {
580 for(Int_t j=0; j<nspec; ++j)
581 {
582 if (includeWeights) {
583 covInv(n, j) += fSData->weight() * pdfvalues[ievt][n] * pdfvalues[ievt][j] / (dsum * dsum);
584 } else {
585 covInv(n, j) += pdfvalues[ievt][n] * pdfvalues[ievt][j] / (dsum * dsum);
586 }
587 }
588 }
589
590 //ADDED WEIGHT ABOVE
591
592 }
593
594 // Covariance inverse should now be computed!
595
596 // Invert to get the covariance matrix
597 if (covInv.Determinant() <=0)
598 {
599 coutE(Eval) << "SPlot Error: covariance matrix is singular; I can't invert it!" << std::endl;
600 covInv.Print();
601 return;
602 }
603
604 TMatrixD covMatrix(TMatrixD::kInverted,covInv);
605
606 //check cov normalization
607 if (RooMsgService::instance().isActive(this, RooFit::Eval, RooFit::DEBUG)) {
608 coutI(Eval) << "Checking Likelihood normalization: " << std::endl;
609 coutI(Eval) << "Yield of specie Sum of Row in Matrix Norm" << std::endl;
610 for(Int_t k=0; k<nspec; ++k)
611 {
612 double covnorm(0) ;
613 for(Int_t m=0; m<nspec; ++m) covnorm += covInv[k][m]*yieldvalues[m] ;
614 double sumrow(0) ;
615 for(Int_t m = 0; m < nspec; ++m) sumrow += covMatrix[k][m] ;
616 coutI(Eval) << yieldvalues[k] << " " << sumrow << " " << covnorm << endl ;
617 }
618 }
619
620 // calculate for each event the sWeight (BAD 509 V2 eq. 21)
621 coutI(Eval) << "Calculating sWeight" << std::endl;
622 std::vector<RooRealVar*> sweightvec ;
623 std::vector<RooRealVar*> pdfvec ;
624 RooArgSet sweightset ;
625
626 // Create and label the variables
627 // used to store the SWeights
628
630
631 for(Int_t k=0; k<nspec; ++k)
632 {
633 std::string wname = std::string(yieldvars[k]->GetName()) + "_sw";
634 RooRealVar* var = new RooRealVar(wname.c_str(),wname.c_str(),0) ;
635 sweightvec.push_back( var) ;
636 sweightset.add(*var) ;
637 fSWeightVars.add(*var);
638
639 wname = "L_" + std::string(yieldvars[k]->GetName());
640 var = new RooRealVar(wname.c_str(),wname.c_str(),0) ;
641 pdfvec.push_back( var) ;
642 sweightset.add(*var) ;
643 }
644
645 // Create and fill a RooDataSet
646 // with the SWeights
647
648 RooDataSet* sWeightData = new RooDataSet("dataset", "dataset with sWeights", sweightset);
649
650 for(Int_t ievt = 0; ievt < numevents; ++ievt)
651 {
652
653 fSData->get(ievt) ;
654
655 // sum for denominator
656 double dsum(0);
657 for(Int_t k = 0; k < nspec; ++k) dsum += pdfvalues[ievt][k] * yieldvalues[k] ;
658 // covariance weighted pdf for each specief
659 for(Int_t n=0; n<nspec; ++n)
660 {
661 double nsum(0) ;
662 for(Int_t j=0; j<nspec; ++j) nsum += covMatrix(n,j) * pdfvalues[ievt][j] ;
663
664
665 //Add the sWeights here!!
666 //Include weights,
667 //ie events weights are absorbed into sWeight
668
669
670 if(includeWeights) sweightvec[n]->setVal(fSData->weight() * nsum/dsum) ;
671 else sweightvec[n]->setVal( nsum/dsum) ;
672
673 pdfvec[n]->setVal( pdfvalues[ievt][n] ) ;
674
675 if( !(std::abs(nsum/dsum)>=0 ) )
676 {
677 coutE(Contents) << "error: " << nsum/dsum << endl ;
678 return;
679 }
680 }
681
682 sWeightData->add(sweightset) ;
683 }
684
685
686 // Add the SWeights to the original data set
687
688 fSData->merge(sWeightData);
689
690 delete sWeightData;
691
692 //Restore yield values
693
694 for(std::size_t i = 0; i < yieldsTmp.size(); i++)
695 static_cast<RooAbsRealLValue*>(yieldsTmp.at(i))->setVal(yieldsHolder.at(i));
696
697 //Make any variables that were forced to constant no longer constant
698
699 for(Int_t i=0; i < (Int_t) constVarHolder.size(); i++)
700 constVarHolder.at(i)->setConstant(false);
701
702 return;
703
704}
#define coutI(a)
#define coutW(a)
#define coutE(a)
int Int_t
Definition RtypesCore.h:45
#define ClassImp(name)
Definition Rtypes.h:382
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:2489
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
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.
bool isConstant() const
Check if the "Constant" attribute is set.
Definition RooAbsArg.h:304
RooFit::OwningPtr< 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...
RooFit::OwningPtr< RooArgSet > getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expression tree)
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...
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.
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.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooFit::OwningPtr< RooFitResult > fitTo(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Fit PDF to given dataset.
Definition RooAbsPdf.h:157
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
void setConstant(bool value=true)
virtual void setVal(double value)=0
Set the current value of the object. Needs to be overridden by implementations.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
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:24
Named container for two doubles, two integers two object points and three string pointers that can be...
Definition RooCmdArg.h:26
Container class to hold unbinned data.
Definition RooDataSet.h:34
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
bool merge(RooDataSet *data1, RooDataSet *data2=nullptr, RooDataSet *data3=nullptr, RooDataSet *data4=nullptr, RooDataSet *data5=nullptr, RooDataSet *data6=nullptr)
void add(const RooArgSet &row, double weight, double weightError)
Add one ore more rows of data.
double weight() const override
Return event weight of current event.
static RooMsgService & instance()
Return reference to singleton instance.
bool isActive(T self, RooFit::MsgTopic topic, RooFit::MsgLevel level)
Check if logging is active for given object/topic/RooFit::MsgLevel combination.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
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:236
void AddSWeight(RooAbsPdf *pdf, const RooArgList &yieldsTmp, const RooArgSet &projDeps=RooArgSet(), bool includeWeights=true, const RooCmdArg &fitToarg5={}, const RooCmdArg &fitToarg6={}, const RooCmdArg &fitToarg7={}, const RooCmdArg &fitToarg8={})
Method which adds the sWeights to the dataset.
Definition SPlot.cxx:404
RooDataSet * fSData
Definition SPlot.h:82
RooArgList GetSWeightVars() const
Return a RooArgList containing all parameters that have s weights.
Definition SPlot.cxx:357
double GetSumOfEventSWeight(Int_t numEvent) const
Sum the SWeights for a particular event.
Definition SPlot.cxx:285
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:371
~SPlot() override
Definition SPlot.cxx:120
RooDataSet * SetSData(RooDataSet *data)
Set dataset (if not passed in constructor).
Definition SPlot.cxx:215
RooDataSet * GetSDataSet() const
Retrieve s-weighted data.
Definition SPlot.cxx:227
double GetYieldFromSWeight(const char *sVariable) const
Sum the SWeights for a particular species over all events.
Definition SPlot.cxx:316
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.
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:199
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:798
RooCmdArg SumW2Error(bool flag)
RooCmdArg PrintEvalErrors(Int_t numErrors)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg Extended(bool flag=true)
const Int_t n
Definition legend1.C:16
@ InputArguments
Namespace for the RooStats classes.
Definition CodegenImpl.h:58
TMarker m
Definition textangle.C:8