[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCTrackGeometry.cxx

// $Id$

//**************************************************************************
//* This file is property of and copyright by the ALICE HLT Project        * 
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
//*                  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.                  *
//**************************************************************************

/// @file   AliHLTTPCTrackGeometry.cxx
/// @author Matthias Richter
/// @date   2011-05-20
/// @brief  Desciption of a track by a sequence of track points
///

#include "AliHLTTPCTrackGeometry.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCSpacePointData.h"
#include "AliHLTTPCClusterDataFormat.h"
#include "AliHLTTPCSpacePointContainer.h"
#include "AliHLTTPCDefinitions.h"
#include "AliHLTComponent.h"
#include "AliHLTGlobalBarrelTrack.h"
#include "TMath.h"
#include "TH2F.h"
#include <memory>

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

AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry()
  : AliHLTTrackGeometry()
  , fRawTrackPoints()
{
  /// standard constructor
}

AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry(const AliHLTTPCTrackGeometry& src)
  : AliHLTTrackGeometry(src)
  , fRawTrackPoints(src.fRawTrackPoints)
{
  /// copy constructor
}

AliHLTTPCTrackGeometry& AliHLTTPCTrackGeometry::operator=(const AliHLTTPCTrackGeometry& src)
{
  /// assignment operator
  AliHLTTrackGeometry::operator=(src);
  fRawTrackPoints.assign(src.fRawTrackPoints.begin(), src.fRawTrackPoints.end());
  return *this;
}

AliHLTTPCTrackGeometry::~AliHLTTPCTrackGeometry()
{
  /// destructor
}

float AliHLTTPCTrackGeometry::GetPlaneAlpha(AliHLTUInt32_t planeId) const
{
  /// alpha of the plane
  UInt_t slice=AliHLTTPCSpacePointData::GetSlice(planeId);
  float alpha=( slice + 0.5 ) * TMath::Pi() / 9.0;
  if (alpha>TMath::TwoPi()) alpha-=TMath::TwoPi();
  return alpha;
}

float AliHLTTPCTrackGeometry::GetPlaneR(AliHLTUInt32_t planeId) const
{
  /// radial distance from global {0,0,0}
  UInt_t partition=AliHLTTPCSpacePointData::GetPatch(planeId);
  UInt_t number=AliHLTTPCSpacePointData::GetNumber(planeId);
  Int_t row=AliHLTTPCTransform::GetFirstRow(partition)+number;
  return AliHLTTPCTransform::Row2X(row);
}

float AliHLTTPCTrackGeometry::GetPlaneTheta(AliHLTUInt32_t /*planeId*/) const
{
  /// theta of the plane
  return 0.0;
}

bool AliHLTTPCTrackGeometry::CheckBounds(AliHLTUInt32_t planeId, float u, float /*v*/) const
{
  /// check bounds in u and v coordinate
  float r=GetPlaneR(planeId);
  if (r<AliHLTTPCTransform::GetFirstRow(0)) return false;

  // TODO: check if the pad width needs to be considered here
  return TMath::Abs(TMath::ASin(u/r))<=TMath::Pi()/18;
}

int AliHLTTPCTrackGeometry::CalculateTrackPoints(const AliHLTExternalTrackParam& track)
{
  /// calculate the track points, expects the global magnetic field to be initialized
  AliHLTGlobalBarrelTrack bt(track);
  return CalculateTrackPoints(bt);
}

int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track)
{
  /// calculate the track points, expects the global magnetic field to be initialized
  int iResult=0;
  int firstpadrow=0;
  for (;
       firstpadrow<AliHLTTPCTransform::GetNRows() && 
	 AliHLTTPCTransform::Row2X(firstpadrow)+AliHLTTPCTransform::GetPadLength(firstpadrow)<track.GetX();
       firstpadrow++);
  if (firstpadrow>=AliHLTTPCTransform::GetNRows()) return 0;
  iResult=CalculateTrackPoints(track, firstpadrow, 1);
  if (iResult>=0 && firstpadrow>0)
    iResult=CalculateTrackPoints(track, firstpadrow-1, -1);
  return iResult;
}

int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track, int firstpadrow, int step)
{
  /// calculate the track points, expects the global magnetic field to be initialized
  float offsetAlpha=0.0;
  for (int padrow=firstpadrow; padrow>=0 && padrow<AliHLTTPCTransform::GetNRows(); padrow+=step) {
    float x=AliHLTTPCTransform::Row2X(padrow);
    float y=0.0;
    float z=0.0;

    int maxshift=9;
    int shift=0;
    int result=0;
    do {
      // start calculation of crossing points with padrow planes in the slice of the first point
      // plane alpha corresponds to alpha of the track, switch to neighboring slice if the result
      // is out of bounds
      if ((result=track.CalculateCrossingPoint(x, track.GetAlpha()-offsetAlpha, y, z))<1) break;
      float pointAlpha=TMath::ATan(y/x);
      if (TMath::Abs(pointAlpha)>TMath::Pi()/18) {
	offsetAlpha+=(pointAlpha>0?-1:1)*TMath::Pi()/9;
      	result=0;
      }
    } while (result==0 && shift++<maxshift);
    if (result<1) continue;
    float planealpha=track.GetAlpha()-offsetAlpha;
    if (planealpha<0) planealpha+=TMath::TwoPi();
    int slice=int(9*planealpha/TMath::Pi());
    //if (z<0) slice+=18;
    int partition=AliHLTTPCTransform::GetPatch(padrow);
    int row=padrow-AliHLTTPCTransform::GetFirstRow(partition);
    UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
    if (TMath::Abs(planealpha-GetPlaneAlpha(id))>0.0001) {
      HLTError("alpha missmatch for plane %08x (slice %d): alpha from id %f (%.0f), expected %f (%.0f)", id, slice, GetPlaneAlpha(id), 180*GetPlaneAlpha(id)/TMath::Pi(), planealpha, 180*planealpha/TMath::Pi());
    }
    if (AddTrackPoint(AliHLTTrackPoint(id, y, z), AliHLTTPCSpacePointData::GetID(slice, partition, 0))>=0) {
      Float_t rpt[3]={0.,y,z}; // row pad time
      AliHLTTPCTransform::LocHLT2Raw(rpt, slice, padrow);
      // FIXME: there is a mismatch in the definition of the pad coordinate
      // should be with respect to middle of pad, that's why the offset of
      // 0.5 has been applied when calling the AliHLTTPCClusterTransformation
      // and for conversion to AliTPCclusterMI in AliHLTTPCClusterAccessHLTOUT
      // AliHLTTPCTransform::LocHLT2Raw seems to define this shift in the
      // opposite direction
      rpt[1]+=1.;
      fRawTrackPoints.push_back(AliHLTTrackPoint(id, rpt[1], rpt[2]));
    }
  }
  return 0;
}

int AliHLTTPCTrackGeometry::FindMatchingTrackPoint(AliHLTUInt32_t spacepointId, float spacepoint[3], AliHLTUInt32_t& planeId)
{
  /// find the track point which can be associated to a spacepoint with coordinates and id
  UInt_t slice=AliHLTTPCSpacePointData::GetSlice(spacepointId);
  UInt_t partition=AliHLTTPCSpacePointData::GetPatch(spacepointId);
  int row=AliHLTTPCTransform::GetPadRow(spacepoint[0]);
  bool bSpecialRow=row==30 || row==90 || row==139;
  if (row<AliHLTTPCTransform::GetFirstRow(partition) || row>AliHLTTPCTransform::GetLastRow(partition)) {
    HLTError("row number %d calculated from x value %f is outside slice %d partition %d", row, spacepoint[0], slice, partition);
    return -EINVAL;
  }

  // find the crossing point of the track with the padrow plane where
  // the spacepoint is
  // 1) calculate plane id from slice, partition and row (within partition)
  row-=AliHLTTPCTransform::GetFirstRow(partition);
  UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
  const AliHLTTrackPoint* point=GetTrackPoint(id);
  // track might be outside the partition and cross the central membrane
  // search in the other half of the TPC
  if (!point && slice<18) {
    // search in the neighboring partition on the C side
    id=AliHLTTPCSpacePointData::GetID(slice+18, partition, row);
    point=GetTrackPoint(id);
  } else if (!point && slice>=18) {
    // search in the neighboring partition on the A side
    id=AliHLTTPCSpacePointData::GetID(slice-18, partition, row);
    point=GetTrackPoint(id);
  }
  
  // search in the neighboring partition, this takes account for rows
  // 30, 90, and 139 which are partly in one and the other partition
  if (!point && bSpecialRow) {
    row+=AliHLTTPCTransform::GetFirstRow(partition);
    row-=AliHLTTPCTransform::GetFirstRow(partition-1);
    id=AliHLTTPCSpacePointData::GetID(slice, partition-1, row);
    point=GetTrackPoint(id);
    if (!point && slice<18) {
      // search in the neighboring partition on the C side
      id=AliHLTTPCSpacePointData::GetID(slice+18, partition-1, row);
      point=GetTrackPoint(id);
    } else if (!point && slice>=18) {
      // search in the neighboring partition on the A side
      id=AliHLTTPCSpacePointData::GetID(slice-18, partition-1, row);
      point=GetTrackPoint(id);
    }
  }

  if (point) {
    planeId=id;
    if (point->HaveAssociatedSpacePoint()) {
      if (GetVerbosity()>2) HLTInfo("descarding spacepoint 0x%08x z=%f y=%f z=%f: track point 0x%08x already occupied", spacepoint[0], spacepoint[1], spacepoint[2], planeId);
      return 0; // already occupied
    }
    float maxdy=2.;
    float maxdz=2.;
    if (TMath::Abs(point->GetU()-spacepoint[1])>maxdy) {
      if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x y %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetU(), maxdy);
      return -ENOENT;
    }
    if (TMath::Abs(point->GetV()-spacepoint[2])>maxdz) {
      if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x z %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetV(), maxdz);
      return -ENOENT;
    }
    return 1;
  }
  return -ENOENT;
}


