commit from Kenneth:
[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( AliHLTTPCDefinitions::fgkDDLPackedRawDataType );
99   else list.push_back( AliHLTTPCDefinitions::fgkUnpackedRawDataType );
100    
101 }
102
103 AliHLTComponentDataType AliHLTTPCClusterFinderComponent::GetOutputDataType()
104 {
105   // see header file for class documentation
106   return AliHLTTPCDefinitions::fgkClustersDataType;
107 }
108
109 void AliHLTTPCClusterFinderComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
110 {
111   // see header file for class documentation
112   // XXX TODO: Find more realistic values.  
113   constBase = 0;
114   if (fPackedSwitch)  inputMultiplier = (6 * 0.4);
115   else  inputMultiplier = 0.4;
116 }
117
118 AliHLTComponent* AliHLTTPCClusterFinderComponent::Spawn()
119 {
120   // see header file for class documentation
121   return new AliHLTTPCClusterFinderComponent(fPackedSwitch);
122 }
123         
124 int AliHLTTPCClusterFinderComponent::DoInit( int argc, const char** argv )
125 {
126   // see header file for class documentation
127   if ( fClusterFinder )
128     return EINPROGRESS;
129
130   fClusterFinder = new AliHLTTPCClusterFinder();
131
132   Int_t rawreadermode =  -1;
133   Int_t sigthresh = -1;
134   Double_t sigmathresh= -1;
135   Float_t occulimit = 1.0;
136   Int_t oldRCUFormat=0;
137   // Data Format version numbers:
138   // 0: RCU Data format as delivered during TPC commissioning, pads/padrows are sorted, RCU trailer is one 32 bit word.
139   // 1: As 0, but pads/padrows are delivered "as is", without sorting
140   // 2: As 0, but RCU trailer is 3 32 bit words.
141   // 3: As 1, but RCU trailer is 3 32 bit words.
142   // -1: use offline raw reader
143
144   Int_t i = 0;
145   Char_t* cpErr;
146
147   while ( i < argc ) {      
148
149     // -- raw reader mode option
150     if ( !strcmp( argv[i], "rawreadermode" ) ) {
151       if ( argc <= i+1 ) {
152         Logging( kHLTLogError, "HLT::TPCClusterFinder::DoInit", "Missing rawreadermode", "Raw Reader Mode not specified" );
153         return ENOTSUP;
154       }
155
156       // Decodes the rawreader mode: either number or string and returns the rawreadermode
157       // -1 on failure, -2 for offline
158       rawreadermode = AliHLTTPCDigitReaderRaw::DecodeMode( argv[i+1] );
159
160       if (rawreadermode == -1 ) {
161         Logging( kHLTLogError, "HLT::TPCClusterFinder::DoInit", "Missing rawreadermode", "Cannot convert rawreadermode specifier '%s'.", argv[i+1] );
162         return EINVAL;
163       }
164
165       i += 2;
166       continue;
167     }
168
169     // -- pp run option
170     if ( !strcmp( argv[i], "pp-run" ) ) {
171       fClusterDeconv = false;
172       i++;
173       continue;
174     }
175
176     // -- zero suppression threshold
177     if ( !strcmp( argv[i], "adc-threshold" ) ) {
178       sigthresh = strtoul( argv[i+1], &cpErr ,0);
179       if ( *cpErr ) {
180         HLTError("Cannot convert threshold specifier '%s'.", argv[i+1]);
181         return EINVAL;
182       }
183       i+=2;
184       continue;
185     }
186
187     // -- pad occupancy limit
188     if ( !strcmp( argv[i], "occupancy-limit" ) ) {
189       occulimit = strtod( argv[i+1], &cpErr);
190       if ( *cpErr ) {
191         HLTError("Cannot convert occupancy specifier '%s'.", argv[i+1]);
192         return EINVAL;
193       }
194       i+=2;
195       continue;
196     }
197
198     // -- number of timebins (default 1024)
199     if ( !strcmp( argv[i], "timebins" ) ) {
200       TString parameter(argv[i+1]);
201       parameter.Remove(TString::kLeading, ' '); // remove all blanks
202       if (parameter.IsDigit()) {
203         AliHLTTPCTransform::SetNTimeBins(parameter.Atoi());
204         HLTInfo("number of timebins set to %d, zbin=%f", AliHLTTPCTransform::GetNTimeBins(), AliHLTTPCTransform::GetZWidth());
205       } else {
206         HLTError("Cannot timebin specifier '%s'.", argv[i+1]);
207         return EINVAL;
208       }
209       i+=2;
210       continue;
211     }
212
213     // -- checking for rcu format
214     if ( !strcmp( argv[i], "oldrcuformat" ) ) {
215       oldRCUFormat = strtoul( argv[i+1], &cpErr ,0);
216       if ( *cpErr ){
217         HLTError("Cannot convert oldrcuformat specifier '%s'. Should  be 0(off) or 1(on), must be integer", argv[i+1]);
218         return EINVAL;
219       }
220       i+=2;
221       continue;
222     }
223       
224     // -- checking for unsorted clusterfinding
225     if ( !strcmp( argv[i], "unsorted" ) ) {
226       fUnsorted = strtoul( argv[i+1], &cpErr ,0);
227       if ( *cpErr ){
228         HLTError("Cannot convert unsorted specifier '%s'. Should  be 0(off) or 1(on), must be integer", argv[i+1]);
229         return EINVAL;
230       }
231       i+=2;
232       continue;
233     }
234       
235     // -- checking for active pads, used in 2007 December run
236     if ( !strcmp( argv[i], "activepads" ) ) {
237       fGetActivePads = strtoul( argv[i+1], &cpErr ,0);
238       if ( *cpErr ){
239         HLTError("Cannot convert activepads specifier '%s'. Should  be 0(off) or 1(on), must be integer", argv[i+1]);
240         return EINVAL;
241       }
242       i+=2;
243       continue;
244     }
245
246     // -- checking for nsigma-threshold, used in 2007 December run in ZeroSuppression
247     if ( !strcmp( argv[i], "nsigma-threshold" ) ) {
248        sigmathresh = strtoul( argv[i+1], &cpErr ,0);
249       if ( *cpErr ){
250         HLTError("Cannot convert nsigma-threshold specifier '%s'. Must be integer", argv[i+1]);
251         return EINVAL;
252       }
253       i+=2;
254       continue;
255     }
256
257     Logging(kHLTLogError, "HLT::TPCClusterFinder::DoInit", "Unknown Option", "Unknown option '%s'", argv[i] );
258     return EINVAL;
259
260   }
261
262   // Choose reader
263
264   if (fPackedSwitch) { 
265     if (rawreadermode == -2) {
266 #if defined(HAVE_ALIRAWDATA) && defined(HAVE_ALITPCRAWSTREAM_H)
267       HLTDebug("using AliHLTTPCDigitReaderPacked");
268       fReader = new AliHLTTPCDigitReaderPacked();
269       if(oldRCUFormat==1){
270         fReader->SetOldRCUFormat(kTRUE);
271       }
272       else if(oldRCUFormat!=0){
273         HLTWarning("Wrong oldrcuformat specifier %d; oldrcuformat set to default(kFALSE)",oldRCUFormat);
274       }
275       if(fUnsorted==1){
276         fReader->SetUnsorted(kTRUE);
277       }
278       fClusterFinder->SetReader(fReader);
279 #else // ! defined(HAVE_ALIRAWDATA) && defined(HAVE_ALITPCRAWSTREAM_H)
280       HLTFatal("DigitReaderPacked not available - check your build");
281       return -ENODEV;
282 #endif //  defined(HAVE_ALIRAWDATA) && defined(HAVE_ALITPCRAWSTREAM_H)
283     } else {
284 #if defined(HAVE_TPC_MAPPING)
285       HLTDebug("using AliHLTTPCDigitReaderRaw mode %d", rawreadermode);
286       fReader = new AliHLTTPCDigitReaderRaw(rawreadermode);
287       fClusterFinder->SetReader(fReader);
288 #else //! defined(HAVE_TPC_MAPPING)
289       HLTFatal("DigitReaderRaw not available - check your build");
290       return -ENODEV;
291 #endif //defined(HAVE_TPC_MAPPING)
292     }
293   }
294   else {
295     HLTDebug("using AliHLTTPCDigitReaderUnpacked");
296     fReader = new AliHLTTPCDigitReaderUnpacked();
297     fClusterFinder->SetReader(fReader);
298   }
299
300   // if pp-run use occupancy limit else set to 1. ==> use all 
301   if ( !fClusterDeconv )
302     fClusterFinder->SetOccupancyLimit(occulimit);
303   else 
304     fClusterFinder->SetOccupancyLimit(1.0);
305       
306   // Variables to setup the Clusterfinder
307   // TODO: this sounds strange and has to be verified; is the cluster finder not working when
308   // fClusterDeconv = false ?
309   fClusterDeconv = true;
310   fXYClusterError = -1;
311   fZClusterError = -1;
312
313  
314   fClusterFinder->SetDeconv( fClusterDeconv );
315   fClusterFinder->SetXYError( fXYClusterError );
316   fClusterFinder->SetZError( fZClusterError );
317   if ( (fXYClusterError>0) && (fZClusterError>0) )
318     fClusterFinder->SetCalcErr( false );
319   fClusterFinder->SetSignalThreshold(sigthresh);
320   fClusterFinder->SetNSigmaThreshold(sigmathresh);
321
322   return 0;
323 }
324
325 int AliHLTTPCClusterFinderComponent::DoDeinit()
326 {
327   // see header file for class documentation
328
329   if ( fClusterFinder )
330     delete fClusterFinder;
331   fClusterFinder = NULL;
332  
333   if ( fReader )
334     delete fReader;
335   fReader = NULL;
336     
337   return 0;
338 }
339
340 int AliHLTTPCClusterFinderComponent::DoEvent( const AliHLTComponentEventData& evtData, 
341                                               const AliHLTComponentBlockData* blocks, 
342                                               AliHLTComponentTriggerData& /*trigData*/, AliHLTUInt8_t* outputPtr, 
343                                               AliHLTUInt32_t& size, 
344                                               vector<AliHLTComponentBlockData>& outputBlocks )
345 {
346   // see header file for class documentation
347
348   //  == init iter (pointer to datablock)
349   const AliHLTComponentBlockData* iter = NULL;
350   unsigned long ndx;
351
352   //  == OUTdatatype pointer
353   AliHLTTPCClusterData* outPtr;
354
355   AliHLTUInt8_t* outBPtr;
356   UInt_t offset, mysize, nSize, tSize = 0;
357
358   outBPtr = outputPtr;
359   outPtr = (AliHLTTPCClusterData*)outBPtr;
360
361   Int_t slice, patch, row[2];
362   unsigned long maxPoints, realPoints = 0;
363
364   for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ )
365     {
366       iter = blocks+ndx;
367       mysize = 0;
368       offset = tSize;
369
370
371       if (fPackedSwitch) {      
372         char tmp1[14], tmp2[14];
373         DataType2Text( iter->fDataType, tmp1 );
374         DataType2Text( AliHLTTPCDefinitions::fgkDDLPackedRawDataType, tmp2 );
375         Logging( kHLTLogDebug, "HLT::TPCClusterFinder::DoEvent", "Event received", 
376                  "Event 0x%08LX (%Lu) received datatype: %s - required datatype: %s",
377                  evtData.fEventID, evtData.fEventID, tmp1, tmp2 );
378
379         if ( iter->fDataType != AliHLTTPCDefinitions::fgkDDLPackedRawDataType ) continue;
380
381       }
382       else {
383         char tmp1[14], tmp2[14];
384         DataType2Text( iter->fDataType, tmp1 );
385         DataType2Text( AliHLTTPCDefinitions::fgkUnpackedRawDataType, tmp2 );
386         Logging( kHLTLogDebug, "HLT::TPCClusterFinder::DoEvent", "Event received", 
387                  "Event 0x%08LX (%Lu) received datatype: %s - required datatype: %s",
388                  evtData.fEventID, evtData.fEventID, tmp1, tmp2 );
389
390         if ( iter->fDataType != AliHLTTPCDefinitions::fgkUnpackedRawDataType ) continue;
391
392       }
393         
394       slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
395       patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
396       row[0] = AliHLTTPCTransform::GetFirstRow( patch );
397       row[1] = AliHLTTPCTransform::GetLastRow( patch );
398         
399       if(fUnsorted){
400         if(fPadArray==NULL){
401           fClusterFinder->SetUnsorted(fUnsorted);
402           fPadArray = new AliHLTTPCPadArray(patch);
403           fPadArray->InitializeVector();
404         }
405         else if(fPadArray->GetPatch()!=patch||fPadArray->GetPatch()==-1){
406           if (GetEventCount()<3) {
407             HLTWarning("pad array not initialized for data of specification 0x%08x, block skipped", iter->fSpecification);
408           } else if ((GetEventCount()%5000)==0) { // assuming 0.5 to 1kHz this gives a message rate of 0.1 to 0.5 Hz
409             HLTWarning("reminder: pad array not initialized for data of specification 0x%08x", iter->fSpecification);
410           }
411           continue;
412         }
413       }
414
415       outPtr = (AliHLTTPCClusterData*)outBPtr;
416
417       maxPoints = (size-tSize-sizeof(AliHLTTPCClusterData))/sizeof(AliHLTTPCSpacePointData);
418
419       fClusterFinder->InitSlice( slice, patch, row[0], row[1], maxPoints );
420       fClusterFinder->SetOutputArray( (AliHLTTPCSpacePointData*)outPtr->fSpacePoints );
421         
422       if(fUnsorted){
423
424
425         fClusterFinder->SetPadArray(fPadArray);
426         /*        
427         double totalT=0;
428         struct timeval startT, endT;
429         gettimeofday( &startT, NULL );
430         */
431         fClusterFinder->ReadDataUnsorted(iter->fPtr, iter->fSize );
432         /*
433         gettimeofday( &endT, NULL );
434         unsigned long long dt;
435         dt = endT.tv_sec-startT.tv_sec;
436         dt *= 1000000ULL;
437         dt += endT.tv_usec-startT.tv_usec;
438         double dtd = ((double)dt);
439         totalT += dtd;
440         //        dtd = dtd / (double)eventIterations;
441         //        if ( iterations<=1 )
442         cout<<endl;
443         printf( "Time needed to read data: %f microsec. / %f millisec. / %f s\n", 
444                 dtd, dtd/1000.0, dtd/1000000.0 );
445           
446         cout<<endl;
447         */
448         fClusterFinder->FindClusters();
449       }
450       else{
451         fClusterFinder->Read(iter->fPtr, iter->fSize );
452         fClusterFinder->ProcessDigits();
453       }
454       realPoints = fClusterFinder->GetNumberOfClusters();
455         
456       outPtr->fSpacePointCnt = realPoints;
457       nSize = sizeof(AliHLTTPCSpacePointData)*realPoints;
458       mysize += nSize+sizeof(AliHLTTPCClusterData);
459
460       Logging( kHLTLogDebug, "HLT::TPCClusterFinder::DoEvent", "Spacepoints", 
461                "Number of spacepoints: %lu Slice/Patch/RowMin/RowMax: %d/%d/%d/%d.",
462                realPoints, slice, patch, row[0], row[1] );
463       AliHLTComponentBlockData bd;
464       FillBlockData( bd );
465       bd.fOffset = offset;
466       bd.fSize = mysize;
467       bd.fSpecification = iter->fSpecification;
468       //AliHLTSubEventDescriptor::FillBlockAttributes( bd.fAttributes );
469       outputBlocks.push_back( bd );
470         
471       tSize += mysize;
472       outBPtr += mysize;
473       outPtr = (AliHLTTPCClusterData*)outBPtr;
474         
475       if(fGetActivePads){
476         AliHLTTPCPadArray::AliHLTTPCActivePads* outPtrActive;
477         UInt_t activePadsSize, activePadsN = 0;
478         outPtrActive = (AliHLTTPCPadArray::AliHLTTPCActivePads*)outBPtr;
479         offset=tSize;
480         Int_t maxActivePads = (size-tSize)/sizeof(AliHLTTPCPadArray::AliHLTTPCActivePads);
481         activePadsSize= fClusterFinder->GetActivePads((AliHLTTPCPadArray::AliHLTTPCActivePads*)outPtrActive,maxActivePads)*sizeof(AliHLTTPCPadArray::AliHLTTPCActivePads);
482         
483         AliHLTComponentBlockData bdActive;
484         FillBlockData( bdActive );
485         bdActive.fOffset = offset;
486         bdActive.fSize = activePadsSize;
487         bdActive.fSpecification = iter->fSpecification;
488         bdActive.fDataType = AliHLTTPCDefinitions::fgkActivePadsDataType;
489         outputBlocks.push_back( bdActive );
490         
491         tSize+=activePadsSize;
492         outBPtr += activePadsSize;
493         outPtrActive = (AliHLTTPCPadArray::AliHLTTPCActivePads*)outBPtr;
494       }
495  
496
497       if ( tSize > size )
498         {
499           Logging( kHLTLogFatal, "HLT::TPCClusterFinder::DoEvent", "Too much data", 
500                    "Data written over allowed buffer. Amount written: %lu, allowed amount: %lu.",
501                    tSize, size );
502           return EMSGSIZE;
503         }
504     }
505     
506   size = tSize;
507
508   return 0;
509 }