EMCAL/AliCaloRawAnalyzerPeakFinder.cxx

   1 /**************************************************************************
   2  * This file is property of and copyright by                              *
   3  * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
   4  *                                                                        *
   5  * Primary Author: Per Thomas Hille  <perthomas.hille@yale.edu>           *
   6  *                                                                        *
   7  * Contributors are mentioned in the code where appropriate.              *
   8  * Please report bugs to p.t.hille@fys.uio.no                             *
   9  *                                                                        *
  10  * Permission to use, copy, modify and distribute this software and its   *
  11  * documentation strictly for non-commercial purposes is hereby granted   *
  12  * without fee, provided that the above copyright notice appears in all   *
  13  * copies and that both the copyright notice and this permission notice   *
  14  * appear in the supporting documentation. The authors make no claims     *
  15  * about the suitability of this software for any purpose. It is          *
  16  * provided "as is" without express or implied warranty.                  *
  17  **************************************************************************/
  18
  19 // The Peak-Finder algorithm
  20 // The amplitude is extracted  as a
  21 // weighted sum of the samples using the
  22 // best possible weights.
  23 // The wights is calculated only once and the
  24 // Actual extraction of amplitude and peak position
  25 // Is done with a simple vector multiplication, allowing for
  26 // Extreemely fast computations.
  27
  28 #include "AliCaloRawAnalyzerPeakFinder.h"
  29 #include "AliCaloBunchInfo.h"
  30 #include "AliCaloFitResults.h"
  31 #include <iostream>
  32 #include "unistd.h"
  33 #include "TMath.h"
  34 #include "AliLog.h"
  35
  36 using namespace std;
  37
  38 ClassImp( AliCaloRawAnalyzerPeakFinder )
  39
  40 AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer("Peak-Finder", "PF")
  41 //    fTof(0),
  42 //                                                            fAmp(0)
  43 {
  44   //comment
  45
  46   fNsampleCut = 5;
  47
  48   for(int i=0; i < MAXSTART; i++)
  49     {
  50       for(int j=0; j < SAMPLERANGE; j++ )
  51         {
  52           fPFAmpVectors[i][j] = new double[100];
  53           fPFTofVectors[i][j] = new double[100];
  54           fPFAmpVectorsCoarse[i][j] = new double[100];
  55           fPFTofVectorsCoarse[i][j] = new double[100];
  56
  57           for(int k=0; k < 100; k++ )
  58             {
  59               fPFAmpVectors[i][j][k] = 0;
  60               fPFTofVectors[i][j][k] = 0;
  61               fPFAmpVectorsCoarse[i][j][k] = 0;
  62               fPFTofVectorsCoarse[i][j][k] = 0;
  63             }
  64         }
  65     }
  66
  67   LoadVectors();
  68
  69 }
  70
  71
  72 AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
  73 {
  74   //comment
  75   for(int i=0; i < MAXSTART; i++)
  76     {
  77       for(int j=0; j < SAMPLERANGE; j++ )
  78         {
  79           delete[] fPFAmpVectors[i][j];
  80           delete[] fPFTofVectors[i][j];
  81           delete[] fPFAmpVectorsCoarse[i][j];
  82           delete[] fPFTofVectorsCoarse[i][j];
  83         }
  84     }
  85 }
  86
  87
  88 Double_t
  89 AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int length ) const
  90 {
  91   Double_t tmpTof = 0;
  92   Double_t tmpAmp= 0;
  93
  94   for(int i=0; i < length; i++)
  95     {
  96       tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i];
  97       tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i];
  98     }
  99
 100   tmpTof = tmpTof / tmpAmp ;
 101   return tmpTof;
 102 }
 103
 104
 105 AliCaloFitResults
 106 AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
 107 {
 108   // Extracting the amplitude using the Peak-Finder algorithm
 109   // The amplitude is a weighted sum of the samples using
 110   // optimum weights.
 111
 112   short maxampindex; //index of maximum amplitude
 113   short maxamp; //Maximum amplitude
 114   //  fAmp = 0;
 115
 116
 117   fAmpA[0] = 0;
 118   fAmpA[1] = 0;
 119   fAmpA[2] = 0;
 120
 121
 122   //  cout << __FILE__ << __LINE__ << "\tendbin = " << bunchvector.at(index).GetEndBin()  <<  "\tstartbin = " << bunchvector.at(index).GetStartBin()  << endl;
 123
 124   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
 125
 126   if( index >= 0)
 127     {
 128       Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
 129       Float_t maxf = TMath::MaxElement(   bunchvector.at(index).GetLength(),  fReversed );
 130       short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
 131       short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1);
 132
 133       if(  maxf < fAmpCut  ||  ( maxamp - ped) > 900  )  // (maxamp - ped) > 900 = Close to saturation (use low gain then)
 134         {
 135           //      cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped  << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex  << endl;
 136           return  AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
 137         }
 138
 139       int first;
 140       int last;
 141
 142       if ( maxf > fAmpCut )
 143         {
 144           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(), maxrev, &first, &last);
 145           int nsamples =  last - first;
 146           if( ( nsamples  )  >= fNsampleCut )
 147             {
 148               int startbin = bunchvector.at(index).GetStartBin();
 149               int n = last - first;
 150               int pfindex = n - fNsampleCut;
 151               pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
 152
 153               int dt =  maxampindex - startbin -2;
 154
 155               //    cout << __FILE__ << __LINE__ <<"\t The coarse estimated t0 is " << ScanCoarse( &fReversed[dt] , n ) << endl;
 156
 157
 158               //     Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
 159
 160               //              cout << __FILE__ << __LINE__ << ",  dt = " << dt << ",\tmaxamindex = " << maxampindex << "\tstartbin = "<< startbin << endl;
 161
 162               for( int i=0; i < SAMPLERANGE; i++ )
 163                 {
 164                   for( int j = 0; j < 3; j++ )
 165                     {
 166                       //    fAmpA[j] += fPFAmpVectors[0][pfindex][i]*tmp[j];
 167                       fAmpA[j] += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i +j -1 ];
 168                     }
 169                 }
 170
 171               double diff = 9999;
 172               int tmpindex = 0;
 173
 174               for(int k=0; k < 3; k ++)
 175                 {
 176                   //              cout << __FILE__ << __LINE__ << "amp[="<< k <<"] = " << fAmpA[k] << endl;
 177                   if(  TMath::Abs(fAmpA[k] - ( maxamp - ped) )  < diff)
 178                     {
 179                       diff = TMath::Abs(fAmpA[k] - ( maxamp - ped));
 180                       tmpindex = k;
 181                     }
 182                 }
 183
 184               Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
 185
 186               if( tmptof < -1 )
 187                 {
 188                   tmpindex = 0;
 189                 }
 190               else
 191                 if( tmptof  > -1 && tmptof < 100 )
 192                   {
 193                     tmpindex =1;
 194                   }
 195                 else
 196                   {
 197                     tmpindex = 2;
 198                   }
 199               double tof = 0;
 200
 201               for(int k=0; k < SAMPLERANGE; k++   )
 202                 {
 203                   tof +=  fPFTofVectors[0][pfindex][k]*fReversed[ dt  +k + tmpindex -1 ];
 204                 }
 205
 206               //              cout << __FILE__ << __LINE__ <<  "tofRaw =   "<< tof /  fAmpA[tmpindex]  << endl;
 207
 208               // tof = tof /  fAmpA[tmpindex] +  (dt + startbin)*100;
 209
 210               if( TMath::Abs(  (maxf - fAmpA[tmpindex])/maxf )  >   0.1 )
 211                 {
 212                   fAmpA[tmpindex] = maxf;
 213                 }
 214
 215               //              timebinOffset
 216
 217                 //      tof = (dt + startbin + tmpindex )*100 - tof/fAmpA[tmpindex];
 218                 // tof = ( timebinOffset )*100 - tof/fAmpA[tmpindex]; // ns
 219               tof = timebinOffset - 0.01*tof/fAmpA[tmpindex]; // clock ticks
 220
 221               //      tof = tof/fAmpA[tmpindex];
 222
 223
 224               return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar, fAmpA[tmpindex], tof, timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy,
 225                                         AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 226             }
 227           else
 228             {
 229               return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kCrude, maxf, timebinOffset);
 230             }
 231         }
 232       else
 233         {
 234           return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
 235         }
 236     }
 237    //  cout << __FILE__ << __LINE__ <<  "WARNING, returning amp = -1 "  <<  endl;
 238   return  AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid);
 239 }
 240
 241
 242 void
 243 AliCaloRawAnalyzerPeakFinder::LoadVectors()
 244 {
 245   //Read in the Peak finder vecors from file
 246   for(int i = 0;  i < MAXSTART ; i++)
 247     {
 248       for( int j=0; j < SAMPLERANGE; j++)
 249         {
 250           char filenameCoarse[256];
 251           char filename[256];
 252
 253           int n = j+fNsampleCut;
 254
 255           //      double start = (double)i+0.5;
 256           double start = (double)i+0;
 257
 258           sprintf(filename,        "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
 259           sprintf(filenameCoarse,  "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
 260
 261           FILE *fp  =  fopen(filename, "r");
 262           FILE *fpc =  fopen(filenameCoarse, "r");
 263
 264           if( fp == 0 )
 265             {
 266               AliFatal( Form( "could not open file: %s", filename ) );
 267             }
 268
 269           if(fpc == 0)
 270             {
 271               AliFatal( Form( "could not open file: %s", filenameCoarse ) );
 272             }
 273           else
 274             {
 275               for(int m = 0; m < n ; m++ )
 276                 {
 277                   // cout << __FILE__ << __LINE__ << "i="<<i <<"\tj=" <<j << "\tm=" << m << endl;
 278
 279                   fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
 280                   //              fPFAmpVectorsCoarse[i][j][m] = 1;
 281                   fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
 282                 }
 283
 284               fscanf(fp,   "\n" );
 285               fscanf(fpc,  "\n" );
 286
 287               for(int m = 0; m < n ; m++ )
 288                 {
 289                   //  fPFTofVectors[i][j][m] = 1;
 290
 291                   fscanf(fp, "%lf\t",   &fPFTofVectors[i][j][m]  );
 292                   fscanf(fpc, "%lf\t",  &fPFTofVectorsCoarse[i][j][m]  );
 293                   //  fPFTofVectorsCoarse[i][j][m] = 1;
 294                 }
 295
 296
 297               fclose (fp);
 298               fclose (fpc);
 299             }
 300         }
 301     }
 302 }
 303
 304
 305
 306 /*
 307 void
 308 AliCaloRawAnalyzerPeakFinder::PolTof( const double rectof ) const
 309 //
 310 {
 311   static Double_t p0 = -55.69;
 312   static Double_t p1 = 3.178;
 313   static Double_t p2 = -0.05587;
 314   static Double_t p3 = 0.0003185;
 315   static Double_t p4 = -7.91E-7;
 316   static Double_t p5 = 7.576E-10;
 317 }
 318 */