EMCAL/AliCaloRawAnalyzerPeakFinder.cxx

   1 /**************************************************************************
   2  * This file is property of and copyright by the Experimental Nuclear     *
   3  * Physics Group, Dep. of Physics                                         *
   4  * University of Oslo, Norway, 2007                                       *
   5  *                                                                        *
   6  * Author: Per Thomas Hille <perthi@fys.uio.no> for the ALICE HLT Project.*
   7  * Contributors are mentioned in the code where appropriate.              *
   8  * Please report bugs to perthi@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 AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer(),
  39                                                                 fTof(0),
  40                                                                 fAmp(0)
  41 {
  42   //comment
  43
  44   fNsampleCut = 5;
  45
  46   for(int i=0; i < MAXSTART; i++)
  47     {
  48       for(int j=0; j < SAMPLERANGE; j++ )
  49         {
  50           fPFAmpVectors[i][j] = new double[100];
  51           fPFTofVectors[i][j] = new double[100];
  52
  53           for(int k=0; k < 100; k++ )
  54             {
  55               fPFAmpVectors[i][j][k] = 0;
  56               fPFTofVectors[i][j][k] = 0;
  57             }
  58         }
  59     }
  60   LoadVectors();
  61 }
  62
  63
  64 AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
  65 {
  66   //comment
  67   for(int i=0; i < MAXSTART; i++)
  68     {
  69       for(int j=0; j < SAMPLERANGE; j++ )
  70         {
  71           delete[] fPFAmpVectors[i][j];
  72           delete[] fPFTofVectors[i][j];
  73         }
  74     }
  75 }
  76
  77
  78 AliCaloFitResults
  79 AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
  80 {
  81   // Extracting the amplitude using the Peak-Finder algorithm
  82   // The amplitude is a weighted sum of the samples using
  83   // optimum weights.
  84
  85   short maxampindex; //index of maximum amplitude
  86   short maxamp; //Maximum amplitude
  87   fAmp = 0;
  88   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
  89
  90   if( index >= 0)
  91     {
  92       Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
  93       Float_t maxf = TMath::MaxElement(   bunchvector.at(index).GetLength(),  fReversed );
  94
  95       if(  maxf < fAmpCut  ||  ( maxamp - ped) > 900  )  // (maxamp - ped) > 900 = Close to saturation (use low gain then)
  96         {
  97           //      cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped  << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex  << endl;
  98           return  AliCaloFitResults( maxamp, ped,  -1, maxf,   maxampindex, -1, -1 );
  99         }
 100
 101       int first;
 102       int last;
 103
 104       if ( maxf > fAmpCut )
 105         {
 106           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxampindex -  bunchvector.at(index).GetStartBin(), &first, &last);
 107           int nsamples =  last - first;
 108           if( ( nsamples  )  >= fNsampleCut )
 109             {
 110               int startbin = bunchvector.at(index).GetStartBin();
 111               int n = last -first;
 112               int pfindex = n - fNsampleCut;
 113               pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
 114
 115               for(int i=0; i < SAMPLERANGE; i++ )
 116                 {
 117                   int dt =  maxampindex - startbin -2;
 118                   double tmp =  fReversed[ dt  +i];
 119                   fAmp += fPFAmpVectors[0][pfindex][i]*tmp;
 120                 }
 121
 122               return AliCaloFitResults( maxamp, ped , -1, fAmp, -1, -1, -1 );
 123             }
 124           else
 125             {
 126               return AliCaloFitResults( maxamp, ped , -1, maxf, -1, -1, -1 );
 127             }
 128         }
 129     }
 130
 131   //  cout << __FILE__ << __LINE__ <<  "WARNING, returning amp = -1 "  <<  endl;
 132
 133   return  AliCaloFitResults(-1, -1);
 134 }
 135
 136
 137 void
 138 AliCaloRawAnalyzerPeakFinder::LoadVectors()
 139 {
 140   //Read in the Peak finder vecors from file
 141   for(int i = 0;  i < MAXSTART ; i++)
 142     {
 143       for( int j=0; j < SAMPLERANGE; j++)
 144         {
 145           char filename[256];
 146           int n = j+fNsampleCut;
 147           double start = (double)i+0.5;
 148           sprintf(filename, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
 149           FILE *fp = fopen(filename, "r");
 150
 151           if(fp == 0)
 152             {
 153               AliFatal( Form( "could not open file: %s", filename ) );
 154             }
 155           else
 156             {
 157               for(int m = 0; m < n ; m++ )
 158                 {
 159                   fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
 160                 }
 161               fscanf(fp, "\n");
 162
 163               for(int m = 0; m < n ; m++ )
 164                 {
 165                   fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] );
 166                 }
 167
 168               fclose (fp);
 169             }
 170         }
 171     }
 172 }