HLT/TPCLib/tracking-ca/AliHLTTPCCATrackerComponent.cxx

   1 // @(#) $Id$
   2 /**************************************************************************
   3  * This file is property of and copyright by the ALICE HLT Project        *
   4  * ALICE Experiment at CERN, All rights reserved.                         *
   5  *                                                                        *
   6  * Primary Authors: Jochen Thaeder <thaeder@kip.uni-heidelberg.de>        *
   7  *                  Ivan Kisel <kisel@kip.uni-heidelberg.de>              *
   8  *                  for The ALICE HLT Project.                            *
   9  *                                                                        *
  10  * Permission to use, copy, modify and distribute this software and its   *
  11  * documentation strictly for non-commercial purposes is hereby granted   *
  12  * without fee, provided that the above copyright notice appears in all   *
  13  * copies and that both the copyright notice and this permission notice   *
  14  * appear in the supporting documentation. The authors make no claims     *
  15  * about the suitability of this software for any purpose. It is          *
  16  * provided "as is" without express or implied warranty.                  *
  17  **************************************************************************/
  18
  19 ///////////////////////////////////////////////////////////////////////////////
  20 //                                                                           //
  21 // a TPC tracker processing component for the HLT based on CA by Ivan Kisel  //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 #if __GNUC__>= 3
  26 using namespace std;
  27 #endif
  28
  29 #include "AliHLTTPCCATrackerComponent.h"
  30 #include "AliHLTTPCTransform.h"
  31 #include "AliHLTTPCCATracker.h"
  32 #include "AliHLTTPCCAHit.h"
  33 #include "AliHLTTPCCAOutTrack.h"
  34
  35 #include "AliHLTTPCSpacePointData.h"
  36 #include "AliHLTTPCClusterDataFormat.h"
  37 #include "AliHLTTPCTransform.h"
  38 #include "AliHLTTPCTrackSegmentData.h"
  39 #include "AliHLTTPCTrackArray.h"
  40 #include "AliHLTTPCTrackletDataFormat.h"
  41 #include "AliHLTTPCDefinitions.h"
  42 #include "TStopwatch.h"
  43 #include "TMath.h"
  44
  45
  46 // this is a global object used for automatic component registration, do not use this
  47 AliHLTTPCCATrackerComponent gAliHLTTPCCATrackerComponent;
  48
  49 ClassImp(AliHLTTPCCATrackerComponent)
  50
  51 AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent()
  52   :
  53   fTracker(NULL),
  54   fBField(0)
  55 {
  56   // see header file for class documentation
  57   // or
  58   // refer to README to build package
  59   // or
  60   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  61 }
  62
  63 AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent(const AliHLTTPCCATrackerComponent&)
  64   :
  65   AliHLTProcessor(),
  66   fTracker(NULL),
  67   fBField(0)
  68 {
  69   // see header file for class documentation
  70   HLTFatal("copy constructor untested");
  71 }
  72
  73 AliHLTTPCCATrackerComponent& AliHLTTPCCATrackerComponent::operator=(const AliHLTTPCCATrackerComponent&)
  74 {
  75   // see header file for class documentation
  76   HLTFatal("assignment operator untested");
  77   return *this;
  78 }
  79
  80 AliHLTTPCCATrackerComponent::~AliHLTTPCCATrackerComponent()
  81 {
  82   // see header file for class documentation
  83   delete fTracker;
  84 }
  85
  86 //
  87 // Public functions to implement AliHLTComponent's interface.
  88 // These functions are required for the registration process
  89 //
  90
  91 const char* AliHLTTPCCATrackerComponent::GetComponentID()
  92 {
  93   // see header file for class documentation
  94   return "TPCCATracker";
  95 }
  96
  97 void AliHLTTPCCATrackerComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list)
  98 {
  99   // see header file for class documentation
 100   list.clear();
 101   list.push_back( AliHLTTPCDefinitions::fgkClustersDataType );
 102 }
 103
 104 AliHLTComponentDataType AliHLTTPCCATrackerComponent::GetOutputDataType()
 105 {
 106   // see header file for class documentation
 107   return AliHLTTPCDefinitions::fgkTrackSegmentsDataType;
 108 }
 109
 110 void AliHLTTPCCATrackerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
 111 {
 112   // define guess for the output data size
 113   constBase = 200;       // minimum size
 114   inputMultiplier = 0.5; // size relative to input
 115 }
 116
 117 AliHLTComponent* AliHLTTPCCATrackerComponent::Spawn()
 118 {
 119   // see header file for class documentation
 120   return new AliHLTTPCCATrackerComponent;
 121 }
 122
 123 int AliHLTTPCCATrackerComponent::DoInit( int argc, const char** argv )
 124 {
 125   // Initialize the CA tracker component
 126   //
 127   // arguments could be:
 128   // bfield - the magnetic field value
 129   //
 130
 131   if ( fTracker ) return EINPROGRESS;
 132
 133   fTracker = new AliHLTTPCCATracker();
 134
 135   // read command line
 136
 137   int i = 0;
 138   char* cpErr;
 139   while ( i < argc ){
 140     if ( !strcmp( argv[i], "bfield" ) ){
 141       if ( i+1 >= argc )
 142         {
 143           Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing B-field", "Missing B-field specifier." );
 144           return ENOTSUP;
 145         }
 146       fBField = strtod( argv[i+1], &cpErr );
 147       if ( *cpErr )
 148         {
 149           Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert B-field specifier '%s'.", argv[i+1] );
 150           return EINVAL;
 151         }
 152
 153       Logging( kHLTLogDebug, "HLT::TPCCATracker::DoInit", "Reading command line",
 154                "Magnetic field value is set to %f kG", fBField );
 155
 156       i += 2;
 157       continue;
 158     }
 159
 160     Logging(kHLTLogError, "HLT::TPCCATracker::DoInit", "Unknown Option", "Unknown option '%s'", argv[i] );
 161     return EINVAL;
 162   }
 163
 164   return 0;
 165 }
 166
 167 int AliHLTTPCCATrackerComponent::DoDeinit()
 168 {
 169   // see header file for class documentation
 170   if ( fTracker ) delete fTracker;
 171   fTracker = NULL;
 172   return 0;
 173 }
 174
 175 int AliHLTTPCCATrackerComponent::DoEvent
 176 (
 177  const AliHLTComponentEventData& evtData,
 178  const AliHLTComponentBlockData* blocks,
 179  AliHLTComponentTriggerData& /*trigData*/,
 180  AliHLTUInt8_t* outputPtr,
 181  AliHLTUInt32_t& size,
 182  vector<AliHLTComponentBlockData>& outputBlocks )
 183 {
 184
 185   AliHLTUInt32_t MaxBufferSize = size;
 186   size = 0; // output size
 187
 188   TStopwatch timer;
 189
 190   // Event reconstruction in one TPC slice with CA Tracker
 191
 192   Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "DoEvent", "DoEvent()" );
 193   if ( evtData.fBlockCnt<=0 )
 194     {
 195       Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "no blocks in event" );
 196       return 0;
 197     }
 198
 199   const AliHLTComponentBlockData* iter = NULL;
 200   unsigned long ndx;
 201   AliHLTTPCClusterData* inPtrSP;
 202
 203   // Determine the slice number
 204
 205   Int_t slice=-1;
 206   {
 207     std::vector<Int_t> slices;
 208     std::vector<Int_t>::iterator slIter;
 209     std::vector<unsigned> sliceCnts;
 210     std::vector<unsigned>::iterator slCntIter;
 211
 212     for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ ){
 213       iter = blocks+ndx;
 214       if ( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType ) continue;
 215
 216       slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
 217
 218       Bool_t found = 0;
 219       slCntIter = sliceCnts.begin();
 220       for( slIter = slices.begin(); slIter!=slices.end(); slIter++, slCntIter++ ){
 221         if ( *slIter == slice ){
 222           found = kTRUE;
 223           break;
 224         }
 225       }
 226       if ( !found ){
 227         slices.push_back( slice );
 228         sliceCnts.push_back( 1 );
 229       } else *slCntIter++;
 230     }
 231
 232
 233     // Determine slice number to really use.
 234     if ( slices.size()>1 )
 235       {
 236         Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
 237                  "Multiple slice numbers found in event 0x%08lX (%lu). Determining maximum occuring slice number...",
 238                  evtData.fEventID, evtData.fEventID );
 239         unsigned maxCntSlice=0;
 240         slCntIter = sliceCnts.begin();
 241         for( slIter = slices.begin(); slIter != slices.end(); slIter++, slCntIter++ )
 242           {
 243             Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
 244                      "Slice %lu found %lu times.", *slIter, *slCntIter );
 245             if ( maxCntSlice<*slCntIter )
 246               {
 247                 maxCntSlice = *slCntIter;
 248                 slice = *slIter;
 249               }
 250           }
 251         Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
 252                  "Using slice %lu.", slice );
 253       }
 254     else if ( slices.size()>0 )
 255       {
 256         slice = *(slices.begin());
 257       }
 258   }
 259
 260   if( slice<0 ){
 261     Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "no slices found in event" );
 262     return 0;
 263   }
 264
 265
 266   // Initialize the tracker
 267
 268   Double_t Bz = fBField;
 269
 270   {
 271     if( !fTracker ) fTracker = new AliHLTTPCCATracker;
 272     Int_t iSec = slice;
 273     Double_t inRmin = 83.65;
 274     //    Double_t inRmax = 133.3;
 275     //    Double_t outRmin = 133.5;
 276     Double_t outRmax = 247.7;
 277     Double_t plusZmin = 0.0529937;
 278     Double_t plusZmax = 249.778;
 279     Double_t minusZmin = -249.645;
 280     Double_t minusZmax = -0.0799937;
 281     Double_t dalpha = 0.349066;
 282     Double_t alpha = 0.174533 + dalpha*iSec;
 283
 284     Bool_t zPlus = (iSec<18 );
 285     Double_t zMin =  zPlus ?plusZmin :minusZmin;
 286     Double_t zMax =  zPlus ?plusZmax :minusZmax;
 287     //TPCZmin = -249.645, ZMax = 249.778
 288     //    Double_t rMin =  inRmin;
 289     //    Double_t rMax =  outRmax;
 290     Int_t NRows = AliHLTTPCTransform::GetNRows();
 291
 292     Double_t padPitch = 0.4;
 293     Double_t sigmaZ = 0.228808;
 294
 295      Double_t rowX[NRows];
 296     for( Int_t irow=0; irow<NRows; irow++){
 297       rowX[irow] = AliHLTTPCTransform::Row2X( irow );
 298     }
 299
 300     AliHLTTPCCAParam param;
 301     param.Initialize( iSec, NRows, rowX, alpha, dalpha,
 302                       inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, Bz );
 303     param.YErrorCorrection() = 1;
 304     param.ZErrorCorrection() = 2;
 305
 306     fTracker->Initialize( param );
 307
 308   }
 309
 310
 311   // min and max patch numbers and row numbers
 312
 313   Int_t row[2] = {0,0};
 314   Int_t minPatch=INT_MAX, maxPatch = -1;
 315
 316   // total n Hits
 317
 318   Int_t nHitsTotal = 0;
 319
 320   // sort patches
 321
 322   std::vector<unsigned long> patchIndices;
 323
 324   for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ ){
 325     iter = blocks+ndx;
 326     if( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType ) continue;
 327     if( slice!=AliHLTTPCDefinitions::GetMinSliceNr( *iter ) ) continue;
 328     inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
 329     nHitsTotal+=inPtrSP->fSpacePointCnt;
 330     Int_t patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
 331     if ( minPatch>patch ){
 332       minPatch = patch;
 333       row[0] = AliHLTTPCTransform::GetFirstRow( patch );
 334     }
 335     if ( maxPatch<patch ){
 336       maxPatch = patch;
 337       row[1] = AliHLTTPCTransform::GetLastRow( patch );
 338     }
 339     std::vector<unsigned long>::iterator pIter = patchIndices.begin();
 340     while( pIter!=patchIndices.end() && AliHLTTPCDefinitions::GetMinPatchNr( blocks[*pIter] ) < patch ){
 341       pIter++;
 342     }
 343     patchIndices.insert( pIter, ndx );
 344   }
 345
 346
 347   // pass event to CA Tracker
 348
 349   fTracker->StartEvent();
 350
 351   AliHLTTPCCAHit vHits[nHitsTotal]; // CA hit array
 352   Double_t vHitStoreX[nHitsTotal];       // hit X coordinates
 353   Int_t vHitStoreID[nHitsTotal];            // hit ID's
 354   Int_t vHitRowID[nHitsTotal];            // hit ID's
 355
 356   Int_t nHits = 0;
 357
 358   Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits",
 359            "Total %d hits to read for slice %d", nHitsTotal, slice );
 360
 361   Int_t nClusters=0;
 362
 363   for( std::vector<unsigned long>::iterator pIter = patchIndices.begin(); pIter!=patchIndices.end(); pIter++ ){
 364     ndx = *pIter;
 365     iter = blocks+ndx;
 366
 367     Int_t patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
 368     inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
 369
 370     Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits",
 371              "Reading %d hits for slice %d - patch %d", inPtrSP->fSpacePointCnt, slice, patch );
 372
 373     // Read patch hits, row by row
 374
 375     Int_t oldRow = -1;
 376     Int_t nRowHits = 0;
 377     Int_t firstRowHit = 0;
 378     for (UInt_t i=0; i<inPtrSP->fSpacePointCnt; i++ ){
 379       AliHLTTPCSpacePointData* pSP = &(inPtrSP->fSpacePoints[i]);
 380
 381       if( pSP->fPadRow != oldRow ){
 382         if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
 383         oldRow = pSP->fPadRow;
 384         firstRowHit = nHits;
 385         nRowHits = 0;
 386       }
 387       AliHLTTPCCAHit &h = vHits[nHits];
 388       if( TMath::Abs(pSP->fX- fTracker->Rows()[pSP->fPadRow].X() )>1.e-4 ) cout<<"row "<<(Int_t)pSP->fPadRow<<" "<<fTracker->Rows()[pSP->fPadRow].X()-pSP->fX <<endl;
 389
 390       h.Y() = pSP->fY;
 391       h.Z() = pSP->fZ;
 392       h.ErrY() = TMath::Sqrt(TMath::Abs(pSP->fSigmaY2));
 393       h.ErrZ() = TMath::Sqrt(TMath::Abs(pSP->fSigmaZ2));
 394       if( h.ErrY()<.1 ) h.ErrY() = .1;
 395       if( h.ErrZ()<.1 ) h.ErrZ() = .1;
 396       if( h.ErrY()>1. ) h.ErrY() = 1.;
 397       if( h.ErrZ()>1. ) h.ErrZ() = 1.;
 398       h.ID() = nHits;
 399       vHitStoreX[nHits] = pSP->fX;
 400       vHitStoreID[nHits] = pSP->fID;
 401       vHitRowID[nHits] = pSP->fPadRow;
 402       nHits++;
 403       nRowHits++;
 404       nClusters++;
 405     }
 406     if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
 407   }
 408
 409   // reconstruct the event
 410
 411   TStopwatch timerReco;
 412
 413   fTracker->Reconstruct();
 414
 415   timerReco.Stop();
 416
 417   Int_t ret = 0;
 418
 419   Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reconstruct",
 420            "%d tracks found for slice %d",fTracker->NOutTracks(), slice);
 421
 422   // write reconstructed tracks
 423
 424   AliHLTTPCTrackletData* outPtr = (AliHLTTPCTrackletData*)(outputPtr);
 425
 426   AliHLTTPCTrackSegmentData* currOutTracklet = outPtr->fTracklets;
 427
 428   Int_t ntracks = fTracker->NOutTracks();
 429
 430   UInt_t mySize =   ((AliHLTUInt8_t *)currOutTracklet) -  ((AliHLTUInt8_t *)outputPtr);
 431
 432   outPtr->fTrackletCnt = 0;
 433
 434   for( int itr=0; itr<ntracks; itr++ ){
 435
 436     AliHLTTPCCAOutTrack &t = fTracker->OutTracks()[itr];
 437
 438     //Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Wrtite output","track %d with %d hits", itr, t.NHits());
 439
 440     // calculate output track size
 441
 442     UInt_t dSize = sizeof(AliHLTTPCTrackSegmentData) + t.NHits()*sizeof(UInt_t);
 443
 444     if( mySize + dSize > MaxBufferSize ){
 445       Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "Wrtite output","Output buffer size exceed (buffer size %d, current size %d), %d tracks are not stored", MaxBufferSize, mySize, ntracks-itr+1);
 446       ret = -ENOSPC;
 447       break;
 448     }
 449
 450     // convert CA track parameters to HLT Track Segment
 451
 452     Int_t iFirstHit = fTracker->OutTrackHits()[t.FirstHitRef()];
 453     Int_t iLastHit = fTracker->OutTrackHits()[t.FirstHitRef()+t.NHits()-1];
 454
 455     AliHLTTPCCAHit &firstHit = vHits[iFirstHit];
 456     AliHLTTPCCAHit &lastHit = vHits[iLastHit];
 457
 458     t.Param().TransportBz(Bz, vHitStoreX[iFirstHit], firstHit.Y(), firstHit.Z() );
 459
 460     currOutTracklet->fX = t.Param().Par()[0];
 461     currOutTracklet->fY = t.Param().Par()[1];
 462     currOutTracklet->fZ = t.Param().Par()[2];
 463     Double_t qp = t.Param().Par()[6];
 464     Double_t p = TMath::Abs(qp)>1.e-5 ?1./TMath::Abs(qp) :1.e5;
 465     Double_t ex = t.Param().Par()[3];
 466     Double_t ey = t.Param().Par()[4];
 467     Double_t ez = t.Param().Par()[5];
 468     Double_t et = TMath::Sqrt( ex*ex + ey*ey );
 469
 470     currOutTracklet->fCharge = (qp>0) ?+1 :(qp<0 ?-1 :0);
 471     currOutTracklet->fPt = p*et;
 472
 473     Double_t h3 =  TMath::Abs(ex) >1.e-5 ? p*ex/et :0;
 474     Double_t h4 =  TMath::Abs(ey) >1.e-5 ? p*ey/et :0;
 475     Double_t h5;
 476     Double_t h6 =  - currOutTracklet->fCharge * p * currOutTracklet->fPt;
 477
 478     currOutTracklet->fPterr = ( h3*h3*t.Param().Cov()[9] + h4*h4*t.Param().Cov()[14] + h6*h6*t.Param().Cov()[27]
 479                                 + 2.*(h3*h4*t.Param().Cov()[13]+h3*h6*t.Param().Cov()[24]+h4*h6*t.Param().Cov()[25] )
 480                                 );
 481     currOutTracklet->fPsi = TMath::ATan2(ey, ex);
 482
 483     h3 =  ex/(et*et);
 484     h4 = -ey/(et*et);
 485     currOutTracklet->fPsierr = h3*h3*t.Param().Cov()[9] + h4*h4*t.Param().Cov()[14] + 2.*h3*h4*t.Param().Cov()[13];
 486
 487     currOutTracklet->fTgl = TMath::Abs(et)>1.e-5  ? ez/et :1.e5;
 488
 489     if( TMath::Abs(et) >1.e-2 ){
 490       h3 = -ez*ex/et/et;
 491       h4 = -ez*ey/et/et;
 492       h5 = 1.;
 493       currOutTracklet->fTglerr =  ( h3*h3*t.Param().Cov()[9] + h4*h4*t.Param().Cov()[14] + h5*h5*t.Param().Cov()[20]
 494                                     + 2.*(h3*h4*t.Param().Cov()[13]+h3*h5*t.Param().Cov()[18]+h4*h5*t.Param().Cov()[19] )
 495                                     )/et/et;
 496     }else{
 497       currOutTracklet->fTglerr =  1.e-2;
 498     }
 499
 500     currOutTracklet->fCharge = -currOutTracklet->fCharge;
 501
 502     t.Param().TransportBz(Bz, vHitStoreX[iLastHit], lastHit.Y(), lastHit.Z() );
 503
 504     currOutTracklet->fLastX = t.Param().Par()[0];
 505     currOutTracklet->fLastY = t.Param().Par()[1];
 506     currOutTracklet->fLastZ = t.Param().Par()[2];
 507
 508 #ifdef INCLUDE_TPC_HOUGH
 509 #ifdef ROWHOUGHPARAMS
 510     currOutTracklet->fTrackID = 0;
 511     currOutTracklet->fRowRange1 = vHitRowID[iFirstHit];
 512     currOutTracklet->fRowRange2 = vHitRowID[iLastHit];
 513     currOutTracklet->fSector = slice;
 514     currOutTracklet->fPID = 211;
 515 #endif
 516 #endif // INCLUDE_TPC_HOUGH
 517
 518
 519     currOutTracklet->fNPoints = t.NHits();
 520
 521     for( Int_t i=0; i<t.NHits(); i++ ){
 522       currOutTracklet->fPointIDs[i] = vHitStoreID[fTracker->OutTrackHits()[t.FirstHitRef()+i]];
 523     }
 524
 525     currOutTracklet = (AliHLTTPCTrackSegmentData*)( (Byte_t *)currOutTracklet + dSize );
 526     mySize+=dSize;
 527     outPtr->fTrackletCnt++;
 528   }
 529
 530
 531   AliHLTComponentBlockData bd;
 532   FillBlockData( bd );
 533   bd.fOffset = 0;
 534   bd.fSize = mySize;
 535   bd.fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification( slice, slice, minPatch, maxPatch );
 536   outputBlocks.push_back( bd );
 537
 538   size = mySize;
 539
 540   timerReco.Stop();
 541
 542   // Set log level to "Warning" for on-line system monitoring
 543
 544   Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "Tracks",
 545            "CATracker slice %d: output %d tracks;  input %d clusters, patches %d..%d, rows %d..%d; reco time %d/%d us",
 546            slice, ntracks, nClusters, minPatch, maxPatch, row[0], row[1], (Int_t) (timer.RealTime()*1.e6), (Int_t) (timerReco.RealTime()*1.e6) );
 547
 548   return ret;
 549
 550 }
 551
 552