laura is hosted by Hepforge, IPPP Durham
Laura++  v2r1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauSPlot.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2006 - 2013.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
19 /**********************************************************************
20  * *
21  * Copyright (c) 2005 ROOT Foundation, CERN/PH-SFT *
22  * *
23  **********************************************************************/
24 
25 #include <cfloat>
26 #include <iostream>
27 #include <vector>
28 using std::cout;
29 using std::cerr;
30 using std::endl;
31 using std::vector;
32 
33 #include "TEventList.h"
34 #include "TFile.h"
35 #include "TLeaf.h"
36 #include "TMath.h"
37 #include "TObjArray.h"
38 #include "TSystem.h"
39 #include "TTree.h"
40 #include "TVirtualFitter.h"
41 
42 #include "LauSPlot.hh"
43 
44 extern void Yields(Int_t &, Double_t *, Double_t &f, Double_t *x, Int_t iflag);
45 
47 
48 
49 LauSPlot::LauSPlot(const TString& fileName, const TString& treeName,
50  Int_t firstExpt,
51  Int_t nExpt,
52  const NameSet& variableNames,
53  const NumbMap& freeSpecies,
54  const NumbMap& fixdSpecies,
55  const TwoDMap& twodimPDFs,
56  Bool_t sigSplit,
57  Bool_t scfDPSmeared) :
58  fileName_(fileName),
59  inputTreeName_(treeName),
60  cnTreeName_(""),
61  sweightTreeName_(""),
62  file_(0),
63  inputTree_(0),
64  cnTree_(0),
65  sweightTree_(0),
66  eventList_(0),
67  variableNames_(variableNames),
68  freeSpecies_(freeSpecies),
69  fixdSpecies_(fixdSpecies),
70  origFreeSpecies_(freeSpecies),
71  origFixdSpecies_(fixdSpecies),
72  twodimPDFs_(twodimPDFs),
73  signalSplit_(sigSplit),
74  scfDPSmear_(scfDPSmeared),
75  readInput_(kFALSE),
76  definedCNBranches_(kFALSE),
77  definedSWeightBranches_(kFALSE),
78  firstExpt_(firstExpt),
79  nExpt_(nExpt),
80  iExpt_(0),
81  nEvents_(0),
82  nDiscVars_(variableNames.size()),
83  nFreeSpecies_(freeSpecies.size()),
84  nFixdSpecies_(fixdSpecies.size()),
85  nSpecies_(freeSpecies.size()+fixdSpecies.size())
86 {
87  this->openInputFileAndTree();
88  this->readInputInfo();
89  this->createSWeightTree();
90  if (nFixdSpecies_>0) {
91  this->createCNTree();
92  }
93 }
94 
96 {
97  // seems that closing the file deletes the tree
98  // so only delete if the file is still open
99  if (file_ && file_->IsOpen()) {
100  delete inputTree_; inputTree_ = 0;
101  delete sweightTree_; sweightTree_ = 0;
102  delete cnTree_; cnTree_ = 0;
103  }
104  delete file_; file_ = 0;
105 }
106 
108 {
109  // first check whether we've already opened up the file or not
110  if (!file_) {
111  // if not, first check the filename and if all ok create the file
112  if (fileName_ == "") {
113  cerr<<"ERROR in LauSPlot::createFileAndTree : Bad filename supplied, not creating file or tree."<<endl;
114  return;
115  }
116  file_ = TFile::Open(fileName_, "update");
117  if (!file_ || file_->IsZombie() || !file_->IsWritable()) {
118  cerr<<"ERROR in LauSPlot::createFileAndTree : Problem opening file \""<<fileName_<<"\" for updating, can't do anything."<<endl;
119  return;
120  }
121  }
122  // next open the input tree for reading
123  if (!inputTree_) {
124  file_->cd();
125  inputTree_ = dynamic_cast<TTree*>(file_->Get(inputTreeName_));
126  inputTree_->SetDirectory(file_);
127  }
128 }
129 
131 {
132  // Read the tree and then setup the maps so we know which leaves to
133  // read from the tree to get the various PDF values
134 
135  Bool_t inputOK = this->readInputLeaves();
136  if (!inputOK) {
137  this->readInput(inputOK);
138  return;
139  }
140 
141  inputOK = this->checkLeaves();
142  this->readInput(inputOK);
143  return;
144 }
145 
147 {
148  // Read all the leaves in the tree and store them in the leaves map
149 
150  if (!inputTree_) {
151  cerr<<"ERROR in LauSPlot::readInputInfo : Invalid pointer to data tree."<<endl;
152  return kFALSE;
153  }
154 
155  Int_t nBranches = inputTree_ ? static_cast<Int_t>(inputTree_->GetNbranches()) : 0;
156  TObjArray* pLeaves = inputTree_->GetListOfLeaves();
157  if (!pLeaves) {
158  cerr<<"ERROR in LauSPlot::readInputInfo : Problem retrieving leaves from the tree."<<endl;
159  return kFALSE;
160  }
161  TObjArray& leaves = *pLeaves;
162 
163  if (nBranches > leaves.GetSize()) {
164  cerr<<"ERROR in LauSPlot::readInputInfo : List of leaves is smaller than number of branches - this is strange!"<<endl;
165  return kFALSE;
166  }
167 
168  for (Int_t iLeaf = 0; iLeaf < nBranches; ++iLeaf) {
169 
170  TLeaf * leaf = dynamic_cast<TLeaf*>(leaves[iLeaf]);
171 
172  // we can't deal with arrays
173  Int_t size = leaf->GetNdata();
174  if (size != 1) {
175  cerr<<"ERROR in LauSPlot::readInputInfo : Tree has array branches, can't deal with those."<<endl;
176  return kFALSE;
177  }
178 
179  // find the name of the leaf
180  TString name = leaf->GetName();
181 
182  // initialise an entry in the maps to hold the value
183  leaves_[name] = leaf;
184  }
185 
186  return kTRUE;
187 }
188 
189 Bool_t LauSPlot::checkLeaves() const
190 {
191  // Analyse the names of the leaves to check that we have a leaf for
192  // all bits of information we require, i.e. a likelihood value for
193  // each combination of variable and species
194 
195  // If we have 2D PDFs then we have to look for some extra leaves
196  for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
197 
198  const TString& specName = twodim_iter->first;
199  const TString& firstVarName = twodim_iter->second.first;
200  const TString& secondVarName = twodim_iter->second.second;
201 
202  TString expectedName(specName);
203  expectedName += firstVarName;
204  expectedName += secondVarName;
205  expectedName += "Like";
206 
207  Bool_t found(kFALSE);
208  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
209  const TString& leafName = leaf_iter->first;
210  const TLeaf* leaf = leaf_iter->second;
211 
212  if ((leafName == expectedName) && (leaf != 0)) {
213  found = kTRUE;
214  break;
215  }
216  }
217 
218  if (!found) {
219  cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
220  return kFALSE;
221  }
222  }
223 
224  // Search for all other "standard" leaves, i.e. <species><var>Like
225 
226  for (NumbMap::const_iterator fixd_iter = fixdSpecies_.begin(); fixd_iter != fixdSpecies_.end(); ++fixd_iter) {
227 
228  const TString& specName = fixd_iter->first;
229 
230  // if the signal is split we need to do a dedicated search
231  // for sigTM and sigSCF, so skip the signal here
232  if ( specName == "sig" && this->signalSplit() ) {
233  continue;
234  }
235 
236  for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
237  const TString& varName = (*vars_iter);
238  TString expectedName(specName);
239  expectedName += varName;
240  expectedName += "Like";
241  Bool_t found(kFALSE);
242  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
243  const TString& leafName = leaf_iter->first;
244  const TLeaf* leaf = leaf_iter->second;
245 
246  if ((leafName == expectedName) && (leaf != 0)) {
247  found = kTRUE;
248  break;
249  }
250  }
251 
252  if (!found) {
253  cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
254  return kFALSE;
255  }
256  }
257  }
258 
259  for (NumbMap::const_iterator free_iter = freeSpecies_.begin(); free_iter != freeSpecies_.end(); ++free_iter) {
260 
261  const TString& specName = free_iter->first;
262 
263  // if the signal is split we need to do a dedicated search
264  // for sigTM and sigSCF, so skip the signal here
265  if ( specName == "sig" && this->signalSplit() ) {
266  continue;
267  }
268 
269  for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
270  const TString& varName = (*vars_iter);
271  TString expectedName(specName);
272  expectedName += varName;
273  expectedName += "Like";
274  Bool_t found(kFALSE);
275  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
276  const TString& leafName = leaf_iter->first;
277  const TLeaf* leaf = leaf_iter->second;
278 
279  if ((leafName == expectedName) && (leaf != 0)) {
280  found = kTRUE;
281  break;
282  }
283  }
284 
285  if (!found) {
286  cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
287  return kFALSE;
288  }
289  }
290  }
291 
292  if ( this->signalSplit() ) {
293 
294  // now need to search for the sigTM and sigSCF leaves
295 
296  for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
297  const TString& varName = (*vars_iter);
298  TString expectedName("sigTM");
299  expectedName += varName;
300  expectedName += "Like";
301  Bool_t found(kFALSE);
302  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
303  const TString& leafName = leaf_iter->first;
304  const TLeaf* leaf = leaf_iter->second;
305 
306  if ((leafName == expectedName) && (leaf != 0)) {
307  found = kTRUE;
308  break;
309  }
310  }
311 
312  if (!found) {
313  cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
314  return kFALSE;
315  }
316 
317  expectedName = "sigSCF";
318  expectedName += varName;
319  expectedName += "Like";
320  found = kFALSE;
321  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
322  const TString& leafName = leaf_iter->first;
323  const TLeaf* leaf = leaf_iter->second;
324 
325  if ((leafName == expectedName) && (leaf != 0)) {
326  found = kTRUE;
327  break;
328  }
329  }
330 
331  if (!found) {
332  cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
333  return kFALSE;
334  }
335 
336  expectedName = "sigSCFFrac";
337  found = kFALSE;
338  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
339  const TString& leafName = leaf_iter->first;
340  const TLeaf* leaf = leaf_iter->second;
341 
342  if ((leafName == expectedName) && (leaf != 0)) {
343  found = kTRUE;
344  break;
345  }
346  }
347 
348  if (!found) {
349  cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
350  return kFALSE;
351  }
352  }
353  }
354 
355  return kTRUE;
356 }
357 
359 {
360  // check whether we've already created the tree
361  if (!cnTree_) {
362  // if not change to the file's directory and create the tree
364  cnTreeName_ += "_cNCoeffs";
365  file_->cd();
366  cnTree_ = new TTree(cnTreeName_, cnTreeName_);
367  cnTree_->SetDirectory(file_);
368  this->definedCNBranches(kFALSE);
369  }
370 }
371 
373 {
374  // check whether we've already created the tree
375  if (!sweightTree_) {
376  // if not change to the file's directory and create the tree
378  sweightTreeName_ += "_sWeights";
379  file_->cd();
381  sweightTree_->SetDirectory(file_);
382  this->definedSWeightBranches(kFALSE);
383  }
384 }
385 
387 {
388  if (this->definedCNBranches()) {
389  cerr<<"ERROR in LauSPlot::defineCNBranches : Already defined branches, not doing it again."<<endl;
390  return;
391  }
392  if (cN_.empty()) {
393  cerr<<"ERROR in LauSPlot::defineCNBranches : No entries in the cN container, can't define branches."<<endl;
394  return;
395  }
396  if (!cnTree_) {
397  cerr<<"ERROR in LauSPlot::defineCNBranches : Invalid pointer to the tree, can't define branches."<<endl;
398  return;
399  }
400 
401  // In the cN tree there is one entry per experiment, so need to know the experiment number.
402  cnTree_->Branch("iExpt", &iExpt_, "iExpt/I");
403 
404  for (std::map<TString,NumbMap>::iterator excl_iter = cN_.begin(); excl_iter != cN_.end(); ++excl_iter) {
405  const TString& exclName = excl_iter->first;
406  NumbMap& species = excl_iter->second;
407  for (NumbMap::iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
408  const TString& specName = spec_iter->first;
409  Double_t * pointer = &(spec_iter->second);
410  TString name(specName); name += "_cN";
411  if (exclName == "none") {
412  name += "_all";
413  } else {
414  name += "_no";
415  name += exclName;
416  }
417  TString thirdPart(name); thirdPart += "/D";
418  cnTree_->Branch(name, pointer, thirdPart);
419  }
420  }
421  this->definedCNBranches(kTRUE);
422 }
423 
425 {
426  if (this->definedSWeightBranches()) {
427  cerr<<"ERROR in LauSPlot::defineSWeightBranches : Already defined branches, not doing it again."<<endl;
428  return;
429  }
430  if (sWeights_.empty()) {
431  cerr<<"ERROR in LauSPlot::defineSWeightBranches : No entries in the sWeights container, can't define branches."<<endl;
432  return;
433  }
434  if (!sweightTree_) {
435  cerr<<"ERROR in LauSPlot::defineSWeightBranches : Invalid pointer to the tree, can't define branches."<<endl;
436  return;
437  }
438 
439  for (std::map<TString,NumbMap>::iterator excl_iter = sWeightsCurrent_.begin(); excl_iter != sWeightsCurrent_.end(); ++excl_iter) {
440  const TString& exclName = excl_iter->first;
441  NumbMap& species = excl_iter->second;
442  for (NumbMap::iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
443  const TString& specName = spec_iter->first;
444  Double_t * pointer = &(spec_iter->second);
445  TString name(specName); name += "_sWeight";
446  if (exclName == "none") {
447  name += "_all";
448  } else {
449  name += "_no";
450  name += exclName;
451  }
452  TString thirdPart(name); thirdPart += "/D";
453  sweightTree_->Branch(name, pointer, thirdPart);
454  }
455  }
456  this->definedSWeightBranches(kTRUE);
457 }
458 
459 void LauSPlot::setExperiment(Int_t iExpt)
460 {
461  if (!inputTree_) {
462  cerr<<"ERROR in LauSPlot::setExperiment : Invalid pointer to data tree."<<endl;
463  return;
464  }
465 
466  // create the event list if we haven't already done so
467  if (!eventList_) {
468  eventList_ = new TEventList("splotelist","splotelist");
469  eventList_->SetDirectory(file_);
470  }
471 
472  // fill the event list with this experiment's events
473  TString listName(eventList_->GetName());
474  listName.Prepend(">>");
475  TString selection("iExpt==");
476  selection += iExpt;
477  cout<<"LauSPlot::setExperiment : Setting tree to experiment number "<<iExpt<<"."<<endl;
478  inputTree_->Draw(listName,selection);
479 
480  // find how many events there are
481  nEvents_ = eventList_->GetN();
482  cout<<" Found "<<nEvents_<<" events."<<endl;
483 
484  // make sure we have enough space in the per-event value vectors
485  pdfTot_.clear(); pdfTot_.resize(nEvents_);
486  discPdf_.clear(); discPdf_.resize(nEvents_);
487  scfFrac_.clear(); scfFrac_.resize(nEvents_);
488  sWeights_.clear(); sWeights_.resize(nEvents_);
489 
490  // read the info for this experiment
491  this->readExpt();
492 }
493 
495 {
496  for (Int_t iEvt(0); iEvt < nEvents_; ++iEvt) {
497  // Find which entry from the full tree contains the requested event
498  Long64_t iEntry = eventList_ ? eventList_->GetEntry(iEvt) : iEvt;
499  if (iEntry<0) { // this shouldn't happen, but just in case...
500  cerr<<"ERROR in LauSPlot::readExpt : Problem retrieving event."<<endl;
501  gSystem->Exit(EXIT_FAILURE);
502  }
503 
504  // Then retrieve that entry from the tree
505  inputTree_->GetEntry(iEntry);
506 
507  // If needed retrieve the SCF fraction values
508  if ( this->signalSplit() ) {
509  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
510  const TString& leafName = leaf_iter->first;
511  const TLeaf* leaf = leaf_iter->second;
512  if ((leafName == "sigSCFFrac") && (leaf != 0)) {
513  scfFrac_[iEvt] = leaf->GetValue();
514  break;
515  }
516  }
517  }
518 
519  // Copy the leaf values into discPdf_
520  for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
521  const TString& specName = twodim_iter->first;
522  const TString& firstVarName = twodim_iter->second.first;
523  const TString& secondVarName = twodim_iter->second.second;
524 
525  TString varName = firstVarName + secondVarName;
526 
527  TString expectedName(specName);
528  expectedName += varName;
529  expectedName += "Like";
530 
531  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
532  const TString& leafName = leaf_iter->first;
533  const TLeaf* leaf = leaf_iter->second;
534  if ((leafName == expectedName) && (leaf != 0)) {
535  discPdf_[iEvt][specName][varName] = leaf->GetValue();
536  break;
537  }
538  }
539  }
540 
541  Bool_t needSignalSearch(kFALSE);
542  for (NumbMap::const_iterator fixd_iter = fixdSpecies_.begin(); fixd_iter != fixdSpecies_.end(); ++fixd_iter) {
543  const TString& specName = fixd_iter->first;
544 
545  // if the signal is split we need to do a dedicated search
546  // for sigTM and sigSCF, so skip the signal here
547  if ( specName == "sig" && this->signalSplit() ) {
548  needSignalSearch = kTRUE;
549  continue;
550  }
551 
552  for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
553  const TString& varName = (*vars_iter);
554  TString expectedName(specName);
555  expectedName += varName;
556  expectedName += "Like";
557  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
558  const TString& leafName = leaf_iter->first;
559  const TLeaf* leaf = leaf_iter->second;
560  if ((leafName == expectedName) && (leaf != 0)) {
561  discPdf_[iEvt][specName][varName] = leaf->GetValue();
562  break;
563  }
564  }
565  }
566  }
567 
568  for (NumbMap::const_iterator free_iter = freeSpecies_.begin(); free_iter != freeSpecies_.end(); ++free_iter) {
569  const TString& specName = free_iter->first;
570 
571  // if the signal is split we need to do a dedicated search
572  // for sigTM and sigSCF, so skip the signal here
573  if ( specName == "sig" && this->signalSplit() ) {
574  needSignalSearch = kTRUE;
575  continue;
576  }
577 
578  for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
579  const TString& varName = (*vars_iter);
580  TString expectedName(specName);
581  expectedName += varName;
582  expectedName += "Like";
583  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
584  const TString& leafName = leaf_iter->first;
585  const TLeaf* leaf = leaf_iter->second;
586  if ((leafName == expectedName) && (leaf != 0)) {
587  discPdf_[iEvt][specName][varName] = leaf->GetValue();
588  break;
589  }
590  }
591  }
592  }
593 
594  if ( needSignalSearch ) {
595  for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
596  const TString& varName = (*vars_iter);
597  TString specName("sigTM");
598  TString expectedName(specName);
599  expectedName += varName;
600  expectedName += "Like";
601  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
602  const TString& leafName = leaf_iter->first;
603  const TLeaf* leaf = leaf_iter->second;
604  if ((leafName == expectedName) && (leaf != 0)) {
605  discPdf_[iEvt][specName][varName] = leaf->GetValue();
606  break;
607  }
608  }
609  specName = "sigSCF";
610  expectedName = specName;
611  expectedName += varName;
612  expectedName += "Like";
613  for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
614  const TString& leafName = leaf_iter->first;
615  const TLeaf* leaf = leaf_iter->second;
616  if ((leafName == expectedName) && (leaf != 0)) {
617  discPdf_[iEvt][specName][varName] = leaf->GetValue();
618  break;
619  }
620  }
621  }
622  }
623  }
624 }
625 
626 void LauSPlot::runCalculations(const TString& option)
627 {
628  // Calculates the sWeights and cN coeffs
629  // The option controls the print level:
630  // "Q" - no print out (default)
631  // "V" - prints the estimated # of events in species
632  // "VV" - as "V" + the MINUIT printing + sums of weights for control
633 
634  if (!this->readInput()) {
635  cerr<<"ERROR in LauSPlot::runCalculations : The input ntuple has not been successfully read, can't calculate anything."<<endl;
636  return;
637  }
638 
639  if (freeSpecies_.empty()) {
640  cerr<<"ERROR in LauSPlot::runCalculations : Numbers of events in species have not been set."<<endl;
641  return;
642  }
643 
644  TString opt(option);
645  opt.ToUpper();
646  opt.ReplaceAll("VV", "W");
647 
648  // Make sure that global fitter is MINUIT
649  this->checkFitter();
650 
651  // Loop over experiments
652  for (iExpt_ = firstExpt_; iExpt_ < (firstExpt_+nExpt_); ++iExpt_) {
653 
654  this->setExperiment(iExpt_);
655 
656  if (nEvents_ < 1) {
657  cerr<<"ERROR in LauSPlot::runCalculations : Zero events in experiment "<<iExpt_<<", skipping..."<<endl;
658  continue;
659  }
660 
661  // Now loop over the PDFs to exclude, including the case where none are to be excluded.
662  NameSet excludePdf;
663  if (variableNames_.size()<2) {
664  excludePdf.insert("none");
665  } else {
666  excludePdf = variableNames_;
667  excludePdf.insert("none");
668  }
669 
670  for (NameSet::const_iterator excl_iter = excludePdf.begin(); excl_iter != excludePdf.end(); ++excl_iter) {
671 
672  const TString& exclName = (*excl_iter);
673 
674  cout<<"LauSPlot::runCalculations : Calculating sWeights, excluding PDF: "<<exclName<<"."<<endl;
675 
676  // Calculate the per-event total PDF values for each species.
677  this->calcTotPDFValues(exclName);
678 
679  // Reset the fitter
680  this->initialiseFitter(opt);
681 
682  // Set the parameters to their initial values
683  this->setFitParameters();
684 
685  // Run the fit
686  this->runFit();
687 
688  // Get the fitted parameter values back
689  this->retrieveFittedParameters(opt);
690 
691  // Get the covariance matrix
692  Bool_t covMatOK = this->calcCovMatrix();
693  Double_t * covmat(0);
694  if (!covMatOK) {
695  TVirtualFitter * fitter = TVirtualFitter::GetFitter();
696  covmat = fitter->GetCovarianceMatrix();
697  }
698  if (opt.Contains("W")) {
699  this->printCovMatrixElements(covmat);
700  }
701 
702  // calculate the cN and sWeights from the covariance matrix
703  if (nFixdSpecies_ > 0) {
704  this->calcCNCoeffs(exclName, covmat);
705  }
706  this->calcSWeights(exclName, covmat);
707 
708  // print verbose info if required
709  if (opt.Contains("W")) {
710  this->printSumOfWeights(exclName);
711  }
712  }
713 
714  // Finally fill all the branches for this experiment
715  if (nFixdSpecies_ > 0) {
716  this->fillCNBranches();
717  }
718  this->fillSWeightBranches();
719  }
720 }
721 
723 {
724  TString minuitName("TFitter");
725  TVirtualFitter * fitter = TVirtualFitter::GetFitter();
726  if (fitter) {
727  TString fitterName(fitter->IsA()->GetName());
728  if (fitterName != minuitName) {
729  delete fitter; fitter = 0;
730  }
731  }
732  fitter = TVirtualFitter::Fitter(0);
733 }
734 
735 void LauSPlot::initialiseFitter(const TString& opt)
736 {
737  TVirtualFitter * fitter = TVirtualFitter::GetFitter();
738  fitter->Clear();
739  fitter->SetFCN(Yields);
740  fitter->SetObjectFit(this);
741 
742  // Set the print level
743  Double_t arglist[10];
744  if (opt.Contains("Q")) {
745  arglist[0]=-1;
746  }
747  if (opt.Contains("V")) {
748  arglist[0]=0;
749  }
750  if (opt.Contains("W")) {
751  arglist[0]=1;
752  }
753  fitter->ExecuteCommand("SET PRINT", arglist, 1);
754 
755  // Need to set the "SET ERR" command to +0.5 for +/-1 sigma errors
756  // for maximum likelihood fit.
757  arglist[0] = 0.5;
758  fitter->ExecuteCommand("SET ERR", arglist, 1);
759 }
760 
762 {
763  TVirtualFitter * fitter = TVirtualFitter::GetFitter();
764 
765  // must add the parameters in the same order as they are stored in pdfTot_
766  Int_t ispecies(0);
767  const NumbMap& species = pdfTot_.front();
768 
769  for (NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
770  const TString& name(spec_iter->first);
771 
772  // starting parameters should be the original values,
773  // not those that came out of the last fit
774  NumbMap::const_iterator free_iter = origFreeSpecies_.find(name);
775  NumbMap::const_iterator fixd_iter = origFixdSpecies_.find(name);
776  Bool_t fixed = fixd_iter != origFixdSpecies_.end();
777 
778  Double_t value(0.0);
779  if (fixed) {
780  value = fixd_iter->second;
781  } else {
782  value = free_iter->second;
783  }
784  fitter->SetParameter(ispecies, name, value, 1, 0, 0);
785  if (fixed) {
786  fitter->FixParameter(ispecies);
787  }
788  ++ispecies;
789  }
790 }
791 
793 {
794  TVirtualFitter * fitter = TVirtualFitter::GetFitter();
795 
796  Double_t arglist[10];
797  arglist[0] = 1000*nFreeSpecies_; // maximum iterations
798  arglist[1] = 0.05; // tolerance : min EDM = 0.001*tolerance
799 
800  // Execute MIGRAD
801  Int_t fitStatus = fitter->ExecuteCommand("MIGRAD", arglist, 2);
802 
803  if (fitStatus != 0) {
804  cerr<<"ERROR in LauSPlot::runFit : Error during MIGRAD minimisation."<<endl;
805  } else {
806  // Execute HESSE
807  fitStatus = fitter->ExecuteCommand("HESSE", arglist, 1);
808  if (fitStatus != 0) {
809  cerr<<"ERROR in LauSPlot::runFit : Error during HESSE error calculation."<<endl;
810  }
811  }
812 }
813 
814 void LauSPlot::printCovMatrixElements(const Double_t * covmat) const
815 {
816  Int_t a(0);
817  cout<<endl;
818  for (NumbMap::const_iterator iter_a = freeSpecies_.begin(); iter_a != freeSpecies_.end(); ++iter_a) {
819  Int_t b(0);
820  for (NumbMap::const_iterator iter_b = freeSpecies_.begin(); iter_b != freeSpecies_.end(); ++iter_b) {
821  if (covmat) {
822  cout<<"CovMat Elem: ["<<iter_a->first<<","<<iter_b->first<<"] = "<<covmat[a*nSpecies_+b]<<endl;
823  } else {
824  cout<<"CovMat Elem: ["<<iter_a->first<<","<<iter_b->first<<"] = "<<covMat_(a,b)<<endl;
825  }
826  ++b;
827  }
828  ++a;
829  }
830  cout<<endl;
831 
832  if (!covmat) {
833  covMat_.Print();
834  }
835 }
836 
837 void LauSPlot::retrieveFittedParameters(const TString& opt)
838 {
839  TVirtualFitter * fitter = TVirtualFitter::GetFitter();
840 
841  // remember fit parameters are in same order as in pdfTot_
842  Int_t ispecies(0);
843  const NumbMap& species = pdfTot_.front();
844 
845  for (NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
846  const TString& name(spec_iter->first);
847  NumbMap::iterator free_iter = freeSpecies_.find(name);
848  if (free_iter != freeSpecies_.end()) {
849  free_iter->second = fitter->GetParameter(ispecies);
850  if (!opt.Contains("Q")) {
851  cout<<"Estimated # of events in species "<<name<<" = "<<free_iter->second<<endl;
852  }
853  }
854  ++ispecies;
855  }
856  if (!opt.Contains("Q")) {cout<<endl;}
857 }
858 
859 void LauSPlot::printSumOfWeights(const TString& exclName) const
860 {
861  for (NumbMap::const_iterator spec_iter = freeSpecies_.begin(); spec_iter != freeSpecies_.end(); ++spec_iter) {
862  const TString& specName = spec_iter->first;
863  Double_t sumweight(0.0);
864  for (Int_t iEvt(0); iEvt < nEvents_; ++iEvt) {
865  const NumbMap& specWeightMap = sWeights_[iEvt].find(exclName)->second;
866  Double_t weight = specWeightMap.find(specName)->second;
867  sumweight += weight;
868  }
869  cout<<"Sum of sWeights for species \""<<specName<<"\" = "<<sumweight<<endl;
870  }
871  cout<<endl;
872 }
873 
875 {
876  // Calculate our inverse covariance matrix from the various PDFs
877 
878  TMatrixD invMatrix(nFreeSpecies_,nFreeSpecies_);
879 
880  // First calculate the denominator, which is common to all elements
881  vector<Double_t> denominator(nEvents_);
882  for (Int_t iEvt(0); iEvt<nEvents_; ++iEvt) {
883  denominator[iEvt] = 0.0;
884  for (NumbMap::const_iterator spec_iter = freeSpecies_.begin(); spec_iter != freeSpecies_.end(); ++spec_iter) {
885  const TString& specName = spec_iter->first;
886  Double_t specNumEvents = spec_iter->second;
887  denominator[iEvt] += specNumEvents * pdfTot_[iEvt][specName];
888  }
889  for (NumbMap::const_iterator spec_iter = fixdSpecies_.begin(); spec_iter != fixdSpecies_.end(); ++spec_iter) {
890  const TString& specName = spec_iter->first;
891  Double_t specNumEvents = spec_iter->second;
892  denominator[iEvt] += specNumEvents * pdfTot_[iEvt][specName];
893  }
894  // Square to get the final denominator
895  denominator[iEvt] *= denominator[iEvt];
896  }
897 
898  // Then calculate each element
899  Int_t i(0);
900  for (NumbMap::const_iterator spec_iter_i = freeSpecies_.begin(); spec_iter_i != freeSpecies_.end(); ++spec_iter_i) {
901  const TString& specName_i = spec_iter_i->first;
902  Int_t j(0);
903  for (NumbMap::const_iterator spec_iter_j = freeSpecies_.begin(); spec_iter_j != freeSpecies_.end(); ++spec_iter_j) {
904  const TString& specName_j = spec_iter_j->first;
905  invMatrix(i,j) = 0.0;
906  for (Int_t iEvt(0); iEvt<nEvents_; ++iEvt) {
907  Double_t numerator = pdfTot_[iEvt][specName_i] * pdfTot_[iEvt][specName_j];
908  invMatrix(i,j) += numerator/denominator[iEvt];
909  }
910  ++j;
911  }
912  ++i;
913  }
914 
915  // Check for a singular matrix
916  if (invMatrix.Determinant() < 1e-15) {
917  cerr<<"ERROR in LauSPlot::calcCovMatrix : Calculated inverse covariance matrix is singular, can't invert it."<<endl;
918  return kFALSE;
919  }
920 
921  // Invert and store in the covariance matrix
923  covMat_ = TMatrixD(TMatrixD::kInverted, invMatrix);
924 
925  return kTRUE;
926 }
927 
928 void LauSPlot::calcTotPDFValues(const TString& exclName)
929 {
930  for (Int_t iEvt(0); iEvt<nEvents_; iEvt++) {
931 
932  Bool_t needSignalSearch(kFALSE);
933 
934  for (NumbMap::const_iterator spec_iter = fixdSpecies_.begin(); spec_iter != fixdSpecies_.end(); ++spec_iter) {
935  const TString& specName = spec_iter->first;
936 
937  // if the signal is split we need to do a dedicated search
938  // for sigTM and sigSCF, so skip the signal here
939  if ( specName == "sig" && this->signalSplit() ) {
940  needSignalSearch = kTRUE;
941  continue;
942  }
943 
944  pdfTot_[iEvt][specName] = 1.0;
945 
946  // loop through the 2D histo list
947  NameSet skipList;
948  for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
949  // if the entry doesn't refer to this
950  // species then skip on
951  if ( specName != twodim_iter->first ) {
952  continue;
953  }
954 
955  // retrieve the two variable names
956  const TString& firstVarName = twodim_iter->second.first;
957  const TString& secondVarName = twodim_iter->second.second;
958  if ( firstVarName != exclName && secondVarName != exclName ) {
959  // if neither name is the one being excluded then...
960  // add them both to the skip list
961  skipList.insert( firstVarName );
962  skipList.insert( secondVarName );
963  // and multiply the total by the combined PDF value
964  TString varName = firstVarName + secondVarName;
965  pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
966  }
967  }
968 
969  // loop through all the variables
970  for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
971  const TString& varName = (*var_iter);
972  // if the variable isn't the one being excluded
973  if (exclName != varName) {
974  // and it's not involved in a 2D PDF
975  if ( skipList.find(varName) == skipList.end() ) {
976  // multiply the total by its PDF value
977  pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
978  }
979  }
980  }
981  }
982 
983  for (NumbMap::const_iterator spec_iter = freeSpecies_.begin(); spec_iter != freeSpecies_.end(); ++spec_iter) {
984  const TString& specName = spec_iter->first;
985 
986  // if the signal is split we need to do a dedicated search
987  // for sigTM and sigSCF, so skip the signal here
988  if ( specName == "sig" && this->signalSplit() ) {
989  needSignalSearch = kTRUE;
990  continue;
991  }
992 
993  pdfTot_[iEvt][specName] = 1.0;
994 
995  // loop through the 2D histo list
996  NameSet skipList;
997  for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
998  // if the entry doesn't refer to this
999  // species then skip on
1000  if ( specName != twodim_iter->first ) {
1001  continue;
1002  }
1003 
1004  // retrieve the two variable names
1005  const TString& firstVarName = twodim_iter->second.first;
1006  const TString& secondVarName = twodim_iter->second.second;
1007  if ( firstVarName != exclName && secondVarName != exclName ) {
1008  // if neither name is the one being excluded then...
1009  // add them both to the skip list
1010  skipList.insert( firstVarName );
1011  skipList.insert( secondVarName );
1012  // and multiply the total by the combined PDF value
1013  TString varName = firstVarName + secondVarName;
1014  pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
1015  }
1016  }
1017 
1018  // loop through all the variables
1019  for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
1020  const TString& varName = (*var_iter);
1021  // if the variable isn't the one being excluded
1022  if (exclName != varName) {
1023  // and it's not involved in a 2D PDF
1024  if ( skipList.find(varName) == skipList.end() ) {
1025  // multiply the total by its PDF value
1026  pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
1027  }
1028  }
1029  }
1030  }
1031 
1032  if ( needSignalSearch ) {
1033 
1034  // loop through the 2D histo list
1035  TString specName("sigTM");
1036  Double_t tmPDFVal(1.0);
1037  NameSet skipList;
1038  for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
1039  // if the entry doesn't refer to this
1040  // species then skip on
1041  if ( specName != twodim_iter->first ) {
1042  continue;
1043  }
1044 
1045  // retrieve the two variable names
1046  const TString& firstVarName = twodim_iter->second.first;
1047  const TString& secondVarName = twodim_iter->second.second;
1048  if ( firstVarName != exclName && secondVarName != exclName ) {
1049  // if neither name is the one being excluded then...
1050  // add them both to the skip list
1051  skipList.insert( firstVarName );
1052  skipList.insert( secondVarName );
1053  // and multiply the total by the combined PDF value
1054  TString varName = firstVarName + secondVarName;
1055  tmPDFVal *= discPdf_[iEvt][specName][varName];
1056  }
1057  }
1058 
1059  // loop through all the variables
1060  for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
1061  const TString& varName = (*var_iter);
1062  // if the variable isn't the one being excluded
1063  if (exclName != varName) {
1064  // and it's not involved in a 2D PDF
1065  if ( skipList.find(varName) == skipList.end() ) {
1066  // multiply the total by its PDF value
1067  tmPDFVal *= discPdf_[iEvt][specName][varName];
1068  }
1069  }
1070  }
1071 
1072  tmPDFVal *= (1.0 - scfFrac_[iEvt]);
1073 
1074 
1075  // loop through the 2D histo list
1076  specName = "sigSCF";
1077  Double_t scfPDFVal(1.0);
1078  skipList.clear();
1079  for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
1080  // if the entry doesn't refer to this
1081  // species then skip on
1082  if ( specName != twodim_iter->first ) {
1083  continue;
1084  }
1085 
1086  // retrieve the two variable names
1087  const TString& firstVarName = twodim_iter->second.first;
1088  const TString& secondVarName = twodim_iter->second.second;
1089  if ( firstVarName != exclName && secondVarName != exclName ) {
1090  // if neither name is the one being excluded then...
1091  // add them both to the skip list
1092  skipList.insert( firstVarName );
1093  skipList.insert( secondVarName );
1094  // and multiply the total by the combined PDF value
1095  TString varName = firstVarName + secondVarName;
1096  scfPDFVal *= discPdf_[iEvt][specName][varName];
1097  }
1098  }
1099 
1100  // loop through all the variables
1101  for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
1102  const TString& varName = (*var_iter);
1103  // if the variable isn't the one being excluded
1104  if (exclName != varName) {
1105  // and it's not involved in a 2D PDF
1106  if ( skipList.find(varName) == skipList.end() ) {
1107  // multiply the total by its PDF value
1108  scfPDFVal *= discPdf_[iEvt][specName][varName];
1109  }
1110  }
1111  }
1112 
1113  if ( exclName == "DP" || !this->scfDPSmear() ) {
1114  scfPDFVal *= scfFrac_[iEvt];
1115  }
1116 
1117  pdfTot_[iEvt]["sig"] = tmPDFVal + scfPDFVal;
1118  }
1119  }
1120 }
1121 
1122 void LauSPlot::calcCNCoeffs(const TString& exclName, const Double_t *covmat)
1123 {
1124  // Computes the cN for the extended sPlots from the covariance matrix
1125 
1126  if (nFixdSpecies_ <= 0) {
1127  return;
1128  }
1129 
1130  Int_t species_n(0);
1131  for (NumbMap::const_iterator iter_n = freeSpecies_.begin(); iter_n != freeSpecies_.end(); ++iter_n) {
1132  Int_t species_j(0);
1133  const TString& specName = iter_n->first;
1134  Double_t value = iter_n->second;
1135  cN_[exclName][specName] = value;
1136  for (NumbMap::const_iterator iter_j = freeSpecies_.begin(); iter_j != freeSpecies_.end(); ++iter_j) {
1137  if (covmat) {
1138  cN_[exclName][specName] -= covmat[species_n*nSpecies_+species_j];
1139  } else {
1140  cN_[exclName][specName] -= covMat_(species_n,species_j);
1141  }
1142  ++species_j;
1143  }
1144  ++species_n;
1145  }
1146 }
1147 
1148 void LauSPlot::calcSWeights(const TString& exclName, Double_t *covmat)
1149 {
1150  // Computes the sWeights from the covariance matrix
1151  // NB for the extended sPlot the sum in the denominator is still over all species,
1152  // while that in the numerator is only over the free species.
1153  // Similarly the sWeights can only be calculated for the free species.
1154 
1155  Double_t numerator(0.0), denominator(0.0);
1156  for (Int_t iEvent = 0; iEvent < nEvents_; ++iEvent) {
1157  denominator = 0.0;
1158  for (NumbMap::const_iterator free_iter = freeSpecies_.begin(); free_iter != freeSpecies_.end(); ++free_iter) {
1159  denominator += free_iter->second * pdfTot_[iEvent][free_iter->first];
1160  }
1161  for (NumbMap::const_iterator fixd_iter = fixdSpecies_.begin(); fixd_iter != fixdSpecies_.end(); ++fixd_iter) {
1162  denominator += fixd_iter->second * pdfTot_[iEvent][fixd_iter->first];
1163  }
1164  Int_t species_n(0);
1165  for (NumbMap::const_iterator iter_n = freeSpecies_.begin(); iter_n != freeSpecies_.end(); ++iter_n) {
1166  numerator = 0.0;
1167  Int_t species_j(0);
1168  for (NumbMap::const_iterator iter_j = freeSpecies_.begin(); iter_j != freeSpecies_.end(); ++iter_j) {
1169  if (covmat) {
1170  numerator += covmat[species_n*nSpecies_+species_j] * pdfTot_[iEvent][iter_j->first];
1171  } else {
1172  numerator += covMat_(species_n,species_j) * pdfTot_[iEvent][iter_j->first];
1173  }
1174  ++species_j;
1175  }
1176  sWeights_[iEvent][exclName][iter_n->first] = numerator/denominator;
1177  ++species_n;
1178  }
1179  }
1180 }
1181 
1183 {
1184  if (!cnTree_) {
1185  cerr<<"ERROR in LauSPlot::fillCNBranches : Tree not created, cannot fill branches."<<endl;
1186  return;
1187  } else if (!this->definedCNBranches()) {
1188  this->defineCNBranches();
1189  }
1190  cnTree_->Fill();
1191 }
1192 
1194 {
1195  if (!sweightTree_) {
1196  cerr<<"ERROR in LauSPlot::fillSWeightBranches : Tree not created, cannot fill branches."<<endl;
1197  return;
1198  } else if (sWeights_.empty()) {
1199  cerr<<"ERROR in LauSPlot::fillSWeightBranches : No sWeights calculated, can't fill branches."<<endl;
1200  return;
1201  } else if (!this->definedSWeightBranches()) {
1202  this->copyEventWeights(0);
1203  this->defineSWeightBranches();
1204  }
1205 
1206  for (Int_t iEvent = 0; iEvent < nEvents_; ++iEvent) {
1207  this->copyEventWeights(iEvent);
1208  sweightTree_->Fill();
1209  }
1210 }
1211 
1212 void LauSPlot::copyEventWeights(Int_t iEvent)
1213 {
1214  for (std::map<TString,NumbMap>::const_iterator excl_iter = sWeights_[iEvent].begin(); excl_iter != sWeights_[iEvent].end(); ++excl_iter) {
1215  const TString& exclName = excl_iter->first;
1216  const NumbMap& species = excl_iter->second;
1217  for (NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
1218  const TString& specName = spec_iter->first;
1219  sWeightsCurrent_[exclName][specName] = spec_iter->second;
1220  }
1221  }
1222 }
1223 
1225 {
1226  // write out the results
1227 
1228  // remove the transient objects that we don't want (re-)written to the file
1229  if (eventList_) {
1230  delete eventList_; eventList_ = 0;
1231  }
1232  if (inputTree_) {
1233  delete inputTree_; inputTree_ = 0;
1234  }
1235 
1236  // first add the input tree as a friend of the output tree
1237  this->addFriendTree();
1238 
1239  // then write everything to the file and clean up
1240  file_->cd();
1241  file_->Write();
1242  file_->Close();
1243  delete file_; file_ = 0;
1244 }
1245 
1247 {
1248  if (!sweightTree_) {
1249  cerr<<"ERROR in LauSPlot::addFriendTree : Tree not created, cannot add friend."<<endl;
1250  return;
1251  }
1253 }
1254 
1255 void Yields(Int_t &, Double_t *, Double_t &f, Double_t *x, Int_t /*iflag*/)
1256 {
1257  // FCN-function for Minuit
1258 
1259  f = 0.0;
1260 
1261  TVirtualFitter *fitter = TVirtualFitter::GetFitter();
1262  LauSPlot* theModel = dynamic_cast<LauSPlot*>(fitter->GetObjectFit());
1263  const std::vector<LauSPlot::NumbMap>& pdfTot = theModel->totalPdf();
1264 
1265  Double_t ntot(0.0);
1266  for (std::vector<LauSPlot::NumbMap>::const_iterator evt_iter = pdfTot.begin(); evt_iter != pdfTot.end(); ++evt_iter) { // loop over events
1267  const LauSPlot::NumbMap& species = (*evt_iter);
1268  Int_t ispecies(0);
1269  Double_t lik(0.0);
1270  ntot = 0.0;
1271  for (LauSPlot::NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) { // loop over species
1272  lik += x[ispecies] * spec_iter->second;
1273  ntot += x[ispecies];
1274  ++ispecies;
1275  }
1276  if (lik < 0.0) {
1277  // make f the worst possible value to force MINUIT
1278  // out of this region of parameter space
1279  f = -DBL_MAX;
1280  break;
1281  } else {
1282  f += TMath::Log(lik);
1283  }
1284  }
1285 
1286  // extended likelihood
1287  f = (ntot-f);
1288 }
1289 
Class for defining the SPlot technique.
Definition: LauSPlot.hh:55
std::map< TString, NumbMap > cN_
The extended sPlot coefficients (for each species and for each combination of excluded vars) ...
Definition: LauSPlot.hh:349
void runFit()
Perform the minimisation wrt the yields alone.
Definition: LauSPlot.cc:792
Bool_t fixed() const
Check whether the parameter is fixed or floated.
std::vector< NumbMap > pdfTot_
The per-event values of the total PDF for each species.
Definition: LauSPlot.hh:335
void initialiseFitter(const TString &opt)
Initialise Minuit, set the verbosity.
Definition: LauSPlot.cc:735
void setFitParameters() const
Add the species yields as fit parameters and fix them as appropriate.
Definition: LauSPlot.cc:761
void runCalculations(const TString &option="q")
Method to calculate the sWeights and cN coeffs.
Definition: LauSPlot.cc:626
ClassImp(LauAbsCoeffSet)
Int_t nFixdSpecies_
Number of species fixed in the fit.
Definition: LauSPlot.hh:330
void calcTotPDFValues(const TString &exclName)
Calculate the total likelihood for each species by multiply together all the PDFs for that species...
Definition: LauSPlot.cc:928
Int_t nEvents_
Number of events in current experiment.
Definition: LauSPlot.hh:324
void fillCNBranches()
Fill the cN branches.
Definition: LauSPlot.cc:1182
NumbMap freeSpecies_
The names and estimated yields of the free species.
Definition: LauSPlot.hh:295
const TString & name() const
The parameter name.
void readInputInfo()
Read the leaf structure from the tree and check the status of the read (calls LauSPlot::readInputLeav...
Definition: LauSPlot.cc:130
std::vector< std::map< TString, NumbMap > > discPdf_
The per-event values of the PDFs for each species for each disc variable.
Definition: LauSPlot.hh:337
TEventList * eventList_
Pointer to an event list, that is used to loop through the experiments.
Definition: LauSPlot.hh:287
void printSumOfWeights(const TString &exclName) const
Print the sum of sWeights for all species.
Definition: LauSPlot.cc:859
std::multimap< TString, std::pair< TString, TString > > TwoDMap
Type to associate the name of the species that have 2D PDFs with the names of the two variables invol...
Definition: LauSPlot.hh:68
void Yields(Int_t &, Double_t *, Double_t &f, Double_t *x, Int_t iflag)
Definition: LauSPlot.cc:1255
void createCNTree()
Create (if not already done) the tree for storing the cN coeffs.
Definition: LauSPlot.cc:358
TString inputTreeName_
The name of the input tree (containing the per-event llhds)
Definition: LauSPlot.hh:272
TFile * file_
Pointer to the data file object.
Definition: LauSPlot.hh:278
void fillSWeightBranches()
Fill the sWeights branches.
Definition: LauSPlot.cc:1193
void readExpt()
Reads the values of each PDF likelihood for every event in the experiment.
Definition: LauSPlot.cc:494
TwoDMap twodimPDFs_
The names of the species that have 2D PDFs and the names of the variables involved.
Definition: LauSPlot.hh:303
std::vector< Double_t > scfFrac_
The per-event values of the SCF fraction.
Definition: LauSPlot.hh:339
TString cnTreeName_
The name of the cn tree (containing the cN coefficients)
Definition: LauSPlot.hh:274
Int_t iExpt_
The current experiment.
Definition: LauSPlot.hh:322
Bool_t signalSplit() const
Check whether the signal is split into Truth Matched and Self Cross Feed.
Definition: LauSPlot.hh:131
void writeOutResults()
Save the sWeight results as a friend tree to the input tree (in the same file)
Definition: LauSPlot.cc:1224
Int_t nExpt_
Number of experiments.
Definition: LauSPlot.hh:320
void checkFitter() const
Make sure that we&#39;re using Minuit.
Definition: LauSPlot.cc:722
void calcSWeights(const TString &exclName, Double_t *covmat=0)
Computes the sWeights from the PDFs and covariance matrix.
Definition: LauSPlot.cc:1148
void defineCNBranches()
Create the branches for each cN coefficient.
Definition: LauSPlot.cc:386
NumbMap origFreeSpecies_
The names and estimated yields of the free species - need to keep the original values.
Definition: LauSPlot.hh:299
LeafMap leaves_
Collection to hold pointers to the leaves of the input tree.
Definition: LauSPlot.hh:290
TTree * cnTree_
Pointer to the output tree containing the cN coefficients.
Definition: LauSPlot.hh:282
Int_t firstExpt_
First experiment.
Definition: LauSPlot.hh:318
NameSet variableNames_
The names of the discriminating variables.
Definition: LauSPlot.hh:293
void copyEventWeights(Int_t iEvent)
Copy the sWeight of a given event into LauSPlot::sWeightsCurrent_, from which they can be stored in t...
Definition: LauSPlot.cc:1212
void openInputFileAndTree()
Method to open the file in &quot;update&quot; mode and grab the input tree for reading.
Definition: LauSPlot.cc:107
void setExperiment(Int_t iExpt)
Set the event list to contain only events from the given experiment.
Definition: LauSPlot.cc:459
std::map< TString, Double_t > NumbMap
Type to associate a category name with a double precision number, e.g. a yield or PDF value for a giv...
Definition: LauSPlot.hh:62
virtual ~LauSPlot()
Destructor.
Definition: LauSPlot.cc:95
File containing declaration of LauSPlot class.
TString sweightTreeName_
The name of the sweight tree (containing the sWeights)
Definition: LauSPlot.hh:276
void printCovMatrixElements(const Double_t *covmat=0) const
Print the supplied covariance matrix or, if pointer is null, the one previously calculated.
Definition: LauSPlot.cc:814
void createSWeightTree()
Create (if not already done) the tree for storing the sWeights.
Definition: LauSPlot.cc:372
Bool_t readInput() const
Check whether the input tree has been successfully read.
Definition: LauSPlot.hh:119
void calcCNCoeffs(const TString &exclName, const Double_t *covmat=0)
Computes the cN for the extended sPlots from the covariance matrix.
Definition: LauSPlot.cc:1122
const std::vector< LauSPlot::NumbMap > & totalPdf() const
Access the per-event total PDF values for each species.
Definition: LauSPlot.hh:112
TTree * sweightTree_
Pointer to the output tree containing the sWeights.
Definition: LauSPlot.hh:284
std::map< TString, NumbMap > sWeightsCurrent_
The current-event values of the computed sWeights.
Definition: LauSPlot.hh:347
Bool_t readInputLeaves()
Read the leaf structure from the tree and setup the leaf map.
Definition: LauSPlot.cc:146
Int_t nSpecies_
Total number of species (free + fixed)
Definition: LauSPlot.hh:332
NumbMap fixdSpecies_
The names and estimated yields of the fixed species.
Definition: LauSPlot.hh:297
Bool_t calcCovMatrix()
Definition: LauSPlot.cc:874
void retrieveFittedParameters(const TString &opt)
Update the yields with the newly fitted values and print them (unless print option is &quot;Q&quot;)...
Definition: LauSPlot.cc:837
TString fileName_
The name of the data file.
Definition: LauSPlot.hh:270
TMatrixD covMat_
The calculated covariance matrix.
Definition: LauSPlot.hh:342
Double_t value() const
The value of the parameter.
void addFriendTree()
Add the sWeightTree as a friend tree of the input tree.
Definition: LauSPlot.cc:1246
Bool_t definedSWeightBranches() const
Check whether the sWeights branches have been already created.
Definition: LauSPlot.hh:155
Bool_t definedCNBranches() const
Check whether the cN branches have been already created.
Definition: LauSPlot.hh:143
std::vector< std::map< TString, NumbMap > > sWeights_
The per-event values of the computed sWeights (for each species and for each combination of excluded ...
Definition: LauSPlot.hh:345
Bool_t checkLeaves() const
Check whether the leaf structure makes sense given the PDFs we are expecting.
Definition: LauSPlot.cc:189
Int_t nFreeSpecies_
Number of species free to float in the fit.
Definition: LauSPlot.hh:328
std::set< TString > NameSet
Type to store names, e.g. of the discriminating/control variables.
Definition: LauSPlot.hh:59
void defineSWeightBranches()
Create the branches to store the sWeights.
Definition: LauSPlot.cc:424
Bool_t scfDPSmear() const
Check whether the Self Cross Feed is smeared in the DP.
Definition: LauSPlot.hh:137
TTree * inputTree_
Pointer to the input tree.
Definition: LauSPlot.hh:280
NumbMap origFixdSpecies_
The names and estimated yields of the fixed species - need to keep the original values.
Definition: LauSPlot.hh:301