Changeset 4568 for trunk/MagicSoft/Mars/mcalib

Timestamp:

08/10/04 17:53:21 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mcalib

Files:

• MCalibrationTestCalc.cc (modified) (14 diffs)
• MCalibrationTestCalc.h (modified) (4 diffs)

trunk/MagicSoft/Mars/mcalib/MCalibrationTestCalc.cc

@@ -85 +85 @@
 //
 MCalibrationTestCalc::MCalibrationTestCalc(const char *name, const char *title)
-    : fBadPixels(NULL), fTestCam(NULL), fGeom(NULL)
+    : fMaxNumBadPixelsCluster(-1), 
+      fBadPixels(NULL), fTestCam(NULL), fGeom(NULL)
 {
         
@@ -169 +170 @@
 {
 
-  if (GetNumExecutions()==0)
-    return kFALSE;
+  //  if (GetNumExecutions()==0)
+  //    return kFALSE;
 
   //
@@ -186 +187 @@
   *fLog << GetDescriptor() << ": Not interpolated pixels statistics:" << endl;  
 
-  PrintNotInterpolated();
-  
+  FinalizeNotInterpolated();
+  CalcMaxNumBadPixelsCluster();
+
+  *fLog << endl;
+  *fLog << " " << setw(7) << fMaxNumBadPixelsCluster 
+        << " pixels in biggest not-interpolateable cluster " << endl;
+  *fLog << endl;  
+
   SetLogStream(&gLog);
 
@@ -267 +274 @@
         }
 
-      areavars[aidx]    = (areavars[aidx] - areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx]) / (numareavalid[aidx]-1.);
-      areaphotons[aidx] = areaphotons[aidx] / numareavalid[aidx];
-
+      if (numareavalid[aidx] == 1)
+        areavars[aidx]  = 0.;
+      else if (numareavalid[aidx] == 0)
+        {
+          areaphotons[aidx] = -1.;
+          areavars[aidx]    = -1.;
+        }
+      else
+        {
+          areavars[aidx]    = (areavars[aidx] - areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx]) 
+                            / (numareavalid[aidx]-1.);
+          areaphotons[aidx] = areaphotons[aidx] / numareavalid[aidx];
+        }
+      
       if (areavars[aidx] < 0.)
         {
@@ -318 +336 @@
       upplim  [aidx] = mean  + fPhotErrLimit*sigma;
 
-      *fLog << " Fitted number of photons in area index " << aidx 
-            << ": "  << Form("%6.2f +- %6.2f",fittedmean[aidx],fittedsigma[aidx]) << endl;
+      *fLog << inf << GetDescriptor() 
+            << ": Fitted number of calib. equiv. Cher. photons in area index " << aidx 
+            << ": "  << Form("%7.2f +- %6.2f",fittedmean[aidx],fittedsigma[aidx]) << endl;
 
       delete hist;
     }
+
+  *fLog << endl;
 
   memset(numareavalid,0,nareas*sizeof(Int_t));
@@ -336 +357 @@
       
       MHCalibrationTestPix &pix = (MHCalibrationTestPix&)(*fTestCam)[i];
-      MBadPixelsPix        &bad = (*fBadPixels)[i];
 
       const Int_t    aidx   = (*fGeom)[i].GetAidx();
@@ -345 +365 @@
       if ( nphot < lowlim[aidx] || nphot > upplim[aidx] )
         {
-          *fLog << warn << "Deviating number of calibrated photons: " 
-                << Form("%4.2f",nphot) << " out of accepted limits: [" 
-                << Form("%4.2f%s%4.2f",lowlim[aidx],",",upplim[aidx]) << "] in pixel " << i << endl;
+          *fLog << warn << GetDescriptor() << ": Number of photons: " 
+                << Form("%8.2f out of %3.1f sigma limit: ",nphot,fPhotErrLimit)
+                << Form("[%8.2f,%8.2f] pixel%4i",lowlim[aidx],upplim[aidx],i) << endl;
+          MBadPixelsPix &bad = (*fBadPixels)[i];
           bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhots );
           bad.SetUnsuitable  ( MBadPixelsPix::kUnsuitableRun     );
@@ -362 +383 @@
     } 
 
+  *fLog << endl;
+
   for (UInt_t aidx=0; aidx<nareas; aidx++)
     {
 
-      areavars   [aidx] -=  areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx];
-      areaphotons[aidx] /=  numareavalid[aidx];
-      areavars   [aidx] /=  numareavalid[aidx]-1.;
-      
-      if (areavars[aidx] <= 0. || areaphotons[aidx] <= 0.)
-        {
-          *fLog << warn << " Mean number of photons per area in area index " 
-                << aidx << " cannot be calculated! Mean: " << areaphotons[aidx] 
+      if (numareavalid[aidx] == 1)
+        areavars[aidx] = 0.;
+      else if (numareavalid[aidx] == 0)
+        {
+          areaphotons[aidx] = -1.;
+          areavars[aidx]    = -1.;
+        }
+      else
+        {
+          areavars[aidx] = (areavars[aidx] - areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx]) 
+                         / (numareavalid[aidx]-1.);
+          areaphotons[aidx] /=  numareavalid[aidx];
+        }
+      
+      if (areavars[aidx] < 0. || areaphotons[aidx] <= 0.)
+        {
+          *fLog << warn << GetDescriptor() 
+                << ": Mean number of photons per area in area index " 
+                << aidx << " could not be calculated! Mean: " << areaphotons[aidx] 
                 << " Variance: " << areavars[aidx] << endl;
           continue;
         }
       
