source: trunk/MagicSoft/Mars/macros/fluxunfold.C@ 5483

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