f3a3e179e853caa4939ae2fa2bac04e133ed040e
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCClusterFinderComponent.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  *                  Timm Steinbeck <timm@kip.uni-heidelberg.de>           *
9  *                  Jochen Thaeder <thaeder@kip.uni-heidelberg.de>        *
10  *                  for The ALICE HLT Project.                            *
11  *                                                                        *
12  * Permission to use, copy, modify and distribute this software and its   *
13  * documentation strictly for non-commercial purposes is hereby granted   *
14  * without fee, provided that the above copyright notice appears in all   *
15  * copies and that both the copyright notice and this permission notice   *
16  * appear in the supporting documentation. The authors make no claims     *
17  * about the suitability of this software for any purpose. It is          *
18  * provided "as is" without express or implied warranty.                  *
19  **************************************************************************/
20
21 /** @file   AliHLTTPCClusterFinderComponent.cxx
22     @author Timm Steinbeck, Matthias Richter, Jochen Thaeder, Kenneth Aamodt
23     @date   
24     @brief  The TPC cluster finder processing component
25 */
26
27 // see header file for class documentation                                   //
28 // or                                                                        //
29 // refer to README to build package                                          //
30 // or                                                                        //
31 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt                          //
32
33 #if __GNUC__>= 3
34 using namespace std;
35 #endif
36 #include "AliHLTTPCClusterFinderComponent.h"
37 #include "AliHLTTPCDigitReaderPacked.h"
38 #include "AliHLTTPCDigitReaderUnpacked.h"
39 #include "AliHLTTPCDigitReaderRaw.h"
40 #include "AliHLTTPCClusterFinder.h"
41 #include "AliHLTTPCSpacePointData.h"
42 #include "AliHLTTPCClusterDataFormat.h"
43 #include "AliHLTTPCTransform.h"
44 #include "AliHLTTPCClusters.h"
45 #include "AliHLTTPCDefinitions.h"
46 #include <cstdlib>
47 #include <cerrno>
48 #include "TString.h"
49 #include <sys/time.h>
50
51 // this is a global object used for automatic component registration, do not use this
52 // use fPackedSwitch = true for packed inputtype "gkDDLPackedRawDataType"
53 // use fPackedSwitch = false for unpacked inputtype "gkUnpackedRawDataType"
54 AliHLTTPCClusterFinderComponent gAliHLTTPCClusterFinderComponentPacked(true);
55 AliHLTTPCClusterFinderComponent gAliHLTTPCClusterFinderComponentUnpacked(false);
56
57 ClassImp(AliHLTTPCClusterFinderComponent)
58
59 AliHLTTPCClusterFinderComponent::AliHLTTPCClusterFinderComponent(bool packed)
60   :
61   fClusterFinder(NULL),
62   fReader(NULL),
63   fClusterDeconv(true),
64   fXYClusterError(-1),
65   fZClusterError(-1),
66   fPackedSwitch(packed),
67   fUnsorted(0),
68   fPatch(0),
69   fPadArray(NULL),
70   fGetActivePads(0)
71 {
72   // see header file for class documentation
73   // or
74   // refer to README to build package
75   // or
76   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
77 }
78
79 AliHLTTPCClusterFinderComponent::~AliHLTTPCClusterFinderComponent()
80 {
81   // see header file for class documentation
82 }
83
84 // Public functions to implement AliHLTComponent's interface.
85 // These functions are required for the registration process
86
87 const char* AliHLTTPCClusterFinderComponent::GetComponentID()
88 {
89   // see header file for class documentation
90   if (fPackedSwitch) return "TPCClusterFinderPacked";
91   else return "TPCClusterFinderUnpacked";
92 }
93
94 void AliHLTTPCClusterFinderComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list)
95 {
96   // see header file for class documentation
97   list.clear(); 
98   if (fPackedSwitch) list.push_back( kAliHLTDataTypeDDLRaw | kAliHLTDataOriginTPC );
99   else list.push_back( AliHLTTPCDefinitions::fgkUnpackedRawDataType );
100    
101 }
102
103 AliHLTComponentDataType AliHLTTPCClusterFinderComponent::GetOutputDataType()
104 {
105   // see header file for class documentation
106   return kAliHLTMultipleDataType;
107 }
108
109 int AliHLTTPCClusterFinderComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)
110
111 {
112   // see header file for class documentation
113   tgtList.clear();
114   tgtList.push_back(AliHLTTPCDefinitions::fgkClustersDataType);
115   tgtList.push_back(AliHLTTPCDefinitions::fgkActivePadsDataType);
116   return tgtList.size();
117 }
118
119 void AliHLTTPCClusterFinderComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
120 {
121   // see header file for class documentation
122   // XXX TODO: Find more realistic values.  
123   constBase = 0;
124   if (fPackedSwitch)  inputMultiplier = (6 * 0.4);
125   else  inputMultiplier = 0.4;
126 }
127
128 AliHLTComponent* AliHLTTPCClusterFinderComponent::Spawn()
129 {
130   // see header file for class documentation
131   return new AliHLTTPCClusterFinderComponent(fPackedSwitch);
132 }
133         
134 int AliHLTTPCClusterFinderComponent::DoInit( int argc, const char** argv )
135 {
136   // see header file for class documentation
137   if ( fClusterFinder )
138     return EINPROGRESS;
139
140   fClusterFinder = new AliHLTTPCClusterFinder();
141
142   Int_t rawreadermode =  -1;
143   Int_t sigthresh = -1;
144   Double_t sigmathresh= -1;
145   Float_t occulimit = 1.0;
146   Int_t oldRCUFormat=0;
147   // Data Format version numbers:
148   // 0: RCU Data format as delivered during TPC commissioning, pads/padrows are sorted, RCU trailer is one 32 bit word.
149   // 1: As 0, but pads/padrows are delivered "as is", without sorting
150   // 2: As 0, but RCU trailer is 3 32 bit words.
151   // 3: As 1, but RCU trailer is 3 32 bit words.
152   // -1: use offline raw reader
153
154   Int_t i = 0;
155   Char_t* cpErr;
156
157   while ( i < argc ) {      
158
159     // -- raw reader mode option
160     if ( !strcmp( argv[i], "rawreadermode" ) ) {
161       if ( argc <= i+1 ) {
162         Logging( kHLTLogError, "HLT::TPCClusterFinder::DoInit", "Missing rawreadermode", "Raw Reader Mode not specified" );
163         return ENOTSUP;
164       }
165
166       // Decodes the rawreader mode: either number or string and returns the rawreadermode
167       // -1 on failure, -2 for offline
168       rawreadermode = AliHLTTPCDigitReaderRaw::DecodeMode( argv[i+1] );
169
170       if (rawreadermode == -1 ) {
171         Logging( kHLTLogError, "HLT::TPCClusterFinder::DoInit", "Missing rawreadermode", "Cannot convert rawreadermode specifier '%s'.", argv[i+1] );
172         return EINVAL;
173       }
174
175       i += 2;
176       continue;
177     }
178
179     // -- pp run option
180     if ( !strcmp( argv[i], "pp-run" ) ) {
181       fClusterDeconv = false;
182       i++;
183       continue;
184     }
185
186     // -- zero suppression threshold
187     if ( !strcmp( argv[i], "adc-threshold" ) ) {
188       sigthresh = strtoul( argv[i+1], &cpErr ,0);
189       if ( *cpErr ) {
190         HLTError("Cannot convert threshold specifier '%s'.", argv[i+1]);
191         return EINVAL;
192       }
193       i+=2;
194       continue;
195     }
196
197     // -- pad occupancy limit
198     if ( !strcmp( argv[i], "occupancy-limit" ) ) {
199       occulimit = strtod( argv[i+1], &cpErr);
200       if ( *cpErr ) {
201         HLTError("Cannot convert occupancy specifier '%s'.", argv[i+1]);
202         return EINVAL;
203       }
204       i+=2;
205       continue;
206     }
207
208     // -- number of timebins (default 1024)
209     if ( !strcmp( argv[i], "timebins" ) ) {
210       TString parameter(argv[i+1]);
211       parameter.Remove(TString::kLeading, ' '); // remove all blanks
212       if (parameter.IsDigit()) {
213         AliHLTTPCTransform::SetNTimeBins(parameter.Atoi());
214         HLTInfo("number of timebins set to %d, zbin=%f", AliHLTTPCTransform::GetNTimeBins(), AliHLTTPCTransform::GetZWidth());
215       } else {
216         HLTError("Cannot timebin specifier '%s'.", argv[i+1]);
217         return EINVAL;
218       }
219       i+=2;
220       continue;
221     }
222
223     // -- checking for rcu format
224     if ( !strcmp( argv[i], "oldrcuformat" ) ) {
225       oldRCUFormat = strtoul( argv[i+1], &cpErr ,0);
226       if ( *cpErr ){
227         HLTError("Cannot convert oldrcuformat specifier '%s'. Should  be 0(off) or 1(on), must be integer", argv[i+1]);
228         return EINVAL;
229       }
230       i+=2;
231       continue;
232     }
233       
234     // -- checking for unsorted clusterfinding
235     if ( !strcmp( argv[i], "unsorted" ) ) {
236       fUnsorted = strtoul( argv[i+1], &cpErr ,0);
237       if ( *cpErr ){
238         HLTError("Cannot convert unsorted specifier '%s'. Should  be 0(off) or 1(on), must be integer", argv[i+1]);
239         return EINVAL;
240       }
241       i+=2;
242       continue;
243     }
244       
245     // -- checking for active pads, used in 2007 December run
246     if ( !strcmp( argv[i], "activepads" ) ) {
247       fGetActivePads = strtoul( argv[i+1], &cpErr ,0);
248       if ( *cpErr ){
249         HLTError("Cannot convert activepads specifier '%s'. Should  be 0(off) or 1(on), must be integer", argv[i+1]);
250         return EINVAL;
251       }
252       i+=2;
253       continue;
254     }
255
256     // -- checking for nsigma-threshold, used in 2007 December run in ZeroSuppression
257     if ( !strcmp( argv[i], "nsigma-threshold" ) ) {
258        sigmathresh = strtoul( argv[i+1], &cpErr ,0);
259       if ( *cpErr ){
260         HLTError("Cannot convert nsigma-threshold specifier '%s'. Must be integer", argv[i+1]);
261         return EINVAL;
262       }
263       i+=2;
264       continue;
265     }
266
267     Logging(kHLTLogError, "HLT::TPCClusterFinder::DoInit", "Unknown Option", "Unknown option '%s'", argv[i] );
268     return EINVAL;
269
270   }
271
272   // Choose reader
273
274   if (fPackedSwitch) { 
275     if (rawreadermode == -2) {
276       HLTDebug("using AliHLTTPCDigitReaderPacked");
277       fReader = new AliHLTTPCDigitReaderPacked();
278       if(oldRCUFormat==1){
279         fReader->SetOldRCUFormat(kTRUE);
280       }
281       else if(oldRCUFormat!=0){
282         HLTWarning("Wrong oldrcuformat specifier %d; oldrcuformat set to default(kFALSE)",oldRCUFormat);
283       }
284       if(fUnsorted==1){
285         fReader->SetUnsorted(kTRUE);
286       }
287       fClusterFinder->SetReader(fReader);
288     } else {
289 #if defined(HAVE_TPC_MAPPING)
290       HLTDebug("using AliHLTTPCDigitReaderRaw mode %d", rawreadermode);
291       fReader = new AliHLTTPCDigitReaderRaw(rawreadermode);
292       fClusterFinder->SetReader(fReader);
293 #else //! defined(HAVE_TPC_MAPPING)
294       HLTFatal("DigitReaderRaw not available - check your build");
295       return -ENODEV;
296 #endif //defined(HAVE_TPC_MAPPING)
297     }
298   }
299   else {
300     HLTDebug("using AliHLTTPCDigitReaderUnpacked");
301     fReader = new AliHLTTPCDigitReaderUnpacked();
302     fClusterFinder->SetReader(fReader);
303   }
304
305   // if pp-run use occupancy limit else set to 1. ==> use all 
306   if ( !fClusterDeconv )
307     fClusterFinder->SetOccupancyLimit(occulimit);
308   else 
309     fClusterFinder->SetOccupancyLimit(1.0);
310       
311   // Variables to setup the Clusterfinder
312   // TODO: this sounds strange and has to be verified; is the cluster finder not working when
313   // fClusterDeconv = false ?
314   fClusterDeconv = true;
315   fXYClusterError = -1;
316   fZClusterError = -1;
317
318  
319   fClusterFinder->SetDeconv( fClusterDeconv );
320   fClusterFinder->SetXYError( fXYClusterError );
321   fClusterFinder->SetZError( fZClusterError );
322   if ( (fXYClusterError>0) && (fZClusterError>0) )
323     fClusterFinder->SetCalcErr( false );
324   fClusterFinder->SetSignalThreshold(sigthresh);
325   fClusterFinder->SetNSigmaThreshold(sigmathresh);
326
327   return 0;
328 }
329
330 int AliHLTTPCClusterFinderComponent::DoDeinit()
331 {
332   // see header file for class documentation
333
334   if ( fClusterFinder )
335     delete fClusterFinder;
336   fClusterFinder = NULL;
337  
338   if ( fReader )
339     delete fReader;
340   fReader = NULL;
341     
342   return 0;
343 }
344
345 int AliHLTTPCClusterFinderComponent::DoEvent( const AliHLTComponentEventData& evtData, 
346                                               const AliHLTComponentBlockData* blocks, 
347                                               AliHLTComponentTriggerData& /*trigData*/, AliHLTUInt8_t* outputPtr, 
348                                               AliHLTUInt32_t& size, 
349                                               vector<AliHLTComponentBlockData>& outputBlocks )
350 {
351   // see header file for class documentation
352
353   //  == init iter (pointer to datablock)
354   const AliHLTComponentBlockData* iter = NULL;
355   unsigned long ndx;
356
357   //  == OUTdatatype pointer
358   AliHLTTPCClusterData* outPtr;
359
360   AliHLTUInt8_t* outBPtr;
361   UInt_t offset, mysize, nSize, tSize = 0;
362
363   outBPtr = outputPtr;
364   outPtr = (AliHLTTPCClusterData*)outBPtr;
365
366   Int_t slice, patch, row[2];
367   unsigned long maxPoints, realPoints = 0;
368
369   for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ )
370     {
371       iter = blocks+ndx;
372       mysize = 0;
373       offset = tSize;
374
375
376       if (fPackedSwitch) {
377         HLTDebug("Event 0x%08LX (%Lu) received datatype: %s - required datatype: %s",
378                  evtData.fEventID, evtData.fEventID, 
379                  DataType2Text( iter->fDataType).c_str(), 
380                  DataType2Text(kAliHLTDataTypeDDLRaw | kAliHLTDataOriginTPC).c_str());
381
382         if (iter->fDataType == AliHLTTPCDefinitions::fgkDDLPackedRawDataType &&
383             GetEventCount()<2) {
384           HLTWarning("data type %s is depricated, use %s (kAliHLTDataTypeDDLRaw)!",
385                      DataType2Text(AliHLTTPCDefinitions::fgkDDLPackedRawDataType).c_str(),
386                      DataType2Text(kAliHLTDataTypeDDLRaw | kAliHLTDataOriginTPC).c_str());
387           }
388
389         if ( iter->fDataType != (kAliHLTDataTypeDDLRaw | kAliHLTDataOriginTPC) &&
390              iter->fDataType != AliHLTTPCDefinitions::fgkDDLPackedRawDataType ) continue;
391
392       }
393       else {
394         HLTDebug("Event 0x%08LX (%Lu) received datatype: %s - required datatype: %s",
395                  evtData.fEventID, evtData.fEventID, 
396                  DataType2Text( iter->fDataType).c_str(), 
397                  DataType2Text(AliHLTTPCDefinitions::fgkUnpackedRawDataType).c_str());
398
399         if ( iter->fDataType != AliHLTTPCDefinitions::fgkUnpackedRawDataType ) continue;
400
401       }
402         
403       slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
404       patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
405       row[0] = AliHLTTPCTransform::GetFirstRow( patch );
406       row[1] = AliHLTTPCTransform::GetLastRow( patch );
407         
408       if(fUnsorted){
409         if(fPadArray==NULL){
410           fClusterFinder->SetUnsorted(fUnsorted);
411           fPadArray = new AliHLTTPCPadArray(patch);
412           fPadArray->InitializeVector();
413         }
414         else if(fPadArray->GetPatch()!=patch||fPadArray->GetPatch()==-1){
415           if (GetEventCount()<3) {
416             HLTWarning("pad array not initialized for data of specification 0x%08x, block skipped", iter->fSpecification);
417           } else if ((GetEventCount()%5000)==0) { // assuming 0.5 to 1kHz this gives a message rate of 0.1 to 0.5 Hz
418             HLTWarning("reminder: pad array not initialized for data of specification 0x%08x", iter->fSpecification);
419           }
420           continue;
421         }
422       }
423
424       outPtr = (AliHLTTPCClusterData*)outBPtr;
425
426       maxPoints = (size-tSize-sizeof(AliHLTTPCClusterData))/sizeof(AliHLTTPCSpacePointData);
427
428       fClusterFinder->InitSlice( slice, patch, row[0], row[1], maxPoints );
429       fClusterFinder->SetOutputArray( (AliHLTTPCSpacePointData*)outPtr->fSpacePoints );
430         
431       if(fUnsorted){
432
433
434         fClusterFinder->SetPadArray(fPadArray);
435         /*        
436         double totalT=0;
437         struct timeval startT, endT;
438         gettimeofday( &startT, NULL );
439         */
440         fClusterFinder->ReadDataUnsorted(iter->fPtr, iter->fSize );
441         /*
442         gettimeofday( &endT, NULL );
443         unsigned long long dt;
444         dt = endT.tv_sec-startT.tv_sec;
445         dt *= 1000000ULL;
446         dt += endT.tv_usec-startT.tv_usec;
447         double dtd = ((double)dt);
448         totalT += dtd;
449         //        dtd = dtd / (double)eventIterations;
450         //        if ( iterations<=1 )
451         cout<<endl;
452         printf( "Time needed to read data: %f microsec. / %f millisec. / %f s\n", 
453                 dtd, dtd/1000.0, dtd/1000000.0 );
454           
455         cout<<endl;
456         */
457         fClusterFinder->FindClusters();
458       }
459       else{
460         fClusterFinder->Read(iter->fPtr, iter->fSize );
461         fClusterFinder->ProcessDigits();
462       }
463       realPoints = fClusterFinder->GetNumberOfClusters();
464         
465       outPtr->fSpacePointCnt = realPoints;
466       nSize = sizeof(AliHLTTPCSpacePointData)*realPoints;
467       mysize += nSize+sizeof(AliHLTTPCClusterData);
468
469       Logging( kHLTLogDebug, "HLT::TPCClusterFinder::DoEvent", "Spacepoints", 
470                "Number of spacepoints: %lu Slice/Patch/RowMin/RowMax: %d/%d/%d/%d.",
471                realPoints, slice, patch, row[0], row[1] );
472       AliHLTComponentBlockData bd;
473       FillBlockData( bd );
474       bd.fOffset = offset;
475       bd.fSize = mysize;
476       bd.fSpecification = iter->fSpecification;
477       bd.fDataType = AliHLTTPCDefinitions::fgkClustersDataType;
478       //AliHLTSubEventDescriptor::FillBlockAttributes( bd.fAttributes );
479       outputBlocks.push_back( bd );
480         
481       tSize += mysize;
482       outBPtr += mysize;
483       outPtr = (AliHLTTPCClusterData*)outBPtr;
484         
485       if(fGetActivePads){
486         AliHLTTPCPadArray::AliHLTTPCActivePads* outPtrActive;
487         UInt_t activePadsSize, activePadsN = 0;
488         outPtrActive = (AliHLTTPCPadArray::AliHLTTPCActivePads*)outBPtr;
489         offset=tSize;
490         Int_t maxActivePads = (size-tSize)/sizeof(AliHLTTPCPadArray::AliHLTTPCActivePads);
491         activePadsSize= fClusterFinder->GetActivePads((AliHLTTPCPadArray::AliHLTTPCActivePads*)outPtrActive,maxActivePads)*sizeof(AliHLTTPCPadArray::AliHLTTPCActivePads);
492         
493         AliHLTComponentBlockData bdActive;
494         FillBlockData( bdActive );
495         bdActive.fOffset = offset;
496         bdActive.fSize = activePadsSize;
497         bdActive.fSpecification = iter->fSpecification;
498         bdActive.fDataType = AliHLTTPCDefinitions::fgkActivePadsDataType;
499         outputBlocks.push_back( bdActive );
500         
501         tSize+=activePadsSize;
502         outBPtr += activePadsSize;
503         outPtrActive = (AliHLTTPCPadArray::AliHLTTPCActivePads*)outBPtr;
504       }
505  
506
507       if ( tSize > size )
508         {
509           Logging( kHLTLogFatal, "HLT::TPCClusterFinder::DoEvent", "Too much data", 
510                    "Data written over allowed buffer. Amount written: %lu, allowed amount: %lu.",
511                    tSize, size );
512           return EMSGSIZE;
513         }
514     }
515     
516   size = tSize;
517
518   return 0;
519 }