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

//$Id$

// Author: Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
//*-- Copyright & copy CL

#include "AliL3StandardIncludes.h"

#include "AliL3Logging.h"
#include "AliL3HoughTransformerLUT.h"
#include "AliL3Transform.h"
#include "AliL3MemHandler.h"
#include "AliL3DigitData.h"
#include "AliL3Histogram.h"

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

//_____________________________________________________________
// AliL3HoughTransformerLUT
//
// Hough transformation class using Look-UP-Tables
//

ClassImp(AliL3HoughTransformerLUT)

AliL3HoughTransformerLUT::AliL3HoughTransformerLUT() : AliL3HoughBaseTransformer()
{
  fMinRow=0;
  fMaxRow=0;

  fNRows=0;
  fNEtas=0;
  fNPhi0=0;
  fSector=0;
  fSlice=0;
  fSectorRow=0;
  fZSign=0;
  fZLengthPlusOff=0;
  fTimeWidth=0;
  fPadPitch=0;
  fEtaSlice=0;

  fLUTX=0;
  fLUTY=0;
  fLUTEta=0;
  fLUTphi0=0;
  fLUT2sinphi0=0;
  fLUT2cosphi0=0;
  fLUTKappa=0;
  
  fLastPad=0;
  fLastIndex=0;
}

AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
{
  fMinRow=0;
  fMaxRow=0;

  fNRows=0;
  fNEtas=0;
  fNPhi0=0;
  fSector=0;
  fSectorRow=0;
  fZSign=0;
  fZLengthPlusOff=0;
  fTimeWidth=0;
  fPadPitch=0;
  fEtaSlice=0;

  fLUTX=0;
  fLUTY=0;
  fLUTEta=0;
  fLUTphi0=0;
  fLUT2sinphi0=0;
  fLUT2cosphi0=0;
  fLUTKappa=0;

  fLastPad=0;
  fLastIndex=0;

  Init(slice,patch,n_eta_segments);
}

AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
{
  DeleteHistograms();

#ifdef do_mc
#if 0
  if(fTrackID)
    {
      for(Int_t i=0; i<fNEtas; i++)
	{
	  if(!fTrackID[i]) continue;
	  delete fTrackID[i];
	}
      delete [] fTrackID;
    }
#endif
#endif

  if(fNRows){
    delete[] fLUTX;
    delete[] fLUTY;
    fNRows=0;
  }

  if(fNEtas){
    delete[] fLUTEta;
    fNEtas=0;
  }
}

void AliL3HoughTransformerLUT::DeleteHistograms()
{
  if(fNPhi0){
    delete[] fLUT2sinphi0;
    delete[] fLUT2cosphi0;
    delete[] fLUTphi0;
    delete[] fLUTKappa;
    fNPhi0=0;
    fLUTphi0=0;
    fLUT2sinphi0=0;
    fLUT2cosphi0=0;
    fLUTKappa=0;
  }

  if(!fParamSpace){
    for(Int_t i=0; i<fNEtas; i++)
      {
	if(!fParamSpace[i]) continue;
	delete fParamSpace[i];
      }
    delete [] fParamSpace;
  }
}

void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t n_eta_segments) 
{
  AliL3HoughBaseTransformer::Init(slice,patch,n_eta_segments);

  //delete old LUT tables
  if(fNRows){
    fNRows=0;
    delete[] fLUTX;
    delete[] fLUTY;
  }
  if(fNEtas){
    delete[] fLUTEta;
    fNEtas=0;
  }

  //member values
  fMinRow=AliL3Transform::GetFirstRow(patch);;
  fMaxRow=AliL3Transform::GetLastRow(patch);;
  fNRows=AliL3Transform::GetNRows(patch);;
  fNEtas=n_eta_segments;
  if(fNEtas!=GetNEtaSegments()) {
    cerr << "AliL3HoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
    exit(1);
  }

  AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
  fZSign = slice < 18 ? 1:-1;
  fSlice = slice;
  fZLengthPlusOff=AliL3Transform::GetZLength()+AliL3Transform::GetZOffset();
  fTimeWidth=AliL3Transform::GetZWidth();

  fPadPitch=0;
  if(fSector<AliL3Transform::GetNSectorLow())
    fPadPitch=AliL3Transform::GetPadPitchWidthLow();
  else
    fPadPitch=AliL3Transform::GetPadPitchWidthUp();  

  Float_t etamax_=GetEtaMax();
  Float_t etamin_=GetEtaMin();
  fEtaSlice=(etamax_-etamin_)/n_eta_segments;

  //lookup tables for X and Y
  fLUTX=new Float_t[fNRows];
  fLUTY=new Float_t[fNRows]; 
  for(Int_t rr=0;rr<fNRows;rr++){
    fLUTX[rr]=Float_t(AliL3Transform::Row2X(rr+fMinRow));
    fLUTY[rr]=Float_t(0.5*(AliL3Transform::GetNPads(rr+fMinRow)-1)*fPadPitch);
    //cout << rr << ": " << (Float_t)fLUTX[rr] << " " << (Float_t)fLUTY[rr] << endl;
  }

  //lookup tables for rz2s <=> etas
  fLUTEta=new Float_t[fNEtas];
  for(Int_t rr=0;rr<fNEtas;rr++){
    Float_t eta=etamin_+(rr+1)*fEtaSlice;
    fLUTEta[rr]=CalcRoverZ2(eta);
    //cout << rr << ": " << eta << " " << fLUTEta[rr] << " " << GetEta(rr,fSlice) << endl;
  }
}

