Changeset 8478 for trunk/MagicSoft/Mars/mcalib

Timestamp:

05/09/07 13:15:53 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mcalib

Files:

• MCalibrationChargeCalc.cc (modified) (13 diffs)
• MCalibrationChargeCalc.h (modified) (3 diffs)
• MCalibrationRelTimeCalc.cc (modified) (7 diffs)

trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc

@@ -1 +1 @@
 /* ======================================================================== *\
-! $Name: not supported by cvs2svn $:$Id: MCalibrationChargeCalc.cc,v 1.177 2007-04-27 10:04:46 tbretz Exp $
+! $Name: not supported by cvs2svn $:$Id: MCalibrationChargeCalc.cc,v 1.178 2007-05-09 12:15:52 tbretz Exp $
 ! --------------------------------------------------------------------------
 !
@@ -19 +19 @@
 !
 !   Author(s): Markus Gaug  02/2004 <mailto:markus@ifae.es>
+!   Author(s): Thomas Bretz <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Copyright: MAGIC Software Development, 2000-2007
 !
 !
@@ -213 +214 @@
 #include "MLogManip.h"
 
+#include "MMath.h"
+
 #include "MParList.h"
 
@@ -255 +258 @@
 const Float_t MCalibrationChargeCalc::fgPheErrUpperLimit       = 5.5;
 const Float_t MCalibrationChargeCalc::fgFFactorErrLimit        = 4.5;
-const Float_t MCalibrationChargeCalc::fgArrTimeRmsLimit        = 3.5;
+const Float_t MCalibrationChargeCalc::fgArrTimeRmsLimit        = 5.0;
 const Float_t MCalibrationChargeCalc::fgUnsuitablesLimit       = 0.1;
 const Float_t MCalibrationChargeCalc::fgUnreliablesLimit       = 0.3;
@@ -672 +675 @@
 }
 
+// ----------------------------------------------------------------------------------------------------
+//
+// Check for outliers. They are marked with
+// MBadPixelsPix::kFluctuatingArrivalTimes
+//
+// see also MCalibrationRelTimeCalc::FinalizeRelTimes
+//
+void MCalibrationChargeCalc::FinalizeAbsTimes()
+{
+    const Int_t npixels = fGeom->GetNumPixels();
+    const Int_t nareas  = fGeom->GetNumAreas();
+
+    // Create an array capable of holding all pixels
+    TArrayF arr(npixels);
+
+    for (Int_t aidx=0; aidx<nareas; aidx++)
+    {
+        Int_t n = 0;
+        for (Int_t i=0; i<npixels; i++)
+        {
+            // Check for this aidx only
+            if ((*fGeom)[i].GetAidx()!=aidx)
+                continue;
+
+            // check if pixel may not contain a valid value
+            if ((*fBadPixels)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
+                continue;
+
+            // check if it was excluded for some reason
+            const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
+            if (pix.IsExcluded())
+                continue;
+
+            // if TimePrecision is valid fill it into array
+            if (pix.GetAbsTimeRms()>0)
+                arr[n++] = pix.GetAbsTimeRms();
+        }
+
+        // Check the ratio of valid entries to the ratio of pixels
+        const Float_t ratio = 100*n/fGeom->GetNumPixWithAidx(aidx);
+        if (3*ratio<2)
+            *fLog << warn << "Area   " << setw(4) << aidx << ": Only " << ratio << "% pixels with valid time resolution found." << endl;
+
+        // Calculate median and median deviation
+        Double_t med;
+        const Double_t dev = MMath::MedianDev(n, arr.GetArray(), med);
+
+        // Calculate upper and lower limit
+        const Float_t lolim = TMath::Max(med-fArrTimeRmsLimit*dev, 0.);
+        const Float_t hilim = TMath::Max(med+fArrTimeRmsLimit*dev, 0.);
+
+        // Now find the outliers
+        for (Int_t i=0; i<npixels; i++)
+        {
+            // Search only within this aidx
+            if ((*fGeom)[i].GetAidx()!=aidx)
+                continue;
+
+            // skip pixels already known to be unsuitable
+            if ((*fBadPixels)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
+                continue;
+
+            // check if a pixel has been excluded. This
+            const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
+
+            // Check if time precision is valid (might be invalid
+            // for example in cae of empty histograms)
+            const Float_t res = pix.GetAbsTimeRms();
+            if (res<0) //FIXME!!! How does this happen?
+            {
+                *fLog << warn << "Pixel  " << setw(4) << i << ": Abs-time resolution could not be calculated." << endl;
+                (*fBadPixels)[i].SetUncalibrated(MBadPixelsPix::kDeviatingTimeResolution);
+                continue;
+            }
+
+            // Now compare to a lower and upper limit
+            if (res<=lolim || res>=hilim)
+            {
+                *fLog << warn << "Pixel  " << setw(4) << i << ": Deviating abs-time resolution: "
+                    << Form("%5.2f", res) << " out of range "
+                    << Form("[%4.2f,%4.2f]", lolim, hilim) << endl;
+
+                (*fBadPixels)[i].SetUncalibrated(MBadPixelsPix::kFluctuatingArrivalTimes);
+            }
+        }
+    }
+}
+
 // -----------------------------------------------------------------------
 //
@@ -723 +814 @@
   // First loop over pixels, call FinalizePedestals and FinalizeCharges
   //
-  Int_t   nvalid      = 0;
-
+  Int_t nvalid = 0;
   for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
-    {
-
-      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[pixid];
+  {
       //
       // Check if the pixel has been excluded from the fits
       //
+      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[pixid];
       if (pix.IsExcluded())
-        continue;
-
-      MPedestalPix  &ped = (*fPedestals)[pixid];
-      MBadPixelsPix &bad = (*fBadPixels)[pixid];
-
-      const Int_t aidx    = (*fGeom)[pixid].GetAidx();
-
-      FinalizePedestals(ped,pix,aidx);
-
-      if (FinalizeCharges(pix,bad,"pixel  "))
-        nvalid++;
-
-      FinalizeArrivalTimes(pix,bad,"pixel  ");
-    }
+          continue;
+
+      FinalizePedestals((*fPedestals)[pixid], pix, (*fGeom)[pixid].GetAidx());
+
+      if (FinalizeCharges(pix, (*fBadPixels)[pixid], "Pixel  "))
+          nvalid++;
+  }
+
+  FinalizeAbsTimes();
 
   *fLog << endl;  
@@ -779 +863 @@
       // (already done for calibration sigma in MHCalibrationCam::CalcAverageSigma()
       //
-      pix.SetPedRMS(pix.GetPedRms()*TMath::Sqrt((Float_t)arr[aidx]),
-                    pix.GetPedRmsErr()*TMath::Sqrt((Float_t)arr[aidx]));
+      const Double_t sqrtnum = TMath::Sqrt((Double_t)arr[aidx]);
+
+      pix.SetPedRMS(pix.GetPedRms()*sqrtnum, pix.GetPedRmsErr()*sqrtnum);
       pix.SetSigma (pix.GetSigma()/pix.GetFFactorFADC2Phe());
 
       FinalizeCharges(pix, fCam->GetAverageBadArea(aidx),"area id");
-      FinalizeArrivalTimes(pix, fCam->GetAverageBadArea(aidx), "area id");
     }
   
@@ -850 +934 @@
   PrintUncalibrated(MBadPixelsPix::kLoGainSaturation,    
                     "Low Gain Saturation:                              ");
-  PrintUncalibrated(MBadPixelsPix::kMeanTimeInFirstBin, 
-                    Form("%s%2.1f%s","Mean Abs. Arr. Time in First ",1.," Bin(s):          "));
-  PrintUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins, 
-                    Form("%s%2.1f%s","Mean Abs. Arr. Time in Last ",2.," Bin(s):           "));
+//  PrintUncalibrated(MBadPixelsPix::kMeanTimeInFirstBin,
+//                    Form("%s%2.1f%s","Mean Abs. Arr. Time in First ",1.," Bin(s):          "));
+//  PrintUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins,
+//                    Form("%s%2.1f%s","Mean Abs. Arr. Time in Last ",2.," Bin(s):           "));
   PrintUncalibrated(MBadPixelsPix::kHiGainOverFlow,    
                     "Pixels with High Gain Overflow:                   ");
@@ -974 +1058 @@
   const TString desc = Form("%7s%4d: ", what, cal.GetPixId());
 
+  if (cal.GetMean()<0)
+  {
+      *fLog << warn << desc << "Charge not fitted." << endl;
+      bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal);
+      return kFALSE;
+  }
+
+  if (cal.GetSigma()<0)
+  {
+      *fLog << warn << desc << "Charge Sigma invalid." << endl;
+      bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal);
+      return kFALSE;
+  }
+
   if (cal.GetMean() < fChargeLimit*cal.GetPedRms())
     {
@@ -1084 +1182 @@
 
   return kTRUE;
-}
-
-// -----------------------------------------------------------------------------------
-//
-// Test the arrival Times RMS of the pixel and set the bit
-// - MBadPixelsPix::kFluctuatingArrivalTimes
-//
-void MCalibrationChargeCalc::FinalizeArrivalTimes(MCalibrationChargePix &cal, MBadPixelsPix &bad, const char* what)
-{
-    if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
-        return;
-
-    const TString desc = Form("%7s%4d: ", what, cal.GetPixId());
-
-    if (cal.GetAbsTimeRms() > fArrTimeRmsLimit)
-    {
-        *fLog << warn;
-        *fLog << desc << "RMS of pulse arrival times: " << Form("%2.1f", cal.GetAbsTimeRms());
-        *fLog << " FADC sl. < " << Form("%2.1f", fArrTimeRmsLimit) << endl;
-        bad.SetUncalibrated( MBadPixelsPix::kFluctuatingArrivalTimes);
-    }
 }
 
@@ -1134 +1211 @@
             bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun   );
         }
-      
+/*
       if (IsCheckExtractionWindow())
         {
@@ -1143 +1220 @@
             bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );
         }
-
+ */
       if (IsCheckDeviatingBehavior())
         {
@@ -1275 +1352 @@
           continue;
         }
