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

// @(#) $Id$

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

#include "AliL3StandardIncludes.h"

#include "AliL3Logging.h"
#include "AliL3HoughTransformerNew.h"
#include "AliL3MemHandler.h"
#include "AliL3Transform.h"
#include "AliL3DigitData.h"
#include "AliL3HistogramAdaptive.h"
#include "AliL3TrackArray.h"
#include "AliL3HoughTrack.h"
#include "AliL3SpacePointData.h"
#include "AliL3Vertex.h"
#include "AliL3Fitter.h"

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

//_____________________________________________________________
// AliL3HoughTransformerNew
//
// Hough transformation class
//

ClassImp(AliL3HoughTransformerNew)

AliL3HoughTransformerNew::AliL3HoughTransformerNew()
{
  //default ctor
  fParamSpace3D=0;
}

AliL3HoughTransformerNew::AliL3HoughTransformerNew(Int_t slice,Int_t patch,Int_t netasegments) : AliL3HoughTransformer(slice,patch,netasegments)
{
  //normal ctor
  fParamSpace3D=0;
}
AliL3HoughTransformerNew::~AliL3HoughTransformerNew()
{
  //dtor
  if(fParamSpace3D)
    delete fParamSpace3D;
}

void AliL3HoughTransformerNew::CreateHistograms(Int_t nxbins,Float_t xlow,Float_t xup,
						Int_t nybins,Float_t ylow,Float_t yup,
						Int_t nzbins,Float_t zlow,Float_t zup)
{
  //Create the histogram which contain the hough space
  Char_t name[1024];
  sprintf(name,"paramspace_%p",(void*)this);
  fParamSpace3D = new TH3F(name,"",nxbins,xlow,xup,nybins,ylow,yup,nzbins,zlow,zup);
}

void AliL3HoughTransformerNew::Reset()
{
  //Reset Hough space
  fParamSpace3D->Reset();
}


void AliL3HoughTransformerNew::TransformLine(Int_t *rowrange,Float_t *phirange)
{
  //Hough Transform
  AliL3DigitRowData *tempPt = GetDataPointer();
  if(!tempPt)
    {
      LOG(AliL3Log::kError,"AliL3HoughTransformerNew::TransformLine","Data")
	<<"No input data "<<ENDLOG;
      return;
    }
  
  TH3 *hist = fParamSpace3D;
  for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
    {
      AliL3DigitData *digPt = tempPt->fDigitData;
      if(i < rowrange[0])
	{
	  AliL3MemHandler::UpdateRowPointer(tempPt);
	  continue;
	}
      else if(i > rowrange[1])
	break;
      if(i != (Int_t)tempPt->fRow)
	{
	  printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
	  continue;
	}
      for(UInt_t j=0; j<tempPt->fNDigit; j++)
	{
	  UShort_t charge = digPt[j].fCharge;
	  UChar_t pad = digPt[j].fPad;
	  UShort_t time = digPt[j].fTime;
	  if(charge < GetLowerThreshold() || charge > GetUpperThreshold()) continue;
	  
	  Int_t sector,row;
	  Float_t xyz[3];
	  AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
	  AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
	  
	  Float_t phi = AliL3Transform::GetPhi(xyz);
	  if(phi < phirange[0] || phi > phirange[1])
	    continue;

	  xyz[0] = xyz[0] - AliL3Transform::Row2X(rowrange[0]);
	  Float_t x = xyz[0] + AliL3Transform::Row2X(rowrange[0]);
	  Float_t r = sqrt(x*x + xyz[1]*xyz[1]);
	  Float_t delta = atan(xyz[2]/r);
	  for(Int_t xbin=hist->GetXaxis()->GetFirst(); xbin<=hist->GetXaxis()->GetLast(); xbin++)
	    {
	      Float_t theta = hist->GetXaxis()->GetBinCenter(xbin);
	      Float_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
	      hist->Fill(theta,rho,delta,charge);
	    }
	}
      AliL3MemHandler::UpdateRowPointer(tempPt);
    }
  
}

