laura is hosted by Hepforge, IPPP Durham
Laura++  v3r5
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauFitNtuple.cc
Go to the documentation of this file.
1 
2 /*
3 Copyright 2004 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 
29 #include <iostream>
30 
31 #include "TFile.h"
32 #include "TMath.h"
33 #include "TMatrixD.h"
34 #include "TSystem.h"
35 #include "TTree.h"
36 
37 #include "LauAbsFitter.hh"
38 #include "LauFitNtuple.hh"
39 #include "LauFitter.hh"
40 #include "LauParameter.hh"
41 #include "LauParamFixed.hh"
42 
44 
45 
46 LauFitNtuple::LauFitNtuple(const TString& fileName, Bool_t storeAsymErrors) :
47  rootFileName_(fileName),
48  rootFile_(0),
49  fitResults_(0),
50  definedFitTree_(kFALSE),
51  storeAsymErrors_(storeAsymErrors),
52  fitStatus_({-1,0.0,0.0}),
53  nFitPars_(0),
54  nFreePars_(0),
55  nExtraPars_(0),
56  iExpt_(0)
57 {
58  rootFile_ = TFile::Open(rootFileName_, "recreate");
59  rootFile_->cd();
60  fitResults_ = new TTree("fitResults", "fitResults");
61  fitResults_->SetDirectory(rootFile_);
62 
63  fitVars_.clear(); extraVars_.clear();
64 }
65 
67 {
68  if (rootFile_ && rootFile_->IsOpen()) {
69  delete fitResults_; fitResults_ = 0;
70  }
71  delete rootFile_; rootFile_ = 0;
72 }
73 
74 void LauFitNtuple::storeCorrMatrix(const UInt_t iExpt, const LauAbsFitter::FitStatus& fitStatus, const TMatrixD& covMatrix)
75 {
76  // store the minimised NLL value, correlation matrix status and experiment number
77  iExpt_ = iExpt;
78  fitStatus_ = fitStatus;
79 
80  // make the correlation matrix the correct dimensions
81  if (definedFitTree_ == kFALSE) {
82  corrMatrix_.clear();
83  corrMatrix_.resize(nFitPars_);
84  for (UInt_t i = 0; i < nFitPars_; ++i) {corrMatrix_[i].resize(nFitPars_);}
85  }
86 
87  // under certain circumstances, e.g. if the fit has failed in the first
88  // stage of a two-stage fit, the covariance matrix might not have the
89  // expected dimensions, or it might even be empty
90  Bool_t needsPadding = kFALSE;
91  const UInt_t nElements = covMatrix.GetNoElements();
92  if ( nElements == 0 ) {
93 
94  // if it's empty we can just make a diagonal matrix
95  std::cerr << "WARNING in LauFitNtuple::storeCorrMatrix : received empty covariance matrix - will store diagonal correlation matrix" << std::endl;
96  for (UInt_t i = 0; i < nFitPars_; ++i) {
97  for (UInt_t j = 0; j < nFitPars_; ++j) {
98  if (i == j) {
99  corrMatrix_[i][j] = 1.0;
100  } else {
101  corrMatrix_[i][j] = 0.0;
102  }
103  }
104  }
105  return;
106 
107  } else if ( nElements != nFreePars_*nFreePars_ ) {
108 
109  UInt_t dimension = covMatrix.GetNrows();
110  UInt_t nSecondStage = 0;
111  for (UInt_t i = 0; i < nFitPars_; ++i) {
112  if ( fitVars_[i]->secondStage() ) {
113  ++nSecondStage;
114  }
115  }
116 
117  if ( (dimension + nSecondStage) == nFreePars_ ) {
118  needsPadding = kTRUE;
119  std::cerr << "WARNING in LauFitNtuple::storeCorrMatrix : received smaller covariance matrix than expected, likely due to failure of first stage fit - will pad the correlation matrix" << std::endl;
120  std::cerr << " : nFitPars_ = " << nFitPars_ << std::endl;
121  std::cerr << " : nFreePars_ = " << nFreePars_ << std::endl;
122  std::cerr << " : nSecondStage = " << nSecondStage << std::endl;
123  std::cerr << " : covMatrix size = " << dimension << std::endl;
124  } else {
125  std::cerr << "WARNING in LauFitNtuple::storeCorrMatrix : received smaller covariance matrix than expected, for unknown reasons - will store diagonal correlation matrix" << std::endl;
126  std::cerr << " : nFitPars_ = " << nFitPars_ << std::endl;
127  std::cerr << " : nFreePars_ = " << nFreePars_ << std::endl;
128  std::cerr << " : nSecondStage = " << nSecondStage << std::endl;
129  std::cerr << " : covMatrix size = " << dimension << std::endl;
130  for (UInt_t i = 0; i < nFitPars_; ++i) {
131  for (UInt_t j = 0; j < nFitPars_; ++j) {
132  if (i == j) {
133  corrMatrix_[i][j] = 1.0;
134  } else {
135  corrMatrix_[i][j] = 0.0;
136  }
137  }
138  }
139  return;
140  }
141 
142  }
143 
144  // calculate the correlation matrix information from the fitter covariance matrix
145  Bool_t iFixed(kFALSE);
146  Bool_t jFixed(kFALSE);
147  Bool_t iSecondStage(kFALSE);
148  Bool_t jSecondStage(kFALSE);
149  UInt_t iFree(0);
150  UInt_t jFree(0);
151 
152  for (UInt_t i = 0; i < nFitPars_; ++i) {
153 
154  iFixed = fitVars_[i]->fixed();
155  iSecondStage = fitVars_[i]->secondStage();
156 
157  // reset the "j" free parameter counter
158  jFree = 0;
159 
160  // NB the supplied covariance matrix is of dimension nFreePars_ x nFreePars_
161  for (UInt_t j = 0; j < nFitPars_; ++j) {
162 
163  jFixed = fitVars_[j]->fixed();
164  jSecondStage = fitVars_[j]->secondStage();
165 
166  if (i == j) {
167  corrMatrix_[i][j] = 1.0;
168  } else if (iFixed == kTRUE || jFixed == kTRUE) {
169  corrMatrix_[i][j] = 0.0;
170  } else if ( needsPadding && ( iSecondStage || jSecondStage ) ) {
171  corrMatrix_[i][j] = 0.0;
172  } else {
173  Double_t r_ij = covMatrix(iFree,jFree);
174  Double_t r_ii = covMatrix(iFree,iFree);
175  Double_t r_jj = covMatrix(jFree,jFree);
176  Double_t denom = r_ii * r_jj;
177  if (denom < 0.0) {
178  r_ij = 0.0;
179  denom = 1.0;
180  }
181  denom = TMath::Sqrt(denom);
182  if (denom < 1e-30) {
183  r_ij = 0.0;
184  denom = 1.0;
185  }
186  corrMatrix_[i][j] = r_ij / denom;
187  }
188 
189  if ( !jFixed && !(needsPadding && jSecondStage) ) {
190  ++jFree;
191  }
192  }
193 
194  if ( !iFixed && !(needsPadding && iSecondStage) ) {
195  ++iFree;
196  }
197  }
198 }
199 
200 void LauFitNtuple::storeParsAndErrors(const std::vector<LauParameter*>& fitVars, const std::vector<LauParameter>& extraVars)
201 {
202  fitVars_ = fitVars;
203  UInt_t nFitPars = fitVars_.size();
204 
205  // the number of parameters being given to us should be the same as the number from the last fit
206  // OR it's the first time so the "last" number is zero
207  if (nFitPars_ != 0 && nFitPars_ != nFitPars) {
208  std::cerr << "ERROR in LauFitNtuple::storeParsAndErrors : expected total number of parameters (" << nFitPars_
209  << ") not the same as the number provided (" << nFitPars << ")." << std::endl;
210  gSystem->Exit(EXIT_FAILURE);
211  }
212 
213  LauParamFixed pred;
214  UInt_t nFreePars = nFitPars - std::count_if(fitVars_.begin(),fitVars_.end(),pred);
215 
216  // the number of free parameters being given to us should be the same as the number from the last fit
217  // OR it's the first time so the "last" number is zero
218  // (NB we check whether nFitPars_ is zero for this since it is possible to have zero free parameters, albeit rather daft)
219  if (nFitPars_ != 0 && nFreePars_ != nFreePars) {
220  std::cerr << "ERROR in LauFitNtuple::storeParsAndErrors : expected number of free parameters (" << nFreePars_
221  << ") not the same as the number provided (" << nFreePars << ")." << std::endl;
222  gSystem->Exit(EXIT_FAILURE);
223  }
224 
225  nFitPars_ = nFitPars;
226  nFreePars_ = nFreePars;
227 
228  extraVars_ = extraVars;
229  nExtraPars_ = extraVars_.size();
230 }
231 
233 {
234  // Now create and fill the stored fit results into an ntuple (TTree)
235  if (definedFitTree_ == kFALSE) {
236 
237  std::cout << "INFO in LauFitNtuple::updateFitNtuple : totNoPars = " << nFitPars_ << std::endl;
238 
239  // Add experiment number as a branch
240  fitResults_->Branch("iExpt", &iExpt_, "iExpt/I");
241  fitResults_->Branch("fitStatus", &fitStatus_.status, "fitStatus/I");
242 
243  // Add NLL (negative log-likelihood) and EDM values from fit
244  fitResults_->Branch("NLL", &fitStatus_.NLL, "NLL/D");
245  fitResults_->Branch("EDM", &fitStatus_.EDM, "EDM/D");
246 
247  for (UInt_t i = 0; i < nFitPars_; i++) {
248 
249  TString parName = fitVars_[i]->name();
250  TString parNameD(parName); parNameD += "/D";
251  fitResults_->Branch(parName.Data(), &fitVars_[i]->value_, parNameD.Data());
252 
253  TString parInitName(parName); parInitName += "_True";
254  TString parInitNameD(parInitName); parInitNameD += "/D";
255  fitResults_->Branch(parInitName.Data(), &fitVars_[i]->genValue_, parInitNameD.Data());
256 
257  if (!fitVars_[i]->fixed()) {
258  TString parErrName(parName); parErrName += "_Error";
259  TString parErrNameD(parErrName); parErrNameD += "/D";
260  fitResults_->Branch(parErrName.Data(), &fitVars_[i]->error_, parErrNameD.Data());
261 
262  if ( storeAsymErrors_ ) {
263  TString parNegErrName(parName); parNegErrName += "_NegError";
264  TString parNegErrNameD(parNegErrName); parNegErrNameD += "/D";
265  fitResults_->Branch(parNegErrName.Data(), &fitVars_[i]->negError_, parNegErrNameD.Data());
266 
267  TString parPosErrName(parName); parPosErrName += "_PosError";
268  TString parPosErrNameD(parPosErrName); parPosErrNameD += "/D";
269  fitResults_->Branch(parPosErrName.Data(), &fitVars_[i]->posError_, parPosErrNameD.Data());
270  }
271 
272  TString parPullName(parName); parPullName += "_Pull";
273  TString parPullNameD(parPullName); parPullNameD += "/D";
274  fitResults_->Branch(parPullName.Data(), &fitVars_[i]->pull_, parPullNameD.Data());
275  }
276 
277  // Now add in the correlation matrix values (only for floating parameters)
278  if (!fitVars_[i]->fixed()) {
279  // First the global correlation coeffs
280  TString parGCCName(parName); parGCCName += "_GCC";
281  TString parGCCNameD(parGCCName); parGCCNameD += "/D";
282  fitResults_->Branch(parGCCName.Data(), &fitVars_[i]->gcc_, parGCCNameD.Data());
283 
284  if ( ! corrMatrix_.empty() ) {
285  // Then the rest
286  for (UInt_t j = 0; j < nFitPars_; j++) {
287  if (!fitVars_[j]->fixed() && i!=j) {
288  TString parName2 = fitVars_[j]->name();
289  TString corrName("corr__");
290  corrName += parName; corrName += "__"; corrName += parName2;
291 
292  TString corrNameD(corrName); corrNameD += "/D";
293  fitResults_->Branch(corrName.Data(), &corrMatrix_[i][j], corrNameD.Data());
294  }
295  }
296  }
297  }
298  }
299 
300  // Update extra parameter values...
301  for (UInt_t i = 0; i < nExtraPars_; i++) {
302 
303  TString parName = extraVars_[i].name();
304  TString parNameD(parName); parNameD += "/D";
305  fitResults_->Branch(parName.Data(), &extraVars_[i].value_, parNameD.Data());
306 
307  TString parInitName(parName); parInitName += "_True";
308  TString parInitNameD(parInitName); parInitNameD += "/D";
309  fitResults_->Branch(parInitName.Data(), &extraVars_[i].genValue_, parInitNameD.Data());
310 
311  //TString parErrName(parName); parErrName += "_Error";
312  //TString parErrNameD(parErrName); parErrNameD += "/D";
313  //fitResults_->Branch(parErrName.Data(), &extraVars_[i].error_, parErrNameD.Data());
314 
315  // Also find the fit fraction pull and store it
316  //TString pullName(parName); pullName += "_Pull";
317  //TString pullNameD(pullName); pullNameD += "/D";
318  //fitResults_->Branch(pullName.Data(), &extraVars_[i].pull_, pullNameD.Data());
319 
320  }
321 
322  definedFitTree_ = kTRUE;
323 
324  }
325 
326  std::cout << "INFO in LauFitNtuple::updateFitNtuple : Stored fit values for experiment " << iExpt_ << "\n";
327  std::cout << " : fitStatus = " << fitStatus_.status << std::endl;
328  std::cout << " : NLL = " << fitStatus_.NLL << std::endl;
329  std::cout << " : EDM = " << fitStatus_.EDM << std::endl;
330 
331  fitResults_->Fill();
332 }
333 
335 {
336  // Write out the fit ntuple to the appropriate root file
337  rootFile_->cd();
338  fitResults_->Write("",TObject::kOverwrite);
339  rootFile_->Close();
340  delete rootFile_; rootFile_ = 0;
341 }
342 
Bool_t definedFitTree_
Flags whether the fit tree has been defined.
std::vector< LauParameter * > fitVars_
Fit variables.
Double_t EDM
The estimated distance to the minimum.
Definition: LauAbsFitter.hh:61
Bool_t fixed() const
Check whether the parameter is fixed or floated.
File containing declaration of LauParamFixed class.
nExtraPars_(0)
Int_t iExpt_
Experiment number.
ClassImp(LauAbsCoeffSet)
UInt_t nExtraPars_
Number of extra parameters.
File containing declaration of LauAbsFitter class.
LauAbsFitter::FitStatus fitStatus_
Status of fit.
nFitPars_(0)
std::vector< std::vector< Double_t > > corrMatrix_
Correlation matrix.
std::vector< LauParameter > extraVars_
Extra variables.
void updateFitNtuple()
Update the fit ntuple.
TFile * rootFile_
Root file.
File containing declaration of LauFitter class.
UInt_t nFitPars_
Number of fit parameters.
void storeCorrMatrix(const UInt_t iExpt, const LauAbsFitter::FitStatus &fitStatus, const TMatrixD &covMatrix)
Store the correlation matrix and other fit information.
Definition: LauFitNtuple.cc:74
void writeOutFitResults()
Write out fit results.
nFreePars_(0)
Double_t NLL
The negative log-likelihood.
Definition: LauAbsFitter.hh:59
Int_t status
The status code of the fit.
Definition: LauAbsFitter.hh:57
File containing declaration of LauParameter class.
Struct to store fit status information.
Definition: LauAbsFitter.hh:55
Predicate to allow counting of the number of fixed parameters.
Bool_t secondStage() const
Check whether the parameter should be floated only in the second stage of a two stage fit...
virtual ~LauFitNtuple()
Destructor.
Definition: LauFitNtuple.cc:66
void storeParsAndErrors(const std::vector< LauParameter * > &fitVars, const std::vector< LauParameter > &extraVars)
Store parameters and their errors.
TTree * fitResults_
Fit results.
Class to store the results from the fit into an ntuple.
Definition: LauFitNtuple.hh:58
iExpt_(0)
Definition: LauFitNtuple.cc:56
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.