Changeset 5155 for trunk/MagicSoft/Mars/msignal

Timestamp:

10/01/04 11:55:29 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/msignal

Files:

• MExtractTimeAndChargeDigitalFilter.cc (modified) (17 diffs)
• MExtractTimeAndChargeDigitalFilter.h (modified) (3 diffs)

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc

@@ -60 +60 @@
 #include "TString.h"
 
+#include <fstream>
+
 ClassImp(MExtractTimeAndChargeDigitalFilter);
 
 using namespace std;
 
-const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst  =  0;
-const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast   = 14;
-const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst  =  3;
-const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast   = 14;
-const Int_t  MExtractTimeAndChargeDigitalFilter::fgHiGainWindowSize  = 6;
-const Int_t  MExtractTimeAndChargeDigitalFilter::fgLoGainWindowSize  = 6;
-const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolution = 10;
+const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst       =  0;
+const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast        = 14;
+const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst       =  3;
+const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast        = 14;
+const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain  = 6;
+const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain  = 6;
+const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
+const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
 // --------------------------------------------------------------------------
 //
@@ -92 +95 @@
   SetBinningResolution();
 
-  SetWeightsFile("");
+  ReadWeightsFile("");
 }
 
@@ -108 +111 @@
 {
 
+  if (windowh != fgWindowSizeHiGain)
+    *fLog << warn << GetDescriptor() 
+          << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default." << endl;
+  if (windowl != fgWindowSizeLoGain)
+    *fLog << warn << GetDescriptor() 
+          << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default" << endl;
+
   fWindowSizeHiGain = windowh;
   fWindowSizeLoGain = windowl;
@@ -162 +172 @@
 // If filenname is empty, then all weights will be set to 1.
 //
-void MExtractTimeAndChargeDigitalFilter::SetWeightsFile(TString filename)
-{
-
-  fAmpWeightsHiGain .Set(fBinningResolution*fWindowSizeHiGain); 
-  fAmpWeightsLoGain .Set(fBinningResolution*fWindowSizeLoGain); 
-  fTimeWeightsHiGain.Set(fBinningResolution*fWindowSizeHiGain); 
-  fTimeWeightsLoGain.Set(fBinningResolution*fWindowSizeLoGain); 
+Bool_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename)
+{
+
+  fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain); 
+  fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain); 
+  fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain); 
+  fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain); 
   
   if (filename.IsNull())
@@ -182 +192 @@
           fTimeWeightsLoGain[i] = 1.;
         }
