source: trunk/MagicSoft/Mars/mmuon/MHSingleMuon.cc@ 7327

Last change on this file since 7327 was 7225, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 14.7 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 express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Keiichi Mase, 10/2004
19! Author(s): Markus Meyer, 02/2005 <mailto:meyer@astro.uni-wuerzburg.de>
20! Author(s): Thomas Bretz, 04/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2005
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28//
29// MHSingleMuon
30//
31// This class is a histogram class for displaying the radial (fHistWidth)
32// and the azimuthal (fHistPhi) intensity distribution for one muon.
33// You can retrieve the histogram (TH1F) using the function GetHistPhi()
34// or GetHistWidth().
35// From these histograms the fraction of the ring segment (ArcPhi) and the
36// Width of the muon ring (ArcWidth) is calculated.
37//
38// First, the radius and center of the ring has to be calculted by
39// MMuonSearchParCalc
40// After that the histograms has to be filled in the following way:
41//
42// MFillH fillmuon("MHSingleMuon", "", "FillMuon");
43//
44// The allowed region to estimate ArcPhi is a certain margin around the
45// radius. The default value is 0.2 deg (60mm). If the estimated radius
46// of the arc is 1.0 deg, the pixel contents in the radius range from
47// 0.8 deg to 1.2 deg are fill in the histogram.
48//
49// For ArcPhi only bins over a certain threshold are supposed to be part
50// of the ring.
51// For ArcWidth, the same algorithm is used to determine the fit region
52// for a gaussian fit to the radial intensity distribution. The ArcWidth
53// is defined as the sigma value of the gaussian fit.
54//
55// The binning of the histograms can be changed in the following way:
56//
57// MBinning bins1("BinningMuonWidth");
58// MBinning bins2("BinningArcPhi");
59// bins1.SetEdges(28, 0.3, 1.7);
60// bins2.SetEdges(20, -180,180);
61// plist.AddToList(&bins1);
62// plist.AddToList(&bins2);
63//
64// The values for the thresholds and the margin are saved in MMuonSetup.
65// They can be easily changed in star.rc.
66//
67// Please have in mind, that changes in this basic parameters will change
68// your results!!
69//
70// Inputcontainer:
71// - MGeomCam
72// - MMuonSearchPar
73//
74//
75////////////////////////////////////////////////////////////////////////////
76#include "MHSingleMuon.h"
77
78#include <TF1.h>
79#include <TMinuit.h>
80#include <TPad.h>
81#include <TCanvas.h>
82
83#include "MLog.h"
84#include "MLogManip.h"
85
86#include "MBinning.h"
87#include "MParList.h"
88
89#include "MGeomCam.h"
90#include "MGeomPix.h"
91
92#include "MSignalCam.h"
93#include "MSignalPix.h"
94
95#include "MMuonSetup.h"
96#include "MMuonCalibPar.h"
97#include "MMuonSearchPar.h"
98
99ClassImp(MHSingleMuon);
100
101using namespace std;
102
103// --------------------------------------------------------------------------
104//
105// Setup histograms
106//
107MHSingleMuon::MHSingleMuon(const char *name, const char *title) :
108 fSignalCam(0), fMuonSearchPar(0), fGeomCam(0), fMargin(0)
109{
110 fName = name ? name : "MHSingleMuon";
111 fTitle = title ? title : "Histograms of muon parameters";
112
113 fHistPhi.SetName("HistPhi");
114 fHistPhi.SetTitle("HistPhi");
115 fHistPhi.SetXTitle("\\phi [\\circ]");
116 fHistPhi.SetYTitle("sum of ADC");
117 fHistPhi.SetDirectory(NULL);
118 fHistPhi.SetFillStyle(4000);
119 fHistPhi.UseCurrentStyle();
120
121 fHistWidth.SetName("HistWidth");
122 fHistWidth.SetTitle("HistWidth");
123 fHistWidth.SetXTitle("distance from the ring center [\\circ]");
124 fHistWidth.SetYTitle("sum of ADC");
125 fHistWidth.SetDirectory(NULL);
126 fHistWidth.SetFillStyle(4000);
127 fHistWidth.UseCurrentStyle();
128
129 MBinning bins;
130 bins.SetEdges(20, -180, 180);
131 bins.Apply(fHistPhi);
132
133 bins.SetEdges(28, 0.3, 1.7);
134 bins.Apply(fHistWidth);
135}
136
137// --------------------------------------------------------------------------
138//
139// Setup the Binning for the histograms automatically if the correct
140// instances of MBinning
141//
142Bool_t MHSingleMuon::SetupFill(const MParList *plist)
143{
144 fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
145 if (!fGeomCam)
146 {
147 *fLog << warn << "MGeomCam not found... abort." << endl;
148 return kFALSE;
149 }
150 fMuonSearchPar = (MMuonSearchPar*)plist->FindObject("MMuonSearchPar");
151 if (!fMuonSearchPar)
152 {
153 *fLog << warn << "MMuonSearchPar not found... abort." << endl;
154 return kFALSE;
155 }
156 fSignalCam = (MSignalCam*)plist->FindObject("MSignalCam");
157 if (!fSignalCam)
158 {
159 *fLog << warn << "MSignalCam not found... abort." << endl;
160 return kFALSE;
161 }
162
163 MMuonSetup *setup = (MMuonSetup*)const_cast<MParList*>(plist)->FindCreateObj("MMuonSetup");
164 if (!setup)
165 return kFALSE;
166
167 fMargin = setup->GetMargin()/fGeomCam->GetConvMm2Deg();
168
169 ApplyBinning(*plist, "ArcPhi", &fHistPhi);
170 ApplyBinning(*plist, "MuonWidth", &fHistWidth);
171
172 return kTRUE;
173}
174
175// --------------------------------------------------------------------------
176//
177// Fill the histograms with data from a MMuonCalibPar and
178// MMuonSearchPar container.
179//
180Bool_t MHSingleMuon::Fill(const MParContainer *par, const Stat_t w)
181{
182 fHistPhi.Reset();
183 fHistWidth.Reset();
184
185 const Int_t entries = fSignalCam->GetNumPixels();
186
187 // the position of the center of a muon ring
188 const Float_t cenx = fMuonSearchPar->GetCenterX();
189 const Float_t ceny = fMuonSearchPar->GetCenterY();
190
191 for (Int_t i=0; i<entries; i++)
192 {
193 const MSignalPix &pix = (*fSignalCam)[i];
194 const MGeomPix &gpix = (*fGeomCam)[i];
195
196 const Float_t dx = gpix.GetX() - cenx;
197 const Float_t dy = gpix.GetY() - ceny;
198
199 const Float_t dist = TMath::Hypot(dx, dy);
200
201 // if the signal is not near the estimated circle, it is ignored.
202 if (dist < fMuonSearchPar->GetRadius() + fMargin &&
203 dist > fMuonSearchPar->GetRadius() - fMargin)
204 fHistPhi.Fill(TMath::ATan2(dx, dy)*TMath::RadToDeg(), pix.GetNumPhotons());
205
206 // use only the inner pixles. FIXME: This is geometry dependent
207 if(i>397)
208 continue;
209
210 fHistWidth.Fill(dist*fGeomCam->GetConvMm2Deg(), pix.GetNumPhotons());
211 }
212
213 // error estimation (temporaly)
214 // The error is estimated from the signal. In order to do so, we have to
215 // once convert the signal from ADC to photo-electron. Then we can get
216 // the fluctuation such as F-factor*sqrt(phe).
217 // Up to now, the error of pedestal is not taken into accout. This is not
218 // of course correct. We will include this soon.
219 const Double_t Ffactor = 1.4;
220 for (Int_t i=0; i<fHistPhi.GetNbinsX()+1; i++)
221 {
222 const Float_t abs = TMath::Abs(fHistPhi.GetBinContent(i));
223 const Float_t rougherr = TMath::Sqrt(abs)*Ffactor;
224
225 fHistPhi.SetBinError(i, rougherr);
226 }
227
228 for (Int_t i=0; i<fHistWidth.GetNbinsX()+1; i++)
229 {
230 const Float_t abs = TMath::Abs(fHistWidth.GetBinContent(i));
231 const Float_t rougherr = TMath::Sqrt(abs)*Ffactor;
232
233 fHistWidth.SetBinError(i, rougherr);
234 }
235
236 return kTRUE;
237}
238
239// --------------------------------------------------------------------------
240//
241// Find the first bins starting at the bin with maximum content in both
242// directions which are below threshold.
243// If in a range of half the histogram size in both directions no bin
244// below the threshold is found, kFALSE is returned.
245//
246Bool_t MHSingleMuon::FindRangeAboveThreshold(const TH1 &h, Float_t thres, Int_t &first, Int_t &last) const
247{
248 const Int_t n = h.GetNbinsX();
249 const Int_t maxbin = h.GetMaximumBin();
250 const Int_t edge = maxbin+n/2;
251
252 // Search from the peak to the right
253 last = -1;
254 for (Int_t i=maxbin; i<edge; i++)
255 {
256 const Float_t val = h.GetBinContent(i%n + 1);
257 if (val<thres)
258 {
259 last = i%n+1;
260 break;
261 }
262 }
263
264 // Search from the peak to the left
265 first = -1;
266 for (Int_t i=maxbin+n-1; i>=edge; i--)
267 {
268 const Float_t val = h.GetBinContent(i%n + 1);
269 if (val<thres)
270 {
271 first = i%n+1;
272 break;
273 }
274 }
275
276 return first>=0 && last>=0;
277}
278
279// --------------------------------------------------------------------------
280//
281// Photon distribution along the estimated circle is fitted with theoritical
282// function in order to get some more information such as Arc Phi and Arc
283// Length.
284//
285Bool_t MHSingleMuon::CalcPhi(Double_t thres, Double_t &peakphi, Double_t &arcphi) const
286{
287 if (fHistPhi.GetMaximum()<thres)
288 return kFALSE;
289
290 peakphi = 180.-fHistPhi.GetBinCenter(fHistPhi.GetMaximumBin());
291
292 // Now find the position at which the peak edges crosses the threshold
293 Int_t first, last;
294
295 FindRangeAboveThreshold(fHistPhi, thres, first, last);
296
297 const Int_t n = fHistPhi.GetNbinsX();
298 const Int_t edge = fHistPhi.GetMaximumBin()+n/2;
299 if (first<0)
300 first = (edge-1)%n+1;
301 if (last<0)
302 last = edge%n+1;;
303
304 const Float_t startfitval = fHistPhi.GetBinLowEdge(first+1);
305 const Float_t endfitval = fHistPhi.GetBinLowEdge(last);
306
307 arcphi = last-1<first ? 360+(endfitval-startfitval) : endfitval-startfitval;
308
309 //if (fEnableImpactCalc)
310 // CalcImpact(effbinnum, startfitval, endfitval);
311
312 return kTRUE;
313}
314
315// --------------------------------------------------------------------------
316//
317// Photon distribution of distance from the center of estimated ring is
318// fitted in order to get some more information such as ARC WIDTH which
319// can represent to the PSF of our reflector.
320//
321// thres: Threshold above zero to determin the edges of the peak which
322// is used as fit range
323// width: ArcWidth returned in deg
324// chi: Chi^2/NDF of the fit
325//
326Bool_t MHSingleMuon::CalcWidth(Double_t thres, Double_t &width, Double_t &chi)
327{
328 Int_t first, last;
329
330 if (!FindRangeAboveThreshold(fHistWidth, thres, first, last))
331 return kFALSE;
332
333 // This happens in some cases
334 const Int_t n = fHistWidth.GetNbinsX()/2;
335 const Int_t m = fHistWidth.GetMaximumBin();
336 if (first>last)
337 if (m>n) // If maximum is on the right side of histogram
338 last = n;
339 else
340 first = 0; // If maximum is on the left side of histogram
341
342 if (last-first<=3)
343 return kFALSE;
344
345 // Now get the fit range
346 const Float_t startfitval = fHistWidth.GetBinLowEdge(first+1);
347 const Float_t endfitval = fHistWidth.GetBinLowEdge(last);
348
349 // Setup the function and perform the fit
350 TF1 f1("f1", "gaus + [3]", startfitval, endfitval);
351 f1.SetLineColor(kBlue);
352
353 // Choose starting values as accurate as possible
354 f1.SetParameter(0, fHistWidth.GetMaximum());
355 f1.SetParameter(1, fHistWidth.GetBinCenter(m));
356// f1.SetParameter(2, (endfitval-startfitval)/2);
357 f1.SetParameter(2, 0.1);
358 f1.SetParameter(3, 1.8);
359
360 // options : N do not store the function, do not draw
361 // I use integral of function in bin rather than value at bin center
362 // R use the range specified in the function range
363 // Q quiet mode
364 fHistWidth.Fit(&f1, "QR0");
365
366 chi = f1.GetChisquare()/f1.GetNDF();
367
368 Double_t err;
369 gMinuit->GetParameter(2, width, err); // get the sigma value
370
371 return kTRUE;
372}
373
374/*
375// --------------------------------------------------------------------------
376//
377// An impact parameter is calculated by fitting the histogram of photon
378// distribution along the circle with a theoritical model.
379// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.
380// The function (6) is used here.)
381//
382// By default this calculation is suppressed because this calculation is
383// very time consuming. If you want to calculate an impact parameter,
384// you can call the function of EnableImpactCalc().
385//
386void MMuonCalibParCalc::CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval)
387{
388 // Fit the distribution with Vacanti function. The function is different
389 // for the impact parameter of inside or outside of our reflector.
390 // Then two different functions are applied to the photon distribution,
391 // and the one which give us smaller chisquare value is taken as a
392 // proper one.
393
394 Double_t val1,err1,val2,err2;
395 // impact parameter inside mirror radius (8.5m)
396 TString func1;
397 Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
398 tmpval = sin(2.*tmpval)*8.5;
399 func1 += "[0]*";
400 func1 += tmpval;
401 func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
402
403 TF1 f1("f1",func1,startfitval,endfitval);
404 f1.SetParameters(2000,3,0);
405 f1.SetParLimits(1,0,8.5);
406 f1.SetParLimits(2,-180.,180.);
407
408 fMuonCalibPar->fHistPhi->Fit("f1","RQEM");
409
410 Float_t chi1 = -1;
411 Float_t chi2 = -1;
412 if(effbinnum>3)
413 chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
414
415 gMinuit->GetParameter(1,val1,err1); // get the estimated IP
416 Float_t estip1 = val1;
417
418 // impact parameter beyond mirror area (8.5m)
419 TString func2;
420 Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
421 tmpval2 = sin(2.*tmpval2)*8.5*2.;
422 func2 += "[0]*";
423 func2 += tmpval2;
424 func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
425
426 TF1 f2("f2",func2,startfitval,endfitval);
427 f2.SetParameters(2000,20,0);
428 f2.SetParLimits(1,8.5,300.);
429 f2.SetParLimits(2,-180.,180.);
430
431 fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");
432
433 if(effbinnum>3)
434 chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
435
436 gMinuit->GetParameter(1,val2,err2); // get the estimated IP
437 Float_t estip2 = val2;
438
439 if(effbinnum<=3)
440 {
441 fMuonCalibPar->SetEstImpact(-1.);
442 }
443 if(chi2 > chi1)
444 {
445 fMuonCalibPar->SetEstImpact(estip1);
446 fMuonCalibPar->SetChiArcPhi(chi1);
447 }
448 else
449 {
450 fMuonCalibPar->SetEstImpact(estip2);
451 fMuonCalibPar->SetChiArcPhi(chi2);
452 }
453}
454*/
455
456void MHSingleMuon::Paint(Option_t *o)
457{
458 TF1 *f = fHistWidth.GetFunction("f1");
459 if (f)
460 f->ResetBit(1<<9);
461}
462
463void MHSingleMuon::Draw(Option_t *o)
464{
465 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
466 pad->SetBorderMode(0);
467
468 AppendPad("");
469
470 pad->Divide(1,2);
471
472 pad->cd(1);
473 gPad->SetBorderMode(0);
474 fHistPhi.Draw();
475
476 pad->cd(2);
477 gPad->SetBorderMode(0);
478 fHistWidth.Draw();
479
480}
Note: See TracBrowser for help on using the repository browser.