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

// @(#) $Id: AliHLTTPCCANeighboursFinder1.cxx 27042 2008-07-02 12:06:02Z richterm $
// **************************************************************************
// 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.                    *
//                                                                          *
//***************************************************************************

#include "AliHLTTPCCANeighboursFinder.h"
#include "AliHLTTPCCAMath.h"
#include "AliHLTTPCCAHitArea.h"
#include "AliHLTTPCCATracker.h"
#include "AliHLTTPCCAHit.h"

//#define DRAW

#ifdef DRAW
#include "AliHLTTPCCADisplay.h"
#endif //DRAW

GPUd() void AliHLTTPCCANeighboursFinder::Thread
( int /*nBlocks*/, int nThreads, int iBlock, int iThread, int iSync,
  AliHLTTPCCASharedMemory &s, AliHLTTPCCATracker &tracker )
{
  //* find neighbours

  if ( iSync == 0 ) {
#ifdef HLTCA_GPUCODE
	for (int i = iThread;i < sizeof(AliHLTTPCCARow) / sizeof(int);i += nThreads)
	{
		reinterpret_cast<int*>(&s.fRow)[i] = reinterpret_cast<int*>(&tracker.SliceDataRows()[iBlock])[i];
		if (iBlock >= 2 && iBlock <= tracker.Param().NRows() - 3)
		{
			reinterpret_cast<int*>(&s.fRowUp)[i] = reinterpret_cast<int*>(&tracker.SliceDataRows()[iBlock + 2])[i];
			reinterpret_cast<int*>(&s.fRowDown)[i] = reinterpret_cast<int*>(&tracker.SliceDataRows()[iBlock - 2])[i];
		}
	}
	__syncthreads();
#endif
    if ( iThread == 0 ) {
      s.fNRows = tracker.Param().NRows();
      s.fIRow = iBlock;
      if ( s.fIRow < s.fNRows ) {
#ifdef HLTCA_GPUCODE
		const AliHLTTPCCARow &row = s.fRow;
#else
		const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
#endif
        s.fNHits = row.NHits();

        if ( ( s.fIRow >= 2 ) && ( s.fIRow <= s.fNRows - 3 ) ) {
          s.fIRowUp = s.fIRow + 2;
          s.fIRowDn = s.fIRow - 2;

          // references to the rows above and below

#ifdef HLTCA_GPUCODE
          const AliHLTTPCCARow &rowUp = s.fRowUp;
          const AliHLTTPCCARow &rowDn = s.fRowDown;
#else
          const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
          const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
#endif
          // the axis perpendicular to the rows
          const float xDn = rowDn.X();
          const float x   = row.X();
          const float xUp = rowUp.X();

          // number of hits in rows above and below
          s.fUpNHits = tracker.Row( s.fIRowUp ).NHits();
          s.fDnNHits = tracker.Row( s.fIRowDn ).NHits();

          // distance of the rows (absolute and relative)
          s.fUpDx = xUp - x;
          s.fDnDx = xDn - x;
          s.fUpTx = xUp / x;
          s.fDnTx = xDn / x;
          // UpTx/DnTx is used to move the HitArea such that central events are preferred (i.e. vertices
          // coming from y = 0, z = 0).

          //s.fGridUp = tracker.Row( s.fIRowUp ).Grid();
          //s.fGridDn = tracker.Row( s.fIRowDn ).Grid();
        }
      }
    }
  } else if ( iSync == 1 ) {
    if ( s.fIRow < s.fNRows ) {
      if ( ( s.fIRow <= 1 ) || ( s.fIRow >= s.fNRows - 2 ) ) {
#ifdef HLTCA_GPUCODE
		const AliHLTTPCCARow &row = s.fRow;
#else
		const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
#endif
        for ( int ih = iThread; ih < s.fNHits; ih += nThreads ) {
          tracker.SetHitLinkUpData( row, ih, -1 );
          tracker.SetHitLinkDownData( row, ih, -1 );
        }
      } else {
/*#ifdef HLTCA_GPUCODE
          const AliHLTTPCCARow &rowUp = s.fRowUp;
          const AliHLTTPCCARow &rowDn = s.fRowDown;
#else
          const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
          const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
#endif

        for ( unsigned int ih = iThread; ih < s.fGridUp.N() + s.fGridUp.Ny() + 2; ih += nThreads ) {
          s.fGridContentUp[ih] = tracker.FirstHitInBin( rowUp, ih );
        }
        for ( unsigned int ih = iThread; ih < s.fGridDn.N() + s.fGridDn.Ny() + 2; ih += nThreads ) {
          s.fGridContentDn[ih] = tracker.FirstHitInBin( rowDn, ih );
        }*/
      }
    }
  } else if ( iSync == 2 ) {
    if ( ( s.fIRow <= 1 ) || ( s.fIRow >= s.fNRows - 2 ) ) return;

    float chi2Cut = 3.*3.*4 * ( s.fUpDx * s.fUpDx + s.fDnDx * s.fDnDx );
    const float kAreaSize = tracker.Param().NeighboursSearchArea();
    //float chi2Cut = 3.*3.*(s.fUpDx*s.fUpDx + s.fDnDx*s.fDnDx ); //SG
#define kMaxN 20

#ifdef HLTCA_GPUCODE
		  const AliHLTTPCCARow &row = s.fRow;
          const AliHLTTPCCARow &rowUp = s.fRowUp;
          const AliHLTTPCCARow &rowDn = s.fRowDown;
#else
		  const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
		  const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
          const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
#endif
    const float y0 = row.Grid().YMin();
    const float z0 = row.Grid().ZMin();
    const float stepY = row.HstepY();
    const float stepZ = row.HstepZ();

	for ( int ih = iThread; ih < s.fNHits; ih += nThreads ) {

      int linkUp = -1;
      int linkDn = -1;

      if ( s.fDnNHits > 0 && s.fUpNHits > 0 ) {

       

        // coordinates of the hit in the current row
#if defined(HLTCA_GPU_TEXTURE_FETCHa)
		ushort2 tmpval = tex1Dfetch(gAliTexRefu2, ((char*) tracker.Data().HitData() - tracker.pData()->GPUTextureBase()) / sizeof(ushort2) + row.HitNumberOffset() + ih);
        const float y = y0 + tmpval.x * stepY;
        const float z = z0 + tmpval.y * stepZ;
#else
        const float y = y0 + tracker.HitDataY( row, ih ) * stepY;
        const float z = z0 + tracker.HitDataZ( row, ih ) * stepZ;
#endif

#ifdef FAST_NEIGHBOURS_FINDER

//#define NFDEBUG

#ifdef NFDEBUG
	printf("\nSearching Neighbours for: %f %f\n", y, z);
#endif

	const float y0Up = rowUp.Grid().YMin();
    const float z0Up = rowUp.Grid().ZMin();
    const float stepYUp = rowUp.HstepY();
    const float stepZUp = rowUp.HstepZ();

    const float y0Dn = rowDn.Grid().YMin();
    const float z0Dn = rowDn.Grid().ZMin();
    const float stepYDn = rowDn.HstepY();
    const float stepZDn = rowDn.HstepZ();

		const float yMinUp = y*s.fUpTx - kAreaSize;
		const float yMaxUp = y*s.fUpTx + kAreaSize;
		const float yMinDn = y*s.fDnTx - kAreaSize;
		const float yMaxDn = y*s.fDnTx + kAreaSize;
		const float zMinUp = z*s.fUpTx - kAreaSize;
		const float zMaxUp = z*s.fUpTx + kAreaSize;
		const float zMinDn = z*s.fDnTx - kAreaSize;
		const float zMaxDn = z*s.fDnTx + kAreaSize;

		int bYUpMin, bZUpMin, bYDnMin, bZDnMin;
		rowUp.Grid().GetBin(yMinUp, zMinUp, &bYUpMin, &bZUpMin);
		rowDn.Grid().GetBin(yMinDn, zMinDn, &bYDnMin, &bZDnMin);

		int bYUpMax, bZUpMax, bYDnMax, bZDnMax;
		rowUp.Grid().GetBin(yMaxUp, zMaxUp, &bYUpMax, &bZUpMax);
		rowDn.Grid().GetBin(yMaxDn, zMaxDn, &bYDnMax, &bZDnMax);

		int nYUp = rowUp.Grid().Ny();
		int nYDn = rowDn.Grid().Ny();

		int ihUp = tracker.Data().FirstHitInBin(rowUp, bZUpMin * nYUp + bYUpMin);
		int ihDn = bZDnMax * nYDn + bYDnMax >= rowDn.Grid().N() ? (rowDn.NHits() - 1) : (tracker.Data().FirstHitInBin(rowDn, bZDnMax * nYDn + bYDnMax + 1) - 1);

		int ihUpMax = tracker.Data().FirstHitInBin(rowUp, bZUpMin * nYUp + bYUpMax + 1) - 1;
		int ihDnMin = tracker.Data().FirstHitInBin(rowDn, bZDnMax * nYDn + bYDnMin);

		float bestD = 1.e10;
		int bestDn, bestUp;

		int lastUp = 0, lastDn = 0;

		while (true)
		{
			float yUp = y0Up + tracker.HitDataY(rowUp, ihUp) * stepYUp;
			float zUp = z0Up + tracker.HitDataZ(rowUp, ihUp) * stepZUp;

			float yDn = y0Dn + tracker.HitDataY(rowDn, ihDn) * stepYDn;
			float zDn = z0Dn + tracker.HitDataZ(rowDn, ihDn) * stepZDn;

			int ihUpNext, ihDnNext;
			if (ihUp >= ihUpMax)
			{
				if (bZUpMin < bZUpMax)
				{
					ihUpNext = tracker.Data().FirstHitInBin(rowUp, (bZUpMin + 1) * nYUp + bYUpMin);
				}
				else
				{
					lastUp = 1;
				}
			}
			else
			{
				ihUpNext = ihUp + 1;
			}

			if (ihDn <= ihDnMin)
			{
				if (bZDnMax > bZDnMin)
				{
					ihDnNext = tracker.Data().FirstHitInBin(rowDn, (bZDnMax - 1) * nYDn + bYDnMax);
				}
				else
				{
					lastDn = 1;
				}
			}
			else
			{
				ihDnNext = ihDn - 1;
			}

			
			

			float dUp, dDn;
			if (!lastUp)
			{
				const float yUpNext = y0Up + tracker.HitDataY(rowUp, ihUpNext) * stepYUp;
				const float zUpNext = z0Up + tracker.HitDataZ(rowUp, ihUpNext) * stepZUp;
				const float dYUp = s.fUpDx * (yUpNext - y) - s.fDnDx * (yDn - y);
				const float dZUp = s.fUpDx * (zUpNext - y) - s.fDnDx * (zDn - z);
#ifdef NFDEBUG
				printf("Checking Up y: %f nexty: %f z: %f nextz: %f\n", yUp, yUpNext, zUp, zUpNext);
#endif
				dUp = dYUp * dYUp + dZUp * dZUp;
			}

			if (!lastDn)
			{
				const float yDnNext = y0Dn + tracker.HitDataY(rowDn, ihDnNext) * stepYDn;
				const float zDnNext = z0Dn + tracker.HitDataZ(rowDn, ihDnNext) * stepZDn;
				const float dYDn = s.fDnDx * (yDnNext - y) - s.fUpDx * (yUp - y);
				const float dZDn = s.fDnDx * (zDnNext - y) - s.fUpDx * (zUp - z);
#ifdef NFDEBUG
				printf("Checking Dn y: %f nexty: %f z: %f nextz: %f\n", yDn, yDnNext, zDn, zDnNext);
#endif
				dDn = dYDn * dYDn + dZDn * dZDn;
			}

			float d;
			if (lastDn || (dUp < dDn && !lastUp))
			{
				if (lastUp) break;
				d = dUp;
				if (ihUp >= ihUpMax)
				{
					bZUpMin++;
					ihUpMax = tracker.Data().FirstHitInBin(rowUp, bZUpMin * nYUp + bYUpMax + 1) - 1;
				}
				ihUp = ihUpNext;
			}
			else
			{
				d = dDn;
				if (ihDn <= ihDnMin)
				{
					bZDnMax--;
					ihDnMin = tracker.Data().FirstHitInBin(rowDn, bZDnMax * nYDn + bYDnMin);
				}
				ihDn = ihDnNext;
			}
			if (d < bestD)
			{
				bestD = d;
				bestUp = ihUp;
				bestDn = ihDn;
			}
		}

		if (bestD < chi2Cut)
		{
			linkUp = bestUp;
			linkDn = bestDn;
		}

#else
		//Old Slow Neighbours finder
      unsigned short *neighUp = s.fB[iThread];
      float2 *yzUp = s.fA[iThread];
#if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
	  unsigned short neighUp2[kMaxN - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP];
	  float2 yzUp2[kMaxN - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP];
#endif
      //unsigned short neighUp[5];
      //float2 yzUp[5];

		int nNeighUp = 0;
        AliHLTTPCCAHitArea areaDn, areaUp;
        areaUp.Init( rowUp, tracker.Data(), y*s.fUpTx, z*s.fUpTx, kAreaSize, kAreaSize );
        areaDn.Init( rowDn, tracker.Data(), y*s.fDnTx, z*s.fDnTx, kAreaSize, kAreaSize );

        do {
          AliHLTTPCCAHit h;
          int i = areaUp.GetNext( tracker, rowUp, tracker.Data(), &h );
          if ( i < 0 ) break;
#if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
		  if (nNeighUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP)
		  {
			neighUp2[nNeighUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] = ( unsigned short ) i;
			yzUp2[nNeighUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] = CAMath::MakeFloat2( s.fDnDx * ( h.Y() - y ), s.fDnDx * ( h.Z() - z ) );
		  }
		  else
#endif
		  {
			neighUp[nNeighUp] = ( unsigned short ) i;
			yzUp[nNeighUp] = CAMath::MakeFloat2( s.fDnDx * ( h.Y() - y ), s.fDnDx * ( h.Z() - z ) );
		  }
          if ( ++nNeighUp >= kMaxN ) break;
        } while ( 1 );

        int nNeighDn = 0;

        if ( nNeighUp > 0 ) {

          int bestDn = -1, bestUp = -1;
          float bestD = 1.e10;

          do {
            AliHLTTPCCAHit h;
            int i = areaDn.GetNext( tracker, rowDn, tracker.Data(), &h );
            if ( i < 0 ) break;

            nNeighDn++;
            float2 yzdn = CAMath::MakeFloat2( s.fUpDx * ( h.Y() - y ), s.fUpDx * ( h.Z() - z ) );

            for ( int iUp = 0; iUp < nNeighUp; iUp++ ) {
#if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
			  float2 yzup = iUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP ? yzUp2[iUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] : yzUp[iUp];
#else
              float2 yzup = yzUp[iUp];
#endif

              float dy = yzdn.x - yzup.x;
              float dz = yzdn.y - yzup.y;
              float d = dy * dy + dz * dz;
              if ( d < bestD ) {
                bestD = d;
                bestDn = i;
                bestUp = iUp;
              }
            }
          } while ( 1 );

          if ( bestD <= chi2Cut ) {
#if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
			linkUp = bestUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP ? neighUp2[bestUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] : neighUp[bestUp];
#else
            linkUp = neighUp[bestUp];
#endif
            linkDn = bestDn;
          }
        }
#ifdef DRAW
        std::cout << "n NeighUp = " << nNeighUp << ", n NeighDn = " << nNeighDn << std::endl;
#endif //DRAW
#endif //FAST_NEIGHBOURS_FINDER
      }

      tracker.SetHitLinkUpData( row, ih, linkUp );
      tracker.SetHitLinkDownData( row, ih, linkDn );
#ifdef DRAW
      std::cout << "Links for row " << s.fIRow << ", hit " << ih << ": " << linkUp << " " << linkDn << std::endl;
      if ( s.fIRow == 22 && ih == 5 ) {
        AliHLTTPCCADisplay::Instance().DrawSliceLink( s.fIRow, ih, -1, -1, 1 );
        AliHLTTPCCADisplay::Instance().DrawSliceHit( s.fIRow, ih, kBlue, 1. );
        AliHLTTPCCADisplay::Instance().Ask();
      }
#endif //DRAW
    }
  }
}
Commit	Line	Data
ce565086	1	// @(#) $Id: AliHLTTPCCANeighboursFinder1.cxx 27042 2008-07-02 12:06:02Z richterm $
ce565086	2	// **************************************************************************
fbb9b71b	3	// This file is property of and copyright by the ALICE HLT Project *
00d07bcd	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. *
ce565086	17	// *
00d07bcd	18	//***************************************************************************
	19
	20	#include "AliHLTTPCCANeighboursFinder.h"
	21	#include "AliHLTTPCCAMath.h"
	22	#include "AliHLTTPCCAHitArea.h"
	23	#include "AliHLTTPCCATracker.h"
