Changeset 3112
- Timestamp:
- 02/12/04 11:34:32 (21 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r3111 r3112 4 4 5 5 -*-*- END OF LINE -*-*- 6 2004/02/12: Thomas Bretz 6 7 8 * mbase/MArray.[h,cc]: 9 - changed argument of CutEdges from pointer to reference 10 11 * mgeom/MGeomPix.h: 12 - fixed calculation of area of pixel. It was too big for a factor 13 or 2 14 15 * mjobs/MJCalibration.cc: 16 - slight change of name of a MHCamera 17 18 * manalysis/MExtractSignal.cc: 19 - slight change to the creation of the satpixels list 20 21 * mcalib/MHCalibrationBlindPixel.cc, mcalib/MHCalibrationPixel.cc, 22 mcalib/MHGausEvent.cc: 23 - corrected call to ProjectArray 24 - corrected call to CutEdges 25 26 * mfilter/MFCosmics.[h,cc]: 27 - small fixes to logging output 28 - small simplification to return statement 29 - declared CosmicsRejection const 30 31 * mhbase/MH.[h,cc]: 32 - changed argument of ProjectArray from pointer to reference 33 - added missing calcualtion of minimum 34 - removed obsolete SetEntries 35 - changed SetDirectory from NULL to gROOT 36 37 38 7 39 2004/02/12: Javier López 8 40 9 41 * macros/pointspreadfunction.C 42 10 43 - added new macro that fits with a 2D gaussian the DC spot for an 11 44 star. It gives you the RMS of the PSF and the center of the star, 12 45 very useful for misspointing studies. 46 47 13 48 14 49 2004/02/11: Hendrik Bartko -
trunk/MagicSoft/Mars/manalysis/MExtractSignal.cc
r3094 r3112 163 163 164 164 if (satlo) 165 165 { 166 166 sat++; 167 satpixels += Form(" %d%s",pixel.GetPixelId()," ");168 167 satpixels += Form(" %d", pixel.GetPixelId()); 168 } 169 169 170 170 … … 186 186 187 187 if (sat) 188 *fLog << warn << "WARNING - Lo Gain saturated in " << sat << " pixels:" << satpixels << endl;188 *fLog << warn << "WARNING - Lo Gain saturated in " << sat << " pixels with Index:" << satpixels << endl; 189 189 190 190 fSignals->SetReadyToSave(); -
trunk/MagicSoft/Mars/mbase/MArray.cc
r3105 r3112 47 47 // Cuts the last entries of an array containing only zeros. 48 48 // 49 void MArray::CutEdges(TArrayD *arr)49 void MArray::CutEdges(TArrayD &arr) 50 50 { 51 52 Int_t i; 53 54 for (i=arr->GetSize()-1;i>=0;i--) 55 if (arr->At(i) != 0) 56 { 57 arr->Set(i+1); 58 break; 59 } 51 const Int_t n = arr.GetSize(); 52 53 for (Int_t i=n-1; i>=0; i--) 54 if (arr[i] != 0) 55 { 56 arr.Set(i+1); 57 break; 58 } 60 59 } 61 60 … … 64 63 // Cuts the last entries of an array containing only zeros. 65 64 // 66 void MArray::CutEdges(TArrayF *arr)65 void MArray::CutEdges(TArrayF &arr) 67 66 { 68 69 Int_t i; 70 71 for (i=arr->GetSize()-1;i>=0;i--) 72 if (arr->At(i) != 0) 73 { 74 arr->Set(i+1); 75 break; 76 } 67 const Int_t n = arr.GetSize(); 68 69 for (Int_t i=n-1; i>=0; i--) 70 if (arr[i] != 0) 71 { 72 arr.Set(i+1); 73 break; 74 } 77 75 } 78 -
trunk/MagicSoft/Mars/mbase/MArray.h
r3105 r3112 32 32 virtual void Set(UInt_t n) = 0; 33 33 34 static void CutEdges(TArrayF *arr); 35 static void CutEdges(TArrayD *arr); 36 34 static void CutEdges(TArrayF &arr); 35 static void CutEdges(TArrayD &arr); 37 36 38 37 ClassDef(MArray, 1) //Abstract array base class for TObject derived Arrays -
trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc
r3108 r3112 268 268 fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains); 269 269 270 fHPSD = MH::ProjectArray(fPSDHiGain, fPSDNbins,271 272 270 fHPSD = ProjectArray(*fPSDHiGain, fPSDNbins, 271 "PSDProjectionBlindPixel", 272 "Power Spectrum Density Projection HiGain Blind Pixel"); 273 273 274 274 // -
trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc
r3108 r3112 269 269 270 270 if (IsUseLoGain()) 271 fHPSD = MH::ProjectArray(fPSDLoGain, fPSDNbins,272 273 271 fHPSD = ProjectArray(*fPSDLoGain, fPSDNbins, 272 Form("%s%d","PSDProjection",fPixId), 273 Form("%s%d","Power Spectrum Density Projection LoGain ",fPixId)); 274 274 else 275 fHPSD = MH::ProjectArray(fPSDHiGain, fPSDNbins,276 277 275 fHPSD = ProjectArray(*fPSDHiGain, fPSDNbins, 276 Form("%s%d","PSDProjection",fPixId), 277 Form("%s%d","Power Spectrum Density Projection HiGain ",fPixId)); 278 278 279 279 // -
trunk/MagicSoft/Mars/mcalib/MHGausEvent.cc
r3108 r3112 177 177 178 178 // This cuts only the non-used zero's, but MFFT will later cut the rest 179 MArray::CutEdges( fEvents);179 MArray::CutEdges(*fEvents); 180 180 181 181 MFFT fourier; 182 182 183 fPowerSpectrum = fourier.PowerSpectrumDensity(fEvents); 184 185 fHPowerProbability = MH::ProjectArray(fPowerSpectrum, fPowerProbabilityBins, 186 "PowerProbability", 187 "Probability of Power occurrance"); 183 fPowerSpectrum = fourier.PowerSpectrumDensity(fEvents); 184 fHPowerProbability = ProjectArray(*fPowerSpectrum, fPowerProbabilityBins, 185 "PowerProbability", 186 "Probability of Power occurrance"); 188 187 // 189 188 // First guesses for the fit (should be as close to reality as possible, -
trunk/MagicSoft/Mars/mfilter/MFCosmics.cc
r3105 r3112 48 48 // 49 49 // Output Containers: 50 // -/- 50 51 // 51 52 ////////////////////////////////////////////////////////////////////////////// … … 77 78 fRawEvt(NULL), fMaxEmptyPixels(230) 78 79 { 79 80 80 fName = name ? name : "MFCosmics"; 81 81 fTitle = title ? title : "Filter to reject cosmics"; … … 91 91 Int_t MFCosmics::PreProcess(MParList *pList) 92 92 { 93 94 93 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); 95 94 if (!fRawEvt) 96 95 { 97 *fLog << err << dbginf <<"MRawEvtData not found... aborting." << endl;96 *fLog << err << "MRawEvtData not found... aborting." << endl; 98 97 return kFALSE; 99 98 } … … 101 100 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); 102 101 if (!fPedestals) 103 104 *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;102 { 103 *fLog << err << "MPedestalCam not found... aborting." << endl; 105 104 return kFALSE; 106 105 } 107 106 108 107 fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam"); 109 108 if (!fSignals) 110 111 *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;109 { 110 *fLog << err << "MExtractedSignalCam not found... aborting." << endl; 112 111 return kFALSE; 113 112 } 114 113 115 114 memset(fCut, 0, sizeof(fCut)); … … 124 123 // - MExtractedSignalCam 125 124 // 126 Bool_t MFCosmics::ReInit(MParList *pList ) 127 { 128 125 Bool_t MFCosmics::ReInit(MParList *pList) 126 { 129 127 fSqrtHiGainSamples = TMath::Sqrt((Float_t) fSignals->GetNumUsedHiGainFADCSlices()); 130 128 … … 140 138 Int_t MFCosmics::Process() 141 139 { 142 143 fResult = CosmicsRejection(); 144 145 fCut[fResult ? 0 : 1]++; 146 return kTRUE; 140 fResult = CosmicsRejection(); 141 142 fCut[fResult ? 0 : 1]++; 143 return kTRUE; 147 144 } 148 145 … … 158 155 // the outer pixels have some defect). 159 156 // 160 Bool_t MFCosmics::CosmicsRejection() 161 { 162 163 MRawEvtPixelIter pixel(fRawEvt); 164 165 Int_t cosmicpix = 0; 166 167 // 168 // Create a first loop to sort out the cosmics ... 169 // 170 while (pixel.Next()) 171 { 172 173 const UInt_t idx = pixel.GetPixelId(); 174 175 MExtractedSignalPix &sig = (*fSignals)[idx]; 176 MPedestalPix &ped = (*fPedestals)[idx]; 177 const Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples; 178 const Float_t sumhi = sig.GetExtractedSignalHiGain(); 179 180 // 181 // We consider a pixel as presumably due to cosmics 182 // if its sum of FADC slices is lower than 3 pedestal RMS 183 // 184 if (sumhi < 3.*pedrms ) 185 cosmicpix++; 186 } 187 188 189 // 190 // If the camera contains more than fMaxEmptyPixels 191 // presumed pixels due to cosmics, then the event is discarted. 192 // 193 if (cosmicpix > fMaxEmptyPixels) 194 return kTRUE; 195 196 return kFALSE; 157 Bool_t MFCosmics::CosmicsRejection() const 158 { 159 MRawEvtPixelIter pixel(fRawEvt); 160 161 Int_t cosmicpix = 0; 162 163 // 164 // Create a first loop to sort out the cosmics ... 165 // 166 while (pixel.Next()) 167 { 168 const UInt_t idx = pixel.GetPixelId(); 169 170 MExtractedSignalPix &sig = (*fSignals)[idx]; 171 MPedestalPix &ped = (*fPedestals)[idx]; 172 173 const Float_t pedrms = ped.GetPedestalRms()*fSqrtHiGainSamples; 174 const Float_t sumhi = sig.GetExtractedSignalHiGain(); 175 176 // 177 // We consider a pixel as presumably due to cosmics 178 // if its sum of FADC slices is lower than 3 pedestal RMS 179 // 180 if (sumhi < 3.*pedrms ) 181 cosmicpix++; 182 } 183 184 // 185 // If the camera contains more than fMaxEmptyPixels 186 // presumed pixels due to cosmics, then the event is discarted. 187 // 188 return cosmicpix > fMaxEmptyPixels; 197 189 } 198 190 199 191 Int_t MFCosmics::PostProcess() 200 192 { 201 202 if (GetNumExecutions()==0) 203 return kTRUE; 204 205 *fLog << inf << endl; 206 *fLog << GetDescriptor() << " execution statistics:" << endl; 207 *fLog << dec << setfill(' '); 208 209 *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ; 210 *fLog << (int)(fCut[1]*100/GetNumExecutions()) ; 211 *fLog << "%) Evts skipped due to: Cosmics Rejection applied " ; 212 *fLog << " (with fMaxEmptyPixels = " << fMaxEmptyPixels << ")" << endl; 213 214 *fLog << " " << setw(7) << fCut[0] << " (" << setw(3) ; 215 *fLog << (int)(fCut[0]*100/GetNumExecutions()) ; 216 *fLog << "%) Evts survived the cosmics rejection!" << endl; 217 *fLog << endl; 218 219 return kTRUE; 220 } 221 193 if (GetNumExecutions()==0) 194 return kTRUE; 195 196 *fLog << inf << endl; 197 *fLog << GetDescriptor() << " execution statistics:" << endl; 198 *fLog << dec << setfill(' '); 199 200 *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ; 201 *fLog << (int)(fCut[1]*100/GetNumExecutions()) ; 202 *fLog << "%) Evts skipped due to: Cosmics Rejection applied " ; 203 *fLog << " (with fMaxEmptyPixels = " << fMaxEmptyPixels << ")" << endl; 204 205 *fLog << " " << setw(7) << fCut[0] << " (" << setw(3) ; 206 *fLog << (int)(fCut[0]*100/GetNumExecutions()) ; 207 *fLog << "%) Evts survived the cosmics rejection!" << endl; 208 *fLog << endl; 209 210 return kTRUE; 211 } 212 -
trunk/MagicSoft/Mars/mfilter/MFCosmics.h
r3084 r3112 1 1 #ifndef MARS_MFCosmics 2 2 #define MARS_MFCosmics 3 4 /////////////////////////////////////////////////////////////////////////////5 // //6 // MFCosmics //7 // //8 // Filter against cosmics (used especially by the calibration //9 // //10 /////////////////////////////////////////////////////////////////////////////11 3 12 4 #ifndef MARS_MFilter … … 22 14 { 23 15 private: 16 MPedestalCam *fPedestals; // Pedestals of all pixels in the camera 17 MExtractedSignalCam *fSignals; // Calibration events of all pixels in the camera 24 18 25 MPedestalCam *fPedestals; // Pedestals of all pixels in the camera 26 MExtractedSignalCam *fSignals; // Calibration events of all pixels in the camera 19 MRawEvtData *fRawEvt; // raw event data (time slices) 27 20 28 MRawEvtData *fRawEvt; // raw event data (time slices) 21 Int_t fCut[2]; 22 Bool_t fResult; 29 23 30 Int_t fCut[2]; 31 Bool_t fResult; 32 Int_t fMaxEmptyPixels; // Maximum number of empty pixels before declaring event as cosmic 33 Float_t fSqrtHiGainSamples; // Square root of the number of used Hi-Gain Samples 34 35 Bool_t ReInit(MParList *pList); 36 Int_t PreProcess(MParList *pList); 37 Int_t Process(); 38 Int_t PostProcess(); 24 Int_t fMaxEmptyPixels; // Maximum number of empty pixels before declaring event as cosmic 25 Float_t fSqrtHiGainSamples; // Square root of the number of used Hi-Gain Samples 39 26 40 Bool_t CosmicsRejection(); 41 42 Bool_t IsExpressionTrue() const { return fResult; } 27 Bool_t ReInit(MParList *pList); 28 Int_t PreProcess(MParList *pList); 29 Int_t Process(); 30 Int_t PostProcess(); 31 32 Bool_t CosmicsRejection() const; 33 34 Bool_t IsExpressionTrue() const { return fResult; } 43 35 44 36 public: 37 MFCosmics(const char *name=NULL, const char *title=NULL); 45 38 46 MFCosmics(const char *name=NULL, const char *title=NULL); 39 void SetMaxEmptyPixels(const Int_t n) { fMaxEmptyPixels = n; } 40 Int_t GetMaxEmptyPixels() const { return fMaxEmptyPixels; } 47 41 48 void SetMaxEmptyPixels(const Int_t n) { fMaxEmptyPixels = n; } 49 Int_t GetMaxEmptyPixels() const { return fMaxEmptyPixels; } 50 51 ClassDef(MFCosmics, 0) // Filter to perform a cosmics rejection 42 ClassDef(MFCosmics, 0) // Filter to perform a cosmics rejection 52 43 }; 53 44 -
trunk/MagicSoft/Mars/mgeom/MGeomPix.h
r2463 r3112 34 34 void Print(Option_t *opt=NULL) const; 35 35 36 void Set(Float_t x, Float_t y, Float_t d, UInt_t s=0) { fX=x; fY=y; fD=d; fA=d*d*gsTan60 ; fSector=s; }36 void Set(Float_t x, Float_t y, Float_t d, UInt_t s=0) { fX=x; fY=y; fD=d; fA=d*d*gsTan60/2; fSector=s; } 37 37 38 38 void SetNeighbors(Short_t i0=-1, Short_t i1=-1, Short_t i2=-1, … … 46 46 UInt_t GetSector() const { return fSector; } 47 47 48 Float_t GetA() const { return fA; /*fD*fD*gsTan60 ;*/ }48 Float_t GetA() const { return fA; /*fD*fD*gsTan60/2;*/ } 49 49 50 50 Byte_t GetNumNeighbors() const { return fNumNeighbors; } -
trunk/MagicSoft/Mars/mhbase/MH.cc
r3081 r3112 1082 1082 } 1083 1083 1084 TH1I* MH::ProjectArray(const TArrayF *array, Int_t nbins, const char* name, const char* title) 1085 { 1086 1087 const Int_t size = array->GetSize(); 1088 TH1I *h1=0; 1089 1090 //check if histogram with identical name exist 1091 TObject *h1obj = gROOT->FindObject(name); 1092 if (h1obj && h1obj->InheritsFrom("TH1I")) { 1093 h1 = (TH1I*)h1obj; 1094 h1->Reset(); 1095 } 1096 1097 Double_t max = 0.; 1098 Double_t min = 0.; 1099 Int_t newbins = 0; 1100 1101 // first loop over array to find the maximum: 1102 for (Int_t i=0; i<size;i++) 1103 if (array->At(i) > max) 1104 max = array->At(i); 1105 1106 FindGoodLimits(nbins, newbins, min, max, kFALSE); 1107 1108 if (!h1) 1109 { 1110 h1 = new TH1I(name, title, newbins, min, max); 1111 h1->SetXTitle(""); 1112 h1->SetYTitle("Counts"); 1113 h1->SetDirectory(NULL); 1114 } 1115 1084 // -------------------------------------------------------------------------- 1085 // 1086 // M.Gaug added this withouz Documentation 1087 // 1088 TH1I* MH::ProjectArray(const TArrayF &array, Int_t nbins, const char* name, const char* title) 1089 { 1090 const Int_t size = array.GetSize(); 1091 1092 TH1I *h1=0; 1093 1094 //check if histogram with identical name exist 1095 TObject *h1obj = gROOT->FindObject(name); 1096 if (h1obj && h1obj->InheritsFrom("TH1I")) 1097 { 1098 h1 = (TH1I*)h1obj; 1099 h1->Reset(); 1100 } 1101 1102 Double_t min = size>0 ? array[0] : 0; 1103 Double_t max = size>0 ? array[0] : 1; 1104 1105 // first loop over array to find the min and max 1106 for (Int_t i=1; i<size;i++) 1107 { 1108 max = TMath::Max((Double_t)array[i], max); 1109 min = TMath::Min((Double_t)array[i], min); 1110 } 1111 1112 Int_t newbins = 0; 1113 FindGoodLimits(nbins, newbins, min, max, kFALSE); 1114 1115 if (!h1) 1116 { 1117 h1 = new TH1I(name, title, newbins, min, max); 1118 h1->SetXTitle(""); 1119 h1->SetYTitle("Counts"); 1120 h1->SetDirectory(gROOT); 1121 } 1122 1116 1123 // Second loop to fill the histogram 1117 for (Int_t i=0;i<size;i++) 1118 h1->Fill(array->At(i)); 1119 1120 h1->SetEntries(size); 1121 1124 for (Int_t i=0;i<size;i++) 1125 h1->Fill(array[i]); 1126 1122 1127 return h1; 1123 1128 } 1124 1129 1125 TH1I* MH::ProjectArray(const TArrayD *array, Int_t nbins, const char* name, const char* title) 1126 { 1127 1128 const Int_t size = array->GetSize(); 1129 1130 TH1I *h1=0; 1131 1132 //check if histogram with identical name exist 1133 TObject *h1obj = gROOT->FindObject(name); 1134 if (h1obj && h1obj->InheritsFrom("TH1I")) { 1135 h1 = (TH1I*)h1obj; 1136 h1->Reset(); 1137 } 1138 1139 Double_t max = 0.; 1140 Double_t min = 0.; 1141 Int_t newbins = 0; 1142 1143 // first loop over array to find the maximum: 1144 for (Int_t i=0; i<size;i++) 1145 if (array->At(i) > max) 1146 max = array->At(i); 1147 1148 FindGoodLimits( nbins, newbins, min, max, kFALSE); 1149 1150 if (!h1) 1151 { 1152 h1 = new TH1I(name, title, newbins, min, max); 1153 h1->SetXTitle(""); 1154 h1->SetYTitle("Counts"); 1155 h1->SetDirectory(NULL); 1156 } 1157 1130 // -------------------------------------------------------------------------- 1131 // 1132 // M.Gaug added this withouz Documentation 1133 // 1134 TH1I* MH::ProjectArray(const TArrayD &array, Int_t nbins, const char* name, const char* title) 1135 { 1136 const Int_t size = array.GetSize(); 1137 TH1I *h1=0; 1138 1139 //check if histogram with identical name exist 1140 TObject *h1obj = gROOT->FindObject(name); 1141 if (h1obj && h1obj->InheritsFrom("TH1I")) 1142 { 1143 h1 = (TH1I*)h1obj; 1144 h1->Reset(); 1145 } 1146 1147 Double_t min = size>0 ? array[0] : 0; 1148 Double_t max = size>0 ? array[0] : 1; 1149 1150 // first loop over array to find the min and max 1151 for (Int_t i=1; i<size;i++) 1152 { 1153 max = TMath::Max(array[i], max); 1154 min = TMath::Min(array[i], min); 1155 } 1156 1157 Int_t newbins = 0; 1158 FindGoodLimits(nbins, newbins, min, max, kFALSE); 1159 1160 if (!h1) 1161 { 1162 h1 = new TH1I(name, title, newbins, min, max); 1163 h1->SetXTitle(""); 1164 h1->SetYTitle("Counts"); 1165 h1->SetDirectory(gROOT); 1166 } 1167 1158 1168 // Second loop to fill the histogram 1159 for (Int_t i=0;i<size;i++) 1160 h1->Fill(array->At(i)); 1161 1162 h1->SetEntries(size); 1163 1169 for (Int_t i=0;i<size;i++) 1170 h1->Fill(array[i]); 1171 1164 1172 return h1; 1165 1173 } -
trunk/MagicSoft/Mars/mhbase/MH.h
r3081 r3112 92 92 static Int_t CutEdges(TH1 *h, Int_t nbins); 93 93 94 static TH1I* ProjectArray(const TArrayF *array, Int_t nbins=30, const char* name="ProjectArray",95 96 static TH1I* ProjectArray(const TArrayD *array, Int_t nbins=30, const char* name="ProjectArray",97 94 static TH1I* ProjectArray(const TArrayF &array, Int_t nbins=30, 95 const char* name="ProjectArray", const char* title="Projected Array"); 96 static TH1I* ProjectArray(const TArrayD &array, Int_t nbins=30, 97 const char* name="ProjectArray", const char* title="Projected Array"); 98 98 99 99 static void ProjectionX(TH1D &dest, const TH2 &src, Int_t firstybin=-1, Int_t lastybin=9999); -
trunk/MagicSoft/Mars/mjobs/MJCalibration.cc
r3085 r3112 223 223 MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted"); 224 224 MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results"); 225 MHCamera disp18 (geomcam, "Cal;Oscillati ng", "Oscillating Pixels");225 MHCamera disp18 (geomcam, "Cal;Oscillation", "Oscillating Pixels"); 226 226 MHCamera disp19 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain"); 227 227
Note:
See TracChangeset
for help on using the changeset viewer.