1 | /* ======================================================================== *\
|
---|
2 | !
|
---|
3 | ! *
|
---|
4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
---|
5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
---|
6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
---|
7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
---|
8 | ! *
|
---|
9 | ! * Permission to use, copy, modify and distribute this software and its
|
---|
10 | ! * documentation for any purpose is hereby granted without fee,
|
---|
11 | ! * provided that the above copyright notice appear in all copies and
|
---|
12 | ! * that both that copyright notice and this permission notice appear
|
---|
13 | ! * in supporting documentation. It is provided "as is" without express
|
---|
14 | ! * or implied warranty.
|
---|
15 | ! *
|
---|
16 | !
|
---|
17 | !
|
---|
18 | ! Author(s): Markus Meyer, 02/2005 <mailto:meyer@astro.uni-wuerzburg.de>
|
---|
19 | ! Author(s): Thomas Bretz, 04/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
|
---|
20 | !
|
---|
21 | ! Copyright: MAGIC Software Development, 2000-2005
|
---|
22 | !
|
---|
23 | !
|
---|
24 | \* ======================================================================== */
|
---|
25 |
|
---|
26 | /////////////////////////////////////////////////////////////////////////////
|
---|
27 | //
|
---|
28 | // MHSingleMuon
|
---|
29 | //
|
---|
30 | // This class is a histogram class for a displaying muonparameters.
|
---|
31 | // The folowing histgrams will be plotted:
|
---|
32 | // - Radius (TH1F)
|
---|
33 | // - ArcWidth (TH1F)
|
---|
34 | //
|
---|
35 | // Inputcontainer:
|
---|
36 | // - MGeomCam
|
---|
37 | // - MMuonSearchPar
|
---|
38 | // - MMuonCalibPar
|
---|
39 | //
|
---|
40 | ////////////////////////////////////////////////////////////////////////////
|
---|
41 | #include "MHSingleMuon.h"
|
---|
42 |
|
---|
43 | #include <TF1.h>
|
---|
44 | #include <TMinuit.h>
|
---|
45 | #include <TPad.h>
|
---|
46 | #include <TCanvas.h>
|
---|
47 |
|
---|
48 | #include "MLog.h"
|
---|
49 | #include "MLogManip.h"
|
---|
50 |
|
---|
51 | #include "MBinning.h"
|
---|
52 | #include "MParList.h"
|
---|
53 |
|
---|
54 | #include "MGeomCam.h"
|
---|
55 | #include "MGeomPix.h"
|
---|
56 |
|
---|
57 | #include "MSignalCam.h"
|
---|
58 | #include "MSignalPix.h"
|
---|
59 |
|
---|
60 | #include "MMuonSetup.h"
|
---|
61 | #include "MMuonCalibPar.h"
|
---|
62 | #include "MMuonSearchPar.h"
|
---|
63 |
|
---|
64 | ClassImp(MHSingleMuon);
|
---|
65 |
|
---|
66 | using namespace std;
|
---|
67 |
|
---|
68 | // --------------------------------------------------------------------------
|
---|
69 | //
|
---|
70 | // Setup histograms
|
---|
71 | //
|
---|
72 | MHSingleMuon::MHSingleMuon(const char *name, const char *title) :
|
---|
73 | fSignalCam(0), fMuonSearchPar(0), fGeomCam(0), fMargin(0)
|
---|
74 | {
|
---|
75 | fName = name ? name : "MHSingleMuon";
|
---|
76 | fTitle = title ? title : "Histograms of muon parameters";
|
---|
77 |
|
---|
78 | fHistPhi.SetName("HistPhi");
|
---|
79 | fHistPhi.SetTitle("HistPhi");
|
---|
80 | fHistPhi.SetXTitle("\\phi [\\circ]");
|
---|
81 | fHistPhi.SetYTitle("sum of ADC");
|
---|
82 | fHistPhi.SetDirectory(NULL);
|
---|
83 | fHistPhi.SetFillStyle(4000);
|
---|
84 | fHistPhi.UseCurrentStyle();
|
---|
85 |
|
---|
86 | fHistWidth.SetName("HistWidth");
|
---|
87 | fHistWidth.SetTitle("HistWidth");
|
---|
88 | fHistWidth.SetXTitle("distance from the ring center [\\circ]");
|
---|
89 | fHistWidth.SetYTitle("sum of ADC");
|
---|
90 | fHistWidth.SetDirectory(NULL);
|
---|
91 | fHistWidth.SetFillStyle(4000);
|
---|
92 | fHistWidth.UseCurrentStyle();
|
---|
93 |
|
---|
94 | MBinning bins;
|
---|
95 | bins.SetEdges(21, -180, 180);
|
---|
96 | bins.Apply(fHistPhi);
|
---|
97 |
|
---|
98 | bins.SetEdges(28, 0.3, 1.7);
|
---|
99 | bins.Apply(fHistWidth);
|
---|
100 | }
|
---|
101 |
|
---|
102 | // --------------------------------------------------------------------------
|
---|
103 | //
|
---|
104 | // Setup the Binning for the histograms automatically if the correct
|
---|
105 | // instances of MBinning
|
---|
106 | //
|
---|
107 | Bool_t MHSingleMuon::SetupFill(const MParList *plist)
|
---|
108 | {
|
---|
109 | fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
|
---|
110 | if (!fGeomCam)
|
---|
111 | {
|
---|
112 | *fLog << warn << "MGeomCam not found... abort." << endl;
|
---|
113 | return kFALSE;
|
---|
114 | }
|
---|
115 | fMuonSearchPar = (MMuonSearchPar*)plist->FindObject("MMuonSearchPar");
|
---|
116 | if (!fMuonSearchPar)
|
---|
117 | {
|
---|
118 | *fLog << warn << "MMuonSearchPar not found... abort." << endl;
|
---|
119 | return kFALSE;
|
---|
120 | }
|
---|
121 | fSignalCam = (MSignalCam*)plist->FindObject("MSignalCam");
|
---|
122 | if (!fSignalCam)
|
---|
123 | {
|
---|
124 | *fLog << warn << "MSignalCam not found... abort." << endl;
|
---|
125 | return kFALSE;
|
---|
126 | }
|
---|
127 |
|
---|
128 | MMuonSetup *setup = (MMuonSetup*)const_cast<MParList*>(plist)->FindCreateObj("MMuonSetup");
|
---|
129 | if (!setup)
|
---|
130 | return kFALSE;
|
---|
131 |
|
---|
132 | fMargin = setup->GetMargin()/fGeomCam->GetConvMm2Deg();
|
---|
133 |
|
---|
134 | ApplyBinning(*plist, "ArcPhi", &fHistPhi);
|
---|
135 | ApplyBinning(*plist, "MuonWidth", &fHistWidth);
|
---|
136 |
|
---|
137 | return kTRUE;
|
---|
138 | }
|
---|
139 |
|
---|
140 | // --------------------------------------------------------------------------
|
---|
141 | //
|
---|
142 | // Fill the histograms with data from a MMuonCalibPar and
|
---|
143 | // MMuonSearchPar container.
|
---|
144 | //
|
---|
145 | Bool_t MHSingleMuon::Fill(const MParContainer *par, const Stat_t w)
|
---|
146 | {
|
---|
147 | fHistPhi.Reset();
|
---|
148 | fHistWidth.Reset();
|
---|
149 |
|
---|
150 | const Int_t entries = fSignalCam->GetNumPixels();
|
---|
151 |
|
---|
152 | // the position of the center of a muon ring
|
---|
153 | const Float_t cenx = fMuonSearchPar->GetCenterX();
|
---|
154 | const Float_t ceny = fMuonSearchPar->GetCenterY();
|
---|
155 |
|
---|
156 | for (Int_t i=0; i<entries; i++)
|
---|
157 | {
|
---|
158 | const MSignalPix &pix = (*fSignalCam)[i];
|
---|
159 | const MGeomPix &gpix = (*fGeomCam)[i];
|
---|
160 |
|
---|
161 | const Float_t dx = gpix.GetX() - cenx;
|
---|
162 | const Float_t dy = gpix.GetY() - ceny;
|
---|
163 |
|
---|
164 | const Float_t dist = TMath::Hypot(dx, dy);
|
---|
165 |
|
---|
166 | // if the signal is not near the estimated circle, it is ignored.
|
---|
167 | if (dist < fMuonSearchPar->GetRadius() + fMargin &&
|
---|
168 | dist > fMuonSearchPar->GetRadius() - fMargin)
|
---|
169 | fHistPhi.Fill(TMath::ATan2(dx, dy)*TMath::RadToDeg(), pix.GetNumPhotons());
|
---|
170 |
|
---|
171 | // use only the inner pixles. This is geometry dependent. This has to
|
---|
172 | // be fixed!
|
---|
173 | if(i>397)
|
---|
174 | continue;
|
---|
175 |
|
---|
176 | fHistWidth.Fill(dist*fGeomCam->GetConvMm2Deg(), pix.GetNumPhotons());
|
---|
177 | }
|
---|
178 |
|
---|
179 | // error estimation (temporaly)
|
---|
180 | // The error is estimated from the signal. In order to do so, we have to
|
---|
181 | // once convert the signal from ADC to photo-electron. Then we can get
|
---|
182 | // the fluctuation such as F-factor*sqrt(phe).
|
---|
183 | // Up to now, the error of pedestal is not taken into accout. This is not
|
---|
184 | // of course correct. We will include this soon.
|
---|
185 | const Double_t Ffactor = 1.4;
|
---|
186 | for (Int_t i=0; i<fHistPhi.GetNbinsX()+1; i++)
|
---|
187 | {
|
---|
188 | const Float_t abs = TMath::Abs(fHistPhi.GetBinContent(i));
|
---|
189 | const Float_t rougherr = TMath::Sqrt(abs)*Ffactor;
|
---|
190 |
|
---|
191 | fHistPhi.SetBinError(i, rougherr);
|
---|
192 | }
|
---|
193 |
|
---|
194 | for (Int_t i=0; i<fHistWidth.GetNbinsX()+1; i++)
|
---|
195 | {
|
---|
196 | const Float_t abs = TMath::Abs(fHistWidth.GetBinContent(i));
|
---|
197 | const Float_t rougherr = TMath::Sqrt(abs)*Ffactor;
|
---|
198 |
|
---|
199 | fHistWidth.SetBinError(i, rougherr);
|
---|
200 | }
|
---|
201 |
|
---|
202 | return kTRUE;
|
---|
203 | }
|
---|
204 |
|
---|
205 | // --------------------------------------------------------------------------
|
---|
206 | //
|
---|
207 | // Find the first bins starting at the bin with maximum content in both
|
---|
208 | // directions which are below threshold.
|
---|
209 | // If in a range of half the histogram size in both directions no bin
|
---|
210 | // below the threshold is found, kFALSE is returned.
|
---|
211 | //
|
---|
212 | Bool_t MHSingleMuon::FindRangeAboveThreshold(const TH1F &h, Float_t thres, Int_t &first, Int_t &last) const
|
---|
213 | {
|
---|
214 | const Int_t n = h.GetNbinsX();
|
---|
215 | const Int_t maxbin = h.GetMaximumBin();
|
---|
216 | const Int_t edge = maxbin+n/2;
|
---|
217 |
|
---|
218 | // Search from the peak to the right
|
---|
219 | last = -1;
|
---|
220 | for (Int_t i=maxbin; i<edge; i++)
|
---|
221 | {
|
---|
222 | const Float_t val = h.GetBinContent(i%n + 1);
|
---|
223 | if (val<thres)
|
---|
224 | {
|
---|
225 | last = i%n+1;
|
---|
226 | break;
|
---|
227 | }
|
---|
228 | }
|
---|
229 |
|
---|
230 | // Search from the peak to the left
|
---|
231 | first = -1;
|
---|
232 | for (Int_t i=maxbin+n-1; i>=edge; i--)
|
---|
233 | {
|
---|
234 | const Float_t val = h.GetBinContent(i%n + 1);
|
---|
235 | if (val<thres)
|
---|
236 | {
|
---|
237 | first = i%n+1;
|
---|
238 | break;
|
---|
239 | }
|
---|
240 | }
|
---|
241 |
|
---|
242 | return first>=0 && last>=0;
|
---|
243 | }
|
---|
244 |
|
---|
245 | // --------------------------------------------------------------------------
|
---|
246 | //
|
---|
247 | // Photon distribution along the estimated circle is fitted with theoritical
|
---|
248 | // function in order to get some more information such as Arc Phi and Arc
|
---|
249 | // Length.
|
---|
250 | //
|
---|
251 | Bool_t MHSingleMuon::CalcPhi(Double_t thres, Double_t &peakphi, Double_t &arcphi) const
|
---|
252 | {
|
---|
253 | if (fHistPhi.GetMaximum()<thres)
|
---|
254 | return kFALSE;
|
---|
255 |
|
---|
256 | peakphi = 180.-fHistPhi.GetBinCenter(fHistPhi.GetMaximumBin());
|
---|
257 |
|
---|
258 | // Now find the position at which the peak edges crosses the threshold
|
---|
259 | Int_t first, last;
|
---|
260 |
|
---|
261 | FindRangeAboveThreshold(fHistPhi, thres, first, last);
|
---|
262 |
|
---|
263 | const Int_t n = fHistPhi.GetNbinsX();
|
---|
264 | const Int_t edge = fHistPhi.GetMaximumBin()+n/2;
|
---|
265 | if (first<0)
|
---|
266 | first = (edge-1)%n+1;
|
---|
267 | if (last<0)
|
---|
268 | last = edge%n+1;;
|
---|
269 |
|
---|
270 | const Float_t startfitval = fHistPhi.GetBinLowEdge(first+1);
|
---|
271 | const Float_t endfitval = fHistPhi.GetBinLowEdge(last);
|
---|
272 |
|
---|
273 | arcphi = last-1<first ? 360+(endfitval-startfitval) : endfitval-startfitval;
|
---|
274 |
|
---|
275 | //if (fEnableImpactCalc)
|
---|
276 | // CalcImpact(effbinnum, startfitval, endfitval);
|
---|
277 |
|
---|
278 | return kTRUE;
|
---|
279 | }
|
---|
280 |
|
---|
281 | // --------------------------------------------------------------------------
|
---|
282 | //
|
---|
283 | // Photon distribution of distance from the center of estimated ring is
|
---|
284 | // fitted in order to get some more information such as ARC WIDTH which
|
---|
285 | // can represent to the PSF of our reflector.
|
---|
286 | //
|
---|
287 | // thres: Threshold above zero to determin the edges of the peak which
|
---|
288 | // is used as fit range
|
---|
289 | // width: ArcWidth returned in deg
|
---|
290 | // chi: Chi^2/NDF of the fit
|
---|
291 | //
|
---|
292 | Bool_t MHSingleMuon::CalcWidth(Double_t thres, Double_t &width, Double_t &chi)
|
---|
293 | {
|
---|
294 | Int_t first, last;
|
---|
295 |
|
---|
296 | if (!FindRangeAboveThreshold(fHistWidth, thres, first, last))
|
---|
297 | return kFALSE;
|
---|
298 |
|
---|
299 | // This happens in some cases
|
---|
300 | const Int_t n = fHistWidth.GetNbinsX()/2;
|
---|
301 | const Int_t m = fHistWidth.GetMaximumBin();
|
---|
302 | if (first>last)
|
---|
303 | if (m>n) // If maximum is on the right side of histogram
|
---|
304 | last = n;
|
---|
305 | else
|
---|
306 | first = 0; // If maximum is on the left side of histogram
|
---|
307 |
|
---|
308 | if (last-first<=3)
|
---|
309 | return kFALSE;
|
---|
310 |
|
---|
311 | // Now get the fit range
|
---|
312 | const Float_t startfitval = fHistWidth.GetBinLowEdge(first+1);
|
---|
313 | const Float_t endfitval = fHistWidth.GetBinLowEdge(last);
|
---|
314 |
|
---|
315 | // Setup the function and perform the fit
|
---|
316 | TF1 f1("f1", "gaus", startfitval, endfitval);
|
---|
317 |
|
---|
318 | // Choose starting values as accurate as possible
|
---|
319 | f1.SetParameter(0, fHistWidth.Integral(first+1, last));
|
---|
320 | f1.SetParameter(1, fHistWidth.GetBinCenter(m));
|
---|
321 | f1.SetParameter(2, (endfitval-startfitval)/2);
|
---|
322 |
|
---|
323 | // options : N do not store the function, do not draw
|
---|
324 | // I use integral of function in bin rather than value at bin center
|
---|
325 | // R use the range specified in the function range
|
---|
326 | // Q quiet mode
|
---|
327 | fHistWidth.Fit(&f1, "NQR");
|
---|
328 |
|
---|
329 | chi = f1.GetChisquare()/f1.GetNDF();
|
---|
330 |
|
---|
331 | Double_t err;
|
---|
332 | gMinuit->GetParameter(2, width, err); // get the sigma value
|
---|
333 |
|
---|
334 | return kTRUE;
|
---|
335 | }
|
---|
336 |
|
---|
337 | /*
|
---|
338 | // --------------------------------------------------------------------------
|
---|
339 | //
|
---|
340 | // An impact parameter is calculated by fitting the histogram of photon
|
---|
341 | // distribution along the circle with a theoritical model.
|
---|
342 | // (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.
|
---|
343 | // The function (6) is used here.)
|
---|
344 | //
|
---|
345 | // By default this calculation is suppressed because this calculation is
|
---|
346 | // very time consuming. If you want to calculate an impact parameter,
|
---|
347 | // you can call the function of EnableImpactCalc().
|
---|
348 | //
|
---|
349 | void MMuonCalibParCalc::CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval)
|
---|
350 | {
|
---|
351 | // Fit the distribution with Vacanti function. The function is different
|
---|
352 | // for the impact parameter of inside or outside of our reflector.
|
---|
353 | // Then two different functions are applied to the photon distribution,
|
---|
354 | // and the one which give us smaller chisquare value is taken as a
|
---|
355 | // proper one.
|
---|
356 |
|
---|
357 | Double_t val1,err1,val2,err2;
|
---|
358 | // impact parameter inside mirror radius (8.5m)
|
---|
359 | TString func1;
|
---|
360 | Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
|
---|
361 | tmpval = sin(2.*tmpval)*8.5;
|
---|
362 | func1 += "[0]*";
|
---|
363 | func1 += tmpval;
|
---|
364 | func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
|
---|
365 |
|
---|
366 | TF1 f1("f1",func1,startfitval,endfitval);
|
---|
367 | f1.SetParameters(2000,3,0);
|
---|
368 | f1.SetParLimits(1,0,8.5);
|
---|
369 | f1.SetParLimits(2,-180.,180.);
|
---|
370 |
|
---|
371 | fMuonCalibPar->fHistPhi->Fit("f1","RQEM");
|
---|
372 |
|
---|
373 | Float_t chi1 = -1;
|
---|
374 | Float_t chi2 = -1;
|
---|
375 | if(effbinnum>3)
|
---|
376 | chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
|
---|
377 |
|
---|
378 | gMinuit->GetParameter(1,val1,err1); // get the estimated IP
|
---|
379 | Float_t estip1 = val1;
|
---|
380 |
|
---|
381 | // impact parameter beyond mirror area (8.5m)
|
---|
382 | TString func2;
|
---|
383 | Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
|
---|
384 | tmpval2 = sin(2.*tmpval2)*8.5*2.;
|
---|
385 | func2 += "[0]*";
|
---|
386 | func2 += tmpval2;
|
---|
387 | func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
|
---|
388 |
|
---|
389 | TF1 f2("f2",func2,startfitval,endfitval);
|
---|
390 | f2.SetParameters(2000,20,0);
|
---|
391 | f2.SetParLimits(1,8.5,300.);
|
---|
392 | f2.SetParLimits(2,-180.,180.);
|
---|
393 |
|
---|
394 | fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");
|
---|
395 |
|
---|
396 | if(effbinnum>3)
|
---|
397 | chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
|
---|
398 |
|
---|
399 | gMinuit->GetParameter(1,val2,err2); // get the estimated IP
|
---|
400 | Float_t estip2 = val2;
|
---|
401 |
|
---|
402 | if(effbinnum<=3)
|
---|
403 | {
|
---|
404 | fMuonCalibPar->SetEstImpact(-1.);
|
---|
405 | }
|
---|
406 | if(chi2 > chi1)
|
---|
407 | {
|
---|
408 | fMuonCalibPar->SetEstImpact(estip1);
|
---|
409 | fMuonCalibPar->SetChiArcPhi(chi1);
|
---|
410 | }
|
---|
411 | else
|
---|
412 | {
|
---|
413 | fMuonCalibPar->SetEstImpact(estip2);
|
---|
414 | fMuonCalibPar->SetChiArcPhi(chi2);
|
---|
415 | }
|
---|
416 | }
|
---|
417 | */
|
---|