a59a784e	24	#include "AliHLTTPCCAHit.h"
00d07bcd	25
ce565086	26	//#define DRAW
	27
	28	#ifdef DRAW
	29	#include "AliHLTTPCCADisplay.h"
31649d4b	30	#endif //DRAW
ce565086	31
00d07bcd	32	GPUd() void AliHLTTPCCANeighboursFinder::Thread
fbb9b71b	33	( int /nBlocks/, int nThreads, int iBlock, int iThread, int iSync,
00d07bcd	34	AliHLTTPCCASharedMemory &s, AliHLTTPCCATracker &tracker )
	35	{
	36	//* find neighbours
fbb9b71b	37
fbb9b71b	38	if ( iSync == 0 ) {
b22af1bf	39	#ifdef HLTCA_GPUCODE
	40	for (int i = iThread;i < sizeof(AliHLTTPCCARow) / sizeof(int);i += nThreads)
	41	{
	42	reinterpret_cast<int>(&s.fRow)[i] = reinterpret_cast<int>(&tracker.SliceDataRows()[iBlock])[i];
	43	if (iBlock >= 2 && iBlock <= tracker.Param().NRows() - 3)
	44	{
	45	reinterpret_cast<int>(&s.fRowUp)[i] = reinterpret_cast<int>(&tracker.SliceDataRows()[iBlock + 2])[i];
	46	reinterpret_cast<int>(&s.fRowDown)[i] = reinterpret_cast<int>(&tracker.SliceDataRows()[iBlock - 2])[i];
	47	}
	48	}
	49	__syncthreads();
	50	#endif
fbb9b71b	51	if ( iThread == 0 ) {
	52	s.fNRows = tracker.Param().NRows();
	53	s.fIRow = iBlock;
	54	if ( s.fIRow < s.fNRows ) {
b22af1bf	55	#ifdef HLTCA_GPUCODE
	56	const AliHLTTPCCARow &row = s.fRow;
	57	#else
	58	const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
	59	#endif
fbb9b71b	60	s.fNHits = row.NHits();
fbb9b71b	61
fbb9b71b	62	if ( ( s.fIRow >= 2 ) && ( s.fIRow <= s.fNRows - 3 ) ) {
	63	s.fIRowUp = s.fIRow + 2;
	64	s.fIRowDn = s.fIRow - 2;
4acc2401	65
4acc2401	66	// references to the rows above and below
b22af1bf	67
	68	#ifdef HLTCA_GPUCODE
	69	const AliHLTTPCCARow &rowUp = s.fRowUp;
	70	const AliHLTTPCCARow &rowDn = s.fRowDown;
	71	#else
4acc2401	72	const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
4acc2401	73	const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
b22af1bf	74	#endif
4acc2401	75	// the axis perpendicular to the rows
	76	const float xDn = rowDn.X();
	77	const float x = row.X();
	78	const float xUp = rowUp.X();
	79
	80	// number of hits in rows above and below
fbb9b71b	81	s.fUpNHits = tracker.Row( s.fIRowUp ).NHits();
4acc2401	82	s.fDnNHits = tracker.Row( s.fIRowDn ).NHits();
	83
	84	// distance of the rows (absolute and relative)
fbb9b71b	85	s.fUpDx = xUp - x;
	86	s.fDnDx = xDn - x;
	87	s.fUpTx = xUp / x;
	88	s.fDnTx = xDn / x;
4acc2401	89	// UpTx/DnTx is used to move the HitArea such that central events are preferred (i.e. vertices
	90	// coming from y = 0, z = 0).
	91
b22af1bf	92	//s.fGridUp = tracker.Row( s.fIRowUp ).Grid();
b22af1bf	93	//s.fGridDn = tracker.Row( s.fIRowDn ).Grid();
fbb9b71b	94	}
00d07bcd	95	}
fbb9b71b	96	}
	97	} else if ( iSync == 1 ) {
	98	if ( s.fIRow < s.fNRows ) {
b22af1bf	99	if ( ( s.fIRow <= 1 ) \|\| ( s.fIRow >= s.fNRows - 2 ) ) {
	100	#ifdef HLTCA_GPUCODE
	101	const AliHLTTPCCARow &row = s.fRow;
	102	#else
	103	const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
	104	#endif
fbb9b71b	105	for ( int ih = iThread; ih < s.fNHits; ih += nThreads ) {
4acc2401	106	tracker.SetHitLinkUpData( row, ih, -1 );
4acc2401	107	tracker.SetHitLinkDownData( row, ih, -1 );
fbb9b71b	108	}
fbb9b71b	109	} else {
b22af1bf	110	/*#ifdef HLTCA_GPUCODE
	111	const AliHLTTPCCARow &rowUp = s.fRowUp;
	112	const AliHLTTPCCARow &rowDn = s.fRowDown;
	113	#else
	114	const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
	115	const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
	116	#endif
fbb9b71b	117
fbb9b71b	118	for ( unsigned int ih = iThread; ih < s.fGridUp.N() + s.fGridUp.Ny() + 2; ih += nThreads ) {
4acc2401	119	s.fGridContentUp[ih] = tracker.FirstHitInBin( rowUp, ih );
fbb9b71b	120	}
fbb9b71b	121	for ( unsigned int ih = iThread; ih < s.fGridDn.N() + s.fGridDn.Ny() + 2; ih += nThreads ) {
4acc2401	122	s.fGridContentDn[ih] = tracker.FirstHitInBin( rowDn, ih );
b22af1bf	123	}*/
00d07bcd	124	}
00d07bcd	125	}
fbb9b71b	126	} else if ( iSync == 2 ) {
	127	if ( ( s.fIRow <= 1 ) \|\| ( s.fIRow >= s.fNRows - 2 ) ) return;
	128
fbb9b71b	129	float chi2Cut = 3.3.4 * ( s.fUpDx * s.fUpDx + s.fDnDx * s.fDnDx );
f0fb467d	130	const float kAreaSize = tracker.Param().NeighboursSearchArea();
fbb9b71b	131	//float chi2Cut = 3.3.(s.fUpDxs.fUpDx + s.fDnDxs.fDnDx ); //SG
b22af1bf	132	#define kMaxN 20
	133
	134	#ifdef HLTCA_GPUCODE
	135	const AliHLTTPCCARow &row = s.fRow;
	136	const AliHLTTPCCARow &rowUp = s.fRowUp;
	137	const AliHLTTPCCARow &rowDn = s.fRowDown;
	138	#else
	139	const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
	140	const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
	141	const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
	142	#endif
4acc2401	143	const float y0 = row.Grid().YMin();
	144	const float z0 = row.Grid().ZMin();
	145	const float stepY = row.HstepY();
	146	const float stepZ = row.HstepZ();
fbb9b71b	147
b22af1bf	148	for ( int ih = iThread; ih < s.fNHits; ih += nThreads ) {
fbb9b71b	149
fbb9b71b	150	int linkUp = -1;
	151	int linkDn = -1;
	152
	153	if ( s.fDnNHits > 0 && s.fUpNHits > 0 ) {
	154
31649d4b	155
fbb9b71b	156
4acc2401	157	// coordinates of the hit in the current row
b22af1bf	158	#if defined(HLTCA_GPU_TEXTURE_FETCHa)
	159	ushort2 tmpval = tex1Dfetch(gAliTexRefu2, ((char*) tracker.Data().HitData() - tracker.pData()->GPUTextureBase()) / sizeof(ushort2) + row.HitNumberOffset() + ih);
	160	const float y = y0 + tmpval.x * stepY;
	161	const float z = z0 + tmpval.y * stepZ;
	162	#else
4acc2401	163	const float y = y0 + tracker.HitDataY( row, ih ) * stepY;
4acc2401	164	const float z = z0 + tracker.HitDataZ( row, ih ) * stepZ;
b22af1bf	165	#endif
fbb9b71b	166
31649d4b	167	#ifdef FAST_NEIGHBOURS_FINDER
	168
	169	//#define NFDEBUG
	170
	171	#ifdef NFDEBUG
	172	printf("\nSearching Neighbours for: %f %f\n", y, z);
	173	#endif
	174
	175	const float y0Up = rowUp.Grid().YMin();
	176	const float z0Up = rowUp.Grid().ZMin();
	177	const float stepYUp = rowUp.HstepY();
	178	const float stepZUp = rowUp.HstepZ();
	179
	180	const float y0Dn = rowDn.Grid().YMin();
	181	const float z0Dn = rowDn.Grid().ZMin();
	182	const float stepYDn = rowDn.HstepY();
	183	const float stepZDn = rowDn.HstepZ();
	184
	185	const float yMinUp = y*s.fUpTx - kAreaSize;
	186	const float yMaxUp = y*s.fUpTx + kAreaSize;
	187	const float yMinDn = y*s.fDnTx - kAreaSize;
	188	const float yMaxDn = y*s.fDnTx + kAreaSize;
	189	const float zMinUp = z*s.fUpTx - kAreaSize;
	190	const float zMaxUp = z*s.fUpTx + kAreaSize;
	191	const float zMinDn = z*s.fDnTx - kAreaSize;
	192	const float zMaxDn = z*s.fDnTx + kAreaSize;
	193
	194	int bYUpMin, bZUpMin, bYDnMin, bZDnMin;
	195	rowUp.Grid().GetBin(yMinUp, zMinUp, &bYUpMin, &bZUpMin);
	196	rowDn.Grid().GetBin(yMinDn, zMinDn, &bYDnMin, &bZDnMin);
	197
	198	int bYUpMax, bZUpMax, bYDnMax, bZDnMax;
	199	rowUp.Grid().GetBin(yMaxUp, zMaxUp, &bYUpMax, &bZUpMax);
	200	rowDn.Grid().GetBin(yMaxDn, zMaxDn, &bYDnMax, &bZDnMax);
	201
	202	int nYUp = rowUp.Grid().Ny();
	203	int nYDn = rowDn.Grid().Ny();
	204
	205	int ihUp = tracker.Data().FirstHitInBin(rowUp, bZUpMin * nYUp + bYUpMin);
	206	int ihDn = bZDnMax * nYDn + bYDnMax >= rowDn.Grid().N() ? (rowDn.NHits() - 1) : (tracker.Data().FirstHitInBin(rowDn, bZDnMax * nYDn + bYDnMax + 1) - 1);
	207
	208	int ihUpMax = tracker.Data().FirstHitInBin(rowUp, bZUpMin * nYUp + bYUpMax + 1) - 1;
	209	int ihDnMin = tracker.Data().FirstHitInBin(rowDn, bZDnMax * nYDn + bYDnMin);
	210
	211	float bestD = 1.e10;
	212	int bestDn, bestUp;
	213
	214	int lastUp = 0, lastDn = 0;
	215
	216	while (true)
	217	{
	218	float yUp = y0Up + tracker.HitDataY(rowUp, ihUp) * stepYUp;
	219	float zUp = z0Up + tracker.HitDataZ(rowUp, ihUp) * stepZUp;
	220
	221	float yDn = y0Dn + tracker.HitDataY(rowDn, ihDn) * stepYDn;
	222	float zDn = z0Dn + tracker.HitDataZ(rowDn, ihDn) * stepZDn;
	223
	224	int ihUpNext, ihDnNext;
	225	if (ihUp >= ihUpMax)
	226	{
	227	if (bZUpMin < bZUpMax)
	228	{
	229	ihUpNext = tracker.Data().FirstHitInBin(rowUp, (bZUpMin + 1) * nYUp + bYUpMin);
	230	}
231	else
232	{
233	lastUp = 1;
234	}
235	}
236	else
237	{
238	ihUpNext = ihUp + 1;
239	}
240
241	if (ihDn <= ihDnMin)
242	{
243	if (bZDnMax > bZDnMin)
244	{
245	ihDnNext = tracker.Data().FirstHitInBin(rowDn, (bZDnMax - 1) * nYDn + bYDnMax);
246	}
247	else
248	{
249	lastDn = 1;
250	}
251	}
252	else
253	{
254	ihDnNext = ihDn - 1;
255	}
256
257
258
259
260	float dUp, dDn;
261	if (!lastUp)
262	{
263	const float yUpNext = y0Up + tracker.HitDataY(rowUp, ihUpNext) * stepYUp;
264	const float zUpNext = z0Up + tracker.HitDataZ(rowUp, ihUpNext) * stepZUp;
265	const float dYUp = s.fUpDx * (yUpNext - y) - s.fDnDx * (yDn - y);
266	const float dZUp = s.fUpDx * (zUpNext - y) - s.fDnDx * (zDn - z);
267	#ifdef NFDEBUG
268	printf("Checking Up y: %f nexty: %f z: %f nextz: %f\n", yUp, yUpNext, zUp, zUpNext);
269	#endif
270	dUp = dYUp * dYUp + dZUp * dZUp;
271	}
272
273	if (!lastDn)
274	{
275	const float yDnNext = y0Dn + tracker.HitDataY(rowDn, ihDnNext) * stepYDn;
276	const float zDnNext = z0Dn + tracker.HitDataZ(rowDn, ihDnNext) * stepZDn;
277	const float dYDn = s.fDnDx * (yDnNext - y) - s.fUpDx * (yUp - y);
278	const float dZDn = s.fDnDx * (zDnNext - y) - s.fUpDx * (zUp - z);
279	#ifdef NFDEBUG
280	printf("Checking Dn y: %f nexty: %f z: %f nextz: %f\n", yDn, yDnNext, zDn, zDnNext);
281	#endif
282	dDn = dYDn * dYDn + dZDn * dZDn;
283	}
284
285	float d;
286	if (lastDn \|\| (dUp < dDn && !lastUp))
287	{
288	if (lastUp) break;
289	d = dUp;
290	if (ihUp >= ihUpMax)
291	{
292	bZUpMin++;
293	ihUpMax = tracker.Data().FirstHitInBin(rowUp, bZUpMin * nYUp + bYUpMax + 1) - 1;
294	}
295	ihUp = ihUpNext;
296	}
297	else
298	{
299	d = dDn;
300	if (ihDn <= ihDnMin)
301	{
302	bZDnMax--;
303	ihDnMin = tracker.Data().FirstHitInBin(rowDn, bZDnMax * nYDn + bYDnMin);
304	}
305	ihDn = ihDnNext;
306	}
307	if (d < bestD)
308	{
309	bestD = d;
310	bestUp = ihUp;
311	bestDn = ihDn;
312	}
313	}
314
315	if (bestD < chi2Cut)
316	{
317	linkUp = bestUp;
318	linkDn = bestDn;
319	}
320
321	#else
322	//Old Slow Neighbours finder
323	unsigned short *neighUp = s.fB[iThread];
324	float2 *yzUp = s.fA[iThread];
325	#if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
326	unsigned short neighUp2[kMaxN - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP];
327	float2 yzUp2[kMaxN - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP];
328	#endif
329	//unsigned short neighUp[5];
330	//float2 yzUp[5];
331
332	int nNeighUp = 0;
fbb9b71b	333	AliHLTTPCCAHitArea areaDn, areaUp;
4acc2401	334	areaUp.Init( rowUp, tracker.Data(), ys.fUpTx, zs.fUpTx, kAreaSize, kAreaSize );
4acc2401	335	areaDn.Init( rowDn, tracker.Data(), ys.fDnTx, zs.fDnTx, kAreaSize, kAreaSize );
fbb9b71b	336
	337	do {
	338	AliHLTTPCCAHit h;
4acc2401	339	int i = areaUp.GetNext( tracker, rowUp, tracker.Data(), &h );
fbb9b71b	340	if ( i < 0 ) break;
b22af1bf	341	#if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
	342	if (nNeighUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP)
	343	{
	344	neighUp2[nNeighUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] = ( unsigned short ) i;
	345	yzUp2[nNeighUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] = CAMath::MakeFloat2( s.fDnDx * ( h.Y() - y ), s.fDnDx * ( h.Z() - z ) );
	346	}
	347	else
	348	#endif
	349	{
	350	neighUp[nNeighUp] = ( unsigned short ) i;
	351	yzUp[nNeighUp] = CAMath::MakeFloat2( s.fDnDx * ( h.Y() - y ), s.fDnDx * ( h.Z() - z ) );
	352	}
fbb9b71b	353	if ( ++nNeighUp >= kMaxN ) break;
	354	} while ( 1 );
	355
	356	int nNeighDn = 0;
	357
	358	if ( nNeighUp > 0 ) {
	359
	360	int bestDn = -1, bestUp = -1;
	361	float bestD = 1.e10;
	362
	363	do {
	364	AliHLTTPCCAHit h;
4acc2401	365	int i = areaDn.GetNext( tracker, rowDn, tracker.Data(), &h );
fbb9b71b	366	if ( i < 0 ) break;
	367
	368	nNeighDn++;
	369	float2 yzdn = CAMath::MakeFloat2( s.fUpDx * ( h.Y() - y ), s.fUpDx * ( h.Z() - z ) );
	370
	371	for ( int iUp = 0; iUp < nNeighUp; iUp++ ) {
b22af1bf	372	#if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
	373	float2 yzup = iUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP ? yzUp2[iUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] : yzUp[iUp];
	374	#else
fbb9b71b	375	float2 yzup = yzUp[iUp];
b22af1bf	376	#endif
b22af1bf	377
fbb9b71b	378	float dy = yzdn.x - yzup.x;
	379	float dz = yzdn.y - yzup.y;
	380	float d = dy * dy + dz * dz;
	381	if ( d < bestD ) {
	382	bestD = d;
	383	bestDn = i;
	384	bestUp = iUp;
	385	}
	386	}
	387	} while ( 1 );
	388
	389	if ( bestD <= chi2Cut ) {
b22af1bf	390	#if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
	391	linkUp = bestUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP ? neighUp2[bestUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] : neighUp[bestUp];
	392	#else
fbb9b71b	393	linkUp = neighUp[bestUp];
b22af1bf	394	#endif
fbb9b71b	395	linkDn = bestDn;
	396	}
	397	}
ce565086	398	#ifdef DRAW
fbb9b71b	399	std::cout << "n NeighUp = " << nNeighUp << ", n NeighDn = " << nNeighDn << std::endl;
31649d4b	400	#endif //DRAW
31649d4b	401	#endif //FAST_NEIGHBOURS_FINDER
fbb9b71b	402	}
fbb9b71b	403
4acc2401	404	tracker.SetHitLinkUpData( row, ih, linkUp );
4acc2401	405	tracker.SetHitLinkDownData( row, ih, linkDn );
ce565086	406	#ifdef DRAW
fbb9b71b	407	std::cout << "Links for row " << s.fIRow << ", hit " << ih << ": " << linkUp << " " << linkDn << std::endl;
	408	if ( s.fIRow == 22 && ih == 5 ) {
	409	AliHLTTPCCADisplay::Instance().DrawSliceLink( s.fIRow, ih, -1, -1, 1 );
	410	AliHLTTPCCADisplay::Instance().DrawSliceHit( s.fIRow, ih, kBlue, 1. );
	411	AliHLTTPCCADisplay::Instance().Ask();
00d07bcd	412	}
31649d4b	413	#endif //DRAW
00d07bcd	414	}
fbb9b71b	415	}
00d07bcd	416	}
00d07bcd	417