Index: trunk/MagicSoft/Mars/mmuon/MHMuonPar.cc =================================================================== --- trunk/MagicSoft/Mars/mmuon/MHMuonPar.cc (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MHMuonPar.cc (revision 6979) @@ -16,9 +16,8 @@ ! ! -! Author(s): Wolfgang Wittek, 03/2003 -! Author(s): Thomas Bretz, 04/2003 ! Author(s): Markus Meyer, 02/2005 -! -! Copyright: MAGIC Software Development, 2000-2004 +! Author(s): Thomas Bretz, 04/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -28,4 +27,5 @@ // // MHMuonPar +// // This class is a histogram class for displaying muonparameters. // The folowing histgrams will be plotted: @@ -37,11 +37,10 @@ // // Inputcontainer: -// MMuonSearchPar -// MMuonCalibPar +// - MGeomCam +// - MMuonSearchPar +// - MMuonCalibPar // //////////////////////////////////////////////////////////////////////////// #include "MHMuonPar.h" - -#include #include @@ -49,4 +48,5 @@ #include #include + #include "MLog.h" #include "MLogManip.h" @@ -56,6 +56,4 @@ #include "MParList.h" -#include "MHillas.h" -#include "MHMuonPar.h" #include "MMuonSearchPar.h" #include "MMuonCalibPar.h" @@ -69,6 +67,6 @@ // Setup histograms // -MHMuonPar::MHMuonPar(const char *name, const char *title):fGeomCam(NULL), fMuonSearchPar(NULL), - fMuonCalibPar(NULL) +MHMuonPar::MHMuonPar(const char *name, const char *title) : + fMuonSearchPar(NULL), fMuonCalibPar(NULL) { fName = name ? name : "MHMuonPar"; @@ -76,6 +74,6 @@ fHistRadius.SetName("Radius"); - fHistRadius.SetTitle("Radius"); - fHistRadius.SetXTitle("Radius[deg]"); + fHistRadius.SetTitle("Distribution of Radius'"); + fHistRadius.SetXTitle("R [\\circ]"); fHistRadius.SetYTitle("Counts"); fHistRadius.SetDirectory(NULL); @@ -84,6 +82,6 @@ fHistArcWidth.SetName("ArcWidth"); - fHistArcWidth.SetTitle("ArcWidth"); - fHistArcWidth.SetXTitle("ArcWidth[deg]"); + fHistArcWidth.SetTitle("Distribution of ArcWidth"); + fHistArcWidth.SetXTitle("W [\\circ]"); fHistArcWidth.SetYTitle("Counts"); fHistArcWidth.SetDirectory(NULL); @@ -91,8 +89,8 @@ fHistArcWidth.SetFillStyle(4000); - fHistBroad.SetName("Ringbroadening"); - fHistBroad.SetTitle("Ringbroadening"); - fHistBroad.SetXTitle("Radius[deg]"); - fHistBroad.SetYTitle("ArcWidth/Radius"); + fHistBroad.SetName("RingBroadening"); + fHistBroad.SetTitle("Profile of ArcWidth/Radius versus Radius"); + fHistBroad.SetXTitle("R [\\circ]"); + fHistBroad.SetYTitle("W/R [1]"); fHistBroad.SetDirectory(NULL); fHistBroad.UseCurrentStyle(); @@ -100,6 +98,7 @@ fHistSize.SetName("SizeVsRadius"); - fHistSize.SetXTitle("Radius"); - fHistSize.SetYTitle("MuonSize"); + fHistSize.SetTitle("Profile of Muon Size vs. Radius"); + fHistSize.SetXTitle("Rs [[\circ]"); + fHistSize.SetYTitle("S [phe]"); fHistSize.SetDirectory(NULL); fHistSize.UseCurrentStyle(); @@ -119,5 +118,4 @@ bins.SetEdges(20, 0.5, 1.5); bins.Apply(fHistSize); - } @@ -129,10 +127,12 @@ Bool_t MHMuonPar::SetupFill(const MParList *plist) { - fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam"); - if (!fGeomCam) + MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam"); + if (!geom) { *fLog << warn << "MGeomCam not found... abort." << endl; return kFALSE; } + fMm2Deg = geom->GetConvMm2Deg(); + fMuonSearchPar = (MMuonSearchPar*)plist->FindObject("MMuonSearchPar"); if (!fMuonSearchPar) @@ -148,13 +148,8 @@ } - fMm2Deg = fGeomCam->GetConvMm2Deg(); - - ApplyBinning(*plist, "Radius", &fHistRadius); - - ApplyBinning(*plist, "ArcWidth", &fHistArcWidth); - - ApplyBinning(*plist, "Ringbroadening", &fHistBroad); - - ApplyBinning(*plist, "SizeVsRadius", &fHistSize); + ApplyBinning(*plist, "Radius", &fHistRadius); + ApplyBinning(*plist, "ArcWidth", &fHistArcWidth); + ApplyBinning(*plist, "RingBroadening", &fHistBroad); + ApplyBinning(*plist, "SizeVsRadius", &fHistSize); return kTRUE; @@ -211,5 +206,5 @@ fHistBroad.Draw(); } - +/* TH1 *MHMuonPar::GetHistByName(const TString name) { @@ -220,3 +215,3 @@ return NULL; } - +*/ Index: trunk/MagicSoft/Mars/mmuon/MHMuonPar.h =================================================================== --- trunk/MagicSoft/Mars/mmuon/MHMuonPar.h (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MHMuonPar.h (revision 6979) @@ -12,5 +12,4 @@ #endif -class MHillas; class MMuonSearchPar; class MMuonCalibPar; @@ -20,17 +19,12 @@ { private: - TH1F fHistRadius; // + TH1F fHistRadius; // + TH1F fHistArcWidth; // - TH1F fHistArcWidth; // - - TProfile fHistBroad; // Area of used pixels - + TProfile fHistBroad; // Area of used pixels TProfile fHistSize; // [ratio] concentration ratio: sum of the two highest pixels / fSize - MGeomCam *fGeomCam; //! Camera geometry for plots (for the moment this is a feature for a loop only!) - - MMuonSearchPar *fMuonSearchPar;//! - - MMuonCalibPar *fMuonCalibPar;//! + MMuonSearchPar *fMuonSearchPar; //! + MMuonCalibPar *fMuonCalibPar; //! Float_t fMm2Deg; @@ -42,13 +36,10 @@ Bool_t Fill(const MParContainer *par, const Stat_t w=1); - TH1 *GetHistByName(const TString name); + //TH1 *GetHistByName(const TString name); - TH1F &GetHistRadius() { return fHistRadius; } - - TH1F &GetHistArcWidth() { return fHistArcWidth; } - - TProfile &GetHistBroad() { return fHistBroad; } - - TProfile &GetHistSize() { return fHistSize; } + const TH1F& GetHistRadius() const { return fHistRadius; } + const TH1F& GetHistArcWidth() const { return fHistArcWidth; } + const TProfile& GetHistBroad() const { return fHistBroad; } + const TProfile& GetHistSize() const { return fHistSize; } void Draw(Option_t *opt=""); Index: trunk/MagicSoft/Mars/mmuon/MHSingleMuon.cc =================================================================== --- trunk/MagicSoft/Mars/mmuon/MHSingleMuon.cc (revision 6979) +++ trunk/MagicSoft/Mars/mmuon/MHSingleMuon.cc (revision 6979) @@ -0,0 +1,407 @@ +/* ======================================================================== *\ +! +! * +! * 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): Markus Meyer, 02/2005 +! Author(s): Thomas Bretz, 04/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHSingleMuon +// +// This class is a histogram class for a displaying muonparameters. +// The folowing histgrams will be plotted: +// - Radius (TH1F) +// - ArcWidth (TH1F) +// +// Inputcontainer: +// - MGeomCam +// - MMuonSearchPar +// - MMuonCalibPar +// +//////////////////////////////////////////////////////////////////////////// +#include "MHSingleMuon.h" + +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MSignalCam.h" +#include "MSignalPix.h" + +#include "MMuonSetup.h" +#include "MMuonCalibPar.h" +#include "MMuonSearchPar.h" + +ClassImp(MHSingleMuon); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Setup histograms +// +MHSingleMuon::MHSingleMuon(const char *name, const char *title) : + fSignalCam(0), fMuonSearchPar(0), fGeomCam(0), fMargin(0) +{ + fName = name ? name : "MHSingleMuon"; + fTitle = title ? title : "Histograms of muon parameters"; + + fHistPhi.SetName("HistPhi"); + fHistPhi.SetTitle("HistPhi"); + fHistPhi.SetXTitle("\\phi [\\circ]"); + 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 [\\circ]"); + fHistWidth.SetYTitle("sum of ADC"); + fHistWidth.SetDirectory(NULL); + fHistWidth.SetFillStyle(4000); + fHistWidth.UseCurrentStyle(); + + MBinning bins; + bins.SetEdges(21, -180, 180); + bins.Apply(fHistPhi); + + bins.SetEdges(28, 0.3, 1.7); + bins.Apply(fHistWidth); +} + +// -------------------------------------------------------------------------- +// +// Setup the Binning for the histograms automatically if the correct +// instances of MBinning +// +Bool_t MHSingleMuon::SetupFill(const MParList *plist) +{ + fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam"); + if (!fGeomCam) + { + *fLog << warn << "MGeomCam not found... abort." << endl; + return kFALSE; + } + fMuonSearchPar = (MMuonSearchPar*)plist->FindObject("MMuonSearchPar"); + if (!fMuonSearchPar) + { + *fLog << warn << "MMuonSearchPar not found... abort." << endl; + return kFALSE; + } + fSignalCam = (MSignalCam*)plist->FindObject("MSignalCam"); + if (!fSignalCam) + { + *fLog << warn << "MSignalCam not found... abort." << endl; + return kFALSE; + } + + MMuonSetup *setup = (MMuonSetup*)plist->FindObject("MMuonSetup"); + if (!setup) + { + *fLog << warn << "MMuonSetup not found... abort." << endl; + return kFALSE; + } + fMargin = setup->GetMargin()/fGeomCam->GetConvMm2Deg(); + + ApplyBinning(*plist, "ArcPhi", &fHistPhi); + ApplyBinning(*plist, "MuonWidth", &fHistWidth); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the histograms with data from a MMuonCalibPar and +// MMuonSearchPar container. +// +Bool_t MHSingleMuon::Fill(const MParContainer *par, const Stat_t w) +{ + fHistPhi.Reset(); + fHistWidth.Reset(); + + const Int_t entries = fSignalCam->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; iGetRadius() + fMargin && + dist > fMuonSearchPar->GetRadius() - fMargin) + fHistPhi.Fill(TMath::ATan2(dx, dy)*TMath::RadToDeg(), pix.GetNumPhotons()); + + // use only the inner pixles. This is geometry dependent. This has to + // be fixed! + if(i>397) + continue; + + fHistWidth.Fill(dist*fGeomCam->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; i=edge; i--) + { + const Float_t val = h.GetBinContent(i%n + 1); + if (val=0 && last>=0; +} + +// -------------------------------------------------------------------------- +// +// 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. +// +Bool_t MHSingleMuon::CalcPhi(Double_t thres, Double_t &peakphi, Double_t &arcphi) const +{ + peakphi = 180.-fHistPhi.GetBinCenter(fHistPhi.GetMaximumBin()); + + // Now find the position at which the peak edges crosses the threshold + Int_t first, last; + + if (!FindRangeAboveThreshold(fHistPhi, thres, first, last)) + return kFALSE; + + //cout << "P: " << first << " " << last << " " << last-first << endl; + + const Float_t startfitval = fHistPhi.GetBinLowEdge(first+1); + const Float_t endfitval = fHistPhi.GetBinLowEdge(last); + + //Int_t effbinnum = TMath::Nint((endfitval-startfitval)/convbin2val); + + arcphi = lastlast) + if (m>n) // If maximum is on the right side of histogram + last = n; + else + first = 0; // If maximum is on the left side of histogram + + // Now get the fit range + const Float_t startfitval = fHistWidth.GetBinLowEdge(first+1); + const Float_t endfitval = fHistWidth.GetBinLowEdge(last); + + // Setup the function and perform the fit + TF1 f1("f1", "gaus", startfitval, endfitval); + + // Choose starting values as accurate as possible + f1.SetParameter(0, fHistWidth.Integral(first+1, last)); + f1.SetParameter(1, fHistWidth.GetBinCenter(m)); + f1.SetParameter(2, (endfitval-startfitval)/2); + + // options : N do not store the function, do not draw + // I use integral of function in bin rather than value at bin center + // R use the range specified in the function range + // Q quiet mode + fHistWidth.Fit(&f1, "NQR"); + + if (last-first>3) + chi = f1.GetChisquare()/(last-first-3); + + Double_t err; + gMinuit->GetParameter(2, width, err); // get the sigma value + + return kTRUE; +} + +/* +// -------------------------------------------------------------------------- +// +// 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); + } +} +*/ Index: trunk/MagicSoft/Mars/mmuon/MHSingleMuon.h =================================================================== --- trunk/MagicSoft/Mars/mmuon/MHSingleMuon.h (revision 6979) +++ trunk/MagicSoft/Mars/mmuon/MHSingleMuon.h (revision 6979) @@ -0,0 +1,47 @@ +#ifndef MARS_MHSingleMuon +#define MARS_MHSingleMuon + +#ifndef MARS_MH +#include "MH.h" +#endif +#ifndef ROOT_TH1 +#include +#endif + +class MSignalCam; +class MMuonSearchPar; +class MMuonCalibPar; +class MMuonSetup; +class MGeomCam; + +class MHSingleMuon : public MH +{ +private: + MSignalCam *fSignalCam; //! + MMuonSearchPar *fMuonSearchPar; //! + MGeomCam *fGeomCam; //! + + Double_t fMargin; //! + + TH1F fHistPhi; // Histogram of photon distribution along the arc. + TH1F fHistWidth; // Histogram of radial photon distribution of the arc. + + Bool_t FindRangeAboveThreshold(const TH1F &h, Float_t thres, Int_t &first, Int_t &last) const; + +public: + MHSingleMuon(const char *name=NULL, const char *title=NULL); + + Bool_t SetupFill(const MParList *plist); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + + Bool_t CalcPhi(Double_t, Double_t &, Double_t &) const; + Bool_t CalcWidth(Double_t, Double_t &, Double_t &); + + const TH1F &GetHistPhi() const { return fHistPhi; } + + ClassDef(MHSingleMuon, 1) // Histograms of new image parameters +}; + +#endif + + Index: trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc (revision 6979) @@ -17,7 +17,7 @@ ! ! Author(s): Keiichi Mase 10/2004 -! Markus Meyer 10/2004 +! Author(s): Markus Meyer 10/2004 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -30,26 +30,14 @@ // Storage Container for muon // -// This class holds some information for a calibraion using muons. Muons +// This class holds some information for a calibration using muons. Muons // are identified by using the class of the MMuonSearchParCalc. You can fill // these information by using the MMuonCalibParCalc. See also these class // manuals. // -// -// Input Containers: -// [MGeomCam] -// [MSignalCam] -// [MMuonSearchPar] -// ///////////////////////////////////////////////////////////////////////////// #include "MMuonCalibPar.h" -#include - #include "MLog.h" #include "MLogManip.h" -#include "MSignalCam.h" -#include "MSignalPix.h" -#include "MMuonSearchPar.h" -#include "MBinning.h" using namespace std; @@ -66,31 +54,5 @@ fTitle = title ? title : "Muon calibration parameters"; - 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(); - -} - -// -------------------------------------------------------------------------- -// -MMuonCalibPar::~MMuonCalibPar() -{ - delete fHistPhi; - delete fHistWidth; + Reset(); } @@ -99,16 +61,12 @@ void MMuonCalibPar::Reset() { - fArcLength = -1.; - fArcPhi = 0.; - fArcWidth = -1.; - fChiArcPhi = -1.; - fChiArcWidth = -1.; - fMuonSize = 0.; - fEstImpact = -1.; - fUseUnmap = kFALSE; - fPeakPhi = 0.; - - fHistPhi->Reset(); - fHistWidth->Reset(); + fArcLength = -1.; + fArcPhi = 0.; + fArcWidth = -1.; + fChiArcPhi = -1.; + fChiArcWidth = -1.; + fMuonSize = 0.; + fEstImpact = -1.; + fPeakPhi = 0.; } @@ -116,14 +74,13 @@ { *fLog << all; - *fLog << "Muon Parameters (" << GetName() << ")" << 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; + *fLog << "Muon Parameters (" << GetName() << ")" << 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 [phe] = " << fMuonSize << endl; + *fLog << " - Peak Phi [deg] = " << fPeakPhi << endl; } Index: trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.h =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.h (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.h (revision 6979) @@ -16,48 +16,45 @@ { private: - 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 - 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 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 + //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 public: - MMuonCalibPar(const char *name=NULL, const char *title=NULL); - ~MMuonCalibPar(); - - TH1F *fHistPhi; // Histogram of photon distribution along the arc. - TH1F *fHistWidth; // Histogram of radial photon distribution of the arc. - - void Reset(); - - 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 GetMuonSize() const { return fMuonSize; } - Float_t GetEstImpact() const { return fEstImpact; } - Bool_t IsUseUnmap() const { return fUseUnmap; } - Float_t GetPeakPhi() const { return fPeakPhi; } - TH1F *GetHistPhi() { return fHistPhi; } - TH1F *GetHistWidth() { return fHistWidth; } - - void SetArcLength(Float_t len) { fArcLength = len; } - void SetArcPhi(Float_t phi) { fArcPhi = phi; } - void SetArcWidth(Float_t wid) { fArcWidth = wid; } - void SetChiArcPhi(Float_t chi) { fChiArcPhi = chi; } - void SetChiArcWidth(Float_t chi) { fChiArcWidth = chi; } - void SetMuonSize(Float_t size) { fMuonSize = size; } - void SetEstImpact(Float_t impact) { fEstImpact = impact; } - void SetUseUnmap() { fUseUnmap = kTRUE; } - void SetPeakPhi(Float_t phi) { fPeakPhi = phi; } - - void Print(Option_t *opt=NULL) const; + MMuonCalibPar(const char *name=NULL, const char *title=NULL); + //~MMuonCalibPar(); - ClassDef(MMuonCalibPar, 1) // Container to hold muon calibration parameters + void Reset(); + + 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 GetMuonSize() const { return fMuonSize; } + Float_t GetEstImpact() const { return fEstImpact; } + //Bool_t IsUseUnmap() const { return fUseUnmap; } + Float_t GetPeakPhi() const { return fPeakPhi; } + // TH1F *GetHistPhi() { return fHistPhi; } + // TH1F *GetHistWidth() { return fHistWidth; } + + void SetArcLength(Float_t len) { fArcLength = len; } + void SetArcPhi(Float_t phi) { fArcPhi = phi; } + void SetArcWidth(Float_t wid) { fArcWidth = wid; } + void SetChiArcPhi(Float_t chi) { fChiArcPhi = chi; } + void SetChiArcWidth(Float_t chi) { fChiArcWidth = chi; } + void SetMuonSize(Float_t size) { fMuonSize = size; } + void SetEstImpact(Float_t impact) { fEstImpact = impact; } + //void SetUseUnmap() { fUseUnmap = kTRUE; } + void SetPeakPhi(Float_t phi) { fPeakPhi = phi; } + + void Print(Option_t *opt=NULL) const; + + ClassDef(MMuonCalibPar, 1) // Container to hold muon calibration parameters }; Index: trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc (revision 6979) @@ -67,10 +67,10 @@ // // -// 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();. +// // 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 @@ -98,16 +98,14 @@ // // Input Containers: -// [MGeomCam] -// [MSignalCam] -// [MMuonSearchPar] +// MGeomCam +// MSignalCam +// MMuonSearchPar // // Output Containers: -// [MMuonCalibPar] +// MMuonCalibPar // ////////////////////////////////////////////////////////////////////////////// #include "MMuonCalibParCalc.h" - -#include #include @@ -115,16 +113,18 @@ #include +#include "MLog.h" +#include "MLogManip.h" + #include "MParList.h" #include "MGeomCam.h" #include "MGeomPix.h" -#include "MSrcPosCam.h" + #include "MSignalCam.h" -#include "MMuonSearchPar.h" + #include "MMuonCalibPar.h" #include "MMuonSetup.h" -#include "MLog.h" -#include "MLogManip.h" -#include "MBinning.h" +#include "MMuonSearchPar.h" +#include "MHSingleMuon.h" using namespace std; @@ -140,9 +140,8 @@ // MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title) +// : fEnableImpactCalc(kFALSE) { fName = name ? name : gsDefName.Data(); fTitle = title ? title : gsDefTitle.Data(); - - fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped. } @@ -151,379 +150,35 @@ Int_t MMuonCalibParCalc::PreProcess(MParList *pList) { - fSignalCam = (MSignalCam*)pList->FindObject("MSignalCam"); - if (!fSignalCam) + fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fGeomCam) { - *fLog << dbginf << "MSignalCam not found... aborting." << endl; - return kFALSE; + *fLog << err << "MGeomCam not found... abort." << endl; + return kFALSE; } - - fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fGeomCam) + + fHist = (MHSingleMuon*)pList->FindObject("MHSingleMuon"); + if (!fHist) { - *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; - return kFALSE; + *fLog << err << "MHSingleMuon not found... abort." << endl; + return kFALSE; } - - fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar", "MMuonCalibPar"); - if (!fMuonCalibPar) - return kFALSE; - - fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", "MMuonSearchPar"); - if (!fMuonSearchPar) - return kFALSE; - - fMuonSetup = (MMuonSetup*)pList->FindCreateObj("MMuonSetup", "MMuonSetup"); - if (!fMuonSetup) - return kFALSE; - - return kTRUE; + + fMuonSetup = (MMuonSetup*)pList->FindObject("MMuonSetup"); + if (!fMuonSetup) + { + *fLog << err << "MMuonSetup not found... abort." << endl; + return kFALSE; + } + + fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar"); + if (!fMuonCalibPar) + return kFALSE; + + fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar"); + if (!fMuonSearchPar) + 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 = fMuonSetup->fArcPhiBinNum; - Int_t binnumwid = fMuonSetup->fArcWidthBinNum; - - // preparation for a histgram - MBinning binsphi; - binsphi.SetEdges(binnumphi, - fMuonSetup->fArcPhiHistStartVal, - fMuonSetup->fArcPhiHistEndVal); - binsphi.Apply(*(fMuonCalibPar->fHistPhi)); - - MBinning binswid; - binswid.SetEdges(binnumwid, - fMuonSetup->fArcWidthHistStartVal, - fMuonSetup->fArcWidthHistEndVal); - binswid.Apply(*(fMuonCalibPar->fHistWidth)); - - const Int_t entries = (*fSignalCam).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() + fMuonSetup->GetMargin() - && dist > (*fMuonSearchPar).GetRadius() - fMuonSetup->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(i>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 = fMuonSetup->GetArcPhiThres(); - Float_t startval = fMuonSetup->fArcPhiHistStartVal; - Float_t endval = fMuonSetup->fArcPhiHistEndVal; - Int_t binnum = fMuonSetup->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 = fMuonSetup->fArcWidthHistStartVal; - Float_t endval = fMuonSetup->fArcWidthHistEndVal; - Int_t binnum = fMuonSetup->fArcWidthBinNum; - Float_t thres = fMuonSetup->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); - - // options : N do not store the function, do not draw - // I use integral of function in bin rather than value at bin center - // R use the range specified in the function range - // Q quiet mode - 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() -{ - // initialization - (*fMuonCalibPar).Reset(); - - // Fill signals to histograms - FillHist(); - - // Calculation of Arc Phi etc... - CalcPhi(); - - // Calculation of Arc Width etc... - fMuonCalibPar->SetArcWidth(CalcWidth()); - - if(fMuonCalibPar->GetArcPhi()>160 && fMuonSearchPar->GetRadius()>170 && - fMuonSearchPar->GetRadius()<400 && fMuonSearchPar->GetDeviation()<50) - fMuonCalibPar->SetReadyToSave(); - - return kTRUE; -} // ------------------------------------------------------------------------- @@ -531,7 +186,36 @@ Int_t MMuonCalibParCalc::Process() { - - Calc(); - - return kTRUE; + // Calculation of Arc Phi etc... + const Float_t thresphi = fMuonSetup->GetThresholdArcPhi() /fGeomCam->GetConvMm2Deg(); + const Float_t threswidth = fMuonSetup->GetThresholdArcWidth()/fGeomCam->GetConvMm2Deg(); + + Double_t peakphi, arcphi; + Double_t width, chi; + + if (!fHist->CalcPhi(thresphi, peakphi, arcphi)) + return kTRUE; + + if (!fHist->CalcWidth(threswidth, width, chi)) + return kTRUE; + + // Get Muon Size + fMuonCalibPar->SetMuonSize(fHist->GetHistPhi().Integral()); + + // Calculate Arc Length + fMuonCalibPar->SetPeakPhi(peakphi); + fMuonCalibPar->SetArcPhi(arcphi); + + const Float_t conv = TMath::RadToDeg()*fGeomCam->GetConvMm2Deg(); + fMuonCalibPar->SetArcLength(fMuonCalibPar->GetArcPhi()*fMuonSearchPar->GetRadius()*conv); + + // Calculation of Arc Width etc... + fMuonCalibPar->SetChiArcWidth(chi); + fMuonCalibPar->SetArcWidth(width); + + // Check if this is a 'Write-Out' candidate + if (fMuonCalibPar->GetArcPhi()>160 && fMuonSearchPar->GetRadius()<400 && + fMuonSearchPar->GetDeviation()<50 && fMuonSearchPar->GetRadius()>170) + fMuonCalibPar->SetReadyToSave(); + + return kTRUE; } Index: trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h (revision 6979) @@ -8,33 +8,31 @@ class MMuonSearchPar; class MMuonCalibPar; -class MSrcPosCam; class MGeomCam; -class MSignalCam; class MMuonSetup; +class MHSingleMuon; class MMuonCalibParCalc : public MTask { private: - MGeomCam *fGeomCam; - MSignalCam *fSignalCam; - MMuonCalibPar *fMuonCalibPar; - MMuonSearchPar *fMuonSearchPar; - MMuonSetup *fMuonSetup; + MGeomCam *fGeomCam; //! + MMuonCalibPar *fMuonCalibPar; //! + MMuonSearchPar *fMuonSearchPar; //! + MMuonSetup *fMuonSetup; //! + MHSingleMuon *fHist; //! - Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time. + //Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time. - Int_t PreProcess(MParList *plist); - Int_t Process(); + Int_t PreProcess(MParList *plist); + Int_t Process(); + + void FillHist(); + void CalcPhi(); + void CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval); + Float_t CalcWidth(); public: MMuonCalibParCalc(const char *name=NULL, const char *title=NULL); - void EnableImpactCalc() { fEnableImpactCalc = kTRUE; } - - void FillHist(); - void CalcPhi(); - void CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval); - Float_t CalcWidth(); - Int_t Calc(); + //void EnableImpactCalc(Bool_t b=kTRUE) { fEnableImpactCalc = b; } ClassDef(MMuonCalibParCalc, 0) // task to calculate muon parameters Index: trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc (revision 6979) @@ -17,7 +17,7 @@ ! ! Author(s): Keiichi Mase 10/2004 -! Markus Meyer 10/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 +! Author(s): Markus Meyer 10/2004 +! +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -40,5 +40,4 @@ // MMuonCalibPar. The information will be available by using the task of // MMuonCalibParCalc. -// // // --- How to search muons --- @@ -47,40 +46,39 @@ // // 1. A temporal center position of a circle is determined by using -// the Hillas parameters. Assumed that the center position will be on the -// line which is perpendicular to the longitudinal image axis and the -// distance from the gravity center of the image to the center position of -// a ring is approximately 1 deg. (corresponding to the Cherenkov angle.). -// Therefore, we will have two candidates of the center positions. +// the Hillas parameters. Assumed that the center position will be on the +// line which is perpendicular to the longitudinal image axis and the +// distance from the gravity center of the image to the center position of +// a ring is approximately 1 deg. (corresponding to the Cherenkov angle.). +// Therefore, we will have two candidates of the center positions. +// // 2. Find the ring radius which gives the minimum RMS between the camera -// images and the estimated circle. +// images and the estimated circle. +// // 3. Select one temporal position which gives smaller RMS as a true temporal -// center position. +// center position. +// // 4. Changing the center position of a circle on the camera plane from the -// determined temporal center position, find the position which gives the -// minimum RMS of the fit. -// -// -// --- Remark --- -// This method to search for muons is not fully optimized yet. However, -// it is good idea to use the temporal estimated center position from -// the Hillas parameters in order to reduce time to search. In addition, -// This method is faster than the MINUIT. +// determined temporal center position, find the position which gives the +// minimum RMS of the fit. // // // Input Containers: -// [MGeomCam] -// [MHillas] -// [MSignalCam] +// MGeomCam +// MHillas +// MSignalCam // ///////////////////////////////////////////////////////////////////////////// #include "MMuonSearchPar.h" -#include +#include #include "MLog.h" #include "MLogManip.h" + #include "MHillas.h" + #include "MGeomCam.h" #include "MGeomPix.h" + #include "MSignalPix.h" #include "MSignalCam.h" @@ -96,6 +94,6 @@ MMuonSearchPar::MMuonSearchPar(const char *name, const char *title) { - fName = name ? name : "MMuonSearchPar"; - fTitle = title ? title : "Muon search parameters"; + fName = name ? name : "MMuonSearchPar"; + fTitle = title ? title : "Muon search parameters"; } @@ -104,30 +102,19 @@ void MMuonSearchPar::Reset() { - fRadius = -1.; - fDeviation = -1.; - fCenterX = 0.; - fCenterY = 0.; -} - -// -------------------------------------------------------------------------- -// -// Get the tempolary center of a ring from the Hillas parameters. -// Two candidates of the position is returened. -// -void MMuonSearchPar::CalcTempCenter(const MHillas &hillas, - Float_t &xtmp1, Float_t &ytmp1, Float_t &xtmp2, Float_t &ytmp2) -{ - Float_t a,dx,dy; - Float_t tmp_r = 300.; // assume that the temporal cherenkov angle is 1 deg. (300 mm). - - 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; - ytmp1 = hillas.GetMeanY() + dy; - xtmp2 = hillas.GetMeanX() - dx; - ytmp2 = hillas.GetMeanY() - dy; + fRadius = -1; + fDeviation = -1; + fCenterX = 0; + fCenterY = 0; +} + +// -------------------------------------------------------------------------- +// +// This is a wrapper function to have direct access to the data members +// in the function calculating the minimization value. +// +void MMuonSearchPar::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + const MMuonSearchPar *optim = (MMuonSearchPar*)gMinuit->GetObjectFit(); + f = optim->Fcn(par); } @@ -137,39 +124,28 @@ // and its RMS for the input position. // -Bool_t MMuonSearchPar::CalcRadius(const MGeomCam &geom, const MSignalCam &evt, - Float_t x, Float_t y, Float_t &r, Float_t &sigma) -{ - Float_t mean_r=0., dev_r=0., sums=0., tmp=0.; - - const Int_t entries = evt.GetNumPixels(); - - for (Int_t i=0; i 360000.) - { - new_sigma=sigma; - opt_rad=rad_tmp; - newx=x; - newy=y; - break; - } - if(CalcRadius(geom,evt,x,y+delta,rad_tmp,sig_tmp)) - { - if(sig_tmp0 ? TMath::Sqrt(sigma) : 0; + + gMinuit = minsave; } @@ -250,36 +226,51 @@ // Calculation of muon parameters // -void MMuonSearchPar::Calc -(const MGeomCam &geom, const MSignalCam &evt, const MHillas &hillas) -{ - Reset(); - - Float_t xtmp1,xtmp2,ytmp1,ytmp2; - Float_t rad,dev,rad2,dev2; - Float_t opt_rad,new_sigma,newx,newy; - - // gets temporaly center - CalcTempCenter(hillas,xtmp1,ytmp1,xtmp2,ytmp2); - - // determine which position will be the true position. Here mainly - // the curvature of a muon arc is relied on. - CalcRadius(geom, evt, xtmp1,ytmp1,rad,dev); - CalcRadius(geom, evt, xtmp2,ytmp2,rad2,dev2); - if(dev2 -! Markus Meyer 09/2004 +! Author(s): Markus Meyer 09/2004 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -41,22 +41,22 @@ // // Input Containers: -// [MGeomCam] -// [MHillas] -// [MSignalCam] +// MGeomCam +// MHillas +// MSignalCam // // Output Containers: -// [MMuonSearchPar] +// MMuonSearchPar // ////////////////////////////////////////////////////////////////////////////// #include "MMuonSearchParCalc.h" -#include +#include "MParList.h" -#include "MParList.h" +#include "MLog.h" +#include "MLogManip.h" + #include "MGeomCam.h" #include "MSignalCam.h" #include "MMuonSearchPar.h" -#include "MLog.h" -#include "MLogManip.h" using namespace std; @@ -71,13 +71,9 @@ // Default constructor. // -MMuonSearchParCalc::MMuonSearchParCalc -(const char *mupar, const char *name, const char *title) +MMuonSearchParCalc::MMuonSearchParCalc(const char *name, const char *title) : fHillas(NULL), fMuonPar(NULL) { fName = name ? name : gsDefName.Data(); fTitle = title ? title : gsDefTitle.Data(); - - fMuparName = mupar; - fHillasInput = "MHillas"; } @@ -86,8 +82,8 @@ Int_t MMuonSearchParCalc::PreProcess(MParList *pList) { - fHillas = (MHillas*)pList->FindObject(fHillasInput, "MHillas"); + fHillas = (MHillas*)pList->FindObject("MHillas"); if (!fHillas) { - *fLog << err << fHillasInput << " [MHillas] not found... aborting." << endl; + *fLog << err << "MHillas not found... aborting." << endl; return kFALSE; } @@ -107,5 +103,5 @@ } - fMuonPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", fMuparName); + fMuonPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar"); if (!fMuonPar) return kFALSE; @@ -118,8 +114,5 @@ Int_t MMuonSearchParCalc::Process() { - - fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas); - - return kTRUE; + fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas); + return kTRUE; } - Index: trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.h (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.h (revision 6979) @@ -14,11 +14,8 @@ { private: - MGeomCam *fGeomCam; - MSignalCam *fSignalCam; - 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; + MGeomCam *fGeomCam; //! + MSignalCam *fSignalCam; //! + MHillas *fHillas; //! Pointer to the source independent hillas parameters + MMuonSearchPar *fMuonPar; //! Pointer to the output container for the new image parameters Int_t PreProcess(MParList *plist); @@ -26,8 +23,5 @@ public: - MMuonSearchParCalc(const char *mupar="MMuonSearchPar", - const char *name=NULL, const char *title=NULL); - - void SetInput(TString hilname) { fHillasInput = hilname; } + MMuonSearchParCalc(const char *name=NULL, const char *title=NULL); ClassDef(MMuonSearchParCalc, 0) // task to calculate muon parameters Index: trunk/MagicSoft/Mars/mmuon/MMuonSetup.cc =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonSetup.cc (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MMuonSetup.cc (revision 6979) @@ -16,8 +16,7 @@ ! ! -! Author(s): Keiichi Mase 10/2004 -! Markus Meyer 10/2004 +! Author(s): Markus Meyer 04/2005 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -28,9 +27,5 @@ // MMuonSetup // -// Storage Container for parameter for the muon analysis -// -// -// Input Containers: -// [MSignalCam] +// Storage Container for setup parameter for the muon analysis // ///////////////////////////////////////////////////////////////////////////// @@ -50,22 +45,36 @@ // MMuonSetup::MMuonSetup(const char *name, const char *title) + : fMargin(0.2), fThresholdArcPhi(0.1), fThresholdArcWidth(0.1) { fName = name ? name : "MMuonSetup"; - fTitle = title ? title : "Muon calibration parameters"; - - fMargin = 60.; // in mm - fArcPhiThres = 30.; - fArcWidthThres = 30.; - fArcPhiBinNum = 20; - fArcPhiHistStartVal = -180.; // deg. - fArcPhiHistEndVal = 180.; // deg. - fArcWidthBinNum = 28; - fArcWidthHistStartVal = 0.3; // deg. - fArcWidthHistEndVal = 1.7; // deg. + fTitle = title ? title : "Muon calibration setup parameters"; } // -------------------------------------------------------------------------- // -MMuonSetup::~MMuonSetup() +// Read the setup from a TEnv, eg: +// MMuonSetup.Margin: 0.2 +// MMuonSetup.ThresholdArcPhi: 0.1 +// MMuonSetup.ThresholdArcWidth: 0.1 +// +Int_t MMuonSetup::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { + Bool_t rc = kFALSE; + if (IsEnvDefined(env, prefix, "Margin", print)) + { + rc = kTRUE; + fMargin = GetEnvValue(env, prefix, "Margin", fMargin); + } + if (IsEnvDefined(env, prefix, "ThresholdArcPhi", print)) + { + rc = kTRUE; + fThresholdArcPhi = GetEnvValue(env, prefix, "ThresholdArcPhi", fThresholdArcPhi); + } + if (IsEnvDefined(env, prefix, "ThresholdArcWidth", print)) + { + rc = kTRUE; + fThresholdArcWidth = GetEnvValue(env, prefix, "ThresholdArcWidth", fThresholdArcWidth); + } + + return rc; } Index: trunk/MagicSoft/Mars/mmuon/MMuonSetup.h =================================================================== --- trunk/MagicSoft/Mars/mmuon/MMuonSetup.h (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MMuonSetup.h (revision 6979) @@ -9,33 +9,24 @@ { private: - 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 + Float_t fMargin; // [deg] margin to evaluate muons. The defaut value is 60 mm, corresponding to 0.2 deg. This value can be changed by using the function of SetMargin + Float_t fThresholdArcPhi; // [deg] The threshold value to define arc phi + Float_t fThresholdArcWidth; // [deg] The threshold value to define arc width + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); public: - MMuonSetup(const char *name=NULL, const char *title=NULL); - ~MMuonSetup(); - - 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 + MMuonSetup(const char *name=NULL, const char *title=NULL); - Float_t GetMargin() const { return fMargin; } - Float_t GetArcPhiThres() const { return fArcPhiThres; } - Float_t GetArcWidthThres() const { return fArcWidthThres; } - Float_t GetArcPhiBinNum() const { return fArcPhiBinNum; } - Float_t GetArcWidthBinNum() const { return fArcWidthBinNum; } - - 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; } - - ClassDef(MMuonSetup, 1) // Container to hold muon parameters for the whole sample + // Getter + Float_t GetMargin() const { return fMargin; } + Float_t GetThresholdArcPhi() const { return fThresholdArcPhi; } + Float_t GetThresholdArcWidth() const { return fThresholdArcWidth; } + + // Setter + void SetMargin(Float_t margin) { fMargin = margin; } + void SetThresholdArcPhi(Float_t thres) { fThresholdArcPhi = thres; } + void SetThresholdArcWidth(Float_t thres) { fThresholdArcWidth = thres; } + + ClassDef(MMuonSetup, 1) // Container to hold setup for muon analysis }; Index: trunk/MagicSoft/Mars/mmuon/Makefile =================================================================== --- trunk/MagicSoft/Mars/mmuon/Makefile (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/Makefile (revision 6979) @@ -30,4 +30,5 @@ MMuonCalibParCalc.cc \ MHMuonPar.cc \ + MHSingleMuon.cc \ MMuonSetup.cc Index: trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h (revision 6973) +++ trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h (revision 6979) @@ -6,16 +6,12 @@ #pragma link C++ class MMuonSearchPar+; -#pragma link C++ class MMuonSearchParCalc+; #pragma link C++ class MMuonCalibPar+; -#pragma link C++ class MMuonCalibParCalc+; -#pragma link C++ class MHMuonPar+; #pragma link C++ class MMuonSetup+; +#pragma link C++ class MMuonSearchParCalc+; +#pragma link C++ class MMuonCalibParCalc+; + +#pragma link C++ class MHMuonPar+; +#pragma link C++ class MHSingleMuon+; + #endif - - - - - - -