source: trunk/MagicSoft/Mars/mtemp/mmpi/MUnfold.cc@ 8782

Last change on this file since 8782 was 5554, checked in by rwagner, 21 years ago
*** empty log message ***
File size: 93.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17! Author(s) : T. Bretz 02/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
18! W. Wittek 09/2002 <mailto:wittek@mppmu.mpg.de>
19! R. Wagner, 11/2004 <mailto:rwagner@mppmu.mpg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2005
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MUnfold
29//
30/////////////////////////////////////////////////////////////////////////////
31#include "MUnfold.h"
32
33#include <TMath.h>
34#include <TRandom3.h>
35#include <TVector.h>
36#include <TMatrixD.h>
37#include <TMatrix.h>
38#include <TH1.h>
39#include <TH2.h>
40#include <TH3.h>
41#include <TProfile.h>
42#include <TF1.h>
43#include <TMinuit.h>
44#include <TCanvas.h>
45#include <TMarker.h>
46
47#include "MLog.h"
48#include "MLogManip.h"
49
50ClassImp(MUnfold);
51using namespace std;
52
53// -----------------------------------------------------------------------
54//
55// fcnSmooth (used by SmoothMigrationMatrix)
56//
57// is called by MINUIT
58// for given values of the parameters it calculates the function
59// to be minimized
60//
61static void fcnSmooth(Int_t &npar, Double_t *gin, Double_t &f,
62 Double_t *par, Int_t iflag)
63{
64 MUnfold &gUnfold = *(MUnfold*)gMinuit->GetObjectFit();
65
66 Double_t a0 = par[0];
67 Double_t a1 = par[1];
68 Double_t a2 = par[2];
69
70 Double_t b0 = par[3];
71 Double_t b1 = par[4];
72 Double_t b2 = par[5];
73
74 // loop over bins of log10(E-true)
75 Double_t chi2 = 0.0;
76 Int_t npoints = 0;
77 Double_t func[20];
78
79 for (UInt_t j=0; j<gUnfold.fNb; j++)
80 {
81 Double_t yj = ((double)j) + 0.5;
82 Double_t mean = a0 + a1*yj + a2*yj*yj + yj;
83 Double_t RMS = b0 + b1*yj + b2*yj*yj;
84
85 if (RMS <= 0.0)
86 {
87 chi2 = 1.e20;
88 break;
89 }
90
91 // loop over bins of log10(E-est)
92
93 //.......................................
94 Double_t function;
95 Double_t sum=0.0;
96 for (UInt_t i=0; i<gUnfold.fNa; i++)
97 {
98 Double_t xlow = (double)i;
99 Double_t xup = xlow + 1.0;
100 Double_t xl = (xlow- mean) / RMS;
101 Double_t xu = (xup - mean) / RMS;
102 function = (TMath::Freq(xu) - TMath::Freq(xl));
103
104 //*fLog << "i, xl, xu, function = " << i << ", " << xl << ", "
105 // << xu << ", " << function << endl;
106
107 if (function < 1.e-10)
108 function = 0.0;
109
110 func[i] = function;
111 sum += function;
112 }
113
114 // *fLog << "mean, RMS = " << mean << ", " << RMS
115 // << ", j , sum of function = " << j << ", " << sum << endl;
116
117 //.......................................
118
119 for (UInt_t i=0; i<gUnfold.fNa; i++)
120 {
121 if (sum != 0.0)
122 func[i] /= sum;
123
124 gUnfold.fMigSmoo(i,j) = func[i];
125 gUnfold.fMigChi2(i,j) = 0.0;
126
127 // if relative error is greater than 30 % ignore the point
128
129 if (gUnfold.fMigOrig(i,j) != 0 &&
130 gUnfold.fMigOrigerr2(i,j) != 0 &&
131 func[i] != 0 )
132 {
133 if (gUnfold.fMigOrigerr2(i,j)/
134 (gUnfold.fMigOrig(i,j)*gUnfold.fMigOrig(i,j)) <= 0.09)
135 {
136 gUnfold.fMigChi2(i,j) = ( gUnfold.fMigOrig(i,j) - func[i] )
137 * ( gUnfold.fMigOrig(i,j) - func[i] )
138 / gUnfold.fMigOrigerr2(i,j);
139 chi2 += gUnfold.fMigChi2(i,j);
140 npoints += 1;
141 }
142 }
143 }
144 //.......................................
145
146 }
147 f = chi2;
148
149 //*fLog << "fcnSmooth : f = " << f << endl;
150
151 //--------------------------------------------------------------------
152 // final calculations
153 if (iflag == 3)
154 {
155 Int_t NDF = npoints - npar;
156 Double_t prob = TMath::Prob(chi2, NDF);
157
158 gLog << "fcnSmooth : npoints, chi2, NDF, prob = " << npoints << ", ";
159 gLog << chi2 << ", " << NDF << ", " << prob << endl;
160 gLog << "=======================================" << endl;
161 }
162}
163
164
165// -----------------------------------------------------------------------
166//
167// fcnTikhonov2 (used by Tikhonov2)
168//
169// is called by MINUIT
170// for given values of the parameters it calculates the function F
171// the free parameters are the first (fNb-1) elements
172// of the normalized unfolded distribution
173//
174static void fcnTikhonov2(Int_t &npar, Double_t *gin, Double_t &f,
175 Double_t *par, Int_t iflag)
176{
177 MUnfold &gUnfold = *(MUnfold*)gMinuit->GetObjectFit();
178
179 // (npar+1) is the number of bins of the unfolded distribuition (fNb)
180 // npar is the number of free parameters (fNb-1)
181
182 UInt_t npar1 = npar + 1;
183
184 UInt_t fNa = gUnfold.fNa;
185 UInt_t fNb = gUnfold.fNb;
186 if (npar1 != fNb)
187 {
188 gLog << "fcnTikhonov2 : inconsistency in number of parameters; npar, fNb = ";
189 gLog << npar << ", " << fNb << endl;
190 //return;
191 }
192 npar1 = fNb;
193
194 TMatrixD p(npar1, 1);
195 TMatrixD &fVb = gUnfold.fVb;
196
197 // p is the normalized unfolded distribution
198 // sum(p(i,0)) from i=0 to npar is equal to 1
199 Double_t sum = 0.0;
200 for (Int_t i=0; i<npar; i++)
201 {
202 p(i,0) = par[i];
203 sum += par[i];
204 }
205 p(npar,0) = 1.0 - sum;
206
207
208 // all p(i,0) have to be greater than zero
209 for (UInt_t i=0; i<npar1; i++)
210 if (p(i,0) <= 0.0)
211 {
212 f = 1.e20;
213 return;
214 }
215
216 //.......................
217 // take least squares result for the normaliztion
218 TMatrixD alpha(gUnfold.fMigrat, TMatrixD::kMult, p);
219
220 //TMatrixD v4 (gUnfold.fVa, TMatrixD::kTransposeMult,
221 // gUnfold.fVacovInv);
222 //TMatrixD norma(v4, TMatrixD::kMult, gUnfold.fVa);
223
224 TMatrixD v5 (alpha, TMatrixD::kTransposeMult, gUnfold.fVacovInv);
225 TMatrixD normb(v5, TMatrixD::kMult, alpha);
226
227 TMatrixD normc(v5, TMatrixD::kMult, gUnfold.fVa);
228
229 Double_t norm = normc(0,0)/normb(0,0);
230
231 //.......................
232
233 // b is the unnormalized unfolded distribution
234 // sum(b(i,0)) from i=0 to npar is equal to norm
235 // (the total number of events)
236 for (UInt_t i=0; i<npar1; i++)
237 fVb(i,0) = p(i,0) * norm;
238
239 TMatrixD Gb(gUnfold.fMigrat, TMatrixD::kMult, fVb);
240 TMatrixD v3(fNa, 1);
241 v3 = gUnfold.fVa;
242 v3 -= Gb;
243
244 TMatrixD v1(1,fNa);
245 for (UInt_t i=0; i<fNa; i++)
246 {
247 v1(0,i) = 0;
248 for (UInt_t j=0; j<fNa; j++)
249 v1(0,i) += v3(j,0) * gUnfold.fVacovInv(j,i) ;
250 }
251
252 for (UInt_t i = 0; i<fNa; i++)
253 gUnfold.Chi2(i,0) = v1(0,i) * v3(i,0);
254
255 Double_t chisqsum = 0;
256 for (Int_t i=0; i<gUnfold.Chi2.GetNrows(); i++)
257 chisqsum += gUnfold.Chi2(i, 0);
258 gUnfold.Chisq = chisqsum;
259
260 //-----------------------------------------------------
261 // calculate 2nd derivative squared
262 // regularization term (second derivative squared)
263 gUnfold.SecDeriv = 0;
264 for (UInt_t j=1; j<(fNb-1); j++)
265 {
266 const Double_t temp =
267 + 2.0*(fVb(j+1,0)-fVb(j,0)) / (fVb(j+1,0)+fVb(j,0))
268 - 2.0*(fVb(j,0)-fVb(j-1,0)) / (fVb(j,0)+fVb(j-1,0));
269
270 gUnfold.SecDeriv += temp*temp;
271 }
272
273 gUnfold.ZerDeriv = 0;
274 for (UInt_t j=0; j<fNb; j++)
275 gUnfold.ZerDeriv += fVb(j,0) * fVb(j,0);
276
277 f = gUnfold.Chisq/2 * gUnfold.fW + gUnfold.SecDeriv;
278
279 //*fLog << "F=" << f << " \tSecDeriv=" << gUnfold.SecDeriv
280 // << " \tchi2="
281 // << gUnfold.Chisq << " \tfW=" << gUnfold.fW << endl;
282
283 //--------------------------------------------------------------------
284 // final calculations
285 if (iflag == 3)
286 {
287 //..............................................
288 // calculate external error matrix of the fitted parameters 'val'
289 // extend it with the covariances for y=1-sum(val)
290 Double_t emat[20][20];
291 Int_t ndim = 20;
292 gMinuit->mnemat(&emat[0][0], ndim);
293
294 Double_t covv = 0;
295 for (UInt_t i=0; i<(gUnfold.fNb-1); i++)
296 {
297 Double_t cov = 0;
298 for (UInt_t k=0; k<(gUnfold.fNb-1); k++)
299 cov += emat[i][k];
300
301 emat[i][gUnfold.fNb-1] = -cov;
302 emat[gUnfold.fNb-1][i] = -cov;
303
304 covv += cov;
305 }
306 emat[gUnfold.fNb-1][gUnfold.fNb-1] = covv;
307
308 for (UInt_t i=0; i<gUnfold.fNb; i++)
309 for (UInt_t k=0; k<gUnfold.fNb; k++)
310 gUnfold.fVbcov(i,k) = emat[i][k] *norm*norm;
311
312 //-----------------------------------------------------
313 //..............................................
314 // put unfolded distribution into fResult
315 // fResult(i,0) value in bin i
316 // fResult(i,1) error of value in bin i
317
318 gUnfold.fResult.ResizeTo(gUnfold.fNb, 5);
319
320 Double_t sum = 0;
321 for (UInt_t i=0; i<(gUnfold.fNb-1); i++)
322 {
323 Double_t val;
324 Double_t err;
325 if (!gMinuit->GetParameter(i, val, err))
326 {
327 gLog << "Error getting parameter #" << i << endl;
328 return;
329 }
330
331 Double_t eplus;
332 Double_t eminus;
333 Double_t eparab;
334 Double_t gcc;
335 gMinuit->mnerrs(i, eplus, eminus, eparab, gcc);
336
337 gUnfold.fVb(i, 0) = val * norm;
338
339 gUnfold.fResult(i, 0) = val * norm;
340 gUnfold.fResult(i, 1) = eparab * norm;
341 gUnfold.fResult(i, 2) = eplus * norm;
342 gUnfold.fResult(i, 3) = eminus * norm;
343 gUnfold.fResult(i, 4) = gcc;
344 sum += val;
345 }
346 gUnfold.fVb(gUnfold.fNb-1, 0) = (1.0-sum) * norm;
347
348 gUnfold.fResult(gUnfold.fNb-1, 0) = (1.0-sum) * norm;
349 gUnfold.fResult(gUnfold.fNb-1, 1) =
350 sqrt(gUnfold.fVbcov(gUnfold.fNb-1,gUnfold.fNb-1));
351 gUnfold.fResult(gUnfold.fNb-1, 2) = 0;
352 gUnfold.fResult(gUnfold.fNb-1, 3) = 0;
353 gUnfold.fResult(gUnfold.fNb-1, 4) = 1;
354 //..............................................
355
356 //-----------------------------------------------------
357 // calculate 0th derivative squared
358 gUnfold.ZerDeriv = 0;
359 for (UInt_t j=0; j<fNb; j++)
360 gUnfold.ZerDeriv += fVb(j,0) * fVb(j,0);
361
362 //-----------------------------------------------------
363 // calculate the entropy
364
365 gUnfold.Entropy = 0;
366 for (UInt_t j=0; j<gUnfold.fNb; j++)
367 if (p(j,0) > 0)
368 gUnfold.Entropy += p(j,0) * log( p(j,0) );
369
370
371 //-----------------------------------------------------
372 // calculate SpurSigma
373 gUnfold.SpurSigma = 0.0;
374 for (UInt_t m=0; m<fNb; m++)
375 gUnfold.SpurSigma += gUnfold.fVbcov(m,m);
376 // *fLog << "SpurSigma =" << SpurSigma << endl;
377
378 //-----------------------------------------------------
379 gUnfold.SpurAR = 0;
380 gUnfold.DiffAR2 = 0;
381
382 //-----------------------------------------------------
383 gUnfold.fNdf = gUnfold.fNa;
384 gUnfold.fChisq = gUnfold.Chisq;
385
386 for (UInt_t i=0; i<fNa; i++)
387 {
388 gUnfold.fChi2(i,0) = gUnfold.Chi2(i,0);
389 }
390
391
392 UInt_t iNdf = (UInt_t) (gUnfold.fNdf+0.5);
393
394 //*fLog << "fcnTikhonov2 : fW, chisq (from fcnF) = "
395 // << gUnfold.fW << ", " << gUnfold.fChisq << endl;
396
397 gUnfold.fProb = iNdf>0 ? TMath::Prob(gUnfold.fChisq, iNdf) : 0;
398 }
399}
400
401
402TH1 *MUnfold::DrawMatrixClone(const TMatrixD &m, Option_t *opt="")
403{
404 const Int_t nrows = m.GetNrows();
405 const Int_t ncols = m.GetNcols();
406
407 TMatrix m2(nrows, ncols);
408 for (int i=0; i<nrows; i++)
409 for (int j=0; j<ncols; j++)
410 m2(i, j) = m(i, j);
411
412 TH2F *hist = new TH2F(m2);
413 hist->SetBit(kCanDelete);
414 hist->Draw(opt);
415 hist->SetDirectory(NULL);
416
417 return hist;
418
419}
420
421
422TH1 *MUnfold::DrawMatrixColClone(const TMatrixD &m, Option_t *opt="", Int_t col=0)
423{
424 const Int_t nrows = m.GetNrows();
425
426 TVector vec(nrows);
427 for (int i=0; i<nrows; i++)
428 vec(i) = m(i, col);
429
430 TH1F *hist = new TH1F("TVector","",nrows,0,nrows);
431 for (int i=0; i<nrows; i++)
432 hist->SetBinContent(i+1, vec(i));
433
434 hist->SetBit(kCanDelete);
435 hist->Draw(opt);
436 hist->SetDirectory(NULL);
437
438 return hist;
439}
440
441
442void MUnfold::PrintTH3Content(const TH3 &hist)
443{
444 *fLog << hist.GetName() << ": " << hist.GetTitle() << endl;
445 *fLog << "-----------------------------------------------------" << endl;
446 for (Int_t i=1; i<=hist.GetNbinsX(); i++)
447 {
448 for (Int_t j=1; j<=hist.GetNbinsY(); j++)
449 for (Int_t k=1; k<=hist.GetNbinsZ(); k++)
450 *fLog << hist.GetBinContent(i,j,k) << " \t";
451 *fLog << endl << endl;
452 }
453}
454
455void MUnfold::PrintTH3Error(const TH3 &hist)
456{
457 *fLog << hist.GetName() << ": " << hist.GetTitle() << " <error>" << endl;
458 *fLog << "-----------------------------------------------------" << endl;
459 for (Int_t i=1; i<=hist.GetNbinsX(); i++)
460 {
461 for (Int_t j=1; j<=hist.GetNbinsY(); j++)
462 for (Int_t k=1; k<=hist.GetNbinsZ(); k++)
463 *fLog << hist.GetBinError(i, j, k) << " \t";
464 *fLog << endl << endl;
465 }
466}
467
468void MUnfold::PrintTH2Content(const TH2 &hist)
469{
470 *fLog << hist.GetName() << ": " << hist.GetTitle() << endl;
471 *fLog << "-----------------------------------------------------" << endl;
472 for (Int_t i=1; i<=hist.GetNbinsX(); i++)
473 for (Int_t j=1; j<=hist.GetNbinsY(); j++)
474 *fLog << hist.GetBinContent(i,j) << " \t";
475 *fLog << endl << endl;
476}
477
478void MUnfold::PrintTH2Error(const TH2 &hist)
479{
480 *fLog << hist.GetName() << ": " << hist.GetTitle() << " <error>" << endl;
481 *fLog << "-----------------------------------------------------" << endl;
482 for (Int_t i=1; i<=hist.GetNbinsX(); i++)
483 for (Int_t j=1; j<=hist.GetNbinsY(); j++)
484 *fLog << hist.GetBinError(i, j) << " \t";
485 *fLog << endl << endl;
486}
487
488void MUnfold::PrintTH1Content(const TH1 &hist)
489{
490 *fLog << hist.GetName() << ": " << hist.GetTitle() << endl;
491 *fLog << "-----------------------------------------------------" << endl;
492 for (Int_t i=1; i<=hist.GetNbinsX(); i++)
493 *fLog << hist.GetBinContent(i) << " \t";
494 *fLog << endl << endl;
495}
496
497void MUnfold::PrintTH1Error(const TH1 &hist)
498{
499 *fLog << hist.GetName() << ": " << hist.GetTitle() << " <error>" << endl;
500 *fLog << "-----------------------------------------------------" << endl;
501 for (Int_t i=1; i<=hist.GetNbinsX(); i++)
502 *fLog << hist.GetBinError(i) << " \t";
503 *fLog << endl << endl;
504}
505
506void MUnfold::CopyCol(TMatrixD &m, const TH1 &h, Int_t col=0)
507{
508 const Int_t n = m.GetNrows();
509
510 for (Int_t i=0; i<n; i++)
511 m(i, col) = h.GetBinContent(i+1);
512}
513
514void MUnfold::CopyCol(TH1 &h, const TMatrixD &m, Int_t col=0)
515{
516 const Int_t n = m.GetNrows();
517
518 for (Int_t i=0; i<n; i++)
519 h.SetBinContent(i+1, m(i, col));
520}
521
522void MUnfold::CopyH2M(TMatrixD &m, const TH2 &h)
523{
524 const Int_t nx = m.GetNrows();
525 const Int_t ny = m.GetNcols();
526
527 for (Int_t i=0; i<nx; i++)
528 for (Int_t j=0; j<ny; j++)
529 m(i, j) = h.GetBinContent(i+1, j+1);
530}
531
532void MUnfold::CopySqr(TMatrixD &m, const TH1 &h)
533{
534 const Int_t nx = m.GetNrows();
535 const Int_t ny = m.GetNcols();
536
537 for (Int_t i=0; i<nx; i++)
538 for (Int_t j=0; j<ny; j++)
539 {
540 const Double_t bin = h.GetBinContent(i+1, j+1);
541 m(i, j) = bin*bin;
542 }
543}
544
545Double_t MUnfold::GetMatrixSumRow(const TMatrixD &m, Int_t row)
546{
547 const Int_t n = m.GetNcols();
548
549 Double_t sum = 0;
550 for (Int_t i=0; i<n; i++)
551 sum += m(row, i);
552
553 return sum;
554}
555
556Double_t MUnfold::GetMatrixSumDiag(const TMatrixD &m)
557{
558 const Int_t n = m.GetNcols();
559
560 Double_t sum = 0;
561 for (Int_t i=0; i<n; i++)
562 sum += m(i, i);
563
564 return sum;
565}
566
567Double_t MUnfold::GetMatrixSumCol(const TMatrixD &m, Int_t col=0)
568{
569 const Int_t n = m.GetNrows();
570
571 Double_t sum = 0;
572 for (Int_t i=0; i<n; i++)
573 sum += m(i, col);
574
575 return sum;
576}
577
578Double_t MUnfold::GetMatrixSum(const TMatrixD &m)
579{
580 const Int_t n = m.GetNrows();
581
582 Double_t sum = 0;
583 for (Int_t i=0; i<n; i++)
584 sum += GetMatrixSumRow(m, i);
585
586 return sum;
587}
588
589
590Double_t MUnfold::CalcSpurSigma(TMatrixD &T, Double_t norm)
591{
592 Double_t spursigma = 0;
593
594 for (UInt_t a=0; a<fNb; a++)
595 {
596 for (UInt_t b=0; b<fNb; b++)
597 {
598 fVbcov(a,b) = 0;
599
600 for (UInt_t c=0; c<fNa; c++)
601 for (UInt_t d=0; d<fNa; d++)
602 fVbcov(a,b) += T(a,d)*fVacov(d,c)*T(b,c);
603
604 fVbcov(a,b) *= norm*norm;
605 }
606 spursigma += fVbcov(a,a);
607 }
608
609 return spursigma;
610}
611
612
613MUnfold::MUnfold(TH1D &ha, TH2D &hacov, TH2D &hmig)
614 : fVEps(hmig.GetYaxis()->GetNbins(),1), fVEps0(fVEps, 0)
615 {
616 // ha is the distribution to be unfolded
617 // hacov is the covariance matrix of ha
618 // hmig is the migration matrix;
619 // this matrix will be used in the unfolding
620 // unless SmoothMigrationMatrix(*hmigrat) is called;
621 // in the latter case hmigrat is smoothed
622 // and the smoothed matrix is used in the unfolding
623
624 // Eigen values of the matrix G, which are smaller than EpsLambda
625 // will be considered as being zero
626 EpsLambda = 1.e-10;
627 fW = 0.0;
628
629 fNa = hmig.GetXaxis()->GetNbins();
630 const Double_t alow = hmig.GetXaxis()->GetXmin();
631 const Double_t aup = hmig.GetXaxis()->GetXmax();
632
633 fNb = hmig.GetYaxis()->GetNbins();
634 const Double_t blow = hmig.GetYaxis()->GetXmin();
635 const Double_t bup = hmig.GetYaxis()->GetXmax();
636
637
638 UInt_t Na = ha.GetNbinsX();
639 if (fNa != Na)
640 {
641 *fLog << "MUnfold::MUnfold : dimensions do not match, fNa = ";
642 *fLog << fNa << ", Na = " << Na << endl;
643 }
644
645 *fLog << "MUnfold::MUnfold :" << endl;
646 *fLog << "==================" << endl;
647 *fLog << " fNa = " << fNa << ", fNb = " << fNb << endl;
648
649 // ------------------------
650
651 fVa.ResizeTo(fNa, 1);
652 CopyCol(fVa, ha, 0);
653
654 *fLog << " fVa = ";
655
656 for (UInt_t i=0; i<fNa; i++)
657 *fLog << fVa(i,0) << " \t";
658 *fLog << endl;
659
660 Double_t vaevents = GetMatrixSumCol(fVa, 0);
661 *fLog << " Total number of events in fVa = " << vaevents << endl;
662
663 // ------------------------
664
665 fChi2.ResizeTo(fNa,1);
666 Chi2.ResizeTo(fNa,1);
667
668 // ------------------------
669
670 fVacov.ResizeTo(fNa, fNa);
671 fSpurVacov = 0;
672
673 CopyH2M(fVacov, hacov);
674
675 fVapoints = 0;
676 for (UInt_t i=0; i<fNa; i++)
677 if (fVa(i,0)>0 && fVacov(i,i)<fVa(i,0)*fVa(i,0))
678 fVapoints++;
679
680 fSpurVacov = GetMatrixSumDiag(fVacov);
681
682 //*fLog << "MUnfold::MUnfold : fVacov = " << endl;
683 //*fLog << "==============================" << endl;
684 //fVacov.Print();
685
686 *fLog << " Number of significant points in fVa = ";
687 *fLog << fVapoints << endl;
688
689 *fLog << " Spur of fVacov = ";
690 *fLog << fSpurVacov << endl;
691
692 // ------------------------
693
694 fVacovInv.ResizeTo(fNa, fNa);
695 fVacovInv = fVacov;
696 fVacovInv.InvertPosDef();
697
698 //*fLog << "MUnfold::MUnfold : fVacovInv = " << endl;
699 //*fLog << "==================================" << endl;
700 //fVacovInv.Print();
701
702 // ------------------------
703 // fMigrat is the migration matrix to be used in the unfolding;
704 // fMigrat may be overwritten by SmoothMigrationMatrix
705
706 fMigrat.ResizeTo(fNa, fNb); // row, col
707
708 CopyH2M(fMigrat, hmig);
709
710
711 // ------------------------
712
713 fMigraterr2.ResizeTo(fNa, fNb); // row, col
714 CopySqr(fMigraterr2, hmig);
715
716 // normaxlize
717
718 for (UInt_t j=0; j<fNb; j++)
719 {
720 const Double_t sum = GetMatrixSumCol(fMigrat, j);
721
722 if (sum==0)
723 continue;
724
725 TMatrixDColumn col1(fMigrat, j);
726 col1 *= 1./sum;
727
728 TMatrixDColumn col2(fMigraterr2, j);
729 col2 *= 1./(sum*sum);
730 }
731
732 //*fLog << "MUnfold::MUnfold : fMigrat = " << endl;
733 //*fLog << "===============================" << endl;
734 //fMigrat.Print();
735
736 //*fLog << "MUnfold::MUnfold : fMigraterr2 = " << endl;
737 //*fLog << "===================================" << endl;
738 //fMigraterr2.Print();
739
740 // ------------------------
741 G.ResizeTo(fNa, fNa);
742 EigenValue.ResizeTo(fNa);
743 Eigen.ResizeTo(fNa, fNa);
744
745 fMigOrig.ResizeTo(fNa, fNb);
746 fMigOrigerr2.ResizeTo(fNa, fNb);
747
748 fMigSmoo.ResizeTo (fNa, fNb);
749 fMigSmooerr2.ResizeTo(fNa, fNb);
750 fMigChi2.ResizeTo (fNa, fNb);
751
752 // ------------------------
753
754 fVEps0 = 1./fNb;
755
756 //*fLog << "MUnfold::MUnfold : Default prior distribution fVEps = " << endl;
757 //*fLog << "========================================================" << endl;
758 //fVEps.Print();
759
760 // ------------------------
761
762 fVb.ResizeTo(fNb,1);
763 fVbcov.ResizeTo(fNb,fNb);
764
765 // ----------------------------------------------------
766 // number and range of weights to be scanned
767 Nix = 30;
768 xmin = 1.e-5;
769 xmax = 1.e5;
770 dlogx = (log10(xmax)-log10(xmin)) / Nix;
771
772 SpSig.ResizeTo (Nix);
773 SpAR.ResizeTo (Nix);
774 chisq.ResizeTo (Nix);
775 SecDer.ResizeTo(Nix);
776 ZerDer.ResizeTo(Nix);
777 Entrop.ResizeTo(Nix);
778 DAR2.ResizeTo (Nix);
779 Dsqbar.ResizeTo(Nix);
780
781 //------------------------------------
782 // plots as a function of the iteration number
783
784 hBchisq = new TH1D("Bchisq", "chisq",
785 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
786
787 hBSpAR = new TH1D("BSpAR", "SpurAR",
788 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
789
790 hBDSpAR = new TH1D("BDSpAR", "Delta(SpurAR)",
791 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
792
793 hBSpSig = new TH1D("BSpSig", "SpurSigma/SpurC",
794 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
795
796 hBDSpSig = new TH1D("BDSpSig", "Delta(SpurSigma/SpurC)",
797 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
798
799 hBSecDeriv = new TH1D("BSecDeriv", "Second Derivative squared",
800 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
801
802 hBDSecDeriv = new TH1D("BDSecDeriv", "Delta(Second Derivative squared)",
803 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
804
805 hBZerDeriv = new TH1D("BZerDeriv", "Zero Derivative squared",
806 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
807
808 hBDZerDeriv = new TH1D("BDZerDeriv", "Delta(Zero Derivative squared)",
809 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
810
811 hBEntropy = new TH1D("BEntrop", "Entropy",
812 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
813
814 hBDEntropy = new TH1D("BDEntrop", "Delta(Entropy)",
815 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
816
817 hBDAR2 = new TH1D("BDAR2", "norm(AR-AR+)",
818 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
819
820 hBD2bar = new TH1D("BD2bar", "(b_unfolded-b_ideal)**2",
821 Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );
822
823 //-------------------------------------
824 // original migration matrix
825 fhmig = new TH2D("fMigrat", "Migration matrix",
826 fNa, alow, aup, fNb, blow, bup);
827 fhmig->Sumw2();
828
829 //-------------------------------------
830 // smoothed migration matrix
831 shmig = new TH2D("sMigrat", "Smoothed migration matrix",
832 fNa, alow, aup, fNb, blow, bup);
833 shmig->Sumw2();
834
835 //-------------------------------------
836 // chi2 contributions for smoothing of migration matrix
837 shmigChi2 = new TH2D("sMigratChi2", "Chi2 contr. for smoothing",
838 fNa, alow, aup, fNb, blow, bup);
839
840 //-------------------------------------
841 // eigen values of matrix G = M * M(transposed)
842 hEigen = new TH1D("Eigen", "Eigen values of M*MT",
843 fNa, 0.5, fNa+0.5);
844
845 //------------------------------------
846 // Ideal distribution
847
848 fhb0 = new TH1D("fhb0", "Ideal distribution", fNb, blow, bup);
849 fhb0->Sumw2();
850
851
852 //------------------------------------
853 // Distribution to be unfolded
854 fha = new TH1D("fha", "Distribution to be unfolded", fNa, alow, aup);
855 fha->Sumw2();
856
857 //------------------------------------
858 // Prior distribution
859 hprior = new TH1D("Prior", "Prior distribution", fNb, blow, bup);
860
861 //-----------------------------S-------
862 // Unfolded distribution
863 hb = new TH1D("DataSp", "Unfolded distribution", fNb, blow, bup);
864 hb->Sumw2();
865
866 }
867
868
869// --------------------------------------------------------------------------
870//
871// Default destructor.
872//
873 MUnfold::~MUnfold()
874 {
875// if (hBchisq) delete hBchisq;
876// if (hBSpAR) delete hBSpAR;
877// if (hBDSpAR) delete hBDSpAR;
878// if (hBSpSig) delete hBSpSig;
879// if (hBDSpSig) delete hBDSpSig;
880// if (hBSecDeriv) delete hBSecDeriv;
881// if (hBSecDeriv) delete hBDSecDeriv;
882// if (hBZerDeriv) delete hBZerDeriv;
883// if (hBDZerDeriv) delete hBDZerDeriv;
884// if (hBEntropy) delete hBEntropy;
885// if (hBDEntropy) delete hBDEntropy;
886// if (hBDAR2) delete hBDAR2;
887// if (hBD2bar) delete hBD2bar;
888// if (fhmig) delete fhmig;
889// if (shmig) delete shmig;
890// if (shmigChi2) delete shmigChi2;
891// if (hEigen) delete hEigen;
892// if (fhb0) delete fhb0;
893// if (fha) delete fha;
894// if (hprior) delete hprior;
895// if (hb) delete hb;
896 }
897
898 // -----------------------------------------------------------------------
899 //
900 // Define prior distribution to be a constant
901 //
902 void MUnfold::SetPriorConstant()
903 {
904 fVEps0 = 1./fNb;
905
906 CopyCol(*hprior, fVEps);
907
908 //*fLog << "SetPriorConstant : Prior distribution fVEps = " << endl;
909 //*fLog << "==============================================" << endl;
910 //fVEps.Print();
911 }
912
913 // -----------------------------------------------------------------------
914 //
915 // Take prior distribution from the histogram 'ha'
916 // which may have a different binning than 'hprior'
917 //
918 Bool_t MUnfold::SetPriorRebin(TH1D &ha)
919 {
920
921 // ------------------------------------------------------------------
922 //
923 // fill the contents of histogram 'ha' into the histogram 'hprior';
924 // the histograms need not have the same binning;
925 // if the binnings are different, the bin contents of histogram 'ha'
926 // are distributed properly (linearly) over the bins of 'hprior'
927 //
928
929 const Int_t na = ha.GetNbinsX();
930 const Double_t alow = ha.GetBinLowEdge(1);
931 const Double_t aup = ha.GetBinLowEdge(na+1);
932
933 const Int_t nb = hprior->GetNbinsX();
934 const Double_t blow = hprior->GetBinLowEdge(1);
935 const Double_t bup = hprior->GetBinLowEdge(nb+1);
936
937 // check whether there is an overlap
938 // between the x ranges of the 2 histograms
939 if (alow>bup || aup<blow)
940 {
941 *fLog << "Rebinning not possible because there is no overlap of the x ranges of the two histograms" << endl;
942 return kFALSE;
943 }
944
945 // there is an overlap
946 //********************
947 Double_t sum = 0;
948 for (Int_t j=1; j<=nb; j++)
949 {
950 const Double_t yl = hprior->GetBinLowEdge(j);
951 const Double_t yh = hprior->GetBinLowEdge(j+1);
952
953 // search bins of histogram ha which contribute
954 // to bin j of histogram hb
955 //----------------
956 Int_t il=0;
957 Int_t ih=0;
958 for (Int_t i=2; i<=na+1; i++)
959 {
960 const Double_t xl = ha.GetBinLowEdge(i);
961 if (xl>yl)
962 {
963 il = i-1;
964
965 //.................................
966 ih = 0;
967 for (Int_t k=(il+1); k<=(na+1); k++)
968 {
969 const Double_t xh = ha.GetBinLowEdge(k);
970 if (xh >= yh)
971 {
972 ih = k-1;
973 break;
974 }
975 }
976 //.................................
977 if (ih == 0)
978 ih = na;
979 break;
980 }
981 }
982 //----------------
983 if (il == 0)
984 {
985 *fLog << "Something is wrong " << endl;
986 *fLog << " na, alow, aup = " << na << ", " << alow
987 << ", " << aup << endl;
988 *fLog << " nb, blow, bup = " << nb << ", " << blow
989 << ", " << bup << endl;
990 return kFALSE;
991 }
992
993 Double_t content=0;
994 // sum up the contribution to bin j
995 for (Int_t i=il; i<=ih; i++)
996 {
997 const Double_t xl = ha.GetBinLowEdge(i);
998 const Double_t xh = ha.GetBinLowEdge(i+1);
999 const Double_t bina = xh-xl;
1000
1001 if (xl<yl && xh<yh)
1002 content += ha.GetBinContent(i) * (xh-yl) / bina;
1003 else
1004 if (xl<yl && xh>=yh)
1005 content += ha.GetBinContent(i) * (yh-yl) / bina;
1006 else
1007 if (xl>=yl && xh<yh)
1008 content += ha.GetBinContent(i);
1009 else if (xl>=yl && xh>=yh)
1010 content += ha.GetBinContent(i) * (yh-xl) / bina;
1011 }
1012 hprior->SetBinContent(j, content);
1013 sum += content;
1014 }
1015
1016 // normalize histogram hb
1017 if (sum==0)
1018 {
1019 *fLog << "histogram hb is empty; sum of weights in ha = ";
1020 *fLog << ha.GetSumOfWeights() << endl;
1021 return kFALSE;
1022 }
1023
1024 for (Int_t j=1; j<=nb; j++)
1025 {
1026 const Double_t content = hprior->GetBinContent(j)/sum;
1027 hprior->SetBinContent(j, content);
1028 fVEps0(j-1) = content;
1029 }
1030
1031 //*fLog << "SetPriorRebin : Prior distribution fVEps = " << endl;
1032 //*fLog << "===========================================" << endl;
1033 //fVEps.Print();
1034
1035 return kTRUE;
1036 }
1037
1038
1039 // -----------------------------------------------------------------------
1040 //
1041 // Set prior distribution to a given distribution 'hpr'
1042 //
1043 Bool_t MUnfold::SetPriorInput(TH1D &hpr)
1044 {
1045 CopyCol(fVEps, hpr);
1046
1047 const Double_t sum = GetMatrixSumCol(fVEps, 0);
1048
1049 if (sum<=0)
1050 {
1051 *fLog << "MUnfold::SetPriorInput: invalid prior distribution" << endl;
1052 return kFALSE;
1053 }
1054
1055 // normalize prior distribution
1056 fVEps0 *= 1./sum;
1057
1058 CopyCol(*hprior, fVEps);
1059
1060 //*fLog << "SetPriorInput : Prior distribution fVEps = " << endl;
1061 //*fLog << "===========================================" << endl;
1062 //fVEps.Print();
1063
1064 return kTRUE;
1065 }
1066
1067 // -----------------------------------------------------------------------
1068 //
1069 // Define prior distribution to be a power law
1070 // use input distribution 'hprior' only
1071 // for defining the histogram parameters
1072 //
1073 Bool_t MUnfold::SetPriorPower(Double_t gamma)
1074 {
1075 // generate distribution according to a power law
1076 // dN/dE = E^{-gamma}
1077 // or with y = lo10(E), E = 10^y :
1078 // dN/dy = ln10 * 10^{y*(1-gamma)}
1079 TH1D hpower(*hprior);
1080
1081 const UInt_t nbin = hprior->GetNbinsX();
1082 const Double_t xmin = hprior->GetBinLowEdge(1);
1083 const Double_t xmax = hprior->GetBinLowEdge(nbin+1);
1084
1085 *fLog << "nbin, xmin, xmax = " << nbin << ", ";
1086 *fLog << xmin << ", " << xmax << endl;
1087
1088 TF1* fpow = new TF1("fpow", "pow(10.0, x*(1.0-[0]))", xmin,xmax);
1089 fpow->SetParName (0,"gamma");
1090 fpow->SetParameter(0, gamma );
1091
1092 hpower.FillRandom("fpow", 100000);
1093
1094 // fill prior distribution
1095 CopyCol(fVEps, hpower);
1096
1097 const Double_t sum = GetMatrixSumCol(fVEps, 0);
1098 if (sum <= 0)
1099 {
1100 *fLog << "MUnfold::SetPriorPower : invalid prior distribution" << endl;
1101 return kFALSE;
1102 }
1103
1104 // normalize prior distribution
1105 fVEps0 *= 1./sum;
1106 CopyCol(*hprior, fVEps);
1107
1108 //*fLog << "SetPriorPower : Prior distribution fVEps = " << endl;
1109 //*fLog << "===========================================" << endl;
1110 //fVEps.Print();
1111
1112 return kTRUE;
1113 }
1114
1115
1116 // -----------------------------------------------------------------------
1117 //
1118 // Set the initial weight
1119 //
1120 Bool_t MUnfold::SetInitialWeight(Double_t &weight)
1121 {
1122 if (weight == 0.0)
1123 {
1124 TMatrixD v1(fVa, TMatrixD::kTransposeMult, fVacovInv);
1125 TMatrixD v2(v1, TMatrixD::kMult, fVa);
1126 weight = 1./sqrt(v2(0,0));
1127 }
1128
1129 *fLog << "MUnfold::SetInitialWeight : Initial Weight = "
1130 << weight << endl;
1131
1132 return kTRUE;
1133 }
1134
1135 // -----------------------------------------------------------------------
1136 //
1137 // Print the unfolded distribution
1138 //
1139 void MUnfold::PrintResults()
1140 {
1141 *fLog << bintitle << endl;
1142 *fLog << "PrintResults : Unfolded distribution fResult " << endl;
1143 *fLog << "=============================================" << endl;
1144 //*fLog << "val, eparab, eplus, eminus, gcc = " << endl;
1145
1146 for (UInt_t i=0; i<fNb; i++)
1147 {
1148 // *fLog << fResult(i, 0) << " \t";
1149 // *fLog << fResult(i, 1) << " \t";
1150 // *fLog << fResult(i, 2) << " \t";
1151 // *fLog << fResult(i, 3) << " \t";
1152 // *fLog << fResult(i, 4) << endl;
1153 }
1154 *fLog << "Chisquared, NDF, chi2 probability, ixbest = "
1155 << fChisq << ", "
1156 << fNdf << ", " << fProb << ", " << ixbest << endl;
1157
1158 }
1159
1160
1161 // -----------------------------------------------------------------------
1162 //
1163 // Schmelling : unfolding by minimizing the function Z
1164 // by Gauss-Newton iteration
1165 //
1166 // the weights are scanned between
1167 // 1.e-5*fWinitial and 1.e5*fWinitial
1168 //
1169 Bool_t MUnfold::Schmelling(TH1D &hb0)
1170 {
1171
1172 //======================================================================
1173 // copy ideal distribution
1174 for (UInt_t i=1; i<=fNb; i++)
1175 {
1176 fhb0->SetBinContent(i, hb0.GetBinContent(i));
1177 fhb0->SetBinError (i, hb0.GetBinError(i));
1178 }
1179
1180 //-----------------------------------------------------------------------
1181 // Initialization
1182 // ==============
1183
1184 Int_t numGiteration;
1185 Int_t MaxGiteration = 1000;
1186
1187 TMatrixD alpha;
1188 alpha.ResizeTo(fNa, 1);
1189
1190
1191 //-----------------------------------------------------------------------
1192 // Newton iteration
1193 // ================
1194
1195 Double_t dga2;
1196 Double_t dga2old;
1197 Double_t EpsG = 1.e-12;
1198
1199 TMatrixD wZdp_inv(fNa, fNa);
1200 TMatrixD d(fNb, 1);
1201 TMatrixD p(fNb, 1);
1202
1203 TMatrixD gamma (fNa, 1);
1204 TMatrixD dgamma(fNa, 1);
1205
1206 Double_t fWinitial;
1207 fWinitial = 0.0;
1208 SetInitialWeight(fWinitial);
1209 // for my example this fWinitial was not good; therefore :
1210 fWinitial = 1.0;
1211
1212 Int_t ix;
1213 Double_t xiter;
1214
1215 //-------- start scanning weights --------------------------
1216 // if full == kFALSE only quantities necessary for the Gauss-Newton
1217 // iteration are calculated in SchmellCore
1218 // if full == kTRUE in addition the unfolded distribution,
1219 // its covariance matrix and quantities like
1220 // Chisq, SpurAR, etc. are computed in SchmellCore
1221 //Bool_t full;
1222 //full = kFALSE;
1223 Int_t full;
1224
1225 *fLog << "Schmelling : start loop over weights" << endl;
1226
1227 dga2 = 1.e20;
1228 for (ix=0; ix<Nix; ix++)
1229 {
1230 xiter = pow(10.0,log10(xmin)+ix*dlogx) * fWinitial;
1231
1232 //---------- start Gauss-Newton iteration ----------------------
1233 numGiteration = 0;
1234
1235 // if there was no convergence and the starting gamma was != 0
1236 // redo unfolding for the same weight starting with gamma = 0
1237 //
1238 Int_t gamma0 = 0;
1239 while (1)
1240 {
1241 if (dga2 > EpsG)
1242 {
1243 gamma0 = 1;
1244 gamma.Zero();
1245 }
1246
1247 dga2 = 1.e20;
1248
1249 while (1)
1250 {
1251 dga2old = dga2;
1252
1253 full = 0;
1254 SchmellCore(full, xiter, gamma, dgamma, dga2);
1255
1256 gamma += dgamma;
1257
1258 //*fLog << "Schmelling : ix, numGiteration, dga2, dga2old = "
1259 // << ix << ", " << numGiteration << ", "
1260 // << dga2 << ", " << dga2old << endl;
1261
1262 numGiteration += 1;
1263
1264 // convergence
1265 if (dga2 < EpsG)
1266 break;
1267
1268 // no convergence
1269 if (numGiteration > MaxGiteration)
1270 break;
1271
1272 // gamma doesn't seem to change any more
1273 if (fabs(dga2-dga2old) < EpsG/100.)
1274 break;
1275 }
1276 //---------- end Gauss-Newton iteration ------------------------
1277 if (dga2<EpsG || gamma0 != 0) break;
1278 }
1279
1280 // if Gauss-Newton iteration has not converged
1281 // go to next weight
1282 if (dga2 > EpsG)
1283 {
1284 *fLog << "Schmelling : Gauss-Newton iteration has not converged;"
1285 << " numGiteration = " << numGiteration << endl;
1286 *fLog << " ix, dga2, dga2old = " << ix << ", "
1287 << dga2 << ", " << dga2old << endl;
1288 continue;
1289 }
1290
1291 //*fLog << "Schmelling : Gauss-Newton iteration has converged;" << endl;
1292 //*fLog << "==================================================" << endl;
1293 //*fLog << " numGiteration = " << numGiteration << endl;
1294 //*fLog << " ix, dga2 = " << ix << ", " << dga2 << endl;
1295
1296 // calculated quantities which will be useful for determining
1297 // the best weight (Chisq, SpurAR, ...)
1298 //full = kTRUE;
1299 full = 1;
1300 SchmellCore(full, xiter, gamma, dgamma, dga2);
1301
1302 // calculate difference between ideal and unfolded distribution
1303 Double_t D2bar = 0.0;
1304 for (UInt_t i = 0; i<fNb; i++)
1305 {
1306 Double_t temp = fVb(i,0)-hb0.GetBinContent(i+1,0);
1307 D2bar += temp*temp;
1308 }
1309
1310 SpAR(ix) = SpurAR;
1311 SpSig(ix) = SpurSigma;
1312 chisq(ix) = Chisq;
1313 SecDer(ix) = SecDeriv;
1314 ZerDer(ix) = ZerDeriv;
1315 Entrop(ix) = Entropy;
1316 DAR2(ix) = DiffAR2;
1317 Dsqbar(ix) = D2bar;
1318
1319 }
1320 //---------- end of scanning weights -------------------------------
1321 *fLog << "Schmelling : end of loop over weights" << endl;
1322 // plots ------------------------------
1323 for (ix=0; ix<Nix; ix++)
1324 {
1325 Double_t xbin = log10(xmin)+ix*dlogx;
1326 xiter = pow(10.0,xbin) * fWinitial;
1327
1328 Int_t bin;
1329 bin = hBchisq->FindBin( xbin );
1330 hBchisq->SetBinContent(bin,chisq(ix));
1331 hBSpAR->SetBinContent(bin,SpAR(ix));
1332 hBSpSig->SetBinContent(bin,SpSig(ix)/fSpurVacov);
1333 hBSecDeriv->SetBinContent(bin,SecDer(ix));
1334 hBZerDeriv->SetBinContent(bin,ZerDer(ix));
1335 hBEntropy->SetBinContent(bin,Entrop(ix));
1336 hBDAR2->SetBinContent(bin,DAR2(ix));
1337 hBD2bar->SetBinContent(bin,Dsqbar(ix));
1338
1339 if (ix > 0)
1340 {
1341 Double_t DSpAR = SpAR(ix) - SpAR(ix-1);
1342 hBDSpAR->SetBinContent(bin,DSpAR);
1343
1344 Double_t diff = SpSig(ix) - SpSig(ix-1);
1345 Double_t DSpSig = diff;
1346 hBDSpSig->SetBinContent(bin, DSpSig/fSpurVacov);
1347
1348 Double_t DEntrop = Entrop(ix) - Entrop(ix-1);
1349 hBDEntropy->SetBinContent(bin,DEntrop);
1350
1351 Double_t DSecDer = SecDer(ix) - SecDer(ix-1);
1352 hBDSecDeriv->SetBinContent(bin,DSecDer);
1353
1354 Double_t DZerDer = ZerDer(ix) - ZerDer(ix-1);
1355 hBDZerDeriv->SetBinContent(bin,DZerDer);
1356 }
1357 }
1358
1359 // Select best weight
1360 SelectBestWeight();
1361
1362 if (ixbest < 0.0)
1363 {
1364 *fLog << "Schmelling : no solution found; " << endl;
1365 return kFALSE;
1366 }
1367
1368 // do the unfolding using the best weight
1369 //full = kTRUE;
1370
1371
1372 xiter = pow(10.0,log10(xmin)+ixbest*dlogx) * fWinitial;
1373
1374 //---------- start Gauss-Newton iteration ----------------------
1375 numGiteration = 0;
1376 gamma.Zero();
1377 dga2 = 1.e20;
1378
1379 while (1)
1380 {
1381 full = 1;
1382 SchmellCore(full, xiter, gamma, dgamma, dga2);
1383 gamma += dgamma;
1384
1385 //*fLog << "Schmelling : sum(dgamma^2) = " << dga2 << endl;
1386
1387 numGiteration += 1;
1388
1389 if (numGiteration > MaxGiteration)
1390 break;
1391
1392 if (dga2 < EpsG)
1393 break;
1394 }
1395 //---------- end Gauss-Newton iteration ------------------------
1396
1397
1398 //-----------------------------------------------------------------------
1399 // termination stage
1400 // =================
1401
1402 *fLog << "Schmelling : best solution found; " << endl;
1403 *fLog << "==================================" << endl;
1404 *fLog << " xiter, ixbest, numGiteration, Chisq = "
1405 << xiter << ", " << ixbest << ", "
1406 << numGiteration << ", " << Chisq << endl;
1407
1408 //------------------------------------
1409 //..............................................
1410 // put unfolded distribution into fResult
1411 // fResult(i,0) value in bin i
1412 // fResult(i,1) error of value in bin i
1413
1414 fNdf = SpurAR;
1415 fChisq = Chisq;
1416
1417 for (UInt_t i=0; i<fNa; i++)
1418 {
1419 fChi2(i,0) = Chi2(i,0);
1420 }
1421
1422 UInt_t iNdf = (UInt_t) (fNdf+0.5);
1423 fProb = iNdf>0 ? TMath::Prob(fChisq, iNdf) : 0;
1424
1425 fResult.ResizeTo(fNb, 5);
1426 for (UInt_t i=0; i<fNb; i++)
1427 {
1428 fResult(i, 0) = fVb(i,0);
1429 fResult(i, 1) = sqrt(fVbcov(i,i));
1430 fResult(i, 2) = 0.0;
1431 fResult(i, 3) = 0.0;
1432 fResult(i, 4) = 1.0;
1433 }
1434
1435 //--------------------------------------------------------
1436
1437 *fLog << "Schmelling : gamma = " << endl;
1438 for (UInt_t i=0; i<fNa; i++)
1439 *fLog << gamma(i,0) << " \t";
1440 *fLog << endl;
1441
1442 return kTRUE;
1443 }
1444
1445
1446
1447
1448 // -----------------------------------------------------------------------
1449 //
1450 // SchmellCore main part of Schmellings calculations
1451 //
1452 void MUnfold::SchmellCore(Int_t full, Double_t &xiter, TMatrixD &gamma,
1453 TMatrixD &dgamma, Double_t &dga2)
1454 {
1455 Double_t norm;
1456 TMatrixD wZdp_inv(fNa, fNa);
1457 TMatrixD d(fNb, 1);
1458 TMatrixD p(fNb, 1);
1459
1460 //--------------------------------------------------------
1461 //-- determine the probability vector p
1462
1463
1464 TMatrixD v3(gamma, TMatrixD::kTransposeMult, fMigrat);
1465 TMatrixD v4(TMatrixD::kTransposed, v3);
1466 d = v4;
1467 Double_t dmax = -1.e10;
1468 for (UInt_t j=0; j<fNb; j++)
1469 if (d(j,0)>dmax)
1470 dmax = d(j,0);
1471
1472 Double_t psum = 0.0;
1473 for (UInt_t j=0; j<fNb; j++)
1474 {
1475 d(j,0) -= dmax;
1476 p(j,0) = fVEps0(j)*exp(xiter*d(j,0));
1477 psum += p(j,0);
1478 }
1479
1480 p *= 1.0/psum;
1481
1482 //-- get the vector alpha
1483
1484 TMatrixD alpha(fMigrat, TMatrixD::kMult, p);
1485
1486 //-- determine the current normalization
1487
1488 TMatrixD v2 (alpha, TMatrixD::kTransposeMult, fVacovInv);
1489 TMatrixD normb(v2, TMatrixD::kMult, alpha);
1490
1491 TMatrixD normc(v2, TMatrixD::kMult, fVa);
1492
1493 norm = normc(0,0)/normb(0,0);
1494
1495 //--------------------------------------------------------
1496 //-- determine the scaled slope vector s and Hessian H
1497
1498 TMatrixD Zp(fNa,1);
1499 for (UInt_t i=0; i<fNa; i++)
1500 {
1501 Zp(i,0) = norm*alpha(i,0) - fVa(i,0);
1502 for (UInt_t k=0; k<fNa; k++)
1503 Zp(i,0) += gamma(k,0)*fVacov(k,i);
1504 }
1505
1506
1507 TMatrixD Q (fNa, fNa);
1508 TMatrixD wZdp(fNa, fNa);
1509 for (UInt_t i=0; i<fNa; i++)
1510 {
1511 for (UInt_t j=0; j<fNa; j++)
1512 {
1513 Q(i,j) = - alpha(i,0)*alpha(j,0);
1514 for (UInt_t k=0; k<fNb; k++)
1515 Q(i,j) += fMigrat(i,k)*fMigrat(j,k)*p(k,0);
1516 wZdp(i,j) = xiter*norm*Q(i,j) + fVacov(i,j);
1517 }
1518 }
1519
1520 //-- invert H and calculate the next Newton step
1521
1522 Double_t determ = 1.0;
1523 wZdp_inv = wZdp;
1524 wZdp_inv.Invert(&determ);
1525
1526 if(determ == 0.0)
1527 {
1528 *fLog << "SchmellCore: matrix inversion for H failed" << endl;
1529 return;
1530 }
1531
1532
1533 dga2 = 0.0;
1534 for (UInt_t i=0; i<fNa; i++)
1535 {
1536 dgamma(i,0) = 0.0;
1537 for (UInt_t j=0; j<fNa; j++)
1538 dgamma(i,0) -= wZdp_inv(i,j)*Zp(j,0);
1539 dga2 += dgamma(i,0)*dgamma(i,0);
1540 }
1541
1542 if (full == 0)
1543 return;
1544
1545 //--------------------------------------------------------
1546 //-- determine chi2 and dNdf (#measurements ignored)
1547 Double_t dNdf = 0;
1548 for (UInt_t i=0; i<fNa; i++)
1549 {
1550 Chi2(i,0) = 0;
1551 for (UInt_t j=0; j<fNa; j++)
1552 {
1553 Chi2(i,0) += fVacov(i,j) * gamma(i,0) * gamma(j,0);
1554 dNdf += fVacov(i,j) * wZdp_inv(j,i);
1555 }
1556 }
1557 Chisq = GetMatrixSumCol(Chi2, 0);
1558 SpurAR = fNa - dNdf;
1559
1560 //-----------------------------------------------------
1561 // calculate the norm |AR - AR+|**2
1562
1563 TMatrixD AR(fNa, fNa);
1564 DiffAR2 = 0.0;
1565 for (UInt_t i=0; i<fNa; i++)
1566 {
1567 for (UInt_t j=0; j<fNa; j++)
1568 {
1569 AR(i,j) = 0.0;
1570 for (UInt_t k=0; k<fNa; k++)
1571 AR(i,j) += fVacov(i,k) * wZdp_inv(k,j);
1572 DiffAR2 += AR(i,j) * AR(i,j);
1573 }
1574 }
1575
1576 //--------------------------------------------------------
1577 //-- fill distribution b(*)
1578 fVb = p;
1579 fVb *= norm;
1580
1581 //-- determine the covariance matrix of b (very expensive)
1582
1583 TMatrixD T(fNb,fNa);
1584 for (UInt_t i=0; i<fNb; i++)
1585 {
1586 for (UInt_t j=0; j<fNa; j++)
1587 {
1588 T(i,j) = 0.0;
1589 for (UInt_t k=0; k<fNa; k++)
1590 T(i,j) += xiter*wZdp_inv(k,j)*(fMigrat(k,i)-alpha(k,0))*p(i,0);
1591 }
1592 }
1593
1594 SpurSigma = CalcSpurSigma(T, norm);
1595
1596 //--------------------------------------------------------
1597
1598 //-----------------------------------------------------
1599 // calculate the second derivative squared
1600
1601 SecDeriv = 0;
1602 for (UInt_t j=1; j<(fNb-1); j++)
1603 {
1604 Double_t temp =
1605 + 2.0*(fVb(j+1,0)-fVb(j,0)) / (fVb(j+1,0)+fVb(j,0))
1606 - 2.0*(fVb(j,0)-fVb(j-1,0)) / (fVb(j,0)+fVb(j-1,0));
1607 SecDeriv += temp*temp;
1608 }
1609
1610 ZerDeriv = 0;
1611 for (UInt_t j=0; j<fNb; j++)
1612 ZerDeriv += fVb(j,0) * fVb(j,0);
1613
1614 //-----------------------------------------------------
1615 // calculate the entropy
1616 Entropy = 0;
1617 for (UInt_t j=0; j<fNb; j++)
1618 if (p(j,0) > 0.0)
1619 Entropy += p(j,0) * log( p(j,0) );
1620
1621 //--------------------------------------------------------
1622 }
1623
1624
1625 // -----------------------------------------------------------------------
1626 //
1627 // Smooth migration matrix
1628 // by fitting a function to the migration matrix
1629 //
1630 Bool_t MUnfold::SmoothMigrationMatrix(TH2D &hmigorig)
1631 {
1632 // copy histograms into matrices; the matrices will be used in fcnSmooth
1633 // ------------------------
1634
1635
1636 //*fLog << "MUnfold::SmoothMigrationMatrix : fNa, fNb = " << fNa << ", " << fNb << endl;
1637
1638 //*fLog << "MUnfold::SmoothMigrationMatrix: fMigOrig = " << endl;
1639 //*fLog << "========================================" << endl;
1640 for (UInt_t i=0; i<fNa; i++)
1641 {
1642 for (UInt_t j=0; j<fNb; j++)
1643 {
1644 fMigOrig(i, j) = hmigorig.GetBinContent(i+1, j+1);
1645 //*fLog << fMigOrig(i, j) << " \t";
1646 }
1647 //*fLog << endl;
1648 }
1649
1650
1651 // ------------------------
1652
1653
1654 //*fLog << "MUnfold::SmoothMigrationMatrix : fMigOrigerr2 = " << endl;
1655 //*fLog << "=============================================" << endl;
1656 for (UInt_t i=0; i<fNa; i++)
1657 {
1658 for (UInt_t j=0; j<fNb; j++)
1659 {
1660 fMigOrigerr2(i, j) = hmigorig.GetBinError(i+1, j+1)
1661 * hmigorig.GetBinError(i+1, j+1);
1662
1663 //*fLog << fMigOrigerr2(i, j) << " \t";
1664 }
1665 //*fLog << endl;
1666 }
1667
1668
1669 // ------------------------
1670 // the number of free parameters (npar) is equal to 6:
1671 // a0mean, a1mean, a2mean
1672 // <log10(Eest)> = a0 + a1*log10(Etrue) + a2*SQR(log10(Etrue))
1673 // + log10(Etrue)
1674 // b0RMS, b1RMS, b2RMS
1675 // RMS(log10(Eest)) = b0 + b1*log10(Etrue) + b2*SQR(log10(Etrue))
1676 //
1677 UInt_t npar = 6;
1678
1679 if (npar > 20)
1680 {
1681 *fLog << "MUnfold::SmoothMigrationMatrix : too many parameters, npar = "
1682 << npar << endl;
1683 return kFALSE;
1684 }
1685
1686
1687 //..............................................
1688 // Find reasonable starting values for a0, a1 and b0, b1
1689
1690 Double_t xbar = 0.0;
1691 Double_t xxbar = 0.0;
1692
1693 Double_t ybarm = 0.0;
1694 Double_t xybarm = 0.0;
1695
1696 Double_t ybarR = 0.0;
1697 Double_t xybarR = 0.0;
1698
1699 Double_t Sum = 0.0;
1700 for (UInt_t j=0; j<fNb; j++)
1701 {
1702 Double_t x = (double)j + 0.5;
1703
1704 Double_t meany = 0.0;
1705 Double_t RMSy = 0.0;
1706 Double_t sum = 0.0;
1707 for (UInt_t i=0; i<fNa; i++)
1708 {
1709 Double_t y = (double)i + 0.5;
1710 meany += y * fMigOrig(i, j);
1711 RMSy += y*y * fMigOrig(i, j);
1712 sum += fMigOrig(i, j);
1713 }
1714 if (sum > 0.0)
1715 {
1716 meany = meany / sum;
1717 RMSy = RMSy / sum - meany*meany;
1718 RMSy = sqrt(RMSy);
1719
1720 Sum += sum;
1721 xbar += x * sum;
1722 xxbar += x*x * sum;
1723
1724 ybarm += meany * sum;
1725 xybarm += x*meany * sum;
1726
1727 ybarR += RMSy * sum;
1728 xybarR += x*RMSy * sum;
1729 }
1730 }
1731
1732 if (Sum > 0.0)
1733 {
1734 xbar /= Sum;
1735 xxbar /= Sum;
1736
1737 ybarm /= Sum;
1738 xybarm /= Sum;
1739
1740 ybarR /= Sum;
1741 xybarR /= Sum;
1742 }
1743
1744 Double_t a1start = (xybarm - xbar*ybarm) / (xxbar - xbar*xbar);
1745 Double_t a0start = ybarm - a1start*xbar;
1746 a1start = a1start - 1.0;
1747
1748 Double_t b1start = (xybarR - xbar*ybarR) / (xxbar - xbar*xbar);
1749 Double_t b0start = ybarR - b1start*xbar;
1750
1751 *fLog << "MUnfold::SmoothMigrationMatrix : " << endl;
1752 *fLog << "============================" << endl;
1753 *fLog << "a0start, a1start = " << a0start << ", " << a1start << endl;
1754 *fLog << "b0start, b1start = " << b0start << ", " << b1start << endl;
1755
1756 //..............................................
1757 // Set starting values and step sizes for parameters
1758
1759 char name[20][100];
1760 Double_t vinit[20];
1761 Double_t step[20];
1762 Double_t limlo[20];
1763 Double_t limup[20];
1764 Int_t fix[20];
1765
1766 sprintf(&name[0][0], "a0mean");
1767 vinit[0] = a0start;
1768 //vinit[0] = 1.0;
1769 step[0] = 0.1;
1770 limlo[0] = 0.0;
1771 limup[0] = 0.0;
1772 fix[0] = 0;
1773
1774 sprintf(&name[1][0], "a1mean");
1775 vinit[1] = a1start;
1776 //vinit[1] = 0.0;
1777 step[1] = 0.1;
1778 limlo[1] = 0.0;
1779 limup[1] = 0.0;
1780 fix[1] = 0;
1781
1782 sprintf(&name[2][0], "a2mean");
1783 vinit[2] = 0.0;
1784 step[2] = 0.1;
1785 limlo[2] = 0.0;
1786 limup[2] = 0.0;
1787 fix[2] = 0;
1788
1789 sprintf(&name[3][0], "b0RMS");
1790 vinit[3] = b0start;
1791 //vinit[3] = 0.8;
1792 step[3] = 0.1;
1793 limlo[3] = 1.e-20;
1794 limup[3] = 10.0;
1795 fix[3] = 0;
1796
1797 sprintf(&name[4][0], "b1RMS");
1798 vinit[4] = b1start;
1799 //vinit[4] = 0.0;
1800 step[4] = 0.1;
1801 limlo[4] = 0.0;
1802 limup[4] = 0.0;
1803 fix[4] = 0;
1804
1805 sprintf(&name[5][0], "b2RMS");
1806 vinit[5] = 0.0;
1807 step[5] = 0.1;
1808 limlo[5] = 0.0;
1809 limup[5] = 0.0;
1810 fix[5] = 0;
1811
1812
1813 if ( CallMinuit(fcnSmooth, npar, name, vinit,
1814 step, limlo, limup, fix) )
1815 {
1816
1817 // ------------------------
1818 // fMigrat is the migration matrix to be used in the unfolding;
1819 // fMigrat, as set by the constructor, is overwritten
1820 // by the smoothed migration matrix
1821
1822 for (UInt_t i=0; i<fNa; i++)
1823 for (UInt_t j=0; j<fNb; j++)
1824 fMigrat(i, j) = fMigSmoo(i, j);
1825
1826 // ------------------------
1827
1828 for (UInt_t i=0; i<fNa; i++)
1829 for (UInt_t j=0; j<fNb; j++)
1830 fMigraterr2(i, j) = fMigSmooerr2(i, j);
1831
1832
1833 // normalize
1834 for (UInt_t j=0; j<fNb; j++)
1835 {
1836 Double_t sum = 0.0;
1837 for (UInt_t i=0; i<fNa; i++)
1838 sum += fMigrat(i, j);
1839
1840 //*fLog << "SmoothMigrationMatrix : normalization fMigrat; j, sum + "
1841 // << j << ", " << sum << endl;
1842
1843 if (sum == 0.0)
1844 continue;
1845
1846 for (UInt_t i=0; i<fNa; i++)
1847 {
1848 fMigrat(i, j) /= sum;
1849 fMigraterr2(i, j) /= (sum*sum);
1850 }
1851 }
1852
1853 *fLog << "MUnfold::SmoothMigrationMatrix : fMigrat = " << endl;
1854 *fLog << "========================================" << endl;
1855 for (UInt_t i=0; i<fNa; i++)
1856 {
1857 for (UInt_t j=0; j<fNb; j++)
1858 *fLog << fMigrat(i, j) << " \t";
1859 *fLog << endl;
1860 }
1861
1862 /*
1863 *fLog << "MUnfold::SmoothMigrationMatrix : fMigraterr2 = " << endl;
1864 *fLog << "============================================" << endl;
1865 for (UInt_t i=0; i<fNa; i++)
1866 {
1867 for (UInt_t j=0; j<fNb; j++)
1868 *fLog << fMigraterr2(i, j) << " \t";
1869 *fLog << endl;
1870 }
1871 */
1872
1873 // ------------------------
1874
1875 return kTRUE;
1876 }
1877
1878 return kFALSE;
1879 }
1880
1881 // -----------------------------------------------------------------------
1882 //
1883 // Prepare the call to MINUIT for the minimization of the function
1884 // f = chi2*w/2 + reg, where reg is the regularization term
1885 // reg is the sum the squared 2nd derivatives
1886 // of the unfolded distribution
1887 //
1888 // the corresponding fcn routine is 'fcnTikhonov2'
1889 //
1890 Bool_t MUnfold::Tikhonov2(TH1D &hb0)
1891 {
1892 // the number of free parameters (npar) is equal to
1893 // the number of bins (fNb) of the unfolded distribution minus 1,
1894 // because of the constraint that the total number of events
1895 // is fixed
1896 UInt_t npar = fNb-1;
1897
1898 if (npar > 20)
1899 {
1900 *fLog << "MUnfold::Tikhonov2 : too many parameters, npar = "
1901 << npar << ", fNb = " << fNb << endl;
1902 return kFALSE;
1903 }
1904
1905 // copy ideal distribution
1906
1907 for (UInt_t i=1; i<=fNb; i++)
1908 {
1909 fhb0->SetBinContent(i, hb0.GetBinContent(i));
1910 fhb0->SetBinError (i, hb0.GetBinError(i));
1911 }
1912
1913
1914 //--- start w loop -----------------------------------
1915
1916 *fLog << "Tikhonov2 : start loop over weights" << endl;
1917
1918 Int_t ix;
1919 Double_t xiter;
1920
1921 for (ix=0; ix<Nix; ix++)
1922 {
1923 fW = pow(10.0,log10(xmin)+ix*dlogx);
1924
1925 //..............................................
1926 // Set starting values and step sizes for parameters
1927
1928 char name[20][100];
1929 Double_t vinit[20];
1930 Double_t step[20];
1931 Double_t limlo[20];
1932 Double_t limup[20];
1933 Int_t fix[20];
1934
1935 for (UInt_t i=0; i<npar; i++)
1936 {
1937 sprintf(&name[i][0], "p%d", i+1);
1938 vinit[i] = fVEps0(i);
1939 step[i] = fVEps0(i)/10;
1940
1941 // lower and upper limits (limlo=limup=0: no limits)
1942 //limlo[i] = 1.e-20;
1943 limlo[i] = -1.0;
1944 limup[i] = 1.0;
1945 fix[i] = 0;
1946 }
1947
1948 // calculate solution for the weight fW
1949 // flag non-convergence by chisq(ix) = 0.0
1950 chisq(ix) = 0.0;
1951 if ( CallMinuit(fcnTikhonov2, npar, name, vinit,
1952 step, limlo, limup, fix) )
1953 {
1954 // calculate difference between ideal and unfolded distribution
1955 Double_t D2bar = 0.0;
1956 for (UInt_t i = 0; i<fNb; i++)
1957 {
1958 Double_t temp = fVb(i,0)-hb0.GetBinContent(i+1,0);
1959 D2bar += temp*temp;
1960 }
1961
1962 SpAR(ix) = SpurAR;
1963 SpSig(ix) = SpurSigma;
1964 chisq(ix) = Chisq;
1965 SecDer(ix) = SecDeriv;
1966 ZerDer(ix) = ZerDeriv;
1967 Entrop(ix) = Entropy;
1968 DAR2(ix) = DiffAR2;
1969 Dsqbar(ix) = D2bar;
1970 }
1971 }
1972 *fLog << "Tikhonov2 : end of loop over weights" << endl;
1973
1974 // plots ------------------------------
1975 for (ix=0; ix<Nix; ix++)
1976 {
1977 // test whether minimization has converged
1978 if (chisq(ix) != 0.0)
1979 {
1980 xiter = pow(10.0,log10(xmin)+ix*dlogx);
1981
1982 Int_t bin;
1983 bin = hBchisq->FindBin( log10(xiter) );
1984 hBchisq->SetBinContent(bin,chisq(ix));
1985
1986 //hBSpAR->SetBinContent(bin,SpAR(ix));
1987 hBSpAR->SetBinContent(bin,0.0);
1988
1989 hBSpSig->SetBinContent(bin,SpSig(ix)/fSpurVacov);
1990 hBSecDeriv->SetBinContent(bin,SecDer(ix));
1991 hBZerDeriv->SetBinContent(bin,ZerDer(ix));
1992 hBEntropy->SetBinContent(bin,Entrop(ix));
1993
1994 //hBDAR2->SetBinContent(bin,DAR2(ix));
1995 hBDAR2->SetBinContent(bin,0.0);
1996
1997 hBD2bar->SetBinContent(bin,Dsqbar(ix));
1998
1999 if (ix > 0)
2000 {
2001 //Double_t DSpAR = SpAR(ix) - SpAR(ix-1);
2002 //hBDSpAR->SetBinContent(bin,DSpAR);
2003
2004 Double_t diff = SpSig(ix) - SpSig(ix-1);
2005 Double_t DSpSig = diff;
2006 hBDSpSig->SetBinContent(bin, DSpSig/fSpurVacov);
2007
2008 Double_t DEntrop = Entrop(ix) - Entrop(ix-1);
2009 hBDEntropy->SetBinContent(bin,DEntrop);
2010
2011 Double_t DSecDer = SecDer(ix) - SecDer(ix-1);
2012 hBDSecDeriv->SetBinContent(bin,DSecDer);
2013
2014 Double_t DZerDer = ZerDer(ix) - ZerDer(ix-1);
2015 hBDZerDeriv->SetBinContent(bin,DZerDer);
2016 }
2017 }
2018 }
2019
2020
2021 //--- end w loop -----------------------------------
2022
2023 // Select best weight
2024 SelectBestWeight();
2025
2026 *fLog << " Tikhonov2 : after SelectBestWeight" << endl;
2027
2028 if (ixbest < 0.0)
2029 {
2030 *fLog << "Tikhonov2 : no result found; " << endl;
2031 return kFALSE;
2032 }
2033
2034 *fLog << "Tikhonov2 : best result found; " << endl;
2035 *fLog << "===============================" << endl;
2036 *fLog << " ixbest = " << ixbest << endl;
2037
2038
2039 // do a final unfolding using the best weight
2040
2041 fW = pow(10.0,log10(xmin)+ixbest*dlogx);
2042
2043 //..............................................
2044 // Set starting values and step sizes for parameters
2045
2046 char name[20][100];
2047 Double_t vinit[20];
2048 Double_t step[20];
2049 Double_t limlo[20];
2050 Double_t limup[20];
2051 Int_t fix[20];
2052
2053 for (UInt_t i=0; i<npar; i++)
2054 {
2055 sprintf(&name[i][0], "p%d", i+1);
2056 vinit[i] = fVEps0(i);
2057 step[i] = fVEps0(i)/10;
2058
2059 // lower and upper limits (limlo=limup=0: no limits)
2060 //limlo[i] = 1.e-20;
2061 limlo[i] = -1.0;
2062 limup[i] = 1.0;
2063 fix[i] = 0;
2064 }
2065
2066 // calculate solution for the best weight
2067 CallMinuit(fcnTikhonov2, npar, name, vinit,
2068 step, limlo, limup, fix);
2069
2070
2071 *fLog << "Tikhonov : Values for best weight " << endl;
2072 *fLog << "==================================" << endl;
2073 *fLog << "fW, ixbest, Chisq, SpurAR, SpurSigma, SecDeriv, ZerDeriv, Entrop, DiffAR2, D2bar = " << endl;
2074 *fLog << " " << fW << ", " << ixbest << ", "
2075 << Chisq << ", " << SpurAR << ", "
2076 << SpurSigma << ", " << SecDeriv << ", " << ZerDeriv << ", "
2077 << Entropy << ", " << DiffAR2 << ", "
2078 << Dsqbar(ixbest) << endl;
2079
2080 return kTRUE;
2081
2082 }
2083
2084
2085 // -----------------------------------------------------------------------
2086 //
2087 // Bertero :
2088 //
2089 // the unfolded distribution is calculated iteratively;
2090 // the number of iterations after which the iteration is stopped
2091 // corresponds to the 'weight' in other methods
2092 // a small number of iterations corresponds to strong regularization
2093 // a high number to no regularization
2094 //
2095 // see : M.Bertero, INFN/TC-88/2 (1988)
2096 // V.B.Anykeyev et al., NIM A303 (1991) 350
2097 //
2098 Bool_t MUnfold::Bertero(TH1D &hb0)
2099 {
2100 // copy ideal distribution
2101
2102 for (UInt_t i=1; i<=fNb; i++)
2103 {
2104 fhb0->SetBinContent(i, hb0.GetBinContent(i));
2105 fhb0->SetBinError (i, hb0.GetBinError(i));
2106 }
2107
2108
2109 TMatrixD bold(fNb, 1);
2110 bold.Zero();
2111
2112 //----------------------------------------------------------
2113
2114 Double_t db2 = 1.e20;
2115
2116
2117 TMatrixD aminusaest(fNa, 1);
2118
2119 //------- scan number of iterations -----------------
2120
2121 *fLog << "Bertero : start loop over number of iterations" << endl;
2122
2123 Int_t ix;
2124
2125 for (ix=0; ix<Nix; ix++)
2126 {
2127 Double_t xiter = pow(10.0,log10(xmin)+ix*dlogx);
2128
2129 // calculate solution for the iteration number xiter
2130 BertCore(xiter);
2131
2132 // calculate difference between ideal and unfolded distribution
2133 Double_t D2bar = 0.0;
2134 for (UInt_t i = 0; i<fNb; i++)
2135 {
2136 Double_t temp = fVb(i,0)-hb0.GetBinContent(i+1,0);
2137 D2bar += temp*temp;
2138 }
2139
2140 SpAR(ix) = SpurAR;
2141 SpSig(ix) = SpurSigma;
2142 chisq(ix) = Chisq;
2143 SecDer(ix) = SecDeriv;
2144 ZerDer(ix) = ZerDeriv;
2145 Entrop(ix) = Entropy;
2146 DAR2(ix) = DiffAR2;
2147 Dsqbar(ix) = D2bar;
2148
2149 db2 = 0.0;
2150 for (UInt_t i = 0; i<fNb; i++)
2151 {
2152 Double_t temp = fVb(i,0)-bold(i,0);
2153 db2 += temp*temp;
2154 }
2155 bold = fVb;
2156
2157 //if (db2 < Epsdb2) break;
2158
2159 }
2160
2161 *fLog << "Bertero : end of loop over number of iterations" << endl;
2162
2163 // plots ------------------------------
2164 for (ix=0; ix<Nix; ix++)
2165 {
2166 Double_t xiter = pow(10.0,log10(xmin)+ix*dlogx);
2167
2168 Int_t bin;
2169 bin = hBchisq->FindBin( log10(xiter) );
2170 hBchisq->SetBinContent(bin,chisq(ix));
2171 hBSpAR->SetBinContent(bin,SpAR(ix));
2172 hBSpSig->SetBinContent(bin,SpSig(ix)/fSpurVacov);
2173 hBSecDeriv->SetBinContent(bin,SecDer(ix));
2174 hBZerDeriv->SetBinContent(bin,ZerDer(ix));
2175 hBEntropy->SetBinContent(bin,Entrop(ix));
2176 hBDAR2->SetBinContent(bin,DAR2(ix));
2177 hBD2bar->SetBinContent(bin,Dsqbar(ix));
2178
2179 if (ix > 0)
2180 {
2181 Double_t DSpAR = SpAR(ix) - SpAR(ix-1);
2182 hBDSpAR->SetBinContent(bin,DSpAR);
2183
2184 Double_t diff = SpSig(ix) - SpSig(ix-1);
2185 Double_t DSpSig = diff;
2186 hBDSpSig->SetBinContent(bin, DSpSig/fSpurVacov);
2187
2188 Double_t DEntrop = Entrop(ix) - Entrop(ix-1);
2189 hBDEntropy->SetBinContent(bin,DEntrop);
2190
2191 Double_t DSecDer = SecDer(ix) - SecDer(ix-1);
2192 hBDSecDeriv->SetBinContent(bin,DSecDer);
2193
2194 Double_t DZerDer = ZerDer(ix) - ZerDer(ix-1);
2195 hBDZerDeriv->SetBinContent(bin,DZerDer);
2196 }
2197 }
2198 //------- end of scan of number of iterations -----------------
2199
2200 // Select best weight
2201 SelectBestWeight();
2202
2203
2204 if (ixbest < 0.0)
2205 {
2206 *fLog << "Bertero : weight iteration has NOT converged; " << endl;
2207 return kFALSE;
2208 }
2209
2210 *fLog << "Bertero : weight iteration has converged; " << endl;
2211 *fLog << " ixbest = " << ixbest << endl;
2212
2213
2214 // do a final unfolding using the best weight
2215
2216 // calculate solution for the iteration number xiter
2217 Double_t xiter = pow(10.0,log10(xmin)+ixbest*dlogx);
2218 BertCore(xiter);
2219
2220 *fLog << "Bertero : Values for best weight " << endl;
2221 *fLog << "=================================" << endl;
2222 *fLog << "fW, ixbest, Chisq, SpurAR, SpurSigma, SecDeriv, ZerDeriv, Entrop, DiffAR2, D2bar = " << endl;
2223 *fLog << " " << fW << ", " << ixbest << ", "
2224 << Chisq << ", " << SpurAR << ", "
2225 << SpurSigma << ", " << SecDeriv << ", " << ZerDeriv << ", "
2226 << Entropy << ", " << DiffAR2 << ", "
2227 << Dsqbar(ixbest) << endl;
2228
2229 //----------------
2230
2231 fNdf = SpurAR;
2232 fChisq = Chisq;
2233
2234 for (UInt_t i=0; i<fNa; i++)
2235 {
2236 fChi2(i,0) = Chi2(i,0);
2237 }
2238
2239 UInt_t iNdf = (UInt_t) (fNdf+0.5);
2240 fProb = iNdf>0 ? TMath::Prob(fChisq, iNdf) : 0;
2241
2242
2243 fResult.ResizeTo(fNb, 5);
2244 for (UInt_t i=0; i<fNb; i++)
2245 {
2246 fResult(i, 0) = fVb(i,0);
2247 fResult(i, 1) = sqrt(fVbcov(i,i));
2248 fResult(i, 2) = 0.0;
2249 fResult(i, 3) = 0.0;
2250 fResult(i, 4) = 1.0;
2251 }
2252
2253 return kTRUE;
2254 }
2255
2256 // -----------------------------------------------------------------------
2257 //
2258 // main part of Bertero's calculations
2259 //
2260 Bool_t MUnfold::BertCore(Double_t &xiter)
2261 {
2262 // ignore eigen values which are smaller than EpsLambda
2263 TMatrixD G_inv(fNa, fNa);
2264 TMatrixD Gtil_inv(fNa, fNa);
2265 TMatrixD atil(fNb, fNa);
2266 TMatrixD aminusaest(fNa, 1);
2267
2268 G_inv.Zero();
2269 Gtil_inv.Zero();
2270 SpurAR = 0.0;
2271
2272 // ----- loop over eigen values ------------------
2273 // calculate the approximate inverse of the matrix G
2274 //*fLog << "flaml = " << endl;
2275
2276 UInt_t flagstart = 2;
2277 Double_t flaml=0;
2278
2279 for (UInt_t l=0; l<fNa; l++)
2280 {
2281 if (EigenValue(l) < EpsLambda)
2282 continue;
2283
2284 switch (flagstart)
2285 {
2286 case 1 :
2287 // This is the expression for f(lambda) if the initial C^(0)
2288 // is chosen to be zero
2289 flaml = 1.0 - pow(1.0-tau*EigenValue(l),xiter);
2290 break;
2291
2292 case 2 :
2293 // This is the expression for f(lambda) if the initial C^(0)
2294 // is chosen to be equal to the measured distribution
2295 flaml = 1.0 - pow(1.0-tau*EigenValue(l),xiter)
2296 + EigenValue(l) * pow(1.0-tau*EigenValue(l),xiter);
2297 break;
2298 }
2299
2300 // *fLog << flaml << ", ";
2301
2302 for (UInt_t m=0; m<fNa; m++)
2303 {
2304 for (UInt_t n=0; n<fNa; n++)
2305 {
2306 G_inv(m,n) += 1.0 /EigenValue(l) * Eigen(m,l)*Eigen(n,l);
2307 Gtil_inv(m,n) += flaml/EigenValue(l) * Eigen(m,l)*Eigen(n,l);
2308 }
2309 }
2310 SpurAR += flaml;
2311 }
2312 //*fLog << endl;
2313
2314
2315 //*fLog << "Gtil_inv =" << endl;
2316 //for (Int_t m=0; m<fNa; m++)
2317 //{
2318 // for (Int_t n=0; n<fNa; n++)
2319 // {
2320 // *fLog << Gtil_inv(m,n) << ", ";
2321 // }
2322 // *fLog << endl;
2323 //}
2324
2325 //-----------------------------------------------------
2326 // calculate the unfolded distribution b
2327 TMatrixD v2(fMigrat, TMatrixD::kTransposeMult, Gtil_inv);
2328 atil = v2;
2329 TMatrixD v4(atil, TMatrixD::kMult, fVa);
2330 fVb = v4;
2331
2332 //-----------------------------------------------------
2333 // calculate AR and AR+
2334 TMatrixD AR(v2, TMatrixD::kMult, fMigrat);
2335
2336 TMatrixD v3(fMigrat, TMatrixD::kTransposeMult, G_inv);
2337 TMatrixD ARplus(v3, TMatrixD::kMult, fMigrat);
2338
2339
2340 //-----------------------------------------------------
2341 // calculate the norm |AR - AR+|**2
2342
2343 DiffAR2 = 0.0;
2344 for (UInt_t j=0; j<fNb; j++)
2345 {
2346 for (UInt_t k=0; k<fNb; k++)
2347 {
2348 Double_t tempo = AR(j,k) - ARplus(j,k);
2349 DiffAR2 += tempo*tempo;
2350 }
2351 }
2352
2353 //-----------------------------------------------------
2354 // calculate the second derivative squared
2355
2356 SecDeriv = 0;
2357 for (UInt_t j=1; j<(fNb-1); j++)
2358 {
2359 // temp = ( 2.0*fVb(j,0)-fVb(j-1,0)-fVb(j+1,0) ) / ( 2.0*fVb(j,0) );
2360 Double_t temp = 2.0*(fVb(j+1,0)-fVb(j,0)) / (fVb(j+1,0)+fVb(j,0))
2361 - 2.0*(fVb(j,0)-fVb(j-1,0)) / (fVb(j,0)+fVb(j-1,0));
2362 SecDeriv += temp*temp;
2363 }
2364
2365 ZerDeriv = 0;
2366 for (UInt_t j=0; j<fNb; j++)
2367 ZerDeriv += fVb(j,0) * fVb(j,0);
2368
2369 //-----------------------------------------------------
2370 // calculate the entropy
2371
2372 Double_t sumb = 0.0;
2373 for (UInt_t j=0; j<fNb; j++)
2374 sumb += fVb(j,0);
2375
2376 TMatrixD p(fNb,1);
2377 p = fVb;
2378 if (sumb > 0.0)
2379 p *= 1.0/sumb;
2380
2381 Entropy = 0;
2382 for (UInt_t j=0; j<fNb; j++)
2383 if (p(j,0) > 0.0)
2384 Entropy += p(j,0) * log( p(j,0) );
2385
2386 //-----------------------------------------------------
2387
2388 TMatrixD Gb(fMigrat, TMatrixD::kMult, fVb);
2389 aminusaest = fVa;
2390 aminusaest -= Gb;
2391
2392 TMatrixD v1(1,fNa);
2393 for (UInt_t i=0; i<fNa; i++)
2394 {
2395 v1(0,i) = 0.0;
2396 for (UInt_t j=0; j<fNa; j++)
2397 v1(0,i) += aminusaest(j,0) * fVacovInv(j,i) ;
2398 }
2399
2400 //-----------------------------------------------------
2401 // calculate error matrix fVbcov of unfolded distribution
2402 SpurSigma = CalcSpurSigma(atil);
2403
2404 //-----------------------------------------------------
2405 // calculate the chi squared
2406 for (UInt_t i = 0; i<fNa; i++)
2407 Chi2(i,0) = v1(0,i) * aminusaest(i,0);
2408
2409 Chisq = GetMatrixSumCol(Chi2,0);
2410 return kTRUE;
2411 }
2412
2413
2414 // -----------------------------------------------------------------------
2415 //
2416 // Calculate the matrix G = M * M(transposed)
2417 // and its eigen values and eigen vectors
2418 //
2419 Bool_t MUnfold::CalculateG()
2420 {
2421 // Calculate matrix G = M*M(transposed) (M = migration matrix)
2422 // the matrix Eigen of the eigen vectors of G
2423 // the vector EigenValues of the eigen values of G
2424 // the parameter tau = 1/lambda_max
2425 //
2426 TMatrixD v5(TMatrixD::kTransposed, fMigrat);
2427 //TMatrixD G(fMigrat, TMatrixD::kMult, v5);
2428 G.Mult(fMigrat, v5);
2429
2430 Eigen = G.EigenVectors(EigenValue);
2431
2432 RankG = 0.0;
2433 for (UInt_t l=0; l<fNa; l++)
2434 {
2435 if (EigenValue(l) < EpsLambda) continue;
2436 RankG += 1.0;
2437 }
2438
2439 tau = 1.0 / EigenValue(0);
2440
2441 // *fLog << "eigen values : " << endl;
2442 // for (Int_t i=0; i<fNa; i++)
2443 // {
2444 // *fLog << EigenValue(i) << ", ";
2445 // }
2446 // *fLog << endl;
2447
2448 //*fLog << "eigen vectors : " << endl;
2449 //for (Int_t i=0; i<fNa; i++)
2450 //{
2451 // *fLog << " vector " << i << endl;
2452 // for (Int_t j=0; j<fNa; j++)
2453 // {
2454 // *fLog << Eigen(j,i) << ", ";
2455 // }
2456 // *fLog << endl;
2457 //}
2458 //*fLog << endl;
2459
2460 //*fLog << "G =" << endl;
2461 //for (Int_t m=0; m<fNa; m++)
2462 //{
2463 // for (Int_t n=0; n<fNa; n++)
2464 // {
2465 // *fLog << G(m,n) << ", ";
2466 // }
2467 // *fLog << endl;
2468 //}
2469
2470 return kTRUE;
2471 }
2472
2473 // -----------------------------------------------------------------------
2474 //
2475 // Select the best weight
2476 //
2477 Bool_t MUnfold::SelectBestWeight()
2478 {
2479 //-------------------------------
2480 // select 'best' weight according to some criterion
2481
2482 Int_t ix;
2483
2484 Double_t DiffSpSigmax = -1.e10;
2485 Int_t ixDiffSpSigmax = -1;
2486
2487 Double_t DiffSigpointsmin = 1.e10;
2488 Int_t ixDiffSigpointsmin = -1;
2489
2490 Double_t DiffRankGmin = 1.e10;
2491 Int_t ixDiffRankGmin = -1;
2492
2493 Double_t D2barmin = 1.e10;
2494 Int_t ixD2barmin = -1;
2495
2496 Double_t DiffSpSig1min = 1.e10;
2497 Int_t ixDiffSpSig1min = -1;
2498
2499
2500 Int_t ixmax = -1;
2501
2502 // first loop over all weights :
2503 // find smallest chi2
2504 Double_t chisqmin = 1.e20;
2505 for (ix=0; ix<Nix; ix++)
2506 {
2507 // consider only weights for which
2508 // - unfolding was successful
2509 if (chisq(ix) != 0.0)
2510 {
2511 if (chisq(ix) < chisqmin)
2512 chisqmin = chisq(ix);
2513 }
2514 }
2515 Double_t chisq0 = chisqmin > fVapoints ? chisqmin : fVapoints/2.0;
2516
2517 // second loop over all weights :
2518 // consider only weights for which chisq(ix) < chisq0
2519 ixbest = -1;
2520 for (ix=0; ix<Nix; ix++)
2521 {
2522 if (chisq(ix) != 0.0 && chisq(ix) < 2.0*chisq0)
2523 {
2524 // ixmax = highest weight with successful unfolding
2525 // (Least squares solution)
2526 ixmax = ix;
2527
2528 SpurSigma = SpSig(ix);
2529 SpurAR = SpAR(ix);
2530 Chisq = chisq(ix);
2531 D2bar = Dsqbar(ix);
2532
2533 //----------------------------------
2534 // search weight where SpurSigma changes most
2535 // (as a function of the weight)
2536 if (ix > 0 && chisq(ix-1) != 0.0)
2537 {
2538 Double_t diff = SpSig(ix) - SpSig(ix-1);
2539 if (diff > DiffSpSigmax)
2540 {
2541 DiffSpSigmax = diff;
2542 ixDiffSpSigmax = ix;
2543 }
2544 }
2545
2546 //----------------------------------
2547 // search weight where Chisq is close
2548 // to the number of significant measurements
2549 Double_t DiffSigpoints = fabs(Chisq-fVapoints);
2550
2551 if (DiffSigpoints < DiffSigpointsmin)
2552 {
2553 DiffSigpointsmin = DiffSigpoints;
2554 ixDiffSigpointsmin = ix;
2555 }
2556
2557 //----------------------------------
2558 // search weight where Chisq is close
2559 // to the rank of matrix G
2560 Double_t DiffRankG = fabs(Chisq-RankG);
2561
2562 if (DiffRankG < DiffRankGmin)
2563 {
2564 DiffRankGmin = DiffRankG;
2565 ixDiffRankGmin = ix;
2566 }
2567
2568 //----------------------------------
2569 // search weight where SpurSigma is close to 1.0
2570 Double_t DiffSpSig1 = fabs(SpurSigma/fSpurVacov-1.0);
2571
2572 if (DiffSpSig1 < DiffSpSig1min)
2573 {
2574 DiffSpSig1min = DiffSpSig1;
2575 ixDiffSpSig1min = ix;
2576 }
2577
2578 //----------------------------------
2579 // search weight where D2bar is minimal
2580
2581 if (D2bar < D2barmin)
2582 {
2583 D2barmin = D2bar;
2584 ixD2barmin = ix;
2585 }
2586
2587 //----------------------------------
2588 }
2589 }
2590
2591
2592 // choose solution where increase of SpurSigma is biggest
2593 //if ( DiffSpSigmax > 0.0)
2594 // ixbest = ixDiffSpSigmax;
2595 //else
2596 // ixbest = ixDiffSigpointsmin;
2597
2598 // choose Least Squares Solution
2599 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2600 // ixbest = ixmax;
2601 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2602
2603 // choose weight where chi2 is close to the number of significant
2604 // measurements
2605 // ixbest = ixDiffSigpointsmin;
2606
2607 // choose weight where chi2 is close to the rank of matrix G
2608 // ixbest = ixDiffRankGmin;
2609
2610 // choose weight where chi2 is close to the rank of matrix G
2611 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2612 ixbest = ixDiffSpSig1min;
2613 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2614
2615 *fLog << "SelectBestWeight : ixDiffSpSigmax, DiffSpSigmax = "
2616 << ixDiffSpSigmax << ", " << DiffSpSigmax << endl;
2617 *fLog << "================== ixDiffSigpointsmin, DiffSigpointsmin = "
2618 << ixDiffSigpointsmin << ", " << DiffSigpointsmin << endl;
2619
2620 *fLog << " ixDiffRankGmin, DiffRankGmin = "
2621 << ixDiffRankGmin << ", " << DiffRankGmin << endl;
2622
2623 *fLog << " ixDiffSpSig1min, DiffSpSig1min = "
2624 << ixDiffSpSig1min << ", " << DiffSpSig1min << endl;
2625
2626 *fLog << " ixD2barmin, D2barmin = "
2627 << ixD2barmin << ", " << D2barmin << endl;
2628 *fLog << " ixmax = " << ixmax << endl;
2629 *fLog << " ixbest = " << ixbest << endl;
2630
2631
2632 return kTRUE;
2633 }
2634
2635 // -----------------------------------------------------------------------
2636 //
2637 // Draw the plots
2638 //
2639 Bool_t MUnfold::DrawPlots()
2640 {
2641
2642 // in the plots, mark the weight which has been selected
2643 Double_t xbin = log10(xmin)+ixbest*dlogx;
2644
2645 TMarker *m = new TMarker();
2646 m->SetMarkerSize(1);
2647 m->SetMarkerStyle(20);
2648
2649 //-------------------------------------
2650 // draw the iteration plots
2651 TString ctitle = bintitle;
2652 ctitle += "Plots versus weight";
2653 TCanvas *c = new TCanvas("iter", ctitle, 900, 600);
2654 c->Divide(3,2);
2655
2656 c->cd(1);
2657 hBchisq->Draw();
2658 gPad->SetLogy();
2659 hBchisq->SetXTitle("log10(iteration number)");
2660 hBchisq->SetYTitle("chisq");
2661 m->DrawMarker(xbin, log10(chisq(ixbest)));
2662
2663 c->cd(2);
2664 hBD2bar->Draw();
2665 gPad->SetLogy();
2666 hBD2bar->SetXTitle("log10(iteration number)");
2667 hBD2bar->SetYTitle("(b_unfolded-b_ideal)**2");
2668 m->DrawMarker(xbin, log10(Dsqbar(ixbest)));
2669
2670 /*
2671 c->cd(3);
2672 hBDAR2->Draw();
2673 gPad->SetLogy();
2674 strgx = "log10(iteration number)";
2675 strgy = "norm(AR-AR+)";
2676 hBDAR2->SetXTitle(strgx);
2677 hBDAR2->SetYTitle(strgy);
2678 m->DrawMarker(xbin, log10(DAR2(ixbest)));
2679 */
2680
2681 c->cd(3);
2682 hBSecDeriv->Draw();
2683 hBSecDeriv->SetXTitle("log10(iteration number)");
2684 hBSecDeriv->SetYTitle("Second Derivative squared");
2685 m->DrawMarker(xbin, SecDer(ixbest));
2686
2687 /*
2688 c->cd(8);
2689 hBDSecDeriv->Draw();
2690 strgx = "log10(iteration number)";
2691 strgy = "Delta(Second Derivative squared)";
2692 hBDSecDeriv->SetXTitle(strgx);
2693 hBDSecDeriv->SetYTitle(strgy);
2694 */
2695
2696 /*
2697 c->cd(4);
2698 hBZerDeriv->Draw();
2699 strgx = "log10(iteration number)";
2700 strgy = "Zero Derivative squared";
2701 hBZerDeriv->SetXTitle(strgx);
2702 hBZerDeriv->SetYTitle(strgy);
2703 m->DrawMarker(xbin, ZerDer(ixbest));
2704 */
2705
2706 /*
2707 c->cd(5);
2708 hBDZerDeriv->Draw();
2709 strgx = "log10(iteration number)";
2710 strgy = "Delta(Zero Derivative squared)";
2711 hBDZerDeriv->SetXTitle(strgx);
2712 hBDZerDeriv->SetYTitle(strgy);
2713 */
2714
2715 c->cd(4);
2716 hBSpAR->Draw();
2717 hBSpAR->SetXTitle("log10(iteration number)");
2718 hBSpAR->SetYTitle("SpurAR");
2719 m->DrawMarker(xbin, SpAR(ixbest));
2720
2721
2722 /*
2723 c->cd(11);
2724 hBDSpAR->Draw();
2725 strgx = "log10(iteration number)";
2726 strgy = "Delta(SpurAR)";
2727 hBDSpAR->SetXTitle(strgx);
2728 hBDSpAR->SetYTitle(strgy);
2729 */
2730
2731 c->cd(5);
2732 hBSpSig->Draw();
2733 hBSpSig->SetXTitle("log10(iteration number)");
2734 hBSpSig->SetYTitle("SpurSig/SpurC");
2735 m->DrawMarker(xbin, SpSig(ixbest)/fSpurVacov);
2736
2737 /*
2738 c->cd(14);
2739 hBDSpSig->Draw();
2740 strgx = "log10(iteration number)";
2741 strgy = "Delta(SpurSig/SpurC)";
2742 hBDSpSig->SetXTitle(strgx);
2743 hBDSpSig->SetYTitle(strgy);
2744 */
2745
2746 c->cd(6);
2747 hBEntropy->Draw();
2748 hBEntropy->SetXTitle("log10(iteration number)");
2749 hBEntropy->SetYTitle("Entropy");
2750 m->DrawMarker(xbin, Entrop(ixbest));
2751
2752 /*
2753 c->cd(17);
2754 hBDEntropy->Draw();
2755 strgx = "log10(iteration number)";
2756 strgy = "Delta(Entropy)";
2757 hBDEntropy->SetXTitle(strgx);
2758 hBDEntropy->SetYTitle(strgy);
2759 */
2760
2761 //-------------------------------------
2762
2763 for (UInt_t i=0; i<fNa; i++)
2764 {
2765 fha->SetBinContent(i+1, fVa(i, 0) );
2766 fha->SetBinError (i+1, sqrt(fVacov(i, i)));
2767
2768 for (UInt_t j=0; j<fNb; j++)
2769 {
2770 fhmig->SetBinContent(i+1, j+1, fMigOrig(i, j) );
2771 fhmig->SetBinError (i+1, j+1, sqrt(fMigOrigerr2(i, j)) );
2772
2773 shmig->SetBinContent(i+1, j+1, fMigrat(i, j) );
2774 shmig->SetBinError (i+1, j+1, sqrt(fMigraterr2(i, j)) );
2775 shmigChi2->SetBinContent(i+1, j+1, fMigChi2(i, j) );
2776 }
2777 }
2778
2779 //PrintTH2Content(*shmig);
2780 //PrintTH2Content(*shmigChi2);
2781
2782 //-------------------------------------
2783 CopyCol(*hprior, fVEps);
2784 CopyCol(*hb, fVb);
2785 for (UInt_t i=0; i<fNb; i++)
2786 hb->SetBinError(i+1, sqrt(fVbcov(i, i)));
2787
2788 PrintTH1Content(*hb);
2789 PrintTH1Error(*hb);
2790
2791 //..............................................
2792 for (UInt_t i=0; i<fNa; i++)
2793 hEigen->SetBinContent(i+1, EigenValue(i));
2794
2795 //..............................................
2796 // draw the plots
2797 TString cctitle = bintitle;
2798 cctitle += "Unfolding input";
2799 TCanvas *cc = new TCanvas("input", cctitle, 900, 600);
2800 cc->Divide(3, 2);
2801
2802 // distribution to be unfolded
2803 cc->cd(1);
2804 fha->Draw();
2805 gPad->SetLogy();
2806 fha->SetXTitle("log10(E_{est}/GeV)");
2807 fha->SetYTitle("Counts");
2808
2809 // superimpose unfolded distribution
2810 hb->Draw("*HSAME");
2811 hb->SetMarkerColor(75);
2812 hb->SetLineColor(75);
2813
2814 // prior distribution
2815 cc->cd(2);
2816 hprior->Draw();
2817 gPad->SetLogy();
2818 hprior->SetXTitle("log10(E_{true}/GeV)");
2819 hprior->SetYTitle("Counts");
2820
2821 // migration matrix
2822 cc->cd(3);
2823 fhmig->Draw("box");
2824 fhmig->SetXTitle("log10(E_{est}/GeV)");
2825 fhmig->SetYTitle("log10(E_{true}/GeV)");
2826
2827 // smoothed migration matrix
2828 cc->cd(4);
2829 shmig->Draw("box");
2830 shmig->SetXTitle("log10(E_{est}/GeV)");
2831 shmig->SetYTitle("log10(E_{true}/GeV)");
2832
2833 // chi2 contributions for smoothing
2834 cc->cd(5);
2835 shmigChi2->Draw("box");
2836 shmigChi2->SetXTitle("log10(E_{est}/GeV)");
2837 shmigChi2->SetYTitle("log10(E_{true}/GeV)");
2838
2839 // Eigenvalues of matrix M*M(transposed)
2840 cc->cd(6);
2841 hEigen->Draw();
2842 hEigen->SetXTitle("l");
2843 hEigen->SetYTitle("Eigen values Lambda_l of M*M(transposed)");
2844
2845
2846 //..............................................
2847 // draw the results
2848 TString crtitle = bintitle;
2849 crtitle += "Unfolding results";
2850 TCanvas *cr = new TCanvas("results", crtitle, 600, 600);
2851 cr->Divide(2, 2);
2852
2853 // unfolded distribution
2854 cr->cd(1);
2855 hb->Draw();
2856 gPad->SetLogy();
2857 hb->SetXTitle("log10(E-true/GeV)");
2858 hb->SetYTitle("Counts");
2859
2860
2861 // covariance matrix of unfolded distribution
2862 cr->cd(2);
2863 TH1 *hbcov=DrawMatrixClone(fVbcov, "lego");
2864 hbcov->SetBins(fNb, hb->GetBinLowEdge(1), hb->GetBinLowEdge(fNb+1),
2865 fNb, hb->GetBinLowEdge(1), hb->GetBinLowEdge(fNb+1));
2866
2867 hbcov->SetName("hbcov");
2868 hbcov->SetTitle("Error matrix of distribution hb");
2869 hbcov->SetXTitle("log10(E_{true}/GeV)");
2870 hbcov->SetYTitle("log10(E_{true}/GeV)");
2871
2872
2873 // chi2 contributions
2874 cr->cd(3);
2875 TH1 *hchi2=DrawMatrixColClone(fChi2);
2876 hchi2->SetBins(fNa, fha->GetBinLowEdge(1), fha->GetBinLowEdge(fNa+1));
2877
2878 hchi2->SetName("Chi2");
2879 hchi2->SetTitle("chi2 contributions");
2880 hchi2->SetXTitle("log10(E_{est}/GeV)");
2881 hchi2->SetYTitle("Chisquared");
2882
2883
2884 // ideal distribution
2885
2886 cr->cd(4);
2887 fhb0->Draw();
2888 gPad->SetLogy();
2889 fhb0->SetXTitle("log10(E_{true}/GeV)");
2890 fhb0->SetYTitle("Counts");
2891
2892
2893 // superimpose unfolded distribution
2894 hb->Draw("*Hsame");
2895
2896
2897 return kTRUE;
2898 }
2899
2900
2901 // -----------------------------------------------------------------------
2902 //
2903 // Interface to MINUIT
2904 //
2905 //
2906 Bool_t MUnfold::CallMinuit(
2907 void (*fcnx)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t),
2908 UInt_t npar, char name[20][100],
2909 Double_t vinit[20], Double_t step[20],
2910 Double_t limlo[20], Double_t limup[20], Int_t fix[20])
2911 {
2912 //
2913 // Be carefull: This is not thread safe
2914 //
2915 UInt_t maxpar = 100;
2916
2917 if (npar > maxpar)
2918 {
2919 *fLog << "MUnfold::CallMinuit : too many parameters, npar = " << fNb
2920 << ", maxpar = " << maxpar << endl;
2921 return kFALSE;
2922 }
2923
2924 //..............................................
2925 // Set the maximum number of parameters
2926 TMinuit minuit(maxpar);
2927
2928
2929 //..............................................
2930 // Set the print level
2931 // -1 no output except SHOW comands
2932 // 0 minimum output
2933 // 1 normal output (default)
2934 // 2 additional ouput giving intermediate results
2935 // 3 maximum output, showing progress of minimizations
2936 //
2937 Int_t printLevel = -1;
2938 minuit.SetPrintLevel(printLevel);
2939
2940 //..............................................
2941 // Printout for warnings
2942 // SET WAR print warnings
2943 // SET NOW suppress warnings
2944 Int_t errWarn;
2945 Double_t tmpwar = 0;
2946 minuit.mnexcm("SET NOW", &tmpwar, 0, errWarn);
2947
2948 //..............................................
2949 // Set the address of the minimization function
2950 minuit.SetFCN(fcnx);
2951
2952 //..............................................
2953 // Set starting values and step sizes for parameters
2954 for (UInt_t i=0; i<npar; i++)
2955 {
2956 if (minuit.DefineParameter(i, &name[i][0], vinit[i], step[i],
2957 limlo[i], limup[i]))
2958 {
2959 *fLog << "MUnfold::CallMinuit: Error in defining parameter "
2960 << name << endl;
2961 return kFALSE;
2962 }
2963 }
2964
2965 //..............................................
2966 //Int_t NumPars = minuit.GetNumPars();
2967 //*fLog << "MUnfold::CallMinuit : number of free parameters = "
2968 // << NumPars << endl;
2969
2970 //..............................................
2971 // Minimization
2972 minuit.SetObjectFit(this);
2973
2974 //..............................................
2975 // Error definition :
2976 //
2977 // for2 chisquare function :
2978 // up = 1.0 means calculate 1-standard deviation error
2979 // = 4.0 means calculate 2-standard deviation error
2980 //
2981 // for log(likelihood) function :
2982 // up = 0.5 means calculate 1-standard deviation error
2983 // = 2.0 means calculate 2-standard deviation error
2984 Double_t up = 1.0;
2985 minuit.SetErrorDef(up);
2986
2987
2988
2989 // Int_t errMigrad;
2990 // Double_t tmp = 0;
2991 // minuit.mnexcm("MIGRAD", &tmp, 0, errMigrad);
2992
2993
2994 //..............................................
2995 // fix a parameter
2996 for (UInt_t i=0; i<npar; i++)
2997 {
2998 if (fix[i] > 0)
2999 {
3000 Int_t parNo = i;
3001 minuit.FixParameter(parNo);
3002 }
3003 }
3004
3005 //..............................................
3006 // Set maximum number of iterations (default = 500)
3007 Int_t maxiter = 100000;
3008 minuit.SetMaxIterations(maxiter);
3009
3010 //..............................................
3011 // minimization by the method of Migrad
3012 // Int_t errMigrad;
3013 // Double_t tmp = 0;
3014 // minuit.mnexcm("MIGRAD", &tmp, 0, errMigrad);
3015
3016 //..............................................
3017 // same minimization as by Migrad
3018 // but switches to the SIMPLEX method if MIGRAD fails to converge
3019 Int_t errMinimize;
3020 Double_t tmp = 0;
3021 minuit.mnexcm("MINIMIZE", &tmp, 0, errMinimize);
3022
3023 //..............................................
3024 // check quality of minimization
3025 // istat = 0 covariance matrix not calculated
3026 // 1 diagonal approximation only (not accurate)
3027 // 2 full matrix, but forced positive-definite
3028 // 3 full accurate covariance matrix
3029 // (indication of normal convergence)
3030 Double_t fmin, fedm, errdef;
3031 Int_t npari, nparx, istat;
3032 minuit.mnstat(fmin, fedm, errdef, npari, nparx, istat);
3033
3034 if (errMinimize || istat < 3)
3035 {
3036 *fLog << "MUnfold::CallMinuit : Minimization failed" << endl;
3037 *fLog << " fmin = " << fmin << ", fedm = " << fedm
3038 << ", errdef = " << errdef << ", istat = " << istat
3039 << endl;
3040 return kFALSE;
3041 }
3042
3043 //..............................................
3044 // Minos error analysis
3045 // minuit.mnmnos();
3046
3047 //..............................................
3048 // Print current status of minimization
3049 // if nkode = 0 only function value
3050 // 1 parameter values, errors, limits
3051 // 2 values, errors, step sizes, internal values
3052 // 3 values, errors, step sizes, 1st derivatives
3053 // 4 values, paraboloc errors, MINOS errors
3054
3055 //Int_t nkode = 4;
3056 //minuit.mnprin(nkode, fmin);
3057
3058 //..............................................
3059 // call fcn with IFLAG = 3 (final calculation : calculate p(chi2))
3060 // iflag = 1 initial calculations only
3061 // 2 calculate 1st derivatives and function
3062 // 3 calculate function only
3063 // 4 calculate function + final calculations
3064 const char *command = "CALL";
3065 Double_t iflag = 3;
3066 Int_t errfcn3;
3067 minuit.mnexcm(command, &iflag, 1, errfcn3);
3068
3069 return kTRUE;
3070 }
3071
3072
3073
3074
3075
3076
3077
3078
3079
Note: See TracBrowser for help on using the repository browser.