Index: /trunk/MagicSoft/Mars/mtemp/MSigmabar.cc
===================================================================
--- /trunk/MagicSoft/Mars/mtemp/MSigmabar.cc	(revision 1676)
+++ /trunk/MagicSoft/Mars/mtemp/MSigmabar.cc	(revision 1676)
@@ -0,0 +1,150 @@
+/* ======================================================================== *\
+!
+! *
+! * This file is part of MARS, the MAGIC Analysis and Reconstruction
+! * Software. It is distributed to you in the hope that it can be a useful
+! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
+! * It is distributed WITHOUT ANY WARRANTY.
+! *
+! * Permission to use, copy, modify and distribute this software and its
+! * documentation for any purpose is hereby granted without fee,
+! * provided that the above copyright notice appear in all copies and
+! * that both that copyright notice and this permission notice appear
+! * in supporting documentation. It is provided "as is" without express
+! * or implied warranty.
+! *
+!
+!
+!   Author(s): Robert Wagner  10/2002 <mailto:magicsoft@rwagner.de>
+!
+!   Copyright: MAGIC Software Development, 2002
+!
+!
+\* ======================================================================== */
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// MSigmabar                                                               //
+//                                                                         //
+// This is the storage container to hold information about the mean sigma  //
+// (aka Sigmabar) of all pedestals                                         //
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+#include "MSigmabar.h"
+
+#include <TMath.h>
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MGeomCam.h"
+#include "MPedestalCam.h"
+#include "MGeomPix.h"
+#include "MPedestalPix.h"
+
+ClassImp(MSigmabar);
+
+MSigmabar::MSigmabar(const char *name, const char *title)
+{
+    fName  = name  ? name  : "MSigmabar";
+    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.
+}
+
+
+// --------------------------------------------------------------------------
+//
+// 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
+//
+// Preliminary! Works for CT1 test, for real MAGIC crosschecks still have
+// to be done. Also implementation details will be updated, like 
+// 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 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 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
+  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/=fInnerPixels;
+  if (fSigmabarOuter != 0) fSigmabarOuter/=fOuterPixels; 
+  fSigmabar=sqrt(fSigmabarInner + fSigmabarOuter/( fOuterPixels==0 ? 1 : fRatioA)); 
+  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]);
+  }
+    
+  // Did all our calculations work? fOuterPixels==0 could happen, however.
+  return (fInnerPixels!=0);
+}
+
+
+void MSigmabar::Print(Option_t *) const
+{
+  *fLog << all;
+  *fLog << "Sigmabar     Overall " << fSigmabar;
+  *fLog << " Sectors: ";
+  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;
+  *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;
+}