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 #if defined(HAVE_ALIRAWDATA) && defined(HAVE_ALITPCRAWSTREAM_H)
 277       HLTDebug("using AliHLTTPCDigitReaderPacked");
 278       fReader = new AliHLTTPCDigitReaderPacked();
 279       if(oldRCUFormat==1){
 280         fReader->SetOldRCUFormat(kTRUE);
 281       }
 282       else if(oldRCUFormat!=0){
 283         HLTWarning("Wrong oldrcuformat specifier %d; oldrcuformat set to default(kFALSE)",oldRCUFormat);
 284       }
 285       if(fUnsorted==1){
 286         fReader->SetUnsorted(kTRUE);
 287       }
 288       fClusterFinder->SetReader(fReader);
 289 #else // ! defined(HAVE_ALIRAWDATA) && defined(HAVE_ALITPCRAWSTREAM_H)
 290       HLTFatal("DigitReaderPacked not available - check your build");
 291       return -ENODEV;
 292 #endif //  defined(HAVE_ALIRAWDATA) && defined(HAVE_ALITPCRAWSTREAM_H)
 293     } else {
 294 #if defined(HAVE_TPC_MAPPING)
 295       HLTDebug("using AliHLTTPCDigitReaderRaw mode %d", rawreadermode);
 296       fReader = new AliHLTTPCDigitReaderRaw(rawreadermode);
 297       fClusterFinder->SetReader(fReader);
 298 #else //! defined(HAVE_TPC_MAPPING)
 299       HLTFatal("DigitReaderRaw not available - check your build");
 300       return -ENODEV;
 301 #endif //defined(HAVE_TPC_MAPPING)
 302     }
 303   }
 304   else {
 305     HLTDebug("using AliHLTTPCDigitReaderUnpacked");
 306     fReader = new AliHLTTPCDigitReaderUnpacked();
 307     fClusterFinder->SetReader(fReader);
 308   }
 309
 310   // if pp-run use occupancy limit else set to 1. ==> use all
 311   if ( !fClusterDeconv )
 312     fClusterFinder->SetOccupancyLimit(occulimit);
 313   else
 314     fClusterFinder->SetOccupancyLimit(1.0);
 315
 316   // Variables to setup the Clusterfinder
 317   // TODO: this sounds strange and has to be verified; is the cluster finder not working when
 318   // fClusterDeconv = false ?
 319   fClusterDeconv = true;
 320   fXYClusterError = -1;
 321   fZClusterError = -1;
 322
 323
 324   fClusterFinder->SetDeconv( fClusterDeconv );
 325   fClusterFinder->SetXYError( fXYClusterError );
 326   fClusterFinder->SetZError( fZClusterError );
 327   if ( (fXYClusterError>0) && (fZClusterError>0) )
 328     fClusterFinder->SetCalcErr( false );
 329   fClusterFinder->SetSignalThreshold(sigthresh);
 330   fClusterFinder->SetNSigmaThreshold(sigmathresh);
 331
 332   return 0;
 333 }
 334
 335 int AliHLTTPCClusterFinderComponent::DoDeinit()
 336 {
 337   // see header file for class documentation
 338
 339   if ( fClusterFinder )
 340     delete fClusterFinder;
 341   fClusterFinder = NULL;
 342
 343   if ( fReader )
 344     delete fReader;
 345   fReader = NULL;
 346
 347   return 0;
 348 }
 349
 350 int AliHLTTPCClusterFinderComponent::DoEvent( const AliHLTComponentEventData& evtData,
 351                                               const AliHLTComponentBlockData* blocks,
 352                                               AliHLTComponentTriggerData& /*trigData*/, AliHLTUInt8_t* outputPtr,
 353                                               AliHLTUInt32_t& size,
 354                                               vector<AliHLTComponentBlockData>& outputBlocks )
 355 {
 356   // see header file for class documentation
 357
 358   //  == init iter (pointer to datablock)
 359   const AliHLTComponentBlockData* iter = NULL;
 360   unsigned long ndx;
 361
 362   //  == OUTdatatype pointer
 363   AliHLTTPCClusterData* outPtr;
 364
 365   AliHLTUInt8_t* outBPtr;
 366   UInt_t offset, mysize, nSize, tSize = 0;
 367
 368   outBPtr = outputPtr;
 369   outPtr = (AliHLTTPCClusterData*)outBPtr;
 370
 371   Int_t slice, patch, row[2];
 372   unsigned long maxPoints, realPoints = 0;
 373
 374   for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ )
 375     {
 376       iter = blocks+ndx;
 377       mysize = 0;
 378       offset = tSize;
 379
 380
 381       if (fPackedSwitch) {
 382         HLTDebug("Event 0x%08LX (%Lu) received datatype: %s - required datatype: %s",
 383                  evtData.fEventID, evtData.fEventID,
 384                  DataType2Text( iter->fDataType).c_str(),
 385                  DataType2Text(kAliHLTDataTypeDDLRaw | kAliHLTDataOriginTPC).c_str());
 386
 387         if (iter->fDataType == AliHLTTPCDefinitions::fgkDDLPackedRawDataType &&
 388             GetEventCount()<2) {
 389           HLTWarning("data type %s is depricated, use %s (kAliHLTDataTypeDDLRaw)!",
 390                      DataType2Text(AliHLTTPCDefinitions::fgkDDLPackedRawDataType).c_str(),
 391                      DataType2Text(kAliHLTDataTypeDDLRaw | kAliHLTDataOriginTPC).c_str());
 392           }
 393
 394         if ( iter->fDataType != (kAliHLTDataTypeDDLRaw | kAliHLTDataOriginTPC) &&
 395              iter->fDataType != AliHLTTPCDefinitions::fgkDDLPackedRawDataType ) continue;
 396
 397       }
 398       else {
 399         HLTDebug("Event 0x%08LX (%Lu) received datatype: %s - required datatype: %s",
 400                  evtData.fEventID, evtData.fEventID,
 401                  DataType2Text( iter->fDataType).c_str(),
 402                  DataType2Text(AliHLTTPCDefinitions::fgkUnpackedRawDataType).c_str());
 403
 404         if ( iter->fDataType != AliHLTTPCDefinitions::fgkUnpackedRawDataType ) continue;
 405
 406       }
 407
 408       slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
 409       patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
 410       row[0] = AliHLTTPCTransform::GetFirstRow( patch );
 411       row[1] = AliHLTTPCTransform::GetLastRow( patch );
 412
 413       if(fUnsorted){
 414         if(fPadArray==NULL){
 415           fClusterFinder->SetUnsorted(fUnsorted);
 416           fPadArray = new AliHLTTPCPadArray(patch);
 417           fPadArray->InitializeVector();
 418         }
 419         else if(fPadArray->GetPatch()!=patch||fPadArray->GetPatch()==-1){
 420           if (GetEventCount()<3) {
 421             HLTWarning("pad array not initialized for data of specification 0x%08x, block skipped", iter->fSpecification);
 422           } else if ((GetEventCount()%5000)==0) { // assuming 0.5 to 1kHz this gives a message rate of 0.1 to 0.5 Hz
 423             HLTWarning("reminder: pad array not initialized for data of specification 0x%08x", iter->fSpecification);
 424           }
 425           continue;
 426         }
 427       }
 428
 429       outPtr = (AliHLTTPCClusterData*)outBPtr;
 430
 431       maxPoints = (size-tSize-sizeof(AliHLTTPCClusterData))/sizeof(AliHLTTPCSpacePointData);
 432
 433       fClusterFinder->InitSlice( slice, patch, row[0], row[1], maxPoints );
 434       fClusterFinder->SetOutputArray( (AliHLTTPCSpacePointData*)outPtr->fSpacePoints );
 435
 436       if(fUnsorted){
 437
 438
 439         fClusterFinder->SetPadArray(fPadArray);
 440         /*
 441         double totalT=0;
 442         struct timeval startT, endT;
 443         gettimeofday( &startT, NULL );
 444         */
 445         fClusterFinder->ReadDataUnsorted(iter->fPtr, iter->fSize );
 446         /*
 447         gettimeofday( &endT, NULL );
 448         unsigned long long dt;
 449         dt = endT.tv_sec-startT.tv_sec;
 450         dt *= 1000000ULL;
 451         dt += endT.tv_usec-startT.tv_usec;
 452         double dtd = ((double)dt);
 453         totalT += dtd;
 454         //        dtd = dtd / (double)eventIterations;
 455         //        if ( iterations<=1 )
 456         cout<<endl;
 457         printf( "Time needed to read data: %f microsec. / %f millisec. / %f s\n",
 458                 dtd, dtd/1000.0, dtd/1000000.0 );
 459
 460         cout<<endl;
 461         */
 462         fClusterFinder->FindClusters();
 463       }
 464       else{
 465         fClusterFinder->Read(iter->fPtr, iter->fSize );
 466         fClusterFinder->ProcessDigits();
 467       }
 468       realPoints = fClusterFinder->GetNumberOfClusters();
 469
 470       outPtr->fSpacePointCnt = realPoints;
 471       nSize = sizeof(AliHLTTPCSpacePointData)*realPoints;
 472       mysize += nSize+sizeof(AliHLTTPCClusterData);
 473
 474       Logging( kHLTLogDebug, "HLT::TPCClusterFinder::DoEvent", "Spacepoints",
 475                "Number of spacepoints: %lu Slice/Patch/RowMin/RowMax: %d/%d/%d/%d.",
 476                realPoints, slice, patch, row[0], row[1] );
 477       AliHLTComponentBlockData bd;
 478       FillBlockData( bd );
 479       bd.fOffset = offset;
 480       bd.fSize = mysize;
 481       bd.fSpecification = iter->fSpecification;
 482       bd.fDataType = AliHLTTPCDefinitions::fgkClustersDataType;
 483       //AliHLTSubEventDescriptor::FillBlockAttributes( bd.fAttributes );
 484       outputBlocks.push_back( bd );
 485
 486       tSize += mysize;
 487       outBPtr += mysize;
 488       outPtr = (AliHLTTPCClusterData*)outBPtr;
 489
 490       if(fGetActivePads){
 491         AliHLTTPCPadArray::AliHLTTPCActivePads* outPtrActive;
 492         UInt_t activePadsSize, activePadsN = 0;
 493         outPtrActive = (AliHLTTPCPadArray::AliHLTTPCActivePads*)outBPtr;
 494         offset=tSize;
 495         Int_t maxActivePads = (size-tSize)/sizeof(AliHLTTPCPadArray::AliHLTTPCActivePads);
 496         activePadsSize= fClusterFinder->GetActivePads((AliHLTTPCPadArray::AliHLTTPCActivePads*)outPtrActive,maxActivePads)*sizeof(AliHLTTPCPadArray::AliHLTTPCActivePads);
 497
 498         AliHLTComponentBlockData bdActive;
 499         FillBlockData( bdActive );
 500         bdActive.fOffset = offset;
 501         bdActive.fSize = activePadsSize;
 502         bdActive.fSpecification = iter->fSpecification;
 503         bdActive.fDataType = AliHLTTPCDefinitions::fgkActivePadsDataType;
 504         outputBlocks.push_back( bdActive );
 505
 506         tSize+=activePadsSize;
 507         outBPtr += activePadsSize;
 508         outPtrActive = (AliHLTTPCPadArray::AliHLTTPCActivePads*)outBPtr;
 509       }
 510
 511
 512       if ( tSize > size )
 513         {
 514           Logging( kHLTLogFatal, "HLT::TPCClusterFinder::DoEvent", "Too much data",
 515                    "Data written over allowed buffer. Amount written: %lu, allowed amount: %lu.",
 516                    tSize, size );
 517           return EMSGSIZE;
 518         }
 519     }
 520
 521   size = tSize;
 522
 523   return 0;
 524 }