[u/mrichter/AliRoot.git] / HLT / TPCLib / tracking-ca / AliHLTTPCCATrackerComponent.cxx

// @(#) $Id$
//***************************************************************************
// This file is property of and copyright by the ALICE HLT Project          * 
// ALICE Experiment at CERN, All rights reserved.                           *
//                                                                          *
// Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
//                  Ivan Kisel <kisel@kip.uni-heidelberg.de>                *
//                  for The ALICE HLT Project.                              *
//                                                                          *
// Permission to use, copy, modify and distribute this software and its     *
// documentation strictly for non-commercial purposes is hereby granted     *
// without fee, provided that the above copyright notice appears in all     *
// copies and that both the copyright notice and this permission notice     *
// appear in the supporting documentation. The authors make no claims       *
// about the suitability of this software for any purpose. It is            *
// provided "as is" without express or implied warranty.                    *
//***************************************************************************
 
///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// a TPC tracker processing component for the HLT based on CA by Ivan Kisel  //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#if __GNUC__>= 3
using namespace std;
#endif

#include <climits>
#include "AliHLTTPCCATrackerComponent.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCCATracker.h"
#include "AliHLTTPCCAHit.h"
#include "AliHLTTPCCAOutTrack.h"
#include "AliHLTTPCCAParam.h"

#include "AliHLTTPCSpacePointData.h"
#include "AliHLTTPCClusterDataFormat.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCTrackSegmentData.h"
#include "AliHLTTPCTrackArray.h"
#include "AliHLTTPCTrackletDataFormat.h"
#include "AliHLTTPCDefinitions.h"
#include "AliExternalTrackParam.h"
#include "TStopwatch.h"
#include "TMath.h"

/** ROOT macro for the implementation of ROOT specific class methods */
ClassImp(AliHLTTPCCATrackerComponent)

AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent()
  :
  fTracker(NULL),
  fBField(0),
  fMinNTrackClusters(30),
  fFullTime(0),
  fRecoTime(0),
  fNEvents(0)
{
  // see header file for class documentation
  // or
  // refer to README to build package
  // or
  // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
}

AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent(const AliHLTTPCCATrackerComponent&)
  :
  AliHLTProcessor(),
  fTracker(NULL),
  fBField(0),
  fMinNTrackClusters(30),
  fFullTime(0),
  fRecoTime(0),
  fNEvents(0)
{
  // see header file for class documentation
  HLTFatal("copy constructor untested");
}

AliHLTTPCCATrackerComponent& AliHLTTPCCATrackerComponent::operator=(const AliHLTTPCCATrackerComponent&)
{
  // see header file for class documentation
  HLTFatal("assignment operator untested");
  return *this;
}

AliHLTTPCCATrackerComponent::~AliHLTTPCCATrackerComponent()
{
  // see header file for class documentation
  delete fTracker;
}

//
// Public functions to implement AliHLTComponent's interface.
// These functions are required for the registration process
//

const char* AliHLTTPCCATrackerComponent::GetComponentID() 
{
  // see header file for class documentation
  return "TPCCATracker";
}

void AliHLTTPCCATrackerComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list) 
{
  // see header file for class documentation
  list.clear();
  list.push_back( AliHLTTPCDefinitions::fgkClustersDataType );
}

AliHLTComponentDataType AliHLTTPCCATrackerComponent::GetOutputDataType() 
{
  // see header file for class documentation
  return AliHLTTPCDefinitions::fgkTrackSegmentsDataType;
}

void AliHLTTPCCATrackerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ) 
{
  // define guess for the output data size
  constBase = 200;       // minimum size
  inputMultiplier = 0.5; // size relative to input 
}

AliHLTComponent* AliHLTTPCCATrackerComponent::Spawn() 
{
  // see header file for class documentation
  return new AliHLTTPCCATrackerComponent;
}