-      *fLog << " Mean number of photons per area in area index " << aidx 
-            << ": "  << Form("%5.4f +- %5.4f",areaphotons[aidx],TMath::Sqrt(areavars[aidx])) << endl;
-    }
+      *fLog << inf << GetDescriptor() << ": Mean number of equiv. Cher. photons "
+            << "per area in area idx " << aidx << ": "  
+            << Form("%5.3f+-%5.4f  [ph/mm^2]",areaphotons[aidx],TMath::Sqrt(areavars[aidx])) << endl;
+    }
+
+  *fLog << endl;
 
   for (UInt_t sector=0; sector<nsectors; sector++)
     {
-
-      sectorvars   [sector] -=  sectorphotons[sector]*sectorphotons[sector]/numsectorvalid[sector];
-      sectorphotons[sector] /=  numsectorvalid[sector];
-      sectorvars   [sector] /=  numsectorvalid[sector]-1.;
-
-      if (sectorvars[sector] <= 0. || sectorphotons[sector] <= 0.)
-        {
-          *fLog << warn << " Mean number of calibrated photons per area from sector " 
-                << sector << " cannot be calculated! " << endl;
-          continue;
-        }
-      
-      *fLog << " Mean number of photons per area in sector " << sector
-            << ": "  << Form("%5.4f +- %5.4f",sectorphotons[sector],TMath::Sqrt(sectorvars[sector])) << endl;
+      
+      if (numsectorvalid[sector] == 1)
+        sectorvars[sector] = 0.;
+      else if (numsectorvalid[sector] == 0)
+        {
+          sectorphotons[sector]  = -1.;
+          sectorvars[sector]     = -1.;
+        }
+      else
+        {
+          sectorvars[sector] = (sectorvars[sector] 
+                               - sectorphotons[sector]*sectorphotons[sector]/numsectorvalid[sector]
+                               ) 
+                             / (numsectorvalid[sector]-1.);
+          sectorphotons[sector] /=  numsectorvalid[sector];
+        }
+      
+      if (sectorvars[sector] < 0. || sectorphotons[sector] <= 0.)
+        {
+          *fLog << warn << GetDescriptor() 
+                << ": Mean number of calibrated photons per area in sector " 
+                << sector << " could not be calculated! Mean: " << sectorphotons[sector] 
+                << " Variance: " << sectorvars[sector] << endl;
+          continue;
+        }
+      
+  
+      *fLog << inf << GetDescriptor() << ": Mean number of equiv. Cher. photons "
+            << "per area in sector " << sector << ":   "  
+            << Form("%5.3f+-%5.4f  [ph/mm^2]",sectorphotons[sector],TMath::Sqrt(sectorvars[sector])) << endl;
     }
 
@@ -407 +460 @@
 // Print out statistics about not interpolated pixels
 // 
-void MCalibrationTestCalc::PrintNotInterpolated() const 
+void MCalibrationTestCalc::FinalizeNotInterpolated() 
 {
   
@@ -418 +471 @@
     newarr[aidx] = new TArrayI(0);
 
-  TArrayI numtot;
-  numtot.Set(areas);
+  fNumUninterpolateable.Set(areas);
 
   for (Int_t i=0; i<arr.GetSize(); i++)
@@ -428 +480 @@
       newarr[aidx]->Set(size+1);
       newarr[aidx]->AddAt(id,size);
-      numtot[aidx]++;
+      fNumUninterpolateable[aidx]++;
     }
 
@@ -436 +488 @@
     {
       *fLog << " " << setw(7) 
-            << numtot[aidx] << " not interpolateable pixels in area index " << aidx << endl;
+            << fNumUninterpolateable[aidx] 
+            << " not interpolateable pixels in area index " << aidx << endl;
       *fLog << " " << setw(7)
             << "Pixel software idx: ";
@@ -447 +500 @@
     }
   
-
   *fLog << " " << setw(7) << num << " total not interpolateable pixels " << endl;
   
 }
