]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - HLT/TPCLib/AliHLTTPCPad.cxx
Updated SNM Glauber fit
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCPad.cxx
... / ...
CommitLineData
1// $Id$
2
3//**************************************************************************
4//* This file is property of and copyright by the ALICE HLT Project *
5//* ALICE Experiment at CERN, All rights reserved. *
6//* *
7//* Primary Authors: Timm Steinbeck, Matthias Richter *
8//* Developers: Kenneth Aamodt <kenneth.aamodt@student.uib.no> *
9//* for The ALICE HLT Project. *
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/// @file AliHLTTPCPad.cxx
21/// @author Matthias Richter, Kenneth Aamodt
22/// @date
23/// @brief Container Class for TPC Pads.
24///
25
26#include <cerrno>
27#include "AliHLTTPCPad.h"
28#include "AliHLTStdIncludes.h"
29
30
31//added by kenneth
32#include "AliHLTTPCTransform.h"
33#include "AliHLTTPCClusters.h"
34#include <sys/time.h>
35#include "TMath.h"
36#include "TFile.h"
37//------------------------------
38
39/** margin for the base line be re-avaluated */
40#define ALIHLTPAD_BASELINE_MARGIN (2*fAverage)
41
42/** ROOT macro for the implementation of ROOT specific class methods */
43ClassImp(AliHLTTPCPad)
44
45AliHLTTPCPad::AliHLTTPCPad()
46 :
47 fClusterCandidates(),
48 fUsedClusterCandidates(),
49 fSelectedPad(kFALSE),
50 fHWAddress(0),
51 fRowNo(-1),
52 fPadNo(-1),
53 fThreshold(0),
54 fAverage(-1),
55 fNofEvents(0),
56 fSum(0),
57 fCount(0),
58 fTotal(0),
59 fBLMax(-1),
60 fBLMaxBin(-1),
61 fBLMin(-1),
62 fBLMinBin(-1),
63 fFirstBLBin(0),
64 fNofBins(0),
65 fReadPos(0),
66 fpRawData(NULL),
67 fDataSignals(NULL),
68 fSignalPositionArray(NULL),
69 fSizeOfSignalPositionArray(0),
70 fNGoodSignalsSent(0),
71 fCandidateDigitsVector()
72{
73 // see header file for class documentation
74 // or
75 // refer to README to build package
76 // or
77 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
78 // HLTInfo("Entering default constructor");
79 fDataSignals= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
80 memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
81
82 fSignalPositionArray= new Int_t[AliHLTTPCTransform::GetNTimeBins()];
83 memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
84 fSizeOfSignalPositionArray=0;
85
86}
87
88AliHLTTPCPad::AliHLTTPCPad(Int_t /*dummy*/)
89 :
90 fClusterCandidates(),
91 fUsedClusterCandidates(),
92 fSelectedPad(kFALSE),
93 fHWAddress(0),
94 fRowNo(-1),
95 fPadNo(-1),
96 fThreshold(0),
97 fAverage(-1),
98 fNofEvents(0),
99 fSum(0),
100 fCount(0),
101 fTotal(0),
102 fBLMax(-1),
103 fBLMaxBin(-1),
104 fBLMin(-1),
105 fBLMinBin(-1),
106 fFirstBLBin(0),
107 fNofBins(0),
108 fReadPos(0),
109 fpRawData(NULL),
110 fDataSignals(NULL),
111 fSignalPositionArray(NULL),
112 fSizeOfSignalPositionArray(0),
113 fNGoodSignalsSent(0),
114 fCandidateDigitsVector()
115{
116 // see header file for class documentation
117 // or
118 // refer to README to build package
119 // or
120 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
121}
122
123AliHLTTPCPad::AliHLTTPCPad(Int_t offset, Int_t nofBins)
124 :
125 fClusterCandidates(),
126 fUsedClusterCandidates(),
127 fSelectedPad(kFALSE),
128 fHWAddress(0),
129 fRowNo(-1),
130 fPadNo(-1),
131 fThreshold(0),
132 fAverage(-1),
133 fNofEvents(0),
134 fSum(0),
135 fCount(0),
136 fTotal(0),
137 fBLMax(-1),
138 fBLMaxBin(-1),
139 fBLMin(-1),
140 fBLMinBin(-1),
141 fFirstBLBin(offset),
142 fNofBins(nofBins),
143 fReadPos(0),
144 fpRawData(NULL),
145 fDataSignals(NULL),
146 fSignalPositionArray(NULL),
147 fSizeOfSignalPositionArray(0),
148 fNGoodSignalsSent(0),
149 fCandidateDigitsVector()
150{
151 // see header file for class documentation
152}
153
154AliHLTTPCPad::~AliHLTTPCPad()
155{
156 // see header file for class documentation
157 if (fpRawData) {
158 HLTWarning("event data acquisition not stopped");
159 StopEvent();
160 }
161 if (fDataSignals) {
162 delete [] fDataSignals;
163 fDataSignals=NULL;
164 }
165 if (fSignalPositionArray!=NULL) {
166 delete [] fSignalPositionArray;
167 fSignalPositionArray=NULL;
168 }
169}
170
171Int_t AliHLTTPCPad::SetID(Int_t rowno, Int_t padno)
172{
173 // see header file for class documentation
174 fRowNo=rowno;
175 fPadNo=padno;
176
177 return 0;
178}
179
180Int_t AliHLTTPCPad::StartEvent()
181{
182 // see header file for class documentation
183 Int_t iResult=0;
184 if (fpRawData==NULL) {
185 fBLMax=-1;
186 fBLMaxBin=-1;
187 fBLMin=-1;
188 fBLMinBin=-1;
189 fSum=0;
190 fCount=0;
191 fTotal=0;
192 if (fNofBins>0) {
193 fpRawData=new AliHLTTPCSignal_t[fNofBins];
194 if (fpRawData) {
195 for (int i=0; i<fNofBins; i++) fpRawData[i]=-1;
196 } else {
197 HLTError("memory allocation failed");
198 iResult=-ENOMEM;
199 }
200 }
201 } else {
202 HLTWarning("event data acquisition already started");
203 iResult=-EALREADY;
204 }
205 return iResult;
206}
207
208Int_t AliHLTTPCPad::CalculateBaseLine(Int_t reqMinCount)
209{
210 // see header file for class documentation
211 Int_t iResult=0;
212 AliHLTTPCSignal_t avBackup=fAverage;
213 //HLTDebug("reqMinCount=%d fCount=%d fTotal=%d fSum=%d fBLMax=%d fBLMin=%d", reqMinCount, fCount, fTotal, fSum, fBLMax, fBLMin);
214 if (fCount>=reqMinCount && fCount>=fTotal/2) {
215 fAverage=fCount>0?fSum/fCount:0;
216 if (fAverage>0) {
217 //HLTDebug("average for current event %d (%d - %d)", fAverage, fBLMax, fBLMin);
218 fCount=0;fSum=-1;
219 if (fBLMax>ALIHLTPAD_BASELINE_MARGIN) {
220 // calculate again
221 //HLTDebug("maximum value %d exceeds margin for base line (%d) "
222 // "-> re-evaluate base line", fBLMax, ALIHLTPAD_BASELINE_MARGIN);
223 if (fpRawData) {
224 for (Int_t i=fFirstBLBin; i<fNofBins; i++)
225 if (fpRawData[i]>=0) AddBaseLineValue(i, fpRawData[i]);
226 if (fCount>0 && fCount>=reqMinCount && fCount>=fTotal/2) {
227 fAverage=fSum/fCount;
228 //HLTDebug("new average %d", fAverage);
229 } else {
230 // HLTDebug("baseline re-eveluation skipped because of to few "
231 // "contributing bins: total=%d, contributing=%d, req=%d"
232 // "\ndata might be already zero suppressed"
233 // , fTotal, fCount, reqMinCount);
234 iResult=-ENODATA;
235 }
236 fCount=0;fSum=-1;
237 } else {
238 HLTError("missing raw data for base line calculation");
239 iResult=-ENOBUFS;
240 }
241 }
242 if (iResult>=0) {
243 // calculate average for all events
244 fAverage=((avBackup*fNofEvents)+fAverage)/(fNofEvents+1);
245 //HLTDebug("base line average for %d event(s): %d", fNofEvents+1, fAverage);
246 } else {
247 fAverage=avBackup;
248 }
249 } else {
250 fAverage=avBackup;
251 }
252 } else {
253 // HLTDebug("baseline calculation skipped because of to few contributing "
254 // "bins: total=%d, contributing=%d, required=%d \ndata might be "
255 // "already zero suppressed", fTotal, fCount, reqMinCount);
256 }
257
258 return iResult;
259}
260
261Int_t AliHLTTPCPad::StopEvent()
262{
263 // see header file for class documentation
264 Int_t iResult=0;
265 if (fpRawData) {
266 AliHLTTPCSignal_t* pData=fpRawData;
267 fpRawData=NULL;
268 delete [] pData;
269 fTotal=0;
270 fNofEvents++;
271 Rewind();
272 } else if (fNofBins>0) {
273 HLTError("event data acquisition not started");
274 iResult=-EBADF;
275 }
276 return iResult;
277}
278
279Int_t AliHLTTPCPad::ResetHistory()
280{
281 // see header file for class documentation
282 Int_t iResult=0;
283 fAverage=-1;
284 fNofEvents=0;
285 return iResult;
286}
287
288Int_t AliHLTTPCPad::SetThreshold(AliHLTTPCSignal_t thresh)
289{
290 // see header file for class documentation
291 Int_t iResult=0;
292 fThreshold=thresh;
293 return iResult;
294}
295
296Int_t AliHLTTPCPad::AddBaseLineValue(Int_t bin, AliHLTTPCSignal_t value)
297{
298 // see header file for class documentation
299 Int_t iResult=0;
300 if (bin>=fFirstBLBin) {
301 if (fAverage<0 || value<ALIHLTPAD_BASELINE_MARGIN) {
302 // add to the current sum and count
303 fSum+=value;
304 fCount++;
305 if (fBLMax<value) {
306 // keep the maximum value for later quality control of the base
307 // line calculation
308 fBLMax=value;
309 fBLMaxBin=bin;
310 }
311 if (fBLMin<0 || fBLMin>value) {
312 // keep the minimum value for later quality control of the base
313 // line calculation
314 fBLMin=value;
315 fBLMinBin=bin;
316 }
317 } else {
318 // HLTDebug("ignoring value %d (bin %d) for base line calculation "
319 // "(current average is %d)",
320 // value, bin, fAverage);
321 }
322 }
323 return iResult;
324}
325
326Int_t AliHLTTPCPad::SetRawData(Int_t bin, AliHLTTPCSignal_t value)
327{
328 // see header file for class documentation
329 // printf("Row: %d Pad: %d Time: %d Charge %d", fRowNo, fPadNo, bin, value);
330 Int_t iResult=0;
331 if (fpRawData) {
332 if (bin<fNofBins) {
333 if (value>=0) {
334 if (fpRawData[bin]<0) {
335 AddBaseLineValue(bin, value);
336 fTotal++;
337 } else {
338 // ignore value for average calculation
339 HLTWarning("overriding content of bin %d (%d)", bin, fpRawData[bin]);
340 }
341 fpRawData[bin]=value;
342 } else {
343 HLTWarning("ignoring neg. raw data");
344 }
345 } else {
346 HLTWarning("bin %d out of range (%d)", bin, fNofBins);
347 iResult=-ERANGE;
348 }
349 } else if (fNofBins>0) {
350 HLTError("event cycle not started");
351 iResult=-EBADF;
352 }
353 return iResult;
354}
355
356Int_t AliHLTTPCPad::Next(Int_t bZeroSuppression)
357{
358 // see header file for class documentation
359 if (fpRawData==NULL) return 0;
360 Int_t iResult=fReadPos<fNofBins;
361 if (iResult>0 && (iResult=(++fReadPos<fNofBins))>0) {
362 if (bZeroSuppression) {
363 while ((iResult=(fReadPos<fNofBins))>0 &&
364 GetCorrectedData(fReadPos)<=0)
365 fReadPos++;
366 }
367 }
368 return iResult;
369}
370
371Int_t AliHLTTPCPad::Rewind(Int_t bZeroSuppression)
372{
373 // see header file for class documentation
374 fReadPos=(bZeroSuppression>0?0:fFirstBLBin)-1;
375 return Next(bZeroSuppression);
376}
377
378AliHLTTPCSignal_t AliHLTTPCPad::GetRawData(Int_t bin) const
379{
380 // see header file for class documentation
381 AliHLTTPCSignal_t data=0;
382 if (fpRawData) {
383 if (bin<fNofBins) {
384 data=fpRawData[bin];
385 } else {
386 HLTWarning("requested bin %d out of range (%d)", bin, fNofBins);
387 }
388 } else if (fNofBins>0) {
389 HLTWarning("data only available within event cycle");
390 }
391 return data;
392}
393
394AliHLTTPCSignal_t AliHLTTPCPad::GetCorrectedData(Int_t bin) const
395{
396 // see header file for class documentation
397 AliHLTTPCSignal_t data=GetRawData(bin)-GetBaseLine(bin);
398 AliHLTTPCSignal_t prev=0;
399 if (bin>1) prev=GetRawData(bin-1)-GetBaseLine(bin-1);
400 AliHLTTPCSignal_t succ=0;
401 if (bin+1<GetSize()) succ=GetRawData(bin+1)-GetBaseLine(bin+1);
402 if (fThreshold>0) {
403 data-=fThreshold;
404 prev-=fThreshold;
405 succ-=fThreshold;
406 }
407
408 // case 1:
409 // the signal is below the base-line and threshold
410 if (data<0) data=0;
411
412 //case 2:
413 // the neighboring bins are both below base-line/threshold
414 // a real signal is always more than one bin wide because of the shaper
415 if (prev<=0 && succ<=0) data=0;
416
417 // case 3:
418 // the bin is inside the range of ignored bins
419 if (bin<fFirstBLBin) data=0;
420 //HLTDebug("fReadPos=%d data=%d threshold=%d raw data=%d base line=%d", fReadPos, data, fThreshold, GetRawData(bin), GetBaseLine(bin));
421 return data;
422}
423
424AliHLTTPCSignal_t AliHLTTPCPad::GetBaseLine(Int_t /*bin*/) const //TODO: Why is bin being ignored?
425{
426 // see header file for class documentation
427 AliHLTTPCSignal_t val=0;
428 if (fAverage>0) {
429 // we take the minumum value as the base line if it doesn't differ from
430 // the average to much
431 val=fAverage;
432#ifdef KEEP_NOISE
433 const AliHLTTPCSignal_t kMaxDifference=15;
434 if ((fAverage-fBLMin)<=kMaxDifference) val=fBLMin;
435 else val>kMaxDifference?val-=kMaxDifference:0;
436#endif
437 }
438 if (val<0) {
439 // here we should never get
440 val=0;
441 HLTFatal("wrong base line value");
442 }
443 return val;
444}
445
446AliHLTTPCSignal_t AliHLTTPCPad::GetAverage() const
447{
448 // see header file for class documentation
449 return fAverage>0?fAverage:0;
450}
451
452Float_t AliHLTTPCPad::GetOccupancy() const
453{
454 // see header file for class documentation
455 Float_t occupancy=0;
456 if (fpRawData && fNofBins>0) {
457 for (Int_t i=fFirstBLBin; i<fNofBins; i++) {
458 if (GetCorrectedData(i)>0) occupancy+=1;
459 }
460 if (fNofBins-fFirstBLBin>0)
461 occupancy/=fNofBins-fFirstBLBin;
462 }
463 return occupancy;
464}
465
466Float_t AliHLTTPCPad::GetAveragedOccupancy() const
467{
468 // see header file for class documentation
469
470 // history is not yet implemented
471 return GetOccupancy();
472}
473void AliHLTTPCPad::PrintRawData()
474{
475 // see header file for class documentation
476 for(Int_t bin=0;bin<AliHLTTPCTransform::GetNTimeBins();bin++){
477 if(GetDataSignal(bin)>0)
478 //This cout should be here since using logging produces output that is much more difficult to read
479 cout<<fRowNo<<"\t"<<fPadNo<<"\t"<<bin<<"\t"<<GetDataSignal(bin)<<endl;
480 }
481}
482
483void AliHLTTPCPad::ClearCandidates(){
484 fClusterCandidates.clear();
485 fUsedClusterCandidates.clear();
486 fCandidateDigitsVector.clear();
487}
488
489void AliHLTTPCPad::SetDataToDefault()
490{
491 // see header file for class documentation
492 // if(fDataSignals && fSignalPositionArray){
493 for(Int_t i =0;i<fSizeOfSignalPositionArray;i++){
494 fDataSignals[fSignalPositionArray[i]]=-1;
495 }
496 fSizeOfSignalPositionArray=0;
497 fNGoodSignalsSent = 0;
498 // }
499}
500
501void AliHLTTPCPad::SetDataSignal(Int_t bin,Int_t signal)
502{
503 // see header file for class documentation
504 fDataSignals[bin]=signal;
505 fSignalPositionArray[fSizeOfSignalPositionArray]=bin;
506 fSizeOfSignalPositionArray++;
507}
508
509Bool_t AliHLTTPCPad::GetNextGoodSignal(Int_t &time,Int_t &bunchSize){
510
511 if(fNGoodSignalsSent<fSizeOfSignalPositionArray&&fSizeOfSignalPositionArray>0){
512 time = fSignalPositionArray[fNGoodSignalsSent];
513 bunchSize=1;
514 fNGoodSignalsSent++;
515 while(fNGoodSignalsSent<fSizeOfSignalPositionArray){
516 if(fDataSignals[time+bunchSize+1]>0){
517 bunchSize++;
518 fNGoodSignalsSent++;
519 }
520 else{
521 break;
522 }
523 }
524 // fNGoodSignalsSent++;
525 return kTRUE;
526 }
527 return kFALSE;
528}
529
530Int_t AliHLTTPCPad::GetDataSignal(Int_t bin) const
531{
532 // see header file for class documentation
533 return fDataSignals[bin];
534}
535
536void AliHLTTPCPad::ZeroSuppress(Double_t nRMS, Int_t threshold, Int_t reqMinPoint, Int_t beginTime, Int_t endTime, Int_t timebinsLeft, Int_t timebinsRight, Int_t valueUnderAverage, bool speedup){
537 //see headerfile for documentation
538
539 //HLTDebug("In Pad: nRMS=%d, threshold=%d, reqMinPoint=%d, beginTime=%d, endTime=%d, timebinsLeft=%d timebinsRight=%d valueUnderAverage=%d \n",nRMS,threshold,reqMinPoint,beginTime,endTime,timebinsLeft,timebinsRight,valueUnderAverage);
540
541 Bool_t useRMS= kFALSE;
542 if(nRMS>0){
543 useRMS=kTRUE;
544 if(threshold>0){
545 HLTInfo("Both RMSThreshold and SignalThreshold defined, using RMSThreshold");
546 }
547 }
548 if(threshold<1 && nRMS<=0){
549 //setting the data to -1 for this pad
550 HLTInfo("Neither of RMSThreshold and SignalThreshold set, zerosuppression aborted");
551 return;
552 }
553
554 Int_t fThresholdUsed=threshold;
555
556 Int_t maxVal=0;
557 Int_t nAdded=0;
558 Int_t sumNAdded=0;
559 fSizeOfSignalPositionArray=0;
560 if(useRMS){
561 for(Int_t i=beginTime;i<endTime+1;i++){
562 if(fDataSignals[i]>0){
563 nAdded++;
564 sumNAdded+=fDataSignals[i]*fDataSignals[i];
565 if (maxVal<fDataSignals[i]) maxVal=fDataSignals[i];
566 }
567 }
568 }
569 else if(threshold>0){
570 for(Int_t i=beginTime;i<endTime+1;i++){
571 if(fDataSignals[i]>0){
572 nAdded++;
573 sumNAdded+=fDataSignals[i];
574 if (maxVal<fDataSignals[i]) maxVal=fDataSignals[i];
575 }
576 }
577 }
578 else{
579 HLTFatal("This should never happen because this is tested earlier in the code.(nRMSThreshold<1&&signal-threshold<1)");
580 }
581 if(nAdded<reqMinPoint){
582 HLTInfo("Number of signals is less than required, zero suppression aborted");
583 return;
584 }
585
586 if(nAdded==0){
587 HLTInfo("No signals added for this pad, zerosuppression aborted: pad %d row %d",fPadNo,fRowNo);
588 return;
589 }
590
591 Double_t averageValue=(Double_t)sumNAdded/nAdded;//true average for threshold approach, average of signals squared for rms approach
592
593 if(useRMS){
594 //Calculate the RMS
595 if(averageValue>0){
596 fThresholdUsed =(Int_t)(TMath::Sqrt(averageValue)*nRMS);
597 }
598 else{
599 HLTFatal("average value in ZeroSuppression less than 0, investigation needed. This should never happen");
600 }
601 }
602 else{
603 fThresholdUsed = (Int_t)(averageValue + threshold);
604 }
605 if (maxVal<fThresholdUsed) return;
606
607 // Do zero suppression on the adc values within [beginTime,endTime](add the good values)
608 for(Int_t i=beginTime;i<endTime;i++){
609 if(fDataSignals[i]>fThresholdUsed){
610 Int_t firstSignalTime=i;
611 for(Int_t left=1;left<timebinsLeft;left++){//looking 5 to the left of the signal to add tail
612 if(fDataSignals[i-left]-averageValue+valueUnderAverage>0 && i-left>=beginTime){
613 firstSignalTime--;
614 }
615 else{
616 break;
617 }
618 }
619 Int_t lastSignalTime=i;
620 while(fDataSignals[lastSignalTime+1]>fThresholdUsed && lastSignalTime+1<endTime){
621 lastSignalTime++;
622 }
623 for(Int_t right=1;right<timebinsRight;right++){//looking 5 to the left of the signal to add tail
624 if(fDataSignals[i+right]-averageValue+valueUnderAverage>0&&i+right<endTime){
625 lastSignalTime++;
626 }
627 else{
628 break;
629 }
630 }
631
632 for(Int_t t=firstSignalTime;t<lastSignalTime;t++){
633 fDataSignals[t]=(AliHLTTPCSignal_t)(fDataSignals[t]-averageValue + valueUnderAverage);
634 fSignalPositionArray[fSizeOfSignalPositionArray]=t;
635 fSizeOfSignalPositionArray++;
636 // Matthias Oct 10 2008: trying hard to make the code faster for the
637 // AltroChannelSelection. For that we only need to know there is data
638 if (speedup) return;
639 }
640 i+=lastSignalTime;
641 }
642 }
643 //reset the rest of the data
644 Int_t counterSize=fSizeOfSignalPositionArray;
645
646 for(Int_t d=endTime;d>=beginTime;d--){
647 if(d==fSignalPositionArray[counterSize-1]&&counterSize-1>=0){
648 counterSize--;
649 }
650 else{
651 fDataSignals[d]=-1;
652 }
653 }
654 if(fDataSignals[beginTime+1]<1){
655 fDataSignals[beginTime]=0;
656 }
657}
658
659void AliHLTTPCPad::AddClusterCandidate(const AliHLTTPCClusters& candidate){
660 fClusterCandidates.push_back(candidate);
661 fUsedClusterCandidates.push_back(0);
662}
663
664void AliHLTTPCPad::AddCandidateDigits(const vector<AliHLTTPCDigitData>& candidateDigits){
665 fCandidateDigitsVector.push_back(candidateDigits);
666}
667
668vector<AliHLTTPCDigitData> *AliHLTTPCPad::GetCandidateDigits(Int_t candidateIndex){
669 return (size_t(candidateIndex) < fCandidateDigitsVector.size()) ? &fCandidateDigitsVector.at(candidateIndex) :0;
670}