int AliHLTTPCTrackGeometry::RegisterTrackPoints(AliHLTTrackGrid* pGrid) const
{
  /// register track points in the index grid, at this step the number
  /// of tracks in each cell is counted
  if (!pGrid) return -EINVAL;
  int iResult=0;
  for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
       tp!=TrackPoints().end() && iResult>=0; tp++) {
    AliHLTUInt32_t id=tp->GetId();
    iResult=pGrid->CountSpacePoint(AliHLTTPCSpacePointData::GetSlice(id),
				   AliHLTTPCSpacePointData::GetPatch(id),
				   AliHLTTPCSpacePointData::GetNumber(id));
  }
  return iResult;
}

int AliHLTTPCTrackGeometry::FillTrackPoints(AliHLTTrackGrid* pGrid) const
{
  /// fill track points to index grid
  if (!pGrid) return -EINVAL;
  int iResult=0;
  for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
       tp!=TrackPoints().end() && iResult>=0; tp++) {
    AliHLTUInt32_t id=tp->GetId();
    iResult=pGrid->AddSpacePoint(GetTrackId(),
				 AliHLTTPCSpacePointData::GetSlice(id),
				 AliHLTTPCSpacePointData::GetPatch(id),
				 AliHLTTPCSpacePointData::GetNumber(id));
  }
  return iResult;
}

AliHLTSpacePointContainer* AliHLTTPCTrackGeometry::ConvertToSpacePoints(bool bAssociated) const
{
  /// create a collection of all points
  std::auto_ptr<AliHLTTPCSpacePointContainer> spacepoints(new AliHLTTPCSpacePointContainer);
  if (!spacepoints.get()) return NULL;

  const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
  unsigned i=0;
  while (i<trackPoints.size()) {
    // allocate buffer for all points, even though the buffer might not be filled
    // completely because of a partition change
    int nofPoints=trackPoints.size()-i;
    int blocksize=sizeof(AliHLTTPCClusterData)+nofPoints*sizeof(AliHLTTPCSpacePointData);
    AliHLTUInt8_t* pBuffer=spacepoints->Alloc(blocksize);
    if (!pBuffer) return NULL;
    AliHLTTPCClusterData* pClusterData=reinterpret_cast<AliHLTTPCClusterData*>(pBuffer);
    pClusterData->fSpacePointCnt=0;
    AliHLTTPCSpacePointData* pClusters=pClusterData->fSpacePoints;
    int currentSlice=-1;
    int currentPartition=-1;
    for (; i<trackPoints.size(); i++) {
      if (bAssociated && !trackPoints[i].HaveAssociatedSpacePoint()) continue;
      AliHLTUInt32_t planeId=trackPoints[i].GetId();
      int slice=AliHLTTPCSpacePointData::GetSlice(planeId);
      int partition=AliHLTTPCSpacePointData::GetPatch(planeId);
      int number=AliHLTTPCSpacePointData::GetNumber(planeId);
      if ((currentSlice>=0 && currentSlice!=slice) || (currentPartition>=0 && currentPartition!=partition)) {
	// change of partition or slice, need to go to next block
	// 2011-07-26 currently all spacepoints go into one block, if separated
	// blocks per partition are needed one has to leave the inner loop here
	// and set the data block specification below
	// Caution: not tested, only the last block seems to make it through
	//break;
      }
      currentSlice=slice;
      currentPartition=partition;
      pClusters[pClusterData->fSpacePointCnt].fX=GetPlaneR(planeId);
      pClusters[pClusterData->fSpacePointCnt].fY=trackPoints[i].GetU();
      pClusters[pClusterData->fSpacePointCnt].fZ=trackPoints[i].GetV();
      pClusters[pClusterData->fSpacePointCnt].fID=planeId;
      pClusters[pClusterData->fSpacePointCnt].fPadRow=AliHLTTPCTransform::GetFirstRow(partition)+number;
      pClusters[pClusterData->fSpacePointCnt].fSigmaY2=0.;
      pClusters[pClusterData->fSpacePointCnt].fSigmaZ2=0.;
      pClusters[pClusterData->fSpacePointCnt].fCharge=0;
      pClusters[pClusterData->fSpacePointCnt].fQMax=0;
      pClusters[pClusterData->fSpacePointCnt].fUsed=0;
      pClusters[pClusterData->fSpacePointCnt].fTrackN=0;
      pClusterData->fSpacePointCnt++;
    }
    AliHLTComponentBlockData bd;
    AliHLTComponent::FillBlockData(bd);
    bd.fPtr=pBuffer;
    bd.fSize=sizeof(AliHLTTPCClusterData)+pClusterData->fSpacePointCnt*sizeof(AliHLTTPCSpacePointData);
    AliHLTComponent::SetDataType(bd.fDataType, "CLUSTERS", "TPC ");
    bd.fSpecification=kAliHLTVoidDataSpec;//AliHLTTPCDefinitions::EncodeDataSpecification(currentSlice, currentSlice, currentPartition, currentPartition);
    spacepoints->AddInputBlock(&bd);
  }

  return spacepoints.release();
}

const AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id) const
{
  /// get raw track point of id
  const AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
  if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
  return p;
}

AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id)
{
  /// get raw track point of id
  AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
  if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
  return p;
}

int AliHLTTPCTrackGeometry::FillRawResidual(int coordinate, TH2* histo, AliHLTSpacePointContainer* points) const
{
  // fill residual histogram
  if (!histo || !points) return -EINVAL;
  const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
  for (vector<AliHLTTrackPoint>::const_iterator trackpoint=trackPoints.begin();
       trackpoint!=trackPoints.end(); trackpoint++) {
    if (!trackpoint->HaveAssociatedSpacePoint()) continue;
    for (vector<AliHLTTrackSpacepoint>::const_iterator sp=(trackpoint->GetSpacepoints()).begin();
	 sp!=(trackpoint->GetSpacepoints()).end(); sp++) {
    AliHLTUInt32_t spacepointId=sp->fId;
    vector<AliHLTTrackPoint>::const_iterator rawpoint=find(fRawTrackPoints.begin(), fRawTrackPoints.end(), trackpoint->GetId());
    if (rawpoint==fRawTrackPoints.end()) {
      HLTError("can not find track raw coordinates of track point 0x%08x", trackpoint->GetId());
      continue;
    }
    if (!points->Check(spacepointId)) {
      //HLTError("can not find associated space point 0x%08x of track point 0x%08x", spacepointId, trackpoint->GetId());
      continue;
    }
    float value=0.;
    if (coordinate==0) {
      value=rawpoint->GetU()-points->GetY(spacepointId);
      histo->Fill(GetPlaneR(trackpoint->GetId()), value);
    } else {
      value=rawpoint->GetV()-points->GetZ(spacepointId);
      //histo->Fill(GetPlaneR(trackpoint->GetId()), value);
      histo->Fill(rawpoint->GetV(), value);
    }
    }
  }
  return 0;
}
Commit	Line	Data
c7585a2a	1	// $Id$
	2
	3	//**************************************************************************
	4	//* This file is property of and copyright by the ALICE HLT Project *
	5	//* ALICE Experiment at CERN, All rights reserved. *
	6	//* *
	7	//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
	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	/// @file AliHLTTPCTrackGeometry.cxx
	20	/// @author Matthias Richter
	21	/// @date 2011-05-20
	22	/// @brief Desciption of a track by a sequence of track points
	23	///
	24
	25	#include "AliHLTTPCTrackGeometry.h"
	26	#include "AliHLTTPCTransform.h"
	27	#include "AliHLTTPCSpacePointData.h"
eb50ad6b	28	#include "AliHLTTPCClusterDataFormat.h"
	29	#include "AliHLTTPCSpacePointContainer.h"
	30	#include "AliHLTTPCDefinitions.h"
	31	#include "AliHLTComponent.h"
c7585a2a	32	#include "AliHLTGlobalBarrelTrack.h"
c7585a2a	33	#include "TMath.h"
bbe2c314	34	#include "TH2F.h"
eb50ad6b	35	#include <memory>
c7585a2a	36
	37	/** ROOT macro for the implementation of ROOT specific class methods */
	38	ClassImp(AliHLTTPCTrackGeometry)
	39
	40	AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry()
	41	: AliHLTTrackGeometry()
bbe2c314	42	, fRawTrackPoints()
c7585a2a	43	{
	44	/// standard constructor
	45	}
	46
	47	AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry(const AliHLTTPCTrackGeometry& src)
	48	: AliHLTTrackGeometry(src)
