Index: /trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc (revision 5209) +++ /trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc (revision 5210) @@ -16,6 +16,6 @@ ! ! -! Author(s): Keiichi Mase 09/2004 -! Markus Meyer 09/2004 +! Author(s): Keiichi Mase 10/2004 +! Markus Meyer 10/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 @@ -72,34 +72,44 @@ fTitle = title ? title : "Muon calibration parameters"; - fHistPhi.SetName("HistPhi"); - fHistPhi.SetTitle("HistPhi"); - fHistPhi.SetXTitle("phi [deg.]"); - fHistPhi.SetYTitle("sum of ADC"); - fHistPhi.SetDirectory(NULL); - fHistPhi.SetFillStyle(4000); - fHistPhi.UseCurrentStyle(); - - fHistWid.SetName("HistWid"); - fHistWid.SetTitle("HistWid"); - fHistWid.SetXTitle("distance from the ring center [deg.]"); - fHistWid.SetYTitle("sum of ADC"); - fHistWid.SetDirectory(NULL); - fHistWid.SetFillStyle(4000); - fHistWid.UseCurrentStyle(); - - fKeepHist = kFALSE; // normally histgram is not kept + fHistPhi = new TH1F; + fHistWidth = new TH1F; + + fHistPhi->SetName("HistPhi"); + fHistPhi->SetTitle("HistPhi"); + fHistPhi->SetXTitle("phi [deg.]"); + fHistPhi->SetYTitle("sum of ADC"); + fHistPhi->SetDirectory(NULL); + fHistPhi->SetFillStyle(4000); + fHistPhi->UseCurrentStyle(); + + fHistWidth->SetName("HistWidth"); + fHistWidth->SetTitle("HistWidth"); + fHistWidth->SetXTitle("distance from the ring center [deg.]"); + fHistWidth->SetYTitle("sum of ADC"); + fHistWidth->SetDirectory(NULL); + fHistWidth->SetFillStyle(4000); + fHistWidth->UseCurrentStyle(); + fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped. fDisablePreCuts = kFALSE; // By default the pre cuts will be applied. - fUseCleanForWid = kFALSE; // By default all the pixels will be used for the histogram of arc width. + fUseCleanForWidth = kFALSE; // By default all the pixels will be used for the histogram of arc width. fMargin = 60.; // in mm fArcPhiThres = 100.; - fArcWidThres = 100.; + fArcWidthThres = 100.; fArcPhiBinNum = 20; fArcPhiHistStartVal = -180.; // deg. fArcPhiHistEndVal = 180.; // deg. - fArcWidBinNum = 28; - fArcWidHistStartVal = 0.3; // deg. - fArcWidHistEndVal = 1.7; // deg. + fArcWidthBinNum = 28; + fArcWidthHistStartVal = 0.3; // deg. + fArcWidthHistEndVal = 1.7; // deg. +} + +// -------------------------------------------------------------------------- +// +MMuonCalibPar::~MMuonCalibPar() +{ + delete fHistPhi; + delete fHistWidth; } @@ -108,9 +118,9 @@ void MMuonCalibPar::Reset() { - fArcLen = -1.; + fArcLength = -1.; fArcPhi = 0.; - fArcWid = -1.; + fArcWidth = -1.; fChiArcPhi = -1.; - fChiArcWid = -1.; + fChiArcWidth = -1.; fMuonSize = 0.; fEstImpact = -1.; @@ -122,7 +132,7 @@ // // This function fill the histograms in order to get muon parameters. -// For the evaluation of the Arc Width, we use only the signals in the inner part. -// You can use the image after the cleaning by using the function of UseCleanForWid(). -// See the manual of MMuonCalibParCalc. +// For the evaluation of the Arc Width, we use only the signals in the inner +// part. You can use the image after the cleaning by using the function of +// UseCleanForWidth(). See the manual of MMuonCalibParCalc. // void MMuonCalibPar::FillHist @@ -133,15 +143,15 @@ MBinning binsphi; binsphi.SetEdges(fArcPhiBinNum, fArcPhiHistStartVal, fArcPhiHistEndVal); - binsphi.Apply(fHistPhi); + binsphi.Apply(*fHistPhi); MBinning binswid; - binswid.SetEdges(fArcWidBinNum, fArcWidHistStartVal, fArcWidHistEndVal); - binswid.Apply(fHistWid); + binswid.SetEdges(fArcWidthBinNum, fArcWidthHistStartVal, fArcWidthHistEndVal); + binswid.Apply(*fHistWidth); const Int_t entries = evt.GetNumPixels(); // the position of the center of a muon ring - const Float_t cenx = musearch.GetCenX(); - const Float_t ceny = musearch.GetCenY(); + const Float_t cenx = musearch.GetCenterX(); + const Float_t ceny = musearch.GetCenterY(); for (Int_t i=0; i musearch.GetRad() - fMargin) + if(dist < musearch.GetRadius() + fMargin && dist > musearch.GetRadius() - fMargin) { // check whether ummapped pixel is used or not. @@ -170,5 +180,5 @@ continue; } - fHistPhi.Fill(ang*kRad2Deg, pix.GetNumPhotons()); + fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons()); fMuonSize += pix.GetNumPhotons(); } @@ -177,5 +187,5 @@ continue; // use only the inner pixles - if(fUseCleanForWid) + if(fUseCleanForWidth) { if(!pix.IsPixelUsed()) @@ -183,5 +193,5 @@ } - fHistWid.Fill(dist*geom.GetConvMm2Deg(), pix.GetNumPhotons()); + fHistWidth->Fill(dist*geom.GetConvMm2Deg(), pix.GetNumPhotons()); } @@ -191,20 +201,20 @@ // once convert the signal from ADC to photo-electron. Then we can get // the fluctuation such as F-factor*sqrt(phe). - // Up to now, the error of pedestal is not taken into accout. This is not of - // course correct. We will include this soon. + // Up to now, the error of pedestal is not taken into accout. This is not + // of course correct. We will include this soon. Double_t ADC2PhEl = 0.14; Double_t Ffactor = 1.4; for(Int_t i=0; iGetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor; + { + fHistPhi->SetBinError(i, rougherr); } } - for(Int_t i=0; iGetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor; + { + fHistWidth->SetBinError(i, rougherr); } } @@ -226,6 +236,6 @@ Float_t maxval = 0.; Int_t maxbin = 0; - maxval = fHistPhi.GetMaximum(); - maxbin = fHistPhi.GetMaximumBin(); + maxval = fHistPhi->GetMaximum(); + maxbin = fHistPhi->GetMaximumBin(); fPeakPhi = 180.-(Float_t)(maxbin-1)*convbin2val; TArrayD tmp; @@ -233,5 +243,5 @@ for(Int_t i=1; iGetBinContent(i); } for(Int_t i=1; iSetBinContent(i,tmp[id]); } maxbin = (Int_t)((Float_t)fArcPhiBinNum/2.)+1; @@ -260,6 +270,6 @@ for(Int_t i=1; iGetBinContent(i); + Float_t content_pre = fHistPhi->GetBinContent(i-1); if(content > fArcPhiThres && content_pre < fArcPhiThres) @@ -281,5 +291,5 @@ fArcPhi = effbinnum*convbin2val; - fArcLen = fArcPhi*3.1415926/180.*musearch.GetRad()*geom.GetConvMm2Deg(); + fArcLength = fArcPhi*3.1415926/180.*musearch.GetRadius()*geom.GetConvMm2Deg(); if(fEnableImpactCalc) @@ -310,5 +320,5 @@ // impact parameter inside mirror radius (8.5m) TString func1; - Float_t tmpval = musearch.GetRad()*geom.GetConvMm2Deg()*3.1415926/180.; + Float_t tmpval = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.; tmpval = sin(2.*tmpval)*8.5; func1 += "[0]*"; @@ -321,12 +331,5 @@ f1.SetParLimits(2,-180.,180.); - if(fKeepHist) - { - fHistPhi.Fit("f1","RQEM"); - } - else - { - fHistPhi.Fit("f1","RQEM0"); - } + fHistPhi->Fit("f1","RQEM"); Float_t chi1 = -1; @@ -340,5 +343,5 @@ // impact parameter beyond mirror area (8.5m) TString func2; - Float_t tmpval2 = musearch.GetRad()*geom.GetConvMm2Deg()*3.1415926/180.; + Float_t tmpval2 = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.; tmpval2 = sin(2.*tmpval2)*8.5*2.; func2 += "[0]*"; @@ -351,12 +354,5 @@ f2.SetParLimits(2,-180.,180.); - if(fKeepHist) - { - fHistPhi.Fit("f2","RQEM+"); - } - else - { - fHistPhi.Fit("f2","RQEM0+"); - } + fHistPhi->Fit("f2","RQEM+"); if(effbinnum>3) @@ -388,28 +384,28 @@ // can represent to the PSF of our reflector. // -Float_t MMuonCalibPar::CalcWid +Float_t MMuonCalibPar::CalcWidth (const MGeomCam &geom, const MCerPhotEvt &evt, const MMuonSearchPar &musearch) { - Float_t convbin2val = (fArcWidHistEndVal - fArcWidHistStartVal) - /(Float_t)fArcWidBinNum; + Float_t convbin2val = (fArcWidthHistEndVal - fArcWidthHistStartVal) + /(Float_t)fArcWidthBinNum; // determination of fitting region - Int_t maxbin = fHistWid.GetMaximumBin(); + Int_t maxbin = fHistWidth->GetMaximumBin(); Float_t startfitval = 0.; Float_t endfitval = 0.; Bool_t IsInMaxim = kFALSE; Int_t effbinnum = 0; - for(Int_t i=1; i fArcWidThres) + Float_t content = fHistWidth->GetBinContent(i); + Float_t content_pre = fHistWidth->GetBinContent(i-1); + + if(content > fArcWidthThres) effbinnum++; - if(content > fArcWidThres && content_pre < fArcWidThres) + if(content > fArcWidthThres && content_pre < fArcWidthThres) { - startfitval = (Float_t)(i-4)*convbin2val + fArcWidHistStartVal; + startfitval = (Float_t)(i-4)*convbin2val + fArcWidthHistStartVal; if(startfitval<0.) startfitval = 0.; } @@ -417,11 +413,11 @@ IsInMaxim = kTRUE; - if(content < fArcWidThres && IsInMaxim == kTRUE) + if(content < fArcWidthThres && IsInMaxim == kTRUE) { - endfitval = (Float_t)(i+2)*convbin2val + fArcWidHistStartVal; + endfitval = (Float_t)(i+2)*convbin2val + fArcWidthHistStartVal; if(endfitval>180.) endfitval = 180.; break; } - endfitval = fArcWidHistEndVal; + endfitval = fArcWidthHistEndVal; } effbinnum = (Int_t)((endfitval-startfitval)/convbin2val); @@ -429,15 +425,8 @@ TF1 f1("f1","gaus",startfitval,endfitval); - if(fKeepHist) - { - fHistWid.Fit("f1","QR","",startfitval,endfitval); - } - else - { - fHistWid.Fit("f1","QR0","",startfitval,endfitval); - } + fHistWidth->Fit("f1","QR","",startfitval,endfitval); if(effbinnum>3) - fChiArcWid = f1.GetChisquare()/((Float_t)(effbinnum-3)); + fChiArcWidth = f1.GetChisquare()/((Float_t)(effbinnum-3)); Double_t val,err; @@ -453,5 +442,5 @@ Int_t MMuonCalibPar::Calc (const MGeomCam &geom, const MCerPhotEvt &evt, - const MMuonSearchPar &musearch) + MMuonSearchPar &musearch, const Float_t *cuts) { // sanity check @@ -459,5 +448,5 @@ return kCONTINUE; - // if an event does not seem to be like muon, it will be skipped. + // If an event does not seem to be like muon, the calculation will be skipped. if(musearch.IsNoMuon()) return kCONTINUE; @@ -466,8 +455,14 @@ if(!fDisablePreCuts) { - if(musearch.GetRad() < 180. || musearch.GetRad() > 400.) - return kCONTINUE; - if(musearch.GetDev() > 50.) - return kCONTINUE; + if(musearch.GetRadius() < cuts[0] || musearch.GetRadius() > cuts[1]) + { + musearch.SetNoMuon(); + return kCONTINUE; + } + if(musearch.GetDeviation() > cuts[2]) + { + musearch.SetNoMuon(); + return kCONTINUE; + } } @@ -481,23 +476,14 @@ if(!fDisablePreCuts) { - if(fMuonSize < 2000. || fArcPhi < 150.) - return kCONTINUE; - } - - fArcWid = CalcWid(geom,evt,musearch); - - // Pre Cuts 3 - if(!fDisablePreCuts) - { - if(fArcWid < 0.16) - return kCONTINUE; - } - + if(fMuonSize < cuts[3] || fArcPhi < cuts[4]) + { + musearch.SetNoMuon(); + return kCONTINUE; + } + } + + fArcWidth = CalcWidth(geom,evt,musearch); + SetReadyToSave(); - - if(!fKeepHist){ - fHistPhi.Reset(); - fHistWid.Reset(); - } return kTRUE; @@ -508,16 +494,16 @@ *fLog << all; *fLog << "Muon Parameters (" << GetName() << ")" << endl; - *fLog << " - Arc Len. [deg.] = " << fArcLen << endl; - *fLog << " - Arc Phi [deg.] = " << fArcPhi << endl; - *fLog << " - Arc Wid. [deg.] = " << fArcWid << endl; - *fLog << " - Chi Arc Phi [x2/ndf]= " << fChiArcPhi << endl; - *fLog << " - Chi Arc Wid [x2/ndf]= " << fChiArcWid << endl; - *fLog << " - Est. I. P. [m] = " << fEstImpact << endl; - *fLog << " - Size of muon = " << fMuonSize << endl; - *fLog << " - Peak Phi [deg.] = " << fPeakPhi << endl; - *fLog << " - UseUnmap = " << fUseUnmap << endl; -} - - - - + *fLog << " - Arc Length [deg.] = " << fArcLength << endl; + *fLog << " - Arc Phi [deg.] = " << fArcPhi << endl; + *fLog << " - Arc Width [deg.] = " << fArcWidth << endl; + *fLog << " - Chi Arc Phi [x2/ndf]= " << fChiArcPhi << endl; + *fLog << " - Chi Arc Width [x2/ndf]= " << fChiArcWidth << endl; + *fLog << " - Est. I. P. [m] = " << fEstImpact << endl; + *fLog << " - Size of muon = " << fMuonSize << endl; + *fLog << " - Peak Phi [deg.] = " << fPeakPhi << endl; + *fLog << " - UseUnmap = " << fUseUnmap << endl; +} + + + + Index: /trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.h =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.h (revision 5209) +++ /trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.h (revision 5210) @@ -17,60 +17,57 @@ { private: - Float_t fArcLen; // An arc length of a muon along the arc [deg.] - Float_t fArcPhi; // A opening angle of a muon arc [deg.] - Float_t fArcWid; // A width of a muon [deg.] (1 sigma of gaussian fit) - Float_t fChiArcPhi; // A chisquare value of the cosine fit for arc phi - Float_t fChiArcWid; // A chisquare value of the cosine fit for arc wid - Float_t fMuonSize; // A SIZE of muon which is defined as a SIZE around the estimated circle - Float_t fEstImpact; // An estimated impact parameter from the photon distribution along the arc image - Bool_t fKeepHist; // A flag to keep histgram - Float_t fMargin; // margin to evaluate muons [mm]. The defaut value is 60 mm, corresponding to 0.2 deg. This value can be changed by using the function of SetMargin - Bool_t fUseUnmap; // This is a flag to know the Unmapped pixels are used. Refer to the class of MImgCleanStd - Float_t fPeakPhi; // The angle which indicates the peak position in the estimated circle - Float_t fArcPhiThres; // The threshold value to define arc phi - Float_t fArcWidThres; // The threshold value to define arc width - Int_t fArcPhiBinNum; // The bin number for the histogram of arc phi. You may change this value. However, if you change this, YOU ALSO HAVE TO CHANGE THE THRESHOLD VALUE TO GET ARC LENGTH. - Int_t fArcWidBinNum; // The bin number for the histogram of arc wid - Float_t fArcPhiHistStartVal; // The starting value for the histogram of arc phi - Float_t fArcPhiHistEndVal; // The end value for the histogram of arc phi - Float_t fArcWidHistStartVal; // The starting value for the histogram of arc width - Float_t fArcWidHistEndVal; // The end value for the histogram of arc width - Float_t fEnableImpactCalc; // If true, the impact calculation will be done, which will consume a lot of time. - Float_t fDisablePreCuts; // If true, the pre cuts to select muons for the calibration will be disabled. - Float_t fUseCleanForWid; // If true, the only signal after image cleaningwill be used for the histogram of arc width + Float_t fArcLength; // An arc length of a muon along the arc [deg.] + Float_t fArcPhi; // A opening angle of a muon arc [deg.] + Float_t fArcWidth; // A width of a muon [deg.] (1 sigma of gaussian fit) + Float_t fChiArcPhi; // A chisquare value of the cosine fit for arc phi + Float_t fChiArcWidth; // A chisquare value of the cosine fit for arc wid + Float_t fMuonSize; // A SIZE of muon which is defined as a SIZE around the estimated circle + Float_t fEstImpact; // An estimated impact parameter from the photon distribution along the arc image + Float_t fMargin; // margin to evaluate muons [mm]. The defaut value is 60 mm, corresponding to 0.2 deg. This value can be changed by using the function of SetMargin + Bool_t fUseUnmap; // This is a flag to know the Unmapped pixels are used. Refer to the class of MImgCleanStd + Float_t fPeakPhi; // The angle which indicates the peak position in the estimated circle + Float_t fArcPhiThres; // The threshold value to define arc phi + Float_t fArcWidthThres; // The threshold value to define arc width + Int_t fArcPhiBinNum; // The bin number for the histogram of arc phi. You may change this value. However, if you change this, YOU ALSO HAVE TO CHANGE THE THRESHOLD VALUE TO GET ARC LENGTH. + Int_t fArcWidthBinNum; // The bin number for the histogram of arc wid + Float_t fArcPhiHistStartVal; // The starting value for the histogram of arc phi + Float_t fArcPhiHistEndVal; // The end value for the histogram of arc phi + Float_t fArcWidthHistStartVal; // The starting value for the histogram of arc width + Float_t fArcWidthHistEndVal; // The end value for the histogram of arc width + Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which will consume a lot of time. + Bool_t fDisablePreCuts; // If true, the pre cuts to select muons for the calibration will be disabled. + Bool_t fUseCleanForWidth; // If true, the only signal after image cleaningwill be used for the histogram of arc width - TH1F fHistPhi; // Histogram of photon distribution along the arc. - TH1F fHistWid; // Histogram of radial photon distribution of the arc. + TH1F *fHistPhi; // Histogram of photon distribution along the arc. + TH1F *fHistWidth; // Histogram of radial photon distribution of the arc. public: MMuonCalibPar(const char *name=NULL, const char *title=NULL); + ~MMuonCalibPar(); void Reset(); - Float_t GetArcLen() const { return fArcLen; } - Float_t GetArcPhi() const { return fArcPhi; } - Float_t GetArcWid() const { return fArcWid; } - Float_t GetChiArcPhi() const { return fChiArcPhi; } - Float_t GetChiArcWid() const { return fChiArcWid; } - Float_t GetMargin() const { return fMargin; } - Float_t GetMuonSize() const { return fMuonSize; } - Float_t GetEstImpact() const { return fEstImpact; } - Bool_t IsUseUnmap() const { return fUseUnmap; } - Float_t GetPeakPhi() const { return fPeakPhi; } - Float_t GetArcPhiThres() const { return fArcPhiThres; } - Float_t GetArcWidThres() const { return fArcWidThres; } - Float_t GetArcPhiBinNum() const { return fArcPhiBinNum; } - Float_t GetArcWidBinNum() const { return fArcWidBinNum; } - TH1F &GetHistPhi() { return fHistPhi; } - TH1F &GetHistWid() { return fHistWid; } + Float_t GetArcLength() const { return fArcLength; } + Float_t GetArcPhi() const { return fArcPhi; } + Float_t GetArcWidth() const { return fArcWidth; } + Float_t GetChiArcPhi() const { return fChiArcPhi; } + Float_t GetChiArcWidth() const { return fChiArcWidth; } + Float_t GetMargin() const { return fMargin; } + Float_t GetMuonSize() const { return fMuonSize; } + Float_t GetEstImpact() const { return fEstImpact; } + Bool_t IsUseUnmap() const { return fUseUnmap; } + Float_t GetPeakPhi() const { return fPeakPhi; } + Float_t GetArcPhiThres() const { return fArcPhiThres; } + Float_t GetArcWidthThres() const { return fArcWidthThres; } + Float_t GetArcPhiBinNum() const { return fArcPhiBinNum; } + Float_t GetArcWidthBinNum() const { return fArcWidthBinNum; } + TH1F *GetHistPhi() { return fHistPhi; } + TH1F *GetHistWidth() { return fHistWidth; } - void SetMargin(Float_t margin) { fMargin = margin; } - void SetArcPhiThres(Float_t thres) { fArcPhiThres = thres; } - void SetArcWidThres(Float_t thres) { fArcWidThres = thres; } - void SetArcPhiBinNum(Int_t num) { fArcPhiBinNum = num; } - void SetArcWidBinNum(Int_t num) { fArcWidBinNum = num; } - - void KeepHist() { fKeepHist = kTRUE; } - Bool_t IsKeepHist() { return fKeepHist; } + void SetMargin(Float_t margin) { fMargin = margin; } + void SetArcPhiThres(Float_t thres) { fArcPhiThres = thres; } + void SetArcWidthThres(Float_t thres) { fArcWidthThres = thres; } + void SetArcPhiBinNum(Int_t num) { fArcPhiBinNum = num; } + void SetArcWidthBinNum(Int_t num) { fArcWidthBinNum = num; } void EnableImpactCalc() { fEnableImpactCalc = kTRUE; } @@ -78,6 +75,6 @@ void DisablePreCuts() { fDisablePreCuts = kTRUE; } Bool_t IsDisablePreCuts() { return fDisablePreCuts; } - void UseCleanForWid() { fUseCleanForWid = kTRUE; } - Bool_t IsUseCleanForWid() { return fUseCleanForWid; } + void UseCleanForWidth() { fUseCleanForWidth = kTRUE; } + Bool_t IsUseCleanForWidth() { return fUseCleanForWidth; } void Print(Option_t *opt=NULL) const; @@ -89,10 +86,10 @@ void CalcImpact(const MGeomCam &geom, const MMuonSearchPar &musearch, Int_t effbinnum, Float_t startfitval, Float_t endfitval); - Float_t CalcWid(const MGeomCam &geom, const MCerPhotEvt &evt, + Float_t CalcWidth(const MGeomCam &geom, const MCerPhotEvt &evt, const MMuonSearchPar &musearch); Int_t Calc(const MGeomCam &geom, const MCerPhotEvt &evt, - const MMuonSearchPar &musearch); + MMuonSearchPar &musearch, const Float_t *cuts); - ClassDef(MMuonCalibPar, 3) // Container to hold muon calibration parameters + ClassDef(MMuonCalibPar, 1) // Container to hold muon calibration parameters }; Index: /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc (revision 5210) +++ /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc (revision 5210) @@ -0,0 +1,226 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Keiichi Mase 10/2004 +! Markus Meyer 10/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MMuonCalibParCalc +// +// Task to calculate the muon parameters +// +// This class allows you to get more muon information especially usefule for +// the calibration of our telescope. This class store the information into the +// container of MMuonCalibPar. The actual calculation is done in the class of +// MMuonCalibPar. Please see the manual. +// +// In order to make this class work, we need the information of the arc +// center and the radius. Therefore, we need to use the task of +// MMuonSearchParCalc. +// +// You can use this class such as the followings; +// +// MTaskList tlist; +// MMuonSearchParCalc musearchcalc; +// MMuonCalibParCalc mucalibcalc; +// tlist.AddToList(&musearchcalc); +// tlist.AddToList(&mucalibcalc);. +// +// You may change the allowed region to estimate muon parameters such as +// Muon SIZE and ARC LENGTH. The default value is 60 mm (0.2 deg.). If the +// estimated radius of the arc is 1.0 degree, we take the photons in the +// radius range from 0.8 to 1.2 degrees. You can change this value such as +// the followings; +// +// MParList plist; +// MMuonCalibPar muparcalib; +// plist.AddToList(&muparcalib);. +// muparcalib.SetMargin(60.); +// +// You can retrieve the histogram (TH1F) using the function of GetHistPhi() +// (also GetHistWid()). Therefore, you can draw the histogram such as +// +// muparcalib.GetHistPhi().Draw();. +// +// In order to use another information of muons such as the center position +// of the estimated circle, the radius of the circle. Use the infomation +// stored in MMuonSearchPar. +// +// +// For the faster computation, by default, the calculation of impact +// parameter is suppressed. If you want to calculate the impact parameter +// from the muon image, you can use the function of EnableImpactCalc(), +// namely; +// +// muparcalib.EnableImpactCalc();. +// +// In addition, for the faster comutation, pre cuts to select the candidates +// of muons for the calibration is done. You can set the values using the +// function of SetPreCuts. This function takes 5 variables. They correspond +// to the cur for the Arc Radius (low and high), the deviation of the fit +// (high), the Muon Size (low) and Arc Phi (low). You can set them such as +// +// mucalibcalc.SetPreCuts(180., 400., 50., 2000., 180.); +// +// If you want to disable the pre cuts, you can disable it by using the +// function of DisablePreCuts(), namely; +// +// muparcalib.DisablePreCuts();. +// +// +// ### TODO ### +// Up to now, in the histogram the error of the signal is estimated from +// the signal using a rough conversion factor and a F-factor and this values +// are global for all pixels. This is not the case for the real data. This +// value should be taken from some containers. In addition, the error of +// the pedestal is not taken into accout. The error treatment should be +// done correctly. +// +// Because of this rough treatment of the error, it may be better to use +// the image after the cleaning for the evaluation of the Arc Width. There +// is a possibility to use the image after the cleaning for the Arc Width. +// You can use the function of UseCleanForWid(). However, this treatment +// should be temporal and if the problem of the error is fixed, this function +// will be removed. +// +// +// Input Containers: +// [MGeomCam] +// [MCerPhotEvt] +// [MMuonSearchPar] +// +// Output Containers: +// [MMuonCalibPar] +// +////////////////////////////////////////////////////////////////////////////// + +#include "MMuonCalibParCalc.h" + +#include + +#include "MParList.h" + +#include "MGeomCam.h" +#include "MSrcPosCam.h" +#include "MCerPhotEvt.h" +#include "MMuonSearchPar.h" +#include "MMuonCalibPar.h" +#include "MLog.h" +#include "MLogManip.h" + +using namespace std; + +ClassImp(MMuonCalibParCalc); + +static const TString gsDefName = "MMuonCalibParCalc"; +static const TString gsDefTitle = "Calculate new image parameters"; + +// ------------------------------------------------------------------------- +// +// Default constructor. +// +MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title) + : fCerPhotName("MCerPhotEvt") +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + fPreCuts[0] = 180.; + fPreCuts[1] = 400.; + fPreCuts[2] = 50.; + fPreCuts[3] = 2000.; + fPreCuts[4] = 150.; +} + +// ------------------------------------------------------------------------- +// +Int_t MMuonCalibParCalc::PreProcess(MParList *pList) +{ + fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fCerPhotName), "MCerPhotEvt"); + if (!fCerPhotEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fGeomCam) + { + *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } + + fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar", "MMuonCalibPar"); + if (!fMuonCalibPar) + { + *fLog << dbginf << "MMuonCalibPar missing in Parameter List... aborting." << endl; + return kFALSE; + } + + fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", "MMuonSearchPar"); + if (!fMuonSearchPar) + { + *fLog << dbginf << "MMuonSearchPar missing in Parameter List... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// ------------------------------------------------------------------------- +// +Int_t MMuonCalibParCalc::Process() +{ + + if(!fMuonCalibPar->Calc(*fGeomCam, *fCerPhotEvt, *fMuonSearchPar, fPreCuts)) + return kCONTINUE; + + return kTRUE; +} + +void MMuonCalibParCalc::SetMargin(Float_t margin) +{ + fMuonCalibPar->SetMargin(margin); +} + +void MMuonCalibParCalc::EnableImpactCalc() +{ + fMuonCalibPar->EnableImpactCalc(); +} + +void MMuonCalibParCalc::DisablePreCuts() +{ + fMuonCalibPar->DisablePreCuts(); +} + +void MMuonCalibParCalc::SetPreCuts +(Float_t radcutlow, Float_t radcuthigh, Float_t devcuthigh, + Float_t musizecutlow, Float_t arcphicutlow) +{ + fPreCuts[0] = radcutlow; + fPreCuts[1] = radcuthigh; + fPreCuts[2] = devcuthigh; + fPreCuts[3] = musizecutlow; + fPreCuts[4] = arcphicutlow; +} + Index: /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h (revision 5210) +++ /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h (revision 5210) @@ -0,0 +1,44 @@ +#ifndef MARS_MMuonCalibParCalc +#define MARS_MMuonCalibParCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MMuonSearchPar; +class MMuonCalibPar; +class MSrcPosCam; +class MGeomCam; +class MCerPhotEvt; + +class MMuonCalibParCalc : public MTask +{ +private: + MGeomCam *fGeomCam; + MCerPhotEvt *fCerPhotEvt; + MMuonCalibPar *fMuonCalibPar; + MMuonSearchPar *fMuonSearchPar; + + Float_t fPreCuts[5]; // The values for pre cuts. + + Int_t PreProcess(MParList *plist); + Int_t Process(); + + TString fCerPhotName; + +public: + MMuonCalibParCalc(const char *name=NULL, const char *title=NULL); + + void SetMargin(Float_t margin); + void EnableImpactCalc(); + void DisablePreCuts(); + void SetPreCuts(Float_t radcutlow, Float_t radcuthigh, Float_t devcuthigh, + Float_t musizecutlow, Float_t arcphicutlow); + + void SetNameCerPhotEvt(const char *name) { fCerPhotName = name; } + + ClassDef(MMuonCalibParCalc, 0) // task to calculate muon parameters +}; + +#endif + Index: /trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc (revision 5209) +++ /trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc (revision 5210) @@ -16,6 +16,6 @@ ! ! -! Author(s): Keiichi Mase 09/2004 -! Markus Meyer 09/2004 +! Author(s): Keiichi Mase 10/2004 +! Markus Meyer 10/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 @@ -105,8 +105,8 @@ void MMuonSearchPar::Reset() { - fRad = -1.; - fDev = -1.; - fCenX = 0.; - fCenY = 0.; + fRadius = -1.; + fDeviation = -1.; + fCenterX = 0.; + fCenterY = 0.; fNoMuon = kFALSE; } @@ -123,8 +123,8 @@ Float_t tmp_r = 300.; // assume that the temporal cherenkov angle is 1 deg. (300 mm). - a = tan(hillas.GetDelta()); - - dx = a/sqrt(tmp_r+a*a)/3.; - dy = -tmp_r/sqrt(1+a*a)/3.; + a = TMath::Tan(hillas.GetDelta()); + + dx = a/TMath::Sqrt(tmp_r+a*a)/3.; + dy = -tmp_r/TMath::Sqrt(1+a*a)/3.; *xtmp1 = hillas.GetMeanX() + dx; @@ -180,5 +180,5 @@ // RMS of the fit, changing the center position of the circle. // -void MMuonSearchPar::CalcMinimumDev +void MMuonSearchPar::CalcMinimumDeviation ( const MGeomCam &geom, const MCerPhotEvt &evt, Float_t x, Float_t y, Float_t xcog, Float_t ycog, Float_t sigma, Float_t *opt_rad, @@ -268,12 +268,12 @@ // find the best fit. - CalcMinimumDev(geom, evt, xtmp1,ytmp1,hillas.GetMeanX(),hillas.GetMeanY(), + CalcMinimumDeviation(geom, evt, xtmp1,ytmp1,hillas.GetMeanX(),hillas.GetMeanY(), dev, &opt_rad, &new_sigma, &newx, &newy); - fRad = opt_rad; - fDev = new_sigma; - - fCenX = newx; - fCenY = newy; + fRadius = opt_rad; + fDeviation = new_sigma; + + fCenterX = newx; + fCenterY = newy; SetReadyToSave(); @@ -284,8 +284,8 @@ *fLog << all; *fLog << "Muon Parameters (" << GetName() << ")" << endl; - *fLog << " - Est. Rad. [mm] = " << fRad << endl; - *fLog << " - Dev. [mm] = " << fDev << endl; - *fLog << " - Cen. Pos. X [mm] = " << fCenX << endl; - *fLog << " - Cen. Pos. Y [mm] = " << fCenY << endl; + *fLog << " - Est. Radius [mm] = " << fRadius << endl; + *fLog << " - Deviation [mm] = " << fDeviation << endl; + *fLog << " - Center Pos. X [mm] = " << fCenterX << endl; + *fLog << " - Center Pos. Y [mm] = " << fCenterY << endl; } @@ -294,10 +294,10 @@ *fLog << all; *fLog << "Muon Parameters (" << GetName() << ")" << endl; - *fLog << " - Est. Rad. [deg.] = " << fRad*geom.GetConvMm2Deg() << endl; - *fLog << " - Dev. [deg.] = " << fDev*geom.GetConvMm2Deg() << endl; - *fLog << " - Cen. Pos. X [deg.] = " << fCenX*geom.GetConvMm2Deg() << endl; - *fLog << " - Cen. Pos. Y [deg.] = " << fCenY*geom.GetConvMm2Deg() << endl; -} - - - + *fLog << " - Est. Radius [deg.] = " << fRadius*geom.GetConvMm2Deg() << endl; + *fLog << " - Deviation [deg.] = " << fDeviation*geom.GetConvMm2Deg() << endl; + *fLog << " - Center Pos. X [deg.] = " << fCenterX*geom.GetConvMm2Deg() << endl; + *fLog << " - Center Pos. Y [deg.] = " << fCenterY*geom.GetConvMm2Deg() << endl; +} + + + Index: /trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h (revision 5209) +++ /trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h (revision 5210) @@ -21,9 +21,9 @@ { private: - Float_t fRad; // An estimated radius of the muon ring [mm] - Float_t fDev; // The standard deviation from the estimated ring [mm] - Float_t fCenX; // An estimated center position in X of the muon ring [mm] - Float_t fCenY; // An estimated center position in Y of the muon ring [mm] - Bool_t fNoMuon; // if the radius is estimated above 600 mm (2 deg.), assumed it's not muon. + Float_t fRadius; // An estimated radius of the muon ring [mm] + Float_t fDeviation; // The standard deviation from the estimated ring [mm] + Float_t fCenterX; // An estimated center position in X of the muon ring [mm] + Float_t fCenterY; // An estimated center position in Y of the muon ring [mm] + Bool_t fNoMuon; // if the radius is estimated above 600 mm (2 deg.), assumed it's not muon. Later on, at the stage of MMuonCalibParCalc, this flag will be changed if the task judge the event as no muon. public: @@ -32,9 +32,10 @@ void Reset(); - Float_t GetRad() const { return fRad; } - Float_t GetDev() const { return fDev; } - Float_t GetCenX() const { return fCenX; } - Float_t GetCenY() const { return fCenY; } - Bool_t IsNoMuon() const { return fNoMuon; } + Float_t GetRadius() const { return fRadius; } + Float_t GetDeviation() const { return fDeviation; } + Float_t GetCenterX() const { return fCenterX; } + Float_t GetCenterY() const { return fCenterY; } + Bool_t IsNoMuon() const { return fNoMuon; } + void SetNoMuon() { fNoMuon = kTRUE; } void Print(Option_t *opt=NULL) const; @@ -44,5 +45,5 @@ Bool_t CalcRadius(const MGeomCam &geom, const MCerPhotEvt &evt, Float_t x, Float_t y, Float_t *r, Float_t *sigma); - void CalcMinimumDev(const MGeomCam &geom, const MCerPhotEvt &evt, + void CalcMinimumDeviation(const MGeomCam &geom, const MCerPhotEvt &evt, Float_t x, Float_t y, Float_t xcog, Float_t ycog, Float_t sigma, Float_t *opt_rad, Index: /trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.cc (revision 5210) +++ /trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.cc (revision 5210) @@ -0,0 +1,125 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Keiichi Mase 09/2004 +! Markus Meyer 09/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MMuonSearchParCalc +// +// Task to calculate the muon parameters +// +// This class search for muons and store the information into the container +// of MMuonSearchPar. The actual calculation is done in the class of +// MMuonSearchPar. Please see the manual. +// +// In order to use further information of muons such as the width of arcs, +// the arc length along it, the muons size. Use the infomation stored in +// MMuonCalibPar. The information will be available by using the task of +// MMuonCalibParCalc. +// +// +// Input Containers: +// [MGeomCam] +// [MHillas] +// [MCerPhotEvt] +// +// Output Containers: +// [MMuonSearchPar] +// +////////////////////////////////////////////////////////////////////////////// +#include "MMuonSearchParCalc.h" + +#include + +#include "MParList.h" +#include "MGeomCam.h" +#include "MCerPhotEvt.h" +#include "MMuonSearchPar.h" +#include "MLog.h" +#include "MLogManip.h" + +using namespace std; + +ClassImp(MMuonSearchParCalc); + +static const TString gsDefName = "MMuonSearchParCalc"; +static const TString gsDefTitle = "Calculate muon parameters"; + +// ------------------------------------------------------------------------- +// +// Default constructor. +// +MMuonSearchParCalc::MMuonSearchParCalc +(const char *mupar, const char *name, const char *title) + : fHillas(NULL), fMuonPar(NULL), fCerPhotName("MCerPhotEvt") +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + fMuparName = mupar; + fHillasInput = "MHillas"; +} + +// ------------------------------------------------------------------------- +// +Int_t MMuonSearchParCalc::PreProcess(MParList *pList) +{ + fHillas = (MHillas*)pList->FindObject(fHillasInput, "MHillas"); + if (!fHillas) + { + *fLog << err << dbginf << "MHillas not found... aborting." << endl; + return kFALSE; + } + + fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fCerPhotName), "MCerPhotEvt"); + if (!fCerPhotEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fGeomCam) + { + *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } + + fMuonPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", fMuparName); + if (!fMuonPar) + return kFALSE; + + return kTRUE; +} + +// ------------------------------------------------------------------------- +// +Int_t MMuonSearchParCalc::Process() +{ + + fMuonPar->Calc(*fGeomCam, *fCerPhotEvt, *fHillas); + + return kTRUE; +} + Index: /trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.h (revision 5210) +++ /trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.h (revision 5210) @@ -0,0 +1,40 @@ +#ifndef MARS_MMuonSearchParCalc +#define MARS_MMuonSearchParCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MHillas; +class MMuonSearchPar; +class MGeomCam; +class MCerPhotEvt; + +class MMuonSearchParCalc : public MTask +{ +private: + MGeomCam *fGeomCam; + MCerPhotEvt *fCerPhotEvt; + MHillas *fHillas; //! Pointer to the source independent hillas parameters + MMuonSearchPar *fMuonPar; //! Pointer to the output container for the new image parameters + + TString fMuparName; + TString fHillasInput; + TString fCerPhotName; + + Int_t PreProcess(MParList *plist); + Int_t Process(); + +public: + MMuonSearchParCalc(const char *mupar="MMuonSearchPar", + const char *name=NULL, const char *title=NULL); + + void SetInput(TString hilname) { fHillasInput = hilname; } + + void SetNameCerPhotEvt(const char *name) { fCerPhotName = name; } + + ClassDef(MMuonSearchParCalc, 0) // task to calculate muon parameters +}; + +#endif +