-      
+
+      // FIXME: WHAT IS THIS FOR? It is overwritten!
       lowlim[i] = areaphes[i] - fPheErrLowerLimit*TMath::Sqrt(areavars[i]);
       upplim[i] = areaphes[i] + fPheErrUpperLimit*TMath::Sqrt(areavars[i]);

trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.h

@@ -83 +83 @@
 
   // Pointers
-//  MBadPixelsIntensityCam         *fIntensBad;      //!  Bad Pixels
   MBadPixelsCam                  *fBadPixels;      //!  Bad Pixels
   MCalibrationChargeCam          *fCam;            //!  Calibrated Charges results of all pixels
@@ -128 +127 @@
 
   // functions
-  void   FinalizeArrivalTimes    ( MCalibrationChargePix &cal, MBadPixelsPix &bad, const char* what);
   void   FinalizeBadPixels       ();
   Bool_t FinalizeBlindCam        ();  
@@ -141 +139 @@
   Bool_t FinalizeUnsuitablePixels();
 
+  void FinalizeAbsTimes();
+
   const char* GetOutputFile();
 

trunk/MagicSoft/Mars/mcalib/MCalibrationRelTimeCalc.cc

@@ -17 +17 @@
 !
 !   Author(s): Markus Gaug  04/2004 <mailto:markus@ifae.es>
