[u/mrichter/AliRoot.git] / HLT / hough / AliHLTHoughClusterTransformer.cxx

// @(#) $Id$

// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
//*-- Copyright &copy ALICE HLT Group

#include "AliHLTStandardIncludes.h"

#include "AliHLTLogging.h"
#include "AliHLTHoughClusterTransformer.h"
#include "AliHLTMemHandler.h"
#include "AliHLTSpacePointData.h"
#include "AliHLTTransform.h"
#include "AliHLTDigitData.h"
#include "AliHLTHistogram.h"
#include "AliHLTClustFinderNew.h"

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

//_____________________________________________________________
// AliHLTHoughClusterTransformer
//
// Hough transformation class. 
//
// This class performs the Hough transform on _clusters_, and not raw data.
// Therefore, it first finds the clusters using the HLT cluster finder, and
// then uses these as input for the transform.

ClassImp(AliHLTHoughClusterTransformer)

AliHLTHoughClusterTransformer::AliHLTHoughClusterTransformer()
{
  //Default constructor
  fParamSpace = 0;
  fNClusters = 0;
  fClusters = 0;
  fMemHandler = 0;
#ifdef do_mc
  fTrackID = 0;
#endif
}

AliHLTHoughClusterTransformer::AliHLTHoughClusterTransformer(Int_t slice,Int_t patch,Int_t netasegments) : AliHLTHoughBaseTransformer(slice,patch,netasegments)
{
  //Normal constructor
  fParamSpace = 0;
  fNClusters=0;
  fClusters=0;
  fMemHandler=0;
#ifdef do_mc
  fTrackID=0;
#endif
}

AliHLTHoughClusterTransformer::~AliHLTHoughClusterTransformer()
{
  //dtor
  DeleteHistograms();
  if(fMemHandler)
    delete fMemHandler;
#ifdef do_mc
  if(fTrackID)
    {
      for(Int_t i=0; i<GetNEtaSegments(); i++)
	{
	  if(!fTrackID[i]) continue;
	  delete fTrackID[i];
	}
      delete [] fTrackID;
    }
#endif
}

//void AliHLTHoughClusterTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t netasegments=100){}

void AliHLTHoughClusterTransformer::DeleteHistograms()
{
  //Delete hough space histograms
  if(fParamSpace)
    {
      for(Int_t i=0; i<GetNEtaSegments(); i++)
	{
	  if(!fParamSpace[i]) continue;
	  delete fParamSpace[i];
	}
      delete [] fParamSpace;
    }
}

void AliHLTHoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t ptmin,
					     Int_t nybin,Float_t phimin,Float_t phimax)
{
  //Create the histograms (parameter space).
  //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
  //The arguments give the range and binning; 
  //nxbin = #bins in kappa
  //nybin = #bins in phi0
  //ptmin = mimium Pt of track (corresponding to maximum kappa)
  //phimin = mimimum phi0 (degrees)
  //phimax = maximum phi0 (degrees)
    
  Double_t x = AliHLTTransform::GetBFact()*AliHLTTransform::GetBField()/ptmin;
  Double_t torad = AliHLTTransform::Pi()/180;
  CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
}

void AliHLTHoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
						    Int_t nybin,Float_t ymin,Float_t ymax)
{
  //Create histograms which contain hough space
  Char_t histname[256];
  
  fParamSpace = new AliHLTHistogram*[GetNEtaSegments()];
  for(Int_t i=0; i<GetNEtaSegments(); i++)
    {
      sprintf(histname,"paramspace_%d",i);
      fParamSpace[i] = new AliHLTHistogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
    }
#ifdef do_mc
  Int_t ncells = (nxbin+2)*(nybin+2);
  cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(AliHLTTrackIndex)<<" bytes to fTrackID"<<endl;
  fTrackID = new AliHLTTrackIndex*[GetNEtaSegments()];
  for(Int_t i=0; i<GetNEtaSegments(); i++)
    fTrackID[i] = new AliHLTTrackIndex[ncells];
#endif
  
}

