source: trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc@ 5892

Last change on this file since 5892 was 5333, checked in by mase, 21 years ago
*** empty log message ***
File size: 17.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// MMuonCalibParCalc
29//
30// Task to calculate the muon parameters
31//
32// This class allows you to get more muon information especially useful for
33// the calibration of our telescope. This class store the information into the
34// container of MMuonCalibPar.
35//
36// In order to make this class work, we need the information of the arc
37// center and the radius. Therefore, we need to use the task of
38// MMuonSearchParCalc.
39//
40// You can use this class such as the followings;
41//
42// MTaskList tlist;
43// MMuonSearchParCalc musearchcalc;
44// MMuonCalibParCalc mucalibcalc;
45// tlist.AddToList(&musearchcalc);
46// tlist.AddToList(&mucalibcalc);.
47//
48// You may change the allowed region to estimate muon parameters such as
49// Muon SIZE and ARC LENGTH. The default value is 60 mm (0.2 deg.). If the
50// estimated radius of the arc is 1.0 degree, we take the photons in the
51// radius range from 0.8 to 1.2 degrees. You can change this value such as
52// the followings;
53//
54// mucalibcalc.SetMargin(60.);
55//
56// You can retrieve the histogram (TH1F) using the function of GetHistPhi()
57// (also GetHistWid()). Therefore, you can draw the histogram such as
58//
59// MParList plist;
60// MMuonCalibPar muparcalib;
61// plist.AddToList(&muparcalib);.
62// muparcalib.GetHistPhi().Draw();.
63//
64// In order to use another information of muons such as the center position
65// of the estimated circle, the radius of the circle. Use the infomation
66// stored in MMuonSearchPar.
67//
68//
69// For the faster computation, by default, the calculation of impact
70// parameter is suppressed. If you want to calculate the impact parameter
71// from the muon image, you can use the function of EnableImpactCalc(),
72// namely;
73//
74// mucalibcalc.EnableImpactCalc();.
75//
76// In addition, for the faster comutation, pre cuts to select the candidates
77// of muons for the calibration is done. You can set the values using the
78// function of SetPreCuts. This function takes 5 variables. They correspond
79// to the cur for the Arc Radius (low and high), the deviation of the fit
80// (high), the Muon Size (low) and Arc Phi (low). You can set them such as
81//
82// mucalibcalc.SetPreCuts(180., 400., 50., 2000., 180.);
83//
84// If you want to disable the pre cuts, you can disable it by using the
85// function of DisablePreCuts(), namely;
86//
87// mucalibcalc.DisablePreCuts();.
88//
89//
90// ### TODO ###
91// Up to now, in the histogram the error of the signal is estimated from
92// the signal using a rough conversion factor and a F-factor and this values
93// are global for all pixels. This is not the case for the real data. This
94// value should be taken from some containers. In addition, the error of
95// the pedestal is not taken into accout. The error treatment should be
96// done correctly.
97//
98//
99// Input Containers:
100// [MGeomCam]
101// [MCerPhotEvt]
102// [MMuonSearchPar]
103//
104// Output Containers:
105// [MMuonCalibPar]
106//
107//////////////////////////////////////////////////////////////////////////////
108
109#include "MMuonCalibParCalc.h"
110
111#include <fstream>
112
113#include <TH1.h>
114#include <TF1.h>
115#include <TMinuit.h>
116
117#include "MParList.h"
118
119#include "MGeomCam.h"
120#include "MGeomPix.h"
121#include "MSrcPosCam.h"
122#include "MCerPhotEvt.h"
123#include "MMuonSearchPar.h"
124#include "MMuonCalibPar.h"
125#include "MLog.h"
126#include "MLogManip.h"
127#include "MBinning.h"
128
129using namespace std;
130
131ClassImp(MMuonCalibParCalc);
132
133static const TString gsDefName = "MMuonCalibParCalc";
134static const TString gsDefTitle = "Calculate new image parameters";
135
136// -------------------------------------------------------------------------
137//
138// Default constructor.
139//
140MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title)
141 : fNameCerPhot("MCerPhotEvt")
142{
143 fName = name ? name : gsDefName.Data();
144 fTitle = title ? title : gsDefTitle.Data();
145
146 fPreCuts[0] = 180.;
147 fPreCuts[1] = 400.;
148 fPreCuts[2] = 50.;
149 fPreCuts[3] = 2000.;
150 fPreCuts[4] = 150.;
151
152 fMargin = 60.;
153 fArcPhiThres = 100.;
154 fArcWidthThres = 100.;
155
156 fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
157 fDisablePreCuts = kFALSE; // By default the pre cuts will be applied.
158}
159
160// -------------------------------------------------------------------------
161//
162Int_t MMuonCalibParCalc::PreProcess(MParList *pList)
163{
164 fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fNameCerPhot), "MCerPhotEvt");
165 if (!fCerPhotEvt)
166 {
167 *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
168 return kFALSE;
169 }
170
171 fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
172 if (!fGeomCam)
173 {
174 *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
175 return kFALSE;
176 }
177
178 fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar", "MMuonCalibPar");
179 if (!fMuonCalibPar)
180 {
181 *fLog << dbginf << "MMuonCalibPar missing in Parameter List... aborting." << endl;
182 return kFALSE;
183 }
184
185 fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", "MMuonSearchPar");
186 if (!fMuonSearchPar)
187 {
188 *fLog << dbginf << "MMuonSearchPar missing in Parameter List... aborting." << endl;
189 return kFALSE;
190 }
191
192 return kTRUE;
193}
194
195// --------------------------------------------------------------------------
196//
197// This function fill the histograms in order to get muon parameters.
198// For the evaluation of the Arc Width, we use only the signals in the inner
199// part.
200//
201void MMuonCalibParCalc::FillHist()
202{
203 Float_t MuonSize = 0.;
204
205 Int_t binnumphi = fMuonCalibPar->fArcPhiBinNum;
206 Int_t binnumwid = fMuonCalibPar->fArcWidthBinNum;
207
208 // preparation for a histgram
209 MBinning binsphi;
210 binsphi.SetEdges(binnumphi,
211 fMuonCalibPar->fArcPhiHistStartVal,
212 fMuonCalibPar->fArcPhiHistEndVal);
213 binsphi.Apply(*(fMuonCalibPar->fHistPhi));
214
215 MBinning binswid;
216 binswid.SetEdges(binnumwid,
217 fMuonCalibPar->fArcWidthHistStartVal,
218 fMuonCalibPar->fArcWidthHistEndVal);
219 binswid.Apply(*(fMuonCalibPar->fHistWidth));
220
221 const Int_t entries = (*fCerPhotEvt).GetNumPixels();
222
223 // the position of the center of a muon ring
224 const Float_t cenx = (*fMuonSearchPar).GetCenterX();
225 const Float_t ceny = (*fMuonSearchPar).GetCenterY();
226
227 for (Int_t i=0; i<entries; i++ )
228 {
229 MCerPhotPix &pix = (*fCerPhotEvt)[i];
230
231 const MGeomPix &gpix = (*fGeomCam)[pix.GetPixId()];
232
233 const Float_t dx = gpix.GetX() - cenx;
234 const Float_t dy = gpix.GetY() - ceny;
235
236 const Float_t dist = TMath::Sqrt(dx*dx+dy*dy);
237
238 Float_t ang = TMath::ACos(dx/dist);
239 if(dy>0)
240 ang *= -1.0;
241
242 // if the signal is not near the estimated circle, it is ignored.
243 if(dist < (*fMuonSearchPar).GetRadius() + fMuonCalibPar->GetMargin()
244 && dist > (*fMuonSearchPar).GetRadius() - fMuonCalibPar->GetMargin())
245 {
246 // check whether ummapped pixel is used or not.
247 // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information.
248 if(pix.IsPixelUnmapped())
249 {
250 fMuonCalibPar->SetUseUnmap();
251 continue;
252 }
253 fMuonCalibPar->fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons());
254 MuonSize += pix.GetNumPhotons();
255 }
256
257 // use only the inner pixles. This is geometry dependent. This has to
258 // be fixed!
259 if(pix.GetPixId()>397)
260 continue;
261
262 fMuonCalibPar->fHistWidth
263 ->Fill(dist*(*fGeomCam).GetConvMm2Deg(), pix.GetNumPhotons());
264 }
265
266 fMuonCalibPar->SetMuonSize(MuonSize);
267
268 // error estimation (temporaly)
269 // The error is estimated from the signal. In order to do so, we have to
270 // once convert the signal from ADC to photo-electron. Then we can get
271 // the fluctuation such as F-factor*sqrt(phe).
272 // Up to now, the error of pedestal is not taken into accout. This is not
273 // of course correct. We will include this soon.
274 Double_t ADC2PhEl = 0.14;
275 Double_t Ffactor = 1.4;
276 for(Int_t i=0; i<binnumphi+1; i++)
277 {
278 Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
279 {
280 fMuonCalibPar->fHistPhi->SetBinError(i, rougherr);
281 }
282 }
283 for(Int_t i=0; i<binnumwid+1; i++)
284 {
285 Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
286 {
287 fMuonCalibPar->fHistWidth->SetBinError(i, rougherr);
288 }
289 }
290}
291
292// --------------------------------------------------------------------------
293//
294// Photon distribution along the estimated circle is fitted with theoritical
295// function in order to get some more information such as Arc Phi and Arc
296// Length.
297//
298void MMuonCalibParCalc::CalcPhi()
299{
300 Float_t thres = fMuonCalibPar->GetArcPhiThres();
301 Float_t startval = fMuonCalibPar->fArcPhiHistStartVal;
302 Float_t endval = fMuonCalibPar->fArcPhiHistEndVal;
303 Int_t binnum = fMuonCalibPar->fArcPhiBinNum;
304
305 Float_t convbin2val = (endval-startval)/(Float_t)binnum;
306
307 // adjust the peak to 0
308 Float_t maxval = 0.;
309 Int_t maxbin = 0;
310 maxval = fMuonCalibPar->fHistPhi->GetMaximum();
311 maxbin = fMuonCalibPar->fHistPhi->GetMaximumBin();
312 fMuonCalibPar->SetPeakPhi(180.-(Float_t)(maxbin-1.)*convbin2val);
313 TArrayD tmp;
314 tmp.Set(binnum+1);
315 for(Int_t i=1; i<binnum+1; i++)
316 {
317 tmp[i] = fMuonCalibPar->fHistPhi->GetBinContent(i);
318 }
319 for(Int_t i=1; i<binnum+1; i++)
320 {
321 Int_t id;
322 id = i + (maxbin-(Int_t)((Float_t)binnum/2.)-1);
323 if(id>binnum)
324 {
325 id-=(binnum);
326 }
327 if(id<=0)
328 {
329 id+=(binnum);
330 }
331 fMuonCalibPar->fHistPhi->SetBinContent(i,tmp[id]);
332 }
333 maxbin = (Int_t)((Float_t)binnum/2.)+1;
334
335 // Determination of fitting region
336 // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore,
337 // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!!
338 Float_t startfitval = 0.;
339 Float_t endfitval = 0.;
340 Bool_t IsInMaxim = kFALSE;
341 Int_t effbinnum = 0;
342 for(Int_t i=1; i<binnum+1; i++)
343 {
344 Float_t content = fMuonCalibPar->fHistPhi->GetBinContent(i);
345 Float_t content_pre = fMuonCalibPar->fHistPhi->GetBinContent(i-1);
346
347 if(content > thres && content_pre < thres)
348 {
349 startfitval = (Float_t)(i-1)*convbin2val+startval;
350 }
351 if(i==maxbin)
352 IsInMaxim = kTRUE;
353
354 if(content < thres && IsInMaxim == kTRUE)
355 {
356 endfitval = (Float_t)(i-1)*convbin2val+startval;
357 break;
358 }
359 endfitval = endval;
360 }
361
362 effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
363
364 fMuonCalibPar->SetArcPhi(endfitval-startfitval);
365
366 fMuonCalibPar->SetArcLength( fMuonCalibPar->GetArcPhi()*TMath::DegToRad()*(*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg());
367
368 if(fEnableImpactCalc)
369 CalcImpact(effbinnum, startfitval, endfitval);
370}
371
372// --------------------------------------------------------------------------
373//
374// An impact parameter is calculated by fitting the histogram of photon
375// distribution along the circle with a theoritical model.
376// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.
377// The function (6) is used here.)
378//
379// By default this calculation is suppressed because this calculation is
380// very time consuming. If you want to calculate an impact parameter,
381// you can call the function of EnableImpactCalc().
382//
383void MMuonCalibParCalc::CalcImpact
384(Int_t effbinnum, Float_t startfitval, Float_t endfitval)
385{
386 // Fit the distribution with Vacanti function. The function is different
387 // for the impact parameter of inside or outside of our reflector.
388 // Then two different functions are applied to the photon distribution,
389 // and the one which give us smaller chisquare value is taken as a
390 // proper one.
391 Double_t val1,err1,val2,err2;
392 // impact parameter inside mirror radius (8.5m)
393 TString func1;
394 Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
395 tmpval = sin(2.*tmpval)*8.5;
396 func1 += "[0]*";
397 func1 += tmpval;
398 func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
399
400 TF1 f1("f1",func1,startfitval,endfitval);
401 f1.SetParameters(2000,3,0);
402 f1.SetParLimits(1,0,8.5);
403 f1.SetParLimits(2,-180.,180.);
404
405 fMuonCalibPar->fHistPhi->Fit("f1","RQEM");
406
407 Float_t chi1 = -1;
408 Float_t chi2 = -1;
409 if(effbinnum>3)
410 chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
411
412 gMinuit->GetParameter(1,val1,err1); // get the estimated IP
413 Float_t estip1 = val1;
414
415 // impact parameter beyond mirror area (8.5m)
416 TString func2;
417 Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
418 tmpval2 = sin(2.*tmpval2)*8.5*2.;
419 func2 += "[0]*";
420 func2 += tmpval2;
421 func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
422
423 TF1 f2("f2",func2,startfitval,endfitval);
424 f2.SetParameters(2000,20,0);
425 f2.SetParLimits(1,8.5,300.);
426 f2.SetParLimits(2,-180.,180.);
427
428 fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");
429
430 if(effbinnum>3)
431 chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
432
433 gMinuit->GetParameter(1,val2,err2); // get the estimated IP
434 Float_t estip2 = val2;
435
436 if(effbinnum<=3)
437 {
438 fMuonCalibPar->SetEstImpact(-1.);
439 }
440 if(chi2 > chi1)
441 {
442 fMuonCalibPar->SetEstImpact(estip1);
443 fMuonCalibPar->SetChiArcPhi(chi1);
444 }
445 else
446 {
447 fMuonCalibPar->SetEstImpact(estip2);
448 fMuonCalibPar->SetChiArcPhi(chi2);
449 }
450}
451
452// --------------------------------------------------------------------------
453//
454// Photon distribution of distance from the center of estimated ring is
455// fitted in order to get some more information such as ARC WIDTH which
456// can represent to the PSF of our reflector.
457//
458Float_t MMuonCalibParCalc::CalcWidth()
459{
460 Float_t startval = fMuonCalibPar->fArcWidthHistStartVal;
461 Float_t endval = fMuonCalibPar->fArcWidthHistEndVal;
462 Int_t binnum = fMuonCalibPar->fArcWidthBinNum;
463 Float_t thres = fMuonCalibPar->GetArcWidthThres();
464
465 Float_t convbin2val = (endval - startval)
466 /(Float_t)binnum;
467
468 // determination of fitting region
469 Int_t maxbin = fMuonCalibPar->fHistWidth->GetMaximumBin();
470 Float_t startfitval = 0.;
471 Float_t endfitval = 0.;
472 Bool_t IsInMaxim = kFALSE;
473 Int_t effbinnum = 0;
474 for(Int_t i=1; i<binnum+1; i++)
475 {
476 Float_t content = fMuonCalibPar->fHistWidth->GetBinContent(i);
477 Float_t content_pre = fMuonCalibPar->fHistWidth->GetBinContent(i-1);
478
479 if(content > thres)
480 effbinnum++;
481
482 if(content > thres && content_pre < thres)
483 {
484 startfitval = (Float_t)(i-4)*convbin2val + startval;
485 if(startfitval<0.) startfitval = 0.;
486 }
487 if(i==maxbin)
488 IsInMaxim = kTRUE;
489
490 if(content < thres && IsInMaxim == kTRUE)
491 {
492 endfitval = (Float_t)(i+2)*convbin2val + startval;
493 if(endfitval>180.) endfitval = 180.;
494 break;
495 }
496 endfitval = endval;
497 }
498 effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
499
500 TF1 f1("f1","gaus",startfitval,endfitval);
501
502 fMuonCalibPar->fHistWidth->Fit("f1","QR","",startfitval,endfitval);
503
504 if(effbinnum>3)
505 fMuonCalibPar->SetChiArcWidth(f1.GetChisquare()/((Float_t)(effbinnum-3)));
506
507 Double_t val,err;
508 gMinuit->GetParameter(2,val,err); // get the sigma value
509
510 return val;
511}
512
513// --------------------------------------------------------------------------
514//
515// Calculation of muon parameters
516//
517Int_t MMuonCalibParCalc::Calc(const Float_t *cuts)
518{
519 // sanity check
520 if((*fCerPhotEvt).GetNumPixels() < 3)
521 return kCONTINUE;
522
523 // If an event does not seem to be like muon, the calculation will be skipped.
524 if((*fMuonSearchPar).IsNoMuon())
525 return kCONTINUE;
526
527 // Pre Cuts 1
528 if(!fDisablePreCuts)
529 {
530 if((*fMuonSearchPar).GetRadius() < cuts[0] || (*fMuonSearchPar).GetRadius() > cuts[1])
531 {
532 (*fMuonSearchPar).SetNoMuon();
533 return kCONTINUE;
534 }
535 if((*fMuonSearchPar).GetDeviation() > cuts[2])
536 {
537 (*fMuonSearchPar).SetNoMuon();
538 return kCONTINUE;
539 }
540 }
541
542 // initialization
543 (*fMuonCalibPar).Reset();
544
545 // Fill signals to histograms
546 FillHist();
547
548 // Calculation of Arc Phi etc...
549 CalcPhi();
550
551 // Pre Cuts 2
552 if(!fDisablePreCuts)
553 {
554 if(fMuonCalibPar->GetMuonSize() < cuts[3]
555 || fMuonCalibPar->GetArcPhi() < cuts[4])
556 {
557 (*fMuonSearchPar).SetNoMuon();
558 return kCONTINUE;
559 }
560 }
561
562 // Calculation of Arc Width etc...
563 fMuonCalibPar->SetArcWidth(CalcWidth());
564
565 return kTRUE;
566}
567
568
569// -------------------------------------------------------------------------
570//
571Int_t MMuonCalibParCalc::Process()
572{
573 fMuonCalibPar->SetMargin(fMargin);
574 fMuonCalibPar->SetArcPhiThres(fArcPhiThres);
575 fMuonCalibPar->SetArcWidthThres(fArcWidthThres);
576
577 if(!Calc(fPreCuts))
578 return kCONTINUE;
579
580 return kTRUE;
581}
582
583void MMuonCalibParCalc::SetPreCuts
584(Float_t radcutlow, Float_t radcuthigh, Float_t devcuthigh,
585 Float_t musizecutlow, Float_t arcphicutlow)
586{
587 fPreCuts[0] = radcutlow;
588 fPreCuts[1] = radcuthigh;
589 fPreCuts[2] = devcuthigh;
590 fPreCuts[3] = musizecutlow;
591 fPreCuts[4] = arcphicutlow;
592}
593
Note: See TracBrowser for help on using the repository browser.