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