struct AliL3Digit {
  Int_t fRow;//padrow index
  Int_t fCharge;//digits charge
  Float_t fY;//Y coordinate of the digit
  Float_t fZ;//Z coordinate of the digit
  AliL3Digit *fNext;//pointer to the next digit
};

struct AliL3EtaContainer {
  void *fFirst;//pointer to the first digit in a sequence
  void *fLast;//pointer to the last digit in a sequence
};

void AliL3HoughTransformerNew::TransformLineC(Int_t *rowrange,Float_t *phirange)
{
  //Hough Transform
  AliL3DigitRowData *tempPt = GetDataPointer();
  if(!tempPt)
    LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
      <<"No input data "<<ENDLOG;
  
  Int_t counter=0;
  for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
    {
      counter += tempPt->fNDigit;
      AliL3MemHandler::UpdateRowPointer(tempPt);
    }
  
  Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
  AliL3EtaContainer *etaPt = new AliL3EtaContainer[bound];
  memset(etaPt,0,bound*sizeof(AliL3EtaContainer));  
  
  AliL3Digit *digits = new AliL3Digit[counter];
  cout<<"Allocating "<<counter*sizeof(AliL3Digit)<<" bytes to digitsarray"<<endl;
  memset(digits,0,counter*sizeof(AliL3Digit));

  Int_t sector,row;
  Float_t xyz[3];
  
  counter=0;
  tempPt = GetDataPointer();
  
  cout<<"Calculating digits in patch "<<GetPatch()<<endl;
  for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
    {
      AliL3DigitData *digPt = tempPt->fDigitData;
      for(UInt_t di=0; di<tempPt->fNDigit; di++)
	{
	  Int_t charge = digPt[di].fCharge;
	  Int_t pad = digPt[di].fPad;
	  Int_t time = digPt[di].fTime;
	  AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
	  AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
	  Double_t eta = AliL3Transform::GetEta(xyz);
	  
	  Float_t phi = atan2(xyz[1],xyz[0]);
	  if(phi < phirange[0] || phi > phirange[1]) continue;
	  
	  digits[counter].fRow = i;
	  digits[counter].fY = xyz[1];
	  digits[counter].fZ = xyz[2];
	  digits[counter].fCharge = charge;
	  
	  Int_t etaindex = GetEtaIndex(eta);
	  Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
	  
	  if(index > 0 && index < bound) 
	    {
	      if(etaPt[index].fFirst == 0)
		etaPt[index].fFirst = (void*)(&digits[counter]);
	      else
		((AliL3Digit*)(etaPt[index].fLast))->fNext = &digits[counter];
	      etaPt[index].fLast = (void*)(&digits[counter]);
	    }
	  counter++;
	}
      AliL3MemHandler::UpdateRowPointer(tempPt);
    }
  
  cout<<"Doing the combinatorics"<<endl;
  
  AliL3Digit *dPt1,*dPt2;
  TH3 *hist = fParamSpace3D;
  for(Int_t e=0; e<GetNEtaSegments(); e++)
    {
      for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
	{
	  Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
	  
	  for(dPt1 = (AliL3Digit*)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3Digit*)dPt1->fNext)
	    {
	      for(Int_t j=i+1; j<=rowrange[1]; j++)
		{
		  Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
		  
		  for(dPt2 = (AliL3Digit*)etaPt[index2].fFirst; dPt2 != 0; dPt2 = (AliL3Digit*)dPt2->fNext)
		    {
		      if(dPt1->fRow == dPt2->fRow)
			{
			  cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
			  exit(5);
			}
		      
		      //Do the transform:
		      Float_t x1 = AliL3Transform::Row2X(dPt1->fRow) - AliL3Transform::Row2X(rowrange[0]);
		      Float_t x2 = AliL3Transform::Row2X(dPt2->fRow) - AliL3Transform::Row2X(rowrange[0]);
		      Float_t y1 = dPt1->fY;
		      Float_t y2 = dPt2->fY;
		      Float_t theta = atan2(x2-x1,y1-y2);
		      Float_t rho = x1*cos(theta)+y1*sin(theta);
		      Float_t r1 = sqrt(pow(AliL3Transform::Row2X(dPt1->fRow),2) + pow(y1,2));
		      Float_t delta1 = atan2(dPt1->fZ,r1);
		      Float_t r2 = sqrt(pow(AliL3Transform::Row2X(dPt2->fRow),2) + pow(y2,2));
		      Float_t delta2 = atan2(dPt2->fZ,r2);
		      Float_t delta = (dPt1->fCharge*delta1+dPt2->fCharge*delta2)/(dPt1->fCharge+dPt2->fCharge);
		      hist->Fill(theta,rho,delta,dPt1->fCharge+dPt2->fCharge);
		    }
		}
	    }
	}
    }
  
  cout<<"done"<<endl;
  delete [] etaPt;
  delete [] digits;
  
}
Commit	Line	Data
3e87ef69	1	// @(#) $Id$
	2
	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
	4	//*-- Copyright &copy ALICE HLT Group
	5
	6	#include "AliL3StandardIncludes.h"
	7
	8	#include "AliL3Logging.h"
	9	#include "AliL3HoughTransformerNew.h"
	10	#include "AliL3MemHandler.h"
	11	#include "AliL3Transform.h"
	12	#include "AliL3DigitData.h"
	13	#include "AliL3HistogramAdaptive.h"
	14	#include "AliL3TrackArray.h"
	15	#include "AliL3HoughTrack.h"
	16	#include "AliL3SpacePointData.h"
	17	#include "AliL3Vertex.h"
	18	#include "AliL3Fitter.h"
	19