int AliHLTTPCCATrackerComponent::DoInit( int argc, const char** argv )
{
  // Initialize the CA tracker component 
  //
  // arguments could be:
  // bfield - the magnetic field value
  //

  if ( fTracker ) return EINPROGRESS;

  fFullTime = 0;
  fRecoTime = 0;
  fNEvents = 0;

  fTracker = new AliHLTTPCCATracker();
  
  // read command line

  int i = 0;
  char* cpErr;
  while ( i < argc ){
    if ( !strcmp( argv[i], "bfield" ) ){
      if ( i+1 >= argc )
	{
	  Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing B-field", "Missing B-field specifier." );
	  return ENOTSUP;
	}
      fBField = strtod( argv[i+1], &cpErr );
      if ( *cpErr )
	{
	  Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert B-field specifier '%s'.", argv[i+1] );
	  return EINVAL;
	}

      Logging( kHLTLogInfo, "HLT::TPCCATracker::DoInit", "Reading command line",
	       "Magnetic field value is set to %f kG", fBField );

      i += 2;
      continue;
    }

    if ( !strcmp( argv[i], "MinNTrackClusters" ) ){
      if ( i+1 >= argc )
	{
	  Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing MinNTrackClusters", "Missing MinNTrackClusters specifier." );
	  return ENOTSUP;
	}
      fMinNTrackClusters = (Int_t ) strtod( argv[i+1], &cpErr );
      if ( *cpErr )
	{
	  Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert MinNTrackClusters '%s'.", argv[i+1] );
	  return EINVAL;
	}

      Logging( kHLTLogInfo, "HLT::TPCCATracker::DoInit", "Reading command line",
	       "MinNTrackClusters is set to %i ", fMinNTrackClusters );

      i += 2;
      continue;
    }
    
    Logging(kHLTLogError, "HLT::TPCCATracker::DoInit", "Unknown Option", "Unknown option '%s'", argv[i] );
    return EINVAL;
  }
  
  return 0;
}

int AliHLTTPCCATrackerComponent::DoDeinit()
{
  // see header file for class documentation
  if ( fTracker ) delete fTracker;
  fTracker = NULL;  
  return 0;
}

