Updated version.
[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzerLMS.cxx
CommitLineData
d655d7dd 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 <p.t.hille@fys.uio.no> *
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
20// Extraction of amplitude and peak position
57839add 21// FRom CALO raw data using
d655d7dd 22// least square fit for the
23// Moment assuming identical and
24// independent errors (equivalent with chi square)
25//
26
57839add 27#include "AliCaloRawAnalyzerLMS.h"
28#include "AliCaloBunchInfo.h"
29#include "AliCaloFitResults.h"
d655d7dd 30#include "AliLog.h"
31#include "TMath.h"
507751ce 32#include <stdexcept>
d655d7dd 33#include <iostream>
34#include "TF1.h"
35#include "TGraph.h"
36
37using namespace std;
38
39
40#define BAD 4 //CRAP PTH
e37e3c84 41
42ClassImp( AliCaloRawAnalyzerLMS )
d655d7dd 43
44
48a2e3eb 45AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square Fit", "LMS"),
57839add 46 fkEulerSquared(7.389056098930650227),
f57baa2d 47 fTf1(0),
48 fTau(2.35),
49 fFixTau(kTRUE)
d655d7dd 50{
51 //comment
52 for(int i=0; i < MAXSAMPLES; i++)
53 {
54 fXaxis[i] = i;
55 }
56
f57baa2d 57 fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^2*exp(-2*(x -[1])/[2])", 0, 30 );
58 if (fFixTau) {
59 fTf1->FixParameter(2, fTau);
60 }
61 else {
62 fTf1->ReleaseParameter(2); // allow par. to vary
63 fTf1->SetParameter(2, fTau);
64 }
65
d655d7dd 66}
67
68
57839add 69AliCaloRawAnalyzerLMS::~AliCaloRawAnalyzerLMS()
d655d7dd 70{
71 delete fTf1;
72}
73
74
75
57839add 76AliCaloFitResults
77AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
d655d7dd 78{
79 // Extracting signal parameters using fitting
80 short maxampindex; //index of maximum amplitude
81 short maxamp; //Maximum amplitude
82 int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
83
84 if( index >= 0)
85 {
86 Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
d655d7dd 87 Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
f57baa2d 88 short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
507751ce 89 // timebinOffset is timebin value at maximum (maxrev)
f57baa2d 90 short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
2cd0ffda 91 if( maxf < fAmpCut || ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
d655d7dd 92 {
2cd0ffda 93 return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
94 }
95 else if ( maxf >= fAmpCut )
96 {
97 int first = 0;
98 int last = 0;
3b8fd9fe 99 SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last);
100 int nsamples = last - first + 1;
d655d7dd 101
f57baa2d 102 if( ( nsamples ) >= fNsampleCut )
d655d7dd 103 {
507751ce 104 Float_t tmax = (maxrev - first); // local tmax estimate
3b8fd9fe 105 TGraph *graph = new TGraph( nsamples, fXaxis, &fReversed[first] );
d655d7dd 106 fTf1->SetParameter(0, maxf*fkEulerSquared );
507751ce 107 fTf1->SetParameter(1, tmax - fTau);
108 // set rather loose parameter limits
109 fTf1->SetParLimits(0, 0.5*maxf*fkEulerSquared, 2*maxf*fkEulerSquared );
110 fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4);
f57baa2d 111
112 if (fFixTau) {
113 fTf1->FixParameter(2, fTau);
114 }
115 else {
116 fTf1->ReleaseParameter(2); // allow par. to vary
117 fTf1->SetParameter(2, fTau);
118 }
d655d7dd 119
507751ce 120 Short_t tmpStatus = 0;
121 try {
122 tmpStatus = graph->Fit(fTf1, "Q0RW");
123 }
124 catch (const std::exception & e) {
125 AliError( Form("TGraph Fit exception %s", e.what()) );
2cd0ffda 126 return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, timebinOffset,
127 timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
507751ce 128 }
129
d655d7dd 130 if( fVerbose == true )
131 {
57839add 132 AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
d655d7dd 133 PrintFitResult( fTf1 ) ;
134 }
507751ce 135 // global tmax
136 tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0
137 + fTf1->GetParameter(2); // +tau, makes sum tmax
d655d7dd 138
139 delete graph;
507751ce 140 return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar,
f57baa2d 141 fTf1->GetParameter(0)/fkEulerSquared,
507751ce 142 tmax,
143 timebinOffset,
f57baa2d 144 fTf1->GetChisquare(),
145 fTf1->GetNDF(),
146 AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
147
d655d7dd 148 // delete graph;
149
150 }
151 else
152 {
2cd0ffda 153 Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
154 Int_t ndf = last - first - 1; // nsamples - 2
155 return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset,
156 timebinOffset, chi2, ndf, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
d655d7dd 157 }
2cd0ffda 158 } // ampcut
d655d7dd 159 }
f57baa2d 160 return AliCaloFitResults( AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid );
d655d7dd 161
162}
163
164
165void
57839add 166AliCaloRawAnalyzerLMS::PrintFitResult(const TF1 *f) const
d655d7dd 167{
168 //comment
169 cout << endl;
170 cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
171 cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
172 cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
173 cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
174 // cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus() << ",.. !!!!" << endl << endl;
175 cout << endl << endl;
176}
f57baa2d 177