Ignore:
Timestamp:
10/25/06 19:36:26 (18 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc

    r8154 r8165  
    11/* ======================================================================== *\
    2 ! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeDigitalFilter.cc,v 1.73 2006-10-24 08:24:52 tbretz Exp $
     2! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeDigitalFilter.cc,v 1.74 2006-10-25 18:36:26 tbretz Exp $
    33! --------------------------------------------------------------------------
    44!
    … …  
    6262#include <fstream>
    6363
    64 #include <TH1.h>
    65 #include <TH2.h>
    66 #include <TMatrix.h>
     64#include <TRandom.h>
    6765
    6866#include "MLog.h"
    … …  
    7674#include "MExtralgoDigitalFilter.h"
    7775
    78 #include "MPedestalPix.h"
    79 
    8076ClassImp(MExtractTimeAndChargeDigitalFilter);
    8177
    8278using namespace std;
    8379
    84 const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst             =  0;
    85 const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast              = 15;
    86 const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst             =  1;
    87 const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast              = 14;
    88 //const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain        =  4;
    89 //const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain        =  6;
    90 const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
    91 const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
    92 const Int_t  MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain    =  4;
    93 const Int_t  MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain    =  4;
    94 const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain           =  0.95;
    95 //const Float_t MExtractTimeAndChargeDigitalFilter::fgLoGainStartShift       = -1.8;
     80const Byte_t  MExtractTimeAndChargeDigitalFilter::fgHiGainFirst             =  0;
     81const Byte_t  MExtractTimeAndChargeDigitalFilter::fgHiGainLast              = 15;
     82const Byte_t  MExtractTimeAndChargeDigitalFilter::fgLoGainFirst             =  1;
     83const Byte_t  MExtractTimeAndChargeDigitalFilter::fgLoGainLast              = 14;
     84const Int_t   MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
     85const Int_t   MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
     86const Int_t   MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain    =  4;
     87const Int_t   MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain    =  4;
     88const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain            =  0.95;
    9689
    9790// --------------------------------------------------------------------------
    … …  
    109102    : fBinningResolutionHiGain(fgBinningResolutionHiGain),
    110103    fBinningResolutionLoGain(fgBinningResolutionLoGain),
    111     fAutomaticWeights(kTRUE), fRandomIter(0)
     104    fAutomaticWeights(kTRUE)
    112105{
    113106    fName  = name  ? name  : "MExtractTimeAndChargeDigitalFilter";
    … …  
    116109    SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
    117110    SetWindowSize(3, 5);
    118 //    SetBinningResolution();
    119 //    SetSignalStartBin();
    120 
    121     SetOffsetLoGain(fgOffsetLoGain);          // Order between both
    122 //    SetLoGainStartShift(fgLoGainStartShift);  // is important
     111    SetOffsetLoGain(fgOffsetLoGain);
    123112}
    124113
    … …  
    189178  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
    190179  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
    191 
    192180}
    193181
    … …  
    256244// Gets called in the ReInit() and initialized the arrays
    257245//
    258 Bool_t MExtractTimeAndChargeDigitalFilter::InitArrays()
     246Bool_t MExtractTimeAndChargeDigitalFilter::InitArrays(Int_t n)
    259247{
    260248    if (!fRunHeader)
    261249        return kFALSE;
    262 
    263     //const Int_t rangehi = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
    264     //const Int_t rangelo = (Int_t)(fLoGainLast - fLoGainFirst + 1);
    265 
    266     //cout << "ARRAYS INITIALIZED TO : " << rangehi << " " << rangelo << endl;
    267 
    268     //fHiGainSignal.Set(rangehi);
    269     //fLoGainSignal.Set(rangelo);
    270250
    271251    return GetWeights();
    … …  
    296276                                                                  Float_t &sum, Float_t &dsum,
    297277                                                                  Float_t &time, Float_t &dtime,
    298                                                                   Byte_t sat, Int_t maxpos)
     278                                                                  Byte_t sat, Int_t maxpos) const
    299279{
    300280    // Do some handling if maxpos is last slice!
    … …  
    308288    if (IsNoiseCalculation())
    309289    {
    310         if (fRandomIter == fBinningResolutionHiGain)
    311             fRandomIter = 0;
    312 
    313         sum = df.ExtractNoise(fRandomIter);
    314         fRandomIter++;
     290        sum = df.ExtractNoise(gRandom->Integer(fBinningResolutionHiGain));
    315291        return;
    316292    }
    … …  
    324300                                                                  Float_t &sum, Float_t &dsum,
    325301                                                                  Float_t &time, Float_t &dtime,
    326                                                                   Byte_t sat, Int_t maxpos)
     302                                                                  Byte_t sat, Int_t maxpos) const
    327303{
    328304    MExtralgoDigitalFilter df(fBinningResolutionLoGain, fWindowSizeLoGain,
    … …  
    335311    if (IsNoiseCalculation())
    336312    {
    337         if (fRandomIter == fBinningResolutionLoGain)
    338             fRandomIter = 0;
    339 
    340         sum = df.ExtractNoise(fRandomIter);
     313        sum = df.ExtractNoise(gRandom->Integer(fBinningResolutionHiGain));
    341314        return;
    342315    }
    … …  
    347320}
    348321
    349 /*
     322
    350323// --------------------------------------------------------------------------
    351324//
    352 // Apply the digital filter algorithm to the high-gain slices.
    353 //
    354 void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
    355                                                                  Float_t &time, Float_t &dtime,
    356                                                                  Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
    357 {
    358     Int_t range = fHiGainLast - fHiGainFirst + 1;
    359     const Byte_t *end = first + range;
    360     Byte_t       *p  = first;
    361 
    362     const Float_t pedes  = ped.GetPedestal();
    363     const Float_t ABoffs = ped.GetPedestalABoffset();
    364 
    365     const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
    366 
    367     Byte_t maxcont = 0;
    368     Int_t  maxpos  = 0;  // obsolete for Digital Filter (used in spline!)
    369     sat = 0;
    370 
    371     //
    372     // Check for saturation in all other slices
    373     //
    374     Int_t ids = fHiGainFirst;
    375     Float_t *sample = fHiGainSignal.GetArray();
    376 
    377     while (p<end)
    378     {
    379         *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
    380 
    381         if (*p > maxcont)
    382         {
    383             maxpos = ids-fHiGainFirst-1;
    384             // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
    385             //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
    386             if (maxpos > 1 && maxpos < (range-1*//* - fWindowSizeHiGain + 1*//*))
    387                 maxcont = *p;
    388         }
    389 
    390         if (*p++ >= fSaturationLimit)
    391             if (!sat)
    392                 sat = ids-fHiGainFirst;
    393     }
    394 
    395     if (fHiLoLast != 0)
    396     {
    397         end = logain + fHiLoLast;
    398 
    399         while (logain<end)
    400         {
    401             *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
    402 
    403             if (*logain > maxcont)
    404             {
    405                 maxpos = ids-fHiGainFirst-1;
    406 
    407                 // FIXME!!! MUST BE THERE!
    408                 // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
    409                 //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
    410                 //    maxcont = *logain;
    411             }
    412 
    413             if (*logain++ >= fSaturationLimit)
    414                 if (!sat)
    415                     sat = ids-fHiGainFirst;
    416 
    417             range++;
    418         }
    419     }
    420 
    421     FindTimeAndChargeHiGain2(fHiGainSignal.GetArray(), range,
    422                              sum, dsum, time, dtime,
    423                              sat, 0);
    424 }
    425 
    426 // --------------------------------------------------------------------------
    427 //
    428 // Apply the digital filter algorithm to the low-gain slices.
    429 //
    430 void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum,
    431                                                                  Float_t &time, Float_t &dtime,
    432                                                                  Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
    433 {
    434     const Int_t range = fLoGainLast - fLoGainFirst + 1;
    435 
    436     const Byte_t *end = ptr + range;
    437     Byte_t *p     = ptr;
    438 
    439     //
    440     // Prepare the low-gain pedestal
    441     //
    442     const Float_t pedes  = ped.GetPedestal();
    443     const Float_t ABoffs = ped.GetPedestalABoffset();
    444 
    445     const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
    446 
    447     //
    448     // Check for saturation in all other slices
    449     //
    450     Float_t *sample = fLoGainSignal.GetArray();
    451     Int_t    ids    = fLoGainFirst;
    452 
    453     while (p<end)
    454     {
    455         *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
    456 
    457         if (*p++ >= fSaturationLimit)
    458             sat++;
    459     }
    460 
    461     FindTimeAndChargeLoGain2(fLoGainSignal.GetArray(), range,
    462                              sum, dsum, time, dtime,
    463                              sat, 0);
    464 
    465 }
    466 */
    467 // --------------------------------------------------------------------------
    468 //
    469325// Read the setup from a TEnv, eg:
    470 //   MJPedestal.MExtractor.WindowSizeHiGain: 6
    471 //   MJPedestal.MExtractor.WindowSizeLoGain: 6
    472 //   MJPedestal.MExtractor.BinningResolutionHiGain: 10
    473 //   MJPedestal.MExtractor.BinningResolutionLoGain: 10
    474326//   MJPedestal.MExtractor.WeightsFile: filename
    475327//   MJPedestal.MExtractor.AutomaticWeights: off
    … …  
    485337      rc = kTRUE;
    486338  }
    487    /*
    488   Byte_t hw = fWindowSizeHiGain;
    489   Byte_t lw = fWindowSizeLoGain;
    490 
    491   if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
    492     {
    493       hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
    494       rc = kTRUE;
    495     }
    496   if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
    497     {
    498       lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
    499       rc = kTRUE;
    500     }
    501  
    502   if (rc)
    503     SetWindowSize(hw, lw);
    504 
    505   Bool_t rc2 = kFALSE;
    506   Int_t brh = fBinningResolutionHiGain;
    507   Int_t brl = fBinningResolutionLoGain;
    508  
    509   if (IsEnvDefined(env, prefix, "BinningResolutionHiGain", print))
    510     {
    511       brh = GetEnvValue(env, prefix, brh);
    512       rc2 = kTRUE;
    513     }
    514   if (IsEnvDefined(env, prefix, "BinningResolutionLoGain", print))
    515     {
    516       brl = GetEnvValue(env, prefix, brl);
    517       rc2 = kTRUE;
    518     }
    519  
    520   if (rc2)
    521     {
    522       SetBinningResolution(brh, brl);
    523       rc = kTRUE;
    524     }
    525    */
     339
    526340  if (IsEnvDefined(env, prefix, "WeightsFile", print))
    527341    {
    … …  
    859673    return ReadWeightsFile(name, path);
    860674}
    861 /*
     675
    862676//----------------------------------------------------------------------------
    863677//
    864 // Create the weights file
    865 // Beware that the shape-histogram has to contain the pulse starting at bin 1
    866 //
    867 Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi,
    868                                                             TH1F *shapelo, TH2F *autocorrlo )
    869 {
    870 
    871   const Int_t nbinshi = shapehi->GetNbinsX();
    872   Float_t binwidth    = shapehi->GetBinWidth(1);
    873 
    874   TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"),
    875                                 Form("%s%s",shapehi->GetTitle()," derivative"),
    876                                 nbinshi,
    877                                 shapehi->GetBinLowEdge(1),
    878                                 shapehi->GetBinLowEdge(nbinshi)+binwidth);
    879 
    880   //
    881   // Calculate the derivative of shapehi
    882   //
    883   for (Int_t i = 1; i<nbinshi+1;i++)
    884     {
    885       derivativehi->SetBinContent(i,
    886                                 ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth));
    887       derivativehi->SetBinError(i,
    888                               (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1)
    889                                     +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth));
    890     }
    891  
    892   //
    893   // normalize the shapehi, such that the integral for fWindowSize slices is one!
    894   //
    895   Float_t sum    = 0;
    896   Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain);
    897   lasttemp       = lasttemp > nbinshi ? nbinshi : lasttemp;
    898  
    899   for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; i<lasttemp; i++) {
    900     sum += shapehi->GetBinContent(i);
    901   }
    902   sum /= fBinningResolutionHiGain;
    903 
    904   shapehi->Scale(1./sum);
    905   derivativehi->Scale(1./sum);
    906 
    907   //
    908   // read in the noise auto-correlation function:
    909   //
    910   TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain);
    911 
    912   for (Int_t i=0; i<fWindowSizeHiGain; i++){
    913     for (Int_t j=0; j<fWindowSizeHiGain; j++){
    914       Bhi[i][j]=autocorrhi->GetBinContent(i+1,j+1); //+fSignalStartBinHiGain +fSignalStartBinHiGain
    915     }
    916   } 
    917   Bhi.Invert();
    918 
    919   const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain;
    920   fAmpWeightsHiGain.Set(nsizehi);
    921   fTimeWeightsHiGain.Set(nsizehi); 
    922  
    923   //
    924   // Loop over relative time in one BinningResolution interval
    925   //
    926   Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1);
    927 
    928   for (Int_t i = -fBinningResolutionHiGain/2+1; i<=fBinningResolutionHiGain/2; i++)
    929     {
    930  
    931       TMatrix g(fWindowSizeHiGain,1);
    932       TMatrix gT(1,fWindowSizeHiGain);
    933       TMatrix d(fWindowSizeHiGain,1);
    934       TMatrix dT(1,fWindowSizeHiGain);
    935    
    936       for (Int_t count=0; count < fWindowSizeHiGain; count++){
    937      
    938         g[count][0]=shapehi->GetBinContent(start
    939                                          +fBinningResolutionHiGain*count+i);
    940         gT[0][count]=shapehi->GetBinContent(start
    941                                           +fBinningResolutionHiGain*count+i);
    942         d[count][0]=derivativehi->GetBinContent(start
    943                                               +fBinningResolutionHiGain*count+i);
    944         dT[0][count]=derivativehi->GetBinContent(start
    945                                                +fBinningResolutionHiGain*count+i);
    946       }
    947    
    948       TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g));
    949       Float_t   denom = m_denom[0][0];  // ROOT thinks, m_denom is still a matrix
    950      
    951       TMatrix m_first = dT*(Bhi*d);       // ROOT thinks, m_first is still a matrix
    952       Float_t   first = m_first[0][0]/denom;
    953      
    954       TMatrix m_last  = gT*(Bhi*d);       // ROOT thinks, m_last  is still a matrix
    955       Float_t   last  = m_last[0][0]/denom;
    956      
    957       TMatrix m1 = gT*Bhi;
    958       m1 *= first;
    959      
    960       TMatrix m2 = dT*Bhi;
    961       m2 *=last;
    962      
    963       TMatrix w_amp = m1 - m2;
    964      
    965       TMatrix m_first1 = gT*(Bhi*g);
    966       Float_t   first1 = m_first1[0][0]/denom;
    967      
    968       TMatrix m_last1  = gT*(Bhi*d);
    969       Float_t   last1  = m_last1 [0][0]/denom;
    970      
    971       TMatrix m11 = dT*Bhi;
    972       m11 *=first1;
    973      
    974       TMatrix m21 = gT*Bhi;
    975       m21 *=last1;
    976      
    977       TMatrix w_time= m11 - m21;
    978      
    979       for (Int_t count=0; count < fWindowSizeHiGain; count++)
    980         {
    981           const Int_t idx = i+fBinningResolutionHiGain/2+fBinningResolutionHiGain*count-1;
    982           fAmpWeightsHiGain [idx] = w_amp [0][count];
    983           fTimeWeightsHiGain[idx] = w_time[0][count];
    984         }
    985      
    986     } // end loop over i
    987 
    988   //
    989   // Low Gain histograms
    990   //
    991   TH1F *derivativelo = NULL;
    992   if (shapelo)
    993     {
    994       const Int_t nbinslo  = shapelo->GetNbinsX();
    995       binwidth = shapelo->GetBinWidth(1);
    996      
    997       derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"),
    998                               Form("%s%s",shapelo->GetTitle()," derivative"),
    999                               nbinslo,
    1000                               shapelo->GetBinLowEdge(1),
    1001                               shapelo->GetBinLowEdge(nbinslo)+binwidth);
    1002      
    1003       //
    1004       // Calculate the derivative of shapelo
    1005       //
    1006       for (Int_t i = 1; i<nbinslo+1;i++)
    1007         {
    1008           derivativelo->SetBinContent(i,
    1009                                       ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth));
    1010           derivativelo->SetBinError(i,
    1011                                     (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1)
    1012                                           +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth));
    1013         }
    1014      
    1015       //
    1016       // normalize the shapelo, such that the integral for fWindowSize slices is one!
    1017       //
    1018       sum      = 0;
    1019       lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain);
    1020       lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp;
    1021      
    1022       for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; i<lasttemp; i++)
    1023         sum += shapelo->GetBinContent(i);
    1024 
    1025       sum /= fBinningResolutionLoGain;
    1026      
    1027       shapelo->Scale(1./sum);
    1028       derivativelo->Scale(1./sum);
    1029      
    1030       //
    1031       // read in the noise auto-correlation function:
    1032       //
    1033       TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain);
    1034      
    1035       for (Int_t i=0; i<fWindowSizeLoGain; i++){
    1036         for (Int_t j=0; j<fWindowSizeLoGain; j++){
    1037           Blo[i][j]=autocorrlo->GetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain);
    1038         }
    1039       } 
    1040       Blo.Invert();
    1041 
    1042       const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain;
    1043       fAmpWeightsLoGain.Set(nsizelo);
    1044       fTimeWeightsLoGain.Set(nsizelo); 
    1045  
    1046       //
    1047       // Loop over relative time in one BinningResolution interval
    1048       //
    1049       Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionLoGain/2;
    1050      
    1051       for (Int_t i = -fBinningResolutionLoGain/2+1; i<=fBinningResolutionLoGain/2; i++)
    1052         {
    1053          
    1054           TMatrix g(fWindowSizeLoGain,1);
    1055           TMatrix gT(1,fWindowSizeLoGain);
    1056           TMatrix d(fWindowSizeLoGain,1);
    1057           TMatrix dT(1,fWindowSizeLoGain);
    1058          
    1059           for (Int_t count=0; count < fWindowSizeLoGain; count++){
    1060            
    1061             g[count][0] = shapelo->GetBinContent(start
    1062                                              +fBinningResolutionLoGain*count+i);
    1063             gT[0][count]= shapelo->GetBinContent(start
    1064                                               +fBinningResolutionLoGain*count+i);
    1065             d[count][0] = derivativelo->GetBinContent(start
    1066                                                   +fBinningResolutionLoGain*count+i);
    1067             dT[0][count]= derivativelo->GetBinContent(start
    1068                                                    +fBinningResolutionLoGain*count+i);
    1069           }
    1070          
    1071           TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g));
    1072           Float_t   denom = m_denom[0][0];  // ROOT thinks, m_denom is still a matrix
    1073          
    1074           TMatrix m_first = dT*(Blo*d);       // ROOT thinks, m_first is still a matrix
    1075           Float_t   first = m_first[0][0]/denom;
    1076          
    1077           TMatrix m_last  = gT*(Blo*d);       // ROOT thinks, m_last  is still a matrix
    1078           Float_t   last  = m_last[0][0]/denom;
    1079          
    1080           TMatrix m1 = gT*Blo;
    1081           m1 *= first;
    1082          
    1083           TMatrix m2 = dT*Blo;
    1084           m2 *=last;
    1085          
    1086           TMatrix w_amp = m1 - m2;
    1087          
    1088           TMatrix m_first1 = gT*(Blo*g);
    1089           Float_t   first1 = m_first1[0][0]/denom;
    1090          
    1091           TMatrix m_last1  = gT*(Blo*d);
    1092           Float_t   last1  = m_last1 [0][0]/denom;
    1093          
    1094           TMatrix m11 = dT*Blo;
    1095           m11 *=first1;
    1096          
    1097           TMatrix m21 = gT*Blo;
    1098           m21 *=last1;
    1099          
    1100           TMatrix w_time= m11 - m21;
    1101          
    1102           for (Int_t count=0; count < fWindowSizeLoGain; count++)
    1103             {
    1104               const Int_t idx = i+fBinningResolutionLoGain/2+fBinningResolutionLoGain*count-1;
    1105               fAmpWeightsLoGain [idx] = w_amp [0][count];
    1106               fTimeWeightsLoGain[idx] = w_time[0][count];
    1107             }
    1108      
    1109         } // end loop over i
    1110     }
    1111  
    1112   ofstream fn(filename.Data());
    1113 
    1114   fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl;
    1115   fn << "# (Amplitude)  (Time) " << endl;
    1116 
    1117   for (Int_t i=0; i<nsizehi; i++)
    1118     fn << "\t" << fAmpWeightsHiGain[i] << "\t" << fTimeWeightsHiGain[i] << endl;
    1119 
    1120   fn << "# Low Gain Weights: " << fWindowSizeLoGain << " " << fBinningResolutionLoGain << endl;
    1121   fn << "# (Amplitude)  (Time) " << endl;
    1122 
    1123   for (Int_t i=0; i<nsizehi; i++)
    1124     fn << "\t" << fAmpWeightsLoGain[i] << "\t" << fTimeWeightsLoGain[i] << endl;
    1125 
    1126   delete derivativehi;
    1127   if (derivativelo)
    1128     delete derivativelo;
    1129 
    1130   return kTRUE;
    1131 }
    1132 */
     678// Print the setup of the digital filter extraction used. Use
     679//  the option "weights" if you want to print also all weights.
     680//
    1133681void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const
    1134682{
    … …  
    1137685
    1138686    MExtractTimeAndCharge::Print(o);
    1139     *fLog << " Window Size HiGain: " << fWindowSizeHiGain        << "  LoGain: " << fWindowSizeLoGain << endl;
    1140     *fLog << " Binning Res HiGain: " << fBinningResolutionHiGain << "  LoGain: " << fBinningResolutionHiGain << endl;
     687    *fLog << " Window Size HiGain: " << setw(2) << fWindowSizeHiGain        << "  LoGain: " << setw(2) << fWindowSizeLoGain << endl;
     688    *fLog << " Binning Res HiGain: " << setw(2) << fBinningResolutionHiGain << "  LoGain: " << setw(2) << fBinningResolutionHiGain << endl;
    1141689    *fLog << " Weights File desired: " << (fNameWeightsFile.IsNull()?"-":fNameWeightsFile.Data()) << endl;
    1142690    if (!fNameWeightsFileSet.IsNull())
Note: See TracChangeset for help on using the changeset viewer.