source: trunk/MagicSoft/Mars/mtemp/mmpi/MUnfoldSpectrum.cc@ 9463

Last change on this file since 9463 was 6452, checked in by rwagner, 20 years ago
*** empty log message ***
File size: 9.4 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without expressed
14! * or implied warranty.
15! *
16!
17! Author(s) : R. Wagner, 02/2004 <mailto:rwagner@mppmu.mpg.de>
18!
19! Copyright: MAGIC Software Development, 2000-2005
20!
21!
22\* ======================================================================== */
23
24/////////////////////////////////////////////////////////////////////////////
25//
26// MUnfoldSpectrum
27//
28// Unfolds a gamma spectrum using the algorithms given in the MUnfold class
29//
30/////////////////////////////////////////////////////////////////////////////
31#include "MUnfoldSpectrum.h"
32
33#include "TH1D.h"
34#include "TH2D.h"
35#include "TH3D.h"
36
37#include "MLog.h"
38#include "MLogManip.h"
39#include "MUnfold.h"
40
41ClassImp(MUnfoldSpectrum);
42using namespace std;
43
44MUnfoldSpectrum::MUnfoldSpectrum()
45 : fUnfoldingMethod(2), fPrior(0)
46{
47
48}
49
50void MUnfoldSpectrum::Calc()
51{
52 // Unfold # Excess events vs. Energy and Theta
53
54 // TH2D* tobeunfolded = hex->GetHist();
55 // TH2D* unfolded = new TH2D(*tobeunfolded);
56 // TH3D* migration = migm->GetHist();
57
58 //const Int_t energyBins = fToBeUnfolded->GetXaxis()->GetNbins();
59 //const Int_t thetaBins = fToBeUnfolded->GetYaxis()->GetNbins();
60 //const TAxis* axisEnergy = fToBeUnfolded->GetXaxis();
61 //const TAxis* axisTheta = fToBeUnfolded->GetYaxis();
62 //cout << "Distribution to be unfolded has " << energyBins
63 // << ", " << thetaBins << " bins" <<endl;
64
65 TAxis &taxis = *fToBeUnfolded->GetYaxis();
66 Int_t numybins = taxis.GetNbins();
67
68 cout << "Processing a total number of " << numybins << " bins. " <<endl;
69 cout << "-------------------" << endl << endl;
70
71 for (Int_t m=1; m<=numybins; m++) {
72 TString bintitle = "Bin ";
73 bintitle += m;
74 bintitle += ": ";
75
76 cout << "Processing " << bintitle << endl;
77
78 // -----------------------------------------
79 // ha : distribution to be unfolded
80
81 TH1D &ha = *fToBeUnfolded->ProjectionX("", m, m, "e");
82 TString title = bintitle;
83 title += "E-est distr. to be unfolded";
84 ha.SetNameTitle("ha", title);
85 TAxis &aaxis = *ha.GetXaxis();
86 Int_t na = aaxis.GetNbins();
87 Double_t alow = aaxis.GetBinLowEdge(1);
88 Double_t aup = aaxis.GetBinLowEdge(na+1);
89
90 cout << ha.GetName() << ": " << ha.GetTitle() << endl;
91 cout << "-----------------------------------------------------" << endl;
92 for (Int_t i=1; i<=ha.GetNbinsX(); i++)
93 cout << ha.GetBinContent(i) << " \t";
94 cout << endl << endl;
95
96 // -----------------------------------------
97 // covariance matrix of the distribution ha
98
99 title = bintitle;
100 title += "Error matrix of distribution ha";
101 TH2D hacov("hacov", title, na, alow, aup, na, alow, aup);
102
103 Double_t errmin = 3.0;
104 for (Int_t i=1; i<=na; i++) {
105 for (Int_t j=1; j<=na; j++)
106 hacov.SetBinContent(i, j, 0.0);
107 const Double_t content = ha.GetBinContent(i);
108 const Double_t error2 = (ha.GetBinError(i))*(ha.GetBinError(i));
109 if (content <= errmin && error2 < errmin)
110 hacov.SetBinContent(i, i, errmin);
111 else
112 hacov.SetBinContent(i, i, error2);
113 }
114
115 //PrintTH2Content(hacov);
116
117 // -----------------------------------------
118 // migration matrix :
119 // x corresponds to measured quantity
120 // y corresponds to true quantity
121
122 // The projection is made for the selected bins only.
123 // To select a bin range along an axis, use TAxis::SetRange, eg
124 // h3.GetYaxis()->SetRange(23,56);
125
126 // taxis->SetRange(m,m);
127 TH2D &hmig = *(TH2D*)fMigrationMatrix->Project3D("yxe");
128 title = bintitle;
129 title += "Migration Matrix";
130 hmig.SetNameTitle("Migrat", title);
131
132 TAxis &aaxismig = *hmig.GetXaxis();
133 Int_t namig = aaxismig.GetNbins();
134
135 if (na != namig) {
136 cout << "doUnfolding : binnings are incompatible; na, namig = "
137 << na << ", " << namig << endl;
138 return;
139 }
140
141 TAxis &baxismig = *hmig.GetYaxis();
142 Int_t nbmig = baxismig.GetNbins();
143 Double_t blow = baxismig.GetBinLowEdge(1);
144 Double_t bup = baxismig.GetBinLowEdge(nbmig+1);
145
146 // -----------------------------------------
147 // dummy ideal distribution
148
149 Int_t nb = nbmig;
150
151 title = bintitle;
152 title += "Dummy Ideal distribution";
153 TH1D hb0("dummyhb0", title, nb, blow, bup);
154 //MH::SetBinning(&hb0, &baxismig);
155 hb0.Sumw2();
156
157 for (Int_t k=1; k<=nb; k++) {
158 hb0.SetBinContent(k, 1.0/nb);
159 hb0.SetBinError (k, 0.1/nb);
160 }
161
162 //PrintTH1Content(hb0);
163
164 // -----------------------------------------
165 // unfolded distribution
166
167 title = bintitle;
168 title += "Unfolded distribution";
169 TH1D hb("hb", title, nb, blow, bup);
170
171 // ha is the distribution to be unfolded
172 // hacov is the covariance matrix of the distribution ha
173 // hmig is the migration matrix;
174 // it is used in the unfolding unless it is overwritten
175 // by SmoothMigrationMatrix by the smoothed migration matrix
176 // hmigor is the migration matrix to be smoothed;
177 // the smoothed migration matrix will be used in the unfolding
178 // hpr the prior distribution
179 // it is only used if SetPriorInput(*hpr) is called
180 // hb unfolded distribution
181
182 // create an MUnfold object;
183 MUnfold unfold(ha, hacov, hmig);
184 unfold.bintitle = bintitle;
185
186 // smooth the migration matrix;
187 // the smoothed migration matrix will be used in the unfolding
188 // hmig is the original (unsmoothed) migration matrix
189 unfold.SmoothMigrationMatrix(hmig);
190
191 // define prior distribution (has always to be defined)
192 // the alternatives are
193 // SetPriorConstant(): isotropic distribution
194 // SetPriorPower(gamma): dN/dE = E^{-gamma}
195 // SetPriorInput(*hpr): the distribution *hpr is used
196 // SetPriorRebin(*ha): use rebinned histogram ha
197 Bool_t errorprior=kTRUE;
198 switch (fPrior)
199 {
200 case 1:
201 unfold.SetPriorConstant();
202 break;
203 case 2:
204 errorprior = unfold.SetPriorPower(fPriorPowerGamma);
205 break;
206 case 3:
207 if (!fPriorInputHist)
208 {
209 cout << "Error: No hpr!" << endl;
210 return;
211 }
212 errorprior = unfold.SetPriorInput(*fPriorInputHist);
213 break;
214 case 4:
215 errorprior = unfold.SetPriorRebin(*fPriorRebinHist);
216 break;
217 }
218 if (!errorprior)
219 {
220 cout << "MUnfoldSpectrum::SetPrior... : failed. fPrior = " ;
221 cout << fPrior << endl;
222 return;
223 }
224
225 // calculate the matrix G = M * M(transposed)
226 // M being the migration matrix
227 unfold.CalculateG();
228
229 switch (fUnfoldingMethod)
230 {
231 case 1: // Schmelling:
232 // minimize the function Z by Gauss-Newton iteration;
233 // the parameters to be fitted are gamma(i) = lambda(i)/w;
234 cout << "Unfolding algorithm : Schmelling" << endl;
235 if (!unfold.Schmelling(hb0)) cout << "MUnfoldSpectrum::Schmelling : failed." << endl;
236 break;
237
238 case 2: // Tikhonov2 :
239 // regularization term is sum of (2nd deriv.)**2 ;
240 // minimization by using MINUIT;
241 // the parameters to be fitted are
242 // the bin contents of the unfolded distribution
243 cout << "Unfolding algorithm : Tikhonov" << endl;
244 if (!unfold.Tikhonov2(hb0)) cout << "MUnfoldSpectrum::Tikhonov2 : failed." << endl;
245 break;
246
247 case 3: // Bertero: minimization by iteration
248 cout << "Unfolding algorithm : Bertero" << endl;
249 if (!unfold.Bertero(hb0)) cout << "MUnfoldSpectrum::Bertero : failed." << endl;
250 break;
251 }
252 unfold.PrintResults();
253 unfold.DrawPlots();
254
255 // get unfolded distribution
256 TMatrixD &Vb = unfold.GetVb();
257 TMatrixD &Vbcov = unfold.GetVbcov();
258 UInt_t fNb = unfold.fNb;
259 for (UInt_t a=0; a<fNb; a++)
260 {
261 hb.SetBinContent(a+1, Vb(a,0));
262 hb.SetBinError(a+1, sqrt(Vbcov(a, a)) );
263 }
264
265 for (Int_t k=1; k<=nb; k++) {
266 Double_t content = hb.GetBinContent(k);
267 Double_t error = hb.GetBinError(k);
268
269 fUnfolded->SetBinContent(k, m, content);
270 fUnfolded->SetBinError(k, m, error);
271
272 //hex->FillBinContent(k, m, content, error);
273 }
274 }
275}
276
277
278// --------------------------------------------------------------------------
279//
280// Default destructor.
281//
282MUnfoldSpectrum::~MUnfoldSpectrum()
283{
284}
285
286
287// -----------------------------------------------------------------------
288//
289// Define prior distribution to be a constant
290//
291void MUnfoldSpectrum::SetPriorConstant()
292{
293 fPrior=1;
294}
295
296void MUnfoldSpectrum::SetPriorRebin(TH1D *ha)
297{
298 fPriorRebinHist=ha;
299 fPrior=2;
300}
301
302void MUnfoldSpectrum::SetPriorInput(TH1D *hpr)
303{
304 fPriorInputHist=hpr;
305 fPrior=3;
306}
307
308void MUnfoldSpectrum::SetPriorPower(Double_t gamma)
309{
310 fPriorPowerGamma=gamma;
311 fPrior=4;
312}
313
Note: See TracBrowser for help on using the repository browser.