+++ /dev/null
-/**************************************************************************
- * This file is property of and copyright by *
- * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 *
- * *
- * Primary Author: Per Thomas Hille <perthomas.hille@yale.edu> *
- * *
- * Contributors are mentioned in the code where appropriate. *
- * Please report bugs to p.t.hille@fys.uio.no *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-// The Peak-Finder algorithm
-// The amplitude is extracted as a
-// weighted sum of the samples using the
-// best possible weights.
-// The wights is calculated only once and the
-// Actual extraction of amplitude and peak position
-// Is done with a simple vector multiplication, allowing for
-// Extreemely fast computations.
-
-#include "AliCaloRawAnalyzerPeakFinderV2.h"
-#include "AliCaloBunchInfo.h"
-#include "AliCaloFitResults.h"
-#include <iostream>
-#include "unistd.h"
-#include "TMath.h"
-#include "AliLog.h"
-
-using namespace std;
-
-ClassImp( AliCaloRawAnalyzerPeakFinderV2 )
-
-AliCaloRawAnalyzerPeakFinderV2::AliCaloRawAnalyzerPeakFinderV2() :AliCaloRawAnalyzer("Peak-FinderV2", "PF")
-// fTof(0),
-// fAmp(0)
-{
- //comment
-
- fNsampleCut = 5;
-
- for(int i=0; i < MAXSTART; i++)
- {
- for(int j=0; j < SAMPLERANGE; j++ )
- {
- fPFAmpVectors[i][j] = new double[100];
- fPFTofVectors[i][j] = new double[100];
- fPFAmpVectorsCoarse[i][j] = new double[100];
- fPFTofVectorsCoarse[i][j] = new double[100];
-
-
- for(int k=0; k < 100; k++ )
- {
- fPFAmpVectors[i][j][k] = 0;
- fPFTofVectors[i][j][k] = 0;
- fPFAmpVectorsCoarse[i][j][k] = 0;
- fPFTofVectorsCoarse[i][j][k] = 0;
- }
- }
- }
-
- LoadVectors();
-
-}
-
-
-AliCaloRawAnalyzerPeakFinderV2::~AliCaloRawAnalyzerPeakFinderV2()
-{
- //comment
- for(int i=0; i < MAXSTART; i++)
- {
- for(int j=0; j < SAMPLERANGE; j++ )
- {
- delete[] fPFAmpVectors[i][j];
- delete[] fPFTofVectors[i][j];
- delete[] fPFAmpVectorsCoarse[i][j];
- delete[] fPFTofVectorsCoarse[i][j];
- }
- }
-}
-
-
-Double_t
-AliCaloRawAnalyzerPeakFinderV2::ScanCoarse(const Double_t *const array, const int length ) const
-{
- Double_t tmpTof = 0;
- Double_t tmpAmp= 0;
-
- for(int i=0; i < length; i++)
- {
- tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i];
- tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i];
- }
-
- tmpTof = tmpTof / tmpAmp ;
- return tmpTof;
-}
-
-
-AliCaloFitResults
-AliCaloRawAnalyzerPeakFinderV2::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
-{
- // Extracting the amplitude using the Peak-FinderV2 algorithm
- // The amplitude is a weighted sum of the samples using
- // optimum weights.
-
- short maxampindex; //index of maximum amplitude
- short maxamp; //Maximum amplitude
- // fAmp = 0;
- fAmpA[0] = 0;
- fAmpA[1] = 0;
- fAmpA[2] = 0;
-
- int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
-
- if( index >= 0)
- {
- Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
- Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
-
- if( maxf < fAmpCut || ( maxamp - ped) > 900 ) // (maxamp - ped) > 900 = Close to saturation (use low gain then)
- {
- // cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex << endl;
- return AliCaloFitResults( maxamp, ped, -1, maxf, maxampindex, -1, -1 );
- }
-
- int first;
- int last;
-
- if ( maxf > fAmpCut )
- {
- SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxampindex - bunchvector.at(index).GetStartBin(), &first, &last);
- int nsamples = last - first;
- if( ( nsamples ) >= fNsampleCut )
- {
- int startbin = bunchvector.at(index).GetStartBin();
- int n = last - first;
- int pfindex = n - fNsampleCut;
- pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
-
- int dt = maxampindex - startbin -1;
-
-
- // cout << __FILE__ << __LINE__ <<"\t The coarse estimated t0 is " << ScanCoarse( &fReversed[dt] , n ) << endl;
-
-
-
- // Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
- // cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl;
-
-
- // cout << __FILE__ << __LINE__ << ", dt = " << dt << ",\tmaxamindex = " << maxampindex << "\tstartbin = "<< startbin << endl;
-
- for( int i=0; i < SAMPLERANGE; i++ )
- {
- for( int j = 0; j < 3; j++ )
- {
- // fAmpA[j] += fPFAmpVectors[0][pfindex][i]*tmp[j];
- fAmpA[j] += fPFAmpVectors[0][pfindex][i]*fReversed[ dt +i +j -1 ];
- }
- }
-
- double diff = 9999;
- int tmpindex = 0;
-
- for(int k=0; k < 3; k ++)
- {
- // cout << __FILE__ << __LINE__ << "amp[="<< k <<"] = " << fAmpA[k] << endl;
- if( TMath::Abs(fAmpA[k] - ( maxamp - ped) ) < diff)
- {
- diff = TMath::Abs(fAmpA[k] - ( maxamp - ped));
- tmpindex = k;
- }
- }
-
- Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
-
- // Double_t tofnew = PolTof(tmptof) + ( dt + startbin )*100 ;
- // int dt = maxampindex - startbin -1;
-
- // Double_t tofnew = PolTof(tmptof) + startbin*100 ;
-
- // Double_t tofnew = PolTof(tmptof) + maxampindex*100 ;
-
-
- // Double_t tofnew = PolTof(tmptof) + (dt + startbin + tmpindex )*100 ;
- Double_t tofnew = PolTof(tmptof) ;
-
- // tmptof= tofnew;
-
- // Double_t tofnew = PolTof(tmptof) + maxampindex ;
-
-
- if(tmptof < 0 )
- {
- cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl;
- }
-
-
- if( tmptof < -1 )
- {
- tmpindex = 0;
- }
- else
- if( tmptof > -1 && tmptof < 100 )
- {
- tmpindex =1;
- }
- else
- {
- tmpindex = 2;
- }
-
-
-
- double tof = 0;
-
- for(int k=0; k < SAMPLERANGE; k++ )
- {
- tof += fPFTofVectors[0][pfindex][k]*fReversed[ dt +k + tmpindex -1 ];
- }
-
- // cout << __FILE__ << __LINE__ << "tofRaw = "<< tof / fAmpA[tmpindex] << endl;
-
- // tof = tof / fAmpA[tmpindex] + (dt + startbin)*100;
-
- if( TMath::Abs( (maxf - fAmpA[tmpindex])/maxf ) > 0.1 )
- {
- fAmpA[tmpindex] = maxf;
- }
-
- // tof = (dt + startbin + tmpindex )*100 - tof/fAmpA[tmpindex]; // ns
- tof = dt + startbin + tmpindex - 0.01*tof/fAmpA[tmpindex]; // clock ticks
-
-
- return AliCaloFitResults( maxamp, ped , -1, fAmpA[tmpindex], tof , -2, -3 );
- }
- else
- {
- return AliCaloFitResults( maxamp, ped , -5, maxf, -6, -7, -8 );
- }
- }
- }
- // cout << __FILE__ << __LINE__ << "WARNING, returning amp = -1 " << endl;
- return AliCaloFitResults(-1, -1);
-}
-
-
-Double_t
-AliCaloRawAnalyzerPeakFinderV2::PolTof( const double fx1 ) const
-{
-
- //Newtons method
- Double_t tolerance = 0.01;
- // cout << "************************" << endl;
- Double_t fx2 = PolValue( fx1 );
- Double_t tmpfx1 = fx1;
-
- while( TMath::Abs( fx2 - fx1 ) > tolerance )
- {
- Double_t der = PolDerivative( tmpfx1 );
- // tmpx = der*( x - val) +x;
- tmpfx1 = ( fx1 - fx2)/der +tmpfx1;
-
- // tmpx = der*( val - tmpx ) +tmpx;
- fx2 = PolValue( tmpfx1 );
- // cout << __FILE__ << __LINE__ << "Der =\t" << der << " x=\t"<<x<<"\tval="<<val << endl;
- // cout << __FILE__ << __LINE__ << "Der =\t" << der << " tmpx=\t"<< tmpfx1 <<"\tval="<< fx2 << endl;
- }
-
- // cout << __FILE__ << __LINE__ << "CONVERGED !! fx1 = "<< fx1 <<" tmpfx1 = "<< tmpfx1 <<" f(tmpfx1) = "<< fx2 << endl;
-
- // cout << "************************" << endl;
-
- return tmpfx1;
-
-}
-
-
-Double_t
-AliCaloRawAnalyzerPeakFinderV2::PolValue(const Double_t x) const
-{
- static Double_t p0 = -55.69;
- static Double_t p1 = 4.718;
- static Double_t p2 = -0.05587;
- static Double_t p3 = 0.0003185;
- static Double_t p4 = -7.91E-7;
- static Double_t p5 = 7.576E-10;
-
- // return p0 + p1*x + p2*TMath::Power(x, 2) + p3*TMath::Power(x, 3) + p4*TMath::Power(x, 4) + p5*TMath::Power(x, 5);
-
- return p0 + p1*x + p2*x*x + p3*x*x*x + p4*x*x*x*x + p5*x*x*x*x*x;
-
-}
-
-
-Double_t
-AliCaloRawAnalyzerPeakFinderV2::PolDerivative(const Double_t x) const
-{
- static Double_t dp0 = 0;
- static Double_t dp1 = 4.718;
- static Double_t dp2 = -0.11174;
- static Double_t dp3 = 0.0009555;
- static Double_t dp4 = -3.164E-6;
- static Double_t dp5 = 3.788E-9;
-
- // return dp0 + dp1 + dp2*x + dp3*TMath::Power(x, 2) + dp4*TMath::Power(x, 3) + dp5*TMath::Power(x, 4);
-
-
-
- return dp0 + dp1 + dp2*x + dp3*x*x + dp4*x*x*x + dp5*x*x*x*x;
-
-}
-
-
-void
-AliCaloRawAnalyzerPeakFinderV2::LoadVectors()
-{
- //Read in the Peak finder vecors from file
- for(int i = 0; i < MAXSTART ; i++)
- {
- for( int j=0; j < SAMPLERANGE; j++)
- {
- char filenameCoarse[256];
- char filename[256];
-
- int n = j+fNsampleCut;
-
- // double start = (double)i+0.5;
- double start = (double)i+0;
-
- sprintf(filename, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
- sprintf(filenameCoarse, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
-
- FILE *fp = fopen(filename, "r");
- FILE *fpc = fopen(filenameCoarse, "r");
-
- if( fp == 0 )
- {
- AliFatal( Form( "could not open file: %s", filename ) );
- }
-
- if(fpc == 0)
- {
- AliFatal( Form( "could not open file: %s", filenameCoarse ) );
- }
- else
- {
- for(int m = 0; m < n ; m++ )
- {
- cout << __FILE__ << __LINE__ << "i="<<i <<"\tj=" <<j << "\tm=" << m << endl;
-
- fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
- // fPFAmpVectorsCoarse[i][j][m] = 1;
- fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
- }
-
- fscanf(fp, "\n" );
- fscanf(fpc, "\n" );
-
- for(int m = 0; m < n ; m++ )
- {
- // fPFTofVectors[i][j][m] = 1;
-
- fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] );
- fscanf(fpc, "%lf\t", &fPFTofVectorsCoarse[i][j][m] );
- // fPFTofVectorsCoarse[i][j][m] = 1;
- }
- fclose (fp);
- fclose (fpc);
- }
- }
- }
-}
-
-
+++ /dev/null
-#ifndef ALICALORAWANALYZERPEAKFINDERV2_H
-#define ALICALORAWANALYZERPEAKFINDERV2_H
-
-
-/**************************************************************************
- * This file is property of and copyright by *
- * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 *
- * *
- * Primary Author: Per Thomas Hille <perthomas.hille@yale.edu> *
- * *
- * Contributors are mentioned in the code where appropriate. *
- * Please report bugs to p.t.hille@fys.uio.no *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-
-// The Peak-Finder algorithm
-// The amplitude is extracted as a
-// weighted sum of the samples using the
-// best possible weights.
-
-
-#include "AliCaloRawAnalyzer.h"
-
-#define MAXSTART 3
-#define SAMPLERANGE 15
-#define SHIF 0.5
-
-class AliCaloBunchInfo;
-
-
-class AliCaloRawAnalyzerPeakFinderV2 : public AliCaloRawAnalyzer
-{
- public:
- AliCaloRawAnalyzerPeakFinderV2();
- virtual ~AliCaloRawAnalyzerPeakFinderV2();
- virtual AliCaloFitResults Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 );
-
- private:
- AliCaloRawAnalyzerPeakFinderV2( const AliCaloRawAnalyzerPeakFinderV2 & );
- AliCaloRawAnalyzerPeakFinderV2 & operator = ( const AliCaloRawAnalyzerPeakFinderV2 & );
-
- void LoadVectors();
- Double_t ScanCoarse(const Double_t *const array, const int length ) const ; // Find a rough estimate of peak position and t0
-
-
- Double_t PolTof(const Double_t x) const;
- Double_t PolValue(const Double_t x) const;
- Double_t PolDerivative(const Double_t x) const;
-
-
-
- double *fPFAmpVectorsCoarse[MAXSTART][SAMPLERANGE]; // Vectors for Amplitude extraction, first iteration
- double *fPFTofVectorsCoarse[MAXSTART][SAMPLERANGE]; // Vectors for TOF extraction, first iteration
-
- double *fPFAmpVectors[MAXSTART][SAMPLERANGE]; // Vectors for Amplitude extraction, second iteration
- double *fPFTofVectors[MAXSTART][SAMPLERANGE]; // Vectors for TOF extraction, second iteration
-
- // double fTof;
- // double fAmp;
-
- double fAmpA[3]; // The amplitude of the signal (eveluate 3 times using 3 differtnt phase shifts of the input samples )
- // double fAmp2;
- // double fAmp3;
-
- ClassDef( AliCaloRawAnalyzerPeakFinderV2, 1 )
-
-};
-
-#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff