0d49c603c7c03a43957d2300b85f273ca711efe7
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCPad.cxx
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: Matthias Richter <Matthias.Richter@ift.uib.no>        *
8  *                  Kenneth Aamodt   <Kenneth.aamodt@ift.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 #if __GNUC__>= 3
27 using namespace std;
28 #endif
29
30 #include <cerrno>
31 #include "AliHLTTPCPad.h"
32 #include "AliHLTStdIncludes.h"
33
34
35 //added by kenneth
36 #include "AliHLTTPCTransform.h"
37 #include "AliHLTTPCClusters.h"
38 #include <sys/time.h>
39 //------------------------------
40
41 /** margin for the base line be re-avaluated */
42 #define ALIHLTPAD_BASELINE_MARGIN (2*fAverage)
43
44 /** ROOT macro for the implementation of ROOT specific class methods */
45 ClassImp(AliHLTTPCPad)
46
47 AliHLTTPCPad::AliHLTTPCPad()
48   :
49   fRowNo(-1),
50   fPadNo(-1),
51   fThreshold(0),
52   fAverage(-1),
53   fNofEvents(0),
54   fSum(0),
55   fCount(0),
56   fTotal(0),
57   fBLMax(-1),
58   fBLMaxBin(-1),
59   fBLMin(-1),
60   fBLMinBin(-1),
61   fFirstBLBin(0),
62   fNofBins(0),
63   fReadPos(0),
64   fpRawData(NULL),
65   fClusterCandidates(),
66   fUsedClusterCandidates(0),
67   fDataSignals(NULL),
68   fSignalPositionArray(NULL),
69   fSizeOfSignalPositionArray(0)
70 {
71   // see header file for class documentation
72   // or
73   // refer to README to build package
74   // or
75   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
76   //  HLTInfo("Entering default constructor");
77   fDataSignals= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
78   memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
79
80   fSignalPositionArray= new AliHLTTPCSignal_t[AliHLTTPCTransform::GetNTimeBins()];
81   memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
82   fSizeOfSignalPositionArray=0;
83 }
84
85 AliHLTTPCPad::AliHLTTPCPad(Int_t offset, Int_t nofBins)
86   :
87   fRowNo(-1),
88   fPadNo(-1),
89   fThreshold(0),
90   fAverage(-1),
91   fNofEvents(0),
92   fSum(0),
93   fCount(0),
94   fTotal(0),
95   fBLMax(-1),
96   fBLMaxBin(-1),
97   fBLMin(-1),
98   fBLMinBin(-1),
99   fFirstBLBin(offset),
100   fNofBins(nofBins),
101   fReadPos(0),
102   fpRawData(NULL),
103   fClusterCandidates(),
104   fUsedClusterCandidates(0),
105   fDataSignals(NULL),
106   fSignalPositionArray(NULL),
107   fSizeOfSignalPositionArray(0)
108 {
109   // see header file for class documentation
110 }
111
112 AliHLTTPCPad::AliHLTTPCPad(const AliHLTTPCPad& srcPad)
113   :
114   fRowNo(srcPad.fRowNo),
115   fPadNo(srcPad.fPadNo),
116   fThreshold(0),
117   fAverage(-1),
118   fNofEvents(0),
119   fSum(0),
120   fCount(0),
121   fTotal(0),
122   fBLMax(-1),
123   fBLMaxBin(-1),
124   fBLMin(-1),
125   fBLMinBin(-1),
126   fFirstBLBin(0),
127   fNofBins(0),
128   fReadPos(0),
129   fpRawData(NULL),
130   fClusterCandidates(),
131   fUsedClusterCandidates(0),
132   fDataSignals(NULL),
133   fSignalPositionArray(NULL),
134   fSizeOfSignalPositionArray(0)
135 {
136   // see header file for class documentation
137   HLTFatal("copy constructor not implemented");
138 }
139
140 AliHLTTPCPad& AliHLTTPCPad::operator=(const AliHLTTPCPad&)
141 {
142   // see header file for class documentation
143   HLTFatal("assignment operator not implemented");
144   return (*this);
145 }
146
147 AliHLTTPCPad::~AliHLTTPCPad()
148 {
149   // see header file for class documentation
150   if (fpRawData) {
151     HLTWarning("event data acquisition not stopped");
152     StopEvent();
153   }
154   if (fDataSignals) {
155     AliHLTTPCSignal_t* pData=fDataSignals;
156     fDataSignals=NULL;
157    delete [] pData;
158   }
159   if (fSignalPositionArray) {
160     AliHLTTPCSignal_t* pData=fSignalPositionArray;
161     fSignalPositionArray=NULL;
162     //   delete [] pData;
163   }
164
165 }
166
167 Int_t AliHLTTPCPad::SetID(Int_t rowno, Int_t padno)
168 {
169   // see header file for class documentation
170   fRowNo=rowno;
171   fPadNo=padno;
172   return 0;
173 }
174
175 Int_t AliHLTTPCPad::StartEvent()
176 {
177   // see header file for class documentation
178   Int_t iResult=0;
179   if (fpRawData==NULL) {
180     fBLMax=-1;
181     fBLMaxBin=-1;
182     fBLMin=-1;
183     fBLMinBin=-1;
184     fSum=0;
185     fCount=0;
186     fTotal=0;
187     if (fNofBins>0) {
188       fpRawData=new AliHLTTPCSignal_t[fNofBins];
189       if (fpRawData) {
190         for (int i=0; i<fNofBins; i++) fpRawData[i]=-1;
191       } else {
192         HLTError("memory allocation failed");
193         iResult=-ENOMEM;
194       }
195     }
196   } else {
197     HLTWarning("event data acquisition already started");
198     iResult=-EALREADY;
199   }
200   return iResult;
201 }
202
203 Int_t AliHLTTPCPad::CalculateBaseLine(Int_t reqMinCount)
204 {
205   // see header file for class documentation
206   Int_t iResult=0;
207   AliHLTTPCSignal_t avBackup=fAverage;
208   //HLTDebug("reqMinCount=%d fCount=%d fTotal=%d fSum=%d fBLMax=%d fBLMin=%d", reqMinCount, fCount, fTotal, fSum, fBLMax, fBLMin);
209   if (fCount>=reqMinCount && fCount>=fTotal/2) {
210     fAverage=fCount>0?fSum/fCount:0;
211     if (fAverage>0) {
212       //HLTDebug("average for current event %d (%d - %d)", fAverage, fBLMax, fBLMin);
213       fCount=0;fSum=-1;
214       if (fBLMax>ALIHLTPAD_BASELINE_MARGIN) {
215         // calculate again
216         //HLTDebug("maximum value %d exceeds margin for base line (%d) "
217         //       "-> re-evaluate base line", fBLMax, ALIHLTPAD_BASELINE_MARGIN);
218         if (fpRawData) {
219           for (Int_t i=fFirstBLBin; i<fNofBins; i++)
220             if (fpRawData[i]>=0) AddBaseLineValue(i, fpRawData[i]);
221           if (fCount>0 && fCount>=reqMinCount && fCount>=fTotal/2) {
222             fAverage=fSum/fCount;
223             //HLTDebug("new average %d", fAverage);
224           } else {
225 //          HLTDebug("baseline re-eveluation skipped because of to few "
226 //                     "contributing bins: total=%d, contributing=%d, req=%d"
227 //                     "\ndata might be already zero suppressed"
228 //                     , fTotal, fCount, reqMinCount);
229             iResult=-ENODATA;
230           }
231           fCount=0;fSum=-1;
232         } else {
233           HLTError("missing raw data for base line calculation");
234           iResult=-ENOBUFS;
235         }
236       }
237       if (iResult>=0) {
238         // calculate average for all events
239         fAverage=((avBackup*fNofEvents)+fAverage)/(fNofEvents+1);
240         //HLTDebug("base line average for %d event(s): %d", fNofEvents+1, fAverage);
241       } else {
242         fAverage=avBackup;      
243       }
244     } else {
245       fAverage=avBackup;
246     }
247   } else {
248 //     HLTDebug("baseline calculation skipped because of to few contributing "
249 //             "bins: total=%d, contributing=%d, required=%d \ndata might be "
250 //             "already zero suppressed", fTotal, fCount, reqMinCount);
251   }
252
253   return iResult;
254 }
255
256 Int_t AliHLTTPCPad::StopEvent()
257 {
258   // see header file for class documentation
259   Int_t iResult=0;
260   if (fpRawData) {
261     AliHLTTPCSignal_t* pData=fpRawData;
262     fpRawData=NULL;
263     delete [] pData;
264     fTotal=0;
265     fNofEvents++;
266     Rewind();
267   } else if (fNofBins>0) {
268     HLTError("event data acquisition not started");
269     iResult=-EBADF;
270   }
271   return iResult;
272 }
273
274 Int_t AliHLTTPCPad::ResetHistory()
275 {
276   // see header file for class documentation
277   Int_t iResult=0;
278   fAverage=-1;
279   fNofEvents=0;
280   return iResult;
281 }
282
283 Int_t AliHLTTPCPad::SetThreshold(AliHLTTPCSignal_t thresh)
284 {
285   // see header file for class documentation
286   Int_t iResult=0;
287   fThreshold=thresh;
288   return iResult;
289 }
290
291 Int_t AliHLTTPCPad::AddBaseLineValue(Int_t bin, AliHLTTPCSignal_t value)
292 {
293   // see header file for class documentation
294   Int_t iResult=0;
295   if (bin>=fFirstBLBin) {
296     if (fAverage<0 || value<ALIHLTPAD_BASELINE_MARGIN) {
297       // add to the current sum and count
298       fSum+=value;
299       fCount++;
300       if (fBLMax<value) {
301         // keep the maximum value for later quality control of the base 
302         // line calculation
303         fBLMax=value;
304         fBLMaxBin=bin;
305       }
306       if (fBLMin<0 || fBLMin>value) {
307         // keep the minimum value for later quality control of the base 
308         // line calculation
309         fBLMin=value;
310         fBLMinBin=bin;
311       }
312     } else {
313 //       HLTDebug("ignoring value %d (bin %d) for base line calculation "
314 //             "(current average is %d)",
315 //             value, bin, fAverage);
316     }
317   }
318   return iResult;
319 }
320
321 Int_t AliHLTTPCPad::SetRawData(Int_t bin, AliHLTTPCSignal_t value)
322 {
323   // see header file for class documentation
324   Int_t iResult=0;
325   if (fpRawData) {
326     if (bin<fNofBins) {
327       if (value>=0) {
328         if (fpRawData[bin]<0) {
329           AddBaseLineValue(bin, value);
330           fTotal++;
331         } else {
332           // ignore value for average calculation
333           HLTWarning("overriding content of bin %d (%d)", bin, fpRawData[bin]);
334         }
335         fpRawData[bin]=value;
336       } else {
337         HLTWarning("ignoring neg. raw data");
338       }
339     } else {
340       HLTWarning("bin %d out of range (%d)", bin, fNofBins);
341       iResult=-ERANGE;
342     }
343   } else if (fNofBins>0) {
344     HLTError("event cycle not started");
345     iResult=-EBADF;
346   }
347   return iResult;
348 }
349
350 Int_t AliHLTTPCPad::Next(Int_t bZeroSuppression) 
351 {
352   // see header file for class documentation
353   if (fpRawData==NULL) return 0;
354   Int_t iResult=fReadPos<fNofBins;
355   if (iResult>0 && (iResult=(++fReadPos<fNofBins))>0) {
356     if (bZeroSuppression) {
357       while ((iResult=(fReadPos<fNofBins))>0 &&
358              GetCorrectedData(fReadPos)<=0)
359         fReadPos++;
360     }
361   }
362   return iResult;
363 }
364
365 Int_t AliHLTTPCPad::Rewind(Int_t bZeroSuppression)
366 {
367   // see header file for class documentation
368   fReadPos=(bZeroSuppression>0?0:fFirstBLBin)-1;
369   return Next(bZeroSuppression);
370 }
371
372 AliHLTTPCSignal_t AliHLTTPCPad::GetRawData(Int_t bin) const
373 {
374   // see header file for class documentation
375   AliHLTTPCSignal_t data=0;
376   if (fpRawData) {
377     if (bin<fNofBins) {
378       data=fpRawData[bin];
379     } else {
380       HLTWarning("requested bin %d out of range (%d)", bin, fNofBins);
381     }
382   } else if (fNofBins>0) {
383     HLTWarning("data only available within event cycle");
384   }
385   return data;
386 }
387
388 AliHLTTPCSignal_t AliHLTTPCPad::GetCorrectedData(Int_t bin) const
389 {
390   // see header file for class documentation
391   AliHLTTPCSignal_t data=GetRawData(bin)-GetBaseLine(bin);
392   AliHLTTPCSignal_t prev=0;
393   if (bin>1) prev=GetRawData(bin-1)-GetBaseLine(bin-1);
394   AliHLTTPCSignal_t succ=0;
395   if (bin+1<GetSize()) succ=GetRawData(bin+1)-GetBaseLine(bin+1);
396   if (fThreshold>0) {
397     data-=fThreshold;
398     prev-=fThreshold;
399     succ-=fThreshold;
400   }
401  
402   // case 1:
403   // the signal is below the base-line and threshold
404   if (data<0) data=0;
405
406   //case 2:
407   // the neighboring bins are both below base-line/threshold
408   // a real signal is always more than one bin wide because of the shaper 
409   if (prev<=0 && succ<=0) data=0;
410  
411   // case 3:
412   // the bin is inside the range of ignored bins
413   if (bin<fFirstBLBin) data=0;
414   //HLTDebug("fReadPos=%d data=%d threshold=%d raw data=%d base line=%d", fReadPos, data, fThreshold, GetRawData(bin), GetBaseLine(bin));
415   return data;
416 }
417
418 AliHLTTPCSignal_t AliHLTTPCPad::GetBaseLine(Int_t bin) const
419 {
420   // see header file for class documentation
421   AliHLTTPCSignal_t val=0;
422   if (fAverage>0) {
423     // we take the minumum value as the base line if it doesn't differ from
424     // the average to much
425     val=fAverage;
426 #ifdef KEEP_NOISE
427     const AliHLTTPCSignal_t kMaxDifference=15;
428     if ((fAverage-fBLMin)<=kMaxDifference) val=fBLMin;
429     else val>kMaxDifference?val-=kMaxDifference:0;
430 #endif
431   }
432   if (val<0) {
433     // here we should never get
434     val=0;
435     HLTFatal("wrong base line value");
436   }
437   return val;
438 }
439
440 AliHLTTPCSignal_t AliHLTTPCPad::GetAverage() const
441 {
442   // see header file for class documentation
443   return fAverage>0?fAverage:0;
444 }
445
446 Float_t AliHLTTPCPad::GetOccupancy() const
447 {
448   // see header file for class documentation
449   Float_t occupancy=0;
450   if (fpRawData && fNofBins>0) {
451     for (Int_t i=fFirstBLBin; i<fNofBins; i++) {
452       if (GetCorrectedData(i)>0) occupancy+=1;
453     }
454     if (fNofBins-fFirstBLBin>0)
455       occupancy/=fNofBins-fFirstBLBin;
456   }
457   return occupancy;
458 }
459
460 Float_t AliHLTTPCPad::GetAveragedOccupancy() const
461 {
462   // see header file for class documentation
463
464   // history is not yet implemented
465   return GetOccupancy();
466 }
467 void AliHLTTPCPad::PrintRawData(){
468   for(Int_t bin=0;bin<AliHLTTPCTransform::GetNTimeBins();bin++){
469     if(GetDataSignal(bin)>0)
470       cout<<fRowNo<<"\t"<<fPadNo<<"\t"<<bin<<"\t"<<GetDataSignal(bin)<<endl;;
471   }
472   cout<<"bins: "<<AliHLTTPCTransform::GetNTimeBins()<<endl;
473 }
474
475 void AliHLTTPCPad::SetDataToDefault(){
476   if(fpRawData){
477     memset( fDataSignals, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
478     memset( fSignalPositionArray, 0xFF, sizeof(Int_t)*(AliHLTTPCTransform::GetNTimeBins()));
479     fSizeOfSignalPositionArray=0;
480   }
481 }
482 void AliHLTTPCPad::SetDataSignal(Int_t bin,Int_t signal){
483   fDataSignals[bin]=signal;
484   fSignalPositionArray[fSizeOfSignalPositionArray]=bin;
485   fSizeOfSignalPositionArray++;
486 }
487 Int_t AliHLTTPCPad::GetDataSignal(Int_t bin){
488   return fDataSignals[bin];
489 }
490 void AliHLTTPCPad::FindClusterCandidates(){
491   UInt_t seqcharge=0;
492   UInt_t seqaverage=0;
493   UInt_t seqerror=0;
494   vector<Int_t> tmpPos;
495   vector<Int_t> tmpSig;
496   UInt_t isFalling=0;
497
498   for(Int_t pos=fSizeOfSignalPositionArray-2;pos>=0;pos--){
499
500     if(fSignalPositionArray[pos]==fSignalPositionArray[pos+1]+1){
501       seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];     
502       seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
503       seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
504           
505       tmpPos.push_back(fSignalPositionArray[pos+1]);
506       tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
507
508       if(fDataSignals[fSignalPositionArray[pos+1]]>fDataSignals[fSignalPositionArray[pos]]){
509         isFalling=1;
510       }
511       if(fDataSignals[fSignalPositionArray[pos+1]]<fDataSignals[fSignalPositionArray[pos]]&&isFalling){
512         Int_t seqmean=0;
513         seqmean = seqaverage/seqcharge;
514         
515         //Calculate mean in pad direction:
516         Int_t padmean = seqcharge*fPadNo;
517         Int_t paderror = fPadNo*padmean;
518         AliHLTTPCClusters candidate;
519         candidate.fTotalCharge   = seqcharge;
520         candidate.fPad       = padmean;
521         candidate.fPad2      = paderror;
522         candidate.fTime      = seqaverage;
523         candidate.fTime2     = seqerror;
524         candidate.fMean          = seqmean;
525         candidate.fLastMergedPad = fPadNo;
526         fClusterCandidates.push_back(candidate);
527         fUsedClusterCandidates.push_back(0);
528         isFalling=0;
529         seqcharge=0;
530         seqaverage=0;
531         seqerror=0;
532
533         tmpPos.clear();
534         tmpSig.clear();
535
536         continue;
537       }
538          
539       if(pos<1){
540         seqcharge+=fDataSignals[fSignalPositionArray[0]];       
541         seqaverage += fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
542         seqerror += fSignalPositionArray[0]*fSignalPositionArray[0]*fDataSignals[fSignalPositionArray[0]];
543         tmpPos.push_back(fSignalPositionArray[0]);
544         tmpSig.push_back(fDataSignals[fSignalPositionArray[0]]);
545           
546         //Calculate mean of sequence:
547         Int_t seqmean=0;
548         seqmean = seqaverage/seqcharge;
549           
550         //Calculate mean in pad direction:
551         Int_t padmean = seqcharge*fPadNo;
552         Int_t paderror = fPadNo*padmean;
553         AliHLTTPCClusters candidate;
554         candidate.fTotalCharge   = seqcharge;
555         candidate.fPad       = padmean;
556         candidate.fPad2      = paderror;
557         candidate.fTime      = seqaverage;
558         candidate.fTime2     = seqerror;
559         candidate.fMean          = seqmean;
560         candidate.fLastMergedPad = fPadNo;
561         fClusterCandidates.push_back(candidate);
562         fUsedClusterCandidates.push_back(0);
563         isFalling=0;
564         seqcharge=0;
565         seqaverage=0;
566         seqerror=0;
567
568         tmpPos.clear();
569         tmpSig.clear();
570       }
571     }
572     else if(seqcharge>0){
573       seqcharge+=fDataSignals[fSignalPositionArray[pos+1]];     
574       seqaverage += fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
575         seqerror += fSignalPositionArray[pos+1]*fSignalPositionArray[pos+1]*fDataSignals[fSignalPositionArray[pos+1]];
576         tmpPos.push_back(fSignalPositionArray[pos+1]);
577         tmpSig.push_back(fDataSignals[fSignalPositionArray[pos+1]]);
578
579         //Calculate mean of sequence:
580         Int_t seqmean=0;
581         seqmean = seqaverage/seqcharge;
582         
583         //Calculate mean in pad direction:
584         Int_t padmean = seqcharge*fPadNo;
585         Int_t paderror = fPadNo*padmean;
586         AliHLTTPCClusters candidate;
587         candidate.fTotalCharge   = seqcharge;
588         candidate.fPad       = padmean;
589         candidate.fPad2      = paderror;
590         candidate.fTime      = seqaverage;
591         candidate.fTime2     = seqerror;
592         candidate.fMean          = seqmean;
593         candidate.fLastMergedPad = fPadNo;
594         fClusterCandidates.push_back(candidate);
595         fUsedClusterCandidates.push_back(0);
596         isFalling=0;
597         seqcharge=0;
598         seqaverage=0;
599         seqerror=0;
600
601         tmpPos.clear();
602         tmpSig.clear();
603     }
604   }
605 }