int AliHLTTPCCATrackerComponent::DoEvent
( 
 const AliHLTComponentEventData& evtData, 
 const AliHLTComponentBlockData* blocks, 
 AliHLTComponentTriggerData& /*trigData*/, 
 AliHLTUInt8_t* outputPtr, 
 AliHLTUInt32_t& size, 
 vector<AliHLTComponentBlockData>& outputBlocks )
{

  AliHLTUInt32_t MaxBufferSize = size;
  size = 0; // output size

  TStopwatch timer;

  // Event reconstruction in one TPC slice with CA Tracker

  //Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "CA::DoEvent()" );
  if ( evtData.fBlockCnt<=0 )
    {
      Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "no blocks in event" );
      return 0;
    }
  
  const AliHLTComponentBlockData* iter = NULL;
  unsigned long ndx;
  AliHLTTPCClusterData* inPtrSP; 
 
  // Determine the slice number 
  
  Int_t slice=-1;
  {
    std::vector<Int_t> slices;
    std::vector<Int_t>::iterator slIter;
    std::vector<unsigned> sliceCnts;
    std::vector<unsigned>::iterator slCntIter;
  
    for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ ){
      iter = blocks+ndx;
      if ( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType ) continue;

      slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );

      Bool_t found = 0;
      slCntIter = sliceCnts.begin();
      for( slIter = slices.begin(); slIter!=slices.end(); slIter++, slCntIter++ ){
	if ( *slIter == slice ){
	  found = kTRUE;
	  break;
	}
      }
      if ( !found ){
	slices.push_back( slice );
	sliceCnts.push_back( 1 );
      } else *slCntIter++;     
    }
  
    
    // Determine slice number to really use.
    if ( slices.size()>1 )
      {
	Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
		 "Multiple slice numbers found in event 0x%08lX (%lu). Determining maximum occuring slice number...",
		 evtData.fEventID, evtData.fEventID );
	unsigned maxCntSlice=0;
	slCntIter = sliceCnts.begin();
	for( slIter = slices.begin(); slIter != slices.end(); slIter++, slCntIter++ )
	  {
	    Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
		     "Slice %lu found %lu times.", *slIter, *slCntIter );
	    if ( maxCntSlice<*slCntIter )
	      {
		maxCntSlice = *slCntIter;
		slice = *slIter;
	      }
	  }
	Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
		 "Using slice %lu.", slice );
      }
    else if ( slices.size()>0 )
      {
	slice = *(slices.begin());
      }
  }
  
  if( slice<0 ){
    Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "CA:: no slices found in event" );
    return 0;
  }


  // Initialize the tracker

  Double_t Bz = fBField;
  
  {
    if( !fTracker ) fTracker = new AliHLTTPCCATracker;
    Int_t iSec = slice;
    Double_t inRmin = 83.65; 
    //    Double_t inRmax = 133.3;
    //    Double_t outRmin = 133.5; 
    Double_t outRmax = 247.7;
    Double_t plusZmin = 0.0529937; 
    Double_t plusZmax = 249.778; 
    Double_t minusZmin = -249.645; 
    Double_t minusZmax = -0.0799937; 
    Double_t dalpha = 0.349066;
    Double_t alpha = 0.174533 + dalpha*iSec;
    
    Bool_t zPlus = (iSec<18 );
    Double_t zMin =  zPlus ?plusZmin :minusZmin;
    Double_t zMax =  zPlus ?plusZmax :minusZmax;
    //TPCZmin = -249.645, ZMax = 249.778    
    //    Double_t rMin =  inRmin;
    //    Double_t rMax =  outRmax;
    Int_t NRows = AliHLTTPCTransform::GetNRows();
        
    Double_t padPitch = 0.4;
    Double_t sigmaZ = 0.228808;
    
    Double_t *rowX = new Double_t [NRows];
    for( Int_t irow=0; irow<NRows; irow++){
      rowX[irow] = AliHLTTPCTransform::Row2X( irow );
    }     
     
    AliHLTTPCCAParam param;
    param.Initialize( iSec, NRows, rowX, alpha, dalpha,
		      inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, Bz );
    param.YErrorCorrection() = 1;
    param.ZErrorCorrection() = 2;

    fTracker->Initialize( param ); 
    delete[] rowX;
  }

    
  // min and max patch numbers and row numbers

  Int_t row[2] = {0,0};
  Int_t minPatch=INT_MAX, maxPatch = -1;

  // total n Hits

  Int_t nHitsTotal = 0;

  // sort patches

  std::vector<unsigned long> patchIndices;

  for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ ){
    iter = blocks+ndx;      
    if( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType ) continue;
    if( slice!=AliHLTTPCDefinitions::GetMinSliceNr( *iter ) ) continue;
    inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
    nHitsTotal+=inPtrSP->fSpacePointCnt;
    Int_t patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
    if ( minPatch>patch ){
      minPatch = patch;
      row[0] = AliHLTTPCTransform::GetFirstRow( patch );
    }
    if ( maxPatch<patch ){
      maxPatch = patch;
      row[1] = AliHLTTPCTransform::GetLastRow( patch );
    }
    std::vector<unsigned long>::iterator pIter = patchIndices.begin(); 
    while( pIter!=patchIndices.end() && AliHLTTPCDefinitions::GetMinPatchNr( blocks[*pIter] ) < patch ){
      pIter++;
    }
    patchIndices.insert( pIter, ndx );
  }
           

  // pass event to CA Tracker
  
  fTracker->StartEvent();

  AliHLTTPCCAHit *vHits = new AliHLTTPCCAHit [nHitsTotal]; // CA hit array
  Double_t *vHitStoreX = new Double_t [nHitsTotal];       // hit X coordinates
  Int_t *vHitStoreID = new Int_t [nHitsTotal];            // hit ID's
  Int_t *vHitRowID = new Int_t [nHitsTotal];            // hit ID's

  Int_t nHits = 0;
 
  Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits",
	   "Total %d hits to read for slice %d", nHitsTotal, slice );

  Int_t nClusters=0;

  for( std::vector<unsigned long>::iterator pIter = patchIndices.begin(); pIter!=patchIndices.end(); pIter++ ){
    ndx = *pIter;
    iter = blocks+ndx;
      
    Int_t patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
    inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
      
    Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits",
	     "Reading %d hits for slice %d - patch %d", inPtrSP->fSpacePointCnt, slice, patch );
      
    // Read patch hits, row by row

    Int_t oldRow = -1;
    Int_t nRowHits = 0;
    Int_t firstRowHit = 0;
    for (UInt_t i=0; i<inPtrSP->fSpacePointCnt; i++ ){	
      AliHLTTPCSpacePointData* pSP = &(inPtrSP->fSpacePoints[i]);

      if( pSP->fPadRow != oldRow ){
	if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
	oldRow = pSP->fPadRow;
	firstRowHit = nHits;
	nRowHits = 0;
      }
      AliHLTTPCCAHit &h = vHits[nHits];
      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;

      h.Y() = pSP->fY;
      h.Z() = pSP->fZ;
      if( TMath::Abs(h.Z())>230.) continue;
      h.ErrY() = TMath::Sqrt(TMath::Abs(pSP->fSigmaY2));
      h.ErrZ() = TMath::Sqrt(TMath::Abs(pSP->fSigmaZ2));  
      if( h.ErrY()<.1 ) h.ErrY() = .1;
      if( h.ErrZ()<.1 ) h.ErrZ() = .1;
      if( h.ErrY()>1. ) h.ErrY() = 1.;
      if( h.ErrZ()>1. ) h.ErrZ() = 1.;
      h.ID() = nHits;
      vHitStoreX[nHits] = pSP->fX;
      vHitStoreID[nHits] = pSP->fID;
      vHitRowID[nHits] = pSP->fPadRow;
      nHits++;	
      nRowHits++;
      nClusters++;
    }	
    if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
  }

  // reconstruct the event  

  TStopwatch timerReco;

  fTracker->Reconstruct();

  timerReco.Stop();

  Int_t ret = 0;

  Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reconstruct",
	   "%d tracks found for slice %d",fTracker->NOutTracks(), slice);

  // write reconstructed tracks

  AliHLTTPCTrackletData* outPtr = (AliHLTTPCTrackletData*)(outputPtr);

  AliHLTTPCTrackSegmentData* currOutTracklet = outPtr->fTracklets;

  Int_t ntracks = fTracker->NOutTracks();

  UInt_t mySize =   ((AliHLTUInt8_t *)currOutTracklet) -  ((AliHLTUInt8_t *)outputPtr);

  outPtr->fTrackletCnt = 0; 

  for( int itr=0; itr<ntracks; itr++ ){
    
    AliHLTTPCCAOutTrack &t = fTracker->OutTracks()[itr];    

    //Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Wrtite output","track %d with %d hits", itr, t.NHits());

    if( t.NHits()<fMinNTrackClusters ) continue;

    // calculate output track size

    UInt_t dSize = sizeof(AliHLTTPCTrackSegmentData) + t.NHits()*sizeof(UInt_t);
    
    if( mySize + dSize > MaxBufferSize ){
      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);
      ret = -ENOSPC;
      break;
    }
    
    // convert CA track parameters to HLT Track Segment

    Int_t iFirstHit = fTracker->OutTrackHits()[t.FirstHitRef()];
    Int_t iLastHit = fTracker->OutTrackHits()[t.FirstHitRef()+t.NHits()-1];
    
    AliHLTTPCCATrackParam par = t.StartPoint();

    par.TransportToX( vHitStoreX[iFirstHit] );

    AliExternalTrackParam tp;
    par.GetExtParam( tp, 0, fBField );

    currOutTracklet->fX = tp.GetX();
    currOutTracklet->fY = tp.GetY();
    currOutTracklet->fZ = tp.GetZ();
    currOutTracklet->fCharge = (Int_t ) tp.GetSign();
    currOutTracklet->fPt = TMath::Abs(tp.GetSignedPt());
    Double_t snp =  tp.GetSnp() ;
    if( snp>.999 ) snp=.999;
    if( snp<-.999 ) snp=-.999;
    currOutTracklet->fPsi = TMath::ASin( snp );
    currOutTracklet->fTgl = tp.GetTgl();
    Double_t h = -currOutTracklet->fPt*currOutTracklet->fPt;
    currOutTracklet->fPterr = h*h*tp.GetSigma1Pt2();
    h = 1./TMath::Sqrt(1-snp*snp);
    currOutTracklet->fPsierr = h*h*tp.GetSigmaSnp2();
    currOutTracklet->fTglerr = tp.GetSigmaTgl2();

    par.TransportToX( vHitStoreX[iLastHit] );     
    currOutTracklet->fLastX = par.GetX();
    currOutTracklet->fLastY = par.GetY();
    currOutTracklet->fLastZ = par.GetZ();

