Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc
===================================================================
--- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc	(revision 5161)
+++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc	(revision 5162)
@@ -56,7 +56,9 @@
 #include "MPedestalPix.h"
 
-#include "TFile.h" // remove!!
-#include "TH1F.h"  // remove!!
+#include "TFile.h" 
+#include "TH1F.h" 
+#include "TH2F.h" 
 #include "TString.h"
+#include "TMatrix.h"
 
 #include <fstream>
@@ -74,4 +76,6 @@
 const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
 const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
+const Int_t  MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain = 4;
+const Int_t  MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain = 4;
 // --------------------------------------------------------------------------
 //
@@ -94,4 +98,5 @@
   SetWindowSize();
   SetBinningResolution();
+  SetSignalStartBin();
 
   ReadWeightsFile("");
@@ -165,149 +170,4 @@
 }
 
-//----------------------------------------------------------------------------
-//
-// Read a pre-defined weights file into the class. 
-// This is mandatory for the extraction
-//
-// If filenname is empty, then all weights will be set to 1.
-//
-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())
-    {
-      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 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;
-}
-
 // --------------------------------------------------------------------------
 //
@@ -337,4 +197,5 @@
   return kTRUE;
 }
+
 
 void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Float_t &dsum, 
@@ -653,2 +514,382 @@
   return kTRUE;
 }
