Changeset 1748 for trunk/MagicSoft/Mars

Timestamp:

02/07/03 13:26:50 (22 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• manalysis/MPadding.cc (added)
• manalysis/MPadding.h (added)
• manalysis/MSigmabar.cc (modified) (6 diffs)
• manalysis/MSigmabar.h (modified) (4 diffs)
• manalysis/MSigmabarCalc.cc (modified) (2 diffs)
• manalysis/MSigmabarCalc.h (modified) (2 diffs)
• mhist/MHSigmaTheta.cc (added)
• mhist/MHSigmaTheta.h (added)

trunk/MagicSoft/Mars/manalysis/MSigmabar.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Robert Wagner  10/2002 <mailto:magicsoft@rwagner.de>
+!   Author(s): Robert Wagner   10/2002 <mailto:magicsoft@rwagner.de>
+!              Wolfgang Wittek 01/2003 <mailto:wittek@mppmu.mpg.de>
 !
 !   Copyright: MAGIC Software Development, 2002
@@ -27 +28 @@
 // MSigmabar                                                               //
 //                                                                         //
-// This is the storage container to hold information about the mean sigma  //
-// (aka Sigmabar) of all pedestals                                         //
+// This is the storage container to hold information about                 //
+// "average" pedestal sigmas                                               //
+//                                                                         //
+// In calculating averages all sigmas are normalized to the area of pixel 0//
 //                                                                         //
 /////////////////////////////////////////////////////////////////////////////
@@ -37 +40 @@
 #include "MLog.h"
 #include "MLogManip.h"
-
+#include "MParList.h"
 #include "MGeomCam.h"
 #include "MPedestalCam.h"
+#include "MPedestalPix.h"
 #include "MGeomPix.h"
-#include "MPedestalPix.h"
+#include "MCerPhotEvt.h"
+#include "MCerPhotPix.h"
 
 ClassImp(MSigmabar);
 
+// --------------------------------------------------------------------------
+//
 MSigmabar::MSigmabar(const char *name, const char *title)
 {
@@ -50 +57 @@
     fTitle = title ? title : "Storage container for Sigmabar";
     
-    fSigmabar = 0.0;
-    fSigmabarInner = 0.0;
-    fSigmabarOuter = 0.0;
-    fRatioA = 0.0;
 
     fCalcPixNum=kTRUE;
 }
 
+// --------------------------------------------------------------------------
+//
 MSigmabar::~MSigmabar()
 {
   // do nothing special.
 }
-
 
 // --------------------------------------------------------------------------
@@ -68 +72 @@
 // Actual calculation of sigmabar. This is done for each of the six sectors
 // separately due to their possibly different HV behavior. Also inner and
-// outer pixels are treated separately
+// outer pixels are treated separately.
 //
 // Preliminary! Works for CT1 test, for real MAGIC crosschecks still have
@@ -74 +78 @@
 // determination of sector to which a respective pixel belongs
 //
-Bool_t MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped)
-{
-  const UInt_t npix = ped.GetSize();
+Float_t MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped, 
+                        const MCerPhotEvt &evt)
+{
+  fSigmabar      = 0.0;
+  fSigmabarInner = 0.0;
+  fSigmabarOuter = 0.0;
+
+
   Float_t innerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
   Float_t outerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
-  Int_t innerPixels[6] = {0,0,0,0,0,0};
-  Int_t outerPixels[6] = {0,0,0,0,0,0};
+  Int_t innerPixels[6]       = {0,0,0,0,0,0};
+  Int_t outerPixels[6]       = {0,0,0,0,0,0};
+
   Int_t currentSector;
   Float_t angle;
   
-  for (UInt_t i=0; i<npix;i++)
-    {
-      const MGeomPix gpix = geom[i];
-      angle=((6*atan2(gpix.GetX(),gpix.GetY())/(2*TMath::Pi()))-0.5);
-      if (angle<0) angle+=6;
-      currentSector=(Int_t)angle;          
-      geom.GetPixRatio(i) == 1 ? innerPixels[currentSector]++ : outerPixels[currentSector]++;   
-      
-      const MPedestalPix pix = ped[i]; 
-      geom.GetPixRatio(i) == 1 ? innerSquaredSum[currentSector]+=(pix.GetSigma()*pix.GetSigma()) : outerSquaredSum[currentSector]+=((pix.GetSigma()*pix.GetSigma()) / geom.GetPixRatio(i));
-
-      // Get once and forever the ratio of areas outer/inner pixels
-      if (fRatioA && (geom.GetPixRatio(i) != 1)) fRatioA=geom.GetPixRatio(i);
-    }
-
-  // Overall Sigma
+  // sum up sigma**2 for each sector, separately for inner and outer region;
+  // all pixels are renormalized to the area of pixel 0
+  //
+  // consider all pixels with Cherenkov photon information
+  // and require "Used"
+  //
+
+  const UInt_t npix = evt.GetNumPixels();
+
+  for (UInt_t i=0; i<npix; i++)
+  {
+      MCerPhotPix &cerpix = evt.operator[](i);
+      if (!cerpix.IsPixelUsed())
+        continue;
+  
+      if ( cerpix.GetNumPhotons() == 0 )
+      {
+        *fLog << "MSigmabar::Calc(); no.of photons is 0 for used pixel" 
+              << endl;
+        continue;
+      }
+
+      Int_t j = cerpix.GetPixId();
+      Double_t Area = geom.GetPixRatio(j);
+
+      const MGeomPix    &gpix = geom[j];
+      const MPedestalPix &pix =  ped[j]; 
+
+      //angle = 6.0*atan2(gpix.GetX(),gpix.GetY()) / (2.0*TMath::Pi()) - 0.5;
+      //if (angle<0.0) angle+=6.0;
+
+      angle = 6.0*atan2(gpix.GetY(),gpix.GetX()) / (2.0*TMath::Pi());
+      if (angle<0.0) angle+=6.0;
+      currentSector=(Int_t)angle;
+       
+      // count only those pixels which have a sigma != 0.0 
+      Float_t sigma = pix.GetMeanRms();
+
+      if ( sigma != 0.0 )
+      {  
+        if (Area < 1.5)
+        {
+          innerPixels[currentSector]++;
+          innerSquaredSum[currentSector]+= sigma*sigma / Area;
+        }
+        else
+        {
+          outerPixels[currentSector]++;
+          outerSquaredSum[currentSector]+= sigma*sigma / Area;
+        }
+      }
+  }
+ 
   fSigmabarInner=0; fInnerPixels=0;
   fSigmabarOuter=0; fOuterPixels=0;
   for (UInt_t i=0; i<6; i++) {
-    fSigmabarInner+=innerSquaredSum[i];
-    fInnerPixels  +=innerPixels[i];
-    fSigmabarOuter+=outerSquaredSum[i];
-    fOuterPixels  +=outerPixels[i];
+    fSigmabarInner += innerSquaredSum[i];
+    fInnerPixels   += innerPixels[i];
+    fSigmabarOuter += outerSquaredSum[i];
+    fOuterPixels   += outerPixels[i];
   }
 
-  fSigmabarInner/=fInnerPixels;
-  if (fSigmabarOuter != 0) fSigmabarOuter/=fOuterPixels; 
-  fSigmabar=sqrt(fSigmabarInner + fSigmabarOuter/( fOuterPixels==0 ? 1 : fRatioA)); 
-  fSigmabarInner=sqrt(fSigmabarInner);
-  fSigmabarOuter=sqrt(fSigmabarOuter);
+
+  // this is the sqrt of the average sigma**2;
+  //
+  if ( (fInnerPixels+fOuterPixels) > 0)
+    fSigmabar=sqrt(   (fSigmabarInner + fSigmabarOuter)
+                    / (fInnerPixels   + fOuterPixels) ); 
+
+  if (fInnerPixels > 0) fSigmabarInner /= fInnerPixels;
+  if (fOuterPixels > 0) fSigmabarOuter /= fOuterPixels;
+ 
+  // this is the sqrt of the average sigma**2
+  // for the inner and outer pixels respectively
+  //
+  fSigmabarInner = sqrt( fSigmabarInner );
+  fSigmabarOuter = sqrt( fSigmabarOuter );
   
+
   for (UInt_t i=0; i<6; i++) {
-    fSigmabarInnerSector[i]=innerSquaredSum[i]/innerPixels[i];
-    fSigmabarOuterSector[i]=outerSquaredSum[i]/( outerSquaredSum[i]==0 ? 1: outerPixels[i] );
-
-    fSigmabarSector[i]=sqrt(fSigmabarInnerSector[i] + fSigmabarOuterSector[i]*fRatioA);
-    fSigmabarInnerSector[i]=sqrt(fSigmabarInnerSector[i]);
-    fSigmabarOuterSector[i]=sqrt(fSigmabarOuterSector[i]);
+    fSigmabarSector[i]         = 0.0;
+    fSigmabarInnerSector[i]    = 0.0;
+    fSigmabarOuterSector[i]    = 0.0;
+
+    if ( (innerPixels[i]+outerPixels[i]) > 0)
+      fSigmabarSector[i] = sqrt(  (innerSquaredSum[i] + outerSquaredSum[i])
+                                / (innerPixels[i]     + outerPixels[i]    ) );
+    if ( innerPixels[i] > 0)
+      fSigmabarInnerSector[i] = innerSquaredSum[i] / innerPixels[i];
+    if ( outerPixels[i] > 0)
+      fSigmabarOuterSector[i] = outerSquaredSum[i] / outerPixels[i];
+
+
+    fSigmabarInnerSector[i] = sqrt( fSigmabarInnerSector[i] );
+    fSigmabarOuterSector[i] = sqrt( fSigmabarOuterSector[i] );
   }
     
-  // Did all our calculations work? fOuterPixels==0 could happen, however.
-  return (fInnerPixels!=0);
-}
-
-
+  return (fSigmabar);
+}
+
+// --------------------------------------------------------------------------
+//
 void MSigmabar::Print(Option_t *) const
 {
-  *fLog << all;
-  *fLog << "Sigmabar     Overall " << fSigmabar;
+  *fLog << endl;
+  *fLog << "Total number of inner pixels is " << fInnerPixels << endl;
+  *fLog << "Total number of outer pixels is " << fOuterPixels << endl;
+  *fLog << endl;
+
+  *fLog << "Sigmabar     Overall : " << fSigmabar << "   ";
   *fLog << " Sectors: ";
-  for (Int_t i=0;i<6;i++) *fLog << fSigmabarSector[i] << " ";
-  *fLog << endl;
-  *fLog << "Sigmabar     Inner   " << fSigmabarInner;
+  for (Int_t i=0;i<6;i++) *fLog << fSigmabarSector[i] << ", ";
+  *fLog << endl;
+
+
+  *fLog << "Sigmabar     Inner   : " << fSigmabarInner << "   ";
   *fLog << " Sectors: ";
-  for (Int_t i=0;i<6;i++) *fLog << fSigmabarInnerSector[i] << " ";
-  *fLog << endl;
-  *fLog << "Sigmabar     Outer   " << fSigmabarOuter;
+  for (Int_t i=0;i<6;i++) *fLog << fSigmabarInnerSector[i] << ", ";
+  *fLog << endl;
+
+
+  *fLog << "Sigmabar     Outer   : " << fSigmabarOuter << "   ";
   *fLog << " Sectors: ";
-  for (Int_t i=0;i<6;i++) *fLog << fSigmabarOuterSector[i] << " ";
-  *fLog << endl;
-  *fLog << "Total number of inner pixels found is " << fInnerPixels << endl;
-  *fLog << "Total number of outer pixels found is " << fOuterPixels << endl;
-  *fLog << "Ratio of areas outer/inner pixels found is " << fRatioA << endl;
-  *fLog << endl;
-}
+  for (Int_t i=0;i<6;i++) *fLog << fSigmabarOuterSector[i] << ", ";
+  *fLog << endl;
+
+}
+// --------------------------------------------------------------------------
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MSigmabar.h