bbe2c314	49	, fRawTrackPoints(src.fRawTrackPoints)
c7585a2a	50	{
	51	/// copy constructor
	52	}
	53
	54	AliHLTTPCTrackGeometry& AliHLTTPCTrackGeometry::operator=(const AliHLTTPCTrackGeometry& src)
	55	{
	56	/// assignment operator
	57	AliHLTTrackGeometry::operator=(src);
bbe2c314	58	fRawTrackPoints.assign(src.fRawTrackPoints.begin(), src.fRawTrackPoints.end());
c7585a2a	59	return *this;
	60	}
	61
	62	AliHLTTPCTrackGeometry::~AliHLTTPCTrackGeometry()
	63	{
	64	/// destructor
	65	}
	66
	67	float AliHLTTPCTrackGeometry::GetPlaneAlpha(AliHLTUInt32_t planeId) const
	68	{
	69	/// alpha of the plane
	70	UInt_t slice=AliHLTTPCSpacePointData::GetSlice(planeId);
	71	float alpha=( slice + 0.5 ) * TMath::Pi() / 9.0;
	72	if (alpha>TMath::TwoPi()) alpha-=TMath::TwoPi();
	73	return alpha;
	74	}
	75
	76	float AliHLTTPCTrackGeometry::GetPlaneR(AliHLTUInt32_t planeId) const
	77	{
	78	/// radial distance from global {0,0,0}
	79	UInt_t partition=AliHLTTPCSpacePointData::GetPatch(planeId);
	80	UInt_t number=AliHLTTPCSpacePointData::GetNumber(planeId);
	81	Int_t row=AliHLTTPCTransform::GetFirstRow(partition)+number;
	82	return AliHLTTPCTransform::Row2X(row);
	83	}
	84
	85	float AliHLTTPCTrackGeometry::GetPlaneTheta(AliHLTUInt32_t /planeId/) const
	86	{
	87	/// theta of the plane
	88	return 0.0;
	89	}
	90
	91	bool AliHLTTPCTrackGeometry::CheckBounds(AliHLTUInt32_t planeId, float u, float /v/) const
	92	{
	93	/// check bounds in u and v coordinate
	94	float r=GetPlaneR(planeId);
	95	if (r<AliHLTTPCTransform::GetFirstRow(0)) return false;
	96
	97	// TODO: check if the pad width needs to be considered here
	98	return TMath::Abs(TMath::ASin(u/r))<=TMath::Pi()/18;
	99	}
	100
	101	int AliHLTTPCTrackGeometry::CalculateTrackPoints(const AliHLTExternalTrackParam& track)
	102	{
	103	/// calculate the track points, expects the global magnetic field to be initialized
	104	AliHLTGlobalBarrelTrack bt(track);
	105	return CalculateTrackPoints(bt);
	106	}
	107
	108	int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track)
	109	{
	110	/// calculate the track points, expects the global magnetic field to be initialized
	111	int iResult=0;
	112	int firstpadrow=0;
eb50ad6b	113	for (;
	114	firstpadrow<AliHLTTPCTransform::GetNRows() &&
	115	AliHLTTPCTransform::Row2X(firstpadrow)+AliHLTTPCTransform::GetPadLength(firstpadrow)<track.GetX();
	116	firstpadrow++);
	117	if (firstpadrow>=AliHLTTPCTransform::GetNRows()) return 0;
c7585a2a	118	iResult=CalculateTrackPoints(track, firstpadrow, 1);
	119	if (iResult>=0 && firstpadrow>0)
	120	iResult=CalculateTrackPoints(track, firstpadrow-1, -1);
	121	return iResult;
	122	}
	123
	124	int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track, int firstpadrow, int step)
	125	{
	126	/// calculate the track points, expects the global magnetic field to be initialized
	127	float offsetAlpha=0.0;
	128	for (int padrow=firstpadrow; padrow>=0 && padrow<AliHLTTPCTransform::GetNRows(); padrow+=step) {
	129	float x=AliHLTTPCTransform::Row2X(padrow);
	130	float y=0.0;
	131	float z=0.0;
	132
	133	int maxshift=9;
	134	int shift=0;
	135	int result=0;
	136	do {
	137	// start calculation of crossing points with padrow planes in the slice of the first point
	138	// plane alpha corresponds to alpha of the track, switch to neighboring slice if the result
	139	// is out of bounds
	140	if ((result=track.CalculateCrossingPoint(x, track.GetAlpha()-offsetAlpha, y, z))<1) break;
	141	float pointAlpha=TMath::ATan(y/x);
	142	if (TMath::Abs(pointAlpha)>TMath::Pi()/18) {
	143	offsetAlpha+=(pointAlpha>0?-1:1)*TMath::Pi()/9;
	144	result=0;
	145	}
	146	} while (result==0 && shift++<maxshift);
	147	if (result<1) continue;
	148	float planealpha=track.GetAlpha()-offsetAlpha;
	149	if (planealpha<0) planealpha+=TMath::TwoPi();
102d9fef	150	int slice=int(9*planealpha/TMath::Pi());
c7585a2a	151	//if (z<0) slice+=18;
	152	int partition=AliHLTTPCTransform::GetPatch(padrow);
	153	int row=padrow-AliHLTTPCTransform::GetFirstRow(partition);
	154	UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
	155	if (TMath::Abs(planealpha-GetPlaneAlpha(id))>0.0001) {
	156	HLTError("alpha missmatch for plane %08x (slice %d): alpha from id %f (%.0f), expected %f (%.0f)", id, slice, GetPlaneAlpha(id), 180GetPlaneAlpha(id)/TMath::Pi(), planealpha, 180planealpha/TMath::Pi());
	157	}
bbe2c314	158	if (AddTrackPoint(AliHLTTrackPoint(id, y, z), AliHLTTPCSpacePointData::GetID(slice, partition, 0))>=0) {
	159	Float_t rpt[3]={0.,y,z}; // row pad time
	160	AliHLTTPCTransform::LocHLT2Raw(rpt, slice, padrow);
	161	// FIXME: there is a mismatch in the definition of the pad coordinate
	162	// should be with respect to middle of pad, that's why the offset of
	163	// 0.5 has been applied when calling the AliHLTTPCClusterTransformation
	164	// and for conversion to AliTPCclusterMI in AliHLTTPCClusterAccessHLTOUT
	165	// AliHLTTPCTransform::LocHLT2Raw seems to define this shift in the
	166	// opposite direction
	167	rpt[1]+=1.;
	168	fRawTrackPoints.push_back(AliHLTTrackPoint(id, rpt[1], rpt[2]));
	169	}
c7585a2a	170	}
	171	return 0;
	172	}
	173
	174	int AliHLTTPCTrackGeometry::FindMatchingTrackPoint(AliHLTUInt32_t spacepointId, float spacepoint[3], AliHLTUInt32_t& planeId)
	175	{
	176	/// find the track point which can be associated to a spacepoint with coordinates and id
	177	UInt_t slice=AliHLTTPCSpacePointData::GetSlice(spacepointId);
	178	UInt_t partition=AliHLTTPCSpacePointData::GetPatch(spacepointId);
	179	int row=AliHLTTPCTransform::GetPadRow(spacepoint[0]);
9003c201	180	bool bSpecialRow=row==30 \|\| row==90 \|\| row==139;
c7585a2a	181	if (row<AliHLTTPCTransform::GetFirstRow(partition) \|\| row>AliHLTTPCTransform::GetLastRow(partition)) {
	182	HLTError("row number %d calculated from x value %f is outside slice %d partition %d", row, spacepoint[0], slice, partition);
	183	return -EINVAL;
	184	}
9003c201	185
	186	// find the crossing point of the track with the padrow plane where
	187	// the spacepoint is
	188	// 1) calculate plane id from slice, partition and row (within partition)
c7585a2a	189	row-=AliHLTTPCTransform::GetFirstRow(partition);
	190	UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
	191	const AliHLTTrackPoint* point=GetTrackPoint(id);
9003c201	192	// track might be outside the partition and cross the central membrane
9003c201	193	// search in the other half of the TPC
c7585a2a	194	if (!point && slice<18) {
9003c201	195	// search in the neighboring partition on the C side
c7585a2a	196	id=AliHLTTPCSpacePointData::GetID(slice+18, partition, row);
	197	point=GetTrackPoint(id);
	198	} else if (!point && slice>=18) {
9003c201	199	// search in the neighboring partition on the A side
c7585a2a	200	id=AliHLTTPCSpacePointData::GetID(slice-18, partition, row);
	201	point=GetTrackPoint(id);
	202	}
9003c201	203
	204	// search in the neighboring partition, this takes account for rows
	205	// 30, 90, and 139 which are partly in one and the other partition
	206	if (!point && bSpecialRow) {
	207	row+=AliHLTTPCTransform::GetFirstRow(partition);
	208	row-=AliHLTTPCTransform::GetFirstRow(partition-1);
	209	id=AliHLTTPCSpacePointData::GetID(slice, partition-1, row);
	210	point=GetTrackPoint(id);
	211	if (!point && slice<18) {
	212	// search in the neighboring partition on the C side
	213	id=AliHLTTPCSpacePointData::GetID(slice+18, partition-1, row);
	214	point=GetTrackPoint(id);
	215	} else if (!point && slice>=18) {
	216	// search in the neighboring partition on the A side
	217	id=AliHLTTPCSpacePointData::GetID(slice-18, partition-1, row);
	218	point=GetTrackPoint(id);
	219	}
	220	}
	221
c7585a2a	222	if (point) {
c7585a2a	223	planeId=id;
9003c201	224	if (point->HaveAssociatedSpacePoint()) {
	225	if (GetVerbosity()>2) HLTInfo("descarding spacepoint 0x%08x z=%f y=%f z=%f: track point 0x%08x already occupied", spacepoint[0], spacepoint[1], spacepoint[2], planeId);
	226	return 0; // already occupied
	227	}
	228	float maxdy=2.;
	229	float maxdz=2.;
	230	if (TMath::Abs(point->GetU()-spacepoint[1])>maxdy) {
	231	if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x y %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetU(), maxdy);
	232	return -ENOENT;
	233	}
	234	if (TMath::Abs(point->GetV()-spacepoint[2])>maxdz) {
	235	if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x z %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetV(), maxdz);
	236	return -ENOENT;
	237	}