#ifdef INCLUDE_TPC_HOUGH
#ifdef ROWHOUGHPARAMS
    currOutTracklet->fTrackID = 0;
    currOutTracklet->fRowRange1 = vHitRowID[iFirstHit];
    currOutTracklet->fRowRange2 = vHitRowID[iLastHit];
    currOutTracklet->fSector = slice;
    currOutTracklet->fPID = 211;
#endif
#endif // INCLUDE_TPC_HOUGH


    currOutTracklet->fNPoints = t.NHits();

    for( Int_t i=0; i<t.NHits(); i++ ){
      currOutTracklet->fPointIDs[i] = vHitStoreID[fTracker->OutTrackHits()[t.FirstHitRef()+i]];
    }

    currOutTracklet = (AliHLTTPCTrackSegmentData*)( (Byte_t *)currOutTracklet + dSize );
    mySize+=dSize;
    outPtr->fTrackletCnt++; 
  }

  delete[] vHits;
  delete[] vHitStoreX;
  delete[] vHitStoreID;
  delete[] vHitRowID;
  
  AliHLTComponentBlockData bd;
  FillBlockData( bd );
  bd.fOffset = 0;
  bd.fSize = mySize;
  bd.fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification( slice, slice, minPatch, maxPatch );      
  outputBlocks.push_back( bd );
  
  size = mySize;
  
  timer.Stop();

  fFullTime+= timer.CpuTime();
  fRecoTime+= timerReco.CpuTime();
  fNEvents++;

  // Set log level to "Warning" for on-line system monitoring

  Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Tracks",
 	   "CATracker slice %d: output %d tracks;  input %d clusters, patches %d..%d, rows %d..%d; reco time %d/%d us", 
	   slice, ntracks, nClusters, minPatch, maxPatch, row[0], row[1], (Int_t) (fFullTime/fNEvents*1.e6), (Int_t) (fRecoTime/fNEvents*1.e6) );

  return ret;

}
Commit	Line	Data
cf471b1e	1	// @(#) $Id$
d54804bf	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	//***************************************************************************
cf471b1e	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
b09247a2	29	#include <climits>
cf471b1e	30	#include "AliHLTTPCCATrackerComponent.h"
	31	#include "AliHLTTPCTransform.h"
	32	#include "AliHLTTPCCATracker.h"
	33	#include "AliHLTTPCCAHit.h"
	34	#include "AliHLTTPCCAOutTrack.h"
