[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: Jochen Thaeder <thaeder@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 "AliHLTTPCSpacePointData.h"
#include "AliHLTTPCClusterDataFormat.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCTrackSegmentData.h"
#include "AliHLTTPCTrackArray.h"
#include "AliHLTTPCTrackletDataFormat.h"
#include "AliHLTTPCDefinitions.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)
{
  // 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)
{
  // 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;
  
  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( kHLTLogDebug, "HLT::TPCCATracker::DoInit", "Reading command line",
	       "Magnetic field value is set to %f kG", fBField );

      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( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "DoEvent", "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", "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[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 ); 

  }

    
  // 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[nHitsTotal]; // CA hit array
  Double_t vHitStoreX[nHitsTotal];       // hit X coordinates
  Int_t vHitStoreID[nHitsTotal];            // hit ID's
  Int_t vHitRowID[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;
      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());

    // 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];
    
    AliHLTTPCCAHit &firstHit = vHits[iFirstHit];
    AliHLTTPCCAHit &lastHit = vHits[iLastHit];
   
    t.Param().TransportBz(Bz, vHitStoreX[iFirstHit], firstHit.Y(), firstHit.Z() );

    currOutTracklet->fX = t.Param().Par()[0];
    currOutTracklet->fY = t.Param().Par()[1];
    currOutTracklet->fZ = t.Param().Par()[2];
    Double_t qp = t.Param().Par()[6];
    Double_t p = TMath::Abs(qp)>1.e-5 ?1./TMath::Abs(qp) :1.e5;
    Double_t ex = t.Param().Par()[3];
    Double_t ey = t.Param().Par()[4];
    Double_t ez = t.Param().Par()[5];
    Double_t et = TMath::Sqrt( ex*ex + ey*ey );
    
    currOutTracklet->fCharge = (qp>0) ?+1 :(qp<0 ?-1 :0);
    currOutTracklet->fPt = p*et;

    Double_t h3 =  TMath::Abs(ex) >1.e-5 ? p*ex/et :0;
    Double_t h4 =  TMath::Abs(ey) >1.e-5 ? p*ey/et :0;
    Double_t h5;
    Double_t h6 =  - currOutTracklet->fCharge * p * currOutTracklet->fPt;
  
    currOutTracklet->fPterr = ( h3*h3*t.Param().Cov()[9] + h4*h4*t.Param().Cov()[14] + h6*h6*t.Param().Cov()[27] 
				+ 2.*(h3*h4*t.Param().Cov()[13]+h3*h6*t.Param().Cov()[24]+h4*h6*t.Param().Cov()[25] )
				);    
    currOutTracklet->fPsi = TMath::ATan2(ey, ex);
    
    h3 =  ex/(et*et);
    h4 = -ey/(et*et);
    currOutTracklet->fPsierr = h3*h3*t.Param().Cov()[9] + h4*h4*t.Param().Cov()[14] + 2.*h3*h4*t.Param().Cov()[13];
    
    currOutTracklet->fTgl = TMath::Abs(et)>1.e-5  ? ez/et :1.e5;
      
    if( TMath::Abs(et) >1.e-2 ){
      h3 = -ez*ex/et/et;
      h4 = -ez*ey/et/et;
      h5 = 1.;
      currOutTracklet->fTglerr =  ( h3*h3*t.Param().Cov()[9] + h4*h4*t.Param().Cov()[14] + h5*h5*t.Param().Cov()[20] 
				    + 2.*(h3*h4*t.Param().Cov()[13]+h3*h5*t.Param().Cov()[18]+h4*h5*t.Param().Cov()[19] )
				    )/et/et;
    }else{
      currOutTracklet->fTglerr =  1.e-2;
    }

    currOutTracklet->fCharge = -currOutTracklet->fCharge;
    
    t.Param().TransportBz(Bz, vHitStoreX[iLastHit], lastHit.Y(), lastHit.Z() );
     
    currOutTracklet->fLastX = t.Param().Par()[0];
    currOutTracklet->fLastY = t.Param().Par()[1];
    currOutTracklet->fLastZ = t.Param().Par()[2];

#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++; 
  }
  
  
  AliHLTComponentBlockData bd;
  FillBlockData( bd );
  bd.fOffset = 0;
  bd.fSize = mySize;
  bd.fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification( slice, slice, minPatch, maxPatch );      
  outputBlocks.push_back( bd );
  
  size = mySize;
  
  timerReco.Stop();

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

  Logging( kHLTLogWarning, "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) (timer.RealTime()*1.e6), (Int_t) (timerReco.RealTime()*1.e6) );

  return ret;

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