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")
  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
  55           for(int k=0; k < 100; k++ )
  56             {
  57               fPFAmpVectors[i][j][k] = 0;
  58               fPFTofVectors[i][j][k] = 0;
  59             }
  60         }
  61     }
  62   LoadVectors();
  63 }
  64
  65
  66 AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
  67 {
  68   //comment
  69   for(int i=0; i < MAXSTART; i++)
  70     {
  71       for(int j=0; j < SAMPLERANGE; j++ )
  72         {
  73           delete[] fPFAmpVectors[i][j];
  74           delete[] fPFTofVectors[i][j];
  75         }
  76     }
  77 }
  78
  79
  80 AliCaloFitResults
  81 AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
  82 {
  83   // Extracting the amplitude using the Peak-Finder algorithm
  84   // The amplitude is a weighted sum of the samples using
  85   // optimum weights.
  86
  87   short maxampindex; //index of maximum amplitude
  88   short maxamp; //Maximum amplitude
  89   //  fAmp = 0;
  90   fAmpA[0] = 0;
  91   fAmpA[1] = 0;
  92   fAmpA[2] = 0;
  93
  94   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
  95
  96   if( index >= 0)
  97     {
  98       Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
  99       Float_t maxf = TMath::MaxElement(   bunchvector.at(index).GetLength(),  fReversed );
 100
 101       if(  maxf < fAmpCut  ||  ( maxamp - ped) > 900  )  // (maxamp - ped) > 900 = Close to saturation (use low gain then)
 102         {
 103           //      cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped  << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex  << endl;
 104           return  AliCaloFitResults( maxamp, ped,  -1, maxf,   maxampindex, -1, -1 );
 105         }
 106
 107       int first;
 108       int last;
 109
 110       if ( maxf > fAmpCut )
 111         {
 112           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxampindex -  bunchvector.at(index).GetStartBin(), &first, &last);
 113           int nsamples =  last - first;
 114           if( ( nsamples  )  >= fNsampleCut )
 115             {
 116               int startbin = bunchvector.at(index).GetStartBin();
 117               int n = last -first;
 118               int pfindex = n - fNsampleCut;
 119               pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
 120
 121               int dt =  maxampindex - startbin -2;
 122               for(int i=0; i < SAMPLERANGE; i++ )
 123                 {
 124
 125                   //      int dt =  maxampindex - startbin -2;
 126                   //              double tmp[3];
 127                   //      tmp[0]  =  fReversed[ dt  +i -1];
 128                   //              tmp[1]  =  fReversed[ dt  +i];
 129                   //              tmp[2]  =  fReversed[ dt  +i +1];
 130
 131                   for(int j = 0; j < 3; j++)
 132                     {
 133                       //    fAmpA[j] += fPFAmpVectors[0][pfindex][i]*tmp[j];
 134                       fAmpA[j] += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i +j -1];
 135
 136                     }
 137                 }
 138
 139               double diff = 9999;
 140               int tmpindex = 0;
 141
 142               for(int k=0; k < 3; k ++)
 143                 {
 144                   //              cout << __FILE__ << __LINE__ << "amp[="<< k <<"] = " << fAmpA[k] << endl;
 145                   if(  TMath::Abs(fAmpA[k] - ( maxamp - ped) )  < diff)
 146                     {
 147                       diff = TMath::Abs(fAmpA[k] - ( maxamp - ped));
 148                       tmpindex = k;
 149                     }
 150                 }
 151
 152               double tof = 0;
 153
 154               for(int k=0; k < SAMPLERANGE; k++   )
 155                 {
 156                   tof +=  fPFTofVectors[0][pfindex][k]*fReversed[ dt  +k + tmpindex -1 ];
 157                 }
 158
 159               tof = tof /  fAmpA[tmpindex];
 160
 161               return AliCaloFitResults( maxamp, ped , -1, fAmpA[tmpindex], tof, -2, -3 );
 162             }
 163
 164           else
 165             {
 166               return AliCaloFitResults( maxamp, ped , -5, maxf, -6, -7, -8 );
 167             }
 168         }
 169     }
 170
 171   //  cout << __FILE__ << __LINE__ <<  "WARNING, returning amp = -1 "  <<  endl;
 172
 173   return  AliCaloFitResults(-1, -1);
 174 }
 175
 176
 177 void
 178 AliCaloRawAnalyzerPeakFinder::LoadVectors()
 179 {
 180   //Read in the Peak finder vecors from file
 181   for(int i = 0;  i < MAXSTART ; i++)
 182     {
 183       for( int j=0; j < SAMPLERANGE; j++)
 184         {
 185           char filename[256];
 186           int n = j+fNsampleCut;
 187           double start = (double)i+0.5;
 188           sprintf(filename, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
 189           FILE *fp = fopen(filename, "r");
 190
 191           if(fp == 0)
 192             {
 193               AliFatal( Form( "could not open file: %s", filename ) );
 194             }
 195           else
 196             {
 197               for(int m = 0; m < n ; m++ )
 198                 {
 199                   fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
 200                 }
 201               fscanf(fp, "\n");
 202
 203               for(int m = 0; m < n ; m++ )
 204                 {
 205                   fPFTofVectors[i][j][m] = 1;
 206          //       fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] );
 207                 }
 208
 209               fclose (fp);
 210             }
 211         }
 212     }
 213 }