Line data Source code
1 : // -*- mode: c++ -*-
2 : /**************************************************************************
3 : * This file is property of and copyright by *
4 : * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 *
5 : * *
6 : * Primary Author: Per Thomas Hille <perthomas.hille@yale.edu> *
7 : * *
8 : * Contributors are mentioned in the code where appropriate. *
9 : * Please report bugs to p.t.hille@fys.uio.no *
10 : * *
11 : * Permission to use, copy, modify and distribute this software and its *
12 : * documentation strictly for non-commercial purposes is hereby granted *
13 : * without fee, provided that the above copyright notice appears in all *
14 : * copies and that both the copyright notice and this permission notice *
15 : * appear in the supporting documentation. The authors make no claims *
16 : * about the suitability of this software for any purpose. It is *
17 : * provided "as is" without express or implied warranty. *
18 : **************************************************************************/
19 :
20 :
21 : // Base class for extraction
22 : // of signal amplitude and peak position
23 : // From CALO Calorimeter RAW data (from the RCU)
24 : // Contains some utilities for preparing / selecting
25 : // Signals suitable for signal extraction
26 : // By derived classes
27 :
28 : #include "AliLog.h"
29 : #include "AliCaloRawAnalyzer.h"
30 : #include "AliCaloBunchInfo.h"
31 : #include "AliCaloFitResults.h"
32 : #include "TMath.h"
33 : #include <iostream>
34 : using namespace std;
35 :
36 42 : ClassImp(AliCaloRawAnalyzer)
37 :
38 20 : AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort) : TObject(),
39 20 : fMinTimeIndex(-1),
40 20 : fMaxTimeIndex(-1),
41 20 : fFitArrayCut(5),
42 20 : fAmpCut(4),
43 20 : fNsampleCut(5),
44 20 : fOverflowCut(950),
45 20 : fNsamplePed(3),
46 20 : fIsZerosupressed( false ),
47 20 : fVerbose( false ),
48 20 : fAlgo(Algo::kNONE),
49 20 : fL1Phase(0),
50 20 : fAmp(0),
51 20 : fTof(0),
52 20 : fTau( EMCAL::TAU ),
53 20 : fFixTau( true )
54 60 : {
55 : // Ctor
56 :
57 20 : snprintf(fName, 256, "%s", name);
58 20 : snprintf(fNameShort,256, "%s", nameshort);
59 :
60 40360 : for(int i=0; i < ALTROMAXSAMPLES; i++ )
61 : {
62 20160 : fReversed[i] = 0;
63 : }
64 20 : }
65 :
66 :
67 : void
68 : AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max )
69 : {
70 : //Require that the bin if the maximum ADC value is between min and max (timebin)
71 :
72 0 : if( ( min > max ) || min > ALTROMAXSAMPLES || max > ALTROMAXSAMPLES )
73 : {
74 0 : AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored", min, max ) );
75 0 : }
76 : else
77 : {
78 0 : fMinTimeIndex = min;
79 0 : fMaxTimeIndex = max;
80 : }
81 0 : }
82 :
83 :
84 : UShort_t
85 : AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
86 : {
87 : // Get maximum of array
88 :
89 0 : UShort_t tmpmax = data[0];
90 :
91 0 : for(int i=0; i < length; i++)
92 : {
93 0 : if( tmpmax < data[i] )
94 : {
95 : tmpmax = data[i];
96 0 : }
97 : }
98 0 : return tmpmax;
99 : }
100 :
101 :
102 : void
103 : AliCaloRawAnalyzer::SelectSubarray( const Double_t *data, int length, short maxindex,
104 : int * first, int * last, int cut) const
105 : {
106 : //Selection of subset of data from one bunch that will be used for fitting or
107 : //Peak finding. Go to the left and right of index of the maximum time bin
108 : //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
109 116 : int tmpfirst = maxindex;
110 : int tmplast = maxindex;
111 58 : Double_t prevFirst = data[maxindex];
112 : Double_t prevLast = data[maxindex];
113 : bool firstJump = false;
114 : bool lastJump = false;
115 :
116 651 : while( (tmpfirst >= 0) && (data[tmpfirst] >= cut ) && (!firstJump) )
117 : {
118 : // jump check:
119 121 : if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
120 63 : if ( data[tmpfirst] >= prevFirst) {
121 : firstJump = true;
122 15 : }
123 : }
124 121 : prevFirst = data[tmpfirst];
125 121 : tmpfirst -- ;
126 : }
127 :
128 1208 : while( (tmplast < length) && (data[tmplast] >= cut ) && (!lastJump) )
129 : {
130 : // jump check:
131 273 : if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
132 215 : if ( data[tmplast] >= prevLast) {
133 : lastJump = true;
134 1 : }
135 : }
136 273 : prevLast = data[tmplast];
137 273 : tmplast ++;
138 : }
139 :
140 : // we keep one pre- or post- sample if we can (as in online)
141 : // check though if we ended up on a 'jump', or out of bounds: if so, back up
142 90 : if (firstJump || tmpfirst<0) tmpfirst ++;
143 116 : if (lastJump || tmplast>=length) tmplast --;
144 :
145 58 : *first = tmpfirst;
146 58 : *last = tmplast;
147 :
148 : return;
149 58 : }
150 :
151 :
152 :
153 : Float_t
154 : AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch,
155 : UInt_t /*altrocfg1*/, UInt_t /*altrocfg2*/,
156 : double *outarray ) const
157 : {
158 : //Time sample comes in reversed order, revers them back
159 : //Subtract the baseline based on content of altrocfg1 and altrocfg2.
160 :
161 116 : Int_t length = bunch->GetLength();
162 58 : const UShort_t *sig = bunch->GetData();
163 :
164 58 : double ped = EvaluatePedestal( sig , length);
165 :
166 1196 : for( int i=0; i < length; i++ )
167 : {
168 540 : outarray[i] = sig[length -i -1] - ped;
169 : }
170 :
171 58 : return ped;
172 : }
173 :
174 :
175 :
176 : Float_t
177 : AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * data, int /*length*/ ) const
178 : {
179 : // Pedestal evaluation if not zero suppressed
180 :
181 : double tmp = 0;
182 :
183 116 : if( fIsZerosupressed == false )
184 : {
185 0 : for(int i=0; i < fNsamplePed; i++ )
186 : {
187 0 : tmp += data[i];
188 : }
189 0 : }
190 :
191 58 : return tmp/fNsamplePed;
192 : }
193 :
194 :
195 : short
196 : AliCaloRawAnalyzer::Max( const AliCaloBunchInfo * bunch , int * maxindex ) const
197 : {
198 : // Get maximum in bunch array
199 :
200 : short tmpmax = -1;
201 : int tmpindex = -1;
202 116 : const UShort_t *sig = bunch->GetData();
203 :
204 1196 : for(int i=0; i < bunch->GetLength(); i++ )
205 : {
206 540 : if( sig[i] > tmpmax )
207 : {
208 : tmpmax = sig[i];
209 : tmpindex = i;
210 383 : }
211 : }
212 :
213 58 : if(maxindex != 0 )
214 : {
215 : // *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
216 58 : *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
217 58 : }
218 :
219 58 : return tmpmax;
220 : }
221 :
222 :
223 : bool
224 : AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const
225 : {
226 : // a bunch is considered invalid if the maximum is in the first or last time-bin
227 : short tmpmax = -1;
228 : int tmpindex = -1;
229 116 : const UShort_t *sig = bunch->GetData();
230 :
231 1196 : for(int i=0; i < bunch->GetLength(); i++ )
232 : {
233 540 : if( sig[i] > tmpmax )
234 : {
235 : tmpmax = sig[i];
236 : tmpindex = i;
237 383 : }
238 : }
239 :
240 : bool bunchOK = true;
241 116 : if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) )
242 : {
243 : bunchOK = false;
244 0 : }
245 :
246 58 : return bunchOK;
247 : }
248 :
249 :
250 : int
251 : AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,
252 : short * maxampbin, short * maxamplitude )
253 : {
254 : //We select the bunch with the highest amplitude unless any time constraints is set
255 : short max = -1;
256 : short bunchindex = -1;
257 : short maxall = -1;
258 116 : int indx = -1;
259 :
260 232 : for(unsigned int i=0; i < bunchvector.size(); i++ )
261 : {
262 58 : max = Max( &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches
263 58 : if( IsInTimeRange( indx, fMaxTimeIndex, fMinTimeIndex) )
264 : {
265 58 : if( max > maxall )
266 : {
267 : maxall = max;
268 58 : bunchindex = i;
269 58 : *maxampbin = indx;
270 58 : *maxamplitude = max;
271 58 : }
272 : }
273 : }
274 :
275 58 : if (bunchindex >= 0) {
276 58 : bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
277 58 : if (! bunchOK) {
278 : bunchindex = -1;
279 0 : }
280 58 : }
281 :
282 116 : return bunchindex;
283 58 : }
284 :
285 :
286 : bool
287 : AliCaloRawAnalyzer::IsInTimeRange( const int maxindex, const int maxtindx, const int mintindx ) const
288 : {
289 : // Ckeck if the index of the max ADC vaue is consistent with trigger.
290 116 : if( ( mintindx < 0 && maxtindx < 0) ||maxtindx < 0 )
291 : {
292 58 : return true;
293 : }
294 0 : return ( maxindex < maxtindx ) && ( maxindex > mintindx ) ? true : false;
295 58 : }
296 :
297 :
298 : void
299 : AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr )
300 : {
301 : // Print bunch vector info
302 0 : cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
303 0 : cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;
304 :
305 0 : for(unsigned int i=0; i < bvctr.size() ; i++ )
306 : {
307 0 : PrintBunch( bvctr.at(i) );
308 0 : cout << " bunch = " << i << endl;
309 : }
310 0 : cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl;
311 0 : }
312 :
313 :
314 : void
315 : AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch )
316 : {
317 : // Print bunch infor
318 0 : cout << __FILE__ << ":" << __LINE__ << endl;
319 0 : cout << __FILE__ << __LINE__ << ", startimebin = " << bunch.GetStartBin() << ", length =" << bunch.GetLength() << endl;
320 0 : const UShort_t *sig = bunch.GetData();
321 :
322 0 : for ( Int_t j = 0; j < bunch.GetLength(); j++)
323 : {
324 0 : printf("%d\t", sig[j] );
325 : }
326 0 : cout << endl;
327 0 : }
328 :
329 :
330 : Double_t
331 : AliCaloRawAnalyzer::CalculateChi2( Double_t amp, Double_t time,
332 : Int_t first, Int_t last,
333 : Double_t adcErr, Double_t tau) const
334 : {
335 : // Input:
336 : // amp - max amplitude;
337 : // time - time of max amplitude;
338 : // first, last - sample array indices to be used
339 : // adcErr - nominal error of amplitude measurement (one value for all channels)
340 : // if adcErr<0 that mean adcErr=1.
341 : // tau - filter time response (in timebin units)
342 : // Output:
343 : // chi2 - chi2
344 :
345 0 : if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined..
346 : // or, first is negative, the indices are not valid
347 0 : return Ret::kDummy;
348 : }
349 :
350 0 : int nsamples = last - first + 1;
351 : // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time);
352 :
353 : Int_t x = 0;
354 : Double_t chi2 = 0;
355 : Double_t dy = 0.0, xx = 0.0, f=0.0;
356 :
357 0 : for (Int_t i=0; i<nsamples; i++) {
358 0 : x = first + i; // timebin
359 0 : xx = (x - time + tau) / tau; // help variable
360 : f = 0;
361 0 : if (xx > 0) {
362 0 : f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ;
363 0 : }
364 0 : dy = fReversed[x] - f;
365 0 : chi2 += dy*dy;
366 : // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy);
367 : }
368 :
369 0 : if (adcErr>0.0) { // weight chi2
370 0 : chi2 /= (adcErr*adcErr);
371 0 : }
372 : return chi2;
373 0 : }
374 :
375 :
376 : void
377 : AliCaloRawAnalyzer::CalculateMeanAndRMS(Int_t first, Int_t last,
378 : Double_t & mean, Double_t & rms)
379 : {
380 : // Input:
381 : // first, last - sample array indices to be used
382 : // Output:
383 : // mean and RMS of samples
384 : //
385 : // To possibly be used to differentiate good signals from bad before fitting
386 : //
387 0 : mean = Ret::kDummy;
388 0 : rms = Ret::kDummy;
389 :
390 0 : if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined..
391 : // or, first is negative, the indices are not valid
392 : return;
393 : }
394 :
395 0 : int nsamples = last - first + 1;
396 : // printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples);
397 :
398 : int x = 0;
399 : Double_t sampleSum = 0; // sum of samples
400 : Double_t squaredSampleSum = 0; // sum of samples squared
401 :
402 0 : for (Int_t i=0; i<nsamples; i++) {
403 0 : x = first + i;
404 0 : sampleSum += fReversed[x];
405 0 : squaredSampleSum += (fReversed[x] * fReversed[x]);
406 : }
407 :
408 0 : mean = sampleSum / nsamples;
409 0 : Double_t squaredMean = squaredSampleSum / nsamples;
410 : // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
411 0 : rms = sqrt(squaredMean - mean*mean);
412 :
413 : return;
414 0 : }
415 :
416 :
417 :
418 : int
419 : AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo> &bunchvector,
420 : UInt_t altrocfg1, UInt_t altrocfg2, Int_t & index,
421 : Float_t & maxf, short & maxamp,
422 : short & maxrev, Float_t & ped,
423 : int & first, int & last, int acut )
424 : {
425 : // method to do the selection of what should possibly be fitted
426 :
427 : int nsamples = 0;
428 116 : short maxampindex = 0;
429 58 : index = SelectBunch( bunchvector, &maxampindex, &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
430 :
431 :
432 116 : if( index >= 0 && maxamp >= acut ) // something valid was found, and non-zero amplitude
433 : {
434 : // use more convenient numbering and possibly subtract pedestal
435 58 : ped = ReverseAndSubtractPed( &(bunchvector.at(index)), altrocfg1, altrocfg2, fReversed );
436 58 : maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
437 :
438 58 : if ( maxf >= acut ) // possibly significant signal
439 : {
440 : // select array around max to possibly be used in fit
441 58 : maxrev = maxampindex - bunchvector.at(index).GetStartBin();
442 58 : SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, acut );
443 :
444 : // sanity check: maximum should not be in first or last bin
445 : // if we should do a fit
446 116 : if (first!=maxrev && last!=maxrev)
447 : {
448 : // calculate how many samples we have
449 58 : nsamples = last - first + 1;
450 58 : }
451 : }
452 : }
453 :
454 58 : return nsamples;
455 58 : }
456 :
|