-!
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Author(s): Thomas Bretz <mailto:tbretz@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2007
 !
 !
@@ -242 +243 @@
 // MBadPixelsPix::kDeviatingTimeResolution
 //
+// see also MCalibrationChargeCalc::FinalizeAbsTimes
+//
 void MCalibrationRelTimeCalc::FinalizeRelTimes()
 {
@@ -260 +263 @@
 
             // check if pixel may not contain a valid value
-            if ((*fBadPixels)[i].IsUnsuitable())
+            if ((*fBadPixels)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
                 continue;
 
@@ -281 +284 @@
         Double_t med;
         const Double_t dev = MMath::MedianDev(n, arr.GetArray(), med);
+
+        // Calculate upper and lower limit
+        const Float_t lolim = TMath::Max(med-fRelTimeResolutionLimit*dev, 0.);
+        const Float_t hilim = TMath::Max(med+fRelTimeResolutionLimit*dev, 0.);
 
         // Now find the outliers
@@ -290 +297 @@
 
             // skip pixels already known to be unsuitable
-            if ((*fBadPixels)[i].IsUnsuitable())
+            if ((*fBadPixels)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
                 continue;
 
@@ -301 +308 @@
             if (res<0) //FIXME!!! How does this happen?
             {
-                *fLog << warn << "Pixel  " << setw(4) << i << ": Time resolution could not be calculated." << endl;
+                *fLog << warn << "Pixel  " << setw(4) << i << ": Rel-time resolution could not be calculated." << endl;
                 (*fBadPixels)[i].SetUncalibrated(MBadPixelsPix::kDeviatingTimeResolution);
                 continue;
@@ -307 +314 @@
 
             // Now compare to a lower and upper limit
-            const Float_t lolim = TMath::Max(med-fRelTimeResolutionLimit*dev, 0.);
-            const Float_t hilim = TMath::Max(med+fRelTimeResolutionLimit*dev, 0.);
-
             if (res<=lolim || res>=hilim)
             {
-                *fLog << warn << "Pixel  " << setw(4) << i << ": Deviating time resolution: "
-                    << Form("%4.2f", res) << " out of range "
+                *fLog << warn << "Pixel  " << setw(4) << i << ": Deviating rel-time resolution: "
+                    << Form("%5.2f", res) << " out of range "
                     << Form("[%4.2f,%4.2f]", lolim, hilim) << endl;