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