c6747af6	20	#if __GNUC__ >= 3
3e87ef69	21	using namespace std;
	22	#endif
	23
	24	//_____________________________________________________________
	25	// AliL3HoughTransformerNew
	26	//
	27	// Hough transformation class
	28	//
	29
	30	ClassImp(AliL3HoughTransformerNew)
	31
	32	AliL3HoughTransformerNew::AliL3HoughTransformerNew()
	33	{
c62b480b	34	//default ctor
c62b480b	35	fParamSpace3D=0;
3e87ef69	36	}
	37
	38	AliL3HoughTransformerNew::AliL3HoughTransformerNew(Int_t slice,Int_t patch,Int_t netasegments) : AliL3HoughTransformer(slice,patch,netasegments)
	39	{
c62b480b	40	//normal ctor
c62b480b	41	fParamSpace3D=0;
3e87ef69	42	}
	43	AliL3HoughTransformerNew::~AliL3HoughTransformerNew()
	44	{
c62b480b	45	//dtor
	46	if(fParamSpace3D)
	47	delete fParamSpace3D;
3e87ef69	48	}
	49
	50	void AliL3HoughTransformerNew::CreateHistograms(Int_t nxbins,Float_t xlow,Float_t xup,
	51	Int_t nybins,Float_t ylow,Float_t yup,
	52	Int_t nzbins,Float_t zlow,Float_t zup)
	53	{
c62b480b	54	//Create the histogram which contain the hough space
3e87ef69	55	Char_t name[1024];
e33f3609	56	sprintf(name,"paramspace_%p",(void*)this);
c62b480b	57	fParamSpace3D = new TH3F(name,"",nxbins,xlow,xup,nybins,ylow,yup,nzbins,zlow,zup);
3e87ef69	58	}
	59
	60	void AliL3HoughTransformerNew::Reset()
	61	{
c62b480b	62	//Reset Hough space
c62b480b	63	fParamSpace3D->Reset();
3e87ef69	64	}
	65
	66
	67
	68	void AliL3HoughTransformerNew::TransformLine(Int_t rowrange,Float_t phirange)
	69	{
c62b480b	70	//Hough Transform
3e87ef69	71	AliL3DigitRowData *tempPt = GetDataPointer();
	72	if(!tempPt)
	73	{
	74	LOG(AliL3Log::kError,"AliL3HoughTransformerNew::TransformLine","Data")
	75	<<"No input data "<<ENDLOG;
	76	return;
	77	}
	78
c62b480b	79	TH3 *hist = fParamSpace3D;
3e87ef69	80	for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
	81	{
	82	AliL3DigitData *digPt = tempPt->fDigitData;
	83	if(i < rowrange[0])
	84	{
	85	AliL3MemHandler::UpdateRowPointer(tempPt);
	86	continue;
	87	}
	88	else if(i > rowrange[1])
	89	break;
	90	if(i != (Int_t)tempPt->fRow)
	91	{
	92	printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
	93	continue;
	94	}
	95	for(UInt_t j=0; j<tempPt->fNDigit; j++)
	96	{
	97	UShort_t charge = digPt[j].fCharge;
	98	UChar_t pad = digPt[j].fPad;
	99	UShort_t time = digPt[j].fTime;
	100	if(charge < GetLowerThreshold() \|\| charge > GetUpperThreshold()) continue;
	101
	102	Int_t sector,row;
	103	Float_t xyz[3];
	104	AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
	105	AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
	106
	107	Float_t phi = AliL3Transform::GetPhi(xyz);
	108	if(phi < phirange[0] \|\| phi > phirange[1])
	109	continue;
	110
	111	xyz[0] = xyz[0] - AliL3Transform::Row2X(rowrange[0]);
	112	Float_t x = xyz[0] + AliL3Transform::Row2X(rowrange[0]);
c62b480b	113	Float_t r = sqrt(xx + xyz[1]xyz[1]);
c62b480b	114	Float_t delta = atan(xyz[2]/r);
3e87ef69	115	for(Int_t xbin=hist->GetXaxis()->GetFirst(); xbin<=hist->GetXaxis()->GetLast(); xbin++)
	116	{
	117	Float_t theta = hist->GetXaxis()->GetBinCenter(xbin);
	118	Float_t rho = xyz[0]cos(theta) + xyz[1]sin(theta);
	119	hist->Fill(theta,rho,delta,charge);
	120	}
	121	}
	122	AliL3MemHandler::UpdateRowPointer(tempPt);
	123	}
	124
	125	}
	126