void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t pt_min,Int_t nybin,Double_t phimin,Double_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 AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,Int_t nybin,Double_t ymin,Double_t ymax)
{
  fParamSpace = new AliL3Histogram*[fNEtas];
  
  Char_t histname[256];
  for(Int_t i=0; i<fNEtas; i++)
    {
      sprintf(histname,"paramspace_%d",i);
      fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
    }
  
#ifdef do_mc
  if(fDoMC)
    {
      cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
      exit(1);
#if 0
      AliL3Histogram *hist = fParamSpace[0];
      Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
      cout<<"Allocating "<<fNEtas*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
      fTrackID = new TrackIndex*[fNEtas];
      for(Int_t i=0; i<fNEtas; i++)
	fTrackID[i] = new TrackIndex[ncells];
#endif
    }
#endif

  //create lookup table for sin and cos
  fNPhi0=nybin+1;

  fLUTphi0=new Float_t[fNPhi0];
  fLUT2sinphi0=new Float_t[fNPhi0];
  fLUT2cosphi0=new Float_t[fNPhi0];
  fLUTKappa=new Float_t[fNPhi0];
  Float_t diff=(ymax-ymin)/nybin;
  Float_t phi0=ymin-0.5*diff;
  for(Int_t i=0; i<fNPhi0; i++){
    phi0+=diff;
    fLUTphi0[i]=phi0;
    fLUT2sinphi0[i]=2.*sin(phi0);
    fLUT2cosphi0[i]=2.*cos(phi0);
    fLUTKappa[i]=0.;
    //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
  }  
}

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

  if(!fParamSpace)
    {
      LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
	<<"No histograms to reset"<<ENDLOG;
      return;
    }
  
  for(Int_t i=0; i<fNEtas; i++)
    fParamSpace[i]->Reset();
#ifdef do_mc
  if(fDoMC)
    {
      cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
      exit(1);
#if 0
      AliL3Histogram *hist = fParamSpace[0];
      Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
      for(Int_t i=0; i<fNEtas; i++)
	memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
#endif
    }
#endif
}

Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
{
  //Return the histogram index of the corresponding eta. 
  
  Float_t rz2=CalcRoverZ2(eta);
  return FindIndex(rz2);
}

Int_t AliL3HoughTransformerLUT::FindIndex(Float_t rz2, Int_t start)
{
  //could improve search through devide and conquere strategy
  
  Int_t index=start; 
  if(index==-100){
    index=0;
    while((index<fNEtas)&&(rz2<=fLUTEta[index])){
      index++;
    }
  } else {
    while((index>=0)&&(rz2>fLUTEta[index])){
      index--;
    }
    index++;
  }
  //cout << start << " - " << index << ": " << rz2 << " " << fLUTEta[index] << endl;
  return index;
}

inline AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t eta_index)
{
  if(!fParamSpace || eta_index >= fNEtas || eta_index < 0)
    return 0;
  if(!fParamSpace[eta_index])
    return 0;
  return fParamSpace[eta_index];
}

Float_t AliL3HoughTransformerLUT::CalcRoverZ2(Float_t eta)
{
  Float_t e=exp(2*eta);
  Float_t ret=(e+1)/(e-1);
  ret*=ret;
  return ret;
}

Float_t AliL3HoughTransformerLUT::CalcEta(Float_t roverz2)
{
  Float_t rz=sqrt(roverz2);
  if(fZSign<0) rz=-rz;
  Float_t ret=(1+rz)/(rz-1);
  ret=0.5*log(ret);
  return ret;
}

Float_t AliL3HoughTransformerLUT::CalcX(Int_t row)
{
  return fLUTX[row];
}

Float_t AliL3HoughTransformerLUT::CalcY(Int_t pad,Int_t row)
{
  return pad*fPadPitch-fLUTY[row];
}

Float_t AliL3HoughTransformerLUT::CalcZ(Int_t time)
{
  Float_t ret=time*fTimeWidth;
  if(fZSign>0) ret=fZLengthPlusOff-ret;
  else ret=ret-fZLengthPlusOff;
  return ret;
}

