]>
Commit | Line | Data |
---|---|---|
c0b10ad4 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
7bc6a013 | 16 | //---------------------------------------------------------------------// |
17 | // // | |
18 | // AliCFUnfolding Class // | |
19 | // Class to handle general unfolding procedure // | |
20 | // For the moment only bayesian unfolding is supported // | |
21 | // The next steps are to add chi2 minimisation and weighting methods // | |
22 | // // | |
23 | // // | |
24 | // // | |
25 | // Use : // | |
26 | // ------- // | |
27 | // The Bayesian unfolding consists of several iterations. // | |
28 | // At each iteration, an inverse response matrix is calculated, given // | |
29 | // the measured spectrum, the a priori (guessed) spectrum, // | |
30 | // the efficiency spectrum and the response matrix. // | |
c7f196bd | 31 | // // |
32 | // Then at each iteration, the unfolded spectrum is calculated using // | |
7bc6a013 | 33 | // the inverse response : the goal is to get an unfolded spectrum // |
34 | // similar (according to some criterion) to the a priori one. // | |
35 | // If the difference is too big, another iteration is performed : // | |
36 | // the a priori spectrum is updated to the unfolded one from the // | |
37 | // previous iteration, and so on so forth, until the maximum number // | |
38 | // of iterations or the similarity criterion is reached. // | |
39 | // // | |
40 | // Currently the similarity criterion is the Chi2 between the a priori // | |
41 | // and the unfolded spectrum. // | |
42 | // // | |
43 | // Currently the user has to define the max. number of iterations // | |
44 | // (::SetMaxNumberOfIterations) // | |
45 | // and the chi2 below which the procedure will stop // | |
46 | // (::SetMaxChi2 or ::SetMaxChi2PerDOF) // | |
47 | // // | |
48 | // An optional possibility is to smooth the unfolded spectrum at the // | |
49 | // end of each iteration, either using a fit function // | |
50 | // (only if #dimensions <=3) // | |
51 | // or a simple averaging using the neighbouring bins values. // | |
52 | // This is possible calling the function ::UseSmoothing // | |
53 | // If no argument is passed to this function, then the second option // | |
54 | // is used. // | |
55 | // // | |
c7f196bd | 56 | // IMPORTANT: // |
57 | //----------- // | |
58 | // With this approach, the efficiency map must be calculated // | |
59 | // with *simulated* values only, otherwise the method won't work. // | |
60 | // // | |
61 | // ex: efficiency(bin_pt) = number_rec(bin_pt) / number_sim(bin_pt) // | |
62 | // // | |
63 | // the pt bin "bin_pt" must always be the same in both the efficiency // | |
64 | // numerator and denominator. // | |
65 | // This is why the efficiency map has to be created by a method // | |
66 | // from which both reconstructed and simulated values are accessible // | |
67 | // simultaneously. // | |
68 | // // | |
69 | // // | |
7bc6a013 | 70 | //---------------------------------------------------------------------// |
71 | // Author : renaud.vernet@cern.ch // | |
72 | //---------------------------------------------------------------------// | |
c0b10ad4 | 73 | |
74 | ||
75 | #include "AliCFUnfolding.h" | |
76 | #include "TMath.h" | |
77 | #include "TAxis.h" | |
78 | #include "AliLog.h" | |
85b6bda9 | 79 | #include "TF1.h" |
80 | #include "TH1D.h" | |
81 | #include "TH2D.h" | |
82 | #include "TH3D.h" | |
c0b10ad4 | 83 | |
84 | ClassImp(AliCFUnfolding) | |
85 | ||
86 | //______________________________________________________________ | |
87 | ||
88 | AliCFUnfolding::AliCFUnfolding() : | |
89 | TNamed(), | |
90 | fResponse(0x0), | |
91 | fPrior(0x0), | |
c0b10ad4 | 92 | fEfficiency(0x0), |
93 | fMeasured(0x0), | |
94 | fMaxNumIterations(0), | |
95 | fNVariables(0), | |
96 | fMaxChi2(0), | |
97 | fUseSmoothing(kFALSE), | |
85b6bda9 | 98 | fSmoothFunction(0x0), |
99 | fSmoothOption(""), | |
100 | fOriginalPrior(0x0), | |
c0b10ad4 | 101 | fInverseResponse(0x0), |
102 | fMeasuredEstimate(0x0), | |
103 | fConditional(0x0), | |
104 | fProjResponseInT(0x0), | |
105 | fUnfolded(0x0), | |
106 | fCoordinates2N(0x0), | |
107 | fCoordinatesN_M(0x0), | |
108 | fCoordinatesN_T(0x0) | |
109 | { | |
110 | // | |
111 | // default constructor | |
112 | // | |
113 | } | |
114 | ||
115 | //______________________________________________________________ | |
116 | ||
117 | AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar, | |
118 | const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior) : | |
119 | TNamed(name,title), | |
120 | fResponse((THnSparse*)response->Clone()), | |
121 | fPrior(0x0), | |
c0b10ad4 | 122 | fEfficiency((THnSparse*)efficiency->Clone()), |
123 | fMeasured((THnSparse*)measured->Clone()), | |
124 | fMaxNumIterations(0), | |
125 | fNVariables(nVar), | |
126 | fMaxChi2(0), | |
127 | fUseSmoothing(kFALSE), | |
85b6bda9 | 128 | fSmoothFunction(0x0), |
129 | fSmoothOption(""), | |
130 | fOriginalPrior(0x0), | |
c0b10ad4 | 131 | fInverseResponse(0x0), |
132 | fMeasuredEstimate(0x0), | |
133 | fConditional(0x0), | |
134 | fProjResponseInT(0x0), | |
135 | fUnfolded(0x0), | |
136 | fCoordinates2N(0x0), | |
137 | fCoordinatesN_M(0x0), | |
138 | fCoordinatesN_T(0x0) | |
139 | { | |
140 | // | |
141 | // named constructor | |
142 | // | |
143 | ||
144 | AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions())); | |
145 | ||
146 | if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution | |
147 | else { | |
148 | fPrior = (THnSparse*) prior->Clone(); | |
149 | fOriginalPrior = (THnSparse*)fPrior->Clone(); | |
85b6bda9 | 150 | if (fPrior->GetNdimensions() != fNVariables) |
151 | AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions())); | |
c0b10ad4 | 152 | } |
85b6bda9 | 153 | |
154 | if (fEfficiency->GetNdimensions() != fNVariables) | |
155 | AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions())); | |
156 | if (fMeasured->GetNdimensions() != fNVariables) | |
157 | AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions())); | |
158 | if (fResponse->GetNdimensions() != 2*fNVariables) | |
159 | AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions())); | |
c0b10ad4 | 160 | |
85b6bda9 | 161 | |
c0b10ad4 | 162 | for (Int_t iVar=0; iVar<fNVariables; iVar++) { |
163 | AliInfo(Form("prior matrix has %d bins in dimension %d",fPrior ->GetAxis(iVar)->GetNbins(),iVar)); | |
164 | AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar)); | |
165 | AliInfo(Form("measured matrix has %d bins in dimension %d",fMeasured ->GetAxis(iVar)->GetNbins(),iVar)); | |
166 | } | |
167 | Init(); | |
168 | } | |
169 | ||
170 | ||
171 | //______________________________________________________________ | |
172 | ||
173 | AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) : | |
174 | TNamed(c), | |
175 | fResponse((THnSparse*)c.fResponse->Clone()), | |
176 | fPrior((THnSparse*)c.fPrior->Clone()), | |
c0b10ad4 | 177 | fEfficiency((THnSparse*)c.fEfficiency->Clone()), |
178 | fMeasured((THnSparse*)c.fMeasured->Clone()), | |
179 | fMaxNumIterations(c.fMaxNumIterations), | |
180 | fNVariables(c.fNVariables), | |
181 | fMaxChi2(c.fMaxChi2), | |
182 | fUseSmoothing(c.fUseSmoothing), | |
85b6bda9 | 183 | fSmoothFunction((TF1*)c.fSmoothFunction->Clone()), |
184 | fSmoothOption(fSmoothOption), | |
185 | fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()), | |
c0b10ad4 | 186 | fInverseResponse((THnSparse*)c.fInverseResponse->Clone()), |
187 | fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()), | |
188 | fConditional((THnSparse*)c.fConditional->Clone()), | |
189 | fProjResponseInT((THnSparse*)c.fProjResponseInT->Clone()), | |
190 | fUnfolded((THnSparse*)c.fUnfolded->Clone()), | |
191 | fCoordinates2N(new Int_t(*c.fCoordinates2N)), | |
192 | fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)), | |
193 | fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T)) | |
194 | { | |
195 | // | |
196 | // copy constructor | |
197 | // | |
198 | } | |
199 | ||
200 | //______________________________________________________________ | |
201 | ||
202 | AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) { | |
203 | // | |
204 | // assignment operator | |
205 | // | |
206 | ||
207 | if (this!=&c) { | |
208 | TNamed::operator=(c); | |
209 | fResponse = (THnSparse*)c.fResponse->Clone() ; | |
210 | fPrior = (THnSparse*)c.fPrior->Clone() ; | |
c0b10ad4 | 211 | fEfficiency = (THnSparse*)c.fEfficiency->Clone() ; |
212 | fMeasured = (THnSparse*)c.fMeasured->Clone() ; | |
213 | fMaxNumIterations = c.fMaxNumIterations ; | |
214 | fNVariables = c.fNVariables ; | |
215 | fMaxChi2 = c.fMaxChi2 ; | |
216 | fUseSmoothing = c.fUseSmoothing ; | |
85b6bda9 | 217 | fSmoothFunction = (TF1*)c.fSmoothFunction->Clone(); |
218 | fSmoothOption = c.fSmoothOption ; | |
219 | fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ; | |
c0b10ad4 | 220 | fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ; |
221 | fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ; | |
222 | fConditional = (THnSparse*)c.fConditional->Clone() ; | |
223 | fProjResponseInT = (THnSparse*)c.fProjResponseInT->Clone() ; | |
224 | fUnfolded = (THnSparse*)c.fUnfolded->Clone() ; | |
225 | fCoordinates2N = new Int_t(*c.fCoordinates2N) ; | |
226 | fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ; | |
227 | fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ; | |
228 | } | |
229 | return *this; | |
230 | } | |
231 | ||
232 | //______________________________________________________________ | |
233 | ||
234 | AliCFUnfolding::~AliCFUnfolding() { | |
235 | // | |
236 | // destructor | |
237 | // | |
238 | if (fResponse) delete fResponse; | |
239 | if (fPrior) delete fPrior; | |
c0b10ad4 | 240 | if (fEfficiency) delete fEfficiency; |
241 | if (fMeasured) delete fMeasured; | |
85b6bda9 | 242 | if (fSmoothFunction) delete fSmoothFunction; |
243 | if (fOriginalPrior) delete fOriginalPrior; | |
c0b10ad4 | 244 | if (fInverseResponse) delete fInverseResponse; |
245 | if (fMeasuredEstimate) delete fMeasuredEstimate; | |
246 | if (fConditional) delete fConditional; | |
247 | if (fProjResponseInT) delete fProjResponseInT; | |
248 | if (fCoordinates2N) delete [] fCoordinates2N; | |
249 | if (fCoordinatesN_M) delete [] fCoordinatesN_M; | |
250 | if (fCoordinatesN_T) delete [] fCoordinatesN_T; | |
251 | } | |
252 | ||
253 | //______________________________________________________________ | |
254 | ||
255 | void AliCFUnfolding::Init() { | |
256 | // | |
257 | // initialisation function : creates internal settings | |
258 | // | |
259 | ||
260 | fCoordinates2N = new Int_t[2*fNVariables]; | |
261 | fCoordinatesN_M = new Int_t[fNVariables]; | |
262 | fCoordinatesN_T = new Int_t[fNVariables]; | |
263 | ||
264 | // create the matrix of conditional probabilities P(M|T) | |
265 | CreateConditional(); | |
266 | ||
267 | // create the frame of the inverse response matrix | |
268 | fInverseResponse = (THnSparse*) fResponse->Clone(); | |
269 | // create the frame of the unfolded spectrum | |
270 | fUnfolded = (THnSparse*) fPrior->Clone(); | |
271 | // create the frame of the measurement estimate spectrum | |
272 | fMeasuredEstimate = (THnSparse*) fMeasured->Clone(); | |
273 | } | |
274 | ||
275 | //______________________________________________________________ | |
276 | ||
277 | void AliCFUnfolding::CreateEstMeasured() { | |
278 | // | |
279 | // This function creates a estimate (M) of the reconstructed spectrum | |
7bc6a013 | 280 | // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND) |
c0b10ad4 | 281 | // |
282 | // --> P(M) = SUM { P(M|T) * P(T) } | |
7bc6a013 | 283 | // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)} |
c0b10ad4 | 284 | // |
285 | // This is needed to calculate the inverse response matrix | |
286 | // | |
287 | ||
288 | ||
289 | // clean the measured estimate spectrum | |
85f9f9e1 | 290 | fMeasuredEstimate->Reset(); |
291 | ||
85b6bda9 | 292 | THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone(); |
293 | priorTimesEff->Multiply(fEfficiency); | |
294 | ||
c0b10ad4 | 295 | // fill it |
7bc6a013 | 296 | for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) { |
c0b10ad4 | 297 | Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N); |
85b6bda9 | 298 | Double_t conditionalError = fConditional->GetBinError (iBin); |
c0b10ad4 | 299 | GetCoordinates(); |
85b6bda9 | 300 | Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T); |
301 | Double_t priorTimesEffError = priorTimesEff->GetBinError (fCoordinatesN_T); | |
302 | Double_t fill = conditionalValue * priorTimesEffValue ; | |
303 | ||
304 | if (fill>0.) { | |
85b6bda9 | 305 | // error calculation : gaussian error propagation (may be overestimated...) |
306 | Double_t err2 = TMath::Power(fMeasuredEstimate->GetBinError(fCoordinatesN_M),2) ; | |
307 | err2 += TMath::Power(conditionalValue*priorTimesEffError,2) + TMath::Power(conditionalError*priorTimesEffValue,2) ; | |
308 | Double_t err = TMath::Sqrt(err2); | |
309 | fMeasuredEstimate->SetBinError(fCoordinatesN_M,err); | |
12e419d5 | 310 | |
311 | fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill); | |
85b6bda9 | 312 | } |
c0b10ad4 | 313 | } |
85b6bda9 | 314 | delete priorTimesEff ; |
c0b10ad4 | 315 | } |
316 | ||
317 | //______________________________________________________________ | |
318 | ||
319 | void AliCFUnfolding::CreateInvResponse() { | |
320 | // | |
321 | // Creates the inverse response matrix (INV) with Bayesian method | |
322 | // : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E) | |
323 | // | |
324 | // --> P(T|M) = P(M|T) * P(T) * eff(T) / SUM { P(M|T) * P(T) } | |
325 | // --> INV(i,j) = COND(i,j) * T(j) * E(j) / SUM_k { COND(i,k) * T(k) } | |
326 | // | |
327 | ||
85b6bda9 | 328 | THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone(); |
329 | priorTimesEff->Multiply(fEfficiency); | |
330 | ||
7bc6a013 | 331 | for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) { |
c0b10ad4 | 332 | Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N); |
85b6bda9 | 333 | Double_t conditionalError = fConditional->GetBinError (iBin); |
c0b10ad4 | 334 | GetCoordinates(); |
85b6bda9 | 335 | Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M); |
336 | Double_t estMeasuredError = fMeasuredEstimate->GetBinError (fCoordinatesN_M); | |
337 | Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T); | |
338 | Double_t priorTimesEffError = priorTimesEff ->GetBinError (fCoordinatesN_T); | |
339 | Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ; | |
340 | // error calculation : gaussian error propagation (may be overestimated...) | |
341 | Double_t err = 0. ; | |
342 | if (estMeasuredValue>0.) { | |
343 | err = TMath::Sqrt( TMath::Power(conditionalError * priorTimesEffValue * estMeasuredValue ,2) + | |
344 | TMath::Power(conditionalValue * priorTimesEffError * estMeasuredValue ,2) + | |
345 | TMath::Power(conditionalValue * priorTimesEffValue * estMeasuredError ,2) ) | |
346 | / TMath::Power(estMeasuredValue,2) ; | |
347 | } | |
348 | if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) { | |
349 | fInverseResponse->SetBinContent(fCoordinates2N,fill); | |
350 | fInverseResponse->SetBinError (fCoordinates2N,err ); | |
351 | } | |
352 | } | |
353 | delete priorTimesEff ; | |
c0b10ad4 | 354 | } |
355 | ||
356 | //______________________________________________________________ | |
357 | ||
358 | void AliCFUnfolding::Unfold() { | |
359 | // | |
360 | // Main routine called by the user : | |
361 | // it calculates the unfolded spectrum from the response matrix and the measured spectrum | |
362 | // several iterations are performed until a reasonable chi2 is reached | |
363 | // | |
364 | ||
365 | Int_t iIterBayes=0 ; | |
366 | Double_t chi2=0 ; | |
367 | ||
368 | for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations | |
369 | CreateEstMeasured(); | |
370 | CreateInvResponse(); | |
371 | CreateUnfolded(); | |
372 | chi2 = GetChi2(); | |
85f9f9e1 | 373 | AliDebug(0,Form("Chi2 at iteration %d is %e",iIterBayes,chi2)); |
c0b10ad4 | 374 | if (fMaxChi2>0. && chi2<fMaxChi2) { |
375 | break; | |
376 | } | |
377 | // update the prior distribution | |
85b6bda9 | 378 | if (fUseSmoothing) { |
379 | if (Smooth()) { | |
380 | AliError("Couldn't smooth the unfolded spectrum!!"); | |
7bc6a013 | 381 | AliInfo(Form("\n\n=======================\nFinished at iteration %d : Chi2 is %e and you required it to be < %e\n=======================\n\n",iIterBayes,chi2,fMaxChi2)); |
85b6bda9 | 382 | return; |
383 | } | |
384 | } | |
c0b10ad4 | 385 | fPrior = (THnSparse*)fUnfolded->Clone() ; // this should be changed (memory) |
386 | } | |
85b6bda9 | 387 | AliInfo(Form("\n\n=======================\nFinished at iteration %d : Chi2 is %e and you required it to be < %e\n=======================\n\n",iIterBayes,chi2,fMaxChi2)); |
c0b10ad4 | 388 | } |
389 | ||
390 | //______________________________________________________________ | |
391 | ||
392 | void AliCFUnfolding::CreateUnfolded() { | |
393 | // | |
394 | // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV) | |
395 | // We have P(T) = SUM { P(T|M) * P(M) } | |
396 | // --> T(i) = SUM_k { INV(i,k) * M(k) } | |
397 | // | |
398 | ||
399 | ||
400 | // clear the unfolded spectrum | |
85f9f9e1 | 401 | fUnfolded->Reset(); |
c0b10ad4 | 402 | |
7bc6a013 | 403 | for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) { |
c0b10ad4 | 404 | Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N); |
85b6bda9 | 405 | Double_t invResponseError = fInverseResponse->GetBinError (iBin); |
c0b10ad4 | 406 | GetCoordinates(); |
85b6bda9 | 407 | Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T); |
408 | Double_t effError = fEfficiency->GetBinError (fCoordinatesN_T); | |
409 | Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M); | |
410 | Double_t measuredError = fMeasured ->GetBinError (fCoordinatesN_M); | |
411 | Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ; | |
412 | ||
413 | if (fill>0.) { | |
85b6bda9 | 414 | // error calculation : gaussian error propagation (may be overestimated...) |
415 | Double_t err2 = TMath::Power(fUnfolded->GetBinError(fCoordinatesN_T),2) ; | |
416 | err2 += TMath::Power(invResponseError * measuredValue * effValue,2) / TMath::Power(effValue,4) ; | |
417 | err2 += TMath::Power(invResponseValue * measuredError * effValue,2) / TMath::Power(effValue,4) ; | |
418 | err2 += TMath::Power(invResponseValue * measuredValue * effError,2) / TMath::Power(effValue,4) ; | |
419 | Double_t err = TMath::Sqrt(err2); | |
420 | fUnfolded->SetBinError(fCoordinatesN_T,err); | |
12e419d5 | 421 | |
422 | fUnfolded->AddBinContent(fCoordinatesN_T,fill); | |
85b6bda9 | 423 | } |
c0b10ad4 | 424 | } |
425 | } | |
426 | ||
85b6bda9 | 427 | //______________________________________________________________ |
428 | ||
c0b10ad4 | 429 | void AliCFUnfolding::GetCoordinates() { |
430 | // | |
431 | // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N) | |
432 | // | |
433 | for (Int_t i = 0; i<fNVariables ; i++) { | |
434 | fCoordinatesN_M[i] = fCoordinates2N[i]; | |
435 | fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables]; | |
436 | } | |
437 | } | |
438 | ||
439 | //______________________________________________________________ | |
440 | ||
441 | void AliCFUnfolding::CreateConditional() { | |
442 | // | |
443 | // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R | |
444 | // | |
445 | // --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) } | |
446 | // | |
447 | ||
448 | fConditional = (THnSparse*) fResponse->Clone(); // output of this function | |
449 | fProjResponseInT = (THnSparse*) fPrior->Clone(); // output denominator : | |
450 | // projection of the response matrix on the TRUE axis | |
85b6bda9 | 451 | Int_t* dim = new Int_t [fNVariables]; |
452 | for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1) | |
453 | fProjResponseInT = fConditional->Projection(fNVariables,dim,"E"); //project | |
454 | delete [] dim; | |
c0b10ad4 | 455 | |
456 | // fill the conditional probability matrix | |
7bc6a013 | 457 | for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) { |
c0b10ad4 | 458 | Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N); |
85b6bda9 | 459 | Double_t responseError = fResponse->GetBinError (iBin); |
c0b10ad4 | 460 | GetCoordinates(); |
85b6bda9 | 461 | Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T); |
462 | Double_t projError = fProjResponseInT->GetBinError (fCoordinatesN_T); | |
463 | ||
464 | Double_t fill = responseValue / projValue ; | |
465 | if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) { | |
466 | fConditional->SetBinContent(fCoordinates2N,fill); | |
467 | // gaussian error for the moment | |
468 | Double_t err2 = TMath::Power(responseError*projValue,2) + TMath::Power(responseValue*projError,2) ; | |
469 | Double_t err = TMath::Sqrt(err2); | |
470 | err /= TMath::Power(projValue,2) ; | |
471 | fConditional->SetBinError (fCoordinates2N,err); | |
472 | } | |
c0b10ad4 | 473 | } |
474 | } | |
475 | ||
476 | //______________________________________________________________ | |
477 | ||
478 | Double_t AliCFUnfolding::GetChi2() { | |
479 | // | |
480 | // Returns the chi2 between unfolded and a priori spectrum | |
481 | // | |
482 | ||
483 | Double_t chi2 = 0. ; | |
7bc6a013 | 484 | for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) { |
85f9f9e1 | 485 | Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T); |
486 | Double_t error_unf = fUnfolded->GetBinError(fCoordinatesN_T); | |
487 | chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ; | |
c0b10ad4 | 488 | } |
489 | return chi2; | |
490 | } | |
491 | ||
492 | //______________________________________________________________ | |
493 | ||
494 | void AliCFUnfolding::SetMaxChi2PerDOF(Double_t val) { | |
495 | // | |
496 | // Max. chi2 per degree of freedom : user setting | |
497 | // | |
498 | ||
499 | Int_t nDOF = 1 ; | |
500 | for (Int_t iDim=0; iDim<fNVariables; iDim++) { | |
501 | nDOF *= fPrior->GetAxis(iDim)->GetNbins(); | |
502 | } | |
503 | AliInfo(Form("Number of degrees of freedom = %d",nDOF)); | |
504 | fMaxChi2 = val * nDOF ; | |
505 | } | |
506 | ||
507 | //______________________________________________________________ | |
508 | ||
85b6bda9 | 509 | Short_t AliCFUnfolding::Smooth() { |
c0b10ad4 | 510 | // |
511 | // Smoothes the unfolded spectrum | |
85b6bda9 | 512 | // |
513 | // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins) | |
514 | // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn | |
c0b10ad4 | 515 | // |
516 | ||
85b6bda9 | 517 | if (fSmoothFunction) { |
7bc6a013 | 518 | AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction)); |
85b6bda9 | 519 | return SmoothUsingFunction(); |
520 | } | |
7036630f | 521 | else return SmoothUsingNeighbours(fUnfolded); |
85b6bda9 | 522 | } |
523 | ||
524 | //______________________________________________________________ | |
525 | ||
7036630f | 526 | Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) { |
85b6bda9 | 527 | // |
528 | // Smoothes the unfolded spectrum using neighouring bins | |
529 | // | |
530 | ||
7036630f | 531 | Int_t const nDimensions = hist->GetNdimensions() ; |
532 | Int_t* coordinates = new Int_t[nDimensions]; | |
533 | ||
534 | Int_t* numBins = new Int_t[nDimensions]; | |
535 | for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins(); | |
c0b10ad4 | 536 | |
7036630f | 537 | //need a copy because hist will be updated during the loop, and this creates problems |
538 | THnSparse* copy = (THnSparse*)hist->Clone(); | |
c0b10ad4 | 539 | |
7bc6a013 | 540 | for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins |
7036630f | 541 | Double_t content = copy->GetBinContent(iBin,coordinates); |
7bc6a013 | 542 | Double_t error2 = TMath::Power(copy->GetBinError(iBin),2); |
c0b10ad4 | 543 | |
7036630f | 544 | printf("coord = [%d,%d]\n",coordinates[0],coordinates[1]); |
545 | ||
c0b10ad4 | 546 | // skip the under/overflow bins... |
547 | Bool_t isOutside = kFALSE ; | |
7036630f | 548 | for (Int_t iVar=0; iVar<nDimensions; iVar++) { |
549 | if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) { | |
c0b10ad4 | 550 | isOutside=kTRUE; |
551 | break; | |
552 | } | |
553 | } | |
554 | if (isOutside) continue; | |
555 | ||
556 | Int_t neighbours = 0; // number of neighbours to average with | |
557 | ||
7036630f | 558 | for (Int_t iVar=0; iVar<nDimensions; iVar++) { |
559 | if (coordinates[iVar] > 1) { // must not be on low edge border | |
560 | coordinates[iVar]-- ; //get lower neighbouring bin | |
561 | content += copy->GetBinContent(coordinates); | |
562 | error2 += TMath::Power(copy->GetBinError(coordinates),2); | |
c0b10ad4 | 563 | neighbours++; |
7036630f | 564 | coordinates[iVar]++ ; //back to initial coordinate |
c0b10ad4 | 565 | } |
7036630f | 566 | if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border |
567 | coordinates[iVar]++ ; //get upper neighbouring bin | |
568 | content += copy->GetBinContent(coordinates); | |
569 | error2 += TMath::Power(copy->GetBinError(coordinates),2); | |
c0b10ad4 | 570 | neighbours++; |
7036630f | 571 | coordinates[iVar]-- ; //back to initial coordinate |
c0b10ad4 | 572 | } |
573 | } | |
85b6bda9 | 574 | // make an average |
7036630f | 575 | hist->SetBinContent(coordinates,content/(1.+neighbours)); |
576 | hist->SetBinError (coordinates,TMath::Sqrt(error2)/(1.+neighbours)); | |
c0b10ad4 | 577 | } |
578 | delete [] numBins; | |
7036630f | 579 | delete [] coordinates ; |
c0b10ad4 | 580 | delete copy; |
85b6bda9 | 581 | return 0; |
c0b10ad4 | 582 | } |
583 | ||
85b6bda9 | 584 | //______________________________________________________________ |
585 | ||
586 | Short_t AliCFUnfolding::SmoothUsingFunction() { | |
587 | // | |
588 | // Fits the unfolded spectrum using the function fSmoothFunction | |
589 | // | |
590 | ||
7bc6a013 | 591 | AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar())); |
85b6bda9 | 592 | |
7bc6a013 | 593 | for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar))); |
594 | ||
595 | Int_t fitResult = 0; | |
85b6bda9 | 596 | |
597 | switch (fNVariables) { | |
7bc6a013 | 598 | case 1 : fitResult = fUnfolded->Projection(0) ->Fit(fSmoothFunction,fSmoothOption); break; |
599 | case 2 : fitResult = fUnfolded->Projection(1,0) ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue | |
600 | case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break; | |
601 | default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1; | |
602 | } | |
603 | ||
604 | if (fitResult != 0) { | |
605 | AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult)); | |
606 | return 1; | |
85b6bda9 | 607 | } |
608 | ||
609 | Int_t nDim = fNVariables; | |
610 | Int_t* bins = new Int_t[nDim]; // number of bins for each variable | |
611 | Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse | |
612 | ||
613 | for (Int_t iVar=0; iVar<nDim; iVar++) { | |
614 | bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins(); | |
615 | nBins *= bins[iVar]; | |
616 | } | |
617 | ||
618 | Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates) | |
619 | Double_t x[3] = {0,0,0} ; // value in bin center (max dimension is 3 (TF3)) | |
620 | ||
621 | // loop on the bins and update of fUnfolded | |
622 | // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin | |
623 | for (Long_t iBin=0; iBin<nBins; iBin++) { | |
624 | Long_t bin_tmp = iBin ; | |
625 | for (Int_t iVar=0; iVar<nDim; iVar++) { | |
626 | bin[iVar] = 1 + bin_tmp % bins[iVar] ; | |
627 | bin_tmp /= bins[iVar] ; | |
628 | x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]); | |
629 | } | |
630 | Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ; | |
12e419d5 | 631 | fUnfolded->SetBinError (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin)); |
85b6bda9 | 632 | fUnfolded->SetBinContent(bin,functionValue); |
85b6bda9 | 633 | } |
634 | return 0; | |
635 | } | |
c0b10ad4 | 636 | |
637 | //______________________________________________________________ | |
638 | ||
639 | void AliCFUnfolding::CreateFlatPrior() { | |
640 | // | |
641 | // Creates a flat prior distribution | |
642 | // | |
643 | ||
644 | AliInfo("Creating a flat a priori distribution"); | |
645 | ||
646 | // create the frame of the THnSparse given (for example) the one from the efficiency map | |
647 | fPrior = (THnSparse*) fEfficiency->Clone(); | |
648 | ||
85b6bda9 | 649 | if (fNVariables != fPrior->GetNdimensions()) |
650 | AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions())); | |
651 | ||
c0b10ad4 | 652 | Int_t nDim = fNVariables; |
653 | Int_t* bins = new Int_t[nDim]; // number of bins for each variable | |
654 | Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse | |
655 | ||
656 | for (Int_t iVar=0; iVar<nDim; iVar++) { | |
657 | bins[iVar] = fPrior->GetAxis(iVar)->GetNbins(); | |
658 | nBins *= bins[iVar]; | |
659 | } | |
660 | ||
661 | Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates) | |
662 | ||
663 | // loop that sets 1 in each bin | |
664 | for (Long_t iBin=0; iBin<nBins; iBin++) { | |
665 | Long_t bin_tmp = iBin ; | |
666 | for (Int_t iVar=0; iVar<nDim; iVar++) { | |
667 | bin[iVar] = 1 + bin_tmp % bins[iVar] ; | |
668 | bin_tmp /= bins[iVar] ; | |
669 | } | |
670 | fPrior->SetBinContent(bin,1.); // put 1 everywhere | |
85b6bda9 | 671 | fPrior->SetBinError (bin,0.); // put 0 everywhere |
c0b10ad4 | 672 | } |
673 | ||
674 | fOriginalPrior = (THnSparse*)fPrior->Clone(); | |
675 | ||
676 | delete [] bin; | |
677 | delete [] bins; | |
678 | } |