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