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