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