Double_t AliL3HoughTransformerLUT::GetEta(Int_t eta_index,Int_t slice)
{
  if(eta_index >= fNEtas || eta_index < 0){
    //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Index out of range."<<ENDLOG;
    return 0.;
  }
  if(slice != fSlice){
    //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Given slice does not match internal slice."<<ENDLOG;
    return 0.;
  }

  return (CalcEta(fLUTEta[eta_index])-0.5*fEtaSlice);
}

void AliL3HoughTransformerLUT::TransformCircle()
{
  //Transform the input data with a circle HT.
  //The function loops over all the data, and transforms each pixel with the equation:
  // 
  //kappa = 2/R^2 * (y*cos(phi0) - x*sin(phi0) )
  //
  //where R^2 = x^2 + y^2
  //
  //Each pixel 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).

  AliL3DigitRowData *tempPt = GetDataPointer();
  if(!tempPt)
    {
      LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
	<<"No input data "<<ENDLOG;
      return;
    }
  
  //Loop over the padrows:
  for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
    {
      //Get the data on this padrow:
      AliL3DigitData *digPt = tempPt->fDigitData;
      if(i != (Int_t)tempPt->fRow)
	{
	  LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
	    <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
	  continue;
	}

      //calculate x for this row
      Float_t x = CalcX(row);
      Float_t x2=x*x;
      Float_t r2=0;

      //start a new row
      fLastPad=-1;

      //Loop over the data on this padrow:
      for(UInt_t j=0; j<tempPt->fNDigit; j++)
	{
	  UShort_t charge = digPt[j].fCharge;

	  //check threshold
	  if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
	    continue;

	  UChar_t pad = digPt[j].fPad;
	  UShort_t time = digPt[j].fTime;

	  if(fLastPad!=pad){ //only update if necessary
	    fLastIndex=fNEtas-1;   
  
	    //calculate hough for this pad
	    Float_t y = CalcY(pad,row);
	    Float_t y2 = y*y;
	    r2 = x2 + y2;
	    Float_t R2 = 1. / r2;
	    for(Int_t b=0; b<fNPhi0; b++)
	      fLUTKappa[b]=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);

	    fLastPad=pad;
	  }

	  //find eta slice
	  Float_t z = CalcZ(time);
	  Float_t z2=z*z;
	  Float_t rz2 = 1 + r2/z2;
	  Int_t eta_index = FindIndex(rz2,fLastIndex);
	  //cout << row << " " << (int)pad << " " << (int)time << " " << eta_index <<" " <<fLastIndex<< endl;
	  if(eta_index < 0 || eta_index >= fNEtas){
	    //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
	    continue;
	  }	  
	  //Get the correct histogrampointer:
	  AliL3Histogram *hist = fParamSpace[eta_index];
	  if(!hist){
	    //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
	    continue;
	  }
	  fLastIndex=eta_index;

	  //Fill the histogram along the phirange
	  for(Int_t b=0; b<fNPhi0; b++){
	    hist->Fill(fLUTKappa[b],fLUTphi0[b],charge);
	    //hist->Fill(fLUTKappa[b],fLUTphi0[b],1);
	    //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;

#ifdef do_mcc
	    if(fDoMC)
	      {
		cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
		exit(1);
#if 0
		Int_t bin = hist->FindBin(kappa,phi0);
		for(Int_t t=0; t<3; t++)
		  {
		    Int_t label = digPt[j].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) LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","MCData") <<"Array reached maximum!! "<<c<<endl;
		    fTrackID[eta_index][bin].fLabel[c] = label;
		    fTrackID[eta_index][bin].fNHits[c]++;
		  }
#endif
	      }
#endif
	  }
	}
      
      //Move the data pointer to the next padrow:
      AliL3MemHandler::UpdateRowPointer(tempPt);
    }
}

