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

// @(#) $Id$

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

#include "AliL3StandardIncludes.h"

#include "AliL3Logging.h"
#include "AliL3HoughClusterTransformer.h"
#include "AliL3MemHandler.h"
#include "AliL3SpacePointData.h"
#include "AliL3Transform.h"
#include "AliL3DigitData.h"
#include "AliL3Histogram.h"
#include "AliL3ClustFinderNew.h"

#if GCCVERSION == 3
using namespace std;
#endif

//_____________________________________________________________
// AliL3HoughClusterTransformer
//
// 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(AliL3HoughClusterTransformer)

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

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

AliL3HoughClusterTransformer::~AliL3HoughClusterTransformer()
{
  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 AliL3HoughClusterTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){}

void AliL3HoughClusterTransformer::DeleteHistograms()
{
  if(fParamSpace)
    {
      for(Int_t i=0; i<GetNEtaSegments(); i++)
	{
	  if(!fParamSpace[i]) continue;
	  delete fParamSpace[i];
	}
      delete [] fParamSpace;
    }
}

void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t pt_min,
					     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
  //pt_min = mimium Pt of track (corresponding to maximum kappa)
  //phi_min = mimimum phi0 (degrees)
  //phi_max = maximum phi0 (degrees)
    
  Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
  Double_t torad = AliL3Transform::Pi()/180;
  CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
}

void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
						    Int_t nybin,Float_t ymin,Float_t ymax)
{
  Char_t histname[256];
  
  fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
  for(Int_t i=0; i<GetNEtaSegments(); i++)
    {
      sprintf(histname,"paramspace_%d",i);
      fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
    }
#ifdef do_mc
  Int_t ncells = (nxbin+2)*(nybin+2);
  cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
  fTrackID = new TrackIndex*[GetNEtaSegments()];
  for(Int_t i=0; i<GetNEtaSegments(); i++)
    fTrackID[i] = new TrackIndex[ncells];
#endif
  
}

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

  if(!fParamSpace)
    {
      LOG(AliL3Log::kWarning,"AliL3HoughClusterTransformer::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
  AliL3Histogram *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(TrackIndex));
#endif
}

Int_t AliL3HoughClusterTransformer::GetEtaIndex(Double_t eta)
{
  //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 AliL3Histogram *AliL3HoughClusterTransformer::GetHistogram(Int_t eta_index)
{
  if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
    return 0;
  if(!fParamSpace[eta_index])
    return 0;
  return fParamSpace[eta_index];
}

Double_t AliL3HoughClusterTransformer::GetEta(Int_t eta_index,Int_t slice)
{
  Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
  Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
  if(slice>17) eta*=-1;
  return eta;
}

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

  const Int_t maxpoints=100000;
  const Int_t pointsize = maxpoints * sizeof(AliL3SpacePointData);

  fMemHandler = new AliL3MemHandler();
  fClusters = (AliL3SpacePointData*)fMemHandler->Allocate(pointsize);
  AliL3ClustFinderNew *cf = new AliL3ClustFinderNew();
  cf->InitSlice(0,GetPatch(),AliL3Transform::GetFirstRow(GetPatch()),AliL3Transform::GetLastRow(GetPatch()),maxpoints);
  cf->SetDeconv(kFALSE);
  cf->SetOutputArray(fClusters);
  cf->Read(1,GetDataPointer());
  cf->ProcessDigits();
  fNClusters = cf->GetNumberOfClusters();
  delete cf;
}

void AliL3HoughClusterTransformer::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(AliL3Log::kError,"AliL3HoughClusterTransformer::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 = AliL3Transform::GetEta(xyz);
      Int_t eta_index = GetEtaIndex(eta);
      if(eta_index < 0 || eta_index >= GetNEtaSegments())
	continue;
      
      AliL3Histogram *hist = fParamSpace[eta_index];
      
      //Do the transformation:
      Double_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); 
      Double_t phi = AliL3Transform::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[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
		  break;
	      if(c == MaxTrack-1) cerr<<"AliL3HoughClusterTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
	      fTrackID[eta_index][bin].fLabel[c] = label;
	      fTrackID[eta_index][bin].fNHits[c]++;
	    }
#endif
	}
    }
}

void AliL3HoughClusterTransformer::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(AliL3Log::kError,"AliL3HoughClusterTransformer::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 = AliL3Transform::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 = AliL3Transform::GetEta(xyz);
	  index2 = GetEtaIndex(eta);
	  if(index1 != index2) continue;
	  
	  AliL3Histogram *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<<"AliL3HoughClusterTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
		      fTrackID[index1][bin].fLabel[c] = label;
		      fTrackID[index1][bin].fNHits[c]++;
		    }
		}
	    }
#endif
	}
    }
}


Int_t AliL3HoughClusterTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
{

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