Changeset 4716 for trunk/MagicSoft
- Timestamp:
- 08/23/04 17:49:53 (20 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r4714 r4716 79 79 * mhist/MHEvent.[h,cc]: 80 80 - added new option for island index 81 - added kEvtCleaningData 81 82 82 83 * mimage/MImgCleanStd.[h,cc]: … … 85 86 - added ReadEnv 86 87 - updated StreamPrimitive 88 - added new cleaning method (probability cleaning) 87 89 88 90 * mimage/Makefile: … … 142 144 * msignal/MArrivalTime.h: 143 145 - removed operator() 146 - added operator[] const 147 148 * manalysis/MCameraData.[h,cc]: 149 - added algorithm for 'Probability cleaning' 150 151 * mbase/MMath.[h,cc]: 152 - added GaussProb 153 154 * mjobs/MSequence.h: 155 - added IsValid 144 156 145 157 -
trunk/MagicSoft/Mars/NEWS
r4710 r4716 25 25 - implemented new image parameters displaying the number of islands, 26 26 saturated hi-gain and lo-gain pixels (MImagePar, MHImagePar) 27 28 - event display in executable changed to support also calibrated files 29 (done with MJCalibrateSignal) 30 31 - added a cleaning which takes signal height _and_ arrival time into 32 account: probability cleaning (for more details see MImgCleanStd) 27 33 28 34 -
trunk/MagicSoft/Mars/manalysis/MCameraData.cc
r4454 r4716 34 34 #include "MCameraData.h" 35 35 36 #include "MMath.h" 37 38 #include "MLog.h" 39 #include "MLogManip.h" 40 41 #include "MGeomCam.h" 42 #include "MGeomPix.h" 43 44 #include "MPedPhotCam.h" 45 #include "MPedPhotPix.h" 46 36 47 #include "MCerPhotEvt.h" 37 48 #include "MCerPhotPix.h" 38 49 39 #include "MGeomCam.h"40 #include "MGeomPix.h"41 42 #include "MLog.h"43 #include "MLogManip.h"44 45 #include "MPedPhotCam.h"46 #include "MPedPhotPix.h"47 48 50 #include "MSigmabar.h" 51 #include "MArrivalTime.h" 49 52 50 53 ClassImp(MCameraData); … … 101 104 const Int_t n = geom.GetNumPixels(); 102 105 103 fData.Set(n); 104 fData.Reset(); 105 106 fValidity.Set(n); 107 fValidity.Reset(); 108 109 const Int_t entries = evt.GetNumPixels(); 110 111 // 112 // check the number of all pixels against the noise level and 113 // set them to 'unused' state if necessary 114 // 106 // Reset arrays 107 fData.Set(n); 108 fData.Reset(); 109 110 fValidity.Set(n); 111 fValidity.Reset(); 112 113 const Int_t entries = evt.GetNumPixels(); 114 115 // calculate cleaning levels 115 116 for (Int_t i=0; i<entries; i++) 116 117 { … … 149 150 const Int_t n = geom.GetNumPixels(); 150 151 151 fData.Set(n); 152 fData.Reset(); 153 154 fValidity.Set(n); 155 fValidity.Reset(); 156 157 const Int_t entries = evt.GetNumPixels(); 152 // reset arrays 153 fData.Set(n); 154 fData.Reset(); 155 156 fValidity.Set(n); 157 fValidity.Reset(); 158 159 // check validity of rms with area index 0 158 160 const Float_t anoise0 = cam.GetArea(0).GetRms(); 159 161 if (anoise0<=0) 160 162 return; 161 163 162 // 163 // check the number of all pixels against the noise level and 164 // set them to 'unused' state if necessary 165 // 164 // calculate cleaning levels 165 const Int_t entries = evt.GetNumPixels(); 166 166 for (Int_t i=0; i<entries; i++) 167 167 { … … 200 200 const Int_t n = geom.GetNumPixels(); 201 201 202 fData.Set(n); 203 fData.Reset(); 204 205 fValidity.Set(n); 206 fValidity.Reset(); 207 202 // reset arrays 203 fData.Set(n); 204 fData.Reset(); 205 206 fValidity.Set(n); 207 fValidity.Reset(); 208 209 // check validity of noise 208 210 if (noise<=0) 209 211 return; 210 212 211 const Int_t entries = evt.GetNumPixels(); 212 213 // 214 // check the number of all pixels against the noise level and 215 // set them to 'unused' state if necessary 216 // 213 // calculate cleaning levels 214 const Int_t entries = evt.GetNumPixels(); 217 215 for (Int_t i=0; i<entries; i++) 218 216 { … … 244 242 const Int_t n = geom.GetNumPixels(); 245 243 246 fData.Set(n); 247 fData.Reset(); 248 249 fValidity.Set(n); 250 fValidity.Reset(); 251 252 const Int_t entries = evt.GetNumPixels(); 253 const Float_t noise0 = cam.GetArea(0).GetRms(); 244 // reset arrays 245 fData.Set(n); 246 fData.Reset(); 247 248 fValidity.Set(n); 249 fValidity.Reset(); 250 251 // check validity of rms with area index 0 252 const Float_t noise0 = cam.GetArea(0).GetRms(); 254 253 if (noise0<=0) 255 254 return; 256 255 257 // 258 // check the number of all pixels against the noise level and 259 // set them to 'unused' state if necessary 260 // 256 // calculate cleaning levels 257 const Int_t entries = evt.GetNumPixels(); 261 258 for (Int_t i=0; i<entries; i++) 262 259 { … … 280 277 // -------------------------------------------------------------------------- 281 278 // 279 // Function to calculate the cleaning level for all pixels in a given event. 280 // The level is the probability that the signal is a real signal (taking 281 // the signal height and the fluctuation of the background into account) 282 // times one minus the probability that the signal is a background 283 // fluctuation (calculated from the time spread of arrival times 284 // around the pixel with the highest signal) 285 // 286 // FIXME: Should the check noise<=0 be replaced by MBadPixels? 287 // 288 void MCameraData::CalcCleaningProbability(const MCerPhotEvt &evt, const MPedPhotCam &pcam, 289 const MGeomCam &geom, const MArrivalTime *tcam) 290 { 291 const Int_t n = geom.GetNumPixels(); 292 293 // Reset arrays 294 fData.Set(n); 295 fData.Reset(); 296 297 fValidity.Set(n); 298 fValidity.Reset(); 299 300 // check validity of noise 301 const Float_t anoise0 = pcam.GetArea(0).GetRms(); 302 if (anoise0<=0) 303 return; 304 305 const Int_t entries = evt.GetNumPixels(); 306 307 // Find pixel with max signal 308 Int_t maxidx = 0; 309 if (tcam) 310 { 311 // Find pixel enty with maximum signal 312 for (Int_t i=0; i<entries; i++) 313 { 314 const Double_t s0 = evt[i].GetNumPhotons() * geom.GetPixRatio(i); 315 const Double_t s1 = evt[maxidx].GetNumPhotons() * geom.GetPixRatio(maxidx); 316 if (s0>s1) 317 maxidx = i; 318 } 319 320 // Get software index of pixel 321 maxidx = evt[maxidx].GetPixId(); 322 } 323 324 const Double_t timemean = tcam ? (*tcam)[maxidx] : 0; 325 const Double_t timerms = 0.75; //[slices] rms time spread around highest pixel 326 327 // calculate cleaning levels 328 for (Int_t i=0; i<entries; i++) 329 { 330 const MCerPhotPix &spix = evt[i]; 331 332 const Int_t idx = spix.GetPixId(); 333 const Float_t rms = pcam[idx].GetRms(); 334 if (rms<=0) // fData[idx]=0, fValidity[idx]=0 335 continue; 336 337 fValidity[idx]=1; 338 339 // get signal and arrival time 340 const UInt_t aidx = geom[idx].GetAidx(); 341 const Float_t ratio = pcam.GetArea(aidx).GetRms()/anoise0; 342 343 const Double_t signal = spix.GetNumPhotons() * geom.GetPixRatio(idx) * ratio / rms; 344 const Double_t time = tcam ? (*tcam)[idx] : 1; 345 346 // if signal<0 the probability is equal 0 347 if (signal<0) 348 continue; 349 350 // Just for convinience for easy readybility 351 const Double_t means = 0; 352 const Double_t meant = timemean; 353 354 const Double_t sigmas = rms; 355 const Double_t sigmat = timerms; 356 357 // Get probabilities 358 const Double_t psignal = MMath::GaussProb(signal, sigmas, means); 359 const Double_t pbckgnd = MMath::GaussProb(time, sigmat, meant); 360 361 // Set probability 362 fData[idx] = psignal*(1-pbckgnd); 363 fValidity[idx] = 1; 364 365 // Make sure, that we don't run into trouble because of 366 // a little numerical uncertanty 367 if (fData[idx]>1) 368 fData[idx]=1; 369 if (fData[idx]<0) 370 fData[idx]=0; 371 } 372 } 373 374 // -------------------------------------------------------------------------- 375 // 282 376 // Returns the contents of the pixel. 283 377 // -
trunk/MagicSoft/Mars/manalysis/MCameraData.h
r4452 r4716 19 19 class MCerPhotEvt; 20 20 class MPedPhotCam; 21 class MArrivalTime; 21 22 22 23 class MCameraData : public MParContainer, public MCamEvent … … 38 39 void CalcCleaningLevel(const MCerPhotEvt &evt, Double_t noise, 39 40 const MGeomCam &geom); 40 41 41 void CalcCleaningLevel2(const MCerPhotEvt &evt, const MPedPhotCam &fCam, 42 const MGeomCam &geom); 43 42 const MGeomCam &geom); 44 43 void CalcCleaningLevelDemocratic(const MCerPhotEvt &evt, const MPedPhotCam &cam, 45 44 const MGeomCam &geom); 45 void CalcCleaningProbability(const MCerPhotEvt &evt, const MPedPhotCam &pcam, 46 const MGeomCam &geom, const MArrivalTime *tcam); 46 47 47 48 const TArrayD &GetData() const { return fData; } -
trunk/MagicSoft/Mars/mbase/MMath.cc
r3666 r4716 47 47 const Double_t f = s+k*k*b; 48 48 49 return f==0 ? 0 : (s-b)/ TMath::Sqrt(f);49 return f==0 ? 0 : (s-b)/Sqrt(f); 50 50 } 51 51 … … 86 86 return -1; 87 87 88 const Double_t l = s* TMath::Log(s/sum*(alpha+1)/alpha);89 const Double_t m = b* TMath::Log(b/sum*(alpha+1) );88 const Double_t l = s*Log(s/sum*(alpha+1)/alpha); 89 const Double_t m = b*Log(b/sum*(alpha+1) ); 90 90 91 return l+m<0 ? -1 : TMath::Sqrt((l+m)*2);91 return l+m<0 ? -1 : Sqrt((l+m)*2); 92 92 } 93 93 … … 104 104 return 0; 105 105 106 return TMath::Sign(sig, s-b);106 return Sign(sig, s-b); 107 107 } 108 109 // -------------------------------------------------------------------------- 110 // 111 // Returns: 2/(sigma*sqrt(2))*integral[0,x](exp(-(x-mu)^2/(2*sigma^2))) 112 // 113 Double_t MMath::GaussProb(Double_t x, Double_t sigma, Double_t mean) 114 { 115 static const Double_t sqrt2 = Sqrt(2.); 116 return Erf((x-mean)/(sigma*sqrt2)); 117 } 118 -
trunk/MagicSoft/Mars/mbase/MMath.h
r3666 r4716 9 9 { 10 10 public: 11 static Double_t GaussProb(Double_t x, Double_t sigma, Double_t mean=0); 12 11 13 static Double_t Significance(Double_t s, Double_t b); 12 14 static Double_t SignificanceSym(Double_t s, Double_t b); -
trunk/MagicSoft/Mars/mhist/MHEvent.cc
r4702 r4716 152 152 fHist->SetYTitle("L"); 153 153 break; 154 case kEvtCleaningData: 155 fHist->SetName("CleanData"); 156 fHist->SetYTitle("L"); 157 break; 154 158 case kEvtIdxMax: 155 159 fHist->SetName("Max Slice Idx"); … … 220 224 fHist->SetLevels(lvl); 221 225 } 226 break; 227 case kEvtCleaningData: 228 fHist->SetCamContent(*event, 0); 222 229 break; 223 230 case kEvtIdxMax: -
trunk/MagicSoft/Mars/mhist/MHEvent.h
r4702 r4716 22 22 kEvtSignalRaw, kEvtSignalDensity, kEvtPedestal, 23 23 kEvtPedestalRMS, kEvtRelativeSignal, kEvtCleaningLevels, 24 kEvtCleaningData, 24 25 kEvtIdxMax, kEvtArrTime, kEvtTrigPix, kEvtIslandIndex 25 26 }; -
trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc
r4702 r4716 549 549 fData->CalcCleaningLevelDemocratic(*fEvt, *fPed, *fCam); 550 550 break; 551 /*552 551 case kProbability: 553 552 fData->CalcCleaningProbability(*fEvt, *fPed, *fCam, fTime); 554 break; */553 break; 555 554 default: 556 555 break; -
trunk/MagicSoft/Mars/mjobs/MSequence.h
r4695 r4716 43 43 void Print(Option_t *o="") const; 44 44 45 Bool_t IsValid() const { return fSequence!=(UInt_t)-1; } 46 45 47 void SetupPedRuns(MDirIter &iter, const char *path=0) const; 46 48 void SetupDatRuns(MDirIter &iter, const char *path=0) const; -
trunk/MagicSoft/Mars/msignal/MArrivalTime.h
r4714 r4716 40 40 const TArrayF &GetDataErr() const { return fDataErr; } 41 41 42 Double_t operator[](int i) { return fData[i]; } 42 Double_t operator[](int i) { return fData[i]; } 43 Double_t operator[](int i) const { return fData[i]; } 43 44 44 45 // Do not do such things! It is highly dangerous to use two very similar operators
Note:
See TracChangeset
for help on using the changeset viewer.