c7585a2a	238	return 1;
	239	}
	240	return -ENOENT;
	241	}
eb50ad6b	242
c2bea8b3	243
	244	int AliHLTTPCTrackGeometry::RegisterTrackPoints(AliHLTTrackGrid* pGrid) const
	245	{
	246	/// register track points in the index grid, at this step the number
	247	/// of tracks in each cell is counted
	248	if (!pGrid) return -EINVAL;
	249	int iResult=0;
	250	for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
	251	tp!=TrackPoints().end() && iResult>=0; tp++) {
	252	AliHLTUInt32_t id=tp->GetId();
	253	iResult=pGrid->CountSpacePoint(AliHLTTPCSpacePointData::GetSlice(id),
	254	AliHLTTPCSpacePointData::GetPatch(id),
	255	AliHLTTPCSpacePointData::GetNumber(id));
	256	}
	257	return iResult;
	258	}
	259
	260	int AliHLTTPCTrackGeometry::FillTrackPoints(AliHLTTrackGrid* pGrid) const
	261	{
	262	/// fill track points to index grid
	263	if (!pGrid) return -EINVAL;
	264	int iResult=0;
	265	for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
	266	tp!=TrackPoints().end() && iResult>=0; tp++) {
	267	AliHLTUInt32_t id=tp->GetId();
	268	iResult=pGrid->AddSpacePoint(GetTrackId(),
	269	AliHLTTPCSpacePointData::GetSlice(id),
	270	AliHLTTPCSpacePointData::GetPatch(id),
	271	AliHLTTPCSpacePointData::GetNumber(id));
	272	}
	273	return iResult;
	274	}
	275
