/* ======================================================================== *\ ! ! * ! * 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 useful for // the calibration of our telescope. This class store the information into the // container of MMuonCalibPar. // // 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; // // 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.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; // // mucalibcalc.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; // // mucalibcalc.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. // // // Input Containers: // [MGeomCam] // [MCerPhotEvt] // [MMuonSearchPar] // // Output Containers: // [MMuonCalibPar] // ////////////////////////////////////////////////////////////////////////////// #include "MMuonCalibParCalc.h" #include #include #include #include #include "MParList.h" #include "MGeomCam.h" #include "MGeomPix.h" #include "MSrcPosCam.h" #include "MCerPhotEvt.h" #include "MMuonSearchPar.h" #include "MMuonCalibPar.h" #include "MLog.h" #include "MLogManip.h" #include "MBinning.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) : fNameCerPhot("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.; 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. } // ------------------------------------------------------------------------- // Int_t MMuonCalibParCalc::PreProcess(MParList *pList) { 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(); } // use only the inner pixles. This is geometry dependent. This has to // be fixed! if(pix.GetPixId()>397) continue; 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; } // ------------------------------------------------------------------------- // Int_t MMuonCalibParCalc::Process() { fMuonCalibPar->SetMargin(fMargin); fMuonCalibPar->SetArcPhiThres(fArcPhiThres); fMuonCalibPar->SetArcWidthThres(fArcWidthThres); if(!Calc(fPreCuts)) return kCONTINUE; return kTRUE; } 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; }