[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 "AliHLTTPCCATrackerComponent.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCCATracker.h"
#include "AliHLTTPCCAHit.h"
#include "AliHLTTPCCAOutTrack.h"

#include "AliHLTTPCVertex.h"
#include "AliHLTTPCSpacePointData.h"
#include "AliHLTTPCVertexData.h"
#include "AliHLTTPCClusterDataFormat.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCTrackSegmentData.h"
#include "AliHLTTPCTrackArray.h"
#include "AliHLTTPCTrackletDataFormat.h"
#include "AliHLTTPCDefinitions.h"
#include "TMath.h"
#include "AliTPC.h"
#include "AliTPCParam.h"
#include "AliRun.h"
#include <stdlib.h>
#include <iostream>
#include <errno.h>


// this is a global object used for automatic component registration, do not use this
AliHLTTPCCATrackerComponent gAliHLTTPCCATrackerComponent;

ClassImp(AliHLTTPCCATrackerComponent)

AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent()
  :
  fTracker(NULL),
  fVertex(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&)
  :
  fTracker(NULL),
  fVertex(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
    }

// 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 );
    list.push_back( AliHLTTPCDefinitions::fgkVertexDataType );
    }

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

void AliHLTTPCCATrackerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ) 
    {
    // see header file for class documentation
    // XXX TODO: Find more realistic values.
    constBase = 0;
    inputMultiplier = 0.2;
    }

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