9003c201	276	AliHLTSpacePointContainer* AliHLTTPCTrackGeometry::ConvertToSpacePoints(bool bAssociated) const
eb50ad6b	277	{
	278	/// create a collection of all points
	279	std::auto_ptr<AliHLTTPCSpacePointContainer> spacepoints(new AliHLTTPCSpacePointContainer);
	280	if (!spacepoints.get()) return NULL;
	281
dd2011bf	282	const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
eb50ad6b	283	unsigned i=0;
	284	while (i<trackPoints.size()) {
	285	// allocate buffer for all points, even though the buffer might not be filled
	286	// completely because of a partition change
	287	int nofPoints=trackPoints.size()-i;
	288	int blocksize=sizeof(AliHLTTPCClusterData)+nofPoints*sizeof(AliHLTTPCSpacePointData);
	289	AliHLTUInt8_t* pBuffer=spacepoints->Alloc(blocksize);
	290	if (!pBuffer) return NULL;
	291	AliHLTTPCClusterData* pClusterData=reinterpret_cast<AliHLTTPCClusterData*>(pBuffer);
	292	pClusterData->fSpacePointCnt=0;
	293	AliHLTTPCSpacePointData* pClusters=pClusterData->fSpacePoints;
	294	int currentSlice=-1;
	295	int currentPartition=-1;
	296	for (; i<trackPoints.size(); i++) {
9003c201	297	if (bAssociated && !trackPoints[i].HaveAssociatedSpacePoint()) continue;
eb50ad6b	298	AliHLTUInt32_t planeId=trackPoints[i].GetId();
	299	int slice=AliHLTTPCSpacePointData::GetSlice(planeId);
	300	int partition=AliHLTTPCSpacePointData::GetPatch(planeId);
	301	int number=AliHLTTPCSpacePointData::GetNumber(planeId);
	302	if ((currentSlice>=0 && currentSlice!=slice) \|\| (currentPartition>=0 && currentPartition!=partition)) {
	303	// change of partition or slice, need to go to next block
	304	// 2011-07-26 currently all spacepoints go into one block, if separated
	305	// blocks per partition are needed one has to leave the inner loop here
	306	// and set the data block specification below
	307	// Caution: not tested, only the last block seems to make it through
	308	//break;
	309	}
	310	currentSlice=slice;
	311	currentPartition=partition;
	312	pClusters[pClusterData->fSpacePointCnt].fX=GetPlaneR(planeId);
	313	pClusters[pClusterData->fSpacePointCnt].fY=trackPoints[i].GetU();
	314	pClusters[pClusterData->fSpacePointCnt].fZ=trackPoints[i].GetV();
	315	pClusters[pClusterData->fSpacePointCnt].fID=planeId;
	316	pClusters[pClusterData->fSpacePointCnt].fPadRow=AliHLTTPCTransform::GetFirstRow(partition)+number;
	317	pClusters[pClusterData->fSpacePointCnt].fSigmaY2=0.;
	318	pClusters[pClusterData->fSpacePointCnt].fSigmaZ2=0.;
	319	pClusters[pClusterData->fSpacePointCnt].fCharge=0;
	320	pClusters[pClusterData->fSpacePointCnt].fQMax=0;
	321	pClusters[pClusterData->fSpacePointCnt].fUsed=0;
	322	pClusters[pClusterData->fSpacePointCnt].fTrackN=0;
	323	pClusterData->fSpacePointCnt++;
	324	}
	325	AliHLTComponentBlockData bd;
	326	AliHLTComponent::FillBlockData(bd);
	327	bd.fPtr=pBuffer;
9003c201	328	bd.fSize=sizeof(AliHLTTPCClusterData)+pClusterData->fSpacePointCnt*sizeof(AliHLTTPCSpacePointData);
eb50ad6b	329	AliHLTComponent::SetDataType(bd.fDataType, "CLUSTERS", "TPC ");
bbe2c314	330	bd.fSpecification=kAliHLTVoidDataSpec;//AliHLTTPCDefinitions::EncodeDataSpecification(currentSlice, currentSlice, currentPartition, currentPartition);
eb50ad6b	331	spacepoints->AddInputBlock(&bd);
	332	}
	333
	334	return spacepoints.release();
	335	}
