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