int AliHLTTPCCATrackerComponent::DoInit( int argc, const char** argv )
    {
  // see header file for class documentation

    if ( fTracker || fVertex )
	return EINPROGRESS;

    fTracker = new AliHLTTPCCATracker();
    fVertex = new AliHLTTPCVertex();


/* ---------------------------------------------------------------------------------
 * cmdline arguments not needed so far

    int i = 0;
    char* cpErr;

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

	Logging(kHLTLogError, "HLT::TPCSliceTracker::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;
    if ( fVertex )
	delete fVertex;
    fVertex = NULL;
    return 0;
    }

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

  // see header file for class documentation

  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;
  AliHLTTPCVertexData* inPtrV = NULL;
  const AliHLTComponentBlockData* vertexIter=NULL;
 
 
  AliHLTUInt32_t vSize = 0;
  UInt_t offset=0, tSize = 0;

  // ------------------------------------------
  
  Int_t slice=-1, patch=-1, row[2];
  Int_t minPatch=INT_MAX, maxPatch = 0;
  offset = 0;
  std::vector<Int_t> slices;
  std::vector<Int_t>::iterator slIter, slEnd;
  std::vector<unsigned> sliceCnts;
  std::vector<unsigned>::iterator slCntIter;
  Int_t vertexSlice=-1;
  
  // Find min/max rows used in total and find and read out vertex if it is present
  // also determine correct slice number, if multiple slice numbers are present in event
  // (which should not happen in the first place) we use the one that occurs the most times
  row[0] = 0;
  row[1] = 0;
  bool found;
  for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ )
    {
      iter = blocks+ndx;
      slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
      found = false;
      slIter = slices.begin();
      slEnd = slices.end();
      slCntIter = sliceCnts.begin();
      while ( slIter != slEnd )
	{
	  if ( *slIter == slice )
	    {
	      found = true;
	      break;
	    }
	  slIter++;
	  slCntIter++;
	}
      if ( !found )
	{
	  slices.insert( slices.end(), slice );
	  sliceCnts.insert( sliceCnts.end(), 1 );
	}
      else
	*slCntIter++;
      
      if ( iter->fDataType == AliHLTTPCDefinitions::fgkVertexDataType )
	{
	  inPtrV = (AliHLTTPCVertexData*)(iter->fPtr);
	  vertexIter = iter;
	  vSize = iter->fSize;
	  fVertex->Read( inPtrV );
	  vertexSlice = slice;
	}
      if ( iter->fDataType == AliHLTTPCDefinitions::fgkClustersDataType )
	{
	  patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
	  if ( minPatch>patch )
	    {
	      minPatch = patch;
	      row[0] = AliHLTTPCTransform::GetFirstRow( patch );
	    }
	  if ( maxPatch<patch )
	    {
	      maxPatch = patch;
	      row[1] = AliHLTTPCTransform::GetLastRow( patch );
	    }
	}
    }
  
  // 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;
      slIter = slices.begin();
      slEnd = slices.end();
      slCntIter = sliceCnts.begin();
      while ( slIter != slEnd )
	{
	  Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
		   "Slice %lu found %lu times.", *slIter, *slCntIter );
	  if ( maxCntSlice<*slCntIter )
	    {
	      maxCntSlice = *slCntIter;
	      slice = *slIter;
	    }
	  slIter++;
	  slCntIter++;
	}
      Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
	       "Using slice %lu.", slice );
    }
  else if ( slices.size()>0 )
    {
      slice = *(slices.begin());
    }
  else
    {
      slice = -1;
    }
  
    
  if ( vertexSlice != slice )
    {
      // multiple vertex blocks in event and we used the wrong one...
      found = false;
      for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ )
	{
	  iter = blocks+ndx;
	  if ( iter->fDataType == AliHLTTPCDefinitions::fgkVertexDataType && slice==AliHLTTPCDefinitions::GetMinSliceNr( *iter ) )
	    {
	      inPtrV = (AliHLTTPCVertexData*)(iter->fPtr);
	      vertexIter = iter;
	      vSize = iter->fSize;
	      fVertex->Read( inPtrV );
	      break;
	    }
	}
    }
    
  // read in all hits
  std::vector<unsigned long> patchIndices;
  std::vector<unsigned long>::iterator pIter, pEnd;
  for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ )
    {
      iter = blocks+ndx;
      
      if ( iter->fDataType == AliHLTTPCDefinitions::fgkClustersDataType && slice==AliHLTTPCDefinitions::GetMinSliceNr( *iter ) )
	{
	  patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
	  pIter = patchIndices.begin();
	  pEnd = patchIndices.end();
	  while ( pIter!=pEnd && AliHLTTPCDefinitions::GetMinSliceNr( blocks[*pIter] ) < patch )
	    pIter++;
	  patchIndices.insert( pIter, ndx );
	}
    }
       
  // Initialize tracker
  Double_t Bz = -5;

  {
    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|| (iSec>=36&&iSec<54) );
    Bool_t rInner = (iSec<36);
    Double_t zMin =  zPlus ?plusZmin :minusZmin;
    Double_t zMax =  zPlus ?plusZmax :minusZmax;
    Double_t rMin =  rInner ?inRmin :outRmin;
    Double_t rMax =  rInner ?inRmax :outRmax;
    Int_t inNRows = 63;
    Int_t outNRows = 96;
    Double_t inRowXFirst = 85.225;
    Double_t outRowXFirst =135.1;
    Double_t inRowXStep = 0.75;
    Double_t outRowXStep = 1.;
    Int_t nRows = rInner ?inNRows :outNRows;
    Double_t rowXFirst = rInner ?inRowXFirst :outRowXFirst;
    Double_t rowXStep = rInner ?inRowXStep :outRowXStep;
    
    Int_t nSectors = 72/2;
    
    Double_t padPitch = 0.4;
    Double_t sigmaZ = 0.228808;
    
    //TPCZmin = -249.645, ZMax = 249.778
    
    if(0){
      if( !gAlice ) return 0;
      AliTPC *tpc = (AliTPC*) gAlice->GetDetector("TPC");      
      AliTPCParam *param = tpc->GetParam();
      cout<<" R inner = "<<param->GetInnerRadiusLow()<<" "<<param->GetInnerRadiusUp()<<endl;
      cout<<" R outer = "<<param->GetOuterRadiusLow()<<" "<<param->GetOuterRadiusUp()<<endl;
      cout<<" Pitch = "<<param->GetPadPitchWidth(0)<<endl;
      cout<<" Sigma Z = "<<param->GetZSigma()<<endl;
      nSectors = param->GetNSector();      
      if( iSec<0 || iSec >= nSectors ) return 0;      
      
      padPitch = param->GetPadPitchWidth(iSec);
      sigmaZ = param->GetZSigma();      
      alpha = param->GetAngle(iSec);      
      
      if( iSec<param->GetNInnerSector() ){
	dalpha = param->GetInnerAngle();
	rMin = param->GetInnerRadiusLow();
	rMax = param->GetInnerRadiusUp();
      } else {  
	dalpha = param->GetOuterAngle();
	rMin = param->GetOuterRadiusLow();
	rMax = param->GetOuterRadiusUp();
      }
      
      TGeoHMatrix  *mat = param->GetClusterMatrix(iSec);
      Double_t p0[3]={0, 0, 0 };
      Double_t p1[3]={0, 0, param->GetZLength(iSec) };
      Double_t p0C[3], p1C[3];
      mat->LocalToMaster(p0,p0C);
      mat->LocalToMaster(p1,p1C);
      Int_t iZ = (iSec%36)/18;
      if( iZ==0 ){
	zMin = p0C[2]; // plus Z
	zMax = p1C[2];
      } else {
	zMin = -p1C[2]; // minus Z
	zMax = p0C[2];      
      }
    }
     
    AliHLTTPCCAParam param;
    param.Initialize( iSec, inNRows+outNRows, inRowXFirst, inRowXStep,alpha, dalpha,
		      inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, Bz );
      
    fTracker->Initialize( param );
 
    for( Int_t irow=0; irow<outNRows; irow++){
      fTracker->Rows()[inNRows+irow].X() = outRowXFirst + irow*outRowXStep;
    }
     
  }

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

  Int_t nHitsTotal = 0;
  pIter = patchIndices.begin();
  pEnd = patchIndices.end();
  while ( pIter!=pEnd ){
    ndx = *pIter;
    iter = blocks+ndx;	
    inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
    nHitsTotal+=inPtrSP->fSpacePointCnt;
    pIter++;
  }
 
  AliHLTTPCCAHit *vHits = new AliHLTTPCCAHit[nHitsTotal];
  Double_t *vHitStore = new Double_t [nHitsTotal];
  Int_t nHits = 0;
    
  pIter = patchIndices.begin();
  while ( pIter!=pEnd )
    {
      ndx = *pIter;
      iter = blocks+ndx;
      
      patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
      inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
      
      Logging( kHLTLogDebug, "HLT::TPCSliceTracker::DoEvent", "Reading hits",
	       "Reading hits for slice %d - patch %d", slice, patch );
      
      // Read patch hits

      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]);
	  //Logging( kHLTLogDebug, "HLT::TPCSliceTracker::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) );
	  if( pSP->fPadRow != oldRow ){
	    if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
	    oldRow = pSP->fPadRow;
	    firstRowHit = nHits;
	    nRowHits = 0;
	  }
	  AliHLTTPCCAHit &h = vHits[nHits];
	  h.Y() = pSP->fY;
	  h.Z() = pSP->fZ;
	  h.ErrY() = TMath::Sqrt(pSP->fSigmaY2);
	  h.ErrZ() = TMath::Sqrt(pSP->fSigmaZ2);  
	  h.ID() = pSP->fID;
	  vHitStore[nHits] = pSP->fX;
	  nHits++;	
	  nRowHits++;
	}	
      if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
      pIter++;
    }
  
  // reconstruct the event  

  fTracker->Reconstruct();


  // write reconstructed tracks

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

  AliHLTTPCTrackSegmentData* currOutTracklet = outPtr->fTracklets;

  Int_t ntracks = fTracker->NOutTracks();
  
  for( int itr=0; itr<ntracks; itr++ ){
    
    AliHLTTPCCAOutTrack &t = fTracker->OutTracks()[itr];    
    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, vHitStore[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(ex)>1.e-5  ? ez/ex :1.e5;
      
    h3 = (TMath::Abs(ex) >1.e-5) ? -ez/ex/ex :0;
    h5 = (TMath::Abs(ex) >1.e-5) ? 1./ex :0;
    currOutTracklet->fTglerr =  h3*h3*t.Param().Cov()[9] + h5*h5*t.Param().Cov()[20] + 2.*h3*h5*t.Param().Cov()[18]; 
    
    currOutTracklet->fCharge = -currOutTracklet->fCharge;
    t.Param().TransportBz(Bz, vHitStore[iLastHit], lastHit.Y(), lastHit.Z() );
    currOutTracklet->fLastX = t.Param().Par()[0];
    currOutTracklet->fLastY = t.Param().Par()[1];
    currOutTracklet->fLastZ = t.Param().Par()[2];

    currOutTracklet->fNPoints = t.NHits();

    for( Int_t i=0; i<t.NHits(); i++ ){
      currOutTracklet->fPointIDs[i] = fTracker->OutTrackHits()[t.FirstHitRef()+i];
    }
    
    Byte_t *tmpP = (Byte_t *)currOutTracklet;
    
    tmpP += sizeof(AliHLTTPCTrackSegmentData) + currOutTracklet->fNPoints*sizeof(UInt_t);
    currOutTracklet = (AliHLTTPCTrackSegmentData*)tmpP;
  }
  
  outPtr->fTrackletCnt = ntracks; 
  
  delete[] vHits;
  delete[] vHitStore;
  
  Logging( kHLTLogDebug, "HLT::TPCSliceTracker::DoEvent", "Tracks",
	   "Input: Number of tracks: %lu Slice/MinPatch/MaxPatch/RowMin/RowMax: %lu/%lu/%lu/%lu/%lu.", 
	   ntracks, slice, minPatch, maxPatch, row[0], row[1] );
  
  AliHLTUInt8_t *pbeg = (AliHLTUInt8_t *)outputPtr;
  AliHLTUInt8_t *pend = (AliHLTUInt8_t *)currOutTracklet;
  UInt_t mySize = pend - pbeg;
  
  AliHLTComponentBlockData bd;
  FillBlockData( bd );
  bd.fOffset = offset;
  bd.fSize = mySize;
  bd.fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification( slice, slice, minPatch, maxPatch );      
  outputBlocks.push_back( bd );
  
#ifdef FORWARD_VERTEX_BLOCK
  if ( vertexIter )
    {
      // Copy the descriptor block for the vertex information.
      //bd = *vertexIter;
      //outputBlocks.push_back( bd );
    }
#endif // FORWARD_VERTEX_BLOCK
  
  size = mySize;
  
  return 0;
}
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
	29	#include "AliHLTTPCCATrackerComponent.h"
	30	#include "AliHLTTPCTransform.h"
	31	#include "AliHLTTPCCATracker.h"
	32	#include "AliHLTTPCCAHit.h"
	33	#include "AliHLTTPCCAOutTrack.h"
	34
	35	#include "AliHLTTPCVertex.h"
	36	#include "AliHLTTPCSpacePointData.h"
	37	#include "AliHLTTPCVertexData.h"
	38	#include "AliHLTTPCClusterDataFormat.h"
	39	#include "AliHLTTPCTransform.h"
	40	#include "AliHLTTPCTrackSegmentData.h"
	41	#include "AliHLTTPCTrackArray.h"
	42	#include "AliHLTTPCTrackletDataFormat.h"
	43	#include "AliHLTTPCDefinitions.h"
	44	#include "TMath.h"
	45	#include "AliTPC.h"
	46	#include "AliTPCParam.h"
	47	#include "AliRun.h"
	48	#include <stdlib.h>
	49	#include <iostream>
	50	#include <errno.h>
	51
	52
	53	// this is a global object used for automatic component registration, do not use this
	54	AliHLTTPCCATrackerComponent gAliHLTTPCCATrackerComponent;
	55
	56	ClassImp(AliHLTTPCCATrackerComponent)
	57
	58	AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent()
	59	:
	60	fTracker(NULL),
	61	fVertex(NULL),
	62	fBField(0)
	63	{
	64	// see header file for class documentation
65	// or
66	// refer to README to build package
67	// or
68	// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
69	}
70
71	AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent(const AliHLTTPCCATrackerComponent&)
72	:
73	fTracker(NULL),
74	fVertex(NULL),
75	fBField(0)
76	{
77	// see header file for class documentation
78	HLTFatal("copy constructor untested");
79	}
80
81	AliHLTTPCCATrackerComponent& AliHLTTPCCATrackerComponent::operator=(const AliHLTTPCCATrackerComponent&)
82	{
83	// see header file for class documentation
84	HLTFatal("assignment operator untested");
85	return *this;
86	}
87
88	AliHLTTPCCATrackerComponent::~AliHLTTPCCATrackerComponent()
89	{
90	// see header file for class documentation
91	}
92
93	// Public functions to implement AliHLTComponent's interface.
94	// These functions are required for the registration process
95
96	const char* AliHLTTPCCATrackerComponent::GetComponentID()
97	{
98	// see header file for class documentation
99	return "TPCCATracker";
100	}
101
102	void AliHLTTPCCATrackerComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list)
103	{
104	// see header file for class documentation
105	list.clear();
106	list.push_back( AliHLTTPCDefinitions::fgkClustersDataType );
107	list.push_back( AliHLTTPCDefinitions::fgkVertexDataType );
108	}
109
110	AliHLTComponentDataType AliHLTTPCCATrackerComponent::GetOutputDataType()
111	{
112	// see header file for class documentation
113	return AliHLTTPCDefinitions::fgkTrackSegmentsDataType;
114	}
115
116	void AliHLTTPCCATrackerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
117	{
118	// see header file for class documentation
119	// XXX TODO: Find more realistic values.
120	constBase = 0;
121	inputMultiplier = 0.2;
122	}
123
124	AliHLTComponent* AliHLTTPCCATrackerComponent::Spawn()
125	{
126	// see header file for class documentation
127	return new AliHLTTPCCATrackerComponent;
128	}
129
130	int AliHLTTPCCATrackerComponent::DoInit( int argc, const char** argv )
131	{
132	// see header file for class documentation
133
134	if ( fTracker \|\| fVertex )
135	return EINPROGRESS;
136
137	fTracker = new AliHLTTPCCATracker();
138	fVertex = new AliHLTTPCVertex();
139
140
141	/* ---------------------------------------------------------------------------------
142	* cmdline arguments not needed so far
143
144	int i = 0;
145	char* cpErr;
146
147	while ( i < argc )
148	{
149	if ( !strcmp( argv[i], "bfield" ) )
150	{
151	if ( argc <= i+1 )
152	{
153	Logging( kHLTLogError, "HLT::TPCSliceTracker::DoInit", "Missing B-field", "Missing B-field specifier." );
154	return ENOTSUP;
155	}
156	fBField = strtod( argv[i+1], &cpErr );
157	if ( *cpErr )
158	{
159	Logging( kHLTLogError, "HLT::TPCSliceTracker::DoInit", "Missing multiplicity", "Cannot convert B-field specifier '%s'.", argv[i+1] );
160	return EINVAL;
161	}
162	i += 2;
163	continue;
164	}
165
166	Logging(kHLTLogError, "HLT::TPCSliceTracker::DoInit", "Unknown Option", "Unknown option '%s'", argv[i] );
167	return EINVAL;
168	}
169	--------------------------------------------------------------------------------- */
170
171	return 0;
172	}
173
174	int AliHLTTPCCATrackerComponent::DoDeinit()
175	{
176	// see header file for class documentation
177	if ( fTracker )
178	delete fTracker;
179	fTracker = NULL;
180	if ( fVertex )
181	delete fVertex;
182	fVertex = NULL;
183	return 0;
184	}
185
186	int AliHLTTPCCATrackerComponent::DoEvent( const AliHLTComponentEventData& evtData, const AliHLTComponentBlockData* blocks,
187	AliHLTComponentTriggerData& trigData, AliHLTUInt8_t* outputPtr,
188	AliHLTUInt32_t& size, vector<AliHLTComponentBlockData>& outputBlocks )
189	{
190
191	// see header file for class documentation
192
193	Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "DoEvent", "DoEvent()" );
194	if ( evtData.fBlockCnt<=0 )
195	{
196	Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "no blocks in event" );
197	return 0;
198	}
199
200	const AliHLTComponentBlockData* iter = NULL;
201	unsigned long ndx;
202	AliHLTTPCClusterData* inPtrSP;
203	AliHLTTPCVertexData* inPtrV = NULL;
204	const AliHLTComponentBlockData* vertexIter=NULL;
205
206
207	AliHLTUInt32_t vSize = 0;
208	UInt_t offset=0, tSize = 0;
209
210	// ------------------------------------------
211
212	Int_t slice=-1, patch=-1, row[2];
213	Int_t minPatch=INT_MAX, maxPatch = 0;
214	offset = 0;
215	std::vector<Int_t> slices;
216	std::vector<Int_t>::iterator slIter, slEnd;
217	std::vector<unsigned> sliceCnts;
218	std::vector<unsigned>::iterator slCntIter;
219	Int_t vertexSlice=-1;
220
221	// Find min/max rows used in total and find and read out vertex if it is present
222	// also determine correct slice number, if multiple slice numbers are present in event
223	// (which should not happen in the first place) we use the one that occurs the most times
224	row[0] = 0;
225	row[1] = 0;
226	bool found;
227	for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ )
228	{
229	iter = blocks+ndx;
230	slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
231	found = false;
232	slIter = slices.begin();
233	slEnd = slices.end();
234	slCntIter = sliceCnts.begin();
235	while ( slIter != slEnd )
236	{
237	if ( *slIter == slice )
238	{
239	found = true;
240	break;
241	}
242	slIter++;
243	slCntIter++;
244	}
245	if ( !found )
246	{
247	slices.insert( slices.end(), slice );
248	sliceCnts.insert( sliceCnts.end(), 1 );
249	}
250	else
251	*slCntIter++;
252
253	if ( iter->fDataType == AliHLTTPCDefinitions::fgkVertexDataType )
254	{
255	inPtrV = (AliHLTTPCVertexData*)(iter->fPtr);
256	vertexIter = iter;
257	vSize = iter->fSize;
258	fVertex->Read( inPtrV );
259	vertexSlice = slice;
260	}
261	if ( iter->fDataType == AliHLTTPCDefinitions::fgkClustersDataType )
262	{
263	patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
264	if ( minPatch>patch )
265	{
266	minPatch = patch;
267	row[0] = AliHLTTPCTransform::GetFirstRow( patch );
268	}
269	if ( maxPatch<patch )
270	{
271	maxPatch = patch;
272	row[1] = AliHLTTPCTransform::GetLastRow( patch );
273	}
274	}
275	}
276
277	// Determine slice number to really use.
278	if ( slices.size()>1 )
279	{
280	Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
281	"Multiple slice numbers found in event 0x%08lX (%lu). Determining maximum occuring slice number...",
282	evtData.fEventID, evtData.fEventID );
283	unsigned maxCntSlice=0;
284	slIter = slices.begin();
285	slEnd = slices.end();
286	slCntIter = sliceCnts.begin();
287	while ( slIter != slEnd )
288	{
289	Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
290	"Slice %lu found %lu times.", slIter, slCntIter );
291	if ( maxCntSlice<*slCntIter )
292	{
293	maxCntSlice = *slCntIter;
294	slice = *slIter;
295	}
296	slIter++;
297	slCntIter++;
298	}
299	Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
300	"Using slice %lu.", slice );
301	}
302	else if ( slices.size()>0 )
303	{
304	slice = *(slices.begin());
305	}
306	else
307	{
308	slice = -1;
309	}
310
311
312	if ( vertexSlice != slice )
313	{
314	// multiple vertex blocks in event and we used the wrong one...
315	found = false;
316	for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ )
317	{
318	iter = blocks+ndx;
319	if ( iter->fDataType == AliHLTTPCDefinitions::fgkVertexDataType && slice==AliHLTTPCDefinitions::GetMinSliceNr( *iter ) )
320	{
321	inPtrV = (AliHLTTPCVertexData*)(iter->fPtr);
322	vertexIter = iter;
323	vSize = iter->fSize;
324	fVertex->Read( inPtrV );
325	break;
326	}
327	}
328	}
329
330	// read in all hits
331	std::vector<unsigned long> patchIndices;
332	std::vector<unsigned long>::iterator pIter, pEnd;
333	for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ )
334	{
335	iter = blocks+ndx;
336
337	if ( iter->fDataType == AliHLTTPCDefinitions::fgkClustersDataType && slice==AliHLTTPCDefinitions::GetMinSliceNr( *iter ) )
338	{
339	patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
340	pIter = patchIndices.begin();
341	pEnd = patchIndices.end();
342	while ( pIter!=pEnd && AliHLTTPCDefinitions::GetMinSliceNr( blocks[*pIter] ) < patch )
343	pIter++;
344	patchIndices.insert( pIter, ndx );
345	}
346	}
347
348	// Initialize tracker
349	Double_t Bz = -5;
350
351	{
352	Int_t iSec = slice;
353	Double_t inRmin = 83.65;
354	Double_t inRmax = 133.3;
355	Double_t outRmin = 133.5;
356	Double_t outRmax = 247.7;
357	Double_t plusZmin = 0.0529937;
358	Double_t plusZmax = 249.778;
359	Double_t minusZmin = -249.645;
360	Double_t minusZmax = -0.0799937;
361	Double_t dalpha = 0.349066;
362	Double_t alpha = 0.174533 + dalpha*iSec;
363
364	Bool_t zPlus = (iSec<18\|\| (iSec>=36&&iSec<54) );
365	Bool_t rInner = (iSec<36);
366	Double_t zMin = zPlus ?plusZmin :minusZmin;
367	Double_t zMax = zPlus ?plusZmax :minusZmax;
368	Double_t rMin = rInner ?inRmin :outRmin;
369	Double_t rMax = rInner ?inRmax :outRmax;
370	Int_t inNRows = 63;
371	Int_t outNRows = 96;
372	Double_t inRowXFirst = 85.225;
373	Double_t outRowXFirst =135.1;
374	Double_t inRowXStep = 0.75;
375	Double_t outRowXStep = 1.;
376	Int_t nRows = rInner ?inNRows :outNRows;
377	Double_t rowXFirst = rInner ?inRowXFirst :outRowXFirst;
378	Double_t rowXStep = rInner ?inRowXStep :outRowXStep;
379
380	Int_t nSectors = 72/2;
381
382	Double_t padPitch = 0.4;
383	Double_t sigmaZ = 0.228808;
384
385	//TPCZmin = -249.645, ZMax = 249.778
386
387	if(0){
388	if( !gAlice ) return 0;
389	AliTPC tpc = (AliTPC) gAlice->GetDetector("TPC");
390	AliTPCParam *param = tpc->GetParam();
391	cout<<" R inner = "<<param->GetInnerRadiusLow()<<" "<<param->GetInnerRadiusUp()<<endl;
392	cout<<" R outer = "<<param->GetOuterRadiusLow()<<" "<<param->GetOuterRadiusUp()<<endl;
393	cout<<" Pitch = "<<param->GetPadPitchWidth(0)<<endl;
394	cout<<" Sigma Z = "<<param->GetZSigma()<<endl;
395	nSectors = param->GetNSector();
396	if( iSec<0 \|\| iSec >= nSectors ) return 0;
397
398	padPitch = param->GetPadPitchWidth(iSec);
399	sigmaZ = param->GetZSigma();
400	alpha = param->GetAngle(iSec);
401
402	if( iSec<param->GetNInnerSector() ){
403	dalpha = param->GetInnerAngle();
404	rMin = param->GetInnerRadiusLow();
405	rMax = param->GetInnerRadiusUp();
406	} else {
407	dalpha = param->GetOuterAngle();
408	rMin = param->GetOuterRadiusLow();
409	rMax = param->GetOuterRadiusUp();
410	}
411
412	TGeoHMatrix *mat = param->GetClusterMatrix(iSec);
413	Double_t p0[3]={0, 0, 0 };
414	Double_t p1[3]={0, 0, param->GetZLength(iSec) };
415	Double_t p0C[3], p1C[3];
416	mat->LocalToMaster(p0,p0C);
417	mat->LocalToMaster(p1,p1C);
418	Int_t iZ = (iSec%36)/18;
419	if( iZ==0 ){
420	zMin = p0C[2]; // plus Z
421	zMax = p1C[2];
422	} else {
423	zMin = -p1C[2]; // minus Z
424	zMax = p0C[2];
425	}
426	}
427
428	AliHLTTPCCAParam param;
429	param.Initialize( iSec, inNRows+outNRows, inRowXFirst, inRowXStep,alpha, dalpha,
430	inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, Bz );
431
432	fTracker->Initialize( param );
433
434	for( Int_t irow=0; irow<outNRows; irow++){
435	fTracker->Rows()[inNRows+irow].X() = outRowXFirst + irow*outRowXStep;
436	}
437
438	}
439
440	// pass event to CA Tracker
441
442	fTracker->StartEvent();
443
444	Int_t nHitsTotal = 0;
445	pIter = patchIndices.begin();
446	pEnd = patchIndices.end();
447	while ( pIter!=pEnd ){
448	ndx = *pIter;
449	iter = blocks+ndx;
450	inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
451	nHitsTotal+=inPtrSP->fSpacePointCnt;
452	pIter++;
453	}
454
455	AliHLTTPCCAHit *vHits = new AliHLTTPCCAHit[nHitsTotal];
456	Double_t *vHitStore = new Double_t [nHitsTotal];
457	Int_t nHits = 0;
458
459	pIter = patchIndices.begin();
460	while ( pIter!=pEnd )
461	{
462	ndx = *pIter;
463	iter = blocks+ndx;
464
465	patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
466	inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
467
468	Logging( kHLTLogDebug, "HLT::TPCSliceTracker::DoEvent", "Reading hits",
469	"Reading hits for slice %d - patch %d", slice, patch );
470
471	// Read patch hits
472
473	Int_t oldRow = -1;
474	Int_t nRowHits = 0;
475	Int_t firstRowHit = 0;
476	for (UInt_t i=0; i<inPtrSP->fSpacePointCnt; i++ )
477	{
478	AliHLTTPCSpacePointData* pSP = &(inPtrSP->fSpacePoints[i]);
479	//Logging( kHLTLogDebug, "HLT::TPCSliceTracker::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) );
480	if( pSP->fPadRow != oldRow ){
481	if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
482	oldRow = pSP->fPadRow;
483	firstRowHit = nHits;
484	nRowHits = 0;
485	}
486	AliHLTTPCCAHit &h = vHits[nHits];
487	h.Y() = pSP->fY;
488	h.Z() = pSP->fZ;
489	h.ErrY() = TMath::Sqrt(pSP->fSigmaY2);
490	h.ErrZ() = TMath::Sqrt(pSP->fSigmaZ2);
491	h.ID() = pSP->fID;
492	vHitStore[nHits] = pSP->fX;
493	nHits++;
494	nRowHits++;
495	}
496	if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
497	pIter++;
498	}
499
500	// reconstruct the event
501
502	fTracker->Reconstruct();
503
504
505	// write reconstructed tracks
506
507	AliHLTTPCTrackletData* outPtr = (AliHLTTPCTrackletData*)(outputPtr);
508
509	AliHLTTPCTrackSegmentData* currOutTracklet = outPtr->fTracklets;
510
511	Int_t ntracks = fTracker->NOutTracks();
512
513	for( int itr=0; itr<ntracks; itr++ ){
514
515	AliHLTTPCCAOutTrack &t = fTracker->OutTracks()[itr];
516	Int_t iFirstHit = fTracker->OutTrackHits()[t.FirstHitRef()];
517	Int_t iLastHit = fTracker->OutTrackHits()[t.FirstHitRef()+t.NHits()-1];
518	AliHLTTPCCAHit &firstHit = vHits[iFirstHit];
519	AliHLTTPCCAHit &lastHit = vHits[iLastHit];
520
521	t.Param().TransportBz(Bz, vHitStore[iFirstHit], firstHit.Y(), firstHit.Z() );
522	currOutTracklet->fX = t.Param().Par()[0];
523	currOutTracklet->fY = t.Param().Par()[1];
524	currOutTracklet->fZ = t.Param().Par()[2];
525	Double_t qp = t.Param().Par()[6];
526	Double_t p = TMath::Abs(qp)>1.e-5 ?1./TMath::Abs(qp) :1.e5;
527	Double_t ex = t.Param().Par()[3];
528	Double_t ey = t.Param().Par()[4];
529	Double_t ez = t.Param().Par()[5];
530	Double_t et = TMath::Sqrt( exex + eyey );
531	currOutTracklet->fCharge = (qp>0) ?+1 :(qp<0 ?-1 :0);
532	currOutTracklet->fPt = p*et;
533
534	Double_t h3 = TMath::Abs(ex) >1.e-5 ? p*ex/et :0;
535	Double_t h4 = TMath::Abs(ey) >1.e-5 ? p*ey/et :0;
536	Double_t h5;
537	Double_t h6 = - currOutTracklet->fCharge * p * currOutTracklet->fPt;
538
539	currOutTracklet->fPterr = ( h3h3t.Param().Cov()[9] + h4h4t.Param().Cov()[14] + h6h6t.Param().Cov()[27]
540	+ 2.(h3h4t.Param().Cov()[13]+h3h6t.Param().Cov()[24]+h4h6*t.Param().Cov()[25] )
541	);
542
543	currOutTracklet->fPsi = TMath::ATan2(ey, ex);
544
545	h3 = ex/(et*et);
546	h4 = -ey/(et*et);
547	currOutTracklet->fPsierr = h3h3t.Param().Cov()[9] + h4h4t.Param().Cov()[14] + 2.h3h4*t.Param().Cov()[13];
548
549	currOutTracklet->fTgl = TMath::Abs(ex)>1.e-5 ? ez/ex :1.e5;
550
551	h3 = (TMath::Abs(ex) >1.e-5) ? -ez/ex/ex :0;
552	h5 = (TMath::Abs(ex) >1.e-5) ? 1./ex :0;
553	currOutTracklet->fTglerr = h3h3t.Param().Cov()[9] + h5h5t.Param().Cov()[20] + 2.h3h5*t.Param().Cov()[18];
554
555	currOutTracklet->fCharge = -currOutTracklet->fCharge;
556	t.Param().TransportBz(Bz, vHitStore[iLastHit], lastHit.Y(), lastHit.Z() );
557	currOutTracklet->fLastX = t.Param().Par()[0];
558	currOutTracklet->fLastY = t.Param().Par()[1];
559	currOutTracklet->fLastZ = t.Param().Par()[2];
560
561	currOutTracklet->fNPoints = t.NHits();
562
563	for( Int_t i=0; i<t.NHits(); i++ ){
564	currOutTracklet->fPointIDs[i] = fTracker->OutTrackHits()[t.FirstHitRef()+i];
565	}
566
567	Byte_t tmpP = (Byte_t )currOutTracklet;
568
569	tmpP += sizeof(AliHLTTPCTrackSegmentData) + currOutTracklet->fNPoints*sizeof(UInt_t);
570	currOutTracklet = (AliHLTTPCTrackSegmentData*)tmpP;
571	}
572
573	outPtr->fTrackletCnt = ntracks;
574
575	delete[] vHits;
576	delete[] vHitStore;
577
578	Logging( kHLTLogDebug, "HLT::TPCSliceTracker::DoEvent", "Tracks",
579	"Input: Number of tracks: %lu Slice/MinPatch/MaxPatch/RowMin/RowMax: %lu/%lu/%lu/%lu/%lu.",
580	ntracks, slice, minPatch, maxPatch, row[0], row[1] );
581
582	AliHLTUInt8_t pbeg = (AliHLTUInt8_t )outputPtr;
583	AliHLTUInt8_t pend = (AliHLTUInt8_t )currOutTracklet;
584	UInt_t mySize = pend - pbeg;
585
586	AliHLTComponentBlockData bd;
587	FillBlockData( bd );
588	bd.fOffset = offset;
589	bd.fSize = mySize;
590	bd.fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification( slice, slice, minPatch, maxPatch );
591	outputBlocks.push_back( bd );
592
593	#ifdef FORWARD_VERTEX_BLOCK
594	if ( vertexIter )
595	{
596	// Copy the descriptor block for the vertex information.
597	//bd = *vertexIter;
598	//outputBlocks.push_back( bd );
599	}
600	#endif // FORWARD_VERTEX_BLOCK
601
602	size = mySize;
603
604	return 0;
605	}
606
607