Int_t AliL3HoughTransformerLUT::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
{
  if(!fDoMC)
    {
      cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
      return -1;
    }
  
#ifdef do_mcc
  cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
  exit(1);
#if 0
  if(eta_index < 0 || eta_index > fNEtas)
    {
      cerr<<"AliL3HoughTransformer::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
#endif
  cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
  return -1;
}

void AliL3HoughTransformerLUT::Print()
{
  cout << "fSlice: " << GetSlice() << endl;
  cout << "fPatch: " << GetPatch() << endl;
  cout << "fSector: " << fSector << endl;
  cout << "fSectorRow: " << fSectorRow << endl;
  cout << "fMinRow: " << fMinRow << endl;
  cout << "fMaxRow: " << fMaxRow << endl;
  cout << "fNRows: " << fNRows << endl;
  cout << "fNEtas: " << fNEtas << endl;
  cout << "fNPhi0: " << fNPhi0 << endl;
  cout << "fZSign: " << fZSign << endl;
  cout << "fZLengthPlusOff: " << fZLengthPlusOff << endl;
  cout << "fPadPitch: " << fPadPitch << endl;
  cout << "fTimeWidth: " << fTimeWidth << endl;

  if(!fNRows) return;
  cout << "fLUTX " << fNRows << endl;
  for(Int_t i=0;i<fNRows;i++) cout << "fLUTX[" << i << "]=" << (Float_t)fLUTX[i] << endl;
  cout << "fLUTY " << fNRows << endl;
  for(Int_t i=0;i<fNRows;i++) cout << "fLUTY[" << i << "]=" << fLUTY[i] << endl;
  if(!fNEtas) return;
  cout << "fLUTEta " << fNEtas << endl;
  for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEta[" << i << "]=" << fLUTEta[i] << endl;
  if(!fNPhi0) return;
  cout << "fLUTphi0 " << fNPhi0 << endl;
  for(Int_t i=0;i<fNPhi0;i++) cout << "fLUTPhi0[" << i << "]=" << fLUTphi0[i] << endl;
  cout << "fLUT2sinphi0 " << fNPhi0 << endl;
  for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2sinphi0[" << i << "]=" << fLUT2sinphi0[i] << endl;
  cout << "fLUT2cosphi0 " << fNPhi0 << endl;
  for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2cosphi0[" << i << "]=" << fLUT2cosphi0[i] << endl;
}
Commit	Line	Data
d5d754e5	1	//$Id$
	2
	3	// Author: Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
	4	//*-- Copyright & copy CL
	5
	6	#include "AliL3StandardIncludes.h"
	7
	8	#include "AliL3Logging.h"
	9	#include "AliL3HoughTransformerLUT.h"
	10	#include "AliL3Transform.h"
	11	#include "AliL3MemHandler.h"
	12	#include "AliL3DigitData.h"
	13	#include "AliL3Histogram.h"
	14
	15	#if GCCVERSION == 3
	16	using namespace std;
	17	#endif
	18
	19	//_____________________________________________________________
	20	// AliL3HoughTransformerLUT
	21	//
	22	// Hough transformation class using Look-UP-Tables
	23	//
	24
	25	ClassImp(AliL3HoughTransformerLUT)
	26
	27	AliL3HoughTransformerLUT::AliL3HoughTransformerLUT() : AliL3HoughBaseTransformer()
	28	{
	29	fMinRow=0;
	30	fMaxRow=0;
	31
	32	fNRows=0;
	33	fNEtas=0;
	34	fNPhi0=0;
	35	fSector=0;
3bcf971b	36	fSlice=0;
d5d754e5	37	fSectorRow=0;
	38	fZSign=0;
	39	fZLengthPlusOff=0;
	40	fTimeWidth=0;
	41	fPadPitch=0;
	42	fEtaSlice=0;
	43
	44	fLUTX=0;
	45	fLUTY=0;
	46	fLUTEta=0;
	47	fLUTphi0=0;
	48	fLUT2sinphi0=0;
	49	fLUT2cosphi0=0;
18659a6b	50	fLUTKappa=0;
700db259	51
18659a6b	52	fLastPad=0;
18659a6b	53	fLastIndex=0;
d5d754e5	54	}
	55
	56	AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
	57	{
3bcf971b	58	fMinRow=0;
	59	fMaxRow=0;
	60
	61	fNRows=0;
	62	fNEtas=0;
	63	fNPhi0=0;
	64	fSector=0;
	65	fSectorRow=0;
	66	fZSign=0;
	67	fZLengthPlusOff=0;
	68	fTimeWidth=0;
	69	fPadPitch=0;
	70	fEtaSlice=0;
	71
	72	fLUTX=0;
	73	fLUTY=0;
	74	fLUTEta=0;
	75	fLUTphi0=0;
	76	fLUT2sinphi0=0;
	77	fLUT2cosphi0=0;
18659a6b	78	fLUTKappa=0;
	79
	80	fLastPad=0;
	81	fLastIndex=0;
d5d754e5	82
	83	Init(slice,patch,n_eta_segments);
	84	}
	85
	86	AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
	87	{
	88	DeleteHistograms();
	89
	90	#ifdef do_mc
efb8e06f	91	#if 0
d5d754e5	92	if(fTrackID)
d5d754e5	93	{
3bcf971b	94	for(Int_t i=0; i<fNEtas; i++)
d5d754e5	95	{
	96	if(!fTrackID[i]) continue;
	97	delete fTrackID[i];
	98	}
	99	delete [] fTrackID;
	100	}
efb8e06f	101	#endif
d5d754e5	102	#endif
	103
	104	if(fNRows){
	105	delete[] fLUTX;
	106	delete[] fLUTY;
	107	fNRows=0;
	108	}
3bcf971b	109
d5d754e5	110	if(fNEtas){
	111	delete[] fLUTEta;
	112	fNEtas=0;
	113	}
	114	}
	115
	116	void AliL3HoughTransformerLUT::DeleteHistograms()
	117	{
	118	if(fNPhi0){
	119	delete[] fLUT2sinphi0;
	120	delete[] fLUT2cosphi0;
	121	delete[] fLUTphi0;
18659a6b	122	delete[] fLUTKappa;
d5d754e5	123	fNPhi0=0;
	124	fLUTphi0=0;
	125	fLUT2sinphi0=0;
	126	fLUT2cosphi0=0;
18659a6b	127	fLUTKappa=0;
d5d754e5	128	}
	129
	130	if(!fParamSpace){
3bcf971b	131	for(Int_t i=0; i<fNEtas; i++)
d5d754e5	132	{
	133	if(!fParamSpace[i]) continue;
	134	delete fParamSpace[i];
	135	}
	136	delete [] fParamSpace;
	137	}
	138	}
	139
	140	void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t n_eta_segments)
	141	{
	142	AliL3HoughBaseTransformer::Init(slice,patch,n_eta_segments);
	143
	144	//delete old LUT tables
	145	if(fNRows){
	146	fNRows=0;
	147	delete[] fLUTX;
	148	delete[] fLUTY;
	149	}
	150	if(fNEtas){
	151	delete[] fLUTEta;
	152	fNEtas=0;
	153	}
	154
	155	//member values
	156	fMinRow=AliL3Transform::GetFirstRow(patch);;
	157	fMaxRow=AliL3Transform::GetLastRow(patch);;
	158	fNRows=AliL3Transform::GetNRows(patch);;
	159	fNEtas=n_eta_segments;
3bcf971b	160	if(fNEtas!=GetNEtaSegments()) {
	161	cerr << "AliL3HoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
	162	exit(1);
	163	}
d5d754e5	164
	165	AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
	166	fZSign = slice < 18 ? 1:-1;
3bcf971b	167	fSlice = slice;
d5d754e5	168	fZLengthPlusOff=AliL3Transform::GetZLength()+AliL3Transform::GetZOffset();
	169	fTimeWidth=AliL3Transform::GetZWidth();
	170
	171	fPadPitch=0;
	172	if(fSector<AliL3Transform::GetNSectorLow())
	173	fPadPitch=AliL3Transform::GetPadPitchWidthLow();
	174	else
	175	fPadPitch=AliL3Transform::GetPadPitchWidthUp();
	176
	177	Float_t etamax_=GetEtaMax();
	178	Float_t etamin_=GetEtaMin();
	179	fEtaSlice=(etamax_-etamin_)/n_eta_segments;
	180
	181	//lookup tables for X and Y
	182	fLUTX=new Float_t[fNRows];
	183	fLUTY=new Float_t[fNRows];
	184	for(Int_t rr=0;rr<fNRows;rr++){
	185	fLUTX[rr]=Float_t(AliL3Transform::Row2X(rr+fMinRow));
	186	fLUTY[rr]=Float_t(0.5(AliL3Transform::GetNPads(rr+fMinRow)-1)fPadPitch);
	187	//cout << rr << ": " << (Float_t)fLUTX[rr] << " " << (Float_t)fLUTY[rr] << endl;
	188	}
	189
	190	//lookup tables for rz2s <=> etas
	191	fLUTEta=new Float_t[fNEtas];
	192	for(Int_t rr=0;rr<fNEtas;rr++){
3bcf971b	193	Float_t eta=etamin_+(rr+1)*fEtaSlice;
	194	fLUTEta[rr]=CalcRoverZ2(eta);
	195	//cout << rr << ": " << eta << " " << fLUTEta[rr] << " " << GetEta(rr,fSlice) << endl;
d5d754e5	196	}
	197	}
	198
	199	void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t pt_min,Int_t nybin,Double_t phimin,Double_t phimax)
	200	{
	201	//Create the histograms (parameter space).
	202	//These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
	203	//The arguments give the range and binning;
	204	//nxbin = #bins in kappa
	205	//nybin = #bins in phi0
	206	//pt_min = mimium Pt of track (corresponding to maximum kappa)
	207	//phi_min = mimimum phi0 (degrees)
	208	//phi_max = maximum phi0 (degrees)
	209
	210	Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
	211	Double_t torad = AliL3Transform::Pi()/180;
	212	CreateHistograms(nxbin,-1.x,x,nybin,phimintorad,phimax*torad);
	213	}
	214
	215	void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,Int_t nybin,Double_t ymin,Double_t ymax)
	216	{
	217	fParamSpace = new AliL3Histogram*[fNEtas];
	218
	219	Char_t histname[256];
	220	for(Int_t i=0; i<fNEtas; i++)
	221	{
	222	sprintf(histname,"paramspace_%d",i);
	223	fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
	224	}
	225
	226	#ifdef do_mc
	227	if(fDoMC)
	228	{
efb8e06f	229	cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
	230	exit(1);
	231	#if 0
d5d754e5	232	AliL3Histogram *hist = fParamSpace[0];
d5d754e5	233	Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
3bcf971b	234	cout<<"Allocating "<<fNEtasncellssizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
	235	fTrackID = new TrackIndex*[fNEtas];
	236	for(Int_t i=0; i<fNEtas; i++)
d5d754e5	237	fTrackID[i] = new TrackIndex[ncells];
efb8e06f	238	#endif
d5d754e5	239	}
	240	#endif
	241
	242	//create lookup table for sin and cos
	243	fNPhi0=nybin+1;
	244
	245	fLUTphi0=new Float_t[fNPhi0];
	246	fLUT2sinphi0=new Float_t[fNPhi0];
	247	fLUT2cosphi0=new Float_t[fNPhi0];
18659a6b	248	fLUTKappa=new Float_t[fNPhi0];
d5d754e5	249	Float_t diff=(ymax-ymin)/nybin;
	250	Float_t phi0=ymin-0.5*diff;
	251	for(Int_t i=0; i<fNPhi0; i++){
	252	phi0+=diff;
	253	fLUTphi0[i]=phi0;
3bcf971b	254	fLUT2sinphi0[i]=2.*sin(phi0);
3bcf971b	255	fLUT2cosphi0[i]=2.*cos(phi0);
18659a6b	256	fLUTKappa[i]=0.;
3bcf971b	257	//cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
d5d754e5	258	}
	259	}
	260
	261	void AliL3HoughTransformerLUT::Reset()
	262	{
	263	//Reset all the histograms
	264
	265	if(!fParamSpace)
	266	{
	267	LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
	268	<<"No histograms to reset"<<ENDLOG;
	269	return;
	270	}
	271
3bcf971b	272	for(Int_t i=0; i<fNEtas; i++)
d5d754e5	273	fParamSpace[i]->Reset();
	274	#ifdef do_mc
	275	if(fDoMC)
	276	{
efb8e06f	277	cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
	278	exit(1);
	279	#if 0
d5d754e5	280	AliL3Histogram *hist = fParamSpace[0];
d5d754e5	281	Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
3bcf971b	282	for(Int_t i=0; i<fNEtas; i++)
d5d754e5	283	memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
efb8e06f	284	#endif
d5d754e5	285	}
	286	#endif
	287	}
	288