void AliHLTHoughClusterTransformer::Reset()
{
  //Reset all the histograms

  if(!fParamSpace)
    {
      LOG(AliHLTLog::kWarning,"AliHLTHoughClusterTransformer::Reset","Histograms")
	<<"No histograms to reset"<<ENDLOG;
      return;
    }
  
  for(Int_t i=0; i<GetNEtaSegments(); i++)
    fParamSpace[i]->Reset();

  if(fMemHandler)
    delete fMemHandler;
  fNClusters=0;
  fClusters=0;
#ifdef do_mc
  AliHLTHistogram *hist = fParamSpace[0];
  Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
  for(Int_t i=0; i<GetNEtaSegments(); i++)
    memset(fTrackID[i],0,ncells*sizeof(AliHLTTrackIndex));
#endif
}

Int_t AliHLTHoughClusterTransformer::GetEtaIndex(Double_t eta) const
{
  //Return the histogram index of the corresponding eta. 
  
  Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
  Double_t index = (eta-GetEtaMin())/etaslice;
  return (Int_t)index;
}

inline AliHLTHistogram *AliHLTHoughClusterTransformer::GetHistogram(Int_t etaindex)
{
  //Returns the histogram which correspond to a given eta slice
  if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
    return 0;
  if(!fParamSpace[etaindex])
    return 0;
  return fParamSpace[etaindex];
}

Double_t AliHLTHoughClusterTransformer::GetEta(Int_t etaindex,Int_t slice) const
{
  //Returns eta associated with given eta slice
  Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
  Double_t eta=(Double_t)((etaindex+0.5)*etaslice);
  if(slice>17) eta*=-1;
  return eta;
}

void AliHLTHoughClusterTransformer::FindClusters()
{
  //Find the clusters
  if(!GetDataPointer())
    {
      cerr<<"AliHLTHoughClusterTransformer::FindClusters : Zero data pointer"<<endl;
      return;
    }

  const Int_t kMaxpoints=100000;
  const Int_t kPointsize = kMaxpoints * sizeof(AliHLTSpacePointData);

  fMemHandler = new AliHLTMemHandler();
  fClusters = (AliHLTSpacePointData*)fMemHandler->Allocate(kPointsize);
  AliHLTClustFinderNew *cf = new AliHLTClustFinderNew();
  cf->InitSlice(0,GetPatch(),AliHLTTransform::GetFirstRow(GetPatch()),AliHLTTransform::GetLastRow(GetPatch()),kMaxpoints);
  cf->SetDeconv(kFALSE);
  cf->SetOutputArray(fClusters);
  cf->Read(1,GetDataPointer());
  cf->ProcessDigits();
  fNClusters = cf->GetNumberOfClusters();
  delete cf;
}

void AliHLTHoughClusterTransformer::TransformCircle()
{
  //Transform the input data with a circle HT.
  //The function loops over all the data, and transforms each cluster with the equations:
  // 
  //kappa = 2/r*sin(phi - phi0)
  //
  //where r = sqrt(x*x +y*y), and phi = arctan(y/x)
  //
  //Each cluster then transforms into a curve in the (kappa,phi0)-space. In order to find
  //which histogram in which the pixel should be transformed, the eta-value is calcluated
  //and the proper histogram index is found by GetEtaIndex(eta).

  FindClusters();
  if(!fClusters)
    {
      LOG(AliHLTLog::kError,"AliHLTHoughClusterTransformer::TransformCircle","Data")
	<<"No input data "<<ENDLOG;
      return;
    }
  
  Float_t xyz[3];
  for(Int_t i=0; i<fNClusters; i++)
    {
      xyz[0] = fClusters[i].fX;
      xyz[1] = fClusters[i].fY;
      xyz[2] = fClusters[i].fZ;
      Double_t eta = AliHLTTransform::GetEta(xyz);
      Int_t etaindex = GetEtaIndex(eta);
      if(etaindex < 0 || etaindex >= GetNEtaSegments())
	continue;
      
      AliHLTHistogram *hist = fParamSpace[etaindex];
      
      //Do the transformation:
      Double_t r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); 
      Double_t phi = AliHLTTransform::GetPhi(xyz);
      
      //Fill the histogram along the phirange
      for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
	{
	  Double_t phi0 = hist->GetBinCenterY(b);
	  Double_t kappa = 2*sin(phi - phi0)/r;
	  hist->Fill(kappa,phi0,1);
#ifdef do_mc
	  Int_t bin = hist->FindBin(kappa,phi0);
	  for(Int_t t=0; t<3; t++)
	    {
	      Int_t label = fClusters[i].fTrackID[t];
	      if(label < 0) break;
	      UInt_t c;
	      for(c=0; c<MaxTrack; c++)
		if(fTrackID[etaindex][bin].fLabel[c] == label || fTrackID[etaindex][bin].fNHits[c] == 0)
		  break;
	      if(c == MaxTrack-1) cerr<<"AliHLTHoughClusterTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
	      fTrackID[etaindex][bin].fLabel[c] = label;
	      fTrackID[etaindex][bin].fNHits[c]++;
	    }
#endif
	}
    }
}

