source: trunk/MagicSoft/Mars/macros/unfold.C@ 5868

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