coverity fix for TFF geometry method
[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{
168c7b3c 51
52 fAlgo = Algo::kLMS;
d655d7dd 53 //comment
54 for(int i=0; i < MAXSAMPLES; i++)
55 {
56 fXaxis[i] = i;
57 }
58
f57baa2d 59 fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^2*exp(-2*(x -[1])/[2])", 0, 30 );
168c7b3c 60 if (fFixTau)
61 {
f57baa2d 62 fTf1->FixParameter(2, fTau);
168c7b3c 63 }
64 else
65 {
66 fTf1->ReleaseParameter(2); // allow par. to vary
67 fTf1->SetParameter(2, fTau);
68 }
d655d7dd 69}
70
71
57839add 72AliCaloRawAnalyzerLMS::~AliCaloRawAnalyzerLMS()
d655d7dd 73{
74 delete fTf1;
75}
76
77
78
57839add 79AliCaloFitResults
80AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
d655d7dd 81{
82 // Extracting signal parameters using fitting
83 short maxampindex; //index of maximum amplitude
84 short maxamp; //Maximum amplitude
85 int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
86
87 if( index >= 0)
88 {
89 Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
d655d7dd 90 Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
f57baa2d 91 short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
507751ce 92 // timebinOffset is timebin value at maximum (maxrev)
f57baa2d 93 short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
2cd0ffda 94 if( maxf < fAmpCut || ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
d655d7dd 95 {
168c7b3c 96 return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
2cd0ffda 97 }
98 else if ( maxf >= fAmpCut )
99 {
100 int first = 0;
101 int last = 0;
3b8fd9fe 102 SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last);
103 int nsamples = last - first + 1;
d655d7dd 104
f57baa2d 105 if( ( nsamples ) >= fNsampleCut )
d655d7dd 106 {
507751ce 107 Float_t tmax = (maxrev - first); // local tmax estimate
3b8fd9fe 108 TGraph *graph = new TGraph( nsamples, fXaxis, &fReversed[first] );
d655d7dd 109 fTf1->SetParameter(0, maxf*fkEulerSquared );
507751ce 110 fTf1->SetParameter(1, tmax - fTau);
111 // set rather loose parameter limits
112 fTf1->SetParLimits(0, 0.5*maxf*fkEulerSquared, 2*maxf*fkEulerSquared );
113 fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4);
f57baa2d 114
115 if (fFixTau) {
116 fTf1->FixParameter(2, fTau);
117 }
118 else {
119 fTf1->ReleaseParameter(2); // allow par. to vary
120 fTf1->SetParameter(2, fTau);
121 }
d655d7dd 122
507751ce 123 Short_t tmpStatus = 0;
124 try {
125 tmpStatus = graph->Fit(fTf1, "Q0RW");
126 }
127 catch (const std::exception & e) {
128 AliError( Form("TGraph Fit exception %s", e.what()) );
168c7b3c 129 return AliCaloFitResults( maxamp, ped, Ret::kNoFit, maxf, timebinOffset,
130 timebinOffset, Ret::kDummy, Ret::kDummy, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
507751ce 131 }
132
d655d7dd 133 if( fVerbose == true )
134 {
57839add 135 AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
d655d7dd 136 PrintFitResult( fTf1 ) ;
137 }
507751ce 138 // global tmax
139 tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0
140 + fTf1->GetParameter(2); // +tau, makes sum tmax
d655d7dd 141
142 delete graph;
168c7b3c 143 return AliCaloFitResults( maxamp, ped , Ret::kFitPar,
f57baa2d 144 fTf1->GetParameter(0)/fkEulerSquared,
507751ce 145 tmax,
146 timebinOffset,
f57baa2d 147 fTf1->GetChisquare(),
148 fTf1->GetNDF(),
168c7b3c 149 Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
f57baa2d 150
d655d7dd 151 // delete graph;
152
153 }
154 else
155 {
2cd0ffda 156 Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
157 Int_t ndf = last - first - 1; // nsamples - 2
168c7b3c 158 return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset,
159 timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
d655d7dd 160 }
2cd0ffda 161 } // ampcut
d655d7dd 162 }
168c7b3c 163 return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
d655d7dd 164
165}
166
167
168void
57839add 169AliCaloRawAnalyzerLMS::PrintFitResult(const TF1 *f) const
d655d7dd 170{
171 //comment
172 cout << endl;
173 cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
174 cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
175 cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
176 cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
177 // cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus() << ",.. !!!!" << endl << endl;
178 cout << endl << endl;
179}
f57baa2d 180