laura is hosted by Hepforge, IPPP Durham
Laura++  v1r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauFitNtuple.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2004 - 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 
15 #include <iostream>
16 using std::cout;
17 using std::cerr;
18 using std::endl;
19 
20 #include "TFile.h"
21 #include "TMath.h"
22 #include "TSystem.h"
23 #include "TTree.h"
24 #include "TVirtualFitter.h"
25 
26 #include "LauFitNtuple.hh"
27 #include "LauFitter.hh"
28 #include "LauParamFixed.hh"
29 
30 ClassImp(LauFitNtuple)
31 
32 
33 LauFitNtuple::LauFitNtuple(const TString& fileName) :
34  rootFileName_(fileName),
35  rootFile_(0),
36  fitResults_(0),
37  definedFitTree_(kFALSE),
38  fitStatus_(0),
39  nFitPars_(0),
40  nFreePars_(0),
41  nExtraPars_(0),
42  NLL_(0.0),
43  iExpt_(0)
44 {
45  rootFile_ = TFile::Open(rootFileName_, "recreate");
46  rootFile_->cd();
47  fitResults_ = new TTree("fitResults", "fitResults");
48  fitResults_->SetDirectory(rootFile_);
49 
50  fitVars_.clear(); extraVars_.clear();
51 }
52 
54 {
55  if (rootFile_ && rootFile_->IsOpen()) {
56  delete fitResults_; fitResults_ = 0;
57  }
58  delete rootFile_; rootFile_ = 0;
59 }
60 
61 void LauFitNtuple::storeCorrMatrix(UInt_t iExpt, Double_t NLL, Int_t fitStatus)
62 {
63  // store the minimised NLL value, correlation matrix status and experiment number
64  iExpt_ = iExpt;
65  NLL_ = NLL;
66  fitStatus_ = fitStatus;
67 
68  // get a pointer to the fitter
69  TVirtualFitter* fitter = LauFitter::fitter();
70 
71  // make a few checks on the number of parameters
72  Int_t nFitPars = fitter->GetNumberTotalParameters();
73  Int_t nFreePars = fitter->GetNumberFreeParameters();
74  if (nFitPars != nFitPars_) {
75  cerr << "WARNING in LauFitNtuple::storeCorrMatrix : total number of parameters from fitter (" << nFitPars
76  << ") not the same as the number held here (" << nFitPars_ << ")." << endl;
77  }
78  if (nFreePars != nFreePars_) {
79  cerr << "WARNING in LauFitNtuple::storeCorrMatrix : number of free parameters from fitter (" << nFreePars
80  << ") not the same as the number held here (" << nFreePars_ << ")." << endl;
81  }
82 
83  // make the two correlation matrices the correct dimensions
84  if (definedFitTree_ == kFALSE) {
85  globalCC_.clear();
86  globalCC_.resize(nFitPars_);
87 
88  corrMatrix_.clear();
89  corrMatrix_.resize(nFitPars_);
90  for (Int_t i = 0; i < nFitPars_; ++i) {corrMatrix_[i].resize(nFitPars_);}
91  }
92 
93  // get correlation matrix information from the fitter
94  Double_t error(0.0);
95  Double_t negError(0.0);
96  Double_t posError(0.0);
97  Double_t globalcc(0.0);
98  Bool_t iFixed(kFALSE);
99  Bool_t jFixed(kFALSE);
100  Int_t iFree(0);
101  Int_t jFree(0);
102 
103  for (Int_t i = 0; i < nFitPars_; i++) {
104 
105  iFixed = fitVars_[i]->fixed();
106 
107  // get the global correlation factor for each parameter
108  if (iFixed) {
109  globalcc = 0.0;
110  } else {
111  fitter->GetErrors(i, posError, negError, error, globalcc);
112  }
113  globalCC_[i] = globalcc;
114 
115  // reset the "j" free parameter counter
116  jFree = 0;
117 
118  // the GetCovarianceMatrixElement method returns elements of the external error matrix,
119  // which is of dimension nFreePars_ x nFreePars_
120  for (Int_t j = 0; j < nFitPars_; j++) {
121 
122  jFixed = fitVars_[j]->fixed();
123 
124  if (i == j) {
125  corrMatrix_[i][j] = 1.0;
126  } else if (iFixed == kTRUE || jFixed == kTRUE) {
127  corrMatrix_[i][j] = 0.0;
128  } else {
129  Double_t r_ij = fitter->GetCovarianceMatrixElement(iFree,jFree);
130  Double_t r_ii = fitter->GetCovarianceMatrixElement(iFree,iFree);
131  Double_t r_jj = fitter->GetCovarianceMatrixElement(jFree,jFree);
132  Double_t denom = r_ii * r_jj;
133  if (denom < 0.0) {
134  r_ij = 0.0;
135  denom = 1.0;
136  }
137  denom = TMath::Sqrt(denom);
138  if (denom < 1e-30) {
139  r_ij = 0.0;
140  denom = 1.0;
141  }
142  corrMatrix_[i][j] = r_ij / denom;
143  }
144 
145  if (!jFixed) {
146  ++jFree;
147  }
148  }
149 
150  if (!iFixed) {
151  ++iFree;
152  }
153  }
154 }
155 
156 void LauFitNtuple::storeParsAndErrors(const std::vector<LauParameter*>& fitVars, const std::vector<LauParameter>& extraVars)
157 {
158  fitVars_ = fitVars;
159  Int_t nFitPars = fitVars_.size();
160  // the number of parameters being given to us should be the same as the number from the last fit
161  // OR it's the first time so the "last" number is zero
162  if (nFitPars_ != 0 && nFitPars_ != nFitPars) {
163  cerr << "ERROR in LauFitNtuple::storeParsAndErrors : expected total number of parameters (" << nFitPars_
164  << ") not the same as the number provided (" << nFitPars << ")." << endl;
165  gSystem->Exit(EXIT_FAILURE);
166  }
167 
168  LauParamFixed pred;
169  Int_t nFreePars = nFitPars - std::count_if(fitVars_.begin(),fitVars_.end(),pred);
170  // the number of free parameters being given to us should be the same as the number from the last fit
171  // OR it's the first time so the "last" number is zero
172  // (NB we check whether nFitPars_ is zero for this since it is possible to have zero free parameters, albeit rather daft)
173  if (nFitPars_ != 0 && nFreePars_ != nFreePars) {
174  cerr << "ERROR in LauFitNtuple::storeParsAndErrors : expected number of free parameters (" << nFreePars_
175  << ") not the same as the number provided (" << nFreePars << ")." << endl;
176  gSystem->Exit(EXIT_FAILURE);
177  }
178 
179  nFitPars_ = nFitPars;
180  nFreePars_ = nFreePars;
181 
182  extraVars_ = extraVars;
183  nExtraPars_ = extraVars_.size();
184 }
185 
187 {
188  // Now create and fill the stored fit results into an ntuple (TTree)
189  if (definedFitTree_ == kFALSE) {
190 
191  cout << "INFO in LauFitNtuple::updateFitNtuple : totNoPars = " << nFitPars_ << endl;
192 
193  // Add experiment number as a branch
194  fitResults_->Branch("iExpt", &iExpt_, "iExpt/I");
195  fitResults_->Branch("fitStatus", &fitStatus_, "fitStatus/I");
196 
197  // Add NLL (negative log-likelihood) value from fit
198  fitResults_->Branch("NLL", &NLL_, "NLL/D");
199 
200  for (Int_t i = 0; i < nFitPars_; i++) {
201 
202  TString parName = fitVars_[i]->name();
203  TString parNameD(parName); parNameD += "/D";
204  fitResults_->Branch(parName.Data(), &fitVars_[i]->value_, parNameD.Data());
205 
206  TString parInitName(parName); parInitName += "_True";
207  TString parInitNameD(parInitName); parInitNameD += "/D";
208  fitResults_->Branch(parInitName.Data(), &fitVars_[i]->genValue_, parInitNameD.Data());
209 
210  if (!fitVars_[i]->fixed()) {
211  TString parErrName(parName); parErrName += "_Error";
212  TString parErrNameD(parErrName); parErrNameD += "/D";
213  fitResults_->Branch(parErrName.Data(), &fitVars_[i]->error_, parErrNameD.Data());
214 
215  TString parPullName(parName); parPullName += "_Pull";
216  TString parPullNameD(parPullName); parPullNameD += "/D";
217  fitResults_->Branch(parPullName.Data(), &fitVars_[i]->pull_, parPullNameD.Data());
218  }
219 
220  // Now add in the correlation matrix values (only for floating parameters)
221  if (!fitVars_[i]->fixed()) {
222  // First the global correlation coeffs
223  TString parGCCName(parName); parGCCName += "_GCC";
224  TString parGCCNameD(parGCCName); parGCCNameD += "/D";
225  fitResults_->Branch(parGCCName.Data(), &globalCC_[i], parGCCNameD.Data());
226 
227  // Then the rest
228  for (Int_t j = 0; j < nFitPars_; j++) {
229  if (!fitVars_[j]->fixed() && i!=j) {
230  TString parName2 = fitVars_[j]->name();
231  TString corrName("corr__");
232  corrName += parName; corrName += "__"; corrName += parName2;
233 
234  TString corrNameD(corrName); corrNameD += "/D";
235  fitResults_->Branch(corrName.Data(), &corrMatrix_[i][j], corrNameD.Data());
236  }
237  }
238  }
239  }
240 
241  // Update extra parameter values...
242  for (Int_t i = 0; i < nExtraPars_; i++) {
243 
244  TString parName = extraVars_[i].name();
245  TString parNameD(parName); parNameD += "/D";
246  fitResults_->Branch(parName.Data(), &extraVars_[i].value_, parNameD.Data());
247 
248  TString parInitName(parName); parInitName += "_True";
249  TString parInitNameD(parInitName); parInitNameD += "/D";
250  fitResults_->Branch(parInitName.Data(), &extraVars_[i].genValue_, parInitNameD.Data());
251 
252  //TString parErrName(parName); parErrName += "_Error";
253  //TString parErrNameD(parErrName); parErrNameD += "/D";
254  //fitResults_->Branch(parErrName.Data(), &extraVars_[i].error_, parErrNameD.Data());
255 
256  // Also find the fit fraction pull and store it
257  //TString pullName(parName); pullName += "_Pull";
258  //TString pullNameD(pullName); pullNameD += "/D";
259  //fitResults_->Branch(pullName.Data(), &extraVars_[i].pull_, pullNameD.Data());
260 
261  }
262 
263  definedFitTree_ = kTRUE;
264 
265  }
266 
267  cout << "NLL = " << NLL_ << endl;
268 
269  fitResults_->Fill();
270 }
271 
273 {
274  // Write out the fit ntuple to the appropriate root file
275  rootFile_->cd();
276  rootFile_->Write();
277  rootFile_->Close();
278  delete rootFile_; rootFile_ = 0;
279 }
280 
Bool_t definedFitTree_
Flags whether the fit tree has been defined.
Definition: LauFitNtuple.hh:93
void storeCorrMatrix(UInt_t iExpt, Double_t NLL, Int_t fitStatus)
Store the correlation matrix and other fit information.
Definition: LauFitNtuple.cc:61
std::vector< LauParameter * > fitVars_
Fit variables.
Definition: LauFitNtuple.hh:83
Bool_t fixed() const
Check whether the parameter is fixed or floated.
File containing declaration of LauParamFixed class.
Int_t iExpt_
Experiment number.
Int_t nFitPars_
Number of fit parameters.
Definition: LauFitNtuple.hh:97
Double_t negError() const
The lower error on the parameter.
std::vector< LauParameter > extraVars_
Extra variables.
Definition: LauFitNtuple.hh:85
void updateFitNtuple()
Update the fit ntuple.
TFile * rootFile_
Root file.
Definition: LauFitNtuple.hh:78
std::vector< std::vector< Double_t > > corrMatrix_
Correlation matrix.
Definition: LauFitNtuple.hh:90
Double_t posError() const
The upper error on the parameter.
File containing LauFitter namespace.
void writeOutFitResults()
Write out fit results.
Predicate to allow counting of the number of fixed parameters.
Double_t error() const
The error on the parameter.
Double_t NLL_
Minimised negative log likelihood.
virtual ~LauFitNtuple()
Destructor.
Definition: LauFitNtuple.cc:53
Int_t nFreePars_
Number of free parameters.
Definition: LauFitNtuple.hh:99
void storeParsAndErrors(const std::vector< LauParameter * > &fitVars, const std::vector< LauParameter > &extraVars)
Store parameters and their errors.
TTree * fitResults_
Fit results.
Definition: LauFitNtuple.hh:80
Class to store the results from the fit into an ntuple.
Definition: LauFitNtuple.hh:41
Int_t nExtraPars_
Number of extra parameters.
Int_t fitStatus_
Status of fit.
Definition: LauFitNtuple.hh:95
File containing declaration of LauFitNtuple class.
std::vector< Double_t > globalCC_
Global correlation coefficients.
Definition: LauFitNtuple.hh:88
TVirtualFitter * fitter(const TString &fitterString="Minuit", Int_t maxPar=100)
Method that provides a pointer to a TVirtualFitter object.
Definition: LauFitter.cc:20