+
+void MCalibrationTestCalc::CalcMaxNumBadPixelsCluster()
+{
+
+  const TArrayI &arr = fTestCam->GetNotInterpolateablePixels();
+  const Int_t   size = arr.GetSize();
+  
+  if (size == 0)
+    {
+      fMaxNumBadPixelsCluster = 0; 
+      return;
+    }
+  
+  if (size == 1)
+    {
+      fMaxNumBadPixelsCluster = 1; 
+      return;
+    }
+
+  TArrayI knownpixels(0);
+  Int_t clustersize    = 1;
+  Int_t oldclustersize = 0;
+  //
+  // Loop over the not-interpolateable pixels:
+  //
+  for (Int_t i=0; i<size; i++)
+    {
+
+      const Int_t id   = arr[i];
+      const Int_t knownsize = knownpixels.GetSize();
+      knownpixels.Set(knownsize+1);
+      knownpixels[knownsize] = id;
+      LoopNeighbours(arr, knownpixels, clustersize, id);
+      if (clustersize > oldclustersize)
+        oldclustersize = clustersize;
+      clustersize = 1;
+    }
+
+  fMaxNumBadPixelsCluster = oldclustersize; 
+
+}
+
+
+void MCalibrationTestCalc::LoopNeighbours( const TArrayI &arr, TArrayI &knownpixels, Int_t &clustersize, const Int_t idx )
+{
+  
+  const MGeomPix &pix = (*fGeom)[idx];
+  const Byte_t neighbours = pix.GetNumNeighbors();
+
+  // 
+  // Loop over the next neighbours: 
+  // - Check if they are already in the list of known pixels
+  // - If not, call loopneighbours for the rest
+  // - grow clustersize for those
+  //
+  for (Int_t i=0;i<neighbours;i++)
+    {
+      const Int_t newid = pix.GetNeighbor(i);
+      Bool_t known = kFALSE;
+
+      for (Int_t j=knownpixels.GetSize()-1;j>=0;j--)
+        if (newid == knownpixels.At(j))
+          {
+            known = kTRUE;
+            break;
+          }
+      if (known)
+        continue;
+
+      for (Int_t k=0;k<arr.GetSize();k++)
+        if (newid == arr.At(k))
+          {
+            // This is an unknown, new pixel in the cluster!!
+            clustersize++;
+            const Int_t knownsize = knownpixels.GetSize();
+            knownpixels.Set(knownsize+1);
+            knownpixels[knownsize] = newid;
+            LoopNeighbours(arr, knownpixels, clustersize, newid);
+          }
+    }
+}
+
+
+
+
 
 // --------------------------------------------------------------------------
@@ -476 +613 @@
   return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
 }
+
+const Int_t MCalibrationTestCalc::GetNumUninterpolateable(const Int_t aidx) const
+{
+  if (aidx < 0)
+    return -1;
+
+  return aidx > fNumUninterpolateable.GetSize() ? -1 : fNumUninterpolateable[aidx]; 
+}

trunk/MagicSoft/Mars/mcalib/MCalibrationTestCalc.h

@@ -15 +15 @@
 #endif
 
+#ifndef ROOT_TArrayI
+#include <TArrayI.h>
+#endif
+
 class MHCalibrationTestCam;
 class MBadPixelsCam;
@@ -26 +30 @@
   // Variables
   Float_t fPhotErrLimit;               // Limit acceptance nr. cal. phots w.r.t. area idx mean (in sigmas)
+  Int_t   fMaxNumBadPixelsCluster;     // Number of not interpolateable pixels in biggest cluster
+
+  TArrayI fNumUninterpolateable;       // Number uninterpolated Pixels per area index 
 
   TString fOutputPath;                 // Path to the output file
@@ -38 +45 @@
   const char* GetOutputFile();
 
-  void   PrintNotInterpolated() const;
+  void   FinalizeNotInterpolated();
   void   FinalizeCalibratedPhotons() const;
-
+  void   CalcMaxNumBadPixelsCluster();
+  void   LoopNeighbours( const TArrayI &arr, TArrayI &known, Int_t &clustersize, const Int_t idx );
+  
   Int_t  PreProcess (MParList *pList);
   Bool_t ReInit     (MParList *pList); 
@@ -50 +59 @@
   MCalibrationTestCalc(const char *name=NULL, const char *title=NULL);
 
+  const Int_t GetMaxNumBadPixelsCluster () const  { return fMaxNumBadPixelsCluster; }
+  const Int_t GetNumUninterpolateable   ( const Int_t aidx ) const;
+
+  void SetOutputFile  ( TString file="TestCalibStat.txt" );
   void SetOutputPath  ( TString path="."                 );
-  void SetOutputFile  ( TString file="TestCalibStat.txt" );
-
-  void SetPhotErrLimit ( const Float_t f=fgPhotErrLimit   ) { fPhotErrLimit = f; }  
+  void SetPhotErrLimit( const Float_t f=fgPhotErrLimit   ) { fPhotErrLimit = f; }  
   
   ClassDef(MCalibrationTestCalc, 1)   // Task retrieving the results of MHCalibrationTestCam