void AliHLTHoughClusterTransformer::TransformCircleC(Int_t */*row_range*/,Int_t /*every*/)
{
  //Circle transform, using combinations of every 2 points lying
  //on different padrows and within the same etaslice.
  
  FindClusters();
  if(!fClusters)
    LOG(AliHLTLog::kError,"AliHLTHoughClusterTransformer::TransformCircleC","Data")
      <<"No input data "<<ENDLOG;
  
  Float_t xyz[3];
  Double_t eta,r1,phi1,r2,phi2,kappa,psi;
  Int_t index1,index2;
  for(Int_t i=0; i<fNClusters; i++)
    {
      xyz[0] = fClusters[i].fX;
      xyz[1] = fClusters[i].fY;
      xyz[2] = fClusters[i].fZ;
      r1 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
      phi1 = atan2(xyz[1],xyz[0]);
      eta = AliHLTTransform::GetEta(xyz);
      index1 = GetEtaIndex(eta);
      if(index1 < 0 || index1 >= GetNEtaSegments()) continue;
      
      for(Int_t j=i+1; j<fNClusters; j++)
	{
	  if(fClusters[j].fPadRow == fClusters[i].fPadRow) continue;
	  xyz[0] = fClusters[j].fX;
	  xyz[1] = fClusters[j].fY;
	  xyz[2] = fClusters[j].fZ;
	  eta = AliHLTTransform::GetEta(xyz);
	  index2 = GetEtaIndex(eta);
	  if(index1 != index2) continue;
	  
	  AliHLTHistogram *hist = fParamSpace[index1];
	  
	  r2 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
	  phi2 = atan2(xyz[1],xyz[0]);
	  
	  psi = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
	  kappa = 2*sin(phi1-psi) / r1;
	  hist->Fill(kappa,psi,1);
#ifdef do_mc
	  Int_t bin = hist->FindBin(kappa,psi);
	  for(Int_t l=0; l<3; l++)
	    {
	      for(Int_t m=0; m<3; m++)
		{
		  if(fClusters[i].fTrackID[l] == fClusters[j].fTrackID[m])
		    {
		      Int_t label = fClusters[i].fTrackID[l];
		      if(label < 0) continue;
		      UInt_t c;
		      for(c=0; c<MaxTrack; c++)
			if(fTrackID[index1][bin].fLabel[c] == label || fTrackID[index1][bin].fNHits[c] == 0)
			  break;
		      if(c == MaxTrack-1) cerr<<"AliHLTHoughClusterTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
		      fTrackID[index1][bin].fLabel[c] = label;
		      fTrackID[index1][bin].fNHits[c]++;
		    }
		}
	    }
#endif
	}
    }
}


