Changeset 2821

Timestamp:

01/15/04 16:44:47 (21 years ago)

Author:

rico

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• macros/pedvsevent.C (added)
• manalysis/MPedCalcPedRun.cc (modified) (6 diffs)
• manalysis/MPedCalcPedRun.h (modified) (4 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -4 +4 @@
 
                                                  -*-*- END OF LINE -*-*-
+ 2004/01/15: Javier Rico
+        
+   * manalysis/MPedCalcPedRun.[h,cc]
+     - optimize the running time
+     - add (some) documentation
+     - correct treatment for the case of several input files 
+
+   * macros/pedvsevent.C
+     - added
+     - draw pedestal mean and rms vs event# for input pixel# and run 
+       file, and compares them to the global pedestal mean and rms 
+
 
  2004/01/15: Raquel de los Reyes

trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc

@@ -18 +18 @@
 !   Author(s): Josep Flix 04/2001 <mailto:jflix@ifae.es>
 !   Author(s): Thomas Bretz 05/2001 <mailto:tbretz@uni-sw.gwdg.de>
+!   Author(s): Sebastian Commichau 12/2003 
+!   Author(s): Javier Rico 01/2004 
 !
-!   Copyright: MAGIC Software Development, 2000-2001
+!   Copyright: MAGIC Software Development, 2000-2004
 !
 !
@@ -28 +30 @@
 //   MPedCalcPedRun                                                        //
 //                                                                         //
+//  This task takes a pedestal run file and fills MPedestalCam during      //
+//  the Process() with the pedestal and rms computed in an event basis.    //
+//  In the PostProcess() MPedestalCam is finally filled with the pedestal  //
+//  mean and rms computed in a run basis.                                  //
+//  More than one run (file) can be merged                                 //
+//                                                                         //
 //  Input Containers:                                                      //
 //   MRawEvtData                                                           //
@@ -33 +41 @@
 //  Output Containers:                                                     //
 //   MPedestalCam                                                          //
+//                                                                         //
 //                                                                         //
 /////////////////////////////////////////////////////////////////////////////
@@ -70 +79 @@
     if (!fRawEvt)
     {
-        *fLog << dbginf << "MRawEvtData not found... aborting." << endl;
+        *fLog << err << "MRawEvtData not found... aborting." << endl;
         return kFALSE;
     }
@@ -78 +87 @@
         return kFALSE;
 
-    MGeomCamMagic magiccam;
-
-    fSumx.Set(magiccam.GetNumPixels());
-    fSumx2.Set(magiccam.GetNumPixels());
-    
-    for(UInt_t i=0;i<magiccam.GetNumPixels();i++)
-    {
-        fSumx.AddAt(0,i);
-        fSumx2.AddAt(0,i);
-    }
+    fNumSamplesTot=0;
 
     return kTRUE;
 }
 
-Bool_t MPedCalcPedRun::ReInit(MParList *pList )   
+Bool_t MPedCalcPedRun::ReInit(MParList *pList)
 {
+    const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
+    if (!runheader)
+    {
+        *fLog << warn << dbginf;
+        *fLog << "Warning - cannot check file type, MRawRunHeader not found." << endl;
+    }
+    else
+        if (runheader->GetRunType() == kRTMonteCarlo)
+            return kTRUE;
 
-    fRunheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
-    if (!fRunheader)
-        {
-        *fLog << warn << dbginf << 
-                "Warning - cannot check file type, MRawRunHeader not found." << endl;
-        }
-    else
-        if (fRunheader->GetRunType() == kRTMonteCarlo)
-        {
-            return kTRUE;
-        }
+    fNumPixels = fPedestals->GetSize();
 
-    fNumHiGainSamples =  fRunheader->GetNumSamplesHiGain();
+    if(fSumx.GetSize()==0)
+      {
+        fSumx.Set(fNumPixels);
+        fSumx2.Set(fNumPixels);
 
-    fPedestals->InitSize(fRunheader->GetNumPixel());
+        memset(fSumx.GetArray(),  0, sizeof(Float_t)*fNumPixels);
+        memset(fSumx2.GetArray(), 0, sizeof(Float_t)*fNumPixels);
+      }
+
+    // Calculate an even number for the hi gain samples to avoid
+    // biases due to the fluctuation in pedestal from one slice to
+    // the other one
+    fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1;
 
     return kTRUE;
-
 }
-
 
 Int_t MPedCalcPedRun::Process()
 {
-
     MRawEvtPixelIter pixel(fRawEvt);
 
     while (pixel.Next())
     {
-        Byte_t shift=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? 0:1;
               Byte_t *ptr = pixel.GetHiGainSamples();
-        const Byte_t *end = ptr + fRawEvt->GetNumHiGainSamples()-shift;
+        const Byte_t *end = ptr + fNumHiGainSamples;
 
-        const Float_t higainped = CalcHiGainMean(ptr, end);
-        const Float_t higainrms = CalcHiGainRms(ptr, end, higainped);
-
-        const UInt_t pixid = pixel.GetPixelId();
-        MPedestalPix &pix = (*fPedestals)[pixid];
+        UInt_t sum = 0;
+        UInt_t sqr = 0;
         
-        // cumulate the sum of pedestals and of pedestal squares
-        fSumx.AddAt(higainped+fSumx.At(pixid),pixid);
-        fSumx2.AddAt(GetSumx2(ptr, end)+fSumx2.At(pixid),pixid);
-
-        // set the value of the pedestal and rms computed from the processed event
-        pix.Set(higainped, higainrms);
+        do
+          {
+            sum += *ptr;
+            sqr += *ptr * *ptr;
+          }
+        while (++ptr != end);
+        
+        const Float_t msum = (Float_t)sum;
+        const Float_t msqr = (Float_t)sqr;
+        
+        const Float_t higainped = msum/fNumHiGainSamples;
+        const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.));
+        
+        const UInt_t idx = pixel.GetPixelId();
+        (*fPedestals)[idx].Set(higainped, higainrms);
+        
+        fSumx[idx]  += msum;
+        fSumx2[idx] += sqr;
     }