d54804bf	35	#include "AliHLTTPCCAParam.h"
cf471b1e	36
cf471b1e	37	#include "AliHLTTPCSpacePointData.h"
cf471b1e	38	#include "AliHLTTPCClusterDataFormat.h"
	39	#include "AliHLTTPCTransform.h"
	40	#include "AliHLTTPCTrackSegmentData.h"
	41	#include "AliHLTTPCTrackArray.h"
	42	#include "AliHLTTPCTrackletDataFormat.h"
	43	#include "AliHLTTPCDefinitions.h"
d54804bf	44	#include "AliExternalTrackParam.h"
dc4788ec	45	#include "TStopwatch.h"
cf471b1e	46	#include "TMath.h"
cf471b1e	47
672f8b8c	48	/** ROOT macro for the implementation of ROOT specific class methods */
cf471b1e	49	ClassImp(AliHLTTPCCATrackerComponent)
	50
	51	AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent()
	52	:
	53	fTracker(NULL),
d54804bf	54	fBField(0),
	55	fMinNTrackClusters(30),
	56	fFullTime(0),
	57	fRecoTime(0),
	58	fNEvents(0)
cf471b1e	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	:
c61a7285	69	AliHLTProcessor(),
cf471b1e	70	fTracker(NULL),
d54804bf	71	fBField(0),
	72	fMinNTrackClusters(30),
	73	fFullTime(0),
	74	fRecoTime(0),
	75	fNEvents(0)
cf471b1e	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()
4c256004	89	{
cf471b1e	90	// see header file for class documentation
4c256004	91	delete fTracker;
4c256004	92	}
cf471b1e	93
4c256004	94	//
cf471b1e	95	// Public functions to implement AliHLTComponent's interface.
cf471b1e	96	// These functions are required for the registration process
4c256004	97	//
cf471b1e	98
cf471b1e	99	const char* AliHLTTPCCATrackerComponent::GetComponentID()
4c256004	100	{
cf471b1e	101	// see header file for class documentation
4c256004	102	return "TPCCATracker";
4c256004	103	}
cf471b1e	104
cf471b1e	105	void AliHLTTPCCATrackerComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list)
4c256004	106	{
cf471b1e	107	// see header file for class documentation
4c256004	108	list.clear();
	109	list.push_back( AliHLTTPCDefinitions::fgkClustersDataType );
	110	}
cf471b1e	111
cf471b1e	112	AliHLTComponentDataType AliHLTTPCCATrackerComponent::GetOutputDataType()
4c256004	113	{
cf471b1e	114	// see header file for class documentation
4c256004	115	return AliHLTTPCDefinitions::fgkTrackSegmentsDataType;
4c256004	116	}
cf471b1e	117
cf471b1e	118	void AliHLTTPCCATrackerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
4c256004	119	{
	120	// define guess for the output data size
	121	constBase = 200; // minimum size
	122	inputMultiplier = 0.5; // size relative to input
	123	}