#ifdef do_mc
Int_t AliHLTHoughClusterTransformer::GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi) const
{
  //Returns the mc label for a given bin in the hough space
  if(etaindex < 0 || etaindex > GetNEtaSegments())
    {
      cerr<<"AliHLTHoughClusterTransformer::GetTrackID : Wrong etaindex "<<etaindex<<endl;
      return -1;
    }
  AliHLTHistogram *hist = fParamSpace[etaindex];
  Int_t bin = hist->FindBin(kappa,psi);
  Int_t label=-1;
  Int_t max=0;
  for(UInt_t i=0; i<MaxTrack; i++)
    {
      Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
      if(nhits == 0) break;
      if(nhits > max)
	{
	  max = nhits;
	  label = fTrackID[etaindex][bin].fLabel[i];
	}
    }
  return label;
#else
  Int_t AliHLTHoughClusterTransformer::GetTrackID(Int_t /*etaindex*/,Double_t /*kappa*/,Double_t /*psi*/) const
{
  // Does nothing if do_mc undefinde
  cout<<"AliHLTHoughClusterTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
  return -1;
#endif
}
Commit	Line	Data
3e87ef69	1	// @(#) $Id$
28ac1a53	2
28ac1a53	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
3e87ef69	4	//*-- Copyright &copy ALICE HLT Group
28ac1a53	5
4aa41877	6	#include "AliHLTStandardIncludes.h"
e06900d5	7
4aa41877	8	#include "AliHLTLogging.h"
	9	#include "AliHLTHoughClusterTransformer.h"
	10	#include "AliHLTMemHandler.h"
	11	#include "AliHLTSpacePointData.h"
	12	#include "AliHLTTransform.h"
	13	#include "AliHLTDigitData.h"
	14	#include "AliHLTHistogram.h"
	15	#include "AliHLTClustFinderNew.h"
28ac1a53	16
c6747af6	17	#if __GNUC__ >= 3
e06900d5	18	using namespace std;
	19	#endif
	20
28ac1a53	21	//_____________________________________________________________
4aa41877	22	// AliHLTHoughClusterTransformer
28ac1a53	23	//
	24	// Hough transformation class.
	25	//
	26	// This class performs the Hough transform on _clusters_, and not raw data.
	27	// Therefore, it first finds the clusters using the HLT cluster finder, and
	28	// then uses these as input for the transform.
	29
4aa41877	30	ClassImp(AliHLTHoughClusterTransformer)
28ac1a53	31
4aa41877	32	AliHLTHoughClusterTransformer::AliHLTHoughClusterTransformer()
28ac1a53	33	{
	34	//Default constructor
	35	fParamSpace = 0;
	36	fNClusters = 0;
	37	fClusters = 0;
	38	fMemHandler = 0;
	39	#ifdef do_mc
	40	fTrackID = 0;
	41	#endif
	42	}
	43
4aa41877	44	AliHLTHoughClusterTransformer::AliHLTHoughClusterTransformer(Int_t slice,Int_t patch,Int_t netasegments) : AliHLTHoughBaseTransformer(slice,patch,netasegments)
28ac1a53	45	{
	46	//Normal constructor
	47	fParamSpace = 0;
	48	fNClusters=0;
	49	fClusters=0;
	50	fMemHandler=0;
	51	#ifdef do_mc
	52	fTrackID=0;
	53	#endif
	54	}
	55
4aa41877	56	AliHLTHoughClusterTransformer::~AliHLTHoughClusterTransformer()
28ac1a53	57	{
c62b480b	58	//dtor
28ac1a53	59	DeleteHistograms();
	60	if(fMemHandler)
	61	delete fMemHandler;
	62	#ifdef do_mc
	63	if(fTrackID)
	64	{
	65	for(Int_t i=0; i<GetNEtaSegments(); i++)
	66	{
	67	if(!fTrackID[i]) continue;
	68	delete fTrackID[i];
	69	}
	70	delete [] fTrackID;
	71	}
	72	#endif
	73	}
	74
4aa41877	75	//void AliHLTHoughClusterTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t netasegments=100){}
28ac1a53	76
4aa41877	77	void AliHLTHoughClusterTransformer::DeleteHistograms()
28ac1a53	78	{
c62b480b	79	//Delete hough space histograms
28ac1a53	80	if(fParamSpace)
	81	{
	82	for(Int_t i=0; i<GetNEtaSegments(); i++)
	83	{
	84	if(!fParamSpace[i]) continue;
	85	delete fParamSpace[i];
	86	}
	87	delete [] fParamSpace;
	88	}
	89	}
	90
4aa41877	91	void AliHLTHoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t ptmin,
b2a02bce	92	Int_t nybin,Float_t phimin,Float_t phimax)
28ac1a53	93	{
	94	//Create the histograms (parameter space).
	95	//These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
	96	//The arguments give the range and binning;
	97	//nxbin = #bins in kappa
	98	//nybin = #bins in phi0
c62b480b	99	//ptmin = mimium Pt of track (corresponding to maximum kappa)
	100	//phimin = mimimum phi0 (degrees)
	101	//phimax = maximum phi0 (degrees)
28ac1a53	102
4aa41877	103	Double_t x = AliHLTTransform::GetBFact()*AliHLTTransform::GetBField()/ptmin;
4aa41877	104	Double_t torad = AliHLTTransform::Pi()/180;
28ac1a53	105	CreateHistograms(nxbin,-1.x,x,nybin,phimintorad,phimax*torad);
	106	}
	107