-
+    
     fPedestals->SetReadyToSave();
-
+    fNumSamplesTot+=fNumHiGainSamples;
+    
     return kTRUE;
 }
@@ -148 +163 @@
 Int_t MPedCalcPedRun::PostProcess()
 {
-// Compute pedestals and rms from the whole run
-
-    MRawEvtPixelIter pixel(fRawEvt);
-
-    while (pixel.Next())
+  // Compute pedestals and rms from the whole run
+  const Int_t n  = fNumSamplesTot;
+  
+  MRawEvtPixelIter pixel(fRawEvt);
+  
+  while (pixel.Next())
     {
-        const UInt_t pixid = pixel.GetPixelId();
-        MPedestalPix &pix = (*fPedestals)[pixid];
-
-        const Int_t N = GetNumExecutions();
-        const Float_t sum = fSumx.At(pixid);
-        const Float_t sum2 = fSumx2.At(pixid);
-        const Float_t higainped = sum/N;
-        const Float_t higainrms = sqrt(1./(N-1.)*(sum2-sum*sum/N));
-        pix.Set(higainped,higainrms);
-
+      const UInt_t pixid = pixel.GetPixelId();
+      
+      const Float_t sum  = fSumx.At(pixid);
+      const Float_t sum2 = fSumx2.At(pixid);
+      
+      const Float_t higainped = sum/n;
+      const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
+      
+      (*fPedestals)[pixid].Set(higainped, higainrms);
     }
-    return kTRUE;
-
+  
+  return kTRUE;
 }
 
 
-Float_t MPedCalcPedRun::CalcHiGainMean(Byte_t *ptr, const Byte_t *end) const
-{
-    Int_t sum=0;
-    Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) 
-                        ? fNumHiGainSamples : fNumHiGainSamples-1;
-
-    do sum += *ptr;
-    while (++ptr != end);
-   
-    return (Float_t)sum/EvenNumSamples;
-}
-
-
-
-Float_t MPedCalcPedRun::GetSumx2(Byte_t *ptr, const Byte_t *end) const
-{
-    Float_t square = 0;
-
-    // Take an even number of time slices to avoid biases due to A/B effect
-    Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? fNumHiGainSamples:fNumHiGainSamples-1;
-
-    do
-    {
-        const Float_t val = (Float_t)(*ptr);
-
-        square += val*val;
-    } while (++ptr != end);
-
-    return square/EvenNumSamples;
-}
-
-
-Float_t MPedCalcPedRun::CalcHiGainRms(Byte_t *ptr, const Byte_t *end, Float_t higainped) const
-{
-    Float_t rms = 0;
-    Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? fNumHiGainSamples:fNumHiGainSamples-1;
-
-    do
-    {
-        const Float_t diff = (Float_t)(*ptr)-higainped;
-
-        rms += diff*diff;
-    } while (++ptr != end);
-
-    return sqrt(rms/(EvenNumSamples-1));
-}
-
-
-
-/*
-Float_t MPedCalcPedRun::CalcHiGainMeanErr(Float_t higainrms) const
-{
-    return higainrms/sqrt((float)fNumHiGainSamples);
-}
-
-Float_t MPedCalcPedRun::CalcHiGainRmsErr(Float_t higainrms) const
-{
-    return higainrms/sqrt(2.*fNumHiGainSamples);
-}
-*/

trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h

@@ -16 +16 @@
 #include <TArrayF.h>
 
-class MRawRunHeader;
 class MRawEvtData;
 class MPedestalCam;
@@ -22 +21 @@
 class MPedCalcPedRun : public MTask
 {
-    Byte_t fNumHiGainSamples;
+    Byte_t   fNumHiGainSamples;
+    UShort_t fNumPixels;
+    ULong_t   fNumSamplesTot;
 
-    MRawRunHeader *fRunheader; // raw event run header
     MRawEvtData  *fRawEvt;     // raw event data (time slices)
     MPedestalCam *fPedestals;  // Pedestals of all pixels in the camera
@@ -30 +30 @@
     TArrayF fSumx;   // sum of values
     TArrayF fSumx2;  // sum of squared values
-
-    Float_t CalcHiGainMean(Byte_t *ptr, const Byte_t *end) const;
-    Float_t CalcHiGainRms(Byte_t *ptr, const Byte_t *end, Float_t higainped) const;
-    Float_t GetSumx2(Byte_t* ptr, const Byte_t* end) const; 
-   //Float_t CalcHiGainMeanErr(Float_t higainrms) const;
-    //Float_t CalcHiGainRmsErr(Float_t higainrms) const;
 
     Bool_t ReInit(MParList *pList);
@@ -44 +38 @@
 
 public:
-
     MPedCalcPedRun(const char *name=NULL, const char *title=NULL);