bbe2c314	336
c2bea8b3	337	const AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id) const
	338	{
	339	/// get raw track point of id
	340	const AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
	341	if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
	342	return p;
	343	}
	344
	345	AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id)
	346	{
	347	/// get raw track point of id
	348	AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
	349	if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
	350	return p;
	351	}
	352
bbe2c314	353	int AliHLTTPCTrackGeometry::FillRawResidual(int coordinate, TH2* histo, AliHLTSpacePointContainer* points) const
	354	{
	355	// fill residual histogram
	356	if (!histo \|\| !points) return -EINVAL;
	357	const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
	358	for (vector<AliHLTTrackPoint>::const_iterator trackpoint=trackPoints.begin();
	359	trackpoint!=trackPoints.end(); trackpoint++) {
6c6615c1	360	if (!trackpoint->HaveAssociatedSpacePoint()) continue;
	361	for (vector<AliHLTTrackSpacepoint>::const_iterator sp=(trackpoint->GetSpacepoints()).begin();
	362	sp!=(trackpoint->GetSpacepoints()).end(); sp++) {
	363	AliHLTUInt32_t spacepointId=sp->fId;
bbe2c314	364	vector<AliHLTTrackPoint>::const_iterator rawpoint=find(fRawTrackPoints.begin(), fRawTrackPoints.end(), trackpoint->GetId());
	365	if (rawpoint==fRawTrackPoints.end()) {
	366	HLTError("can not find track raw coordinates of track point 0x%08x", trackpoint->GetId());
	367	continue;
	368	}
	369	if (!points->Check(spacepointId)) {
	370	//HLTError("can not find associated space point 0x%08x of track point 0x%08x", spacepointId, trackpoint->GetId());
	371	continue;
	372	}
	373	float value=0.;
	374	if (coordinate==0) {
	375	value=rawpoint->GetU()-points->GetY(spacepointId);
	376	histo->Fill(GetPlaneR(trackpoint->GetId()), value);
	377	} else {
	378	value=rawpoint->GetV()-points->GetZ(spacepointId);
	379	//histo->Fill(GetPlaneR(trackpoint->GetId()), value);
	380	histo->Fill(rawpoint->GetV(), value);
	381	}
6c6615c1	382	}
bbe2c314	383	}
	384	return 0;
	385	}