Changeset 5152 for trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc

Timestamp:

09/30/04 15:08:58 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (modified) (10 diffs)

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc

@@ -80 +80 @@
 // - SetBinningResolution();
 //
+// Sets all weights to 1.
+//
 MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title) 
 {
@@ -90 +92 @@
   SetBinningResolution();
 
-  for (int i=0; i<60; i++){
-    fw_amp[i]  = 1;
-    fw_time[i] = 1;
-  }
+  SetWeightsFile("");
 }
 
@@ -161 +160 @@
 // This is mandatory for the extraction
 //
+// 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); 
+  
+  if (filename.IsNull())
+    {
+      for (UInt_t i=0; i<fAmpWeightsHiGain.GetSize(); i++)
+        {
+          fAmpWeightsHiGain [i] = 1.;
+          fTimeWeightsHiGain[i] = 1.;
+        }
+      for (UInt_t i=0; i<fAmpWeightsLoGain.GetSize(); i++)
+        {
+          fAmpWeightsLoGain [i] = 1.;
+          fTimeWeightsLoGain[i] = 1.;
+        }
+      return;
+    }
 
   TFile *f = new TFile(filename);
@@ -182 +203 @@
     }
   
-  for (int i=0; i<fBinningResolution*fWindowSizeHiGain;i++)
-    {
-    fw_amp [i] = hw_amp->GetBinContent(i+1);
-    fw_time[i] = hw_time->GetBinContent(i+1);
-    //    y[i]=hshape->GetBinContent(i+1);
-    //    derivative[i]=hderivative->GetBinContent(i+1);
-  
-    //    cout << "w_amp[" << i << "] = " << fw_amp[i] << endl;
-  }
+  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;
+}
+
+// --------------------------------------------------------------------------
+//
+// ReInit
+//
+// Calls:
+// - MExtractor::ReInit(pList);
+// - Creates new arrays according to the extraction range
+//
+Bool_t MExtractTimeAndChargeDigitalFilter::ReInit(MParList *pList)
+{
+
+  if (!MExtractTimeAndCharge::ReInit(pList))
+    return kFALSE;
+
+  Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
+
+  fHiGainSignal.Set(range);
+
+  range = fLoGainLast - fLoGainFirst + 1;
+
+  fLoGainSignal.Set(range);
+
+  return kTRUE;
 }
 
@@ -199 +247 @@
 {
 
-  const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
+  Int_t range = fHiGainLast - fHiGainFirst + 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)
+    {
+
+      fHiGainSignal[count++] = (Float_t)*p;
+
+      if (*p > max)
+        {
+          max    = *p;
+          maxpos =  p-ptr;
+        }
+
+      if (*p++ >= fSaturationLimit)
+        sat++;
+    }
+  
+  if (fHiLoLast != 0)
+    {
+
+      end = logain + fHiLoLast;
+
+      while (logain<end)
+        {
+
+          fHiGainSignal[count++] = (Float_t)*logain;
+          range++;
+
+          if (*logain > max)
+            {
+              max    = *logain;
+              maxpos =  range;
+            }
+          
+          if (*logain++ >= fSaturationLimit)
+            sat++;
+        }
+    }
+  
+  //
+  // allow one saturated slice 
+  //
+  if (sat > 0)
+    return;
+
+#if 0
+  if (maxpos < 2)
+    {
+      time = (Float_t)fHiGainFirst;
+      return;
+    }
+#endif
   
   Float_t time_sum  = 0.;
   Float_t fmax      = 0.;
   Float_t ftime_max = 0.;
-  
-  Float_t pedes = ped.GetPedestal();
-
+  Int_t   max_p     = 0;
+  
+  const Float_t pedes  = ped.GetPedestal();
   const Float_t ABoffs = ped.GetPedestalABoffset();
         
@@ -216 +323 @@
   // Calculate the sum of the first fWindowSize slices
   //
-  Byte_t *p     = ptr;
-  Byte_t *max_p = ptr;
-
-  while (p < end-fWindowSizeHiGain)
+  for (Int_t i=0;i<range-fWindowSizeHiGain;i++)
     {
       sum      = 0.;
       time_sum = 0.;
-
+      
       //
       // Slide with a window of size fWindowSizeHiGain over the sample 
       // and multiply the entries with the corresponding weights
       //
-      Byte_t *p1 = p;
       for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
       {
         const Int_t   idx = fBinningResolution*sample+fBinningResolutionHalf;
-        const Float_t pex = (Float_t)*p1-PedMean[(Int_t(p1-ptr)+abflag) & 0x1];
-        sum      += fw_amp [idx]*pex; 
-        time_sum += fw_time[idx]*pex;
+        const Float_t pex = fHiGainSignal[sample+i]-PedMean[(sample+i+abflag) & 0x1];
+        sum              += fAmpWeightsHiGain [idx]*pex; 
+        time_sum         += fTimeWeightsHiGain[idx]*pex;
+      }
+
+      if (sum>fmax)
+      {
+        fmax      = sum;
+        ftime_max = time_sum;
+        max_p     = i;
+      }
+    } /*   for (Int_t i=0;i<range-fWindowSizeHiGain;i++) */
+  
+  if (fmax!=0)
+    {
+      ftime_max        /= fmax;
+      Int_t t_iter      = Int_t(ftime_max*fBinningResolution);
+      Int_t sample_iter = 0;
+
+      while ( t_iter > fBinningResolutionHalf-1 || t_iter < -fBinningResolutionHalf )
+        {
+          if (t_iter > fBinningResolutionHalf-1)
+            {
+              t_iter -= fBinningResolution;
+              max_p--; 
+              sample_iter--;
+            }
+          if (t_iter < -fBinningResolutionHalf)
+            {
+              t_iter += fBinningResolution;
+              max_p++; 
+              sample_iter++;
+            }
+        }
+      
+      sum = 0.;
+      //
+      // Slide with a window of size fWindowSizeHiGain over the sample 
+      // and multiply the entries with the corresponding weights
+      //
+      for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
+      {
+        const Int_t   idx = fBinningResolution*sample + fBinningResolutionHalf + t_iter;
+        const Int_t   ids = max_p + sample;
+        const Float_t pex = ids < 0 ? 0. : 
+          ( ids > range ? 0. : fHiGainSignal[ids]-PedMean[(ids+abflag) & 0x1]);
+        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; 
+      else 
+        time = 0.;
+    } /* if (max!=0) */ 
+    else 
+      time=0.;
+    
+  return;
+}
+
+void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum, 
+                                                                 Float_t &time, Float_t &dtime,
+                                                                 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
+{
+
+  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();
+
+  const Float_t ABoffs = ped.GetPedestalABoffset();
+        
+  Float_t PedMean[2];
+  PedMean[0] = pedes + ABoffs;
+  PedMean[1] = pedes - ABoffs;
+
+  //
+  // 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;
+
+      //
+      // 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++;
       }