4aa41877	108	void AliHLTHoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
b2a02bce	109	Int_t nybin,Float_t ymin,Float_t ymax)
28ac1a53	110	{
c62b480b	111	//Create histograms which contain hough space
28ac1a53	112	Char_t histname[256];
28ac1a53	113
4aa41877	114	fParamSpace = new AliHLTHistogram*[GetNEtaSegments()];
28ac1a53	115	for(Int_t i=0; i<GetNEtaSegments(); i++)
	116	{
	117	sprintf(histname,"paramspace_%d",i);
4aa41877	118	fParamSpace[i] = new AliHLTHistogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
28ac1a53	119	}
	120	#ifdef do_mc
	121	Int_t ncells = (nxbin+2)*(nybin+2);
4aa41877	122	cout<<"Allocating "<<GetNEtaSegments()ncellssizeof(AliHLTTrackIndex)<<" bytes to fTrackID"<<endl;
4aa41877	123	fTrackID = new AliHLTTrackIndex*[GetNEtaSegments()];
28ac1a53	124	for(Int_t i=0; i<GetNEtaSegments(); i++)
4aa41877	125	fTrackID[i] = new AliHLTTrackIndex[ncells];
28ac1a53	126	#endif
	127
	128	}
	129
4aa41877	130	void AliHLTHoughClusterTransformer::Reset()
28ac1a53	131	{
	132	//Reset all the histograms
	133
	134	if(!fParamSpace)
	135	{
4aa41877	136	LOG(AliHLTLog::kWarning,"AliHLTHoughClusterTransformer::Reset","Histograms")
28ac1a53	137	<<"No histograms to reset"<<ENDLOG;
	138	return;
	139	}
	140
	141	for(Int_t i=0; i<GetNEtaSegments(); i++)
	142	fParamSpace[i]->Reset();
	143
	144	if(fMemHandler)
	145	delete fMemHandler;
	146	fNClusters=0;
	147	fClusters=0;
	148	#ifdef do_mc
4aa41877	149	AliHLTHistogram *hist = fParamSpace[0];
28ac1a53	150	Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
28ac1a53	151	for(Int_t i=0; i<GetNEtaSegments(); i++)
4aa41877	152	memset(fTrackID[i],0,ncells*sizeof(AliHLTTrackIndex));
28ac1a53	153	#endif
	154	}
	155
4aa41877	156	Int_t AliHLTHoughClusterTransformer::GetEtaIndex(Double_t eta) const
28ac1a53	157	{
	158	//Return the histogram index of the corresponding eta.
	159
	160	Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
	161	Double_t index = (eta-GetEtaMin())/etaslice;
	162	return (Int_t)index;
	163	}
	164
4aa41877	165	inline AliHLTHistogram *AliHLTHoughClusterTransformer::GetHistogram(Int_t etaindex)
28ac1a53	166	{
c62b480b	167	//Returns the histogram which correspond to a given eta slice
c62b480b	168	if(!fParamSpace \|\| etaindex >= GetNEtaSegments() \|\| etaindex < 0)
28ac1a53	169	return 0;
c62b480b	170	if(!fParamSpace[etaindex])
28ac1a53	171	return 0;
c62b480b	172	return fParamSpace[etaindex];
28ac1a53	173	}
28ac1a53	174
4aa41877	175	Double_t AliHLTHoughClusterTransformer::GetEta(Int_t etaindex,Int_t slice) const
28ac1a53	176	{
c62b480b	177	//Returns eta associated with given eta slice
	178	Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
	179	Double_t eta=(Double_t)((etaindex+0.5)*etaslice);
28ac1a53	180	if(slice>17) eta*=-1;
	181	return eta;
	182	}
	183
4aa41877	184	void AliHLTHoughClusterTransformer::FindClusters()
28ac1a53	185	{
	186	//Find the clusters
	187	if(!GetDataPointer())
	188	{
4aa41877	189	cerr<<"AliHLTHoughClusterTransformer::FindClusters : Zero data pointer"<<endl;
28ac1a53	190	return;
	191	}
	192
c62b480b	193	const Int_t kMaxpoints=100000;
4aa41877	194	const Int_t kPointsize = kMaxpoints * sizeof(AliHLTSpacePointData);
28ac1a53	195
4aa41877	196	fMemHandler = new AliHLTMemHandler();
	197	fClusters = (AliHLTSpacePointData*)fMemHandler->Allocate(kPointsize);
	198	AliHLTClustFinderNew *cf = new AliHLTClustFinderNew();
	199	cf->InitSlice(0,GetPatch(),AliHLTTransform::GetFirstRow(GetPatch()),AliHLTTransform::GetLastRow(GetPatch()),kMaxpoints);
28ac1a53	200	cf->SetDeconv(kFALSE);
	201	cf->SetOutputArray(fClusters);
	202	cf->Read(1,GetDataPointer());
	203	cf->ProcessDigits();
	204	fNClusters = cf->GetNumberOfClusters();
	205	delete cf;
	206	}
	207