d5d754e5	289	Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
	290	{
	291	//Return the histogram index of the corresponding eta.
3bcf971b	292
	293	Float_t rz2=CalcRoverZ2(eta);
	294	return FindIndex(rz2);
	295	}
	296
18659a6b	297	Int_t AliL3HoughTransformerLUT::FindIndex(Float_t rz2, Int_t start)
3bcf971b	298	{
18659a6b	299	//could improve search through devide and conquere strategy
	300
	301	Int_t index=start;
	302	if(index==-100){
	303	index=0;
	304	while((index<fNEtas)&&(rz2<=fLUTEta[index])){
	305	index++;
	306	}
	307	} else {
	308	while((index>=0)&&(rz2>fLUTEta[index])){
	309	index--;
	310	}
efb8e06f	311	index++;
3bcf971b	312	}
18659a6b	313	//cout << start << " - " << index << ": " << rz2 << " " << fLUTEta[index] << endl;
3bcf971b	314	return index;
d5d754e5	315	}
	316
	317	inline AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t eta_index)
	318	{
3bcf971b	319	if(!fParamSpace \|\| eta_index >= fNEtas \|\| eta_index < 0)
d5d754e5	320	return 0;
	321	if(!fParamSpace[eta_index])
	322	return 0;
	323	return fParamSpace[eta_index];
	324	}
	325
	326	Float_t AliL3HoughTransformerLUT::CalcRoverZ2(Float_t eta)
	327	{
	328	Float_t e=exp(2*eta);
	329	Float_t ret=(e+1)/(e-1);
	330	ret*=ret;
	331	return ret;
	332	}
	333
	334	Float_t AliL3HoughTransformerLUT::CalcEta(Float_t roverz2)
	335	{
	336	Float_t rz=sqrt(roverz2);
3bcf971b	337	if(fZSign<0) rz=-rz;
d5d754e5	338	Float_t ret=(1+rz)/(rz-1);
	339	ret=0.5*log(ret);
	340	return ret;
	341	}
	342
	343	Float_t AliL3HoughTransformerLUT::CalcX(Int_t row)
	344	{
	345	return fLUTX[row];
	346	}
	347
	348	Float_t AliL3HoughTransformerLUT::CalcY(Int_t pad,Int_t row)
	349	{
	350	return pad*fPadPitch-fLUTY[row];
	351	}
	352
	353	Float_t AliL3HoughTransformerLUT::CalcZ(Int_t time)
	354	{
	355	Float_t ret=time*fTimeWidth;
	356	if(fZSign>0) ret=fZLengthPlusOff-ret;
	357	else ret=ret-fZLengthPlusOff;
	358	return ret;
	359	}
	360
	361	Double_t AliL3HoughTransformerLUT::GetEta(Int_t eta_index,Int_t slice)
	362	{
3bcf971b	363	if(eta_index >= fNEtas \|\| eta_index < 0){
	364	//LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Index out of range."<<ENDLOG;
	365	return 0.;
	366	}
	367	if(slice != fSlice){
	368	//LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Given slice does not match internal slice."<<ENDLOG;
	369	return 0.;
	370	}
	371
	372	return (CalcEta(fLUTEta[eta_index])-0.5*fEtaSlice);
d5d754e5	373	}
	374
	375	void AliL3HoughTransformerLUT::TransformCircle()
	376	{
d5d754e5	377	//Transform the input data with a circle HT.
3bcf971b	378	//The function loops over all the data, and transforms each pixel with the equation:
d5d754e5	379	//
3bcf971b	380	//kappa = 2/R^2 * (ycos(phi0) - xsin(phi0) )
d5d754e5	381	//
3bcf971b	382	//where R^2 = x^2 + y^2
d5d754e5	383	//
	384	//Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
	385	//which histogram in which the pixel should be transformed, the eta-value is calcluated
	386	//and the proper histogram index is found by GetEtaIndex(eta).
	387
d5d754e5	388	AliL3DigitRowData *tempPt = GetDataPointer();
	389	if(!tempPt)
	390	{
3bcf971b	391	LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
d5d754e5	392	<<"No input data "<<ENDLOG;
	393	return;
	394	}
	395
	396	//Loop over the padrows:
3bcf971b	397	for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
d5d754e5	398	{
	399	//Get the data on this padrow:
	400	AliL3DigitData *digPt = tempPt->fDigitData;
	401	if(i != (Int_t)tempPt->fRow)
	402	{
3bcf971b	403	LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
b490510d	404	<<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
d5d754e5	405	continue;
	406	}
	407
efb8e06f	408	//calculate x for this row
18659a6b	409	Float_t x = CalcX(row);
18659a6b	410	Float_t x2=x*x;
18659a6b	411	Float_t r2=0;
efb8e06f	412
efb8e06f	413	//start a new row
18659a6b	414	fLastPad=-1;
18659a6b	415
d5d754e5	416	//Loop over the data on this padrow:
	417	for(UInt_t j=0; j<tempPt->fNDigit; j++)
	418	{
	419	UShort_t charge = digPt[j].fCharge;
3bcf971b	420
3bcf971b	421	//check threshold
d5d754e5	422	if((Int_t)charge <= GetLowerThreshold() \|\| (Int_t)charge > GetUpperThreshold())
d5d754e5	423	continue;
3bcf971b	424
	425	UChar_t pad = digPt[j].fPad;
	426	UShort_t time = digPt[j].fTime;
	427
18659a6b	428	if(fLastPad!=pad){ //only update if necessary
18659a6b	429	fLastIndex=fNEtas-1;
efb8e06f	430
	431	//calculate hough for this pad
	432	Float_t y = CalcY(pad,row);
	433	Float_t y2 = y*y;
18659a6b	434	r2 = x2 + y2;
efb8e06f	435	Float_t R2 = 1. / r2;
18659a6b	436	for(Int_t b=0; b<fNPhi0; b++)
	437	fLUTKappa[b]=R2(yfLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
	438
	439	fLastPad=pad;
	440	}
	441
3bcf971b	442	//find eta slice
efb8e06f	443	Float_t z = CalcZ(time);
	444	Float_t z2=z*z;
	445	Float_t rz2 = 1 + r2/z2;
18659a6b	446	Int_t eta_index = FindIndex(rz2,fLastIndex);
18659a6b	447	//cout << row << " " << (int)pad << " " << (int)time << " " << eta_index <<" " <<fLastIndex<< endl;
3bcf971b	448	if(eta_index < 0 \|\| eta_index >= fNEtas){
3bcf971b	449	//LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
d5d754e5	450	continue;
3bcf971b	451	}
d5d754e5	452	//Get the correct histogrampointer:
d5d754e5	453	AliL3Histogram *hist = fParamSpace[eta_index];
3bcf971b	454	if(!hist){
	455	//LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
	456	continue;
	457	}
efb8e06f	458	fLastIndex=eta_index;
3bcf971b	459
d5d754e5	460	//Fill the histogram along the phirange
18659a6b	461	for(Int_t b=0; b<fNPhi0; b++){
700db259	462	hist->Fill(fLUTKappa[b],fLUTphi0[b],charge);
700db259	463	//hist->Fill(fLUTKappa[b],fLUTphi0[b],1);
3bcf971b	464	//cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
3bcf971b	465
b490510d	466	#ifdef do_mcc
3bcf971b	467	if(fDoMC)
3bcf971b	468	{
efb8e06f	469	cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
	470	exit(1);
	471	#if 0
3bcf971b	472	Int_t bin = hist->FindBin(kappa,phi0);
	473	for(Int_t t=0; t<3; t++)
	474	{
	475	Int_t label = digPt[j].fTrackID[t];
	476	if(label < 0) break;
	477	UInt_t c;
	478	for(c=0; c<MaxTrack; c++)
	479	if(fTrackID[eta_index][bin].fLabel[c] == label \|\| fTrackID[eta_index][bin].fNHits[c] == 0)
	480	break;
	481	if(c == MaxTrack-1) LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","MCData") <<"Array reached maximum!! "<<c<<endl;
	482	fTrackID[eta_index][bin].fLabel[c] = label;
	483	fTrackID[eta_index][bin].fNHits[c]++;
	484	}
efb8e06f	485	#endif
3bcf971b	486	}
d5d754e5	487	#endif
3bcf971b	488	}
d5d754e5	489	}
efb8e06f	490
d5d754e5	491	//Move the data pointer to the next padrow:
	492	AliL3MemHandler::UpdateRowPointer(tempPt);
	493	}
d5d754e5	494	}
d5d754e5	495
868b6166	496	Int_t AliL3HoughTransformerLUT::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
d5d754e5	497	{
	498	if(!fDoMC)
	499	{
	500	cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
	501	return -1;
	502	}
	503
b490510d	504	#ifdef do_mcc
efb8e06f	505	cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
	506	exit(1);
	507	#if 0
3bcf971b	508	if(eta_index < 0 \|\| eta_index > fNEtas)
d5d754e5	509	{
	510	cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
	511	return -1;
	512	}
	513	AliL3Histogram *hist = fParamSpace[eta_index];
	514	Int_t bin = hist->FindBin(kappa,psi);
	515	Int_t label=-1;
	516	Int_t max=0;
	517	for(UInt_t i=0; i<MaxTrack; i++)
	518	{
	519	Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
	520	if(nhits == 0) break;
	521	if(nhits > max)
	522	{
	523	max = nhits;
	524	label = fTrackID[eta_index][bin].fLabel[i];
	525	}
	526	}
	527	return label;
efb8e06f	528	#endif
d5d754e5	529	#endif
	530	cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
	531	return -1;
	532	}
	533
	534	void AliL3HoughTransformerLUT::Print()
	535	{
	536	cout << "fSlice: " << GetSlice() << endl;
	537	cout << "fPatch: " << GetPatch() << endl;
	538	cout << "fSector: " << fSector << endl;
	539	cout << "fSectorRow: " << fSectorRow << endl;
	540	cout << "fMinRow: " << fMinRow << endl;
	541	cout << "fMaxRow: " << fMaxRow << endl;
	542	cout << "fNRows: " << fNRows << endl;
	543	cout << "fNEtas: " << fNEtas << endl;
	544	cout << "fNPhi0: " << fNPhi0 << endl;
	545	cout << "fZSign: " << fZSign << endl;
	546	cout << "fZLengthPlusOff: " << fZLengthPlusOff << endl;
	547	cout << "fPadPitch: " << fPadPitch << endl;
	548	cout << "fTimeWidth: " << fTimeWidth << endl;
	549
	550	if(!fNRows) return;
	551	cout << "fLUTX " << fNRows << endl;
	552	for(Int_t i=0;i<fNRows;i++) cout << "fLUTX[" << i << "]=" << (Float_t)fLUTX[i] << endl;
	553	cout << "fLUTY " << fNRows << endl;
	554	for(Int_t i=0;i<fNRows;i++) cout << "fLUTY[" << i << "]=" << fLUTY[i] << endl;
	555	if(!fNEtas) return;
	556	cout << "fLUTEta " << fNEtas << endl;
	557	for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEta[" << i << "]=" << fLUTEta[i] << endl;
	558	if(!fNPhi0) return;
	559	cout << "fLUTphi0 " << fNPhi0 << endl;
	560	for(Int_t i=0;i<fNPhi0;i++) cout << "fLUTPhi0[" << i << "]=" << fLUTphi0[i] << endl;
	561	cout << "fLUT2sinphi0 " << fNPhi0 << endl;
	562	for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2sinphi0[" << i << "]=" << fLUT2sinphi0[i] << endl;
	563	cout << "fLUT2cosphi0 " << fNPhi0 << endl;
	564	for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2cosphi0[" << i << "]=" << fLUT2cosphi0[i] << endl;
	565	}