source: trunk/MagicSoft/Mars/mtemp/mmpi/macros/fluxMUnfold.C@ 5934

Last change on this file since 5934 was 5554, checked in by rwagner, 21 years ago
*** empty log message ***
File size: 22.1 KB
Line 
1////////////////////////////////////////////////////////////////////////////
2// //
3// This program should be run under root: //
4// root fluxMUnfold.C++ //
5// //
6// Authors: T. Bretz 02/2002 <mailto:tbretz@astro.uni-wuerzburg.de> //
7// W. Wittek 09/2002 <mailto:wittek@mppmu.mpg.de> //
8// R. Wagner 11/2004 <mailto:rwagner@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//
227// SteerUnfold
228//
229void SteerUnfold(TString bintitle,
230 TH1D &ha, TH2D &hacov, TH2D &hmig,
231 TH2D &hmigor, TH1D &hb0, TH1D *hpr,
232 TH1D &hb)
233{
234 // ha is the distribution to be unfolded
235 // hacov is the covariance matrix of the distribution ha
236 // hmig is the migration matrix;
237 // it is used in the unfolding unless it is overwritten
238 // by SmoothMigrationMatrix by the smoothed migration matrix
239 // hmigor is the migration matrix to be smoothed;
240 // the smoothed migration matrix will be used in the unfolding
241 // hpr the prior distribution
242 // it is only used if SetPriorInput(*hpr) is called
243 // hb unfolded distribution
244
245 //..............................................
246 // create an MUnfold object;
247 // fill histograms into vectors and matrices
248
249 MUnfold unfold(ha, hacov, hmig);
250 unfold.bintitle = bintitle;
251
252 //..............................................
253 // smooth the migration matrix;
254 // the smoothed migration matrix will be used in the unfolding
255 // hmig is the original (unsmoothed) migration matrix
256
257 unfold.SmoothMigrationMatrix(hmigor);
258
259 //..............................................
260 // define prior distribution (has always to be defined)
261 // the alternatives are :
262
263 // 1 SetPriorConstant() : isotropic distribution
264 // 2 SetPriorPower(gamma) : dN/dE = E^{-gamma}
265 // 3 SetPriorInput(*hpr): the distribution *hpr is used
266 // 4 SetPriorRebin(*ha) : use rebinned histogram ha
267
268 UInt_t flagprior = 4;
269 cout << "SteerUnfold : flagprior = " << flagprior << endl;
270 cout << "==========================="<< endl;
271
272 Bool_t errorprior=kTRUE;
273 switch (flagprior)
274 {
275 case 1:
276 unfold.SetPriorConstant();
277 break;
278 case 2:
279 errorprior = unfold.SetPriorPower(1.5);
280 break;
281 case 3:
282 if (!hpr)
283 {
284 cout << "Error: No hpr!" << endl;
285 return;
286 }
287 errorprior = unfold.SetPriorInput(*hpr);
288 break;
289 case 4:
290 errorprior = unfold.SetPriorRebin(ha);
291 break;
292 }
293 if (!errorprior)
294 {
295 cout << "MUnfold::SetPrior... : failed. flagprior = " ;
296 cout << flagprior << endl;
297 return;
298 }
299
300 //..............................................
301 // calculate the matrix G = M * M(transposed)
302 // M being the migration matrix
303
304 unfold.CalculateG();
305
306 //..............................................
307 // call steering routine for the actual unfolding;
308 // the alternatives are :
309
310 // 1 Schmelling : minimize the function Z by Gauss-Newton iteration;
311 // the parameters to be fitted are gamma(i) = lambda(i)/w;
312
313 // 2 Tikhonov2 : regularization term is sum of (2nd deriv.)**2 ;
314 // minimization by using MINUIT;
315 // the parameters to be fitted are
316 // the bin contents of the unfolded distribution
317
318 // 3 Bertero: minimization by iteration
319 //
320
321 UInt_t flagunfold = 1;
322 cout << "SteerUnfold : flagunfold = " << flagunfold << endl;
323 cout << "============================" << endl;
324
325
326
327 switch (flagunfold)
328 {
329 case 1: // Schmelling
330 cout << "" << endl;
331 cout << "Unfolding algorithm : Schmelling" << endl;
332 cout << "================================" << endl;
333 if (!unfold.Schmelling(hb0))
334 cout << "MUnfold::Schmelling : failed." << endl;
335 break;
336
337 case 2: // Tikhonov2
338 cout << "" << endl;
339 cout << "Unfolding algorithm : Tikhonov" << endl;
340 cout << "================================" << endl;
341 if (!unfold.Tikhonov2(hb0))
342 cout << "MUnfold::Tikhonov2 : failed." << endl;
343 break;
344
345 case 3: // Bertero
346 cout << "" << endl;
347 cout << "Unfolding algorithm : Bertero" << endl;
348 cout << "================================" << endl;
349 if (!unfold.Bertero(hb0))
350 cout << "MUnfold::Bertero : failed." << endl;
351 break;
352 }
353
354
355 //..............................................
356 // Print fResult
357 unfold.PrintResults();
358
359
360 //..............................................
361 // Draw the plots
362 unfold.DrawPlots();
363
364 //..............................................
365 // get unfolded distribution
366 TMatrixD &Vb = unfold.GetVb();
367 TMatrixD &Vbcov = unfold.GetVbcov();
368
369 UInt_t fNb = unfold.fNb;
370
371 for (UInt_t a=0; a<fNb; a++)
372 {
373 hb.SetBinContent(a+1, Vb(a,0));
374 hb.SetBinError(a+1, sqrt(Vbcov(a, a)) );
375 }
376
377}
378
379//__________________________________________________________________________
380////////////////////////////////////////////////////////////////////////////
381// //
382// doUnfolding (to be called in the analysis) //
383// //
384// arguments : //
385// //
386// INPUT //
387// TH2D &tobeunfolded : no.of excess events and its error //
388// vs. (E-est, Theta) //
389// TH3D &migration : migration matrix (E-est, E_true, Theta) //
390// //
391// OUITPUT //
392// TH2D &unfolded : no.of excess events and its error //
393// vs. (E-true, Theta) //
394// //
395// calls SteerUnfold to do the unfolding //
396// //
397// The "Theta" axis is only used to loop over the bins of theta //
398// and to do the unfolding for each bin of theta. Instead of theta //
399// any other variable (or a dummy variable) may be used. //
400// //
401// //
402////////////////////////////////////////////////////////////////////////////
403
404void doUnfolding(TH2D &tobeunfolded, TH3D &migration, TH2D &unfolded)
405{
406 TAxis &taxis = *tobeunfolded.GetYaxis();
407 Int_t numybins = taxis.GetNbins();
408
409 for (Int_t m=1; m<=numybins; m++)
410 {
411 TString bintitle = "Bin ";
412 bintitle += m;
413 bintitle += ": ";
414
415 // -----------------------------------------
416 // ha : distribution to be unfolded
417
418 TH1D &ha = *tobeunfolded.ProjectionX("", m, m, "e");
419 TString title = bintitle;
420 title += "E-est distr. to be unfolded";
421 ha.SetNameTitle("ha", title);
422 TAxis &aaxis = *ha.GetXaxis();
423 Int_t na = aaxis.GetNbins();
424 Double_t alow = aaxis.GetBinLowEdge(1);
425 Double_t aup = aaxis.GetBinLowEdge(na+1);
426
427 PrintTH1Content(ha);
428 PrintTH1Error(ha);
429
430 // -----------------------------------------
431 // covariance matrix of the distribution ha
432
433 title = bintitle;
434 title += "Error matrix of distribution ha";
435 TH2D hacov("hacov", title, na, alow, aup, na, alow, aup);
436 //MH::SetBinning(&hacov, &aaxis, &aaxis);
437
438 Double_t errmin = 3.0;
439 for (Int_t i=1; i<=na; i++)
440 {
441 for (Int_t j=1; j<=na; j++)
442 {
443 hacov.SetBinContent(i, j, 0.0);
444 }
445 const Double_t content = ha.GetBinContent(i);
446 const Double_t error2 = (ha.GetBinError(i))*(ha.GetBinError(i));
447 if (content <= errmin && error2 < errmin)
448 hacov.SetBinContent(i, i, errmin);
449 else
450 hacov.SetBinContent(i, i, error2);
451 }
452
453 //PrintTH2Content(hacov);
454
455
456 // -----------------------------------------
457 // migration matrix :
458 // x corresponds to measured quantity
459 // y corresponds to true quantity
460 TH2D &hmig = *(TH2D*)migration.Project3D("yxe");
461 title = bintitle;
462 title += "Migration Matrix";
463 hmig.SetNameTitle("Migrat", title);
464
465 TAxis &aaxismig = *hmig.GetXaxis();
466 Int_t namig = aaxismig.GetNbins();
467
468 if (na != namig)
469 {
470 cout << "doUnfolding : binnings are incompatible; na, namig = "
471 << na << ", " << namig << endl;
472 return;
473 }
474
475 TAxis &baxismig = *hmig.GetYaxis();
476 Int_t nbmig = baxismig.GetNbins();
477 Double_t blow = baxismig.GetBinLowEdge(1);
478 Double_t bup = baxismig.GetBinLowEdge(nbmig+1);
479
480 PrintTH2Content(hmig);
481 //PrintTH2Error(hmig);
482
483
484 // -----------------------------------------
485 // dummy ideal distribution
486
487 Int_t nb = nbmig;
488
489 title = bintitle;
490 title += "Dummy Ideal distribution";
491 TH1D hb0("dummyhb0", title, nb, blow, bup);;
492 //MH::SetBinning(&hb0, &baxismig);
493 hb0.Sumw2();
494
495 for (Int_t k=1; k<=nb; k++)
496 {
497 hb0.SetBinContent(k, 1.0/nb);
498 hb0.SetBinError (k, 0.1/nb);
499 }
500
501 //PrintTH1Content(hb0);
502
503 // -----------------------------------------
504 // here the prior distribution can be defined for the call
505 // to SetPriorInput(*hpr)
506
507 title = bintitle;
508 title += "Dummy Prior distribution";
509 TH1D hpr("hpr", title, nb, blow, bup);
510 //MH::SetBinning(&hpr, &baxismig);
511
512 for (Int_t k=1; k<=nb; k++)
513 hpr.SetBinContent(k, 1.0/nb);
514
515 //PrintTH1Content(hpr);
516
517
518 // -----------------------------------------
519 // unfolded distribution
520
521
522 title = bintitle;
523 title += "Unfolded distribution";
524 TH1D hb("hb", title, nb, blow, bup);
525 //MH::SetBinning(&hb, &baxismig);
526
527 // -----------------------------------------
528 SteerUnfold(bintitle, ha, hacov, hmig, hmig, hb0, &hpr, hb);
529
530 for (Int_t k=1; k<=nb; k++)
531 {
532 Double_t content = hb.GetBinContent(k);
533 Double_t error = hb.GetBinError(k);
534
535 unfolded.SetBinContent(k, m, content);
536 unfolded.SetBinError(k, m, error);
537 }
538
539 delete &ha;
540 delete &hmig;
541 }
542
543}
544//========================================================================//
545
546
547////////////////////////////////////////////////////////////////////////////
548// //
549// Main program (for testing purposes) //
550// //
551// defines the ideal distribution (hb0) //
552// defines the migration matrix (hMigrat) //
553// defines the distribution to be unfolded (hVa) //
554// //
555// calls doUnfolding //
556// to do the unfolding //
557// //
558////////////////////////////////////////////////////////////////////////////
559void fluxMUnfold()
560{
561 // -----------------------------------------
562 // migration matrix :
563 // x corresponds to measured quantity
564 // y corresponds to true quantity
565
566 const Int_t na = 13;
567 //const Int_t na = 18;
568 const Axis_t alow = 0.25;
569 const Axis_t aup = 3.50;
570
571 const Int_t nb = 11;
572 //const Int_t nb = 22;
573 const Axis_t blow = 0.50;
574 const Axis_t bup = 3.25;
575
576 const Int_t nc = 1;
577 const Axis_t clow = 0.0;
578 const Axis_t cup = 1.0;
579
580 Int_t m = 1;
581
582 TH3D migration("migration", "Migration Matrix",
583 na, alow, aup, nb, blow, bup, nc, clow, cup);
584 migration.Sumw2();
585
586 // parametrize migration matrix as
587 // <log10(Eest)> = a0 + a1*log10(Etrue) + a2*log10(Etrue)**2
588 // + log10(Etrue)
589 // RMS( log10(Eest) ) = b0 + b1*log10(Etrue) + b2*log10(Etrue)**2
590 Double_t a0 = 0.0;
591 Double_t a1 = 0.0;
592 Double_t a2 = 0.0;
593
594 Double_t b0 = 0.26;
595 Double_t b1 =-0.054;
596 Double_t b2 = 0.0;
597
598 TF1 f2("f2", "gaus(0)", alow, aup);
599 f2.SetParName(0, "ampl");
600 f2.SetParName(1, "mean");
601 f2.SetParName(2, "sigma");
602
603 // loop over log10(Etrue) bins
604 TAxis &yaxis = *migration.GetYaxis();
605 for (Int_t j=1; j<=nb; j++)
606 {
607 Double_t yvalue = yaxis.GetBinCenter(j);
608
609 const Double_t mean = a0 + a1*yvalue + a2*yvalue*yvalue + yvalue;
610 const Double_t sigma = b0 + b1*yvalue + b2*yvalue*yvalue;
611 const Double_t ampl = 1./ ( sigma*TMath::Sqrt(2.0*TMath::Pi()));
612
613 // gaus(0) is a substitute for [0]*exp( -0.5*( (x-[1])/[2] )**2 )
614 f2.SetParameter(0, ampl);
615 f2.SetParameter(1, mean);
616 f2.SetParameter(2, sigma);
617
618 // fill temporary 1-dim histogram with the function
619 // fill the histogram using
620 // - either FillRandom
621 // - or using Freq
622 TH1D htemp("temp", "temp", na, alow, aup);
623 htemp.Sumw2();
624
625 for (Int_t k=0; k<1000000; k++)
626 htemp.Fill(f2.GetRandom());
627
628 // copy it into the migration matrix
629 Double_t sum = 0;
630 for (Int_t i=1; i<=na; i++)
631 {
632 const Stat_t content = htemp.GetBinContent(i);
633 migration.SetBinContent(i, j, m, content);
634 sum += content;
635 }
636
637 // normalize migration matrix
638 if (sum==0)
639 continue;
640
641 for (Int_t i=1; i<=na; i++)
642 {
643 const Stat_t content = migration.GetBinContent(i,j,m);
644 migration.SetBinContent(i,j,m, content/sum);
645 migration.SetBinError (i,j,m, sqrt(content)/sum);
646 }
647 }
648
649 //PrintTH3Content(migration);
650 //PrintTH3Error(migration);
651
652 // -----------------------------------------
653 // ideal distribution
654
655 TH1D hb0("hb0", "Ideal distribution", nb, blow, bup);
656 hb0.Sumw2();
657
658 // fill histogram with random numbers according to
659 // an exponential function dN/dE = E^{-gamma}
660 // or with y = log10(E), E = 10^y :
661 // dN/dy = ln10 * 10^{y*(1-gamma)}
662 TF1 f1("f1", "pow(10.0, x*(1.0-[0]))", blow, bup);
663 f1.SetParName(0,"gamma");
664 f1.SetParameter(0, 2.7);
665
666 // ntimes is the number of entries
667 for (Int_t k=0; k<10000; k++)
668 hb0.Fill(f1.GetRandom());
669
670 // introduce energy threshold at 50 GeV
671
672 const Double_t lgEth = 1.70;
673 const Double_t dlgEth = 0.09;
674
675 for (Int_t j=1; j<=nb; j++)
676 {
677 const Double_t lgE = hb0.GetBinCenter(j);
678 const Double_t c = hb0.GetBinContent(j);
679 const Double_t dc = hb0.GetBinError(j);
680 const Double_t f = 1.0 / (1.0 + exp( -(lgE-lgEth)/dlgEth ));
681
682 hb0.SetBinContent(j, f* c);
683 hb0.SetBinError (j, f*dc);
684 }
685
686 //PrintTH1Content(hb0);
687
688 // -----------------------------------------
689 // generate distribution to be unfolded (ha)
690 // by smearing the ideal distribution (hb0)
691 //
692 TH2D tobeunfolded("tobeunfolded", "Distribution to be unfolded",
693 na, alow, aup, nc, clow, cup);
694 tobeunfolded.Sumw2();
695
696 for (Int_t i=1; i<=na; i++)
697 {
698 Double_t cont = 0;
699 for (Int_t j=1; j<=nb; j++)
700 cont += migration.GetBinContent(i, j, m) * hb0.GetBinContent(j);
701
702 tobeunfolded.SetBinContent(i, m, cont);
703 tobeunfolded.SetBinError(i, m, sqrt(cont));
704 }
705
706 //PrintTH2Content(tobeunfolded);
707 //PrintTH2Error(tobeunfolded);
708
709 // -----------------------------------------
710 // unfolded distribution
711
712 TH2D unfolded("unfolded", "Unfolded distribution",
713 nb, blow, bup, nc, clow, cup);
714 unfolded.Sumw2();
715
716 // -----------------------------------------
717 doUnfolding(tobeunfolded, migration, unfolded);
718
719}
720//========================================================================//
721
722
723
Note: See TracBrowser for help on using the repository browser.