Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 5331) +++ /trunk/MagicSoft/Mars/Changelog (revision 5332) @@ -23,5 +23,5 @@ 2004/11/01: Markus Meyer and Keiichi Mase - * mmuon/MMuonCalibParCalc.[h,cc], + * mmuon/MMuonCalibParCalc.[h,cc] - changed according to Thomas Bretz's suggestion. The main part of the calculation is moved from MMuonCalibPar to @@ -31,11 +31,10 @@ that we should get the Arc Width with uncleaned images. + * mmuon/MMuonCalibPar.[h,cc] + - changed according to Thomas Bretz's suggestion as written above. + * mmuon/MMuonSearchPar.[h,cc] - changed according to Thomas Bretz's suggestion. The pointer of the funciotn is changed to the reference. - - * mmuon/MMuonCalibPar.[h,cc], - - changed according to Thomas Bretz's suggestion as written above. - Index: /trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc (revision 5331) +++ /trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc (revision 5332) @@ -46,12 +46,6 @@ #include -#include -#include -#include - #include "MLog.h" #include "MLogManip.h" -#include "MGeomCam.h" -#include "MGeomPix.h" #include "MCerPhotEvt.h" #include "MCerPhotPix.h" @@ -91,8 +85,4 @@ fHistWidth->UseCurrentStyle(); - fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped. - fDisablePreCuts = kFALSE; // By default the pre cuts will be applied. - fUseCleanForWidth = kFALSE; // By default all the pixels will be used for the histogram of arc width. - fMargin = 60.; // in mm fArcPhiThres = 100.; @@ -102,6 +92,6 @@ fArcPhiHistEndVal = 180.; // deg. fArcWidthBinNum = 28; - fArcWidthHistStartVal = 0.3; // deg. - fArcWidthHistEndVal = 1.7; // deg. + fArcWidthHistStartVal = 0.3; // deg. + fArcWidthHistEndVal = 1.7; // deg. } @@ -127,366 +117,8 @@ fUseUnmap = kFALSE; fPeakPhi = 0.; + + fHistPhi->Reset(); + fHistWidth->Reset(); } - -// -------------------------------------------------------------------------- -// -// 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 -// UseCleanForWidth(). See the manual of MMuonCalibParCalc. -// -void MMuonCalibPar::FillHist -( const MGeomCam &geom, const MCerPhotEvt &evt, - const MMuonSearchPar &musearch) -{ - // preparation for a histgram - MBinning binsphi; - binsphi.SetEdges(fArcPhiBinNum, fArcPhiHistStartVal, fArcPhiHistEndVal); - binsphi.Apply(*fHistPhi); - - MBinning binswid; - 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.GetCenterX(); - const Float_t ceny = musearch.GetCenterY(); - - for (Int_t i=0; i0) - ang *= -1.0; - - // if the signal is not near the estimated circle, it is ignored. - if(dist < musearch.GetRadius() + fMargin && dist > musearch.GetRadius() - fMargin) - { - // check whether ummapped pixel is used or not. - // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information. - if(pix.IsPixelUnmapped()) - { - fUseUnmap = kTRUE; - continue; - } - fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons()); - fMuonSize += pix.GetNumPhotons(); - } - - if(pix.GetPixId()>397) - continue; // use only the inner pixles - - if(fUseCleanForWidth) - { - if(!pix.IsPixelUsed()) - continue; - } - - fHistWidth->Fill(dist*geom.GetConvMm2Deg(), pix.GetNumPhotons()); - } - - - // error estimation (temporaly) - // The error is estimated from the signal. In order to do so, we have to - // 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. - 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); - } - } -} - -// -------------------------------------------------------------------------- -// -// Photon distribution along the estimated circle is fitted with theoritical -// function in order to get some more information such as Arc Phi and Arc Length. -// -void MMuonCalibPar::CalcPhi -(const MGeomCam &geom, const MCerPhotEvt &evt, - const MMuonSearchPar &musearch) -{ - Float_t convbin2val = (fArcPhiHistEndVal-fArcPhiHistStartVal)/ - (Float_t)fArcPhiBinNum; - - // adjust the peak to 0 - Float_t maxval = 0.; - Int_t maxbin = 0; - maxval = fHistPhi->GetMaximum(); - maxbin = fHistPhi->GetMaximumBin(); - fPeakPhi = 180.-(Float_t)(maxbin-1)*convbin2val; - TArrayD tmp; - tmp.Set(fArcPhiBinNum+1); - for(Int_t i=1; iGetBinContent(i); - } - for(Int_t i=1; ifArcPhiBinNum) - { - id-=(fArcPhiBinNum); - } - if(id<=0) - { - id+=(fArcPhiBinNum); - } - fHistPhi->SetBinContent(i,tmp[id]); - } - maxbin = (Int_t)((Float_t)fArcPhiBinNum/2.)+1; - - // Determination of fitting region - // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore, - // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!! - Float_t startfitval = 0.; - Float_t endfitval = 0.; - Bool_t IsInMaxim = kFALSE; - Int_t effbinnum = 0; - for(Int_t i=1; iGetBinContent(i); - Float_t content_pre = fHistPhi->GetBinContent(i-1); - - if(content > fArcPhiThres && content_pre < fArcPhiThres) - { - startfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal; - } - if(i==maxbin) - IsInMaxim = kTRUE; - - if(content < fArcPhiThres && IsInMaxim == kTRUE) - { - endfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal; - break; - } - endfitval = fArcPhiHistEndVal; - } - - effbinnum = (Int_t)((endfitval-startfitval)/convbin2val); - - fArcPhi = effbinnum*convbin2val; - fArcLength = fArcPhi*3.1415926/180.*musearch.GetRadius()*geom.GetConvMm2Deg(); - - if(fEnableImpactCalc) - CalcImpact(geom, musearch, effbinnum, startfitval, endfitval); -} - -// -------------------------------------------------------------------------- -// -// An impact parameter is calculated by fitting the histogram of photon -// distribution along the circle with a theoritical model. -// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11. -// The function (6) is used here.) -// -// By default this calculation is suppressed because this calculation is -// very time consuming. If you want to calculate an impact parameter, -// you can call the function of EnableImpactCalc(). -// -void MMuonCalibPar::CalcImpact -( const MGeomCam &geom, const MMuonSearchPar &musearch, - Int_t effbinnum, Float_t startfitval, Float_t endfitval) -{ - // Fit the distribution with Vacanti function. The function is different - // for the impact parameter of inside or outside of our reflector. - // Then two different functions are applied to the photon distribution, - // and the one which give us smaller chisquare value is taken as a - // proper one. - Double_t val1,err1,val2,err2; - // impact parameter inside mirror radius (8.5m) - TString func1; - Float_t tmpval = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.; - tmpval = sin(2.*tmpval)*8.5; - func1 += "[0]*"; - func1 += tmpval; - func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))"; - - TF1 f1("f1",func1,startfitval,endfitval); - f1.SetParameters(2000,3,0); - f1.SetParLimits(1,0,8.5); - f1.SetParLimits(2,-180.,180.); - - fHistPhi->Fit("f1","RQEM"); - - Float_t chi1 = -1; - Float_t chi2 = -1; - if(effbinnum>3) - chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3)); - - gMinuit->GetParameter(1,val1,err1); // get the estimated IP - Float_t estip1 = val1; - - // impact parameter beyond mirror area (8.5m) - TString func2; - Float_t tmpval2 = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.; - tmpval2 = sin(2.*tmpval2)*8.5*2.; - func2 += "[0]*"; - func2 += tmpval2; - func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)"; - - TF1 f2("f2",func2,startfitval,endfitval); - f2.SetParameters(2000,20,0); - f2.SetParLimits(1,8.5,300.); - f2.SetParLimits(2,-180.,180.); - - fHistPhi->Fit("f2","RQEM+"); - - if(effbinnum>3) - chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3)); - - gMinuit->GetParameter(1,val2,err2); // get the estimated IP - Float_t estip2 = val2; - - if(effbinnum<=3) - { - fEstImpact = -1.; - } - if(chi2 > chi1) - { - fEstImpact = estip1; - fChiArcPhi = chi1; - } - else - { - fEstImpact = estip2; - fChiArcPhi = chi2; - } -} - -// -------------------------------------------------------------------------- -// -// Photon distribution of distance from the center of estimated ring is -// fitted in order to get some more information such as ARC WIDTH which -// can represent to the PSF of our reflector. -// -Float_t MMuonCalibPar::CalcWidth -(const MGeomCam &geom, const MCerPhotEvt &evt, - const MMuonSearchPar &musearch) -{ - Float_t convbin2val = (fArcWidthHistEndVal - fArcWidthHistStartVal) - /(Float_t)fArcWidthBinNum; - - // determination of fitting region - 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; iGetBinContent(i); - Float_t content_pre = fHistWidth->GetBinContent(i-1); - - if(content > fArcWidthThres) - effbinnum++; - - if(content > fArcWidthThres && content_pre < fArcWidthThres) - { - startfitval = (Float_t)(i-4)*convbin2val + fArcWidthHistStartVal; - if(startfitval<0.) startfitval = 0.; - } - if(i==maxbin) - IsInMaxim = kTRUE; - - if(content < fArcWidthThres && IsInMaxim == kTRUE) - { - endfitval = (Float_t)(i+2)*convbin2val + fArcWidthHistStartVal; - if(endfitval>180.) endfitval = 180.; - break; - } - endfitval = fArcWidthHistEndVal; - } - effbinnum = (Int_t)((endfitval-startfitval)/convbin2val); - - TF1 f1("f1","gaus",startfitval,endfitval); - - fHistWidth->Fit("f1","QR","",startfitval,endfitval); - - if(effbinnum>3) - fChiArcWidth = f1.GetChisquare()/((Float_t)(effbinnum-3)); - - Double_t val,err; - gMinuit->GetParameter(2,val,err); // get the sigma value - - return val; -} - -// -------------------------------------------------------------------------- -// -// Calculation of muon parameters -// -Int_t MMuonCalibPar::Calc -(const MGeomCam &geom, const MCerPhotEvt &evt, - MMuonSearchPar &musearch, const Float_t *cuts) -{ - // sanity check - if(evt.GetNumPixels() < 3) - return kCONTINUE; - - // If an event does not seem to be like muon, the calculation will be skipped. - if(musearch.IsNoMuon()) - return kCONTINUE; - - // Pre Cuts 1 - if(!fDisablePreCuts) - { - if(musearch.GetRadius() < cuts[0] || musearch.GetRadius() > cuts[1]) - { - musearch.SetNoMuon(); - return kCONTINUE; - } - if(musearch.GetDeviation() > cuts[2]) - { - musearch.SetNoMuon(); - return kCONTINUE; - } - } - - Reset(); - - FillHist(geom,evt,musearch); - - CalcPhi(geom,evt,musearch); - - // Pre Cuts 2 - if(!fDisablePreCuts) - { - if(fMuonSize < cuts[3] || fArcPhi < cuts[4]) - { - musearch.SetNoMuon(); - return kCONTINUE; - } - } - - fArcWidth = CalcWidth(geom,evt,musearch); - - SetReadyToSave(); - - return kTRUE; -} void MMuonCalibPar::Print(Option_t *) const Index: /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc (revision 5331) +++ /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc (revision 5332) @@ -30,8 +30,7 @@ // Task to calculate the muon parameters // -// This class allows you to get more muon information especially usefule for +// This class allows you to get more muon information especially useful 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. +// container of MMuonCalibPar. // // In order to make this class work, we need the information of the arc @@ -53,12 +52,12 @@ // the followings; // +// mucalibcalc.SetMargin(60.); +// +// You can retrieve the histogram (TH1F) using the function of GetHistPhi() +// (also GetHistWid()). Therefore, you can draw the histogram such as +// // 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();. // @@ -73,5 +72,5 @@ // namely; // -// muparcalib.EnableImpactCalc();. +// mucalibcalc.EnableImpactCalc();. // // In addition, for the faster comutation, pre cuts to select the candidates @@ -86,5 +85,5 @@ // function of DisablePreCuts(), namely; // -// muparcalib.DisablePreCuts();. +// mucalibcalc.DisablePreCuts();. // // @@ -97,11 +96,4 @@ // 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: @@ -119,7 +111,12 @@ #include +#include +#include +#include + #include "MParList.h" #include "MGeomCam.h" +#include "MGeomPix.h" #include "MSrcPosCam.h" #include "MCerPhotEvt.h" @@ -128,4 +125,5 @@ #include "MLog.h" #include "MLogManip.h" +#include "MBinning.h" using namespace std; @@ -141,5 +139,5 @@ // MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title) - : fCerPhotName("MCerPhotEvt") + : fNameCerPhot("MCerPhotEvt") { fName = name ? name : gsDefName.Data(); @@ -151,4 +149,11 @@ fPreCuts[3] = 2000.; fPreCuts[4] = 150.; + + fMargin = 60.; + fArcPhiThres = 100.; + fArcWidthThres = 100.; + + fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped. + fDisablePreCuts = kFALSE; // By default the pre cuts will be applied. } @@ -157,28 +162,406 @@ Int_t MMuonCalibParCalc::PreProcess(MParList *pList) { - fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fCerPhotName), "MCerPhotEvt"); - if (!fCerPhotEvt) - { - *fLog << err << fCerPhotName << " [MCerPhotEvt] not found... aborting." << endl; - return kFALSE; - } - - fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fGeomCam) - { - *fLog << err << "MGeomCam not found... aborting." << endl; - return kFALSE; - } - - fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar"); - if (!fMuonCalibPar) - return kFALSE; - - fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar"); - if (!fMuonSearchPar) - return kFALSE; - - return kTRUE; -} + fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fNameCerPhot), "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; +} + +// -------------------------------------------------------------------------- +// +// 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. +// +void MMuonCalibParCalc::FillHist() +{ + Float_t MuonSize = 0.; + + Int_t binnumphi = fMuonCalibPar->fArcPhiBinNum; + Int_t binnumwid = fMuonCalibPar->fArcWidthBinNum; + + // preparation for a histgram + MBinning binsphi; + binsphi.SetEdges(binnumphi, + fMuonCalibPar->fArcPhiHistStartVal, + fMuonCalibPar->fArcPhiHistEndVal); + binsphi.Apply(*(fMuonCalibPar->fHistPhi)); + + MBinning binswid; + binswid.SetEdges(binnumwid, + fMuonCalibPar->fArcWidthHistStartVal, + fMuonCalibPar->fArcWidthHistEndVal); + binswid.Apply(*(fMuonCalibPar->fHistWidth)); + + const Int_t entries = (*fCerPhotEvt).GetNumPixels(); + + // the position of the center of a muon ring + const Float_t cenx = (*fMuonSearchPar).GetCenterX(); + const Float_t ceny = (*fMuonSearchPar).GetCenterY(); + + for (Int_t i=0; i0) + ang *= -1.0; + + // if the signal is not near the estimated circle, it is ignored. + if(dist < (*fMuonSearchPar).GetRadius() + fMuonCalibPar->GetMargin() + && dist > (*fMuonSearchPar).GetRadius() - fMuonCalibPar->GetMargin()) + { + // check whether ummapped pixel is used or not. + // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information. + if(pix.IsPixelUnmapped()) + { + fMuonCalibPar->SetUseUnmap(); + continue; + } + fMuonCalibPar->fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons()); + MuonSize += pix.GetNumPhotons(); + } + + if(pix.GetPixId()>397) + continue; // use only the inner pixles + + fMuonCalibPar->fHistWidth + ->Fill(dist*(*fGeomCam).GetConvMm2Deg(), pix.GetNumPhotons()); + } + + fMuonCalibPar->SetMuonSize(MuonSize); + + // error estimation (temporaly) + // The error is estimated from the signal. In order to do so, we have to + // 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. + Double_t ADC2PhEl = 0.14; + Double_t Ffactor = 1.4; + for(Int_t i=0; ifHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor; + { + fMuonCalibPar->fHistPhi->SetBinError(i, rougherr); + } + } + for(Int_t i=0; ifHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor; + { + fMuonCalibPar->fHistWidth->SetBinError(i, rougherr); + } + } +} + +// -------------------------------------------------------------------------- +// +// Photon distribution along the estimated circle is fitted with theoritical +// function in order to get some more information such as Arc Phi and Arc +// Length. +// +void MMuonCalibParCalc::CalcPhi() +{ + Float_t thres = fMuonCalibPar->GetArcPhiThres(); + Float_t startval = fMuonCalibPar->fArcPhiHistStartVal; + Float_t endval = fMuonCalibPar->fArcPhiHistEndVal; + Int_t binnum = fMuonCalibPar->fArcPhiBinNum; + + Float_t convbin2val = (endval-startval)/(Float_t)binnum; + + // adjust the peak to 0 + Float_t maxval = 0.; + Int_t maxbin = 0; + maxval = fMuonCalibPar->fHistPhi->GetMaximum(); + maxbin = fMuonCalibPar->fHistPhi->GetMaximumBin(); + fMuonCalibPar->SetPeakPhi(180.-(Float_t)(maxbin-1.)*convbin2val); + TArrayD tmp; + tmp.Set(binnum+1); + for(Int_t i=1; ifHistPhi->GetBinContent(i); + } + for(Int_t i=1; ibinnum) + { + id-=(binnum); + } + if(id<=0) + { + id+=(binnum); + } + fMuonCalibPar->fHistPhi->SetBinContent(i,tmp[id]); + } + maxbin = (Int_t)((Float_t)binnum/2.)+1; + + // Determination of fitting region + // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore, + // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!! + Float_t startfitval = 0.; + Float_t endfitval = 0.; + Bool_t IsInMaxim = kFALSE; + Int_t effbinnum = 0; + for(Int_t i=1; ifHistPhi->GetBinContent(i); + Float_t content_pre = fMuonCalibPar->fHistPhi->GetBinContent(i-1); + + if(content > thres && content_pre < thres) + { + startfitval = (Float_t)(i-1)*convbin2val+startval; + } + if(i==maxbin) + IsInMaxim = kTRUE; + + if(content < thres && IsInMaxim == kTRUE) + { + endfitval = (Float_t)(i-1)*convbin2val+startval; + break; + } + endfitval = endval; + } + + effbinnum = (Int_t)((endfitval-startfitval)/convbin2val); + + fMuonCalibPar->SetArcPhi(endfitval-startfitval); + + fMuonCalibPar->SetArcLength( fMuonCalibPar->GetArcPhi()*TMath::DegToRad()*(*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()); + + if(fEnableImpactCalc) + CalcImpact(effbinnum, startfitval, endfitval); +} + +// -------------------------------------------------------------------------- +// +// An impact parameter is calculated by fitting the histogram of photon +// distribution along the circle with a theoritical model. +// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11. +// The function (6) is used here.) +// +// By default this calculation is suppressed because this calculation is +// very time consuming. If you want to calculate an impact parameter, +// you can call the function of EnableImpactCalc(). +// +void MMuonCalibParCalc::CalcImpact +(Int_t effbinnum, Float_t startfitval, Float_t endfitval) +{ + // Fit the distribution with Vacanti function. The function is different + // for the impact parameter of inside or outside of our reflector. + // Then two different functions are applied to the photon distribution, + // and the one which give us smaller chisquare value is taken as a + // proper one. + Double_t val1,err1,val2,err2; + // impact parameter inside mirror radius (8.5m) + TString func1; + Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad(); + tmpval = sin(2.*tmpval)*8.5; + func1 += "[0]*"; + func1 += tmpval; + func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))"; + + TF1 f1("f1",func1,startfitval,endfitval); + f1.SetParameters(2000,3,0); + f1.SetParLimits(1,0,8.5); + f1.SetParLimits(2,-180.,180.); + + fMuonCalibPar->fHistPhi->Fit("f1","RQEM"); + + Float_t chi1 = -1; + Float_t chi2 = -1; + if(effbinnum>3) + chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3)); + + gMinuit->GetParameter(1,val1,err1); // get the estimated IP + Float_t estip1 = val1; + + // impact parameter beyond mirror area (8.5m) + TString func2; + Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad(); + tmpval2 = sin(2.*tmpval2)*8.5*2.; + func2 += "[0]*"; + func2 += tmpval2; + func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)"; + + TF1 f2("f2",func2,startfitval,endfitval); + f2.SetParameters(2000,20,0); + f2.SetParLimits(1,8.5,300.); + f2.SetParLimits(2,-180.,180.); + + fMuonCalibPar->fHistPhi->Fit("f2","RQEM+"); + + if(effbinnum>3) + chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3)); + + gMinuit->GetParameter(1,val2,err2); // get the estimated IP + Float_t estip2 = val2; + + if(effbinnum<=3) + { + fMuonCalibPar->SetEstImpact(-1.); + } + if(chi2 > chi1) + { + fMuonCalibPar->SetEstImpact(estip1); + fMuonCalibPar->SetChiArcPhi(chi1); + } + else + { + fMuonCalibPar->SetEstImpact(estip2); + fMuonCalibPar->SetChiArcPhi(chi2); + } +} + +// -------------------------------------------------------------------------- +// +// Photon distribution of distance from the center of estimated ring is +// fitted in order to get some more information such as ARC WIDTH which +// can represent to the PSF of our reflector. +// +Float_t MMuonCalibParCalc::CalcWidth() +{ + Float_t startval = fMuonCalibPar->fArcWidthHistStartVal; + Float_t endval = fMuonCalibPar->fArcWidthHistEndVal; + Int_t binnum = fMuonCalibPar->fArcWidthBinNum; + Float_t thres = fMuonCalibPar->GetArcWidthThres(); + + Float_t convbin2val = (endval - startval) + /(Float_t)binnum; + + // determination of fitting region + Int_t maxbin = fMuonCalibPar->fHistWidth->GetMaximumBin(); + Float_t startfitval = 0.; + Float_t endfitval = 0.; + Bool_t IsInMaxim = kFALSE; + Int_t effbinnum = 0; + for(Int_t i=1; ifHistWidth->GetBinContent(i); + Float_t content_pre = fMuonCalibPar->fHistWidth->GetBinContent(i-1); + + if(content > thres) + effbinnum++; + + if(content > thres && content_pre < thres) + { + startfitval = (Float_t)(i-4)*convbin2val + startval; + if(startfitval<0.) startfitval = 0.; + } + if(i==maxbin) + IsInMaxim = kTRUE; + + if(content < thres && IsInMaxim == kTRUE) + { + endfitval = (Float_t)(i+2)*convbin2val + startval; + if(endfitval>180.) endfitval = 180.; + break; + } + endfitval = endval; + } + effbinnum = (Int_t)((endfitval-startfitval)/convbin2val); + + TF1 f1("f1","gaus",startfitval,endfitval); + + fMuonCalibPar->fHistWidth->Fit("f1","QR","",startfitval,endfitval); + + if(effbinnum>3) + fMuonCalibPar->SetChiArcWidth(f1.GetChisquare()/((Float_t)(effbinnum-3))); + + Double_t val,err; + gMinuit->GetParameter(2,val,err); // get the sigma value + + return val; +} + +// -------------------------------------------------------------------------- +// +// Calculation of muon parameters +// +Int_t MMuonCalibParCalc::Calc(const Float_t *cuts) +{ + // sanity check + if((*fCerPhotEvt).GetNumPixels() < 3) + return kCONTINUE; + + // If an event does not seem to be like muon, the calculation will be skipped. + if((*fMuonSearchPar).IsNoMuon()) + return kCONTINUE; + + // Pre Cuts 1 + if(!fDisablePreCuts) + { + if((*fMuonSearchPar).GetRadius() < cuts[0] || (*fMuonSearchPar).GetRadius() > cuts[1]) + { + (*fMuonSearchPar).SetNoMuon(); + return kCONTINUE; + } + if((*fMuonSearchPar).GetDeviation() > cuts[2]) + { + (*fMuonSearchPar).SetNoMuon(); + return kCONTINUE; + } + } + + // initialization + (*fMuonCalibPar).Reset(); + + // Fill signals to histograms + FillHist(); + + // Calculation of Arc Phi etc... + CalcPhi(); + + // Pre Cuts 2 + if(!fDisablePreCuts) + { + if(fMuonCalibPar->GetMuonSize() < cuts[3] + || fMuonCalibPar->GetArcPhi() < cuts[4]) + { + (*fMuonSearchPar).SetNoMuon(); + return kCONTINUE; + } + } + + // Calculation of Arc Width etc... + fMuonCalibPar->SetArcWidth(CalcWidth()); + + return kTRUE; +} + // ------------------------------------------------------------------------- @@ -186,24 +569,12 @@ Int_t MMuonCalibParCalc::Process() { - - if(!fMuonCalibPar->Calc(*fGeomCam, *fCerPhotEvt, *fMuonSearchPar, fPreCuts)) + fMuonCalibPar->SetMargin(fMargin); + fMuonCalibPar->SetArcPhiThres(fArcPhiThres); + fMuonCalibPar->SetArcWidthThres(fArcWidthThres); + + if(!Calc(fPreCuts)) return kCONTINUE; return kTRUE; -} - -void MMuonCalibParCalc::SetMargin(Float_t margin) -{ - fMuonCalibPar->SetMargin(margin); -} - -void MMuonCalibParCalc::EnableImpactCalc() -{ - fMuonCalibPar->EnableImpactCalc(); -} - -void MMuonCalibParCalc::DisablePreCuts() -{ - fMuonCalibPar->DisablePreCuts(); } Index: /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h (revision 5331) +++ /trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h (revision 5332) @@ -20,4 +20,10 @@ MMuonSearchPar *fMuonSearchPar; + 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 + Float_t fArcPhiThres; // The threshold value to define arc phi + Float_t fArcWidthThres; // The threshold value to define arc width + Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time. + Bool_t fDisablePreCuts; // If true, the pre cuts to select muons for the calibration will be disabled. + Float_t fPreCuts[5]; // The values for pre cuts. @@ -25,16 +31,24 @@ Int_t Process(); - TString fCerPhotName; + TString fNameCerPhot; public: MMuonCalibParCalc(const char *name=NULL, const char *title=NULL); - void SetMargin(Float_t margin); - void EnableImpactCalc(); - void DisablePreCuts(); + void SetMargin(Float_t margin) { fMargin = margin; } + void SetArcPhiThres(Float_t thres) { fArcPhiThres = thres; } + void SetArcWidthThres(Float_t thres) { fArcWidthThres = thres; } + void EnableImpactCalc() { fEnableImpactCalc = kTRUE; } + void DisablePreCuts() { fDisablePreCuts = kTRUE; } void SetPreCuts(Float_t radcutlow, Float_t radcuthigh, Float_t devcuthigh, Float_t musizecutlow, Float_t arcphicutlow); - void SetNameCerPhotEvt(const char *name) { fCerPhotName = name; } + void SetNameCerPhotEvt(const char *name) { fNameCerPhot = name; } + + void FillHist(); + void CalcPhi(); + void CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval); + Float_t CalcWidth(); + Int_t Calc(const Float_t *cuts); ClassDef(MMuonCalibParCalc, 0) // task to calculate muon parameters Index: /trunk/MagicSoft/Mars/mmuon/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mmuon/Makefile (revision 5331) +++ /trunk/MagicSoft/Mars/mmuon/Makefile (revision 5332) @@ -26,5 +26,7 @@ SRCFILES = MMuonSearchPar.cc \ - MMuonCalibPar.cc + MMuonSearchParCalc.cc \ + MMuonCalibPar.cc \ + MMuonCalibParCalc.cc ############################################################ @@ -37,2 +39,9 @@ mrproper: clean rmbak + + + + + + + Index: /trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h (revision 5331) +++ /trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h (revision 5332) @@ -6,5 +6,7 @@ #pragma link C++ class MMuonSearchPar+; +#pragma link C++ class MMuonSearchParCalc+; #pragma link C++ class MMuonCalibPar+; +#pragma link C++ class MMuonCalibParCalc+; #endif