EMCAL/AliCaloRawAnalyzerPeakFinderV2.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 "AliCaloRawAnalyzerPeakFinderV2.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( AliCaloRawAnalyzerPeakFinderV2 )
  39
  40 AliCaloRawAnalyzerPeakFinderV2::AliCaloRawAnalyzerPeakFinderV2() :AliCaloRawAnalyzer("Peak-FinderV2", "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
  58           for(int k=0; k < 100; k++ )
  59             {
  60               fPFAmpVectors[i][j][k] = 0;
  61               fPFTofVectors[i][j][k] = 0;
  62               fPFAmpVectorsCoarse[i][j][k] = 0;
  63               fPFTofVectorsCoarse[i][j][k] = 0;
  64             }
  65         }
  66     }
  67
  68   LoadVectors();
  69
  70 }
  71
  72
  73 AliCaloRawAnalyzerPeakFinderV2::~AliCaloRawAnalyzerPeakFinderV2()
  74 {
  75   //comment
  76   for(int i=0; i < MAXSTART; i++)
  77     {
  78       for(int j=0; j < SAMPLERANGE; j++ )
  79         {
  80           delete[] fPFAmpVectors[i][j];
  81           delete[] fPFTofVectors[i][j];
  82           delete[] fPFAmpVectorsCoarse[i][j];
  83           delete[] fPFTofVectorsCoarse[i][j];
  84         }
  85     }
  86 }
  87
  88
  89 Double_t
  90 AliCaloRawAnalyzerPeakFinderV2::ScanCoarse(const Double_t *const array, const int length ) const
  91 {
  92   Double_t tmpTof = 0;
  93   Double_t tmpAmp= 0;
  94
  95   for(int i=0; i < length; i++)
  96     {
  97       tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i];
  98       tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i];
  99     }
 100
 101   tmpTof = tmpTof / tmpAmp ;
 102   return tmpTof;
 103 }
 104
 105
 106 AliCaloFitResults
 107 AliCaloRawAnalyzerPeakFinderV2::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
 108 {
 109   // Extracting the amplitude using the Peak-FinderV2 algorithm
 110   // The amplitude is a weighted sum of the samples using
 111   // optimum weights.
 112
 113   short maxampindex; //index of maximum amplitude
 114   short maxamp; //Maximum amplitude
 115   //  fAmp = 0;
 116   fAmpA[0] = 0;
 117   fAmpA[1] = 0;
 118   fAmpA[2] = 0;
 119
 120   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
 121
 122   if( index >= 0)
 123     {
 124       Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
 125       Float_t maxf = TMath::MaxElement(   bunchvector.at(index).GetLength(),  fReversed );
 126
 127       if(  maxf < fAmpCut  ||  ( maxamp - ped) > 900  )  // (maxamp - ped) > 900 = Close to saturation (use low gain then)
 128         {
 129           //      cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped  << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex  << endl;
 130           return  AliCaloFitResults( maxamp, ped,  -1, maxf,   maxampindex, -1, -1 );
 131         }
 132
 133       int first;
 134       int last;
 135
 136       if ( maxf > fAmpCut )
 137         {
 138           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxampindex -  bunchvector.at(index).GetStartBin(), &first, &last);
 139           int nsamples =  last - first;
 140           if( ( nsamples  )  >= fNsampleCut )
 141             {
 142               int startbin = bunchvector.at(index).GetStartBin();
 143               int n = last - first;
 144               int pfindex = n - fNsampleCut;
 145               pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
 146
 147               int dt =  maxampindex - startbin -1;
 148
 149
 150               //    cout << __FILE__ << __LINE__ <<"\t The coarse estimated t0 is " << ScanCoarse( &fReversed[dt] , n ) << endl;
 151
 152
 153
 154               //  Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
 155               //  cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl;
 156
 157
 158               //              cout << __FILE__ << __LINE__ << ",  dt = " << dt << ",\tmaxamindex = " << maxampindex << "\tstartbin = "<< startbin << endl;
 159
 160               for( int i=0; i < SAMPLERANGE; i++ )
 161                 {
 162                   for( int j = 0; j < 3; j++ )
 163                     {
 164                       //    fAmpA[j] += fPFAmpVectors[0][pfindex][i]*tmp[j];
 165                       fAmpA[j] += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i +j -1 ];
 166                     }
 167                 }
 168
 169               double diff = 9999;
 170               int tmpindex = 0;
 171
 172               for(int k=0; k < 3; k ++)
 173                 {
 174                   //              cout << __FILE__ << __LINE__ << "amp[="<< k <<"] = " << fAmpA[k] << endl;
 175                   if(  TMath::Abs(fAmpA[k] - ( maxamp - ped) )  < diff)
 176                     {
 177                       diff = TMath::Abs(fAmpA[k] - ( maxamp - ped));
 178                       tmpindex = k;
 179                     }
 180                 }
 181
 182               Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
 183
 184               //      Double_t tofnew = PolTof(tmptof) + ( dt + startbin  )*100  ;
 185               //             int dt =  maxampindex - startbin -1;
 186
 187               //              Double_t tofnew = PolTof(tmptof) + startbin*100 ;
 188
 189               //      Double_t tofnew = PolTof(tmptof)  +  maxampindex*100 ;
 190
 191
 192               //              Double_t tofnew = PolTof(tmptof)  + (dt + startbin + tmpindex )*100 ;
 193               Double_t tofnew = PolTof(tmptof)  ;
 194
 195               //              tmptof= tofnew;
 196
 197               //              Double_t tofnew = PolTof(tmptof) + maxampindex  ;
 198
 199
 200               if(tmptof < 0 )
 201                 {
 202                   cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl;
 203                 }
 204
 205
 206               if( tmptof < -1 )
 207                 {
 208                   tmpindex = 0;
 209                 }
 210               else
 211                 if( tmptof  > -1 && tmptof < 100 )
 212                   {
 213                     tmpindex =1;
 214                   }
 215                 else
 216                   {
 217                     tmpindex = 2;
 218                   }
 219
 220
 221
 222               double tof = 0;
 223
 224               for(int k=0; k < SAMPLERANGE; k++   )
 225                 {
 226                   tof +=  fPFTofVectors[0][pfindex][k]*fReversed[ dt  +k + tmpindex -1 ];
 227                 }
 228
 229               //              cout << __FILE__ << __LINE__ <<  "tofRaw =   "<< tof /  fAmpA[tmpindex]  << endl;
 230
 231               // tof = tof /  fAmpA[tmpindex] +  (dt + startbin)*100;
 232
 233               if( TMath::Abs(  (maxf - fAmpA[tmpindex])/maxf )  >   0.1 )
 234                 {
 235                   fAmpA[tmpindex] = maxf;
 236                 }
 237
 238               // tof = (dt + startbin + tmpindex )*100 - tof/fAmpA[tmpindex]; // ns
 239               tof = dt + startbin + tmpindex - 0.01*tof/fAmpA[tmpindex]; // clock ticks
 240
 241
 242               return AliCaloFitResults( maxamp, ped , -1, fAmpA[tmpindex], tof  , -2, -3 );
 243             }
 244           else
 245             {
 246               return AliCaloFitResults( maxamp, ped , -5, maxf, -6, -7, -8 );
 247             }
 248         }
 249     }
 250    //  cout << __FILE__ << __LINE__ <<  "WARNING, returning amp = -1 "  <<  endl;
 251   return  AliCaloFitResults(-1, -1);
 252 }
 253
 254
 255 Double_t
 256 AliCaloRawAnalyzerPeakFinderV2::PolTof( const double fx1 ) const
 257 {
 258
 259   //Newtons method
 260   Double_t tolerance = 0.01;
 261   //  cout << "************************"  << endl;
 262   Double_t fx2 = PolValue( fx1 );
 263   Double_t tmpfx1 = fx1;
 264
 265   while(  TMath::Abs( fx2 - fx1  ) >  tolerance )
 266     {
 267       Double_t der = PolDerivative( tmpfx1 );
 268       //     tmpx = der*( x - val) +x;
 269       tmpfx1 = ( fx1 - fx2)/der  +tmpfx1;
 270
 271       //    tmpx = der*( val - tmpx ) +tmpx;
 272       fx2 = PolValue(  tmpfx1 );
 273       //     cout << __FILE__ << __LINE__ <<  "Der =\t" << der << "  x=\t"<<x<<"\tval="<<val <<  endl;
 274       //     cout << __FILE__ << __LINE__ <<  "Der =\t" << der << "  tmpx=\t"<< tmpfx1 <<"\tval="<< fx2 <<  endl;
 275     }
 276
 277   //  cout << __FILE__ << __LINE__ << "CONVERGED !! fx1 = "<< fx1 <<"  tmpfx1 = "<< tmpfx1  <<"  f(tmpfx1) = "<< fx2 << endl;
 278
 279   //  cout << "************************"  << endl;
 280
 281   return tmpfx1;
 282
 283 }
 284
 285
 286 Double_t
 287 AliCaloRawAnalyzerPeakFinderV2::PolValue(const Double_t x) const
 288 {
 289   static Double_t p0 = -55.69;
 290   static Double_t p1 = 4.718;
 291   static Double_t p2 = -0.05587;
 292   static Double_t p3 = 0.0003185;
 293   static Double_t p4 = -7.91E-7;
 294   static Double_t p5 = 7.576E-10;
 295
 296   //  return  p0 + p1*x + p2*TMath::Power(x, 2) + p3*TMath::Power(x, 3) + p4*TMath::Power(x, 4)  + p5*TMath::Power(x, 5);
 297
 298   return  p0 + p1*x + p2*x*x + p3*x*x*x + p4*x*x*x*x  + p5*x*x*x*x*x;
 299
 300 }
 301
 302
 303 Double_t
 304 AliCaloRawAnalyzerPeakFinderV2::PolDerivative(const Double_t x) const
 305 {
 306   static Double_t dp0 = 0;
 307   static Double_t dp1 = 4.718;
 308   static Double_t dp2 = -0.11174;
 309   static Double_t dp3 = 0.0009555;
 310   static Double_t dp4 = -3.164E-6;
 311   static Double_t dp5 = 3.788E-9;
 312
 313   //  return  dp0 + dp1 + dp2*x + dp3*TMath::Power(x, 2) + dp4*TMath::Power(x, 3)  + dp5*TMath::Power(x, 4);
 314
 315
 316
 317   return  dp0 + dp1 + dp2*x + dp3*x*x + dp4*x*x*x  + dp5*x*x*x*x;
 318
 319 }
 320
 321
 322 void
 323 AliCaloRawAnalyzerPeakFinderV2::LoadVectors()
 324 {
 325   //Read in the Peak finder vecors from file
 326   for(int i = 0;  i < MAXSTART ; i++)
 327     {
 328       for( int j=0; j < SAMPLERANGE; j++)
 329         {
 330           char filenameCoarse[256];
 331           char filename[256];
 332
 333           int n = j+fNsampleCut;
 334
 335           //      double start = (double)i+0.5;
 336           double start = (double)i+0;
 337
 338           sprintf(filename,        "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
 339           sprintf(filenameCoarse,  "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
 340
 341           FILE *fp  =  fopen(filename, "r");
 342           FILE *fpc =  fopen(filenameCoarse, "r");
 343
 344           if( fp == 0 )
 345             {
 346               AliFatal( Form( "could not open file: %s", filename ) );
 347             }
 348
 349           if(fpc == 0)
 350             {
 351               AliFatal( Form( "could not open file: %s", filenameCoarse ) );
 352             }
 353           else
 354             {
 355               for(int m = 0; m < n ; m++ )
 356                 {
 357                   cout << __FILE__ << __LINE__ << "i="<<i <<"\tj=" <<j << "\tm=" << m << endl;
 358
 359                   fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
 360                   //              fPFAmpVectorsCoarse[i][j][m] = 1;
 361                   fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
 362                 }
 363
 364               fscanf(fp,   "\n" );
 365               fscanf(fpc,  "\n" );
 366
 367               for(int m = 0; m < n ; m++ )
 368                 {
 369                   //  fPFTofVectors[i][j][m] = 1;
 370
 371                   fscanf(fp, "%lf\t",   &fPFTofVectors[i][j][m]  );
 372                   fscanf(fpc, "%lf\t",  &fPFTofVectorsCoarse[i][j][m]  );
 373                   //  fPFTofVectorsCoarse[i][j][m] = 1;
 374                 }
 375               fclose (fp);
 376               fclose (fpc);
 377             }
 378         }
 379     }
 380 }
 381
 382