source: branches/MarsMoreSimulationTruth/mmuon/MHSingleMuon.cc@ 19004

Last change on this file since 19004 was 14911, checked in by tbretz, 12 years ago
Reverting to last revision.
File size: 17.2 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// Class Version 2:
76// ----------------
77// + Double_t fRelTimeMean; // Result of the gaus fit to the arrival time
78// + Double_t fRelTimeSigma; // Result of the gaus fit to the arrival time
79//
80////////////////////////////////////////////////////////////////////////////
81#include "MHSingleMuon.h"
82
83#include <TF1.h>
84#include <TPad.h>
85#include <TCanvas.h>
86
87#include "MLog.h"
88#include "MLogManip.h"
89
90#include "MBinning.h"
91#include "MParList.h"
92
93#include "MGeomCam.h"
94#include "MGeomPix.h"
95
96#include "MSignalCam.h"
97#include "MSignalPix.h"
98
99#include "MMuonSetup.h"
100#include "MMuonCalibPar.h"
101#include "MMuonSearchPar.h"
102
103ClassImp(MHSingleMuon);
104
105using namespace std;
106
107// --------------------------------------------------------------------------
108//
109// Setup histograms
110//
111MHSingleMuon::MHSingleMuon(const char *name, const char *title) :
112 fSignalCam(0), fMuonSearchPar(0), fGeomCam(0), fMargin(0)
113{
114 fName = name ? name : "MHSingleMuon";
115 fTitle = title ? title : "Histograms of muon parameters";
116
117 fHistPhi.SetName("HistPhi");
118 fHistPhi.SetTitle("HistPhi");
119 fHistPhi.SetXTitle("\\phi [\\circ]");
120 fHistPhi.SetYTitle("sum of ADC");
121 fHistPhi.SetDirectory(NULL);
122 fHistPhi.SetFillStyle(4000);
123 fHistPhi.UseCurrentStyle();
124
125 fHistWidth.SetName("HistWidth");
126 fHistWidth.SetTitle("HistWidth");
127 fHistWidth.SetXTitle("distance from the ring center [\\circ]");
128 fHistWidth.SetYTitle("sum of ADC");
129 fHistWidth.SetDirectory(NULL);
130 fHistWidth.SetFillStyle(4000);
131 fHistWidth.UseCurrentStyle();
132
133 fHistTime.SetName("HistTime");
134 fHistTime.SetTitle("HistTime");
135 fHistTime.SetXTitle("timing difference");
136 fHistTime.SetYTitle("Counts");
137 fHistTime.SetDirectory(NULL);
138 fHistTime.SetFillStyle(4000);
139 fHistTime.UseCurrentStyle();
140
141 MBinning bins;
142 bins.SetEdges(20, -180, 180);
143 bins.Apply(fHistPhi);
144
145 bins.SetEdges(28, 0.3, 1.7);
146 bins.Apply(fHistWidth);
147
148 bins.SetEdges(101, -33, 33); // +/- 33ns
149 bins.Apply(fHistTime);
150}
151
152// --------------------------------------------------------------------------
153//
154// Setup the Binning for the histograms automatically if the correct
155// instances of MBinning
156//
157Bool_t MHSingleMuon::SetupFill(const MParList *plist)
158{
159 fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
160 if (!fGeomCam)
161 {
162 *fLog << warn << "MGeomCam not found... abort." << endl;
163 return kFALSE;
164 }
165 fMuonSearchPar = (MMuonSearchPar*)plist->FindObject("MMuonSearchPar");
166 if (!fMuonSearchPar)
167 {
168 *fLog << warn << "MMuonSearchPar not found... abort." << endl;
169 return kFALSE;
170 }
171 fSignalCam = (MSignalCam*)plist->FindObject("MSignalCam");
172 if (!fSignalCam)
173 {
174 *fLog << warn << "MSignalCam not found... abort." << endl;
175 return kFALSE;
176 }
177
178 MMuonSetup *setup = (MMuonSetup*)const_cast<MParList*>(plist)->FindCreateObj("MMuonSetup");
179 if (!setup)
180 return kFALSE;
181
182 fMargin = setup->GetMargin()/fGeomCam->GetConvMm2Deg();
183
184 ApplyBinning(*plist, "ArcPhi", fHistPhi);
185 ApplyBinning(*plist, "MuonWidth", fHistWidth);
186 ApplyBinning(*plist, "MuonTime", fHistTime);
187
188 return kTRUE;
189}
190
191// --------------------------------------------------------------------------
192//
193// Fill the histograms with data from a MMuonCalibPar and
194// MMuonSearchPar container.
195//
196Int_t MHSingleMuon::Fill(const MParContainer *par, const Stat_t w)
197{
198 fRelTimeMean = 0;
199 fRelTimeSigma = -1;
200
201 fHistPhi.Reset();
202 fHistWidth.Reset();
203 fHistTime.Reset();
204
205 const Int_t entries = fSignalCam->GetNumPixels();
206
207 // the position of the center of a muon ring
208 const Float_t cenx = fMuonSearchPar->GetCenterX();
209 const Float_t ceny = fMuonSearchPar->GetCenterY();
210
211 for (Int_t i=0; i<entries; i++)
212 {
213 const MSignalPix &pix = (*fSignalCam)[i];
214 const MGeom &gpix = (*fGeomCam)[i];
215
216 const Float_t dx = gpix.GetX() - cenx;
217 const Float_t dy = gpix.GetY() - ceny;
218
219 const Float_t dist = TMath::Hypot(dx, dy);
220
221 // if the signal is not near the estimated circle, it is ignored.
222 if (TMath::Abs(dist-fMuonSearchPar->GetRadius())<fMargin)
223 {
224 // The arrival time is aligned around 0 for smaller
225 // and more stable histogram range
226 fHistTime.Fill(pix.GetArrivalTime()-fMuonSearchPar->GetTime());
227 }
228
229 // use only the inner pixles. FIXME: This is geometry dependent
230 if (gpix.GetAidx()>0)
231 continue;
232
233 fHistWidth.Fill(dist*fGeomCam->GetConvMm2Deg(), pix.GetNumPhotons());
234 }
235
236 // Setup the function and perform the fit
237 TF1 g1("g1", "gaus");//, -fHistTime.GetXmin(), fHistTime.GetXmax());
238
239 // Choose starting values as accurate as possible
240 g1.SetParameter(0, fHistTime.GetMaximum());
241 g1.SetParameter(1, 0);
242 g1.SetParameter(2, 0.7); // FIXME! GetRMS instead???
243
244 // According to fMuonSearchPar->GetTimeRMS() identified muons
245 // do not have an arrival time rms>3
246 g1.SetParLimits(1, -1.7, 1.7);
247 g1.SetParLimits(2, 0, 3.4);
248
249 // options : N do not store the function, do not draw
250 // I use integral of function in bin rather than value at bin center
251 // R use the range specified in the function range
252 // Q quiet mode
253 if (fHistTime.Fit(&g1, "QNB"))
254 return kTRUE;
255
256 fRelTimeMean = g1.GetParameter(1);
257 fRelTimeSigma = g1.GetParameter(2);
258
259 // The mean arrival time which was subtracted before will
260 // be added again, now
261 const Double_t tm0 = fMuonSearchPar->GetTime()+fRelTimeMean;
262
263 for (Int_t i=0; i<entries; i++)
264 {
265 const MSignalPix &pix = (*fSignalCam)[i];
266 const MGeom &gpix = (*fGeomCam)[i];
267
268 const Float_t dx = gpix.GetX() - cenx;
269 const Float_t dy = gpix.GetY() - ceny;
270
271 const Float_t dist = TMath::Hypot(dx, dy);
272
273 // if the signal is not near the estimated circle, it is ignored.
274 if (TMath::Abs(dist-fMuonSearchPar->GetRadius())<fMargin &&
275 TMath::Abs(pix.GetArrivalTime()-tm0) < 2*fRelTimeSigma)
276 {
277 fHistPhi.Fill(TMath::ATan2(dx, dy)*TMath::RadToDeg(), pix.GetNumPhotons());
278 }
279 }
280
281 return kTRUE;
282
283/*
284 // Because the errors (sqrt(content)) are only scaled by a fixed
285 // factor, and the absolute value of the error is nowhere
286 // needed we skip this step
287
288 // error estimation (temporarily)
289 // The error is estimated from the signal. In order to do so, we have to
290 // once convert the signal from ADC to photo-electron. Then we can get
291 // the fluctuation such as F-factor*sqrt(phe).
292 // Up to now, the error of pedestal is not taken into accout. This is not
293 // of course correct. We will include this soon.
294 const Double_t Ffactor = 1.4;
295 for (Int_t i=0; i<fHistPhi.GetNbinsX()+1; i++)
296 fHistPhi.SetBinError(i, fHistPhi.GetBinError(i)*Ffactor);
297
298 for (Int_t i=0; i<fHistWidth.GetNbinsX()+1; i++)
299 fHistWidth.SetBinError(i, fHistWidth.GetBinError(i)*Ffactor);
300
301 return kTRUE;
302 */
303}
304
305// --------------------------------------------------------------------------
306//
307// Find the first bins starting at the bin with maximum content in both
308// directions which are below threshold.
309// If in a range of half the histogram size in both directions no bin
310// below the threshold is found, kFALSE is returned.
311//
312Bool_t MHSingleMuon::FindRangeAboveThreshold(const TProfile &h, Float_t thres, Int_t &first, Int_t &last) const
313{
314 const Int_t n = h.GetNbinsX();
315 const Int_t maxbin = h.GetMaximumBin();
316 const Int_t edge = maxbin+n/2;
317
318 // Search from the peak to the right
319 last = -1;
320 for (Int_t i=maxbin; i<edge; i++)
321 {
322 const Float_t val = h.GetBinContent(i%n + 1);
323 if (val<thres)
324 {
325 last = i%n+1;
326 break;
327 }
328 }
329
330 // Search from the peak to the left
331 first = -1;
332 for (Int_t i=maxbin+n-1; i>=edge; i--)
333 {
334 const Float_t val = h.GetBinContent(i%n + 1);
335 if (val<thres)
336 {
337 first = i%n+1;
338 break;
339 }
340 }
341
342 return first>=0 && last>=0;
343}
344
345// --------------------------------------------------------------------------
346//
347// Photon distribution along the estimated circle is fitted with theoritical
348// function in order to get some more information such as Arc Phi and Arc
349// Length.
350//
351Bool_t MHSingleMuon::CalcPhi(Double_t thres, Double_t &peakphi, Double_t &arcphi) const
352{
353 if (fHistPhi.GetMaximum()<thres)
354 return kFALSE;
355
356 peakphi = 180.-fHistPhi.GetBinCenter(fHistPhi.GetMaximumBin());
357
358 // Now find the position at which the peak edges crosses the threshold
359 Int_t first, last;
360
361 FindRangeAboveThreshold(fHistPhi, thres, first, last);
362
363 const Int_t n = fHistPhi.GetNbinsX();
364 const Int_t edge = fHistPhi.GetMaximumBin()+n/2;
365 if (first<0)
366 first = (edge-1)%n+1;
367 if (last<0)
368 last = edge%n+1;;
369
370 const Float_t startfitval = fHistPhi.GetBinLowEdge(first+1);
371 const Float_t endfitval = fHistPhi.GetBinLowEdge(last);
372
373 arcphi = last-1<first ? 360+(endfitval-startfitval) : endfitval-startfitval;
374
375 //if (fEnableImpactCalc)
376 // CalcImpact(effbinnum, startfitval, endfitval);
377
378 return kTRUE;
379}
380
381// --------------------------------------------------------------------------
382//
383// Photon distribution of distance from the center of estimated ring is
384// fitted in order to get some more information such as ARC WIDTH which
385// can represent to the PSF of our reflector.
386//
387// thres: Threshold above zero to determin the edges of the peak which
388// is used as fit range
389// width: ArcWidth returned in deg
390// chi: Chi^2/NDF of the fit
391//
392Bool_t MHSingleMuon::CalcWidth(Double_t thres, Double_t &width, Double_t &chi)
393{
394 Int_t first, last;
395
396 if (!FindRangeAboveThreshold(fHistWidth, thres, first, last))
397 return kFALSE;
398
399 // This happens in some cases
400 const Int_t n = fHistWidth.GetNbinsX()/2;
401 const Int_t m = fHistWidth.GetMaximumBin();
402 if (first>last)
403 {
404 if (m>n) // If maximum is on the right side of histogram
405 last = n;
406 else
407 first = 0; // If maximum is on the left side of histogram
408 }
409
410 if (last-first<=3)
411 return kFALSE;
412
413 // Now get the fit range
414 const Float_t startfitval = fHistWidth.GetBinLowEdge(first+1);
415 const Float_t endfitval = fHistWidth.GetBinLowEdge(last);
416
417 // Setup the function and perform the fit
418 TF1 f1("f1", "gaus + [3]", startfitval, endfitval);
419 f1.SetLineColor(kBlue);
420
421 // Choose starting values as accurate as possible
422 f1.SetParameter(0, fHistWidth.GetMaximum());
423 f1.SetParameter(1, fHistWidth.GetBinCenter(m));
424// f1.SetParameter(2, (endfitval-startfitval)/2);
425 f1.SetParameter(2, 0.1);
426 f1.SetParameter(3, 1.8);
427
428 // options : N do not store the function, do not draw
429 // I use integral of function in bin rather than value at bin center
430 // R use the range specified in the function range
431 // Q quiet mode
432// fHistWidth.Fit(&f1, "QRO");
433 if (fHistWidth.Fit(&f1, "QRN"))
434 return kFALSE;
435
436 chi = f1.GetChisquare()/f1.GetNDF();
437 width = f1.GetParameter(2);
438
439 return kTRUE;
440}
441
442/*
443// --------------------------------------------------------------------------
444//
445// An impact parameter is calculated by fitting the histogram of photon
446// distribution along the circle with a theoritical model.
447// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.
448// The function (6) is used here.)
449//
450// By default this calculation is suppressed because this calculation is
451// very time consuming. If you want to calculate an impact parameter,
452// you can call the function of EnableImpactCalc().
453//
454void MMuonCalibParCalc::CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval)
455{
456 // Fit the distribution with Vacanti function. The function is different
457 // for the impact parameter of inside or outside of our reflector.
458 // Then two different functions are applied to the photon distribution,
459 // and the one which give us smaller chisquare value is taken as a
460 // proper one.
461
462 Double_t val1,err1,val2,err2;
463 // impact parameter inside mirror radius (8.5m)
464 TString func1;
465 Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
466 tmpval = sin(2.*tmpval)*8.5;
467 func1 += "[0]*";
468 func1 += tmpval;
469 func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
470
471 TF1 f1("f1",func1,startfitval,endfitval);
472 f1.SetParameters(2000,3,0);
473 f1.SetParLimits(1,0,8.5);
474 f1.SetParLimits(2,-180.,180.);
475
476 fMuonCalibPar->fHistPhi->Fit("f1","RQEM");
477
478 Float_t chi1 = -1;
479 Float_t chi2 = -1;
480 if(effbinnum>3)
481 chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
482
483 gMinuit->GetParameter(1,val1,err1); // get the estimated IP
484 Float_t estip1 = val1;
485
486 // impact parameter beyond mirror area (8.5m)
487 TString func2;
488 Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
489 tmpval2 = sin(2.*tmpval2)*8.5*2.;
490 func2 += "[0]*";
491 func2 += tmpval2;
492 func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
493
494 TF1 f2("f2",func2,startfitval,endfitval);
495 f2.SetParameters(2000,20,0);
496 f2.SetParLimits(1,8.5,300.);
497 f2.SetParLimits(2,-180.,180.);
498
499 fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");
500
501 if(effbinnum>3)
502 chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
503
504 gMinuit->GetParameter(1,val2,err2); // get the estimated IP
505 Float_t estip2 = val2;
506
507 if(effbinnum<=3)
508 {
509 fMuonCalibPar->SetEstImpact(-1.);
510 }
511 if(chi2 > chi1)
512 {
513 fMuonCalibPar->SetEstImpact(estip1);
514 fMuonCalibPar->SetChiArcPhi(chi1);
515 }
516 else
517 {
518 fMuonCalibPar->SetEstImpact(estip2);
519 fMuonCalibPar->SetChiArcPhi(chi2);
520 }
521}
522*/
523
524Float_t MHSingleMuon::CalcSize() const
525{
526 const Int_t n = fHistPhi.GetNbinsX();
527
528 Double_t sz=0;
529 for (Int_t i=1; i<=n; i++)
530 sz += fHistPhi.GetBinContent(i)*fHistPhi.GetBinEntries(i);
531
532 return sz;
533}
534
535void MHSingleMuon::Paint(Option_t *o)
536{
537 TF1 *f = fHistWidth.GetFunction("f1");
538 if (f)
539 f->ResetBit(1<<9);
540}
541
542void MHSingleMuon::Draw(Option_t *o)
543{
544 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
545 pad->SetBorderMode(0);
546
547 AppendPad("");
548
549 pad->Divide(1,2);
550
551 pad->cd(1);
552 gPad->SetBorderMode(0);
553 fHistPhi.Draw();
554
555 pad->cd(2);
556 gPad->SetBorderMode(0);
557 fHistWidth.Draw();
558}
Note: See TracBrowser for help on using the repository browser.