4aa41877	208	void AliHLTHoughClusterTransformer::TransformCircle()
28ac1a53	209	{
	210	//Transform the input data with a circle HT.
	211	//The function loops over all the data, and transforms each cluster with the equations:
	212	//
c62b480b	213	//kappa = 2/r*sin(phi - phi0)
28ac1a53	214	//
c62b480b	215	//where r = sqrt(xx +yy), and phi = arctan(y/x)
28ac1a53	216	//
	217	//Each cluster then transforms into a curve in the (kappa,phi0)-space. In order to find
	218	//which histogram in which the pixel should be transformed, the eta-value is calcluated
	219	//and the proper histogram index is found by GetEtaIndex(eta).
	220
	221	FindClusters();
	222	if(!fClusters)
	223	{
4aa41877	224	LOG(AliHLTLog::kError,"AliHLTHoughClusterTransformer::TransformCircle","Data")
28ac1a53	225	<<"No input data "<<ENDLOG;
	226	return;
	227	}
	228
	229	Float_t xyz[3];
	230	for(Int_t i=0; i<fNClusters; i++)
	231	{
	232	xyz[0] = fClusters[i].fX;
	233	xyz[1] = fClusters[i].fY;
	234	xyz[2] = fClusters[i].fZ;
4aa41877	235	Double_t eta = AliHLTTransform::GetEta(xyz);
c62b480b	236	Int_t etaindex = GetEtaIndex(eta);
c62b480b	237	if(etaindex < 0 \|\| etaindex >= GetNEtaSegments())
28ac1a53	238	continue;
28ac1a53	239
4aa41877	240	AliHLTHistogram *hist = fParamSpace[etaindex];
28ac1a53	241
28ac1a53	242	//Do the transformation:
c62b480b	243	Double_t r = sqrt(xyz[0]xyz[0] + xyz[1]xyz[1]);
4aa41877	244	Double_t phi = AliHLTTransform::GetPhi(xyz);
28ac1a53	245
	246	//Fill the histogram along the phirange
	247	for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
	248	{
	249	Double_t phi0 = hist->GetBinCenterY(b);
c62b480b	250	Double_t kappa = 2*sin(phi - phi0)/r;
28ac1a53	251	hist->Fill(kappa,phi0,1);
	252	#ifdef do_mc
	253	Int_t bin = hist->FindBin(kappa,phi0);
	254	for(Int_t t=0; t<3; t++)
	255	{
	256	Int_t label = fClusters[i].fTrackID[t];
	257	if(label < 0) break;
	258	UInt_t c;
	259	for(c=0; c<MaxTrack; c++)
c62b480b	260	if(fTrackID[etaindex][bin].fLabel[c] == label \|\| fTrackID[etaindex][bin].fNHits[c] == 0)
28ac1a53	261	break;
4aa41877	262	if(c == MaxTrack-1) cerr<<"AliHLTHoughClusterTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
c62b480b	263	fTrackID[etaindex][bin].fLabel[c] = label;
c62b480b	264	fTrackID[etaindex][bin].fNHits[c]++;
28ac1a53	265	}
	266	#endif
	267	}
	268	}
	269	}
	270