@@ -4 +4 @@
 #ifndef MARS_MParContainer
 #include "MParContainer.h"
+#endif
+
+#ifndef MARS_MParList
+#include "MParList.h"
 #endif
 
@@ -14 +18 @@
 #endif
 
+class MCerPhotEvt;
+
 class MSigmabar : public MParContainer
 {
 private:
-    Float_t fSigmabar; // Sigmabar (mean standard deviation) of pedestal
+    Float_t fSigmabar;          // Sigmabar ("average" RMS) of pedestal
+    Float_t fSigmabarInner;     // --only for inner pixels
+    Float_t fSigmabarOuter;     // --only for outer pixels  
+
     Float_t fSigmabarSector[6]; // --for the 6 sectors of the camera
     Float_t fSigmabarInnerSector[6];
     Float_t fSigmabarOuterSector[6];
-    Float_t fSigmabarInner; // --only for inner pixels
-    Float_t fSigmabarOuter; // --only for outer pixels  
-    UInt_t  fInnerPixels; // Overall number of inner pixels
-    UInt_t  fOuterPixels; // Overall number of outer pixels
-    Float_t fRatioA; // Ratio of areas (outer/inner pixels)
-    Bool_t fCalcPixNum;
 
+    UInt_t  fInnerPixels;       // Overall number of inner pixels
+    UInt_t  fOuterPixels;       // Overall number of outer pixels
+
+    Bool_t  fCalcPixNum;
+    
 public:
 
@@ -46 +54 @@
     //    void SetSigmabarOuter(Float_t f) { fSigmabarOuter = f; }   
 
-    Bool_t Calc(const MGeomCam &geom, const MPedestalCam &ped);
+    Float_t Calc(const MGeomCam &geom, const MPedestalCam &ped, const MCerPhotEvt &evt);
       
     ClassDef(MSigmabar, 1)  // Storage Container for Sigmabar
@@ -53 +61 @@
 #endif
 
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc

@@ -110 +110 @@
     if (!fMcEvt)
       {
-        *fLog << err << dbginf << "MHSigmabarTheta : MMcEvt not found... aborting." << endl;
+        *fLog << err << dbginf << "MHSigmabaCalc : MMcEvt not found... aborting." << endl;
+        return kFALSE;
+      }
+
+    fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
+    if (!fEvt)
+      {
+        *fLog << err << dbginf << "MHSigmabarCalc : MCerPhotEvt not found... aborting." << endl;
         return kFALSE;
       }
@@ -126 +133 @@
 Bool_t MSigmabarCalc::Process()
 {
-  Bool_t rc = fSig->Calc(*fCam, *fPed);    
-  fSigmabarMax = TMath::Max((Double_t)fSig->GetSigmabar(), fSigmabarMax);
-  fSigmabarMin = TMath::Min((Double_t)fSig->GetSigmabar(), fSigmabarMin);
+  Float_t rc = fSig->Calc(*fCam, *fPed, *fEvt);    
+  fSigmabarMax = TMath::Max((Double_t)rc, fSigmabarMax);
+  fSigmabarMin = TMath::Min((Double_t)rc, fSigmabarMin);
 
-  if (fMcEvt->GetTelescopeTheta()*kRad2Deg < 90)
+  if (fMcEvt->GetTelescopeTheta()*kRad2Deg < 120)
     fThetaMax    = TMath::Max(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMax);
   fThetaMin    = TMath::Min(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMin);
 
-  return rc;
+  return kTRUE;
 }
 

trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.h

@@ -43 +43 @@
     MSigmabarParam *fSigParam;
     MMcEvt *fMcEvt;
+    MCerPhotEvt *fEvt;
     void Reset();
 
@@ -58 +59 @@
 #endif
 
+
+