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