@@ -244 +450 @@
         max_p     = p;
       }
-      p++;
-    } /* while (p<ptr+fWindowSizeHiGain-fHiLoLast) */
-  
-#if 0
-  //
-  // If part of the "low-Gain" slices are used, 
-  // repeat steps one and two for the logain part until fHiLoLast
-  //
-  Byte_t *l = logain;
-  while (l<logain+fHiLoLast)
-    {
-      sum      = 0;
-      time_sum = 0;
-
-      //
-      // Slide with a window of size fWindowSizeHiGain over the sample 
-      // and multiply the entries with the corresponding weights
-      //
-      Byte_t *p1 = l++;
-      for (Int_t sample=0; sample < fHiLoLast; sample++)
-      {
-        const Int_t   idx = fBinningResolution*sample+fBinningResolutionHalf;
-        const Float_t pex = (Float_t)*p1-PedMean[(Int_t(p1-logain)+abflag) & 0x1];
-        sum      += fw_amp [idx]*pex; 
-        time_sum += fw_time[idx]*pex;
-        p1++;
-      }
-
-      if (sum>fmax)
-      {
-        fmax      = sum;
-        ftime_max = time_sum;
-        max_p     = p;
-      }
-
-    }
-#endif  
+    } /* while (p<end-fWindowSizeLoGain) */
   
     if (fmax!=0)
@@ -305 +475 @@
         }
       
-      sum        = 0.;
+      sum       = 0;
       Byte_t *p2 = p;
      
       //
-      // Slide with a window of size fWindowSizeHiGain over the sample 
+      // Slide with a window of size fWindowSizeLoGain over the sample 
       // and multiply the entries with the corresponding weights
       //
-      for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
+      for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
       {
-        const Int_t   idx = fBinningResolution*sample + fBinningResolutionHalf + t_iter;
-        const Float_t pex = (Float_t)*p2-PedMean[(Int_t(p2-ptr)+abflag) & 0x1];
-        sum      += fw_amp [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;
-        time_sum += fw_time[idx]*pex;
+        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++;
       }
@@ -324 +493 @@
       if (sum != 0)
         time = (Float_t)(max_p - ptr) + 1. + sample_iter 
-             - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum ; 
+             - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum 
+             + (Float_t)fLoGainFirst; 
       else 
         time = 0.;
@@ -334 +504 @@
 }
 
-void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum, 
-                                                                 Float_t &time, Float_t &dtime,
-                                                                 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
-{
-
-  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();
-
-  const Float_t ABoffs = ped.GetPedestalABoffset();
-        
-  Float_t PedMean[2];
-  PedMean[0] = pedes + ABoffs;
-  PedMean[1] = pedes - ABoffs;
-
-  //
-  // 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;
-
-      //
-      // 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      += fw_amp [idx]*(*p1-PedMean[adx]); 
-        time_sum += fw_time[idx]*(*p1-PedMean[adx]);
-        p1++;
-      }
-
-      if (sum>fmax)
-      {
-        fmax      = sum;
-        ftime_max = time_sum;
-        max_p     = p;
-      }
-    } /* while (p<end-fWindowSizeLoGain) */
-  
-    if (fmax!=0)
-    {
-      ftime_max        /= fmax;
-      Int_t t_iter      = Int_t(ftime_max*fBinningResolution);
-      Int_t sample_iter = 0;
-      p                 = max_p;
-
-      while((t_iter)>fBinningResolutionHalf-1 || t_iter<-fBinningResolutionHalf)
-        {
-          if (t_iter>fBinningResolutionHalf-1)
-            {
-              t_iter -= fBinningResolution;
-              p--;
-              sample_iter--;
-            }
-          if (t_iter < -fBinningResolutionHalf)
-            {
-              t_iter += fBinningResolution;
-              p++; 
-              sample_iter++;
-            }
-        }
-      
-      sum       = 0;
-      Byte_t *p2 = p;
-     
-      //
-      // Slide with a window of size fWindowSizeLoGain over the sample 
-      // and multiply the entries with the corresponding weights
-      //
-      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      += fw_amp [idx]*(*p2-PedMean[adx]); 
-        time_sum += fw_time[idx]*(*p2-PedMean[adx]);
-        p2++;
-      }
-
-      if (sum != 0)
-        time = (Float_t)(max_p - ptr) + 1. + sample_iter 
-             - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum ; 
-      else 
-        time = 0.;
-    } /* if (max!=0) */ 
-    else 
-      time=0.;
-    
-  return;
-}
-
 Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
 {