+
+//----------------------------------------------------------------------------
+//
+// Read a pre-defined weights file into the class. 
+// This is mandatory for the extraction
+//
+// If filenname is empty, then all weights will be set to 1.
+//
+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())
+    {
+      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 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;
+}
+
+//----------------------------------------------------------------------------
+//
+// Create the weights file
+// Beware that the shape-histogram has to contain the pulse starting at bin 1
+//
+Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi,
+                                                            TH1F *shapelo, TH2F *autocorrlo )
+{
+
+  const Int_t nbinshi = shapehi->GetNbinsX();
+  Float_t binwidth    = shapehi->GetBinWidth(1);
+
+  TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"),
+                                Form("%s%s",shapehi->GetTitle()," derivative"),
+                                nbinshi,
+                                shapehi->GetBinLowEdge(1),
+                                shapehi->GetBinLowEdge(nbinshi)+binwidth);
+
+  //
+  // Calculate the derivative of shapehi
+  // 
+  for (Int_t i = 1; i<nbinshi+1;i++)
+    {
+      derivativehi->SetBinContent(i,
+                                ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth));
+      derivativehi->SetBinError(i,
+                              (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1)
+                                    +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth));
+    }
+  
+  //
+  // normalize the shapehi, such that the integral for fWindowSize slices is one!
+  //
+  Float_t sum    = 0;
+  Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain);
+  lasttemp       = lasttemp > nbinshi ? nbinshi : lasttemp;
+  
+  for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; i<lasttemp; i++) {
+    sum += shapehi->GetBinContent(i);
+  }
+  sum /= fBinningResolutionHiGain;
+
+  shapehi->Scale(1./sum);
+  derivativehi->Scale(1./sum);
+
+  //
+  // read in the noise auto-correlation function:
+  //
+  TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain);
+
+  for (Int_t i=0; i<fWindowSizeHiGain; i++){
+    for (Int_t j=0; j<fWindowSizeHiGain; j++){
+      Bhi[i][j]=autocorrhi->GetBinContent(i+1+fSignalStartBinHiGain,j+1+fSignalStartBinHiGain);
+    }
+  }  
+  Bhi.Invert();
+
+  const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain;
+  fAmpWeightsHiGain.Set(nsizehi);
+  fTimeWeightsHiGain.Set(nsizehi);  
+  
+  //
+  // Loop over relative time in one BinningResolution interval
+  //
+  Int_t start = fBinningResolutionHiGain*fSignalStartBinHiGain + fBinningResolutionHalfHiGain;
+
+  for (Int_t i = -fBinningResolutionHalfHiGain+1; i<=fBinningResolutionHalfHiGain; i++)
+    {
+  
+      TMatrix g(fWindowSizeHiGain,1);
+      TMatrix gT(1,fWindowSizeHiGain);
+      TMatrix d(fWindowSizeHiGain,1);
+      TMatrix dT(1,fWindowSizeHiGain);
+    
+      for (Int_t count=0; count < fWindowSizeHiGain; count++){
+      
+        g[count][0]=shapehi->GetBinContent(start
+                                         +fBinningResolutionHiGain*count+i); 
+        gT[0][count]=shapehi->GetBinContent(start
+                                          +fBinningResolutionHiGain*count+i);
+        d[count][0]=derivativehi->GetBinContent(start
+                                              +fBinningResolutionHiGain*count+i);
+        dT[0][count]=derivativehi->GetBinContent(start
+                                               +fBinningResolutionHiGain*count+i);
+      }
+    
+      TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g));
+      Float_t   denom = m_denom[0][0];  // ROOT thinks, m_denom is still a matrix
+      
+      TMatrix m_first = dT*(Bhi*d);       // ROOT thinks, m_first is still a matrix
+      Float_t   first = m_first[0][0]/denom;
+      
+      TMatrix m_last  = gT*(Bhi*d);       // ROOT thinks, m_last  is still a matrix
+      Float_t   last  = m_last[0][0]/denom;
+      
+      TMatrix m1 = gT*Bhi;
+      m1 *= first;
+      
+      TMatrix m2 = dT*Bhi; 
+      m2 *=last;
+      
+      TMatrix w_amp = m1 - m2;
+      
+      TMatrix m_first1 = gT*(Bhi*g);
+      Float_t   first1 = m_first1[0][0]/denom;
+      
+      TMatrix m_last1  = gT*(Bhi*d);
+      Float_t   last1  = m_last1 [0][0]/denom;
+      
+      TMatrix m11 = dT*Bhi; 
+      m11 *=first1;
+      
+      TMatrix m21 = gT*Bhi;
+      m21 *=last1;
+      
+      TMatrix w_time= m11 - m21; 
+      
+      for (Int_t count=0; count < fWindowSizeHiGain; count++)
+        {
+          const Int_t idx = i+fBinningResolutionHalfHiGain+fBinningResolutionHiGain*count-1;
+          fAmpWeightsHiGain [idx] = w_amp [0][count];
+          fTimeWeightsHiGain[idx] = w_time[0][count];
+        }
+      
+    } // end loop over i
+
+  //
+  // Low Gain histograms
+  //
+  if (shapelo)
+    {
+      const Int_t nbinslo  = shapelo->GetNbinsX();
+      binwidth = shapelo->GetBinWidth(1);
+      
+      TH1F *derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"),
+                                    Form("%s%s",shapelo->GetTitle()," derivative"),
+                                    nbinslo,
+                                    shapelo->GetBinLowEdge(1),
+                                    shapelo->GetBinLowEdge(nbinslo)+binwidth);
+      
+      //
+      // Calculate the derivative of shapelo
+      // 
+      for (Int_t i = 1; i<nbinslo+1;i++)
+        {
+          derivativelo->SetBinContent(i,
+                                      ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth));
+          derivativelo->SetBinError(i,
+                                    (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1)
+                                          +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth));
+        }
+      
+      //
+      // normalize the shapelo, such that the integral for fWindowSize slices is one!
+      //
+      sum      = 0;
+      lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain);
+      lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp;
+      
+      for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; i<lasttemp; i++) 
+        sum += shapelo->GetBinContent(i);
+
+      sum /= fBinningResolutionLoGain;
+      
+      shapelo->Scale(1./sum);
+      derivativelo->Scale(1./sum);
+      
+      //
+      // read in the noise auto-correlation function:
+      //
+      TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain);
+      
+      for (Int_t i=0; i<fWindowSizeLoGain; i++){
+        for (Int_t j=0; j<fWindowSizeLoGain; j++){
+          Blo[i][j]=autocorrlo->GetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain);
+        }
+      }  
+      Blo.Invert();
+
+      const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain;
+      fAmpWeightsLoGain.Set(nsizelo);
+      fTimeWeightsLoGain.Set(nsizelo);  
+  
+      //
+      // Loop over relative time in one BinningResolution interval
+      //
+      Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionHalfLoGain;
+      
+      for (Int_t i = -fBinningResolutionHalfLoGain+1; i<=fBinningResolutionHalfLoGain; i++)
+        {
+          
+          TMatrix g(fWindowSizeLoGain,1);
+          TMatrix gT(1,fWindowSizeLoGain);
+          TMatrix d(fWindowSizeLoGain,1);
+          TMatrix dT(1,fWindowSizeLoGain);
+          
+          for (Int_t count=0; count < fWindowSizeLoGain; count++){
+            
+            g[count][0] = shapelo->GetBinContent(start
+                                             +fBinningResolutionLoGain*count+i); 
+            gT[0][count]= shapelo->GetBinContent(start
+                                              +fBinningResolutionLoGain*count+i);
+            d[count][0] = derivativelo->GetBinContent(start
+                                                  +fBinningResolutionLoGain*count+i);
+            dT[0][count]= derivativelo->GetBinContent(start
+                                                   +fBinningResolutionLoGain*count+i);
+          }
+          
+          TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g));
+          Float_t   denom = m_denom[0][0];  // ROOT thinks, m_denom is still a matrix
+          
+          TMatrix m_first = dT*(Blo*d);       // ROOT thinks, m_first is still a matrix
+          Float_t   first = m_first[0][0]/denom;
+          
+          TMatrix m_last  = gT*(Blo*d);       // ROOT thinks, m_last  is still a matrix
+          Float_t   last  = m_last[0][0]/denom;
+          
+          TMatrix m1 = gT*Blo;
+          m1 *= first;
+          
+          TMatrix m2 = dT*Blo; 
+          m2 *=last;
+          
+          TMatrix w_amp = m1 - m2;
+          
+          TMatrix m_first1 = gT*(Blo*g);
+          Float_t   first1 = m_first1[0][0]/denom;
+          
+          TMatrix m_last1  = gT*(Blo*d);
+          Float_t   last1  = m_last1 [0][0]/denom;
+          
+          TMatrix m11 = dT*Blo; 
+          m11 *=first1;
+          
+          TMatrix m21 = gT*Blo;
+          m21 *=last1;
+          
+          TMatrix w_time= m11 - m21; 
+          
+          for (Int_t count=0; count < fWindowSizeLoGain; count++)
+            {
+              const Int_t idx = i+fBinningResolutionHalfLoGain+fBinningResolutionLoGain*count-1;
+              fAmpWeightsLoGain [idx] = w_amp [0][count];
+              fTimeWeightsLoGain[idx] = w_time[0][count];
+            }
+      
+        } // end loop over i
+    }
+  
+
+  ofstream fn(filename.Data());
+
+  fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl;
+  fn << "# (Amplitude)  (Time) " << endl;
+
+  for (Int_t i=0; i<nsizehi; i++)
+    fn << "\t" << fAmpWeightsHiGain[i] << "\t" << fTimeWeightsHiGain[i] << endl;
+
+  fn << "# Low Gain Weights: " << fWindowSizeLoGain << " " << fBinningResolutionLoGain << endl;
+  fn << "# (Amplitude)  (Time) " << endl;
+
+  for (Int_t i=0; i<nsizehi; i++)
+    fn << "\t" << fAmpWeightsLoGain[i] << "\t" << fTimeWeightsLoGain[i] << endl;
+
+  return kTRUE;
+}
Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h
===================================================================
--- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h	(revision 5161)
+++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h	(revision 5162)
@@ -10,4 +10,6 @@
 #endif
 