c62b480b	127	struct AliL3Digit {
	128	Int_t fRow;//padrow index
	129	Int_t fCharge;//digits charge
	130	Float_t fY;//Y coordinate of the digit
	131	Float_t fZ;//Z coordinate of the digit
	132	AliL3Digit *fNext;//pointer to the next digit
3e87ef69	133	};
1f1942b8	134
c62b480b	135	struct AliL3EtaContainer {
	136	void *fFirst;//pointer to the first digit in a sequence
	137	void *fLast;//pointer to the last digit in a sequence
3e87ef69	138	};
1f1942b8	139
3e87ef69	140	void AliL3HoughTransformerNew::TransformLineC(Int_t rowrange,Float_t phirange)
3e87ef69	141	{
c62b480b	142	//Hough Transform
3e87ef69	143	AliL3DigitRowData *tempPt = GetDataPointer();
	144	if(!tempPt)
	145	LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
	146	<<"No input data "<<ENDLOG;
	147
	148	Int_t counter=0;
	149	for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
	150	{
	151	counter += tempPt->fNDigit;
	152	AliL3MemHandler::UpdateRowPointer(tempPt);
	153	}
	154
	155	Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
c62b480b	156	AliL3EtaContainer *etaPt = new AliL3EtaContainer[bound];
c62b480b	157	memset(etaPt,0,bound*sizeof(AliL3EtaContainer));
3e87ef69	158
c62b480b	159	AliL3Digit *digits = new AliL3Digit[counter];
	160	cout<<"Allocating "<<counter*sizeof(AliL3Digit)<<" bytes to digitsarray"<<endl;
	161	memset(digits,0,counter*sizeof(AliL3Digit));
3e87ef69	162
	163	Int_t sector,row;
	164	Float_t xyz[3];
	165
	166	counter=0;
	167	tempPt = GetDataPointer();
	168
	169	cout<<"Calculating digits in patch "<<GetPatch()<<endl;
	170	for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
	171	{
	172	AliL3DigitData *digPt = tempPt->fDigitData;
	173	for(UInt_t di=0; di<tempPt->fNDigit; di++)
	174	{
	175	Int_t charge = digPt[di].fCharge;
	176	Int_t pad = digPt[di].fPad;
	177	Int_t time = digPt[di].fTime;
	178	AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
	179	AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
	180	Double_t eta = AliL3Transform::GetEta(xyz);
	181
	182	Float_t phi = atan2(xyz[1],xyz[0]);
	183	if(phi < phirange[0] \|\| phi > phirange[1]) continue;
	184
c62b480b	185	digits[counter].fRow = i;
	186	digits[counter].fY = xyz[1];
	187	digits[counter].fZ = xyz[2];
	188	digits[counter].fCharge = charge;
3e87ef69	189
c62b480b	190	Int_t etaindex = GetEtaIndex(eta);
c62b480b	191	Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
3e87ef69	192
	193	if(index > 0 && index < bound)
	194	{
c62b480b	195	if(etaPt[index].fFirst == 0)
c62b480b	196	etaPt[index].fFirst = (void*)(&digits[counter]);
3e87ef69	197	else
c62b480b	198	((AliL3Digit*)(etaPt[index].fLast))->fNext = &digits[counter];
c62b480b	199	etaPt[index].fLast = (void*)(&digits[counter]);
3e87ef69	200	}
	201	counter++;
	202	}
	203	AliL3MemHandler::UpdateRowPointer(tempPt);
	204	}
	205
	206	cout<<"Doing the combinatorics"<<endl;
	207
c62b480b	208	AliL3Digit dPt1,dPt2;
c62b480b	209	TH3 *hist = fParamSpace3D;
3e87ef69	210	for(Int_t e=0; e<GetNEtaSegments(); e++)
	211	{
	212	for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
	213	{
	214	Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
	215
c62b480b	216	for(dPt1 = (AliL3Digit)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3Digit)dPt1->fNext)
3e87ef69	217	{
	218	for(Int_t j=i+1; j<=rowrange[1]; j++)
	219	{
	220	Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
	221
c62b480b	222	for(dPt2 = (AliL3Digit)etaPt[index2].fFirst; dPt2 != 0; dPt2 = (AliL3Digit)dPt2->fNext)
3e87ef69	223	{
c62b480b	224	if(dPt1->fRow == dPt2->fRow)
3e87ef69	225	{
	226	cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
	227	exit(5);
	228	}
	229
	230	//Do the transform:
c62b480b	231	Float_t x1 = AliL3Transform::Row2X(dPt1->fRow) - AliL3Transform::Row2X(rowrange[0]);
	232	Float_t x2 = AliL3Transform::Row2X(dPt2->fRow) - AliL3Transform::Row2X(rowrange[0]);
	233	Float_t y1 = dPt1->fY;
	234	Float_t y2 = dPt2->fY;
3e87ef69	235	Float_t theta = atan2(x2-x1,y1-y2);
3e87ef69	236	Float_t rho = x1cos(theta)+y1sin(theta);
c62b480b	237	Float_t r1 = sqrt(pow(AliL3Transform::Row2X(dPt1->fRow),2) + pow(y1,2));
	238	Float_t delta1 = atan2(dPt1->fZ,r1);
	239	Float_t r2 = sqrt(pow(AliL3Transform::Row2X(dPt2->fRow),2) + pow(y2,2));
	240	Float_t delta2 = atan2(dPt2->fZ,r2);
	241	Float_t delta = (dPt1->fChargedelta1+dPt2->fChargedelta2)/(dPt1->fCharge+dPt2->fCharge);
	242	hist->Fill(theta,rho,delta,dPt1->fCharge+dPt2->fCharge);
3e87ef69	243	}
	244	}
	245	}
	246	}
	247	}
	248
	249	cout<<"done"<<endl;
	250	delete [] etaPt;
	251	delete [] digits;
	252
	253	}
	254