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

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