cf471b1e	124
cf471b1e	125	AliHLTComponent* AliHLTTPCCATrackerComponent::Spawn()
4c256004	126	{
cf471b1e	127	// see header file for class documentation
4c256004	128	return new AliHLTTPCCATrackerComponent;
4c256004	129	}
cf471b1e	130
cf471b1e	131	int AliHLTTPCCATrackerComponent::DoInit( int argc, const char** argv )
4c256004	132	{
	133	// Initialize the CA tracker component
	134	//
	135	// arguments could be:
	136	// bfield - the magnetic field value
	137	//
cf471b1e	138
4c256004	139	if ( fTracker ) return EINPROGRESS;
d54804bf	140
	141	fFullTime = 0;
	142	fRecoTime = 0;
	143	fNEvents = 0;
	144
4c256004	145	fTracker = new AliHLTTPCCATracker();
	146
	147	// read command line
cf471b1e	148
4c256004	149	int i = 0;
	150	char* cpErr;
	151	while ( i < argc ){
	152	if ( !strcmp( argv[i], "bfield" ) ){
	153	if ( i+1 >= argc )
cf471b1e	154	{
4c256004	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;
cf471b1e	163	}
cf471b1e	164
d54804bf	165	Logging( kHLTLogInfo, "HLT::TPCCATracker::DoInit", "Reading command line",
4c256004	166	"Magnetic field value is set to %f kG", fBField );
	167
	168	i += 2;
	169	continue;
cf471b1e	170	}
d54804bf	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	}
4c256004	191
	192	Logging(kHLTLogError, "HLT::TPCCATracker::DoInit", "Unknown Option", "Unknown option '%s'", argv[i] );
	193	return EINVAL;
	194	}
	195
	196	return 0;
	197	}
cf471b1e	198
cf471b1e	199	int AliHLTTPCCATrackerComponent::DoDeinit()
4c256004	200	{
cf471b1e	201	// see header file for class documentation
4c256004	202	if ( fTracker ) delete fTracker;
	203	fTracker = NULL;
	204	return 0;
	205	}
cf471b1e	206
4c256004	207	int AliHLTTPCCATrackerComponent::DoEvent
	208	(
	209	const AliHLTComponentEventData& evtData,
	210	const AliHLTComponentBlockData* blocks,
c61a7285	211	AliHLTComponentTriggerData& /trigData/,
4c256004	212	AliHLTUInt8_t* outputPtr,
	213	AliHLTUInt32_t& size,
	214	vector<AliHLTComponentBlockData>& outputBlocks )
cf471b1e	215	{
dc4788ec	216
	217	AliHLTUInt32_t MaxBufferSize = size;
	218	size = 0; // output size
	219
	220	TStopwatch timer;
	221
4c256004	222	// Event reconstruction in one TPC slice with CA Tracker
cf471b1e	223
d54804bf	224	//Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "CA::DoEvent()" );
cf471b1e	225	if ( evtData.fBlockCnt<=0 )
	226	{
	227	Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "no blocks in event" );
	228	return 0;
	229	}
4c256004	230
cf471b1e	231	const AliHLTComponentBlockData* iter = NULL;
cf471b1e	232	unsigned long ndx;
4c256004	233	AliHLTTPCClusterData* inPtrSP;
cf471b1e	234
4c256004	235	// Determine the slice number
cf471b1e	236
4c256004	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;
cf471b1e	243
4c256004	244	for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ ){
cf471b1e	245	iter = blocks+ndx;
4c256004	246	if ( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType ) continue;
4c256004	247
cf471b1e	248	slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
dc4788ec	249
4c256004	250	Bool_t found = 0;
cf471b1e	251	slCntIter = sliceCnts.begin();
4c256004	252	for( slIter = slices.begin(); slIter!=slices.end(); slIter++, slCntIter++ ){
	253	if ( *slIter == slice ){
	254	found = kTRUE;
	255	break;
cf471b1e	256	}
4c256004	257	}
	258	if ( !found ){
	259	slices.push_back( slice );
	260	sliceCnts.push_back( 1 );
	261	} else *slCntIter++;
cf471b1e	262	}
	263
	264