+class TH1F;
+class TH2F;
 class MPedestalPix;
 class MExtractTimeAndChargeDigitalFilter : public MExtractTimeAndCharge
@@ -23,11 +25,16 @@
   static const Int_t  fgBinningResolutionHiGain;
   static const Int_t  fgBinningResolutionLoGain;
+  static const Int_t  fgSignalStartBinHiGain;
+  static const Int_t  fgSignalStartBinLoGain;
   
-  MArrayF fHiGainSignal;               //! Need fast access to the signals in a float way
-  MArrayF fLoGainSignal;               //! Store them in separate arrays
+  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   fSignalStartBinHiGain;
+  Int_t   fSignalStartBinLoGain;
+
   Int_t   fWindowSizeHiGain;            
   Int_t   fWindowSizeLoGain;            
@@ -61,5 +68,7 @@
   MExtractTimeAndChargeDigitalFilter(const char *name=NULL, const char *title=NULL);  
 
-  Bool_t WriteWeightsFile(TString filename="cosmic_weights.dat");
+  Bool_t WriteWeightsFile(TString filename,
+                          TH1F *shapehi, TH2F *autocorrhi,
+                          TH1F *shapelo=NULL, TH2F *autocorrlo=NULL );
 
   Bool_t ReadWeightsFile(TString filename="cosmic_weights.dat");
@@ -74,4 +83,9 @@
   }
   
+
+  void SetSignalStartBin( const Int_t sh=fgSignalStartBinHiGain, const Int_t sl=fgSignalStartBinLoGain) {
+    fSignalStartBinHiGain = sh;
+    fSignalStartBinLoGain = sl; 
+  }
   
   ClassDef(MExtractTimeAndChargeDigitalFilter, 0)   // Hendrik's digital filter