source: trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc@ 5245

Last change on this file since 5245 was 5210, checked in by mase, 21 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 <mailto:mase@mppmu.mpg.de>
19! Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MMuonCalibPar
29//
30// Storage Container for muon
31//
32// This class holds some information for a calibraion using muons. Muons
33// are identified by using the class of the MMuonSearchParCalc. You can fill
34// these information by using the MMuonCalibParCalc. See also these class
35// manuals.
36//
37//
38// Input Containers:
39// [MGeomCam]
40// [MCerPhotEvt]
41// [MMuonSearchPar]
42//
43/////////////////////////////////////////////////////////////////////////////
44#include "MMuonCalibPar.h"
45
46#include <fstream>
47
48#include <TH1.h>
49#include <TF1.h>
50#include <TMinuit.h>
51
52#include "MLog.h"
53#include "MLogManip.h"
54#include "MGeomCam.h"
55#include "MGeomPix.h"
56#include "MCerPhotEvt.h"
57#include "MCerPhotPix.h"
58#include "MMuonSearchPar.h"
59#include "MBinning.h"
60
61using namespace std;
62
63ClassImp(MMuonCalibPar);
64
65// --------------------------------------------------------------------------
66//
67// Default constructor.
68//
69MMuonCalibPar::MMuonCalibPar(const char *name, const char *title)
70{
71 fName = name ? name : "MMuonCalibPar";
72 fTitle = title ? title : "Muon calibration parameters";
73
74 fHistPhi = new TH1F;
75 fHistWidth = new TH1F;
76
77 fHistPhi->SetName("HistPhi");
78 fHistPhi->SetTitle("HistPhi");
79 fHistPhi->SetXTitle("phi [deg.]");
80 fHistPhi->SetYTitle("sum of ADC");
81 fHistPhi->SetDirectory(NULL);
82 fHistPhi->SetFillStyle(4000);
83 fHistPhi->UseCurrentStyle();
84
85 fHistWidth->SetName("HistWidth");
86 fHistWidth->SetTitle("HistWidth");
87 fHistWidth->SetXTitle("distance from the ring center [deg.]");
88 fHistWidth->SetYTitle("sum of ADC");
89 fHistWidth->SetDirectory(NULL);
90 fHistWidth->SetFillStyle(4000);
91 fHistWidth->UseCurrentStyle();
92
93 fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
94 fDisablePreCuts = kFALSE; // By default the pre cuts will be applied.
95 fUseCleanForWidth = kFALSE; // By default all the pixels will be used for the histogram of arc width.
96
97 fMargin = 60.; // in mm
98 fArcPhiThres = 100.;
99 fArcWidthThres = 100.;
100 fArcPhiBinNum = 20;
101 fArcPhiHistStartVal = -180.; // deg.
102 fArcPhiHistEndVal = 180.; // deg.
103 fArcWidthBinNum = 28;
104 fArcWidthHistStartVal = 0.3; // deg.
105 fArcWidthHistEndVal = 1.7; // deg.
106}
107
108// --------------------------------------------------------------------------
109//
110MMuonCalibPar::~MMuonCalibPar()
111{
112 delete fHistPhi;
113 delete fHistWidth;
114}
115
116// --------------------------------------------------------------------------
117//
118void MMuonCalibPar::Reset()
119{
120 fArcLength = -1.;
121 fArcPhi = 0.;
122 fArcWidth = -1.;
123 fChiArcPhi = -1.;
124 fChiArcWidth = -1.;
125 fMuonSize = 0.;
126 fEstImpact = -1.;
127 fUseUnmap = kFALSE;
128 fPeakPhi = 0.;
129}
130
131// --------------------------------------------------------------------------
132//
133// This function fill the histograms in order to get muon parameters.
134// For the evaluation of the Arc Width, we use only the signals in the inner
135// part. You can use the image after the cleaning by using the function of
136// UseCleanForWidth(). See the manual of MMuonCalibParCalc.
137//
138void MMuonCalibPar::FillHist
139( const MGeomCam &geom, const MCerPhotEvt &evt,
140 const MMuonSearchPar &musearch)
141{
142 // preparation for a histgram
143 MBinning binsphi;
144 binsphi.SetEdges(fArcPhiBinNum, fArcPhiHistStartVal, fArcPhiHistEndVal);
145 binsphi.Apply(*fHistPhi);
146
147 MBinning binswid;
148 binswid.SetEdges(fArcWidthBinNum, fArcWidthHistStartVal, fArcWidthHistEndVal);
149 binswid.Apply(*fHistWidth);
150
151 const Int_t entries = evt.GetNumPixels();
152
153 // the position of the center of a muon ring
154 const Float_t cenx = musearch.GetCenterX();
155 const Float_t ceny = musearch.GetCenterY();
156
157 for (Int_t i=0; i<entries; i++ )
158 {
159 MCerPhotPix &pix = (evt)[i];
160
161 const MGeomPix &gpix = (geom)[pix.GetPixId()];
162
163 const Float_t dx = gpix.GetX() - cenx;
164 const Float_t dy = gpix.GetY() - ceny;
165
166 const Float_t dist = TMath::Sqrt(dx*dx+dy*dy);
167
168 Float_t ang = TMath::ACos(dx/dist);
169 if(dy>0)
170 ang *= -1.0;
171
172 // if the signal is not near the estimated circle, it is ignored.
173 if(dist < musearch.GetRadius() + fMargin && dist > musearch.GetRadius() - fMargin)
174 {
175 // check whether ummapped pixel is used or not.
176 // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information.
177 if(pix.IsPixelUnmapped())
178 {
179 fUseUnmap = kTRUE;
180 continue;
181 }
182 fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons());
183 fMuonSize += pix.GetNumPhotons();
184 }
185
186 if(pix.GetPixId()>397)
187 continue; // use only the inner pixles
188
189 if(fUseCleanForWidth)
190 {
191 if(!pix.IsPixelUsed())
192 continue;
193 }
194
195 fHistWidth->Fill(dist*geom.GetConvMm2Deg(), pix.GetNumPhotons());
196 }
197
198
199 // error estimation (temporaly)
200 // The error is estimated from the signal. In order to do so, we have to
201 // once convert the signal from ADC to photo-electron. Then we can get
202 // the fluctuation such as F-factor*sqrt(phe).
203 // Up to now, the error of pedestal is not taken into accout. This is not
204 // of course correct. We will include this soon.
205 Double_t ADC2PhEl = 0.14;
206 Double_t Ffactor = 1.4;
207 for(Int_t i=0; i<fArcPhiBinNum+1; i++)
208 {
209 Float_t rougherr = TMath::Sqrt(TMath::Abs(fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
210 {
211 fHistPhi->SetBinError(i, rougherr);
212 }
213 }
214 for(Int_t i=0; i<fArcWidthBinNum+1; i++)
215 {
216 Float_t rougherr = TMath::Sqrt(TMath::Abs(fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
217 {
218 fHistWidth->SetBinError(i, rougherr);
219 }
220 }
221}
222
223// --------------------------------------------------------------------------
224//
225// Photon distribution along the estimated circle is fitted with theoritical
226// function in order to get some more information such as Arc Phi and Arc Length.
227//
228void MMuonCalibPar::CalcPhi
229(const MGeomCam &geom, const MCerPhotEvt &evt,
230 const MMuonSearchPar &musearch)
231{
232 Float_t convbin2val = (fArcPhiHistEndVal-fArcPhiHistStartVal)/
233 (Float_t)fArcPhiBinNum;
234
235 // adjust the peak to 0
236 Float_t maxval = 0.;
237 Int_t maxbin = 0;
238 maxval = fHistPhi->GetMaximum();
239 maxbin = fHistPhi->GetMaximumBin();
240 fPeakPhi = 180.-(Float_t)(maxbin-1)*convbin2val;
241 TArrayD tmp;
242 tmp.Set(fArcPhiBinNum+1);
243 for(Int_t i=1; i<fArcPhiBinNum+1; i++)
244 {
245 tmp[i] = fHistPhi->GetBinContent(i);
246 }
247 for(Int_t i=1; i<fArcPhiBinNum+1; i++)
248 {
249 Int_t id;
250 id = i + (maxbin-(Int_t)((Float_t)fArcPhiBinNum/2.)-1);
251 if(id>fArcPhiBinNum)
252 {
253 id-=(fArcPhiBinNum);
254 }
255 if(id<=0)
256 {
257 id+=(fArcPhiBinNum);
258 }
259 fHistPhi->SetBinContent(i,tmp[id]);
260 }
261 maxbin = (Int_t)((Float_t)fArcPhiBinNum/2.)+1;
262
263 // Determination of fitting region
264 // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore,
265 // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!!
266 Float_t startfitval = 0.;
267 Float_t endfitval = 0.;
268 Bool_t IsInMaxim = kFALSE;
269 Int_t effbinnum = 0;
270 for(Int_t i=1; i<fArcPhiBinNum+1; i++)
271 {
272 Float_t content = fHistPhi->GetBinContent(i);
273 Float_t content_pre = fHistPhi->GetBinContent(i-1);
274
275 if(content > fArcPhiThres && content_pre < fArcPhiThres)
276 {
277 startfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal;
278 }
279 if(i==maxbin)
280 IsInMaxim = kTRUE;
281
282 if(content < fArcPhiThres && IsInMaxim == kTRUE)
283 {
284 endfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal;
285 break;
286 }
287 endfitval = fArcPhiHistEndVal;
288 }
289
290 effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
291
292 fArcPhi = effbinnum*convbin2val;
293 fArcLength = fArcPhi*3.1415926/180.*musearch.GetRadius()*geom.GetConvMm2Deg();
294
295 if(fEnableImpactCalc)
296 CalcImpact(geom, musearch, effbinnum, startfitval, endfitval);
297}
298
299// --------------------------------------------------------------------------
300//
301// An impact parameter is calculated by fitting the histogram of photon
302// distribution along the circle with a theoritical model.
303// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.
304// The function (6) is used here.)
305//
306// By default this calculation is suppressed because this calculation is
307// very time consuming. If you want to calculate an impact parameter,
308// you can call the function of EnableImpactCalc().
309//
310void MMuonCalibPar::CalcImpact
311( const MGeomCam &geom, const MMuonSearchPar &musearch,
312 Int_t effbinnum, Float_t startfitval, Float_t endfitval)
313{
314 // Fit the distribution with Vacanti function. The function is different
315 // for the impact parameter of inside or outside of our reflector.
316 // Then two different functions are applied to the photon distribution,
317 // and the one which give us smaller chisquare value is taken as a
318 // proper one.
319 Double_t val1,err1,val2,err2;
320 // impact parameter inside mirror radius (8.5m)
321 TString func1;
322 Float_t tmpval = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.;
323 tmpval = sin(2.*tmpval)*8.5;
324 func1 += "[0]*";
325 func1 += tmpval;
326 func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
327
328 TF1 f1("f1",func1,startfitval,endfitval);
329 f1.SetParameters(2000,3,0);
330 f1.SetParLimits(1,0,8.5);
331 f1.SetParLimits(2,-180.,180.);
332
333 fHistPhi->Fit("f1","RQEM");
334
335 Float_t chi1 = -1;
336 Float_t chi2 = -1;
337 if(effbinnum>3)
338 chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
339
340 gMinuit->GetParameter(1,val1,err1); // get the estimated IP
341 Float_t estip1 = val1;
342
343 // impact parameter beyond mirror area (8.5m)
344 TString func2;
345 Float_t tmpval2 = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.;
346 tmpval2 = sin(2.*tmpval2)*8.5*2.;
347 func2 += "[0]*";
348 func2 += tmpval2;
349 func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
350
351 TF1 f2("f2",func2,startfitval,endfitval);
352 f2.SetParameters(2000,20,0);
353 f2.SetParLimits(1,8.5,300.);
354 f2.SetParLimits(2,-180.,180.);
355
356 fHistPhi->Fit("f2","RQEM+");
357
358 if(effbinnum>3)
359 chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
360
361 gMinuit->GetParameter(1,val2,err2); // get the estimated IP
362 Float_t estip2 = val2;
363
364 if(effbinnum<=3)
365 {
366 fEstImpact = -1.;
367 }
368 if(chi2 > chi1)
369 {
370 fEstImpact = estip1;
371 fChiArcPhi = chi1;
372 }
373 else
374 {
375 fEstImpact = estip2;
376 fChiArcPhi = chi2;
377 }
378}
379
380// --------------------------------------------------------------------------
381//
382// Photon distribution of distance from the center of estimated ring is
383// fitted in order to get some more information such as ARC WIDTH which
384// can represent to the PSF of our reflector.
385//
386Float_t MMuonCalibPar::CalcWidth
387(const MGeomCam &geom, const MCerPhotEvt &evt,
388 const MMuonSearchPar &musearch)
389{
390 Float_t convbin2val = (fArcWidthHistEndVal - fArcWidthHistStartVal)
391 /(Float_t)fArcWidthBinNum;
392
393 // determination of fitting region
394 Int_t maxbin = fHistWidth->GetMaximumBin();
395 Float_t startfitval = 0.;
396 Float_t endfitval = 0.;
397 Bool_t IsInMaxim = kFALSE;
398 Int_t effbinnum = 0;
399 for(Int_t i=1; i<fArcWidthBinNum+1; i++)
400 {
401 Float_t content = fHistWidth->GetBinContent(i);
402 Float_t content_pre = fHistWidth->GetBinContent(i-1);
403
404 if(content > fArcWidthThres)
405 effbinnum++;
406
407 if(content > fArcWidthThres && content_pre < fArcWidthThres)
408 {
409 startfitval = (Float_t)(i-4)*convbin2val + fArcWidthHistStartVal;
410 if(startfitval<0.) startfitval = 0.;
411 }
412 if(i==maxbin)
413 IsInMaxim = kTRUE;
414
415 if(content < fArcWidthThres && IsInMaxim == kTRUE)
416 {
417 endfitval = (Float_t)(i+2)*convbin2val + fArcWidthHistStartVal;
418 if(endfitval>180.) endfitval = 180.;
419 break;
420 }
421 endfitval = fArcWidthHistEndVal;
422 }
423 effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
424
425 TF1 f1("f1","gaus",startfitval,endfitval);
426
427 fHistWidth->Fit("f1","QR","",startfitval,endfitval);
428
429 if(effbinnum>3)
430 fChiArcWidth = f1.GetChisquare()/((Float_t)(effbinnum-3));
431
432 Double_t val,err;
433 gMinuit->GetParameter(2,val,err); // get the sigma value
434
435 return val;
436}
437
438// --------------------------------------------------------------------------
439//
440// Calculation of muon parameters
441//
442Int_t MMuonCalibPar::Calc
443(const MGeomCam &geom, const MCerPhotEvt &evt,
444 MMuonSearchPar &musearch, const Float_t *cuts)
445{
446 // sanity check
447 if(evt.GetNumPixels() < 3)
448 return kCONTINUE;
449
450 // If an event does not seem to be like muon, the calculation will be skipped.
451 if(musearch.IsNoMuon())
452 return kCONTINUE;
453
454 // Pre Cuts 1
455 if(!fDisablePreCuts)
456 {
457 if(musearch.GetRadius() < cuts[0] || musearch.GetRadius() > cuts[1])
458 {
459 musearch.SetNoMuon();
460 return kCONTINUE;
461 }
462 if(musearch.GetDeviation() > cuts[2])
463 {
464 musearch.SetNoMuon();
465 return kCONTINUE;
466 }
467 }
468
469 Reset();
470
471 FillHist(geom,evt,musearch);
472
473 CalcPhi(geom,evt,musearch);
474
475 // Pre Cuts 2
476 if(!fDisablePreCuts)
477 {
478 if(fMuonSize < cuts[3] || fArcPhi < cuts[4])
479 {
480 musearch.SetNoMuon();
481 return kCONTINUE;
482 }
483 }
484
485 fArcWidth = CalcWidth(geom,evt,musearch);
486
487 SetReadyToSave();
488
489 return kTRUE;
490}
491
492void MMuonCalibPar::Print(Option_t *) const
493{
494 *fLog << all;
495 *fLog << "Muon Parameters (" << GetName() << ")" << endl;
496 *fLog << " - Arc Length [deg.] = " << fArcLength << endl;
497 *fLog << " - Arc Phi [deg.] = " << fArcPhi << endl;
498 *fLog << " - Arc Width [deg.] = " << fArcWidth << endl;
499 *fLog << " - Chi Arc Phi [x2/ndf]= " << fChiArcPhi << endl;
500 *fLog << " - Chi Arc Width [x2/ndf]= " << fChiArcWidth << endl;
501 *fLog << " - Est. I. P. [m] = " << fEstImpact << endl;
502 *fLog << " - Size of muon = " << fMuonSize << endl;
503 *fLog << " - Peak Phi [deg.] = " << fPeakPhi << endl;
504 *fLog << " - UseUnmap = " << fUseUnmap << endl;
505}
506
507
508
509
Note: See TracBrowser for help on using the repository browser.