-      return;
-    }
-
-  TFile *f = new TFile(filename);
-  TH1F  *hw_amp  = (TH1F*)f->Get("hw_amp");
-  TH1F  *hw_time = (TH1F*)f->Get("hw_time");
-
-  if (!hw_amp)
-    {
-      *fLog << warn << GetDescriptor()
-            << ": Could not find hw_amp in " << filename << endl;
-      return;
-    }
-  
-  if (!hw_time)
-    {
-      *fLog << warn << GetDescriptor()
-            << ": Could not find hw_time in " << filename << endl;
-      return;
-    }
-  
-  for (UInt_t i=0; i<fAmpWeightsHiGain.GetSize(); i++)
-    {
-      fAmpWeightsHiGain [i] = hw_amp->GetBinContent(i+1);
-      fTimeWeightsHiGain[i] = hw_time->GetBinContent(i+1);
-    }
-  for (UInt_t i=0; i<fAmpWeightsLoGain.GetSize(); i++)
-    {
-      fAmpWeightsLoGain [i] = hw_amp->GetBinContent(i+1);
-      fTimeWeightsLoGain[i] = hw_time->GetBinContent(i+1);
-    }
-  //    cout << "w_amp[" << i << "] = " << fw_amp[i] << endl;
-  delete f;
+      return kTRUE;
+    }
+
+  ifstream fin(filename.Data());
+  
+  if (!fin)
+    {
+      *fLog << err << GetDescriptor()
+            << ": No weights file found: " << filename << endl;
+      return kFALSE;
+    }
+
+  Int_t len = 0;
+  Int_t cnt = 0;
+  Bool_t hi = kFALSE;
+  Bool_t lo = kFALSE;
+  
+  TString str;
+
+  while (1)
+    {
+
+      str.ReadLine(fin);
+      if (!fin)
+        break;
+      
+
+      if (str.Contains("# High Gain Weights:"))
+        {
+          str.ReplaceAll("# High Gain Weights:","");
+          sscanf(str.Data(),"%2i%2i",&fWindowSizeHiGain,&fBinningResolutionHiGain);
+          *fLog << inf << "Found number of High Gain slices: " << fWindowSizeHiGain 
+                       << " and High Gain resolution: " << fBinningResolutionHiGain << endl;
+          len = fBinningResolutionHiGain*fWindowSizeHiGain;
+          fAmpWeightsHiGain .Set(len); 
+          fTimeWeightsHiGain.Set(len); 
+          hi = kTRUE;
+          continue;
+        }
+      
+      if (str.Contains("# Low Gain Weights:"))
+        {
+          str.ReplaceAll("# Low Gain Weights:","");
+          sscanf(str.Data(),"%2i%2i",&fWindowSizeLoGain,&fBinningResolutionLoGain);
+          *fLog << inf << "Found number of Low Gain slices: " << fWindowSizeLoGain 
+                       << " and Low Gain resolution: " << fBinningResolutionLoGain << endl;
+          len = fBinningResolutionLoGain*fWindowSizeHiGain;
+          fAmpWeightsLoGain .Set(len); 
+          fTimeWeightsLoGain.Set(len); 
+          lo = kTRUE;
+          continue;
+        }
+      
+      if (str.Contains("#"))
+        continue;
+
+      if (len == 0)
+        continue;
+
+      sscanf(str.Data(),"\t%f\t%f",lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt],
+                                         lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt]);
+
+      if (++cnt == len)
+        {
+          len = 0;
+          cnt = 0;
+        }
+    }
+
+  if (cnt != len)
+    {
+      *fLog << err << GetDescriptor()
+            << ": Size mismatch in weights file " << filename << endl;
+      return kFALSE;
+    }
+
+  if (!hi)
+    {
+      *fLog << err << GetDescriptor()
+            << ": No correct header found in weights file " << filename << endl;
+      return kFALSE;
+    }
+
+  return kTRUE;
+}
+
+//----------------------------------------------------------------------------
+//
+// Read a pre-defined weights file into the class. 
+// This is mandatory for the extraction
+//
+Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename)
+{
+
+
+  Float_t amp[60];
+  Float_t tim[60];  
+  
+  ofstream fn(filename.Data());
+
+  fn << "# High Gain Weights: 6 10" << endl;
+  fn << "# (Amplitude)  (Time) " << endl;
+
+  for (UInt_t i=0; i<60; i++)
+    {
+      fn << "\t" << amp[i] << "\t" << tim[i] << endl;
+    }
+
+  fn << "# Low Gain Weights: 6 10" << endl;
+  fn << "# (Amplitude)  (Time) " << endl;
+
+  for (UInt_t i=0; i<60; i++)
+    {
+      fn << "\t" << amp[i] << "\t" << tim[i] << endl;
+    }
+  return kTRUE;
 }
 
@@ -238 +331 @@
 
   fLoGainSignal.Set(range);
+
+  fTimeShiftHiGain = (Float_t)fHiGainFirst + 0.5 + 1./fBinningResolutionHiGain;
+  fTimeShiftLoGain = (Float_t)fLoGainFirst + 0.5 + 1./fBinningResolutionLoGain;
 
   return kTRUE;
@@ -300 +396 @@
     return;
 
-#if 0
-  if (maxpos < 2)
-    {
-      time = (Float_t)fHiGainFirst;
-      return;
-    }
-#endif
-  
   Float_t time_sum  = 0.;
   Float_t fmax      = 0.;
