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