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

Last change on this file since 5332 was 5332, checked in by mase, 20 years ago
*** empty log message ***
File size: 17.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 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 if(pix.GetPixId()>397)
258 continue; // use only the inner pixles
259
260 fMuonCalibPar->fHistWidth
261 ->Fill(dist*(*fGeomCam).GetConvMm2Deg(), pix.GetNumPhotons());
262 }
263
264 fMuonCalibPar->SetMuonSize(MuonSize);
265
266 // error estimation (temporaly)
267 // The error is estimated from the signal. In order to do so, we have to
268 // once convert the signal from ADC to photo-electron. Then we can get
269 // the fluctuation such as F-factor*sqrt(phe).
270 // Up to now, the error of pedestal is not taken into accout. This is not
271 // of course correct. We will include this soon.
272 Double_t ADC2PhEl = 0.14;
273 Double_t Ffactor = 1.4;
274 for(Int_t i=0; i<binnumphi+1; i++)
275 {
276 Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
277 {
278 fMuonCalibPar->fHistPhi->SetBinError(i, rougherr);
279 }
280 }
281 for(Int_t i=0; i<binnumwid+1; i++)
282 {
283 Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
284 {
285 fMuonCalibPar->fHistWidth->SetBinError(i, rougherr);
286 }
287 }
288}
289
290// --------------------------------------------------------------------------
291//
292// Photon distribution along the estimated circle is fitted with theoritical
293// function in order to get some more information such as Arc Phi and Arc
294// Length.
295//
296void MMuonCalibParCalc::CalcPhi()
297{
298 Float_t thres = fMuonCalibPar->GetArcPhiThres();
299 Float_t startval = fMuonCalibPar->fArcPhiHistStartVal;
300 Float_t endval = fMuonCalibPar->fArcPhiHistEndVal;
301 Int_t binnum = fMuonCalibPar->fArcPhiBinNum;
302
303 Float_t convbin2val = (endval-startval)/(Float_t)binnum;
304
305 // adjust the peak to 0
306 Float_t maxval = 0.;
307 Int_t maxbin = 0;
308 maxval = fMuonCalibPar->fHistPhi->GetMaximum();
309 maxbin = fMuonCalibPar->fHistPhi->GetMaximumBin();
310 fMuonCalibPar->SetPeakPhi(180.-(Float_t)(maxbin-1.)*convbin2val);
311 TArrayD tmp;
312 tmp.Set(binnum+1);
313 for(Int_t i=1; i<binnum+1; i++)
314 {
315 tmp[i] = fMuonCalibPar->fHistPhi->GetBinContent(i);
316 }
317 for(Int_t i=1; i<binnum+1; i++)
318 {
319 Int_t id;
320 id = i + (maxbin-(Int_t)((Float_t)binnum/2.)-1);
321 if(id>binnum)
322 {
323 id-=(binnum);
324 }
325 if(id<=0)
326 {
327 id+=(binnum);
328 }
329 fMuonCalibPar->fHistPhi->SetBinContent(i,tmp[id]);
330 }
331 maxbin = (Int_t)((Float_t)binnum/2.)+1;
332
333 // Determination of fitting region
334 // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore,
335 // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!!
336 Float_t startfitval = 0.;
337 Float_t endfitval = 0.;
338 Bool_t IsInMaxim = kFALSE;
339 Int_t effbinnum = 0;
340 for(Int_t i=1; i<binnum+1; i++)
341 {
342 Float_t content = fMuonCalibPar->fHistPhi->GetBinContent(i);
343 Float_t content_pre = fMuonCalibPar->fHistPhi->GetBinContent(i-1);
344
345 if(content > thres && content_pre < thres)
346 {
347 startfitval = (Float_t)(i-1)*convbin2val+startval;
348 }
349 if(i==maxbin)
350 IsInMaxim = kTRUE;
351
352 if(content < thres && IsInMaxim == kTRUE)
353 {
354 endfitval = (Float_t)(i-1)*convbin2val+startval;
355 break;
356 }
357 endfitval = endval;
358 }
359
360 effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
361
362 fMuonCalibPar->SetArcPhi(endfitval-startfitval);
363
364 fMuonCalibPar->SetArcLength( fMuonCalibPar->GetArcPhi()*TMath::DegToRad()*(*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg());
365
366 if(fEnableImpactCalc)
367 CalcImpact(effbinnum, startfitval, endfitval);
368}
369
370// --------------------------------------------------------------------------
371//
372// An impact parameter is calculated by fitting the histogram of photon
373// distribution along the circle with a theoritical model.
374// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.
375// The function (6) is used here.)
376//
377// By default this calculation is suppressed because this calculation is
378// very time consuming. If you want to calculate an impact parameter,
379// you can call the function of EnableImpactCalc().
380//
381void MMuonCalibParCalc::CalcImpact
382(Int_t effbinnum, Float_t startfitval, Float_t endfitval)
383{
384 // Fit the distribution with Vacanti function. The function is different
385 // for the impact parameter of inside or outside of our reflector.
386 // Then two different functions are applied to the photon distribution,
387 // and the one which give us smaller chisquare value is taken as a
388 // proper one.
389 Double_t val1,err1,val2,err2;
390 // impact parameter inside mirror radius (8.5m)
391 TString func1;
392 Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
393 tmpval = sin(2.*tmpval)*8.5;
394 func1 += "[0]*";
395 func1 += tmpval;
396 func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
397
398 TF1 f1("f1",func1,startfitval,endfitval);
399 f1.SetParameters(2000,3,0);
400 f1.SetParLimits(1,0,8.5);
401 f1.SetParLimits(2,-180.,180.);
402
403 fMuonCalibPar->fHistPhi->Fit("f1","RQEM");
404
405 Float_t chi1 = -1;
406 Float_t chi2 = -1;
407 if(effbinnum>3)
408 chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
409
410 gMinuit->GetParameter(1,val1,err1); // get the estimated IP
411 Float_t estip1 = val1;
412
413 // impact parameter beyond mirror area (8.5m)
414 TString func2;
415 Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
416 tmpval2 = sin(2.*tmpval2)*8.5*2.;
417 func2 += "[0]*";
418 func2 += tmpval2;
419 func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
420
421 TF1 f2("f2",func2,startfitval,endfitval);
422 f2.SetParameters(2000,20,0);
423 f2.SetParLimits(1,8.5,300.);
424 f2.SetParLimits(2,-180.,180.);
425
426 fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");
427
428 if(effbinnum>3)
429 chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
430
431 gMinuit->GetParameter(1,val2,err2); // get the estimated IP
432 Float_t estip2 = val2;
433
434 if(effbinnum<=3)
435 {
436 fMuonCalibPar->SetEstImpact(-1.);
437 }
438 if(chi2 > chi1)
439 {
440 fMuonCalibPar->SetEstImpact(estip1);
441 fMuonCalibPar->SetChiArcPhi(chi1);
442 }
443 else
444 {
445 fMuonCalibPar->SetEstImpact(estip2);
446 fMuonCalibPar->SetChiArcPhi(chi2);
447 }
448}
449
450// --------------------------------------------------------------------------
451//
452// Photon distribution of distance from the center of estimated ring is
453// fitted in order to get some more information such as ARC WIDTH which
454// can represent to the PSF of our reflector.
455//
456Float_t MMuonCalibParCalc::CalcWidth()
457{
458 Float_t startval = fMuonCalibPar->fArcWidthHistStartVal;
459 Float_t endval = fMuonCalibPar->fArcWidthHistEndVal;
460 Int_t binnum = fMuonCalibPar->fArcWidthBinNum;
461 Float_t thres = fMuonCalibPar->GetArcWidthThres();
462
463 Float_t convbin2val = (endval - startval)
464 /(Float_t)binnum;
465
466 // determination of fitting region
467 Int_t maxbin = fMuonCalibPar->fHistWidth->GetMaximumBin();
468 Float_t startfitval = 0.;
469 Float_t endfitval = 0.;
470 Bool_t IsInMaxim = kFALSE;
471 Int_t effbinnum = 0;
472 for(Int_t i=1; i<binnum+1; i++)
473 {
474 Float_t content = fMuonCalibPar->fHistWidth->GetBinContent(i);
475 Float_t content_pre = fMuonCalibPar->fHistWidth->GetBinContent(i-1);
476
477 if(content > thres)
478 effbinnum++;
479
480 if(content > thres && content_pre < thres)
481 {
482 startfitval = (Float_t)(i-4)*convbin2val + startval;
483 if(startfitval<0.) startfitval = 0.;
484 }
485 if(i==maxbin)
486 IsInMaxim = kTRUE;
487
488 if(content < thres && IsInMaxim == kTRUE)
489 {
490 endfitval = (Float_t)(i+2)*convbin2val + startval;
491 if(endfitval>180.) endfitval = 180.;
492 break;
493 }
494 endfitval = endval;
495 }
496 effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
497
498 TF1 f1("f1","gaus",startfitval,endfitval);
499
500 fMuonCalibPar->fHistWidth->Fit("f1","QR","",startfitval,endfitval);
501
502 if(effbinnum>3)
503 fMuonCalibPar->SetChiArcWidth(f1.GetChisquare()/((Float_t)(effbinnum-3)));
504
505 Double_t val,err;
506 gMinuit->GetParameter(2,val,err); // get the sigma value
507
508 return val;
509}
510
511// --------------------------------------------------------------------------
512//
513// Calculation of muon parameters
514//
515Int_t MMuonCalibParCalc::Calc(const Float_t *cuts)
516{
517 // sanity check
518 if((*fCerPhotEvt).GetNumPixels() < 3)
519 return kCONTINUE;
520
521 // If an event does not seem to be like muon, the calculation will be skipped.
522 if((*fMuonSearchPar).IsNoMuon())
523 return kCONTINUE;
524
525 // Pre Cuts 1
526 if(!fDisablePreCuts)
527 {
528 if((*fMuonSearchPar).GetRadius() < cuts[0] || (*fMuonSearchPar).GetRadius() > cuts[1])
529 {
530 (*fMuonSearchPar).SetNoMuon();
531 return kCONTINUE;
532 }
533 if((*fMuonSearchPar).GetDeviation() > cuts[2])
534 {
535 (*fMuonSearchPar).SetNoMuon();
536 return kCONTINUE;
537 }
538 }
539
540 // initialization
541 (*fMuonCalibPar).Reset();
542
543 // Fill signals to histograms
544 FillHist();
545
546 // Calculation of Arc Phi etc...
547 CalcPhi();
548
549 // Pre Cuts 2
550 if(!fDisablePreCuts)
551 {
552 if(fMuonCalibPar->GetMuonSize() < cuts[3]
553 || fMuonCalibPar->GetArcPhi() < cuts[4])
554 {
555 (*fMuonSearchPar).SetNoMuon();
556 return kCONTINUE;
557 }
558 }
559
560 // Calculation of Arc Width etc...
561 fMuonCalibPar->SetArcWidth(CalcWidth());
562
563 return kTRUE;
564}
565
566
567// -------------------------------------------------------------------------
568//
569Int_t MMuonCalibParCalc::Process()
570{
571 fMuonCalibPar->SetMargin(fMargin);
572 fMuonCalibPar->SetArcPhiThres(fArcPhiThres);
573 fMuonCalibPar->SetArcWidthThres(fArcWidthThres);
574
575 if(!Calc(fPreCuts))
576 return kCONTINUE;
577
578 return kTRUE;
579}
580
581void MMuonCalibParCalc::SetPreCuts
582(Float_t radcutlow, Float_t radcuthigh, Float_t devcuthigh,
583 Float_t musizecutlow, Float_t arcphicutlow)
584{
585 fPreCuts[0] = radcutlow;
586 fPreCuts[1] = radcuthigh;
587 fPreCuts[2] = devcuthigh;
588 fPreCuts[3] = musizecutlow;
589 fPreCuts[4] = arcphicutlow;
590}
591
Note: See TracBrowser for help on using the repository browser.