From d655d7dd7062c413d1837b7fd20121d1632d2dba Mon Sep 17 00:00:00 2001 From: dsilverm Date: Fri, 22 Jan 2010 02:51:04 +0000 Subject: [PATCH] peak finder code from Per Thomas/Yale --- EMCAL/AliEMCALBunchInfo.cxx | 67 +++++++ EMCAL/AliEMCALBunchInfo.h | 54 ++++++ EMCAL/AliEMCALFitResults.cxx | 67 +++++++ EMCAL/AliEMCALFitResults.h | 66 +++++++ EMCAL/AliEMCALRawAnalyzer.cxx | 240 ++++++++++++++++++++++++ EMCAL/AliEMCALRawAnalyzer.h | 74 ++++++++ EMCAL/AliEMCALRawAnalyzerLMS.cxx | 152 +++++++++++++++ EMCAL/AliEMCALRawAnalyzerLMS.h | 53 ++++++ EMCAL/AliEMCALRawAnalyzerPeakFinder.cxx | 172 +++++++++++++++++ EMCAL/AliEMCALRawAnalyzerPeakFinder.h | 55 ++++++ EMCAL/libEMCALrec.pkg | 7 +- 11 files changed, 1006 insertions(+), 1 deletion(-) create mode 100644 EMCAL/AliEMCALBunchInfo.cxx create mode 100644 EMCAL/AliEMCALBunchInfo.h create mode 100644 EMCAL/AliEMCALFitResults.cxx create mode 100644 EMCAL/AliEMCALFitResults.h create mode 100644 EMCAL/AliEMCALRawAnalyzer.cxx create mode 100644 EMCAL/AliEMCALRawAnalyzer.h create mode 100644 EMCAL/AliEMCALRawAnalyzerLMS.cxx create mode 100644 EMCAL/AliEMCALRawAnalyzerLMS.h create mode 100644 EMCAL/AliEMCALRawAnalyzerPeakFinder.cxx create mode 100644 EMCAL/AliEMCALRawAnalyzerPeakFinder.h diff --git a/EMCAL/AliEMCALBunchInfo.cxx b/EMCAL/AliEMCALBunchInfo.cxx new file mode 100644 index 00000000000..585b4dd21f0 --- /dev/null +++ b/EMCAL/AliEMCALBunchInfo.cxx @@ -0,0 +1,67 @@ +/************************************************************************** + * This file is property of and copyright by * + * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 * + * * + * Primary Author: Per Thomas Hille * + * * + * 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. * + **************************************************************************/ + + +#include "AliEMCALBunchInfo.h" + +// Container class to hold information +// about of bunches. +// Used by the +// AliEMCALRawAnalyzer +// classes +AliEMCALBunchInfo::AliEMCALBunchInfo( UInt_t start, Int_t length, const UShort_t * data ) : fStartTimebin(start), + fLength(length), + fkData(data) +{ + + +} + + + +AliEMCALBunchInfo::~AliEMCALBunchInfo() +{ + + +} + + + +AliEMCALBunchInfo::AliEMCALBunchInfo( const AliEMCALBunchInfo & rhs) :fStartTimebin( rhs.fStartTimebin ), + fLength( rhs.fLength ), + fkData( rhs.fkData ) +{ + + +} + + + + +AliEMCALBunchInfo& AliEMCALBunchInfo::operator = ( const AliEMCALBunchInfo & rhs) +{ + //This is just to get of compliation warning. Its not really needed + if(this != & rhs) + { + fStartTimebin = rhs.fStartTimebin; + fLength = rhs.fLength; + fkData = rhs.fkData; + } + return *this; +} + diff --git a/EMCAL/AliEMCALBunchInfo.h b/EMCAL/AliEMCALBunchInfo.h new file mode 100644 index 00000000000..44e44061a6d --- /dev/null +++ b/EMCAL/AliEMCALBunchInfo.h @@ -0,0 +1,54 @@ +#ifndef ALIEMCALBUNCHINFO_H +#define ALIEMCALBUNCHINFO_H + +/************************************************************************** + * This file is property of and copyright by * + * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 * + * * + * Primary Author: Per Thomas Hille * + * * + * 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. * + **************************************************************************/ + +#include "Rtypes.h" + +// Container class to hold +// information about ALTRO +// Bunces from the altro stream. +// Each bunch has a start marker, ( fStartTimebin ) +// the number of ADC samples in the bunch fLength, and a pointer +// to the last (fStartTimebin + fLength ) time bin of the bunch. +// +class AliEMCALBunchInfo +{ + public: + AliEMCALBunchInfo( UInt_t starttimebin, Int_t length, const UShort_t * data ); + virtual ~AliEMCALBunchInfo(); + + AliEMCALBunchInfo( const AliEMCALBunchInfo & rhs); + AliEMCALBunchInfo & operator = ( const AliEMCALBunchInfo & rhs); + + + UInt_t GetStartBin( ) const { return fStartTimebin;}; + Int_t GetLength() const { return fLength; }; + const UShort_t *GetData() const { return fkData; }; + + private: + AliEMCALBunchInfo(); + UInt_t fStartTimebin; //Starttimebin as given by the ALTRO stream + Int_t fLength; //Length of the bunch + const UShort_t *fkData; //Pointer to the last data enetry of the bunch (data is reversed with respect to fStartTimebin) +}; + + + +#endif diff --git a/EMCAL/AliEMCALFitResults.cxx b/EMCAL/AliEMCALFitResults.cxx new file mode 100644 index 00000000000..055ee6ef5aa --- /dev/null +++ b/EMCAL/AliEMCALFitResults.cxx @@ -0,0 +1,67 @@ +/************************************************************************** + * This file is property of and copyright by * + * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 * + * * + * Primary Author: Per Thomas Hille * + * * + * 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. * + **************************************************************************/ + + +#include "AliEMCALFitResults.h" + + +// Container class to hold results from fitting +// as well as other methods for +// raw data signals extraction. The class memebers +// fChi2Sig, fNdfSig is only relevant if a fitting procedure is +// Applied. fStatus holds information on wether or not +// The signal was fitted sucessfully. fStatus might have a different meaning If other +// procedures than A different meaning Fitting is applied +AliEMCALFitResults::AliEMCALFitResults(const UShort_t maxSig, const Float_t ped, + const Short_t fitstatus, const Float_t amp, + const Float_t t0, const Float_t chi, + const UShort_t ndf, UShort_t minSig ) : fMaxSig(maxSig), + fPed(ped), + fStatus(fitstatus), + fAmpSig(amp), + fT0(t0), + fChi2Sig(chi), + fNdfSig(ndf), + fMinSig(minSig) +{ + + +} + + + + +AliEMCALFitResults::AliEMCALFitResults(const UShort_t maxSig, const UShort_t minSig) : fMaxSig(maxSig), + fPed(-98), + fStatus( -1 ), + fAmpSig( -1 ), + fT0(-99), + fChi2Sig( -1 ), + fNdfSig( -1), + fMinSig (minSig) +{ + +} + + + +AliEMCALFitResults::~AliEMCALFitResults() +{ + +} + diff --git a/EMCAL/AliEMCALFitResults.h b/EMCAL/AliEMCALFitResults.h new file mode 100644 index 00000000000..799b66edcee --- /dev/null +++ b/EMCAL/AliEMCALFitResults.h @@ -0,0 +1,66 @@ +// -*- mode: c++ -*- +#ifndef ALIEMCALFITRESULTS_H +#define ALIEMCALFITRESULTS_H +/************************************************************************** + * This file is property of and copyright by * + * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 * + * * + * Primary Author: Per Thomas Hille * + * * + * 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. * + **************************************************************************/ + + +#include "Rtypes.h" + +// Container class to hold results from fitting +// as well as other methods for +// raw data signals extraction +class AliEMCALFitResults +{ + public: + explicit AliEMCALFitResults( const UShort_t maxSig, + const Float_t ped, + const Short_t fitStatus, + const Float_t amp, + const Float_t t0, + const Float_t chi, + const UShort_t ndf, + const UShort_t minSig = -99); + + explicit AliEMCALFitResults( const UShort_t maxSig, const UShort_t minSig ); + //AliEMCALFitResults( const UShort_t maxSig, const UShort_t minSig ); + + + virtual ~AliEMCALFitResults(); + UShort_t GetMaxSig() const { return fMaxSig;}; + Float_t GetPed() const { return fPed;}; + UShort_t GetMinSig() const { return fMinSig;}; + UShort_t GetStatus() const { return fStatus;}; + Float_t GetAmp() const { return fAmpSig; }; + Float_t GetTof() const { return fT0; }; + Float_t GetChisSquare() const { return fChi2Sig;}; + UShort_t GetNdf() const { return fNdfSig; }; + + private: + AliEMCALFitResults(); + UShort_t fMaxSig; //Maximum sample value ( 0 - 1023 ) + Float_t fPed; //Pedestal + UShort_t fStatus; //Sucess or failure of fitting pocedure + Float_t fAmpSig; //Amplitude in entities of ADC counts + Float_t fT0; //Start time of signal in entities of sample intervals + Float_t fChi2Sig; //Chi Square of fit + UShort_t fNdfSig; //Number of degrees of freedom of fit + UShort_t fMinSig; //Pedestal +}; + +#endif diff --git a/EMCAL/AliEMCALRawAnalyzer.cxx b/EMCAL/AliEMCALRawAnalyzer.cxx new file mode 100644 index 00000000000..87385465025 --- /dev/null +++ b/EMCAL/AliEMCALRawAnalyzer.cxx @@ -0,0 +1,240 @@ +/************************************************************************** + * This file is property of and copyright by * + * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 * + * * + * Primary Author: Per Thomas Hille * + * * + * 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. * + **************************************************************************/ + + +// Base class for extraction +// of signal amplitude and peak position +// From EMCAL Calorimeter RAW data (from the RCU) +// Contains some utilities for preparing / selecting +// Signals suitable for signal extraction +// By derived classes + +#include "AliLog.h" +#include "AliEMCALRawAnalyzer.h" +#include "AliEMCALBunchInfo.h" +#include "AliEMCALFitResults.h" +#include +using namespace std; + + +AliEMCALRawAnalyzer::AliEMCALRawAnalyzer() : TObject(), + fMinTimeIndex(-1), + fMaxTimeIndex(-1), + fFitArrayCut(5), + fAmpCut(4), + fNsampleCut(5), + fIsZerosupressed( false ), + fVerbose( false ) +{ + //Comment + for(int i=0; i < MAXSAMPLES; i++ ) + { + fReversed[i] = 0; + } +} + +AliEMCALRawAnalyzer::~AliEMCALRawAnalyzer() +{ + +} + + +void +AliEMCALRawAnalyzer::SetTimeConstraint(const int min, const int max ) +{ + //Require that the bin if the maximum ADC value is between min and max (timebin) + if( ( min > max ) || min > MAXSAMPLES || max > MAXSAMPLES ) + { + AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored", min, max ) ); + } + else + { + fMinTimeIndex = min; + fMaxTimeIndex = max; + } +} + + +UShort_t +AliEMCALRawAnalyzer::Max(const UShort_t *data, const int length ) const +{ + //------------ + UShort_t tmpmax = data[0]; + + for(int i=0; i < length; i++) + { + if( tmpmax < data[i] ) + { + tmpmax = data[i]; + } + } + return tmpmax; +} + + +void +AliEMCALRawAnalyzer::SelectSubarray( const Double_t *fData, const int length, const short maxindex, int *const first, int *const last ) const +{ + //Selection of subset of data from one bunch that will be used for fitting or + //Peak finding. Go to the left and right of index of the maximum time bin + //Untile the ADC value is less that fFitArrayCut + int tmpfirst = maxindex; + int tmplast = maxindex; + + while(( tmpfirst ) > 0 && ( fData[tmpfirst] > fFitArrayCut )) + { + tmpfirst -- ; + } + + while(( tmplast ) < length && ( fData [tmplast] > fFitArrayCut )) + { + tmplast ++; + } + + *first = tmpfirst +1; + *last = tmplast -1; +} + + + +Float_t +AliEMCALRawAnalyzer::ReverseAndSubtractPed( const AliEMCALBunchInfo *bunch, const UInt_t /*altrocfg1*/, const UInt_t /*altrocfg2*/, double *outarray ) const +{ + //Time sample comes in reversed order, revers them back + //Subtract the baseline based on content of altrocfg1 and altrocfg2. + Int_t length = bunch->GetLength(); + const UShort_t *sig = bunch->GetData(); + + double ped = 0; + double tmp = 0; + + if( fIsZerosupressed == false ) + { + for(int i=0; i < 5; i++ ) + { + tmp += sig[i]; + } + } + ped = tmp / 5; + for( int i=0; i < length; i++ ) + { + outarray[i] = sig[length -i -1] - ped; + } + + return ped; +} + + +short +AliEMCALRawAnalyzer::Max( const AliEMCALBunchInfo *const bunch , int *const maxindex ) const +{ + //comment + short tmpmax = -1; + int tmpindex = -1; + const UShort_t *sig = bunch->GetData(); + + for(int i=0; i < bunch->GetLength(); i++ ) + { + if( sig[i] > tmpmax ) + { + tmpmax = sig[i]; + tmpindex = i; + } + } + + if(maxindex != 0 ) + { + *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); + } + + return tmpmax; +} + + +int +AliEMCALRawAnalyzer::SelectBunch( const vector &bunchvector,short *const maxampbin, short *const maxamplitude ) const +{ + //We select the bunch with the highest amplitude unless any time constraints is set + short max = -1; + short bunchindex = -1; + short maxall = -1; + int indx = -1; + + for(unsigned int i=0; i < bunchvector.size(); i++ ) + { + max = Max( &bunchvector.at(i), &indx ); + if( IsInTimeRange( indx) ) + { + if( max > maxall ) + { + maxall = max; + bunchindex = i; + } + } + } + + *maxampbin = indx; + *maxamplitude = max; + return bunchindex; +} + + +bool +AliEMCALRawAnalyzer::IsInTimeRange( const int maxindex ) const +{ + // Ckeck if the index of the max ADC vaue is consistent with trigger. + if( ( fMinTimeIndex < 0 && fMaxTimeIndex < 0) ||fMaxTimeIndex < 0 ) + { + return true; + } + return ( maxindex < fMaxTimeIndex ) && ( maxindex > fMinTimeIndex ) ? true : false; +} + + +void +AliEMCALRawAnalyzer::PrintBunches( const vector &bvctr ) const +{ + //comment + cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl; + cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl; + + for(unsigned int i=0; i < bvctr.size() ; i++ ) + { + PrintBunch( bvctr.at(i) ); + cout << " bunch = " << i << endl; + } + cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl; +} + + +void +AliEMCALRawAnalyzer::PrintBunch( const AliEMCALBunchInfo &bunch ) const +{ + //comment + cout << __FILE__ << ":" << __LINE__ << endl; + cout << __FILE__ << __LINE__ << ", startimebin = " << bunch.GetStartBin() << ", length =" << bunch.GetLength() << endl; + const UShort_t *sig = bunch.GetData(); + + for ( Int_t j = 0; j < bunch.GetLength(); j++) + { + printf("%d\t", sig[j] ); + } + cout << endl; +} + + diff --git a/EMCAL/AliEMCALRawAnalyzer.h b/EMCAL/AliEMCALRawAnalyzer.h new file mode 100644 index 00000000000..b4d103fbddf --- /dev/null +++ b/EMCAL/AliEMCALRawAnalyzer.h @@ -0,0 +1,74 @@ +#ifndef ALIEMCALRAWANALYZER_H +#define ALIEMCALRAWANALYZER_H +/************************************************************************** + * This file is property of and copyright by * + * the Relatvistic Heavy Ion Group (RHIG), Yale University, US, 2009 * + * * + * Primary Author: Per Thomas Hille * + * * + * 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. * + **************************************************************************/ + +//Base class for extraction +//of signal amplitude and peak position +//From EMCAL Calorimeter RAW data + + +#include "Rtypes.h" +#include "TObject.h" +#include +using namespace std; + +#define MAXSAMPLES 1008 //CRAP PTH +#include "AliEMCALRawAnalyzer.h" + +class AliEMCALBunchInfo; +class AliEMCALFitResults; + + +class AliEMCALRawAnalyzer : public TObject +{ + public: + AliEMCALRawAnalyzer(); + virtual ~AliEMCALRawAnalyzer(); + virtual AliEMCALFitResults Evaluate( const vector &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 ) = 0; + + void PrintBunches( const vector &bunchvector ) const; + void PrintBunch( const AliEMCALBunchInfo &bunch ) const ; + + void SetTimeConstraint(const int min, const int max ); + void SetVerbose(bool verbose = true){ fVerbose = verbose; }; + void SetIsZeroSuppressed(const bool iszs = true) { fIsZerosupressed = iszs; } ; + void SetAmpCut(const Float_t cut) { fAmpCut = cut ; } ; + + protected: + short Max( const AliEMCALBunchInfo *const bunch, int *const maxindex) const; + UShort_t Max(const UShort_t *data, const int length ) const; + bool IsInTimeRange( const int maxindex ) const; + Float_t ReverseAndSubtractPed( const AliEMCALBunchInfo *bunch, const UInt_t altrocfg1, const UInt_t altrocfg2, double *outarray ) const; + int SelectBunch( const vector &bunchvector, short *const maxampbin, short *const maxamplitude ) const; + void SelectSubarray( const Double_t *fData, const int length, const short maxindex, int *const first, int *const last ) const; + + Double_t fReversed[MAXSAMPLES]; //Reversed sequence of samples (pedestalsubtracted) + + // private: + int fMinTimeIndex; //The timebin of the max signal value must be between fMinTimeIndex and fMaxTimeIndex + int fMaxTimeIndex; //The timebin of the max signal value must be between fMinTimeIndex and fMaxTimeIndex + int fFitArrayCut; //Cut on ADC value (after ped. subtraction) for signals used for fit + Float_t fAmpCut; //Max ADC - pedestal must be higher than this befor attemting to extract the amplitude + int fNsampleCut; //Minimum number of sample require before attemting to extract signal parameters + bool fIsZerosupressed; //Wether or not the data is zeros supressed, by default its assumed that the baseline is also subtracted if set to true + bool fVerbose; //Print debug information to std out if set to true + +}; + +#endif diff --git a/EMCAL/AliEMCALRawAnalyzerLMS.cxx b/EMCAL/AliEMCALRawAnalyzerLMS.cxx new file mode 100644 index 00000000000..64b342418cb --- /dev/null +++ b/EMCAL/AliEMCALRawAnalyzerLMS.cxx @@ -0,0 +1,152 @@ +/************************************************************************** + * This file is property of and copyright by * + * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 * + * * + * Primary Author: Per Thomas Hille * + * * + * 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. * + **************************************************************************/ + + +// Extraction of amplitude and peak position +// FRom EMCAL raw data using +// least square fit for the +// Moment assuming identical and +// independent errors (equivalent with chi square) +// + +#include "AliEMCALRawAnalyzerLMS.h" +#include "AliEMCALBunchInfo.h" +#include "AliEMCALFitResults.h" +#include "AliLog.h" +#include "TMath.h" +#include +#include "TF1.h" +#include "TGraph.h" + +using namespace std; + + +#define BAD 4 //CRAP PTH + + +AliEMCALRawAnalyzerLMS::AliEMCALRawAnalyzerLMS() : AliEMCALRawAnalyzer(), + fkEulerSquared(7.389056098930650227), + fSig(0), + fTf1(0) +{ + //comment + for(int i=0; i < MAXSAMPLES; i++) + { + fXaxis[i] = i; + } + + fTf1 = new TF1( "myformula", "[0]*((x - [1])/2.35)^2*exp(-2*(x -[1])/2.35)", 0, 30 ); + // fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^[3]*exp(-[3]*(x -[1])/[2])", 0, 30 ); + +} + + +AliEMCALRawAnalyzerLMS::~AliEMCALRawAnalyzerLMS() +{ + delete fTf1; +} + + + +AliEMCALFitResults +AliEMCALRawAnalyzerLMS::Evaluate( const vector &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 ) +{ + // Extracting signal parameters using fitting + short maxampindex; //index of maximum amplitude + short maxamp; //Maximum amplitude + int index = SelectBunch( bunchvector, &maxampindex, &maxamp ); + + if( index >= 0) + { + Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed ); + int first; + int last; + Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed ); + + if ( maxf > fAmpCut ) + { + SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxampindex - bunchvector.at(index).GetStartBin(), &first, &last); + int nsamples = last - first; + + if( ( nsamples ) > fNsampleCut ) + { + + TGraph *graph = new TGraph( last - first, fXaxis, &fReversed[first] ); + fTf1->SetParameter(0, maxf*fkEulerSquared ); + fTf1->SetParameter(1, 0.235); + // fTf1->SetParameter(2, 2); + // fTf1->SetParameter(2, 3); + + // return AliEMCALFitResults( -1 , -1 , -1, -1, -1, -1 , -1); + + Short_t tmpStatus = graph->Fit(fTf1, "Q0RW"); + + // return AliEMCALFitResults( -1 , -1 , -1, -1, -1, -1 , -1); + + if( fVerbose == true ) + { + AliEMCALRawAnalyzer::PrintBunch( bunchvector.at(index) ); + PrintFitResult( fTf1 ) ; + } + + delete graph; + + + + return AliEMCALFitResults( maxamp, ped , tmpStatus, + fTf1->GetParameter(0)/fkEulerSquared, + fTf1->GetParameter(1) + maxampindex, + fTf1->GetChisquare(), + fTf1->GetNDF()); + + + // delete graph; + + } + else + { + return AliEMCALFitResults( maxamp, ped, -1, maxf, -1, -1, -1 ); + } + } + else + { + return AliEMCALFitResults( maxamp , ped, -1, maxf, -1, -1, -1 ); + } + } + else + { + return AliEMCALFitResults( -99, -99 ); + } + + return AliEMCALFitResults( -99, -99 ); + +} + + +void +AliEMCALRawAnalyzerLMS::PrintFitResult(const TF1 *f) const +{ + //comment + cout << endl; + cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl; + cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl; + cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl; + cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl; + // cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus() << ",.. !!!!" << endl << endl; + cout << endl << endl; +} diff --git a/EMCAL/AliEMCALRawAnalyzerLMS.h b/EMCAL/AliEMCALRawAnalyzerLMS.h new file mode 100644 index 00000000000..10c56bf9ca3 --- /dev/null +++ b/EMCAL/AliEMCALRawAnalyzerLMS.h @@ -0,0 +1,53 @@ +#ifndef ALIEMCALRAWANALYZERLMS_H +#define ALIEMCALRAWANALYZERLMS_H +/************************************************************************** + * This file is property of and copyright by * + * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 * + * * + * Primary Author: Per Thomas Hille * + * * + * 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. * + **************************************************************************/ + + +// Extraction of amplitude and peak position +// FRom EMCAL raw data using +// Chi square fit + +#include "AliEMCALRawAnalyzer.h" + + +class TF1; +class TGraph; + +class AliEMCALRawAnalyzerLMS : public AliEMCALRawAnalyzer +{ + public: + AliEMCALRawAnalyzerLMS(); + virtual ~AliEMCALRawAnalyzerLMS(); + virtual AliEMCALFitResults Evaluate( const vector &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 ); + void PrintFitResult(const TF1 *f) const; + + private: + AliEMCALRawAnalyzerLMS(const AliEMCALRawAnalyzerLMS & ); + AliEMCALRawAnalyzerLMS & operator = (const AliEMCALRawAnalyzerLMS &); + + double fXaxis[MAXSAMPLES]; //Axis if time bins, ( used by TGraph ) + const double fkEulerSquared; //e^2 = 7.389056098930650227 + TGraph *fSig; // Signale holding the data to be fitted. + TF1 *fTf1; // Analytical formula of the Semi Gaussian to be fitted + + + +}; + +#endif diff --git a/EMCAL/AliEMCALRawAnalyzerPeakFinder.cxx b/EMCAL/AliEMCALRawAnalyzerPeakFinder.cxx new file mode 100644 index 00000000000..ff1f5c22f21 --- /dev/null +++ b/EMCAL/AliEMCALRawAnalyzerPeakFinder.cxx @@ -0,0 +1,172 @@ +/************************************************************************** + * This file is property of and copyright by the Experimental Nuclear * + * Physics Group, Dep. of Physics * + * University of Oslo, Norway, 2007 * + * * + * Author: Per Thomas Hille for the ALICE HLT Project.* + * Contributors are mentioned in the code where appropriate. * + * Please report bugs to perthi@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 "AliEMCALRawAnalyzerPeakFinder.h" +#include "AliEMCALBunchInfo.h" +#include "AliEMCALFitResults.h" +#include +#include "unistd.h" +#include "TMath.h" +#include "AliLog.h" + +using namespace std; + +AliEMCALRawAnalyzerPeakFinder::AliEMCALRawAnalyzerPeakFinder() :AliEMCALRawAnalyzer(), + 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]; + + for(int k=0; k < 100; k++ ) + { + fPFAmpVectors[i][j][k] = 0; + fPFTofVectors[i][j][k] = 0; + } + } + } + LoadVectors(); +} + + +AliEMCALRawAnalyzerPeakFinder::~AliEMCALRawAnalyzerPeakFinder() +{ + //comment + for(int i=0; i < MAXSTART; i++) + { + for(int j=0; j < SAMPLERANGE; j++ ) + { + delete[] fPFAmpVectors[i][j]; + delete[] fPFTofVectors[i][j]; + } + } +} + + +AliEMCALFitResults +AliEMCALRawAnalyzerPeakFinder::Evaluate( const vector &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 ) +{ + // Extracting the amplitude using the Peak-Finder 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; + 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 AliEMCALFitResults( 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; + + for(int i=0; i < SAMPLERANGE; i++ ) + { + int dt = maxampindex - startbin -2; + double tmp = fReversed[ dt +i]; + fAmp += fPFAmpVectors[0][pfindex][i]*tmp; + } + + return AliEMCALFitResults( maxamp, ped , -1, fAmp, -1, -1, -1 ); + } + else + { + return AliEMCALFitResults( maxamp, ped , -1, maxf, -1, -1, -1 ); + } + } + } + + // cout << __FILE__ << __LINE__ << "WARNING, returning amp = -1 " << endl; + + return AliEMCALFitResults(-1, -1); +} + + +void +AliEMCALRawAnalyzerPeakFinder::LoadVectors() +{ + //Read in the Peak finder vecors from file + for(int i = 0; i < MAXSTART ; i++) + { + for( int j=0; j < SAMPLERANGE; j++) + { + char filename[256]; + int n = j+fNsampleCut; + double start = (double)i+0.5; + sprintf(filename, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n); + FILE *fp = fopen(filename, "r"); + + if(fp == 0) + { + AliFatal( Form( "could not open file: %s", filename ) ); + } + else + { + for(int m = 0; m < n ; m++ ) + { + fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] ); + } + fscanf(fp, "\n"); + + for(int m = 0; m < n ; m++ ) + { + fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] ); + } + + fclose (fp); + } + } + } +} diff --git a/EMCAL/AliEMCALRawAnalyzerPeakFinder.h b/EMCAL/AliEMCALRawAnalyzerPeakFinder.h new file mode 100644 index 00000000000..c7e1d7b5546 --- /dev/null +++ b/EMCAL/AliEMCALRawAnalyzerPeakFinder.h @@ -0,0 +1,55 @@ +#ifndef ALIEMCALRAWANALYZERPEAKFINDER_H +#define ALIEMCALRAWANALYZERPEAKFINDER_H + +/************************************************************************** + * This file is property of and copyright by the Experimental Nuclear * + * Physics Group, Dep. of Physics * + * University of Oslo, Norway, 2007 * + * * + * Author: Per Thomas Hille for the ALICE HLT Project.* + * Contributors are mentioned in the code where appropriate. * + * Please report bugs to perthi@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 "AliEMCALRawAnalyzer.h" + +#define MAXSTART 3 +#define SAMPLERANGE 15 +#define SHIF 0.5 + +class AliEMCALBunchInfo; + +class AliEMCALRawAnalyzerPeakFinder : public AliEMCALRawAnalyzer +{ + public: + AliEMCALRawAnalyzerPeakFinder(); + virtual ~AliEMCALRawAnalyzerPeakFinder(); + virtual AliEMCALFitResults Evaluate( const vector &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 ); + + private: + AliEMCALRawAnalyzerPeakFinder( const AliEMCALRawAnalyzerPeakFinder & ); + AliEMCALRawAnalyzerPeakFinder & operator = ( const AliEMCALRawAnalyzerPeakFinder & ); + + void LoadVectors(); + double *fPFAmpVectors[MAXSTART][SAMPLERANGE]; // Vectors for Amplitude extraction + double *fPFTofVectors[MAXSTART][SAMPLERANGE]; // Vectors for TOF extraction + double fTof; + double fAmp; + +}; + +#endif diff --git a/EMCAL/libEMCALrec.pkg b/EMCAL/libEMCALrec.pkg index 42db9c7de3a..eda174f59c6 100644 --- a/EMCAL/libEMCALrec.pkg +++ b/EMCAL/libEMCALrec.pkg @@ -9,7 +9,12 @@ AliEMCALTracker.cxx \ AliEMCALPID.cxx \ AliEMCALQADataMakerRec.cxx \ AliEMCALAodCluster.cxx \ -AliCaloNeuralFit.cxx +AliCaloNeuralFit.cxx \ +AliEMCALRawAnalyzer.cxx \ +AliEMCALRawAnalyzerLMS.cxx \ +AliEMCALRawAnalyzerPeakFinder.cxx \ +AliEMCALBunchInfo.cxx \ +AliEMCALFitResults.cxx HDRS= $(SRCS:.cxx=.h) -- 2.39.3