4aa41877	271	void AliHLTHoughClusterTransformer::TransformCircleC(Int_t /row_range/,Int_t /every*/)
28ac1a53	272	{
	273	//Circle transform, using combinations of every 2 points lying
	274	//on different padrows and within the same etaslice.
	275
	276	FindClusters();
	277	if(!fClusters)
4aa41877	278	LOG(AliHLTLog::kError,"AliHLTHoughClusterTransformer::TransformCircleC","Data")
28ac1a53	279	<<"No input data "<<ENDLOG;
	280
	281	Float_t xyz[3];
	282	Double_t eta,r1,phi1,r2,phi2,kappa,psi;
	283	Int_t index1,index2;
	284	for(Int_t i=0; i<fNClusters; i++)
	285	{
	286	xyz[0] = fClusters[i].fX;
	287	xyz[1] = fClusters[i].fY;
	288	xyz[2] = fClusters[i].fZ;
	289	r1 = sqrt(xyz[0]xyz[0] + xyz[1]xyz[1]);
	290	phi1 = atan2(xyz[1],xyz[0]);
4aa41877	291	eta = AliHLTTransform::GetEta(xyz);
28ac1a53	292	index1 = GetEtaIndex(eta);
	293	if(index1 < 0 \|\| index1 >= GetNEtaSegments()) continue;
	294
	295	for(Int_t j=i+1; j<fNClusters; j++)
	296	{
	297	if(fClusters[j].fPadRow == fClusters[i].fPadRow) continue;
	298	xyz[0] = fClusters[j].fX;
	299	xyz[1] = fClusters[j].fY;
	300	xyz[2] = fClusters[j].fZ;
4aa41877	301	eta = AliHLTTransform::GetEta(xyz);
28ac1a53	302	index2 = GetEtaIndex(eta);
	303	if(index1 != index2) continue;
	304
4aa41877	305	AliHLTHistogram *hist = fParamSpace[index1];
28ac1a53	306
	307	r2 = sqrt(xyz[0]xyz[0] + xyz[1]xyz[1]);
	308	phi2 = atan2(xyz[1],xyz[0]);
	309
	310	psi = atan( (r2sin(phi1) - r1sin(phi2)) / (r2cos(phi1) - r1cos(phi2)) );
	311	kappa = 2*sin(phi1-psi) / r1;
	312	hist->Fill(kappa,psi,1);
	313	#ifdef do_mc
	314	Int_t bin = hist->FindBin(kappa,psi);
	315	for(Int_t l=0; l<3; l++)
	316	{
	317	for(Int_t m=0; m<3; m++)
	318	{
	319	if(fClusters[i].fTrackID[l] == fClusters[j].fTrackID[m])
	320	{
	321	Int_t label = fClusters[i].fTrackID[l];
	322	if(label < 0) continue;
	323	UInt_t c;
	324	for(c=0; c<MaxTrack; c++)
	325	if(fTrackID[index1][bin].fLabel[c] == label \|\| fTrackID[index1][bin].fNHits[c] == 0)
	326	break;
4aa41877	327	if(c == MaxTrack-1) cerr<<"AliHLTHoughClusterTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
28ac1a53	328	fTrackID[index1][bin].fLabel[c] = label;
	329	fTrackID[index1][bin].fNHits[c]++;
	330	}
	331	}
	332	}
	333	#endif
	334	}
	335	}
	336	}
	337
	338
3cb5c013	339	#ifdef do_mc
4aa41877	340	Int_t AliHLTHoughClusterTransformer::GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi) const
28ac1a53	341	{
c62b480b	342	//Returns the mc label for a given bin in the hough space
c62b480b	343	if(etaindex < 0 \|\| etaindex > GetNEtaSegments())
28ac1a53	344	{
4aa41877	345	cerr<<"AliHLTHoughClusterTransformer::GetTrackID : Wrong etaindex "<<etaindex<<endl;
28ac1a53	346	return -1;
28ac1a53	347	}
4aa41877	348	AliHLTHistogram *hist = fParamSpace[etaindex];
28ac1a53	349	Int_t bin = hist->FindBin(kappa,psi);
	350	Int_t label=-1;
	351	Int_t max=0;
	352	for(UInt_t i=0; i<MaxTrack; i++)
	353	{
c62b480b	354	Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
28ac1a53	355	if(nhits == 0) break;
	356	if(nhits > max)
	357	{
	358	max = nhits;
c62b480b	359	label = fTrackID[etaindex][bin].fLabel[i];
28ac1a53	360	}
	361	}
	362	return label;
3cb5c013	363	#else
4aa41877	364	Int_t AliHLTHoughClusterTransformer::GetTrackID(Int_t /etaindex/,Double_t /kappa/,Double_t /psi/) const
3cb5c013	365	{
3cb5c013	366	// Does nothing if do_mc undefinde
4aa41877	367	cout<<"AliHLTHoughClusterTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
28ac1a53	368	return -1;
3cb5c013	369	#endif
28ac1a53	370	}
28ac1a53	371