4c256004	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 ){
d54804bf	293	Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "CA:: no slices found in event" );
4c256004	294	return 0;
4c256004	295	}
cf471b1e	296
4c256004	297
	298	// Initialize the tracker
	299
	300	Double_t Bz = fBField;
	301
cf471b1e	302	{
4c256004	303	if( !fTracker ) fTracker = new AliHLTTPCCATracker;
cf471b1e	304	Int_t iSec = slice;
cf471b1e	305	Double_t inRmin = 83.65;
c61a7285	306	// Double_t inRmax = 133.3;
c61a7285	307	// Double_t outRmin = 133.5;
cf471b1e	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
dc4788ec	316	Bool_t zPlus = (iSec<18 );
cf471b1e	317	Double_t zMin = zPlus ?plusZmin :minusZmin;
cf471b1e	318	Double_t zMax = zPlus ?plusZmax :minusZmax;
dc4788ec	319	//TPCZmin = -249.645, ZMax = 249.778
c61a7285	320	// Double_t rMin = inRmin;
c61a7285	321	// Double_t rMax = outRmax;
dc4788ec	322	Int_t NRows = AliHLTTPCTransform::GetNRows();
dc4788ec	323
cf471b1e	324	Double_t padPitch = 0.4;
	325	Double_t sigmaZ = 0.228808;
	326
d54804bf	327	Double_t *rowX = new Double_t [NRows];
dc4788ec	328	for( Int_t irow=0; irow<NRows; irow++){
	329	rowX[irow] = AliHLTTPCTransform::Row2X( irow );
	330	}
cf471b1e	331
cf471b1e	332	AliHLTTPCCAParam param;
dc4788ec	333	param.Initialize( iSec, NRows, rowX, alpha, dalpha,
cf471b1e	334	inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, Bz );
dc4788ec	335	param.YErrorCorrection() = 1;
dc4788ec	336	param.ZErrorCorrection() = 2;
4c256004	337
dc4788ec	338	fTracker->Initialize( param );
d54804bf	339	delete[] rowX;
cf471b1e	340	}
cf471b1e	341
4c256004	342
	343	// min and max patch numbers and row numbers
	344
	345	Int_t row[2] = {0,0};
dc4788ec	346	Int_t minPatch=INT_MAX, maxPatch = -1;
4c256004	347
4c256004	348	// total n Hits
cf471b1e	349
cf471b1e	350	Int_t nHitsTotal = 0;
4c256004	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;
cf471b1e	360	inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
cf471b1e	361	nHitsTotal+=inPtrSP->fSpacePointCnt;
4c256004	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 );
cf471b1e	376	}
4c256004	377
	378
	379	// pass event to CA Tracker
	380
	381	fTracker->StartEvent();
	382
d54804bf	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
4c256004	387
cf471b1e	388	Int_t nHits = 0;
4c256004	389
	390	Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits",
	391	"Total %d hits to read for slice %d", nHitsTotal, slice );
	392
dc4788ec	393	Int_t nClusters=0;
dc4788ec	394
4c256004	395	for( std::vector<unsigned long>::iterator pIter = patchIndices.begin(); pIter!=patchIndices.end(); pIter++ ){
	396	ndx = *pIter;
	397	iter = blocks+ndx;
cf471b1e	398
4c256004	399	Int_t patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
4c256004	400	inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
cf471b1e	401
4c256004	402	Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits",
4c256004	403	"Reading %d hits for slice %d - patch %d", inPtrSP->fSpacePointCnt, slice, patch );
cf471b1e	404
4c256004	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
4c256004	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];
dc4788ec	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;
dc4788ec	421
4c256004	422	h.Y() = pSP->fY;
4c256004	423	h.Z() = pSP->fZ;
d54804bf	424	if( TMath::Abs(h.Z())>230.) continue;
dc4788ec	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.;
4c256004	431	h.ID() = nHits;
	432	vHitStoreX[nHits] = pSP->fX;
	433	vHitStoreID[nHits] = pSP->fID;
dc4788ec	434	vHitRowID[nHits] = pSP->fPadRow;
4c256004	435	nHits++;
4c256004	436	nRowHits++;
dc4788ec	437	nClusters++;
4c256004	438	}
	439	if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
	440	}
	441
cf471b1e	442	// reconstruct the event
cf471b1e	443
dc4788ec	444	TStopwatch timerReco;
dc4788ec	445
cf471b1e	446	fTracker->Reconstruct();
cf471b1e	447
dc4788ec	448	timerReco.Stop();
dc4788ec	449
4c256004	450	Int_t ret = 0;
	451
	452	Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reconstruct",
	453	"%d tracks found for slice %d",fTracker->NOutTracks(), slice);