@@ -334 +422 @@
       for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
       {
-        const Int_t   idx = fBinningResolution*sample+fBinningResolutionHalf;
+        const Int_t   idx = fBinningResolutionHiGain*sample+fBinningResolutionHalfHiGain;
         const Float_t pex = fHiGainSignal[sample+i]-PedMean[(sample+i+abflag) & 0x1];
         sum              += fAmpWeightsHiGain [idx]*pex; 
@@ -351 +439 @@
     {
       ftime_max        /= fmax;
-      Int_t t_iter      = Int_t(ftime_max*fBinningResolution);
+      Int_t t_iter      = Int_t(ftime_max*fBinningResolutionHiGain);
       Int_t sample_iter = 0;
 
-      while ( t_iter > fBinningResolutionHalf-1 || t_iter < -fBinningResolutionHalf )
-        {
-          if (t_iter > fBinningResolutionHalf-1)
+      while ( t_iter > fBinningResolutionHalfHiGain-1 || t_iter < -fBinningResolutionHalfHiGain )
+        {
+          if (t_iter > fBinningResolutionHalfHiGain-1)
             {
-              t_iter -= fBinningResolution;
+              t_iter -= fBinningResolutionHiGain;
               max_p--; 
               sample_iter--;
             }
-          if (t_iter < -fBinningResolutionHalf)
+          if (t_iter < -fBinningResolutionHalfHiGain)
             {
-              t_iter += fBinningResolution;
+              t_iter += fBinningResolutionHiGain;
               max_p++; 
               sample_iter++;
@@ -377 +465 @@
       for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
       {
-        const Int_t   idx = fBinningResolution*sample + fBinningResolutionHalf + t_iter;
+        const Int_t   idx = fBinningResolutionHiGain*sample + fBinningResolutionHalfHiGain + t_iter;
         const Int_t   ids = max_p + sample;
         const Float_t pex = ids < 0 ? 0. : 
@@ -383 +471 @@
         sum              += fAmpWeightsHiGain [idx]*pex; 
         time_sum         += fTimeWeightsHiGain[idx]*pex;
-        //        cout << idx << " " << fw_amp[idx] << "  " << pex << "  " << t_iter << "  " << sum << "  " << ((Int_t(p2-ptr)+abflag) & 0x1) << "  " << PedMean[(Int_t(p2-ptr)+abflag) & 0x1] << endl;
       }
 
       if (sum != 0.)
-        time = max_p + 1. + sample_iter 
-             - ((Float_t)t_iter)/fBinningResolution - 0.4 - time_sum/sum 
-             + (Float_t)fHiGainFirst; 
+        time = max_p + fTimeShiftHiGain /* this shifts the time to the start of the rising edge */
+             - ((Float_t)t_iter)/fBinningResolutionHiGain - time_sum/sum;
       else 
         time = 0.;
@@ -404 +490 @@
 {
 
-  const Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1;
-  
-  Float_t time_sum  = 0;
-  Float_t fmax      = 0;
-  Float_t ftime_max = 0;
-  
-  Float_t pedes = ped.GetPedestal();
-
+  Int_t range = fLoGainLast - fLoGainFirst + 1;
+  const Byte_t *end = ptr + range;
+  Byte_t *p     = ptr;
+  Byte_t maxpos = 0;
+  Byte_t max    = 0;
+  Int_t count   = 0;
+
+  //
+  // Check for saturation in all other slices
+  //
+  while (p<end)
+    {
+
+      fLoGainSignal[count++] = (Float_t)*p;
+
+      if (*p > max)
+        {
+          max    = *p;
+          maxpos =  p-ptr;
+        }
+
+      if (*p++ >= fSaturationLimit)
+        sat++;
+    }
+  
+  Float_t time_sum  = 0.;
+  Float_t fmax      = 0.;
+  Float_t ftime_max = 0.;
+  Int_t   max_p     = 0;
+  
+  const Float_t pedes  = ped.GetPedestal();
   const Float_t ABoffs = ped.GetPedestalABoffset();
         
@@ -421 +530 @@
   // Calculate the sum of the first fWindowSize slices
   //
-  sat           = 0;
-  Byte_t *p     = ptr;
-  Byte_t *max_p = ptr;
-
-  while (p < end-fWindowSizeLoGain)
-    {
-      sum      = 0;
-      time_sum = 0;
-
+  for (Int_t i=0;i<range-fWindowSizeLoGain;i++)
+    {
+      sum      = 0.;
+      time_sum = 0.;
+      
       //
       // Slide with a window of size fWindowSizeLoGain over the sample 
       // and multiply the entries with the corresponding weights
       //
-      Byte_t *p1 = p++;
       for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
       {
-        const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf;
-        const Int_t adx = (Int_t(p1-ptr)+abflag) & 0x1;
-        sum      += fAmpWeightsLoGain [idx]*(*p1-PedMean[adx]); 
-        time_sum += fTimeWeightsLoGain[idx]*(*p1-PedMean[adx]);
-        p1++;
+        const Int_t   idx = fBinningResolutionLoGain*sample+fBinningResolutionHalfLoGain;
+        const Float_t pex = fLoGainSignal[sample+i]-PedMean[(sample+i+abflag) & 0x1];
+        sum              += fAmpWeightsLoGain [idx]*pex; 
+        time_sum         += fTimeWeightsLoGain[idx]*pex;
       }
 
@@ -448 +551 @@
         fmax      = sum;
         ftime_max = time_sum;
-        max_p     = p;
+        max_p     = i;
       }
-    } /* while (p<end-fWindowSizeLoGain) */
-  
-    if (fmax!=0)
+    } /*   for (Int_t i=0;i<range-fWindowSizeLoGain;i++) */
+  
+  if (fmax!=0)
     {
       ftime_max        /= fmax;
-      Int_t t_iter      = Int_t(ftime_max*fBinningResolution);
+      Int_t t_iter      = Int_t(ftime_max*fBinningResolutionLoGain);
       Int_t sample_iter = 0;
-      p                 = max_p;
-
-      while((t_iter)>fBinningResolutionHalf-1 || t_iter<-fBinningResolutionHalf)
-        {
-          if (t_iter>fBinningResolutionHalf-1)
+
+      while ( t_iter > fBinningResolutionHalfLoGain-1 || t_iter < -fBinningResolutionHalfLoGain )
+        {
+          if (t_iter > fBinningResolutionHalfLoGain-1)
             {
-              t_iter -= fBinningResolution;
-              p--;
+              t_iter -= fBinningResolutionLoGain;
+              max_p--; 
               sample_iter--;
             }
-          if (t_iter < -fBinningResolutionHalf)
+          if (t_iter < -fBinningResolutionHalfLoGain)
             {
-              t_iter += fBinningResolution;
-              p++; 
+              t_iter += fBinningResolutionLoGain;
+              max_p++; 
               sample_iter++;
             }
         }
       
-      sum       = 0;
-      Byte_t *p2 = p;
-     
+      sum = 0.;
       //
       // Slide with a window of size fWindowSizeLoGain over the sample 
@@ -484 +584 @@
       for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
       {
-        const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf+t_iter;
-        const Int_t adx = (Int_t(p2-ptr)+abflag) & 0x1;
-        sum      += fAmpWeightsLoGain [idx]*(*p2-PedMean[adx]); 
-        time_sum += fTimeWeightsLoGain[idx]*(*p2-PedMean[adx]);
-        p2++;
+        const Int_t   idx = fBinningResolutionLoGain*sample + fBinningResolutionHalfLoGain + t_iter;
+        const Int_t   ids = max_p + sample;
+        const Float_t pex = ids < 0 ? 0. : 
+          ( ids > range ? 0. : fLoGainSignal[ids]-PedMean[(ids+abflag) & 0x1]);
+        sum              += fAmpWeightsLoGain [idx]*pex; 
+        time_sum         += fTimeWeightsLoGain[idx]*pex;
       }
 
-      if (sum != 0)
-        time = (Float_t)(max_p - ptr) + 1. + sample_iter 
-             - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum 
-             + (Float_t)fLoGainFirst; 
+      if (sum != 0.)
+        time = max_p + fTimeShiftLoGain /* this shifts the time to the start of the rising edge */
+             - ((Float_t)t_iter)/fBinningResolutionLoGain - time_sum/sum;
       else 
         time = 0.;
@@ -527 +627 @@
   if (IsEnvDefined(env, prefix, "BinningResolution", print))
     {
-      SetBinningResolution(GetEnvValue(env, prefix, "BinningResolution", fBinningResolution));
+      SetBinningResolution(GetEnvValue(env, prefix, "BinningResolutionHiGain", fBinningResolutionHiGain),
+                           GetEnvValue(env, prefix, "BinningResolutionLoGain", fBinningResolutionLoGain));
       rc = kTRUE;
     }
@@ -536 +637 @@
 }
 
+Int_t  MExtractTimeAndChargeDigitalFilter::PostProcess()
+{
+  
+  *fLog << endl;
+  *fLog << inf << "Used High Gain weights in the extractor: " << endl;
+  
+  for (Int_t i=0;i<fBinningResolutionHiGain*fWindowSizeHiGain; i++)
+    *fLog << inf << fAmpWeightsHiGain[i] << "   " << fTimeWeightsHiGain[i] << endl;
+
+  *fLog << endl;
+  *fLog << inf << "Used Low Gain weights in the extractor: " << endl;
+  for (Int_t i=0;i<fBinningResolutionLoGain*fWindowSizeLoGain; i++)
+    *fLog << inf << fAmpWeightsLoGain[i] << "   " << fTimeWeightsLoGain[i] << endl;
+  
+  return kTRUE;
+}

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h

@@ -19 +19 @@
   static const Byte_t fgLoGainFirst;
   static const Byte_t fgLoGainLast;
-  static const Int_t  fgHiGainWindowSize;
-  static const Int_t  fgLoGainWindowSize;
-  static const Int_t  fgBinningResolution;
+  static const Int_t  fgWindowSizeHiGain;
+  static const Int_t  fgWindowSizeLoGain;
+  static const Int_t  fgBinningResolutionHiGain;
+  static const Int_t  fgBinningResolutionLoGain;
   
   MArrayF fHiGainSignal;               //! Need fast access to the signals in a float way
   MArrayF fLoGainSignal;               //! Store them in separate arrays
 
+  Float_t fTimeShiftHiGain;
+  Float_t fTimeShiftLoGain;
+  
   Int_t   fWindowSizeHiGain;            
   Int_t   fWindowSizeLoGain;            
 
-  Int_t   fBinningResolution;
-  Int_t   fBinningResolutionHalf;
+  Int_t   fBinningResolutionHiGain;
+  Int_t   fBinningResolutionHalfHiGain;
+  Int_t   fBinningResolutionLoGain;
+  Int_t   fBinningResolutionHalfLoGain;
   
   MArrayF fAmpWeightsHiGain; 
@@ -38 +44 @@
 
   Bool_t ReInit( MParList *pList );
-
+  Int_t PostProcess();  
+  
   Int_t  ReadEnv(const TEnv &env, TString prefix, Bool_t print);
   
@@ -54 +61 @@
   MExtractTimeAndChargeDigitalFilter(const char *name=NULL, const char *title=NULL);  
 
-  void SetWeightsFile(TString filename="pulpo_weights.root");
-  void SetWindowSize(Int_t windowh=fgHiGainWindowSize,
-                     Int_t windowl=fgLoGainWindowSize);
+  Bool_t WriteWeightsFile(TString filename="cosmic_weights.dat");
 
-  void SetBinningResolution(Int_t r=fgBinningResolution)  {
-    fBinningResolution     = r & ~1;
-    fBinningResolutionHalf = fBinningResolution/2; }
+  Bool_t ReadWeightsFile(TString filename="cosmic_weights.dat");
+  void SetWindowSize(Int_t windowh=fgWindowSizeHiGain,
+                     Int_t windowl=fgWindowSizeLoGain);
+
+  void SetBinningResolution(const Int_t rh=fgBinningResolutionHiGain, const Int_t rl=fgBinningResolutionLoGain)  {
+    fBinningResolutionHiGain     = rh & ~1;
+    fBinningResolutionHalfHiGain = fBinningResolutionHiGain/2;
+    fBinningResolutionLoGain     = rl & ~1;
+    fBinningResolutionHalfLoGain = fBinningResolutionLoGain/2;
+  }