Changeset 6979 for trunk/MagicSoft/Mars/mmuon
- Timestamp:
- 04/27/05 16:46:24 (20 years ago)
- Location:
- trunk/MagicSoft/Mars/mmuon
- Files:
-
- 2 added
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mmuon/MHMuonPar.cc
r6973 r6979 16 16 ! 17 17 ! 18 ! Author(s): Wolfgang Wittek, 03/2003 <mailto:wittek@mppmu.mpg.de>19 ! Author(s): Thomas Bretz, 04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>20 18 ! Author(s): Markus Meyer, 02/2005 <mailto:meyer@astro.uni-wuerzburg.de> 21 ! 22 ! Copyright: MAGIC Software Development, 2000-2004 19 ! Author(s): Thomas Bretz, 04/2005 <mailto:tbretz@astro.uni-wuerzburg.de> 20 ! 21 ! Copyright: MAGIC Software Development, 2000-2005 23 22 ! 24 23 ! … … 28 27 // 29 28 // MHMuonPar 29 // 30 30 // This class is a histogram class for displaying muonparameters. 31 31 // The folowing histgrams will be plotted: … … 37 37 // 38 38 // Inputcontainer: 39 // MMuonSearchPar 40 // MMuonCalibPar 39 // - MGeomCam 40 // - MMuonSearchPar 41 // - MMuonCalibPar 41 42 // 42 43 //////////////////////////////////////////////////////////////////////////// 43 44 #include "MHMuonPar.h" 44 45 #include <math.h>46 45 47 46 #include <TH1.h> … … 49 48 #include <TCanvas.h> 50 49 #include <TProfile.h> 50 51 51 #include "MLog.h" 52 52 #include "MLogManip.h" … … 56 56 #include "MParList.h" 57 57 58 #include "MHillas.h"59 #include "MHMuonPar.h"60 58 #include "MMuonSearchPar.h" 61 59 #include "MMuonCalibPar.h" … … 69 67 // Setup histograms 70 68 // 71 MHMuonPar::MHMuonPar(const char *name, const char *title) :fGeomCam(NULL), fMuonSearchPar(NULL),72 fMuon CalibPar(NULL)69 MHMuonPar::MHMuonPar(const char *name, const char *title) : 70 fMuonSearchPar(NULL), fMuonCalibPar(NULL) 73 71 { 74 72 fName = name ? name : "MHMuonPar"; … … 76 74 77 75 fHistRadius.SetName("Radius"); 78 fHistRadius.SetTitle(" Radius");79 fHistRadius.SetXTitle("R adius[deg]");76 fHistRadius.SetTitle("Distribution of Radius'"); 77 fHistRadius.SetXTitle("R [\\circ]"); 80 78 fHistRadius.SetYTitle("Counts"); 81 79 fHistRadius.SetDirectory(NULL); … … 84 82 85 83 fHistArcWidth.SetName("ArcWidth"); 86 fHistArcWidth.SetTitle(" ArcWidth");87 fHistArcWidth.SetXTitle(" ArcWidth[deg]");84 fHistArcWidth.SetTitle("Distribution of ArcWidth"); 85 fHistArcWidth.SetXTitle("W [\\circ]"); 88 86 fHistArcWidth.SetYTitle("Counts"); 89 87 fHistArcWidth.SetDirectory(NULL); … … 91 89 fHistArcWidth.SetFillStyle(4000); 92 90 93 fHistBroad.SetName("Ring broadening");94 fHistBroad.SetTitle(" Ringbroadening");95 fHistBroad.SetXTitle("R adius[deg]");96 fHistBroad.SetYTitle(" ArcWidth/Radius");91 fHistBroad.SetName("RingBroadening"); 92 fHistBroad.SetTitle("Profile of ArcWidth/Radius versus Radius"); 93 fHistBroad.SetXTitle("R [\\circ]"); 94 fHistBroad.SetYTitle("W/R [1]"); 97 95 fHistBroad.SetDirectory(NULL); 98 96 fHistBroad.UseCurrentStyle(); … … 100 98 101 99 fHistSize.SetName("SizeVsRadius"); 102 fHistSize.SetXTitle("Radius"); 103 fHistSize.SetYTitle("MuonSize"); 100 fHistSize.SetTitle("Profile of Muon Size vs. Radius"); 101 fHistSize.SetXTitle("Rs [[\circ]"); 102 fHistSize.SetYTitle("S [phe]"); 104 103 fHistSize.SetDirectory(NULL); 105 104 fHistSize.UseCurrentStyle(); … … 119 118 bins.SetEdges(20, 0.5, 1.5); 120 119 bins.Apply(fHistSize); 121 122 120 } 123 121 … … 129 127 Bool_t MHMuonPar::SetupFill(const MParList *plist) 130 128 { 131 fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");132 if (! fGeomCam)129 MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam"); 130 if (!geom) 133 131 { 134 132 *fLog << warn << "MGeomCam not found... abort." << endl; 135 133 return kFALSE; 136 134 } 135 fMm2Deg = geom->GetConvMm2Deg(); 136 137 137 fMuonSearchPar = (MMuonSearchPar*)plist->FindObject("MMuonSearchPar"); 138 138 if (!fMuonSearchPar) … … 148 148 } 149 149 150 fMm2Deg = fGeomCam->GetConvMm2Deg(); 151 152 ApplyBinning(*plist, "Radius", &fHistRadius); 153 154 ApplyBinning(*plist, "ArcWidth", &fHistArcWidth); 155 156 ApplyBinning(*plist, "Ringbroadening", &fHistBroad); 157 158 ApplyBinning(*plist, "SizeVsRadius", &fHistSize); 150 ApplyBinning(*plist, "Radius", &fHistRadius); 151 ApplyBinning(*plist, "ArcWidth", &fHistArcWidth); 152 ApplyBinning(*plist, "RingBroadening", &fHistBroad); 153 ApplyBinning(*plist, "SizeVsRadius", &fHistSize); 159 154 160 155 return kTRUE; … … 211 206 fHistBroad.Draw(); 212 207 } 213 208 /* 214 209 TH1 *MHMuonPar::GetHistByName(const TString name) 215 210 { … … 220 215 return NULL; 221 216 } 222 217 */ -
trunk/MagicSoft/Mars/mmuon/MHMuonPar.h
r6973 r6979 12 12 #endif 13 13 14 class MHillas;15 14 class MMuonSearchPar; 16 15 class MMuonCalibPar; … … 20 19 { 21 20 private: 22 TH1F fHistRadius; // 21 TH1F fHistRadius; // 22 TH1F fHistArcWidth; // 23 23 24 TH1F fHistArcWidth; // 25 26 TProfile fHistBroad; // Area of used pixels 27 24 TProfile fHistBroad; // Area of used pixels 28 25 TProfile fHistSize; // [ratio] concentration ratio: sum of the two highest pixels / fSize 29 26 30 MGeomCam *fGeomCam; //! Camera geometry for plots (for the moment this is a feature for a loop only!) 31 32 MMuonSearchPar *fMuonSearchPar;//! 33 34 MMuonCalibPar *fMuonCalibPar;//! 27 MMuonSearchPar *fMuonSearchPar; //! 28 MMuonCalibPar *fMuonCalibPar; //! 35 29 36 30 Float_t fMm2Deg; … … 42 36 Bool_t Fill(const MParContainer *par, const Stat_t w=1); 43 37 44 TH1 *GetHistByName(const TString name);38 //TH1 *GetHistByName(const TString name); 45 39 46 TH1F &GetHistRadius() { return fHistRadius; } 47 48 TH1F &GetHistArcWidth() { return fHistArcWidth; } 49 50 TProfile &GetHistBroad() { return fHistBroad; } 51 52 TProfile &GetHistSize() { return fHistSize; } 40 const TH1F& GetHistRadius() const { return fHistRadius; } 41 const TH1F& GetHistArcWidth() const { return fHistArcWidth; } 42 const TProfile& GetHistBroad() const { return fHistBroad; } 43 const TProfile& GetHistSize() const { return fHistSize; } 53 44 54 45 void Draw(Option_t *opt=""); -
trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc
r6973 r6979 17 17 ! 18 18 ! Author(s): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de> 19 ! 19 ! Author(s): Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de> 20 20 ! 21 ! Copyright: MAGIC Software Development, 2000-200 421 ! Copyright: MAGIC Software Development, 2000-2005 22 22 ! 23 23 ! … … 30 30 // Storage Container for muon 31 31 // 32 // This class holds some information for a calibra ion using muons. Muons32 // This class holds some information for a calibration using muons. Muons 33 33 // are identified by using the class of the MMuonSearchParCalc. You can fill 34 34 // these information by using the MMuonCalibParCalc. See also these class 35 35 // manuals. 36 36 // 37 //38 // Input Containers:39 // [MGeomCam]40 // [MSignalCam]41 // [MMuonSearchPar]42 //43 37 ///////////////////////////////////////////////////////////////////////////// 44 38 #include "MMuonCalibPar.h" 45 39 46 #include <fstream>47 48 40 #include "MLog.h" 49 41 #include "MLogManip.h" 50 #include "MSignalCam.h"51 #include "MSignalPix.h"52 #include "MMuonSearchPar.h"53 #include "MBinning.h"54 42 55 43 using namespace std; … … 66 54 fTitle = title ? title : "Muon calibration parameters"; 67 55 68 fHistPhi = new TH1F; 69 fHistWidth = new TH1F; 70 71 fHistPhi->SetName("HistPhi"); 72 fHistPhi->SetTitle("HistPhi"); 73 fHistPhi->SetXTitle("phi [deg.]"); 74 fHistPhi->SetYTitle("sum of ADC"); 75 fHistPhi->SetDirectory(NULL); 76 fHistPhi->SetFillStyle(4000); 77 fHistPhi->UseCurrentStyle(); 78 79 fHistWidth->SetName("HistWidth"); 80 fHistWidth->SetTitle("HistWidth"); 81 fHistWidth->SetXTitle("distance from the ring center [deg.]"); 82 fHistWidth->SetYTitle("sum of ADC"); 83 fHistWidth->SetDirectory(NULL); 84 fHistWidth->SetFillStyle(4000); 85 fHistWidth->UseCurrentStyle(); 86 87 } 88 89 // -------------------------------------------------------------------------- 90 // 91 MMuonCalibPar::~MMuonCalibPar() 92 { 93 delete fHistPhi; 94 delete fHistWidth; 56 Reset(); 95 57 } 96 58 … … 99 61 void MMuonCalibPar::Reset() 100 62 { 101 fArcLength = -1.; 102 fArcPhi = 0.; 103 fArcWidth = -1.; 104 fChiArcPhi = -1.; 105 fChiArcWidth = -1.; 106 fMuonSize = 0.; 107 fEstImpact = -1.; 108 fUseUnmap = kFALSE; 109 fPeakPhi = 0.; 110 111 fHistPhi->Reset(); 112 fHistWidth->Reset(); 63 fArcLength = -1.; 64 fArcPhi = 0.; 65 fArcWidth = -1.; 66 fChiArcPhi = -1.; 67 fChiArcWidth = -1.; 68 fMuonSize = 0.; 69 fEstImpact = -1.; 70 fPeakPhi = 0.; 113 71 } 114 72 … … 116 74 { 117 75 *fLog << all; 118 *fLog << "Muon Parameters (" << GetName() << ")" << endl; 119 *fLog << " - Arc Length [deg.] = " << fArcLength << endl; 120 *fLog << " - Arc Phi [deg.] = " << fArcPhi << endl; 121 *fLog << " - Arc Width [deg.] = " << fArcWidth << endl; 122 *fLog << " - Chi Arc Phi [x2/ndf]= " << fChiArcPhi << endl; 123 *fLog << " - Chi Arc Width [x2/ndf]= " << fChiArcWidth << endl; 124 *fLog << " - Est. I. P. [m] = " << fEstImpact << endl; 125 *fLog << " - Size of muon = " << fMuonSize << endl; 126 *fLog << " - Peak Phi [deg.] = " << fPeakPhi << endl; 127 *fLog << " - UseUnmap = " << fUseUnmap << endl; 76 *fLog << "Muon Parameters (" << GetName() << ")" << endl; 77 *fLog << " - Arc Length [deg] = " << fArcLength << endl; 78 *fLog << " - Arc Phi [deg] = " << fArcPhi << endl; 79 *fLog << " - Arc Width [deg] = " << fArcWidth << endl; 80 *fLog << " - Chi Arc Phi [x2/ndf]= " << fChiArcPhi << endl; 81 *fLog << " - Chi Arc Width [x2/ndf]= " << fChiArcWidth << endl; 82 *fLog << " - Est. I. P. [m] = " << fEstImpact << endl; 83 *fLog << " - Size of muon [phe] = " << fMuonSize << endl; 84 *fLog << " - Peak Phi [deg] = " << fPeakPhi << endl; 128 85 } 129 86 -
trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.h
r6973 r6979 16 16 { 17 17 private: 18 Float_t fArcLength; // An arc length of a muon along the arc [deg.]19 Float_t fArcPhi; // A opening angle of a muon arc [deg.]20 Float_t fArcWidth; // A width of a muon [deg.] (1 sigma of gaussian fit)21 Float_t fChiArcPhi; // A chisquare value of the cosine fit for arc phi22 Float_t fChiArcWidth; // A chisquare value of the cosine fit for arc wid23 Float_t fMuonSize; // A SIZE of muon which is defined as a SIZE around the estimated circle24 Float_t fEstImpact; // An estimated impact parameter from the photon distribution along the arc image25 Bool_t fUseUnmap; // This is a flag to know the Unmapped pixels are used. Refer to the class of MImgCleanStd26 Float_t fPeakPhi; // The angle which indicates the peak position in the estimated circle18 Float_t fArcLength; // An arc length of a muon along the arc [deg.] 19 Float_t fArcPhi; // A opening angle of a muon arc [deg.] 20 Float_t fArcWidth; // A width of a muon [deg.] (1 sigma of gaussian fit) 21 Float_t fChiArcPhi; // A chisquare value of the cosine fit for arc phi 22 Float_t fChiArcWidth; // A chisquare value of the cosine fit for arc wid 23 Float_t fMuonSize; // A SIZE of muon which is defined as a SIZE around the estimated circle 24 Float_t fEstImpact; // An estimated impact parameter from the photon distribution along the arc image 25 //Bool_t fUseUnmap; // This is a flag to know the Unmapped pixels are used. Refer to the class of MImgCleanStd 26 Float_t fPeakPhi; // The angle which indicates the peak position in the estimated circle 27 27 28 28 public: 29 MMuonCalibPar(const char *name=NULL, const char *title=NULL); 30 ~MMuonCalibPar(); 31 32 TH1F *fHistPhi; // Histogram of photon distribution along the arc. 33 TH1F *fHistWidth; // Histogram of radial photon distribution of the arc. 34 35 void Reset(); 36 37 Float_t GetArcLength() const { return fArcLength; } 38 Float_t GetArcPhi() const { return fArcPhi; } 39 Float_t GetArcWidth() const { return fArcWidth; } 40 Float_t GetChiArcPhi() const { return fChiArcPhi; } 41 Float_t GetChiArcWidth() const { return fChiArcWidth; } 42 Float_t GetMuonSize() const { return fMuonSize; } 43 Float_t GetEstImpact() const { return fEstImpact; } 44 Bool_t IsUseUnmap() const { return fUseUnmap; } 45 Float_t GetPeakPhi() const { return fPeakPhi; } 46 TH1F *GetHistPhi() { return fHistPhi; } 47 TH1F *GetHistWidth() { return fHistWidth; } 48 49 void SetArcLength(Float_t len) { fArcLength = len; } 50 void SetArcPhi(Float_t phi) { fArcPhi = phi; } 51 void SetArcWidth(Float_t wid) { fArcWidth = wid; } 52 void SetChiArcPhi(Float_t chi) { fChiArcPhi = chi; } 53 void SetChiArcWidth(Float_t chi) { fChiArcWidth = chi; } 54 void SetMuonSize(Float_t size) { fMuonSize = size; } 55 void SetEstImpact(Float_t impact) { fEstImpact = impact; } 56 void SetUseUnmap() { fUseUnmap = kTRUE; } 57 void SetPeakPhi(Float_t phi) { fPeakPhi = phi; } 58 59 void Print(Option_t *opt=NULL) const; 29 MMuonCalibPar(const char *name=NULL, const char *title=NULL); 30 //~MMuonCalibPar(); 60 31 61 ClassDef(MMuonCalibPar, 1) // Container to hold muon calibration parameters 32 void Reset(); 33 34 Float_t GetArcLength() const { return fArcLength; } 35 Float_t GetArcPhi() const { return fArcPhi; } 36 Float_t GetArcWidth() const { return fArcWidth; } 37 Float_t GetChiArcPhi() const { return fChiArcPhi; } 38 Float_t GetChiArcWidth() const { return fChiArcWidth; } 39 Float_t GetMuonSize() const { return fMuonSize; } 40 Float_t GetEstImpact() const { return fEstImpact; } 41 //Bool_t IsUseUnmap() const { return fUseUnmap; } 42 Float_t GetPeakPhi() const { return fPeakPhi; } 43 // TH1F *GetHistPhi() { return fHistPhi; } 44 // TH1F *GetHistWidth() { return fHistWidth; } 45 46 void SetArcLength(Float_t len) { fArcLength = len; } 47 void SetArcPhi(Float_t phi) { fArcPhi = phi; } 48 void SetArcWidth(Float_t wid) { fArcWidth = wid; } 49 void SetChiArcPhi(Float_t chi) { fChiArcPhi = chi; } 50 void SetChiArcWidth(Float_t chi) { fChiArcWidth = chi; } 51 void SetMuonSize(Float_t size) { fMuonSize = size; } 52 void SetEstImpact(Float_t impact) { fEstImpact = impact; } 53 //void SetUseUnmap() { fUseUnmap = kTRUE; } 54 void SetPeakPhi(Float_t phi) { fPeakPhi = phi; } 55 56 void Print(Option_t *opt=NULL) const; 57 58 ClassDef(MMuonCalibPar, 1) // Container to hold muon calibration parameters 62 59 }; 63 60 -
trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc
r6973 r6979 67 67 // 68 68 // 69 // For the faster computation, by default, the calculation of impact70 // parameter is suppressed. If you want to calculate the impact parameter71 // from the muon image, you can use the function of EnableImpactCalc(),72 // namely;73 // 74 // mucalibcalc.EnableImpactCalc();.69 // // For the faster computation, by default, the calculation of impact 70 // // parameter is suppressed. If you want to calculate the impact parameter 71 // // from the muon image, you can use the function of EnableImpactCalc(), 72 // // namely; 73 // // 74 // // mucalibcalc.EnableImpactCalc();. 75 75 // 76 76 // In addition, for the faster comutation, pre cuts to select the candidates … … 98 98 // 99 99 // Input Containers: 100 // [MGeomCam]101 // [MSignalCam]102 // [MMuonSearchPar]100 // MGeomCam 101 // MSignalCam 102 // MMuonSearchPar 103 103 // 104 104 // Output Containers: 105 // [MMuonCalibPar]105 // MMuonCalibPar 106 106 // 107 107 ////////////////////////////////////////////////////////////////////////////// 108 108 109 109 #include "MMuonCalibParCalc.h" 110 111 #include <fstream>112 110 113 111 #include <TH1.h> … … 115 113 #include <TMinuit.h> 116 114 115 #include "MLog.h" 116 #include "MLogManip.h" 117 117 118 #include "MParList.h" 118 119 119 120 #include "MGeomCam.h" 120 121 #include "MGeomPix.h" 121 #include "MSrcPosCam.h" 122 122 123 #include "MSignalCam.h" 123 #include "MMuonSearchPar.h" 124 124 125 #include "MMuonCalibPar.h" 125 126 #include "MMuonSetup.h" 126 #include "MLog.h" 127 #include "MLogManip.h" 128 #include "MBinning.h" 127 #include "MMuonSearchPar.h" 128 #include "MHSingleMuon.h" 129 129 130 130 using namespace std; … … 140 140 // 141 141 MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title) 142 // : fEnableImpactCalc(kFALSE) 142 143 { 143 144 fName = name ? name : gsDefName.Data(); 144 145 fTitle = title ? title : gsDefTitle.Data(); 145 146 fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.147 146 } 148 147 … … 151 150 Int_t MMuonCalibParCalc::PreProcess(MParList *pList) 152 151 { 153 fSignalCam = (MSignalCam*)pList->FindObject("MSignalCam");154 if (!fSignalCam)152 fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam"); 153 if (!fGeomCam) 155 154 { 156 *fLog << dbginf << "MSignalCam not found... aborting." << endl;157 return kFALSE;155 *fLog << err << "MGeomCam not found... abort." << endl; 156 return kFALSE; 158 157 } 159 160 fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");161 if (!fGeomCam)158 159 fHist = (MHSingleMuon*)pList->FindObject("MHSingleMuon"); 160 if (!fHist) 162 161 { 163 *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;164 return kFALSE;162 *fLog << err << "MHSingleMuon not found... abort." << endl; 163 return kFALSE; 165 164 } 166 167 fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar", "MMuonCalibPar"); 168 if (!fMuonCalibPar) 169 return kFALSE; 170 171 fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", "MMuonSearchPar"); 172 if (!fMuonSearchPar) 173 return kFALSE; 174 175 fMuonSetup = (MMuonSetup*)pList->FindCreateObj("MMuonSetup", "MMuonSetup"); 176 if (!fMuonSetup) 177 return kFALSE; 178 179 return kTRUE; 165 166 fMuonSetup = (MMuonSetup*)pList->FindObject("MMuonSetup"); 167 if (!fMuonSetup) 168 { 169 *fLog << err << "MMuonSetup not found... abort." << endl; 170 return kFALSE; 171 } 172 173 fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar"); 174 if (!fMuonCalibPar) 175 return kFALSE; 176 177 fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar"); 178 if (!fMuonSearchPar) 179 return kFALSE; 180 181 return kTRUE; 180 182 } 181 182 // --------------------------------------------------------------------------183 //184 // This function fill the histograms in order to get muon parameters.185 // For the evaluation of the Arc Width, we use only the signals in the inner186 // part.187 //188 void MMuonCalibParCalc::FillHist()189 {190 Float_t MuonSize = 0.;191 192 Int_t binnumphi = fMuonSetup->fArcPhiBinNum;193 Int_t binnumwid = fMuonSetup->fArcWidthBinNum;194 195 // preparation for a histgram196 MBinning binsphi;197 binsphi.SetEdges(binnumphi,198 fMuonSetup->fArcPhiHistStartVal,199 fMuonSetup->fArcPhiHistEndVal);200 binsphi.Apply(*(fMuonCalibPar->fHistPhi));201 202 MBinning binswid;203 binswid.SetEdges(binnumwid,204 fMuonSetup->fArcWidthHistStartVal,205 fMuonSetup->fArcWidthHistEndVal);206 binswid.Apply(*(fMuonCalibPar->fHistWidth));207 208 const Int_t entries = (*fSignalCam).GetNumPixels();209 210 // the position of the center of a muon ring211 const Float_t cenx = (*fMuonSearchPar).GetCenterX();212 const Float_t ceny = (*fMuonSearchPar).GetCenterY();213 214 for (Int_t i=0; i<entries; i++ )215 {216 MSignalPix &pix = (*fSignalCam)[i];217 218 const MGeomPix &gpix = (*fGeomCam)[i/*pix.GetPixId()*/];219 220 const Float_t dx = gpix.GetX() - cenx;221 const Float_t dy = gpix.GetY() - ceny;222 223 const Float_t dist = TMath::Sqrt(dx*dx+dy*dy);224 225 Float_t ang = TMath::ACos(dx/dist);226 if(dy>0)227 ang *= -1.0;228 229 // if the signal is not near the estimated circle, it is ignored.230 if(dist < (*fMuonSearchPar).GetRadius() + fMuonSetup->GetMargin()231 && dist > (*fMuonSearchPar).GetRadius() - fMuonSetup->GetMargin())232 {233 // check whether ummapped pixel is used or not.234 // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information.235 if(pix.IsPixelUnmapped())236 {237 fMuonCalibPar->SetUseUnmap();238 continue;239 }240 fMuonCalibPar->fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons());241 MuonSize += pix.GetNumPhotons();242 }243 244 // use only the inner pixles. This is geometry dependent. This has to245 // be fixed!246 if(i>397)247 continue;248 249 fMuonCalibPar->fHistWidth250 ->Fill(dist*(*fGeomCam).GetConvMm2Deg(), pix.GetNumPhotons());251 }252 253 fMuonCalibPar->SetMuonSize(MuonSize);254 255 // error estimation (temporaly)256 // The error is estimated from the signal. In order to do so, we have to257 // once convert the signal from ADC to photo-electron. Then we can get258 // the fluctuation such as F-factor*sqrt(phe).259 // Up to now, the error of pedestal is not taken into accout. This is not260 // of course correct. We will include this soon.261 Double_t ADC2PhEl = 0.14;262 Double_t Ffactor = 1.4;263 for(Int_t i=0; i<binnumphi+1; i++)264 {265 Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;266 {267 fMuonCalibPar->fHistPhi->SetBinError(i, rougherr);268 }269 }270 for(Int_t i=0; i<binnumwid+1; i++)271 {272 Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;273 {274 fMuonCalibPar->fHistWidth->SetBinError(i, rougherr);275 }276 }277 }278 279 // --------------------------------------------------------------------------280 //281 // Photon distribution along the estimated circle is fitted with theoritical282 // function in order to get some more information such as Arc Phi and Arc283 // Length.284 //285 void MMuonCalibParCalc::CalcPhi()286 {287 Float_t thres = fMuonSetup->GetArcPhiThres();288 Float_t startval = fMuonSetup->fArcPhiHistStartVal;289 Float_t endval = fMuonSetup->fArcPhiHistEndVal;290 Int_t binnum = fMuonSetup->fArcPhiBinNum;291 292 Float_t convbin2val = (endval-startval)/(Float_t)binnum;293 294 // adjust the peak to 0295 Float_t maxval = 0.;296 Int_t maxbin = 0;297 maxval = fMuonCalibPar->fHistPhi->GetMaximum();298 maxbin = fMuonCalibPar->fHistPhi->GetMaximumBin();299 fMuonCalibPar->SetPeakPhi(180.-(Float_t)(maxbin-1.)*convbin2val);300 TArrayD tmp;301 tmp.Set(binnum+1);302 for(Int_t i=1; i<binnum+1; i++)303 {304 tmp[i] = fMuonCalibPar->fHistPhi->GetBinContent(i);305 }306 for(Int_t i=1; i<binnum+1; i++)307 {308 Int_t id;309 id = i + (maxbin-(Int_t)((Float_t)binnum/2.)-1);310 if(id>binnum)311 {312 id-=(binnum);313 }314 if(id<=0)315 {316 id+=(binnum);317 }318 fMuonCalibPar->fHistPhi->SetBinContent(i,tmp[id]);319 }320 maxbin = (Int_t)((Float_t)binnum/2.)+1;321 322 // Determination of fitting region323 // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore,324 // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!!325 Float_t startfitval = 0.;326 Float_t endfitval = 0.;327 Bool_t IsInMaxim = kFALSE;328 Int_t effbinnum = 0;329 for(Int_t i=1; i<binnum+1; i++)330 {331 Float_t content = fMuonCalibPar->fHistPhi->GetBinContent(i);332 Float_t content_pre = fMuonCalibPar->fHistPhi->GetBinContent(i-1);333 334 if(content > thres && content_pre < thres)335 {336 startfitval = (Float_t)(i-1)*convbin2val+startval;337 }338 if(i==maxbin)339 IsInMaxim = kTRUE;340 341 if(content < thres && IsInMaxim == kTRUE)342 {343 endfitval = (Float_t)(i-1)*convbin2val+startval;344 break;345 }346 endfitval = endval;347 }348 349 effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);350 351 fMuonCalibPar->SetArcPhi(endfitval-startfitval);352 353 fMuonCalibPar->SetArcLength( fMuonCalibPar->GetArcPhi()*TMath::DegToRad()*(*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg());354 355 if(fEnableImpactCalc)356 CalcImpact(effbinnum, startfitval, endfitval);357 }358 359 // --------------------------------------------------------------------------360 //361 // An impact parameter is calculated by fitting the histogram of photon362 // distribution along the circle with a theoritical model.363 // (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.364 // The function (6) is used here.)365 //366 // By default this calculation is suppressed because this calculation is367 // very time consuming. If you want to calculate an impact parameter,368 // you can call the function of EnableImpactCalc().369 //370 void MMuonCalibParCalc::CalcImpact371 (Int_t effbinnum, Float_t startfitval, Float_t endfitval)372 {373 // Fit the distribution with Vacanti function. The function is different374 // for the impact parameter of inside or outside of our reflector.375 // Then two different functions are applied to the photon distribution,376 // and the one which give us smaller chisquare value is taken as a377 // proper one.378 Double_t val1,err1,val2,err2;379 // impact parameter inside mirror radius (8.5m)380 TString func1;381 Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();382 tmpval = sin(2.*tmpval)*8.5;383 func1 += "[0]*";384 func1 += tmpval;385 func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";386 387 TF1 f1("f1",func1,startfitval,endfitval);388 f1.SetParameters(2000,3,0);389 f1.SetParLimits(1,0,8.5);390 f1.SetParLimits(2,-180.,180.);391 392 fMuonCalibPar->fHistPhi->Fit("f1","RQEM");393 394 Float_t chi1 = -1;395 Float_t chi2 = -1;396 if(effbinnum>3)397 chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));398 399 gMinuit->GetParameter(1,val1,err1); // get the estimated IP400 Float_t estip1 = val1;401 402 // impact parameter beyond mirror area (8.5m)403 TString func2;404 Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();405 tmpval2 = sin(2.*tmpval2)*8.5*2.;406 func2 += "[0]*";407 func2 += tmpval2;408 func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";409 410 TF1 f2("f2",func2,startfitval,endfitval);411 f2.SetParameters(2000,20,0);412 f2.SetParLimits(1,8.5,300.);413 f2.SetParLimits(2,-180.,180.);414 415 fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");416 417 if(effbinnum>3)418 chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));419 420 gMinuit->GetParameter(1,val2,err2); // get the estimated IP421 Float_t estip2 = val2;422 423 if(effbinnum<=3)424 {425 fMuonCalibPar->SetEstImpact(-1.);426 }427 if(chi2 > chi1)428 {429 fMuonCalibPar->SetEstImpact(estip1);430 fMuonCalibPar->SetChiArcPhi(chi1);431 }432 else433 {434 fMuonCalibPar->SetEstImpact(estip2);435 fMuonCalibPar->SetChiArcPhi(chi2);436 }437 }438 439 // --------------------------------------------------------------------------440 //441 // Photon distribution of distance from the center of estimated ring is442 // fitted in order to get some more information such as ARC WIDTH which443 // can represent to the PSF of our reflector.444 //445 Float_t MMuonCalibParCalc::CalcWidth()446 {447 Float_t startval = fMuonSetup->fArcWidthHistStartVal;448 Float_t endval = fMuonSetup->fArcWidthHistEndVal;449 Int_t binnum = fMuonSetup->fArcWidthBinNum;450 Float_t thres = fMuonSetup->GetArcWidthThres();451 452 Float_t convbin2val = (endval - startval)453 /(Float_t)binnum;454 455 // determination of fitting region456 Int_t maxbin = fMuonCalibPar->fHistWidth->GetMaximumBin();457 Float_t startfitval = 0.;458 Float_t endfitval = 0.;459 Bool_t IsInMaxim = kFALSE;460 Int_t effbinnum = 0;461 for(Int_t i=1; i<binnum+1; i++)462 {463 Float_t content = fMuonCalibPar->fHistWidth->GetBinContent(i);464 Float_t content_pre = fMuonCalibPar->fHistWidth->GetBinContent(i-1);465 466 if(content > thres)467 effbinnum++;468 469 if(content > thres && content_pre < thres)470 {471 startfitval = (Float_t)(i-4)*convbin2val + startval;472 if(startfitval<0.) startfitval = 0.;473 }474 if(i==maxbin)475 IsInMaxim = kTRUE;476 477 if(content < thres && IsInMaxim == kTRUE)478 {479 endfitval = (Float_t)(i+2)*convbin2val + startval;480 if(endfitval>180.) endfitval = 180.;481 break;482 }483 endfitval = endval;484 }485 effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);486 487 TF1 f1("f1","gaus",startfitval,endfitval);488 489 // options : N do not store the function, do not draw490 // I use integral of function in bin rather than value at bin center491 // R use the range specified in the function range492 // Q quiet mode493 fMuonCalibPar->fHistWidth->Fit("f1","QR","",startfitval,endfitval);494 495 if(effbinnum>3)496 fMuonCalibPar->SetChiArcWidth(f1.GetChisquare()/((Float_t)(effbinnum-3)));497 498 Double_t val,err;499 gMinuit->GetParameter(2,val,err); // get the sigma value500 501 return val;502 }503 504 // --------------------------------------------------------------------------505 //506 // Calculation of muon parameters507 //508 Int_t MMuonCalibParCalc::Calc()509 {510 // initialization511 (*fMuonCalibPar).Reset();512 513 // Fill signals to histograms514 FillHist();515 516 // Calculation of Arc Phi etc...517 CalcPhi();518 519 // Calculation of Arc Width etc...520 fMuonCalibPar->SetArcWidth(CalcWidth());521 522 if(fMuonCalibPar->GetArcPhi()>160 && fMuonSearchPar->GetRadius()>170 &&523 fMuonSearchPar->GetRadius()<400 && fMuonSearchPar->GetDeviation()<50)524 fMuonCalibPar->SetReadyToSave();525 526 return kTRUE;527 }528 183 529 184 // ------------------------------------------------------------------------- … … 531 186 Int_t MMuonCalibParCalc::Process() 532 187 { 533 534 Calc(); 535 536 return kTRUE; 188 // Calculation of Arc Phi etc... 189 const Float_t thresphi = fMuonSetup->GetThresholdArcPhi() /fGeomCam->GetConvMm2Deg(); 190 const Float_t threswidth = fMuonSetup->GetThresholdArcWidth()/fGeomCam->GetConvMm2Deg(); 191 192 Double_t peakphi, arcphi; 193 Double_t width, chi; 194 195 if (!fHist->CalcPhi(thresphi, peakphi, arcphi)) 196 return kTRUE; 197 198 if (!fHist->CalcWidth(threswidth, width, chi)) 199 return kTRUE; 200 201 // Get Muon Size 202 fMuonCalibPar->SetMuonSize(fHist->GetHistPhi().Integral()); 203 204 // Calculate Arc Length 205 fMuonCalibPar->SetPeakPhi(peakphi); 206 fMuonCalibPar->SetArcPhi(arcphi); 207 208 const Float_t conv = TMath::RadToDeg()*fGeomCam->GetConvMm2Deg(); 209 fMuonCalibPar->SetArcLength(fMuonCalibPar->GetArcPhi()*fMuonSearchPar->GetRadius()*conv); 210 211 // Calculation of Arc Width etc... 212 fMuonCalibPar->SetChiArcWidth(chi); 213 fMuonCalibPar->SetArcWidth(width); 214 215 // Check if this is a 'Write-Out' candidate 216 if (fMuonCalibPar->GetArcPhi()>160 && fMuonSearchPar->GetRadius()<400 && 217 fMuonSearchPar->GetDeviation()<50 && fMuonSearchPar->GetRadius()>170) 218 fMuonCalibPar->SetReadyToSave(); 219 220 return kTRUE; 537 221 } -
trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h
r6973 r6979 8 8 class MMuonSearchPar; 9 9 class MMuonCalibPar; 10 class MSrcPosCam;11 10 class MGeomCam; 12 class MSignalCam;13 11 class MMuonSetup; 12 class MHSingleMuon; 14 13 15 14 class MMuonCalibParCalc : public MTask 16 15 { 17 16 private: 18 MGeomCam *fGeomCam; 19 M SignalCam *fSignalCam;20 MMuon CalibPar *fMuonCalibPar;21 MMuonSe archPar *fMuonSearchPar;22 M MuonSetup *fMuonSetup;17 MGeomCam *fGeomCam; //! 18 MMuonCalibPar *fMuonCalibPar; //! 19 MMuonSearchPar *fMuonSearchPar; //! 20 MMuonSetup *fMuonSetup; //! 21 MHSingleMuon *fHist; //! 23 22 24 Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time.23 //Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time. 25 24 26 Int_t PreProcess(MParList *plist); 27 Int_t Process(); 25 Int_t PreProcess(MParList *plist); 26 Int_t Process(); 27 28 void FillHist(); 29 void CalcPhi(); 30 void CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval); 31 Float_t CalcWidth(); 28 32 29 33 public: 30 34 MMuonCalibParCalc(const char *name=NULL, const char *title=NULL); 31 35 32 void EnableImpactCalc() { fEnableImpactCalc = kTRUE; } 33 34 void FillHist(); 35 void CalcPhi(); 36 void CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval); 37 Float_t CalcWidth(); 38 Int_t Calc(); 36 //void EnableImpactCalc(Bool_t b=kTRUE) { fEnableImpactCalc = b; } 39 37 40 38 ClassDef(MMuonCalibParCalc, 0) // task to calculate muon parameters -
trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc
r6973 r6979 17 17 ! 18 18 ! Author(s): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de> 19 ! 20 ! 21 ! Copyright: MAGIC Software Development, 2000-200 419 ! Author(s): Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de> 20 ! 21 ! Copyright: MAGIC Software Development, 2000-2005 22 22 ! 23 23 ! … … 40 40 // MMuonCalibPar. The information will be available by using the task of 41 41 // MMuonCalibParCalc. 42 //43 42 // 44 43 // --- How to search muons --- … … 47 46 // 48 47 // 1. A temporal center position of a circle is determined by using 49 // the Hillas parameters. Assumed that the center position will be on the 50 // line which is perpendicular to the longitudinal image axis and the 51 // distance from the gravity center of the image to the center position of 52 // a ring is approximately 1 deg. (corresponding to the Cherenkov angle.). 53 // Therefore, we will have two candidates of the center positions. 48 // the Hillas parameters. Assumed that the center position will be on the 49 // line which is perpendicular to the longitudinal image axis and the 50 // distance from the gravity center of the image to the center position of 51 // a ring is approximately 1 deg. (corresponding to the Cherenkov angle.). 52 // Therefore, we will have two candidates of the center positions. 53 // 54 54 // 2. Find the ring radius which gives the minimum RMS between the camera 55 // images and the estimated circle. 55 // images and the estimated circle. 56 // 56 57 // 3. Select one temporal position which gives smaller RMS as a true temporal 57 // center position. 58 // center position. 59 // 58 60 // 4. Changing the center position of a circle on the camera plane from the 59 // determined temporal center position, find the position which gives the 60 // minimum RMS of the fit. 61 // 62 // 63 // --- Remark --- 64 // This method to search for muons is not fully optimized yet. However, 65 // it is good idea to use the temporal estimated center position from 66 // the Hillas parameters in order to reduce time to search. In addition, 67 // This method is faster than the MINUIT. 61 // determined temporal center position, find the position which gives the 62 // minimum RMS of the fit. 68 63 // 69 64 // 70 65 // Input Containers: 71 // [MGeomCam]72 // [MHillas]73 // [MSignalCam]66 // MGeomCam 67 // MHillas 68 // MSignalCam 74 69 // 75 70 ///////////////////////////////////////////////////////////////////////////// 76 71 #include "MMuonSearchPar.h" 77 72 78 #include < fstream>73 #include <TMinuit.h> 79 74 80 75 #include "MLog.h" 81 76 #include "MLogManip.h" 77 82 78 #include "MHillas.h" 79 83 80 #include "MGeomCam.h" 84 81 #include "MGeomPix.h" 82 85 83 #include "MSignalPix.h" 86 84 #include "MSignalCam.h" … … 96 94 MMuonSearchPar::MMuonSearchPar(const char *name, const char *title) 97 95 { 98 fName = name ? name : "MMuonSearchPar";99 fTitle = title ? title : "Muon search parameters";96 fName = name ? name : "MMuonSearchPar"; 97 fTitle = title ? title : "Muon search parameters"; 100 98 } 101 99 … … 104 102 void MMuonSearchPar::Reset() 105 103 { 106 fRadius = -1.; 107 fDeviation = -1.; 108 fCenterX = 0.; 109 fCenterY = 0.; 110 } 111 112 // -------------------------------------------------------------------------- 113 // 114 // Get the tempolary center of a ring from the Hillas parameters. 115 // Two candidates of the position is returened. 116 // 117 void MMuonSearchPar::CalcTempCenter(const MHillas &hillas, 118 Float_t &xtmp1, Float_t &ytmp1, Float_t &xtmp2, Float_t &ytmp2) 119 { 120 Float_t a,dx,dy; 121 Float_t tmp_r = 300.; // assume that the temporal cherenkov angle is 1 deg. (300 mm). 122 123 a = TMath::Tan(hillas.GetDelta()); 124 125 dx = a/TMath::Sqrt(tmp_r+a*a)/3.; 126 dy = -tmp_r/TMath::Sqrt(1+a*a)/3.; 127 128 xtmp1 = hillas.GetMeanX() + dx; 129 ytmp1 = hillas.GetMeanY() + dy; 130 xtmp2 = hillas.GetMeanX() - dx; 131 ytmp2 = hillas.GetMeanY() - dy; 104 fRadius = -1; 105 fDeviation = -1; 106 fCenterX = 0; 107 fCenterY = 0; 108 } 109 110 // -------------------------------------------------------------------------- 111 // 112 // This is a wrapper function to have direct access to the data members 113 // in the function calculating the minimization value. 114 // 115 void MMuonSearchPar::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) 116 { 117 const MMuonSearchPar *optim = (MMuonSearchPar*)gMinuit->GetObjectFit(); 118 f = optim->Fcn(par); 132 119 } 133 120 … … 137 124 // and its RMS for the input position. 138 125 // 139 Bool_t MMuonSearchPar::CalcRadius(const MGeomCam &geom, const MSignalCam &evt, 140 Float_t x, Float_t y, Float_t &r, Float_t &sigma) 141 { 142 Float_t mean_r=0., dev_r=0., sums=0., tmp=0.; 143 144 const Int_t entries = evt.GetNumPixels(); 145 146 for (Int_t i=0; i<entries; i++ ){ 147 const MSignalPix &pix = evt[i]; 148 149 if (!pix.IsPixelUsed()) 150 continue; 151 152 const MGeomPix &gpix = geom[i/*pix.GetPixId()*/]; 153 154 tmp=TMath::Sqrt((gpix.GetX()-x)*(gpix.GetX()-x) 155 +(gpix.GetY()-y)*(gpix.GetY()-y)); 156 157 mean_r += pix.GetNumPhotons()*tmp; 158 dev_r += pix.GetNumPhotons()*tmp*tmp; 159 sums += pix.GetNumPhotons(); 160 } 161 162 if(sums<1.E-10) 163 return kFALSE; 164 165 r = mean_r/sums; 166 167 if(dev_r/sums-(r)*(r)<1.E-10) 168 return kFALSE; 169 170 sigma = TMath::Sqrt(dev_r/sums-(r)*(r)); 171 172 return kTRUE; 173 } 126 Double_t MMuonSearchPar::Fcn(Double_t *par) const 127 { 128 const Int_t entries = fSignal.GetSize(); 129 130 Double_t meanr=0; 131 Double_t devr =0; 132 Double_t sums =0; 133 134 // It seems that the loop is easy enough for a compiler optimization. 135 // Using pointer arithmetics doesn't improve the speed of the fit. 136 for (Int_t i=0; i<entries; i++ ) 137 { 138 Double_t tmp = TMath::Hypot(fX[i]-par[0], fY[i]-par[1]); 139 140 sums += fSignal[i]; 141 meanr += fSignal[i] * tmp; 142 devr += fSignal[i] * tmp*tmp; 143 } 144 145 par[2] = meanr/sums; 146 par[3] = devr/sums - par[2]*par[2]; 147 148 return par[3]; 149 } 174 150 175 151 // -------------------------------------------------------------------------- … … 178 154 // RMS of the fit, changing the center position of the circle. 179 155 // 180 void MMuonSearchPar::CalcMinimumDeviation 181 ( const MGeomCam &geom, const MSignalCam &evt, Float_t x, Float_t y, 182 Float_t xcog, Float_t ycog, Float_t sigma, Float_t &opt_rad, 183 Float_t &new_sigma, Float_t &newx, Float_t &newy ) 184 { 185 Float_t delta = 3.; // 3 mm (1/10 of an inner pixel size) Step to move. 186 Float_t rad_tmp,sig_tmp; 187 Float_t r2; 188 189 while(1) 190 { 191 r2=(xcog-x)*(xcog-x)+(ycog-y)*(ycog-y); 192 // Exit if the new estimated radius is above 2 deg. (600 mm). 193 if(r2 > 360000.) 194 { 195 new_sigma=sigma; 196 opt_rad=rad_tmp; 197 newx=x; 198 newy=y; 199 break; 200 } 201 if(CalcRadius(geom,evt,x,y+delta,rad_tmp,sig_tmp)) 202 { 203 if(sig_tmp<sigma) 204 { 205 sigma=sig_tmp; 206 opt_rad=rad_tmp; 207 y=y+delta; 208 continue; 209 } 210 } 211 if(CalcRadius(geom,evt,x-delta,y,rad_tmp,sig_tmp)) 212 { 213 if(sig_tmp<sigma) 214 { 215 sigma=sig_tmp; 216 opt_rad=rad_tmp; 217 x=x-delta; 218 continue; 219 } 220 } 221 if(CalcRadius(geom,evt,x+delta,y,rad_tmp,sig_tmp)) 222 { 223 if(sig_tmp<sigma) 224 { 225 sigma=sig_tmp; 226 opt_rad=rad_tmp; 227 x=x+delta; 228 continue; 229 } 230 } 231 if(CalcRadius(geom,evt,x,y-delta,rad_tmp,sig_tmp)) 232 { 233 if(sig_tmp<sigma) 234 { 235 sigma=sig_tmp; 236 opt_rad=rad_tmp; 237 y=y-delta; 238 continue; 239 } 240 } 241 new_sigma=sigma; 242 newx=x; 243 newy=y; 244 break; 156 void MMuonSearchPar::CalcMinimumDeviation(const MGeomCam &geom, 157 const MSignalCam &evt, 158 Double_t &x, Double_t &y, 159 Double_t &sigma, Double_t &radius) 160 { 161 // ------- Make a temporaray copy of the signal --------- 162 const Int_t n = geom.GetNumPixels(); 163 164 fSignal.Set(n); 165 fX.Set(n); 166 fY.Set(n); 167 168 Int_t p=0; 169 for (int i=0; i<n; i++) 170 { 171 const MSignalPix &pix = evt[i]; 172 if (pix.IsPixelUsed()) 173 { 174 fSignal[p] = pix.GetNumPhotons(); 175 176 fX[p] = geom[i].GetX(); 177 fY[p] = geom[i].GetY(); 178 p++; 179 } 245 180 } 181 fSignal.Set(p); 182 183 184 // ----------------- Setup and call minuit ------------------- 185 const Float_t delta = 30.; // 3 mm (1/10 of an inner pixel size) Step to move. 186 const Double_t r = geom.GetMaxRadius()*2; 187 188 // Save gMinuit 189 TMinuit *minsave = gMinuit; 190 191 // Initialize TMinuit with 4 parameters 192 TMinuit minuit; 193 minuit.SetPrintLevel(-1); // Switch off printing 194 minuit.Command("set nowarn"); // Switch off warning 195 // Define Parameters 196 minuit.DefineParameter(0, "x", x, delta, -r, r); 197 minuit.DefineParameter(1, "y", y, delta, -r, r); 198 minuit.DefineParameter(2, "r", 0, 1, 0, 0); 199 minuit.DefineParameter(3, "sigma", 0, 1, 0, 0); 200 // Fix return parameters 201 minuit.FixParameter(2); 202 minuit.FixParameter(3); 203 // Setup Minuit for 'this' and use fit function fcn 204 minuit.SetObjectFit(this); 205 minuit.SetFCN(fcn); 206 207 // Perform Simplex minimization 208 Int_t err; 209 Double_t tmp[2] = { 0, 0 }; 210 minuit.mnexcm("simplex", tmp, 2, err); 211 212 // Get resulting parameters 213 Double_t error; 214 minuit.GetParameter(0, x, error); 215 minuit.GetParameter(1, y, error); 216 minuit.GetParameter(2, radius, error); 217 minuit.GetParameter(3, sigma, error); 218 219 sigma = sigma>0 ? TMath::Sqrt(sigma) : 0; 220 221 gMinuit = minsave; 246 222 } 247 223 … … 250 226 // Calculation of muon parameters 251 227 // 252 void MMuonSearchPar::Calc 253 (const MGeomCam &geom, const MSignalCam &evt, const MHillas &hillas) 254 { 255 Reset(); 256 257 Float_t xtmp1,xtmp2,ytmp1,ytmp2; 258 Float_t rad,dev,rad2,dev2; 259 Float_t opt_rad,new_sigma,newx,newy; 260 261 // gets temporaly center 262 CalcTempCenter(hillas,xtmp1,ytmp1,xtmp2,ytmp2); 263 264 // determine which position will be the true position. Here mainly 265 // the curvature of a muon arc is relied on. 266 CalcRadius(geom, evt, xtmp1,ytmp1,rad,dev); 267 CalcRadius(geom, evt, xtmp2,ytmp2,rad2,dev2); 268 if(dev2<dev){ 269 xtmp1=xtmp2; ytmp1=ytmp2; dev=dev2; rad=rad2; 270 } 271 272 // find the best fit. 273 CalcMinimumDeviation(geom, evt, xtmp1,ytmp1,hillas.GetMeanX(), 274 hillas.GetMeanY(), dev, opt_rad, new_sigma, 275 newx, newy); 276 277 fRadius = opt_rad; 278 fDeviation = new_sigma; 279 280 fCenterX = newx; 281 fCenterY = newy; 282 283 //SetReadyToSave(); 228 void MMuonSearchPar::Calc(const MGeomCam &geom, const MSignalCam &evt, 229 const MHillas &hillas) 230 { 231 Double_t x = hillas.GetMeanX(); 232 Double_t y = hillas.GetMeanY(); 233 234 // ------------------------------------------------- 235 // Keiichi suggested trying to precalculate the Muon 236 // center a bit better, but it neither improves the 237 // fit result nor the speed 238 // 239 // const Float_t tmpr = 300.; // assume that the temporal cherenkov angle is 1 deg. (300 mm). 240 // 241 // const Double_t a = TMath::Tan(hillas.GetDelta()); 242 // 243 // const Double_t dx = a/TMath::Sqrt(tmpr+a*a)/3.; 244 // const Double_t dy = -tmpr/TMath::Sqrt(1+a*a)/3.; 245 // 246 // Double_t par1[] = { x+dx, y+dy, 0, 0 }; 247 // Double_t par2[] = { x-dx, y-dy, 0, 0 }; 248 // 249 // const Double_t dev1 = MMuonSearchPar::Fcn(par1); 250 // const Double_t dev2 = MMuonSearchPar::Fcn(par2); 251 // 252 // if (dev1<dev2) 253 // { 254 // x += dx; 255 // y += dy; 256 // } 257 // else 258 // { 259 // x -= dx; 260 // y -= dy; 261 // } 262 // ------------------------------------------------- 263 264 Double_t sigma, rad; 265 266 // find the best fit. 267 CalcMinimumDeviation(geom, evt, x, y, sigma, rad); 268 269 fCenterX = x; 270 fCenterY = y; 271 fRadius = rad; 272 fDeviation = sigma; 273 274 //SetReadyToSave(); 284 275 } 285 276 … … 298 289 *fLog << all; 299 290 *fLog << "Muon Parameters (" << GetName() << ")" << endl; 300 *fLog << " - Est. Radius [deg.] = " << fRadius*geom.GetConvMm2Deg() << endl; 301 *fLog << " - Deviation [deg.] = " << fDeviation*geom.GetConvMm2Deg() << endl; 302 *fLog << " - Center Pos. X [deg.] = " << fCenterX*geom.GetConvMm2Deg() << endl; 303 *fLog << " - Center Pos. Y [deg.] = " << fCenterY*geom.GetConvMm2Deg() << endl; 304 } 305 306 307 291 *fLog << " - Est. Radius [deg] = " << fRadius*geom.GetConvMm2Deg() << endl; 292 *fLog << " - Deviation [deg] = " << fDeviation*geom.GetConvMm2Deg() << endl; 293 *fLog << " - Center Pos. X [deg] = " << fCenterX*geom.GetConvMm2Deg() << endl; 294 *fLog << " - Center Pos. Y [deg] = " << fCenterY*geom.GetConvMm2Deg() << endl; 295 } -
trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h
r6973 r6979 4 4 #ifndef MARS_MParContainer 5 5 #include "MParContainer.h" 6 #endif 7 8 #ifndef MARS_MArrayF 9 #include "MArrayF.h" 6 10 #endif 7 11 … … 18 22 Float_t fCenterY; // An estimated center position in Y of the muon ring [mm] 19 23 24 MArrayF fSignal; //! Temporary storage for signal 25 MArrayF fX; //! Temporary storage for pixels X-position 26 MArrayF fY; //! Temporary storage for pixels Y-position 27 28 static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); 29 Double_t Fcn(Double_t *par) const; 30 20 31 public: 21 32 MMuonSearchPar(const char *name=NULL, const char *title=NULL); … … 28 39 Float_t GetCenterY() const { return fCenterY; } 29 40 30 void CalcTempCenter(const MHillas &hillas, Float_t &xtmp1,31 Float_t &ytmp1, Float_t &xtmp2, Float_t &ytmp2);32 Bool_t CalcRadius(const MGeomCam &geom, const MSignalCam &evt, Float_t x,33 Float_t y, Float_t &r, Float_t &sigma);34 41 void CalcMinimumDeviation(const MGeomCam &geom, const MSignalCam &evt, 35 Float_t x, Float_t y, Float_t xcog, 36 Float_t ycog, Float_t sigma, Float_t &opt_rad, 37 Float_t &new_sigma, Float_t &newx, Float_t &newy); 42 Double_t &x, Double_t &y, Double_t &sigma, Double_t &rad); 43 38 44 void Calc(const MGeomCam &geom, const MSignalCam &evt, 39 45 const MHillas &hillas); -
trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.cc
r6973 r6979 17 17 ! 18 18 ! Author(s): Keiichi Mase 09/2004 <mailto:mase@mppmu.mpg.de> 19 ! 19 ! Author(s): Markus Meyer 09/2004 <mailto:meyer@astro.uni-wuerzburg.de> 20 20 ! 21 ! Copyright: MAGIC Software Development, 2000-200 421 ! Copyright: MAGIC Software Development, 2000-2005 22 22 ! 23 23 ! … … 41 41 // 42 42 // Input Containers: 43 // [MGeomCam]44 // [MHillas]45 // [MSignalCam]43 // MGeomCam 44 // MHillas 45 // MSignalCam 46 46 // 47 47 // Output Containers: 48 // [MMuonSearchPar]48 // MMuonSearchPar 49 49 // 50 50 ////////////////////////////////////////////////////////////////////////////// 51 51 #include "MMuonSearchParCalc.h" 52 52 53 #include <fstream>53 #include "MParList.h" 54 54 55 #include "MParList.h" 55 #include "MLog.h" 56 #include "MLogManip.h" 57 56 58 #include "MGeomCam.h" 57 59 #include "MSignalCam.h" 58 60 #include "MMuonSearchPar.h" 59 #include "MLog.h"60 #include "MLogManip.h"61 61 62 62 using namespace std; … … 71 71 // Default constructor. 72 72 // 73 MMuonSearchParCalc::MMuonSearchParCalc 74 (const char *mupar, const char *name, const char *title) 73 MMuonSearchParCalc::MMuonSearchParCalc(const char *name, const char *title) 75 74 : fHillas(NULL), fMuonPar(NULL) 76 75 { 77 76 fName = name ? name : gsDefName.Data(); 78 77 fTitle = title ? title : gsDefTitle.Data(); 79 80 fMuparName = mupar;81 fHillasInput = "MHillas";82 78 } 83 79 … … 86 82 Int_t MMuonSearchParCalc::PreProcess(MParList *pList) 87 83 { 88 fHillas = (MHillas*)pList->FindObject( fHillasInput,"MHillas");84 fHillas = (MHillas*)pList->FindObject("MHillas"); 89 85 if (!fHillas) 90 86 { 91 *fLog << err << fHillasInput << " [MHillas]not found... aborting." << endl;87 *fLog << err << "MHillas not found... aborting." << endl; 92 88 return kFALSE; 93 89 } … … 107 103 } 108 104 109 fMuonPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar" , fMuparName);105 fMuonPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar"); 110 106 if (!fMuonPar) 111 107 return kFALSE; … … 118 114 Int_t MMuonSearchParCalc::Process() 119 115 { 120 121 fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas); 122 123 return kTRUE; 116 fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas); 117 return kTRUE; 124 118 } 125 -
trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.h
r6973 r6979 14 14 { 15 15 private: 16 MGeomCam *fGeomCam; 17 MSignalCam *fSignalCam; 18 MHillas *fHillas; //! Pointer to the source independent hillas parameters 19 MMuonSearchPar *fMuonPar; //! Pointer to the output container for the new image parameters 20 21 TString fMuparName; 22 TString fHillasInput; 16 MGeomCam *fGeomCam; //! 17 MSignalCam *fSignalCam; //! 18 MHillas *fHillas; //! Pointer to the source independent hillas parameters 19 MMuonSearchPar *fMuonPar; //! Pointer to the output container for the new image parameters 23 20 24 21 Int_t PreProcess(MParList *plist); … … 26 23 27 24 public: 28 MMuonSearchParCalc(const char *mupar="MMuonSearchPar", 29 const char *name=NULL, const char *title=NULL); 30 31 void SetInput(TString hilname) { fHillasInput = hilname; } 25 MMuonSearchParCalc(const char *name=NULL, const char *title=NULL); 32 26 33 27 ClassDef(MMuonSearchParCalc, 0) // task to calculate muon parameters -
trunk/MagicSoft/Mars/mmuon/MMuonSetup.cc
r6973 r6979 16 16 ! 17 17 ! 18 ! Author(s): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de> 19 ! Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de> 18 ! Author(s): Markus Meyer 04/2005 <mailto:meyer@astro.uni-wuerzburg.de> 20 19 ! 21 ! Copyright: MAGIC Software Development, 2000-200 420 ! Copyright: MAGIC Software Development, 2000-2005 22 21 ! 23 22 ! … … 28 27 // MMuonSetup 29 28 // 30 // Storage Container for parameter for the muon analysis 31 // 32 // 33 // Input Containers: 34 // [MSignalCam] 29 // Storage Container for setup parameter for the muon analysis 35 30 // 36 31 ///////////////////////////////////////////////////////////////////////////// … … 50 45 // 51 46 MMuonSetup::MMuonSetup(const char *name, const char *title) 47 : fMargin(0.2), fThresholdArcPhi(0.1), fThresholdArcWidth(0.1) 52 48 { 53 49 fName = name ? name : "MMuonSetup"; 54 fTitle = title ? title : "Muon calibration parameters"; 55 56 fMargin = 60.; // in mm 57 fArcPhiThres = 30.; 58 fArcWidthThres = 30.; 59 fArcPhiBinNum = 20; 60 fArcPhiHistStartVal = -180.; // deg. 61 fArcPhiHistEndVal = 180.; // deg. 62 fArcWidthBinNum = 28; 63 fArcWidthHistStartVal = 0.3; // deg. 64 fArcWidthHistEndVal = 1.7; // deg. 50 fTitle = title ? title : "Muon calibration setup parameters"; 65 51 } 66 52 67 53 // -------------------------------------------------------------------------- 68 54 // 69 MMuonSetup::~MMuonSetup() 55 // Read the setup from a TEnv, eg: 56 // MMuonSetup.Margin: 0.2 57 // MMuonSetup.ThresholdArcPhi: 0.1 58 // MMuonSetup.ThresholdArcWidth: 0.1 59 // 60 Int_t MMuonSetup::ReadEnv(const TEnv &env, TString prefix, Bool_t print) 70 61 { 62 Bool_t rc = kFALSE; 63 if (IsEnvDefined(env, prefix, "Margin", print)) 64 { 65 rc = kTRUE; 66 fMargin = GetEnvValue(env, prefix, "Margin", fMargin); 67 } 68 if (IsEnvDefined(env, prefix, "ThresholdArcPhi", print)) 69 { 70 rc = kTRUE; 71 fThresholdArcPhi = GetEnvValue(env, prefix, "ThresholdArcPhi", fThresholdArcPhi); 72 } 73 if (IsEnvDefined(env, prefix, "ThresholdArcWidth", print)) 74 { 75 rc = kTRUE; 76 fThresholdArcWidth = GetEnvValue(env, prefix, "ThresholdArcWidth", fThresholdArcWidth); 77 } 78 79 return rc; 71 80 } -
trunk/MagicSoft/Mars/mmuon/MMuonSetup.h
r6973 r6979 9 9 { 10 10 private: 11 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 SetMargin12 Float_t fArcPhiThres; //The threshold value to define arc phi13 Float_t fArcWidthThres; //The threshold value to define arc width11 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 12 Float_t fThresholdArcPhi; // [deg] The threshold value to define arc phi 13 Float_t fThresholdArcWidth; // [deg] The threshold value to define arc width 14 14 15 Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); 15 16 16 17 public: 17 MMuonSetup(const char *name=NULL, const char *title=NULL); 18 ~MMuonSetup(); 19 20 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. 21 Int_t fArcWidthBinNum; // The bin number for the histogram of arc wid 22 Float_t fArcPhiHistStartVal; // The starting value for the histogram of arc phi 23 Float_t fArcPhiHistEndVal; // The end value for the histogram of arc phi 24 Float_t fArcWidthHistStartVal; // The starting value for the histogram of arc width 25 Float_t fArcWidthHistEndVal; // The end value for the histogram of arc width 18 MMuonSetup(const char *name=NULL, const char *title=NULL); 26 19 27 Float_t GetMargin() const { return fMargin; } 28 Float_t GetArcPhiThres() const { return fArcPhiThres; } 29 Float_t GetArcWidthThres() const { return fArcWidthThres; } 30 Float_t GetArcPhiBinNum() const { return fArcPhiBinNum; } 31 Float_t GetArcWidthBinNum() const { return fArcWidthBinNum; } 32 33 void SetMargin(Float_t margin) { fMargin = margin; } 34 void SetArcPhiThres(Float_t thres) { fArcPhiThres = thres; } 35 void SetArcWidthThres(Float_t thres) { fArcWidthThres = thres; } 36 void SetArcPhiBinNum(Int_t num) { fArcPhiBinNum = num; } 37 void SetArcWidthBinNum(Int_t num) { fArcWidthBinNum = num; } 38 39 ClassDef(MMuonSetup, 1) // Container to hold muon parameters for the whole sample 20 // Getter 21 Float_t GetMargin() const { return fMargin; } 22 Float_t GetThresholdArcPhi() const { return fThresholdArcPhi; } 23 Float_t GetThresholdArcWidth() const { return fThresholdArcWidth; } 24 25 // Setter 26 void SetMargin(Float_t margin) { fMargin = margin; } 27 void SetThresholdArcPhi(Float_t thres) { fThresholdArcPhi = thres; } 28 void SetThresholdArcWidth(Float_t thres) { fThresholdArcWidth = thres; } 29 30 ClassDef(MMuonSetup, 1) // Container to hold setup for muon analysis 40 31 }; 41 32 -
trunk/MagicSoft/Mars/mmuon/Makefile
r6973 r6979 30 30 MMuonCalibParCalc.cc \ 31 31 MHMuonPar.cc \ 32 MHSingleMuon.cc \ 32 33 MMuonSetup.cc 33 34 -
trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h
r6973 r6979 6 6 7 7 #pragma link C++ class MMuonSearchPar+; 8 #pragma link C++ class MMuonSearchParCalc+;9 8 #pragma link C++ class MMuonCalibPar+; 10 #pragma link C++ class MMuonCalibParCalc+;11 #pragma link C++ class MHMuonPar+;12 9 #pragma link C++ class MMuonSetup+; 13 10 11 #pragma link C++ class MMuonSearchParCalc+; 12 #pragma link C++ class MMuonCalibParCalc+; 13 14 #pragma link C++ class MHMuonPar+; 15 #pragma link C++ class MHSingleMuon+; 16 14 17 #endif 15 16 17 18 19 20 21
Note:
See TracChangeset
for help on using the changeset viewer.