]>
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 // | |
769f5114 | 41 | // and the unfolded spectrum (OBSOLETE). // |
42 | // // | |
43 | // Chi2 calculation became absolute with the correlated error // | |
44 | // calculation. // | |
45 | // Errors on the unfolded distribution are not known until the end // | |
46 | // Use the convergence criterion instead // | |
7bc6a013 | 47 | // // |
48 | // Currently the user has to define the max. number of iterations // | |
49 | // (::SetMaxNumberOfIterations) // | |
769f5114 | 50 | // and // |
51 | // - the chi2 below which the procedure will stop // | |
52 | // (::SetMaxChi2 or ::SetMaxChi2PerDOF) (OBSOLETE) // | |
53 | // - the convergence criterion below which the procedure will stop // | |
54 | // SetMaxConvergencePerDOF(Double_t val); // | |
55 | // // | |
56 | // Correlated error calculation can be activated by using: // | |
57 | // SetUseCorrelatedErrors(Bool_t b) in combination with convergence // | |
58 | // criterion // | |
59 | // Documentation about correlated error calculation method can be // | |
60 | // found in AliCFUnfolding::CalculateCorrelatedErrors() // | |
61 | // Author: marta.verweij@cern.ch // | |
7bc6a013 | 62 | // // |
63 | // An optional possibility is to smooth the unfolded spectrum at the // | |
64 | // end of each iteration, either using a fit function // | |
65 | // (only if #dimensions <=3) // | |
66 | // or a simple averaging using the neighbouring bins values. // | |
67 | // This is possible calling the function ::UseSmoothing // | |
68 | // If no argument is passed to this function, then the second option // | |
69 | // is used. // | |
70 | // // | |
c7f196bd | 71 | // IMPORTANT: // |
72 | //----------- // | |
73 | // With this approach, the efficiency map must be calculated // | |
74 | // with *simulated* values only, otherwise the method won't work. // | |
75 | // // | |
76 | // ex: efficiency(bin_pt) = number_rec(bin_pt) / number_sim(bin_pt) // | |
77 | // // | |
78 | // the pt bin "bin_pt" must always be the same in both the efficiency // | |
79 | // numerator and denominator. // | |
80 | // This is why the efficiency map has to be created by a method // | |
81 | // from which both reconstructed and simulated values are accessible // | |
82 | // simultaneously. // | |
83 | // // | |
84 | // // | |
7bc6a013 | 85 | //---------------------------------------------------------------------// |
86 | // Author : renaud.vernet@cern.ch // | |
87 | //---------------------------------------------------------------------// | |
c0b10ad4 | 88 | |
89 | ||
90 | #include "AliCFUnfolding.h" | |
91 | #include "TMath.h" | |
92 | #include "TAxis.h" | |
93 | #include "AliLog.h" | |
85b6bda9 | 94 | #include "TF1.h" |
95 | #include "TH1D.h" | |
96 | #include "TH2D.h" | |
97 | #include "TH3D.h" | |
769f5114 | 98 | #include "TProfile.h" |
99 | #include "TRandom3.h" | |
c0b10ad4 | 100 | |
101 | ClassImp(AliCFUnfolding) | |
102 | ||
103 | //______________________________________________________________ | |
104 | ||
105 | AliCFUnfolding::AliCFUnfolding() : | |
106 | TNamed(), | |
107 | fResponse(0x0), | |
108 | fPrior(0x0), | |
c0b10ad4 | 109 | fEfficiency(0x0), |
110 | fMeasured(0x0), | |
769f5114 | 111 | fMeasuredOrig(0x0), |
112 | fMaxNumIterations(20), | |
c0b10ad4 | 113 | fNVariables(0), |
c0b10ad4 | 114 | fUseSmoothing(kFALSE), |
85b6bda9 | 115 | fSmoothFunction(0x0), |
116 | fSmoothOption(""), | |
769f5114 | 117 | fMaxConvergence(0), |
118 | fUseCorrelatedErrors(kTRUE), | |
119 | fNRandomIterations(20), | |
85b6bda9 | 120 | fOriginalPrior(0x0), |
c0b10ad4 | 121 | fInverseResponse(0x0), |
122 | fMeasuredEstimate(0x0), | |
123 | fConditional(0x0), | |
124 | fProjResponseInT(0x0), | |
125 | fUnfolded(0x0), | |
126 | fCoordinates2N(0x0), | |
127 | fCoordinatesN_M(0x0), | |
769f5114 | 128 | fCoordinatesN_T(0x0), |
129 | fRandomizedDist(0x0), | |
130 | fRandom3(0x0), | |
131 | fRandomUnfolded(0x0), | |
132 | fDeltaUnfoldedP(0x0), | |
133 | fNCalcCorrErrors(0), | |
134 | fRandomSeed(0) | |
c0b10ad4 | 135 | { |
136 | // | |
137 | // default constructor | |
138 | // | |
139 | } | |
140 | ||
141 | //______________________________________________________________ | |
142 | ||
143 | AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar, | |
144 | const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior) : | |
145 | TNamed(name,title), | |
146 | fResponse((THnSparse*)response->Clone()), | |
147 | fPrior(0x0), | |
c0b10ad4 | 148 | fEfficiency((THnSparse*)efficiency->Clone()), |
149 | fMeasured((THnSparse*)measured->Clone()), | |
769f5114 | 150 | fMeasuredOrig((THnSparse*)measured->Clone()), |
c0b10ad4 | 151 | fMaxNumIterations(0), |
152 | fNVariables(nVar), | |
c0b10ad4 | 153 | fUseSmoothing(kFALSE), |
85b6bda9 | 154 | fSmoothFunction(0x0), |
155 | fSmoothOption(""), | |
769f5114 | 156 | fMaxConvergence(0), |
157 | fUseCorrelatedErrors(kTRUE), | |
158 | fNRandomIterations(20), | |
85b6bda9 | 159 | fOriginalPrior(0x0), |
c0b10ad4 | 160 | fInverseResponse(0x0), |
161 | fMeasuredEstimate(0x0), | |
162 | fConditional(0x0), | |
163 | fProjResponseInT(0x0), | |
164 | fUnfolded(0x0), | |
165 | fCoordinates2N(0x0), | |
166 | fCoordinatesN_M(0x0), | |
769f5114 | 167 | fCoordinatesN_T(0x0), |
168 | fRandomizedDist(0x0), | |
169 | fRandom3(0x0), | |
170 | fRandomUnfolded(0x0), | |
171 | fDeltaUnfoldedP(0x0), | |
172 | fNCalcCorrErrors(0), | |
173 | fRandomSeed(0) | |
c0b10ad4 | 174 | { |
175 | // | |
176 | // named constructor | |
177 | // | |
178 | ||
179 | AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions())); | |
180 | ||
181 | if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution | |
182 | else { | |
183 | fPrior = (THnSparse*) prior->Clone(); | |
184 | fOriginalPrior = (THnSparse*)fPrior->Clone(); | |
85b6bda9 | 185 | if (fPrior->GetNdimensions() != fNVariables) |
186 | AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions())); | |
c0b10ad4 | 187 | } |
85b6bda9 | 188 | |
189 | if (fEfficiency->GetNdimensions() != fNVariables) | |
190 | AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions())); | |
191 | if (fMeasured->GetNdimensions() != fNVariables) | |
192 | AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions())); | |
193 | if (fResponse->GetNdimensions() != 2*fNVariables) | |
194 | AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions())); | |
c0b10ad4 | 195 | |
85b6bda9 | 196 | |
c0b10ad4 | 197 | for (Int_t iVar=0; iVar<fNVariables; iVar++) { |
198 | AliInfo(Form("prior matrix has %d bins in dimension %d",fPrior ->GetAxis(iVar)->GetNbins(),iVar)); | |
199 | AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar)); | |
200 | AliInfo(Form("measured matrix has %d bins in dimension %d",fMeasured ->GetAxis(iVar)->GetNbins(),iVar)); | |
201 | } | |
769f5114 | 202 | |
c0b10ad4 | 203 | Init(); |
204 | } | |
205 | ||
206 | ||
207 | //______________________________________________________________ | |
208 | ||
209 | AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) : | |
210 | TNamed(c), | |
211 | fResponse((THnSparse*)c.fResponse->Clone()), | |
212 | fPrior((THnSparse*)c.fPrior->Clone()), | |
c0b10ad4 | 213 | fEfficiency((THnSparse*)c.fEfficiency->Clone()), |
214 | fMeasured((THnSparse*)c.fMeasured->Clone()), | |
769f5114 | 215 | fMeasuredOrig((THnSparse*)c.fMeasuredOrig->Clone()), |
c0b10ad4 | 216 | fMaxNumIterations(c.fMaxNumIterations), |
217 | fNVariables(c.fNVariables), | |
c0b10ad4 | 218 | fUseSmoothing(c.fUseSmoothing), |
85b6bda9 | 219 | fSmoothFunction((TF1*)c.fSmoothFunction->Clone()), |
220 | fSmoothOption(fSmoothOption), | |
769f5114 | 221 | fMaxConvergence(c.fMaxConvergence), |
222 | fUseCorrelatedErrors(c.fUseCorrelatedErrors), | |
223 | fNRandomIterations(c.fNRandomIterations), | |
85b6bda9 | 224 | fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()), |
c0b10ad4 | 225 | fInverseResponse((THnSparse*)c.fInverseResponse->Clone()), |
226 | fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()), | |
227 | fConditional((THnSparse*)c.fConditional->Clone()), | |
228 | fProjResponseInT((THnSparse*)c.fProjResponseInT->Clone()), | |
229 | fUnfolded((THnSparse*)c.fUnfolded->Clone()), | |
230 | fCoordinates2N(new Int_t(*c.fCoordinates2N)), | |
231 | fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)), | |
769f5114 | 232 | fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T)), |
233 | fRandomizedDist((THnSparse*)c.fRandomizedDist->Clone()), | |
234 | fRandom3((TRandom3*)c.fRandom3->Clone()), | |
235 | fRandomUnfolded((THnSparse*)c.fRandomUnfolded->Clone()), | |
236 | fDeltaUnfoldedP((TProfile*)c.fDeltaUnfoldedP), | |
237 | fNCalcCorrErrors(c.fNCalcCorrErrors), | |
238 | fRandomSeed(c.fRandomSeed) | |
c0b10ad4 | 239 | { |
240 | // | |
241 | // copy constructor | |
242 | // | |
243 | } | |
244 | ||
245 | //______________________________________________________________ | |
246 | ||
247 | AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) { | |
248 | // | |
249 | // assignment operator | |
250 | // | |
251 | ||
252 | if (this!=&c) { | |
253 | TNamed::operator=(c); | |
254 | fResponse = (THnSparse*)c.fResponse->Clone() ; | |
255 | fPrior = (THnSparse*)c.fPrior->Clone() ; | |
c0b10ad4 | 256 | fEfficiency = (THnSparse*)c.fEfficiency->Clone() ; |
257 | fMeasured = (THnSparse*)c.fMeasured->Clone() ; | |
769f5114 | 258 | fMeasuredOrig = ((THnSparse*)c.fMeasuredOrig->Clone()), |
c0b10ad4 | 259 | fMaxNumIterations = c.fMaxNumIterations ; |
260 | fNVariables = c.fNVariables ; | |
769f5114 | 261 | fMaxConvergence = c.fMaxConvergence ; |
c0b10ad4 | 262 | fUseSmoothing = c.fUseSmoothing ; |
85b6bda9 | 263 | fSmoothFunction = (TF1*)c.fSmoothFunction->Clone(); |
264 | fSmoothOption = c.fSmoothOption ; | |
769f5114 | 265 | fUseCorrelatedErrors = c.fUseCorrelatedErrors; |
266 | fNRandomIterations = c.fNRandomIterations; | |
85b6bda9 | 267 | fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ; |
c0b10ad4 | 268 | fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ; |
269 | fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ; | |
270 | fConditional = (THnSparse*)c.fConditional->Clone() ; | |
271 | fProjResponseInT = (THnSparse*)c.fProjResponseInT->Clone() ; | |
272 | fUnfolded = (THnSparse*)c.fUnfolded->Clone() ; | |
273 | fCoordinates2N = new Int_t(*c.fCoordinates2N) ; | |
274 | fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ; | |
275 | fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ; | |
769f5114 | 276 | fRandomizedDist = (THnSparse*)c.fRandomizedDist->Clone(); |
277 | fRandom3 = (TRandom3*)c.fRandom3->Clone(); | |
278 | fRandomUnfolded = (THnSparse*)c.fRandomUnfolded->Clone(); | |
279 | fDeltaUnfoldedP = (TProfile*)c.fDeltaUnfoldedP; | |
280 | fNCalcCorrErrors = c.fNCalcCorrErrors; | |
281 | fRandomSeed = c.fRandomSeed ; | |
c0b10ad4 | 282 | } |
283 | return *this; | |
284 | } | |
285 | ||
286 | //______________________________________________________________ | |
287 | ||
288 | AliCFUnfolding::~AliCFUnfolding() { | |
289 | // | |
290 | // destructor | |
291 | // | |
292 | if (fResponse) delete fResponse; | |
293 | if (fPrior) delete fPrior; | |
c0b10ad4 | 294 | if (fEfficiency) delete fEfficiency; |
295 | if (fMeasured) delete fMeasured; | |
769f5114 | 296 | if (fMeasuredOrig) delete fMeasuredOrig; |
85b6bda9 | 297 | if (fSmoothFunction) delete fSmoothFunction; |
298 | if (fOriginalPrior) delete fOriginalPrior; | |
c0b10ad4 | 299 | if (fInverseResponse) delete fInverseResponse; |
300 | if (fMeasuredEstimate) delete fMeasuredEstimate; | |
301 | if (fConditional) delete fConditional; | |
302 | if (fProjResponseInT) delete fProjResponseInT; | |
303 | if (fCoordinates2N) delete [] fCoordinates2N; | |
304 | if (fCoordinatesN_M) delete [] fCoordinatesN_M; | |
305 | if (fCoordinatesN_T) delete [] fCoordinatesN_T; | |
769f5114 | 306 | if (fRandomizedDist) delete fRandomizedDist; |
307 | if (fRandom3) delete fRandom3; | |
308 | if (fRandomUnfolded) delete fRandomUnfolded; | |
309 | if (fDeltaUnfoldedP) delete fDeltaUnfoldedP; | |
c0b10ad4 | 310 | } |
311 | ||
312 | //______________________________________________________________ | |
313 | ||
314 | void AliCFUnfolding::Init() { | |
315 | // | |
316 | // initialisation function : creates internal settings | |
317 | // | |
318 | ||
769f5114 | 319 | fRandom3 = new TRandom3(fRandomSeed); |
320 | ||
c0b10ad4 | 321 | fCoordinates2N = new Int_t[2*fNVariables]; |
322 | fCoordinatesN_M = new Int_t[fNVariables]; | |
323 | fCoordinatesN_T = new Int_t[fNVariables]; | |
324 | ||
325 | // create the matrix of conditional probabilities P(M|T) | |
326 | CreateConditional(); | |
327 | ||
328 | // create the frame of the inverse response matrix | |
329 | fInverseResponse = (THnSparse*) fResponse->Clone(); | |
330 | // create the frame of the unfolded spectrum | |
331 | fUnfolded = (THnSparse*) fPrior->Clone(); | |
769f5114 | 332 | // create the frame of the random unfolded spectrum |
333 | fRandomUnfolded = (THnSparse*) fPrior->Clone(); | |
c0b10ad4 | 334 | // create the frame of the measurement estimate spectrum |
335 | fMeasuredEstimate = (THnSparse*) fMeasured->Clone(); | |
769f5114 | 336 | // create the frame of the original measurement spectrum |
337 | fMeasuredOrig = (THnSparse*) fMeasured->Clone(); | |
338 | ||
339 | InitDeltaUnfoldedProfile(); | |
340 | ||
c0b10ad4 | 341 | } |
342 | ||
343 | //______________________________________________________________ | |
344 | ||
345 | void AliCFUnfolding::CreateEstMeasured() { | |
346 | // | |
347 | // This function creates a estimate (M) of the reconstructed spectrum | |
7bc6a013 | 348 | // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND) |
c0b10ad4 | 349 | // |
350 | // --> P(M) = SUM { P(M|T) * P(T) } | |
7bc6a013 | 351 | // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)} |
c0b10ad4 | 352 | // |
353 | // This is needed to calculate the inverse response matrix | |
354 | // | |
355 | ||
356 | ||
357 | // clean the measured estimate spectrum | |
85f9f9e1 | 358 | fMeasuredEstimate->Reset(); |
359 | ||
85b6bda9 | 360 | THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone(); |
361 | priorTimesEff->Multiply(fEfficiency); | |
362 | ||
c0b10ad4 | 363 | // fill it |
7bc6a013 | 364 | for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) { |
c0b10ad4 | 365 | Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N); |
366 | GetCoordinates(); | |
85b6bda9 | 367 | Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T); |
85b6bda9 | 368 | Double_t fill = conditionalValue * priorTimesEffValue ; |
369 | ||
370 | if (fill>0.) { | |
12e419d5 | 371 | fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill); |
769f5114 | 372 | fMeasuredEstimate->SetBinError(fCoordinatesN_M,0.); |
85b6bda9 | 373 | } |
c0b10ad4 | 374 | } |
85b6bda9 | 375 | delete priorTimesEff ; |
c0b10ad4 | 376 | } |
377 | ||
378 | //______________________________________________________________ | |
379 | ||
380 | void AliCFUnfolding::CreateInvResponse() { | |
381 | // | |
382 | // Creates the inverse response matrix (INV) with Bayesian method | |
383 | // : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E) | |
384 | // | |
385 | // --> P(T|M) = P(M|T) * P(T) * eff(T) / SUM { P(M|T) * P(T) } | |
386 | // --> INV(i,j) = COND(i,j) * T(j) * E(j) / SUM_k { COND(i,k) * T(k) } | |
387 | // | |
388 | ||
85b6bda9 | 389 | THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone(); |
390 | priorTimesEff->Multiply(fEfficiency); | |
391 | ||
7bc6a013 | 392 | for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) { |
c0b10ad4 | 393 | Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N); |
394 | GetCoordinates(); | |
85b6bda9 | 395 | Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M); |
85b6bda9 | 396 | Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T); |
85b6bda9 | 397 | Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ; |
85b6bda9 | 398 | if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) { |
399 | fInverseResponse->SetBinContent(fCoordinates2N,fill); | |
769f5114 | 400 | fInverseResponse->SetBinError (fCoordinates2N,0.); |
85b6bda9 | 401 | } |
402 | } | |
403 | delete priorTimesEff ; | |
c0b10ad4 | 404 | } |
405 | ||
406 | //______________________________________________________________ | |
407 | ||
408 | void AliCFUnfolding::Unfold() { | |
409 | // | |
410 | // Main routine called by the user : | |
769f5114 | 411 | // it calculates the unfolded spectrum from the response matrix, measured spectrum and efficiency |
412 | // several iterations are performed until a reasonable chi2 or convergence criterion is reached | |
c0b10ad4 | 413 | // |
414 | ||
769f5114 | 415 | Int_t iIterBayes = 0 ; |
416 | Double_t convergence = 0.; | |
c0b10ad4 | 417 | |
418 | for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations | |
419 | CreateEstMeasured(); | |
420 | CreateInvResponse(); | |
421 | CreateUnfolded(); | |
769f5114 | 422 | if (fUseCorrelatedErrors) { |
423 | convergence = GetConvergence(); | |
424 | AliDebug(0,Form("convergence at iteration %d is %e",iIterBayes,convergence)); | |
425 | } | |
426 | else AliWarning("No errors will be calculated. Activate SetUseCorrelatedErrors(kTRUE)\n"); | |
427 | ||
428 | if (fMaxConvergence>0. && convergence<fMaxConvergence && fNCalcCorrErrors<1) { | |
429 | fNRandomIterations=iIterBayes; | |
c0b10ad4 | 430 | break; |
431 | } | |
769f5114 | 432 | |
c0b10ad4 | 433 | // update the prior distribution |
85b6bda9 | 434 | if (fUseSmoothing) { |
435 | if (Smooth()) { | |
436 | AliError("Couldn't smooth the unfolded spectrum!!"); | |
769f5114 | 437 | if (fUseCorrelatedErrors) { |
438 | if (fNCalcCorrErrors>0) { | |
439 | AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence)); | |
440 | } | |
441 | else { | |
442 | AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence)); | |
443 | } | |
444 | } | |
85b6bda9 | 445 | return; |
446 | } | |
447 | } | |
769f5114 | 448 | if(fNCalcCorrErrors>0) { |
449 | if (fPrior) delete fPrior ; | |
450 | fPrior = (THnSparse*)fRandomUnfolded->Clone() ; | |
451 | } | |
452 | else { | |
453 | if (fPrior) delete fPrior ; | |
454 | fPrior = (THnSparse*)fUnfolded->Clone() ; | |
455 | } | |
456 | } | |
457 | ||
458 | if (fUseCorrelatedErrors && fNCalcCorrErrors==0) { | |
459 | fNCalcCorrErrors=1; | |
460 | CalculateCorrelatedErrors(); | |
461 | } | |
462 | ||
463 | if (fUseCorrelatedErrors) { | |
464 | if (fNCalcCorrErrors>1) { | |
465 | AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence)); | |
466 | } | |
467 | else if(fNCalcCorrErrors>0) { | |
468 | AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence)); | |
469 | } | |
c0b10ad4 | 470 | } |
c0b10ad4 | 471 | } |
472 | ||
473 | //______________________________________________________________ | |
474 | ||
475 | void AliCFUnfolding::CreateUnfolded() { | |
476 | // | |
477 | // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV) | |
478 | // We have P(T) = SUM { P(T|M) * P(M) } | |
479 | // --> T(i) = SUM_k { INV(i,k) * M(k) } | |
480 | // | |
481 | ||
482 | ||
483 | // clear the unfolded spectrum | |
769f5114 | 484 | if(fNCalcCorrErrors>0) { |
485 | //unfold randomized dist | |
486 | fRandomUnfolded->Reset(); | |
487 | } | |
488 | else { | |
489 | //unfold measured dist | |
490 | fUnfolded->Reset(); | |
491 | } | |
c0b10ad4 | 492 | |
7bc6a013 | 493 | for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) { |
c0b10ad4 | 494 | Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N); |
495 | GetCoordinates(); | |
85b6bda9 | 496 | Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T); |
85b6bda9 | 497 | Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M); |
85b6bda9 | 498 | Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ; |
769f5114 | 499 | |
85b6bda9 | 500 | if (fill>0.) { |
769f5114 | 501 | |
502 | Double_t err = 0.; | |
503 | ||
504 | if(fNCalcCorrErrors>0) { | |
505 | fRandomUnfolded->SetBinError(fCoordinatesN_T,err); | |
506 | fRandomUnfolded->AddBinContent(fCoordinatesN_T,fill); | |
507 | } | |
508 | else { | |
509 | fUnfolded->SetBinError(fCoordinatesN_T,err); | |
510 | fUnfolded->AddBinContent(fCoordinatesN_T,fill); | |
511 | } | |
85b6bda9 | 512 | } |
c0b10ad4 | 513 | } |
514 | } | |
769f5114 | 515 | |
516 | //______________________________________________________________ | |
517 | ||
518 | void AliCFUnfolding::CalculateCorrelatedErrors() { | |
519 | ||
520 | fRandomizedDist = (THnSparse*) fMeasuredOrig->Clone(); | |
521 | fPrior = (THnSparse*) fOriginalPrior->Clone(); | |
522 | ||
523 | // Step 1: Create randomized distribution (fRandomizedDist) of each bin of the measured spectrum to calculate correlated errors. Poisson statistics: mean = measured value of bin | |
524 | // Step 2: Unfold randomized distribution | |
525 | // Step 3: Store difference of unfolded spectrum from measured distribution and unfolded distribution from randomized distribution -> fDeltaUnfoldedP (TProfile with option "S") | |
526 | // Step 4: Repeat Step 1-3 several times (fNRandomIterations) | |
527 | // Step 5: The spread of fDeltaUnfoldedP for each bin is the error on the unfolded spectrum of that specific bin | |
528 | ||
529 | //Do fNRandomIterations = bayes iterations performed | |
530 | for(int i=0; i<fNRandomIterations; i++) { | |
531 | if (fPrior) delete fPrior ; | |
532 | if (fRandomizedDist) delete fRandomizedDist ; | |
533 | fPrior = (THnSparse*) fOriginalPrior->Clone(); | |
534 | fRandomizedDist = (THnSparse*) fMeasuredOrig->Clone(); | |
535 | CreateRandomizedDist(); | |
536 | if (fMeasured) delete fMeasured ; | |
537 | fMeasured = (THnSparse*) fRandomizedDist->Clone(); | |
538 | //Unfold fRandomizedDist | |
539 | Unfold(); | |
540 | FillDeltaUnfoldedProfile(); | |
541 | } | |
542 | ||
543 | // Get statistical errors for final unfolded spectrum | |
544 | // ie. spread of each pt bin in fDeltaUnfoldedP | |
545 | Double_t sigma = 0.; | |
546 | Double_t dummy = 0.; | |
547 | for (Long_t iBin=0; iBin<fRandomUnfolded->GetNbins(); iBin++) { | |
548 | dummy = fUnfolded->GetBinContent(iBin,fCoordinatesN_M); | |
549 | sigma = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M[0]); | |
550 | fUnfolded->SetBinError(fCoordinatesN_M,sigma); | |
551 | fNCalcCorrErrors = 2; | |
552 | } | |
553 | ||
554 | } | |
555 | ||
556 | //______________________________________________________________ | |
557 | void AliCFUnfolding::InitDeltaUnfoldedProfile() { | |
558 | // | |
559 | //Initialize the fDeltaUnfoldedP profile | |
560 | //Errors will be filled with spread between unfolded measured and unfolded randomized spectra | |
561 | // | |
562 | ||
563 | Int_t nbinsx = fResponse->GetAxis(0)->GetNbins(); | |
564 | Double_t xbins[nbinsx]; | |
565 | for(int ix=0; ix<nbinsx; ix++) { | |
566 | xbins[ix] = fResponse->GetAxis(0)->GetBinLowEdge(ix+1); | |
567 | } | |
568 | xbins[nbinsx] = fResponse->GetAxis(0)->GetBinUpEdge(nbinsx); | |
569 | fDeltaUnfoldedP = new TProfile("fDeltaUnfoldedP","Profile of pTUnfolded with spread in error",nbinsx,xbins,"S"); | |
570 | } | |
571 | //______________________________________________________________ | |
572 | void AliCFUnfolding::CreateRandomizedDist() { | |
573 | // | |
574 | // Create randomized dist from measured distribution | |
575 | // | |
576 | ||
577 | Double_t random = 0.; | |
578 | Double_t measuredValue = 0.; | |
579 | Double_t measuredError = 0.; | |
580 | for (Long_t iBin=0; iBin<fRandomizedDist->GetNbins(); iBin++) { | |
581 | measuredValue = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean | |
582 | measuredError = fMeasuredOrig->GetBinError(fCoordinatesN_M); //used as sigma | |
583 | // random = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10) | |
584 | random = fRandom3->Gaus(measuredValue,measuredError); | |
585 | fRandomizedDist->SetBinContent(iBin,random); | |
586 | } | |
587 | } | |
588 | ||
589 | //______________________________________________________________ | |
590 | void AliCFUnfolding::FillDeltaUnfoldedProfile() { | |
591 | // | |
592 | // Store difference of unfolded spectrum from measured distribution and unfolded distribution from randomized distribution | |
593 | // | |
594 | ||
595 | for (Long_t iBin2=0; iBin2<fRandomUnfolded->GetNbins(); iBin2++) { | |
596 | Double_t delta = fUnfolded->GetBinContent(iBin2,fCoordinatesN_M)-fRandomUnfolded->GetBinContent(iBin2,fCoordinatesN_M); | |
597 | fDeltaUnfoldedP->Fill(fDeltaUnfoldedP->GetBinCenter(fCoordinatesN_M[0]),delta); | |
598 | } | |
599 | } | |
600 | ||
85b6bda9 | 601 | //______________________________________________________________ |
602 | ||
c0b10ad4 | 603 | void AliCFUnfolding::GetCoordinates() { |
604 | // | |
605 | // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N) | |
606 | // | |
607 | for (Int_t i = 0; i<fNVariables ; i++) { | |
608 | fCoordinatesN_M[i] = fCoordinates2N[i]; | |
609 | fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables]; | |
610 | } | |
611 | } | |
612 | ||
613 | //______________________________________________________________ | |
614 | ||
615 | void AliCFUnfolding::CreateConditional() { | |
616 | // | |
617 | // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R | |
618 | // | |
619 | // --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) } | |
620 | // | |
621 | ||
622 | fConditional = (THnSparse*) fResponse->Clone(); // output of this function | |
623 | fProjResponseInT = (THnSparse*) fPrior->Clone(); // output denominator : | |
624 | // projection of the response matrix on the TRUE axis | |
85b6bda9 | 625 | Int_t* dim = new Int_t [fNVariables]; |
626 | for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1) | |
627 | fProjResponseInT = fConditional->Projection(fNVariables,dim,"E"); //project | |
628 | delete [] dim; | |
769f5114 | 629 | |
c0b10ad4 | 630 | // fill the conditional probability matrix |
7bc6a013 | 631 | for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) { |
c0b10ad4 | 632 | Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N); |
633 | GetCoordinates(); | |
85b6bda9 | 634 | Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T); |
769f5114 | 635 | |
85b6bda9 | 636 | Double_t fill = responseValue / projValue ; |
637 | if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) { | |
638 | fConditional->SetBinContent(fCoordinates2N,fill); | |
769f5114 | 639 | Double_t err = 0.; |
85b6bda9 | 640 | fConditional->SetBinError (fCoordinates2N,err); |
641 | } | |
c0b10ad4 | 642 | } |
643 | } | |
769f5114 | 644 | //______________________________________________________________ |
645 | ||
646 | Int_t AliCFUnfolding::GetDOF() { | |
647 | // | |
648 | // number of dof = number of bins | |
649 | // | |
650 | ||
651 | Int_t nDOF = 1 ; | |
652 | for (Int_t iDim=0; iDim<fNVariables; iDim++) { | |
653 | nDOF *= fPrior->GetAxis(iDim)->GetNbins(); | |
654 | } | |
655 | AliDebug(0,Form("Number of degrees of freedom = %d",nDOF)); | |
656 | return nDOF; | |
657 | } | |
c0b10ad4 | 658 | |
659 | //______________________________________________________________ | |
660 | ||
661 | Double_t AliCFUnfolding::GetChi2() { | |
662 | // | |
663 | // Returns the chi2 between unfolded and a priori spectrum | |
769f5114 | 664 | // This function became absolute with the correlated error calculation. |
665 | // Errors on the unfolded distribution are not known until the end | |
666 | // Use the convergence criterion instead | |
c0b10ad4 | 667 | // |
668 | ||
769f5114 | 669 | Double_t chi2 = 0. ; |
670 | Double_t error_unf = 0.; | |
7bc6a013 | 671 | for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) { |
85f9f9e1 | 672 | Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T); |
769f5114 | 673 | error_unf = fUnfolded->GetBinError(fCoordinatesN_T); |
85f9f9e1 | 674 | chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ; |
c0b10ad4 | 675 | } |
676 | return chi2; | |
677 | } | |
678 | ||
679 | //______________________________________________________________ | |
680 | ||
769f5114 | 681 | Double_t AliCFUnfolding::GetConvergence() { |
682 | // | |
683 | // Returns convergence criterion = \sum_t ((U_t^{n-1}-U_t^n)/U_t^{n-1})^2 | |
684 | // U is unfolded spectrum, t is the bin, n = current, n-1 = previous | |
685 | // | |
686 | Double_t convergence = 0.; | |
687 | Double_t priorValue = 0.; | |
688 | Double_t currentValue = 0.; | |
689 | for (Long_t iBin=0; iBin < fPrior->GetNbins(); iBin++) { | |
690 | priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T); | |
691 | if (fNCalcCorrErrors > 0) | |
692 | currentValue = fRandomUnfolded->GetBinContent(fCoordinatesN_T); | |
693 | else | |
694 | currentValue = fUnfolded->GetBinContent(fCoordinatesN_T); | |
695 | ||
696 | if (priorValue > 0.) | |
697 | convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue); | |
698 | else { | |
699 | AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue)); | |
700 | convergence += 0.; | |
701 | } | |
702 | } | |
703 | return convergence; | |
704 | } | |
705 | ||
706 | //______________________________________________________________ | |
707 | ||
708 | void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) { | |
c0b10ad4 | 709 | // |
769f5114 | 710 | // Max. convergence criterion per degree of freedom : user setting |
711 | // convergence criterion = DOF*val; DOF = number of bins | |
712 | // In Jan-Fiete's multiplicity note: Convergence criterion = DOF*0.001^2 | |
c0b10ad4 | 713 | // |
714 | ||
769f5114 | 715 | fUseCorrelatedErrors = kTRUE; |
716 | Int_t nDOF = GetDOF() ; | |
717 | fMaxConvergence = val * nDOF ; | |
718 | AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF)); | |
c0b10ad4 | 719 | } |
720 | ||
721 | //______________________________________________________________ | |
722 | ||
85b6bda9 | 723 | Short_t AliCFUnfolding::Smooth() { |
c0b10ad4 | 724 | // |
725 | // Smoothes the unfolded spectrum | |
85b6bda9 | 726 | // |
727 | // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins) | |
728 | // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn | |
c0b10ad4 | 729 | // |
730 | ||
85b6bda9 | 731 | if (fSmoothFunction) { |
7bc6a013 | 732 | AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction)); |
85b6bda9 | 733 | return SmoothUsingFunction(); |
734 | } | |
7036630f | 735 | else return SmoothUsingNeighbours(fUnfolded); |
85b6bda9 | 736 | } |
737 | ||
738 | //______________________________________________________________ | |
739 | ||
7036630f | 740 | Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) { |
85b6bda9 | 741 | // |
742 | // Smoothes the unfolded spectrum using neighouring bins | |
743 | // | |
744 | ||
7036630f | 745 | Int_t const nDimensions = hist->GetNdimensions() ; |
746 | Int_t* coordinates = new Int_t[nDimensions]; | |
747 | ||
748 | Int_t* numBins = new Int_t[nDimensions]; | |
749 | for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins(); | |
c0b10ad4 | 750 | |
7036630f | 751 | //need a copy because hist will be updated during the loop, and this creates problems |
752 | THnSparse* copy = (THnSparse*)hist->Clone(); | |
c0b10ad4 | 753 | |
7bc6a013 | 754 | for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins |
7036630f | 755 | Double_t content = copy->GetBinContent(iBin,coordinates); |
7bc6a013 | 756 | Double_t error2 = TMath::Power(copy->GetBinError(iBin),2); |
c0b10ad4 | 757 | |
758 | // skip the under/overflow bins... | |
759 | Bool_t isOutside = kFALSE ; | |
7036630f | 760 | for (Int_t iVar=0; iVar<nDimensions; iVar++) { |
761 | if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) { | |
c0b10ad4 | 762 | isOutside=kTRUE; |
763 | break; | |
764 | } | |
765 | } | |
766 | if (isOutside) continue; | |
767 | ||
768 | Int_t neighbours = 0; // number of neighbours to average with | |
769 | ||
7036630f | 770 | for (Int_t iVar=0; iVar<nDimensions; iVar++) { |
771 | if (coordinates[iVar] > 1) { // must not be on low edge border | |
772 | coordinates[iVar]-- ; //get lower neighbouring bin | |
773 | content += copy->GetBinContent(coordinates); | |
774 | error2 += TMath::Power(copy->GetBinError(coordinates),2); | |
c0b10ad4 | 775 | neighbours++; |
7036630f | 776 | coordinates[iVar]++ ; //back to initial coordinate |
c0b10ad4 | 777 | } |
7036630f | 778 | if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border |
779 | coordinates[iVar]++ ; //get upper neighbouring bin | |
780 | content += copy->GetBinContent(coordinates); | |
781 | error2 += TMath::Power(copy->GetBinError(coordinates),2); | |
c0b10ad4 | 782 | neighbours++; |
7036630f | 783 | coordinates[iVar]-- ; //back to initial coordinate |
c0b10ad4 | 784 | } |
785 | } | |
85b6bda9 | 786 | // make an average |
7036630f | 787 | hist->SetBinContent(coordinates,content/(1.+neighbours)); |
788 | hist->SetBinError (coordinates,TMath::Sqrt(error2)/(1.+neighbours)); | |
c0b10ad4 | 789 | } |
790 | delete [] numBins; | |
7036630f | 791 | delete [] coordinates ; |
c0b10ad4 | 792 | delete copy; |
85b6bda9 | 793 | return 0; |
c0b10ad4 | 794 | } |
795 | ||
85b6bda9 | 796 | //______________________________________________________________ |
797 | ||
798 | Short_t AliCFUnfolding::SmoothUsingFunction() { | |
799 | // | |
800 | // Fits the unfolded spectrum using the function fSmoothFunction | |
801 | // | |
802 | ||
7bc6a013 | 803 | AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar())); |
85b6bda9 | 804 | |
7bc6a013 | 805 | for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar))); |
806 | ||
807 | Int_t fitResult = 0; | |
85b6bda9 | 808 | |
809 | switch (fNVariables) { | |
7bc6a013 | 810 | case 1 : fitResult = fUnfolded->Projection(0) ->Fit(fSmoothFunction,fSmoothOption); break; |
811 | case 2 : fitResult = fUnfolded->Projection(1,0) ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue | |
812 | case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break; | |
813 | default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1; | |
814 | } | |
815 | ||
816 | if (fitResult != 0) { | |
817 | AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult)); | |
818 | return 1; | |
85b6bda9 | 819 | } |
820 | ||
821 | Int_t nDim = fNVariables; | |
822 | Int_t* bins = new Int_t[nDim]; // number of bins for each variable | |
823 | Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse | |
824 | ||
825 | for (Int_t iVar=0; iVar<nDim; iVar++) { | |
826 | bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins(); | |
827 | nBins *= bins[iVar]; | |
828 | } | |
829 | ||
830 | Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates) | |
831 | Double_t x[3] = {0,0,0} ; // value in bin center (max dimension is 3 (TF3)) | |
832 | ||
833 | // loop on the bins and update of fUnfolded | |
834 | // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin | |
835 | for (Long_t iBin=0; iBin<nBins; iBin++) { | |
836 | Long_t bin_tmp = iBin ; | |
837 | for (Int_t iVar=0; iVar<nDim; iVar++) { | |
838 | bin[iVar] = 1 + bin_tmp % bins[iVar] ; | |
839 | bin_tmp /= bins[iVar] ; | |
840 | x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]); | |
841 | } | |
842 | Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ; | |
12e419d5 | 843 | fUnfolded->SetBinError (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin)); |
85b6bda9 | 844 | fUnfolded->SetBinContent(bin,functionValue); |
85b6bda9 | 845 | } |
846 | return 0; | |
847 | } | |
c0b10ad4 | 848 | |
849 | //______________________________________________________________ | |
850 | ||
851 | void AliCFUnfolding::CreateFlatPrior() { | |
852 | // | |
853 | // Creates a flat prior distribution | |
854 | // | |
855 | ||
856 | AliInfo("Creating a flat a priori distribution"); | |
857 | ||
858 | // create the frame of the THnSparse given (for example) the one from the efficiency map | |
859 | fPrior = (THnSparse*) fEfficiency->Clone(); | |
860 | ||
85b6bda9 | 861 | if (fNVariables != fPrior->GetNdimensions()) |
862 | AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions())); | |
863 | ||
c0b10ad4 | 864 | Int_t nDim = fNVariables; |
865 | Int_t* bins = new Int_t[nDim]; // number of bins for each variable | |
866 | Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse | |
867 | ||
868 | for (Int_t iVar=0; iVar<nDim; iVar++) { | |
869 | bins[iVar] = fPrior->GetAxis(iVar)->GetNbins(); | |
870 | nBins *= bins[iVar]; | |
871 | } | |
872 | ||
873 | Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates) | |
874 | ||
875 | // loop that sets 1 in each bin | |
876 | for (Long_t iBin=0; iBin<nBins; iBin++) { | |
877 | Long_t bin_tmp = iBin ; | |
878 | for (Int_t iVar=0; iVar<nDim; iVar++) { | |
879 | bin[iVar] = 1 + bin_tmp % bins[iVar] ; | |
880 | bin_tmp /= bins[iVar] ; | |
881 | } | |
882 | fPrior->SetBinContent(bin,1.); // put 1 everywhere | |
85b6bda9 | 883 | fPrior->SetBinError (bin,0.); // put 0 everywhere |
c0b10ad4 | 884 | } |
885 | ||
886 | fOriginalPrior = (THnSparse*)fPrior->Clone(); | |
887 | ||
888 | delete [] bin; | |
889 | delete [] bins; | |
890 | } |