HLT/TPCLib/AliHLTTPCdEdxComponent.cxx

   1 // **************************************************************************
   2 // This file is property of and copyright by the ALICE HLT Project          *
   3 // ALICE Experiment at CERN, All rights reserved.                           *
   4 //                                                                          *
   5 // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
   6 //                  for The ALICE HLT Project.                              *
   7 //                                                                          *
   8 // Permission to use, copy, modify and distribute this software and its     *
   9 // documentation strictly for non-commercial purposes is hereby granted     *
  10 // without fee, provided that the above copyright notice appears in all     *
  11 // copies and that both the copyright notice and this permission notice     *
  12 // appear in the supporting documentation. The authors make no claims       *
  13 // about the suitability of this software for any purpose. It is            *
  14 // provided "as is" without express or implied warranty.                    *
  15 //                                                                          *
  16 //***************************************************************************
  17
  18 ///  @file   AliHLTTPCdEdxComponent.cxx
  19 ///  @author Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de>
  20 ///  @date   June 2009
  21 ///  @brief  An ITS tracker processing component for the HLT
  22
  23
  24 /////////////////////////////////////////////////////
  25 //                                                 //
  26 // dEdx calculation component for the HLT TPC       //
  27 //                                                 //
  28 /////////////////////////////////////////////////////
  29
  30 #include "AliHLTTPCdEdxComponent.h"
  31 #include "TStopwatch.h"
  32 #include "TMath.h"
  33 #include "AliCDBEntry.h"
  34 #include "AliCDBManager.h"
  35 #include "TObjString.h"
  36 #include "TObjArray.h"
  37 #include "AliHLTDataTypes.h"
  38 #include "AliHLTExternalTrackParam.h"
  39 #include "AliHLTGlobalBarrelTrack.h"
  40 #include "AliHLTTPCTransform.h"
  41 #include "AliHLTTPCSpacePointData.h"
  42 #include "AliHLTTPCClusterDataFormat.h"
  43 #include "AliHLTTPCTrackletDataFormat.h"
  44 #include "AliHLTTPCDefinitions.h"
  45 #include "AliTPCseed.h"
  46 #include "AliTPCclusterMI.h"
  47 #include "TGeoGlobalMagField.h"
  48
  49
  50 /** ROOT macro for the implementation of ROOT specific class methods */
  51 ClassImp( AliHLTTPCdEdxComponent )
  52 AliHLTTPCdEdxComponent::AliHLTTPCdEdxComponent()
  53     :
  54     fSolenoidBz( 0 ),
  55     fStatTime( 0 ),
  56     fStatNEvents( 0 )
  57 {
  58   // see header file for class documentation
  59   // or
  60   // refer to README to build package
  61   // or
  62   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  63   for( int i=0; i<fkNPatches; i++ ){
  64     fPatchClusters[i] = 0;
  65     fNPatchClusters[i] = 0;
  66   }
  67 }
  68
  69 AliHLTTPCdEdxComponent::AliHLTTPCdEdxComponent( const AliHLTTPCdEdxComponent& )
  70     :
  71     AliHLTProcessor(),
  72     fSolenoidBz( 0 ),
  73     fStatTime( 0 ),
  74     fStatNEvents( 0 )
  75 {
  76   // dummy
  77   for( int i=0; i<fkNPatches; i++ ){
  78     fPatchClusters[i] = 0;
  79     fNPatchClusters[i] = 0;
  80   }
  81   HLTFatal( "copy constructor untested" );
  82 }
  83
  84 AliHLTTPCdEdxComponent& AliHLTTPCdEdxComponent::operator=( const AliHLTTPCdEdxComponent& )
  85 {
  86   // see header file for class documentation
  87   HLTFatal( "assignment operator untested" );
  88   for( int i=0; i<fkNPatches; i++ ){
  89     fPatchClusters[i] = 0;
  90     fNPatchClusters[i] = 0;
  91   }
  92   return *this;
  93 }
  94
  95 AliHLTTPCdEdxComponent::~AliHLTTPCdEdxComponent()
  96 {
  97   // see header file for class documentation
  98   for( int i=0; i<fkNPatches; i++ ){
  99     delete[] fPatchClusters[i];
 100   }
 101 }
 102
 103 //
 104 // Public functions to implement AliHLTComponent's interface.
 105 // These functions are required for the registration process
 106 //
 107
 108 const char* AliHLTTPCdEdxComponent::GetComponentID()
 109 {
 110   // see header file for class documentation
 111   return "TPCdEdx";
 112 }
 113
 114 void AliHLTTPCdEdxComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list )
 115 {
 116   // see header file for class documentation
 117   list.clear();
 118   list.push_back( kAliHLTDataTypeTrack|kAliHLTDataOriginTPC );
 119   list.push_back( AliHLTTPCDefinitions::fgkClustersDataType );
 120 }
 121
 122 AliHLTComponentDataType AliHLTTPCdEdxComponent::GetOutputDataType()
 123 {
 124   // see header file for class documentation
 125   return kAliHLTDataTypedEdx|kAliHLTDataOriginTPC;
 126 }
 127
 128 void AliHLTTPCdEdxComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
 129 {
 130   // define guess for the output data size
 131   constBase = 200;       // minimum size
 132   inputMultiplier = 0.1; // size relative to input
 133 }
 134
 135 AliHLTComponent* AliHLTTPCdEdxComponent::Spawn()
 136 {
 137   // see header file for class documentation
 138   return new AliHLTTPCdEdxComponent;
 139 }
 140
 141 void AliHLTTPCdEdxComponent::SetDefaultConfiguration()
 142 {
 143   // Set default configuration for the CA tracker component
 144   // Some parameters can be later overwritten from the OCDB
 145
 146   fSolenoidBz = -5.00668;
 147   fStatTime = 0;
 148   fStatNEvents = 0;
 149 }
 150
 151 int AliHLTTPCdEdxComponent::ReadConfigurationString(  const char* arguments )
 152 {
 153   // Set configuration parameters for the CA tracker component from the string
 154
 155   int iResult = 0;
 156   if ( !arguments ) return iResult;
 157
 158   TString allArgs = arguments;
 159   TString argument;
 160   int bMissingParam = 0;
 161
 162   TObjArray* pTokens = allArgs.Tokenize( " " );
 163
 164   int nArgs =  pTokens ? pTokens->GetEntries() : 0;
 165
 166   for ( int i = 0; i < nArgs; i++ ) {
 167     argument = ( ( TObjString* )pTokens->At( i ) )->GetString();
 168     if ( argument.IsNull() ) continue;
 169
 170     if ( argument.CompareTo( "-solenoidBz" ) == 0 ) {
 171       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
 172       HLTWarning("argument -solenoidBz is deprecated, magnetic field set up globally (%f)", GetBz());
 173       continue;
 174     }
 175
 176     HLTError( "Unknown option \"%s\"", argument.Data() );
 177     iResult = -EINVAL;
 178   }
 179   delete pTokens;
 180
 181   if ( bMissingParam ) {
 182     HLTError( "Specifier missed for parameter \"%s\"", argument.Data() );
 183     iResult = -EINVAL;
 184   }
 185
 186   return iResult;
 187 }
 188
 189
 190 int AliHLTTPCdEdxComponent::ReadCDBEntry( const char* cdbEntry, const char* chainId )
 191 {
 192   // see header file for class documentation
 193
 194   const char* defaultNotify = "";
 195
 196   if ( !cdbEntry ) {
 197     return 0;// need to add the HLT/ConfigTPC/TPCdEdx directory to cdb SG!!!
 198     //cdbEntry = "HLT/ConfigTPC/TPCdEdx";
 199     //defaultNotify = " (default)";
 200     //chainId = 0;
 201   }
 202
 203   HLTInfo( "configure from entry \"%s\"%s, chain id %s", cdbEntry, defaultNotify, ( chainId != NULL && chainId[0] != 0 ) ? chainId : "<none>" );
 204   AliCDBEntry *pEntry = AliCDBManager::Instance()->Get( cdbEntry );//,GetRunNo());
 205
 206   if ( !pEntry ) {
 207     HLTError( "cannot fetch object \"%s\" from CDB", cdbEntry );
 208     return -EINVAL;
 209   }
 210
 211   TObjString* pString = dynamic_cast<TObjString*>( pEntry->GetObject() );
 212
 213   if ( !pString ) {
 214     HLTError( "configuration object \"%s\" has wrong type, required TObjString", cdbEntry );
 215     return -EINVAL;
 216   }
 217
 218   HLTInfo( "received configuration object string: \"%s\"", pString->GetString().Data() );
 219
 220   return  ReadConfigurationString( pString->GetString().Data() );
 221 }
 222
 223
 224 int AliHLTTPCdEdxComponent::Configure( const char* cdbEntry, const char* chainId, const char *commandLine )
 225 {
 226   // Configure the component
 227   // There are few levels of configuration,
 228   // parameters which are set on one step can be overwritten on the next step
 229
 230   //* read hard-coded values
 231
 232   SetDefaultConfiguration();
 233
 234   //* read the default CDB entry
 235
 236   int iResult1 = ReadCDBEntry( NULL, chainId );
 237
 238   //* read magnetic field
 239
 240   fSolenoidBz=GetBz();
 241
 242   //* read the actual CDB entry if required
 243
 244   int iResult2 = ( cdbEntry ) ? ReadCDBEntry( cdbEntry, chainId ) : 0;
 245
 246   //* read extra parameters from input (if they are)
 247
 248   int iResult3 = 0;
 249
 250   if ( commandLine && commandLine[0] != '\0' ) {
 251     HLTInfo( "received configuration string from HLT framework: \"%s\"", commandLine );
 252     iResult3 = ReadConfigurationString( commandLine );
 253   }
 254
 255   // Initialise the tracker here
 256
 257   return iResult1 ? iResult1 : ( iResult2 ? iResult2 : iResult3 );
 258 }
 259
 260
 261
 262 int AliHLTTPCdEdxComponent::DoInit( int argc, const char** argv )
 263 {
 264   // Configure the component
 265
 266   TString arguments = "";
 267   for ( int i = 0; i < argc; i++ ) {
 268     if ( !arguments.IsNull() ) arguments += " ";
 269     arguments += argv[i];
 270   }
 271
 272   int ret = Configure( NULL, NULL, arguments.Data() );
 273
 274   // Check field
 275   if (!TGeoGlobalMagField::Instance()) {
 276     HLTError("magnetic field not initialized, please set up TGeoGlobalMagField and AliMagF");
 277     return -ENODEV;
 278   }
 279
 280   //AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform();
 281   //if( transform ){
 282   //AliTPCRecoParam *reco = transform->GetCurrentRecoParam();
 283   //if( ! reco ){
 284   //reco = new AliTPCRecoParam;
 285   //reco->SetUseTotalCharge(0);
 286   //transform->SetCurrentRecoParam( reco );
 287   //}
 288
 289   return ret;
 290 }
 291
 292
 293 int AliHLTTPCdEdxComponent::DoDeinit()
 294 {
 295   // see header file for class documentation
 296   return 0;
 297 }
 298
 299
 300
 301 int AliHLTTPCdEdxComponent::Reconfigure( const char* cdbEntry, const char* chainId )
 302 {
 303   // Reconfigure the component from OCDB .
 304
 305   return Configure( cdbEntry, chainId, NULL );
 306 }
 307
 308
 309
 310 int AliHLTTPCdEdxComponent::DoEvent
 311 (
 312   const AliHLTComponentEventData& evtData,
 313   const AliHLTComponentBlockData* blocks,
 314   AliHLTComponentTriggerData& /*trigData*/,
 315   AliHLTUInt8_t* outputPtr,
 316   AliHLTUInt32_t& size,
 317   vector<AliHLTComponentBlockData>& outputBlocks )
 318 {
 319   //* process the event
 320
 321   AliHLTUInt32_t maxBufferSize = size;
 322   size = 0; // output size
 323
 324   if (!IsDataEvent()) return 0;
 325
 326   if ( evtData.fBlockCnt <= 0 ) {
 327     HLTWarning( "no blocks in event" );
 328     return 0;
 329   }
 330   if ( !outputPtr ) {
 331     return -ENOSPC;
 332   }
 333
 334   int iResult=0;
 335
 336   TStopwatch timer;
 337
 338   // Initialise the data pointers
 339
 340   for( int i=0; i<fkNPatches; i++ ){
 341     delete[] fPatchClusters[i];
 342     fPatchClusters[i] = 0;
 343     fNPatchClusters[i] = 0;
 344   }
 345
 346   int nBlocks = (int)evtData.fBlockCnt;
 347
 348   int nInputClusters = 0;
 349   int nInputTracks = 0;
 350
 351   // first read all the clusters
 352
 353   for (int ndx=0; ndx<nBlocks && iResult>=0; ndx++) {
 354     const AliHLTComponentBlockData* iter = blocks+ndx;
 355     if ( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType ) continue;
 356     Int_t slice=AliHLTTPCDefinitions::GetMinSliceNr(iter->fSpecification);
 357     Int_t patch=AliHLTTPCDefinitions::GetMinPatchNr(iter->fSpecification);
 358     Int_t slicepatch=slice*6+patch;
 359     if( slicepatch >= fkNPatches ){
 360       HLTWarning("Wrong header of TPC cluster data, slice %d, patch %d",
 361                  slice, patch );
 362       continue;
 363     }
 364     AliHLTTPCClusterData* inPtrSP = ( AliHLTTPCClusterData* )( iter->fPtr );
 365     nInputClusters += inPtrSP->fSpacePointCnt;
 366
 367     delete[] fPatchClusters[slicepatch];
 368     fPatchClusters[slicepatch] = new AliTPCclusterMI[inPtrSP->fSpacePointCnt];
 369     fNPatchClusters[slicepatch] = inPtrSP->fSpacePointCnt;
 370
 371     // create  off-line clusters out of the HLT clusters
 372     // todo: check which cluster information is really needed for the dEdx
 373
 374     for ( unsigned int i = 0; i < inPtrSP->fSpacePointCnt; i++ ) {
 375       AliHLTTPCSpacePointData *chlt = &( inPtrSP->fSpacePoints[i] );
 376       AliTPCclusterMI *c = fPatchClusters[slicepatch]+i;
 377       c->SetX(chlt->fX);
 378       c->SetY(chlt->fY);
 379       c->SetZ(chlt->fZ);
 380       c->SetSigmaY2(chlt->fSigmaY2);
 381       c->SetSigmaYZ( 0 );
 382       c->SetSigmaZ2(chlt->fSigmaZ2);
 383       c->SetQ( chlt->fCharge );
 384       c->SetMax( chlt->fQMax );
 385       Int_t sector, row;
 386       Float_t padtime[3]={0,chlt->fY,chlt->fZ};
 387       AliHLTTPCTransform::Slice2Sector(slice,chlt->fPadRow, sector, row);
 388       AliHLTTPCTransform::Local2Raw( padtime, sector, row);
 389       c->SetDetector( sector );
 390       c->SetRow( row );
 391       c->SetPad( (Int_t) padtime[1] );
 392       c->SetTimeBin( (Int_t) padtime[2] );
 393     }
 394   }
 395
 396
 397   // loop over the input tracks: calculate dEdx and write output
 398
 399   unsigned int outSize = 0;
 400   AliHLTFloat32_t *outPtr = ( AliHLTFloat32_t * )( outputPtr );
 401
 402   for (const AliHLTComponentBlockData* pBlock=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
 403        pBlock!=NULL; pBlock=GetNextInputBlock()) {
 404
 405     AliHLTComponentBlockData outBlock;
 406     FillBlockData( outBlock );
 407     outBlock.fOffset = outSize;
 408     outBlock.fSize = 0;
 409     outBlock.fDataType = kAliHLTDataTypedEdx|kAliHLTDataOriginTPC;
 410     outBlock.fSpecification = pBlock->fSpecification;
 411
 412     AliHLTTracksData* dataPtr = ( AliHLTTracksData* ) pBlock->fPtr;
 413     int nTracks = dataPtr->fCount;
 414     AliHLTExternalTrackParam* currTrack = dataPtr->fTracklets;
 415     nInputTracks+=nTracks;
 416
 417     for( int itr=0;
 418          itr<nTracks && ( (AliHLTUInt8_t *)currTrack < ((AliHLTUInt8_t *) pBlock->fPtr)+pBlock->fSize);
 419          itr++ ){
 420
 421       // create an off-line track
 422       AliHLTGlobalBarrelTrack gb(*currTrack);
 423       AliTPCseed tTPC;
 424       tTPC.Set( gb.GetX(), gb.GetAlpha(),
 425                 gb.GetParameter(), gb.GetCovariance() );
 426
 427       // set the clusters
 428
 429       for( UInt_t ic=0; ic<currTrack->fNPoints; ic++){
 430         UInt_t id = currTrack->fPointIDs[ic];
 431         int iSlice = AliHLTTPCSpacePointData::GetSlice(id);
 432         int iPatch = AliHLTTPCSpacePointData::GetPatch(id);
 433         int iCluster = AliHLTTPCSpacePointData::GetNumber(id);
 434         if( iSlice<0 || iSlice>36 || iPatch<0 || iPatch>5 ){
 435           HLTError("Corrupted TPC cluster Id: slice %d, patch %d, cluster %d",
 436                    iSlice, iPatch,iCluster );
 437           continue;
 438         }
 439         AliTPCclusterMI *patchClusters = fPatchClusters[iSlice*6 + iPatch];
 440         if( !patchClusters ){
 441           HLTError("Clusters are missed for slice %d, patch %d",
 442                    iSlice, iPatch );
 443           continue;
 444         }
 445         if( iCluster >= fNPatchClusters[iSlice*6 + iPatch] ){
 446           HLTError("TPC slice %d, patch %d: ClusterID==%d >= N Cluaters==%d ",
 447                    iSlice, iPatch,iCluster, fNPatchClusters[iSlice*6 + iPatch] );
 448           continue;
 449         }
 450         AliTPCclusterMI *c = &(patchClusters[iCluster]);
 451
 452         int sec = c->GetDetector();
 453         int row = c->GetRow();
 454         if ( sec >= 36 ) {
 455           row = row + AliHLTTPCTransform::GetNRowLow();
 456         }
 457         tTPC.SetClusterPointer( row, c );
 458
 459         AliTPCTrackerPoint &point = *( tTPC.GetTrackPoint( row ) );
 460         tTPC.Propagate( TMath::DegToRad()*(sec%18*20.+10.), c->GetX(), fSolenoidBz );
 461         Double_t angle2 = tTPC.GetSnp()*tTPC.GetSnp();
 462         angle2 = (angle2<1) ?TMath::Sqrt(angle2/(1-angle2)) :10.;
 463         point.SetAngleY( angle2 );
 464         point.SetAngleZ( tTPC.GetTgl() );
 465       }
 466
 467       // Cook dEdx
 468
 469       if( outSize+3*sizeof( AliHLTFloat32_t ) > maxBufferSize ){
 470         HLTWarning( "Output buffer size %d exceed", maxBufferSize );
 471         iResult = -ENOSPC;
 472         break;
 473       }
 474
 475       tTPC.CookdEdx( 0.02, 0.6 );
 476       outPtr[0] = tTPC.GetdEdx();
 477       outPtr[1] = tTPC.GetSDEDX(0);
 478       outPtr[2] = tTPC.GetNCDEDX(0);
 479       outPtr+=3;
 480       outSize+=3*sizeof( AliHLTFloat32_t );
 481       outBlock.fSize+=3*sizeof( AliHLTFloat32_t );
 482
 483       unsigned int step = sizeof( AliHLTExternalTrackParam ) + currTrack->fNPoints * sizeof( unsigned int );
 484       currTrack = ( AliHLTExternalTrackParam* )( (( Byte_t * )currTrack) + step );
 485     }
 486     if( iResult<0 ) break;
 487     outputBlocks.push_back( outBlock );
 488     size = outSize;
 489   }
 490
 491   timer.Stop();
 492   fStatTime += timer.RealTime();
 493   fStatNEvents++;
 494
 495   // Set log level to "Warning" for on-line system monitoring
 496   int hz = ( int ) ( fStatTime > 1.e-10 ? fStatNEvents / fStatTime : 100000 );
 497
 498   HLTInfo( "TPCdEdx : %d clusters and %d tracks processed; average time %d Hz",
 499            nInputClusters, nInputTracks, hz );
 500
 501   return iResult;
 502 }