cf471b1e	454
	455	// write reconstructed tracks
	456
	457	AliHLTTPCTrackletData* outPtr = (AliHLTTPCTrackletData*)(outputPtr);
	458
	459	AliHLTTPCTrackSegmentData* currOutTracklet = outPtr->fTracklets;
	460
	461	Int_t ntracks = fTracker->NOutTracks();
4c256004	462
	463	UInt_t mySize = ((AliHLTUInt8_t )currOutTracklet) - ((AliHLTUInt8_t )outputPtr);
	464
	465	outPtr->fTrackletCnt = 0;
	466
cf471b1e	467	for( int itr=0; itr<ntracks; itr++ ){
	468
	469	AliHLTTPCCAOutTrack &t = fTracker->OutTracks()[itr];
4c256004	470
	471	//Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Wrtite output","track %d with %d hits", itr, t.NHits());
	472
d54804bf	473	if( t.NHits()<fMinNTrackClusters ) continue;
d54804bf	474
4c256004	475	// calculate output track size
	476
	477	UInt_t dSize = sizeof(AliHLTTPCTrackSegmentData) + t.NHits()*sizeof(UInt_t);
	478
dc4788ec	479	if( mySize + dSize > MaxBufferSize ){
dc4788ec	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);
4c256004	481	ret = -ENOSPC;
	482	break;
	483	}
	484
	485	// convert CA track parameters to HLT Track Segment
	486
cf471b1e	487	Int_t iFirstHit = fTracker->OutTrackHits()[t.FirstHitRef()];
cf471b1e	488	Int_t iLastHit = fTracker->OutTrackHits()[t.FirstHitRef()+t.NHits()-1];
4c256004	489
d54804bf	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;
5b41cd63	504	if( snp<-.999 ) snp=-.999;
d54804bf	505	currOutTracklet->fPsi = TMath::ASin( snp );
	506	currOutTracklet->fTgl = tp.GetTgl();
	507	Double_t h = -currOutTracklet->fPt*currOutTracklet->fPt;
	508	currOutTracklet->fPterr = hhtp.GetSigma1Pt2();
	509	h = 1./TMath::Sqrt(1-snp*snp);
	510	currOutTracklet->fPsierr = hhtp.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();
cf471b1e	517
dc4788ec	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
cf471b1e	529	currOutTracklet->fNPoints = t.NHits();
	530
	531	for( Int_t i=0; i<t.NHits(); i++ ){
4c256004	532	currOutTracklet->fPointIDs[i] = vHitStoreID[fTracker->OutTrackHits()[t.FirstHitRef()+i]];
cf471b1e	533	}
4c256004	534
	535	currOutTracklet = (AliHLTTPCTrackSegmentData)( (Byte_t )currOutTracklet + dSize );
	536	mySize+=dSize;
	537	outPtr->fTrackletCnt++;
cf471b1e	538	}
d54804bf	539
	540	delete[] vHits;
	541	delete[] vHitStoreX;
	542	delete[] vHitStoreID;
	543	delete[] vHitRowID;
cf471b1e	544
cf471b1e	545	AliHLTComponentBlockData bd;
cf471b1e	546	FillBlockData( bd );
4c256004	547	bd.fOffset = 0;
cf471b1e	548	bd.fSize = mySize;
	549	bd.fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification( slice, slice, minPatch, maxPatch );
	550	outputBlocks.push_back( bd );
	551
cf471b1e	552	size = mySize;
cf471b1e	553
d54804bf	554	timer.Stop();
	555
	556	fFullTime+= timer.CpuTime();
	557	fRecoTime+= timerReco.CpuTime();
	558	fNEvents++;
dc4788ec	559
	560	// Set log level to "Warning" for on-line system monitoring
	561
d54804bf	562	Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Tracks",
dc4788ec	563	"CATracker slice %d: output %d tracks; input %d clusters, patches %d..%d, rows %d..%d; reco time %d/%d us",
d54804bf	564	slice, ntracks, nClusters, minPatch, maxPatch, row[0], row[1], (Int_t) (fFullTime/fNEvents1.e6), (Int_t) (fRecoTime/fNEvents1.e6) );
dc4788ec	565
4c256004	566	return ret;
4c256004	567
cf471b1e	568	}
	569
	570