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

// @(#) $Id$

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

#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

//#define FULLLUT

//_____________________________________________________________
// 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;
  fLUTEtaReal=0;
  fLUTphi0=0;
  fLUT2sinphi0=0;
  fLUT2cosphi0=0;
  fLUTKappa=0;
  
  fLastPad=0;
  fLastIndex=0;
  fAccCharge=0;

  fDoMC = kFALSE;
}

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;
  fLUTEtaReal=0;
  fLUTphi0=0;
  fLUT2sinphi0=0;
  fLUT2cosphi0=0;
  fLUTKappa=0;

  fLastPad=0;
  fLastIndex=0;
  fAccCharge=0;
  fDoMC=kFALSE;

  Init(slice,patch,n_eta_segments);
}

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

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

  if(fNEtas){
    delete[] fLUTEta;
    delete[] fLUTEtaReal;
    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,Int_t n_seqs) 
{
  AliL3HoughBaseTransformer::Init(slice,patch,n_eta_segments);

  //delete old LUT tables
  if(fNRows){
    fNRows=0;
    delete[] fLUTX;
    delete[] fLUTY;
  }
  if(fNEtas){
    delete[] fLUTEta;
    delete[] fLUTEtaReal;
    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];
  fLUTEtaReal=new Float_t[fNEtas];
  for(Int_t rr=0;rr<fNEtas;rr++){
    Float_t eta=etamin_+(rr+1)*fEtaSlice;
    fLUTEta[rr]=CalcRoverZ2(eta);
    fLUTEtaReal[rr]=eta-0.5*fEtaSlice;
    //cout << rr << ": " << eta << " " << fLUTEtaReal[rr] << " " << GetEta(rr,fSlice) << " - " << fLUTEta[rr] << 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);
    }
  
  //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();
}

Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
{
  //Return the histogram index of the corresponding eta. 
  
#ifdef FULLLUT 
  /* try to imitate a circuit -> should 
     go into the VHDL implementation of transformer */
  Float_t rz2=CalcRoverZ2(eta);
  return FindIndex(rz2);
#else /* optimize for speed on the computer */
  Double_t index = (eta-GetEtaMin())/fEtaSlice;
  return (Int_t)index;
#endif
}

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);
  return(fLUTEtaReal[eta_index]);
}

void AliL3HoughTransformerLUT::AddCurveToHistogram(Int_t new_eta_index)
{
  //Get the correct histogrampointer:
  AliL3Histogram *hist = fParamSpace[fLastIndex];
  if(!hist){
    //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"Error getting histogram in index "<<fLastIndex<<"."<<ENDLOG;
    return;
  }

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

  fAccCharge=0;
  fLastIndex=new_eta_index;
}

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;

      if(fAccCharge>0) cerr << "Big error " << endl;

      //accumulate charge per histogram
      fAccCharge=0;
      fLastIndex=-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

	    //if there is accumalated charge, put it in the right histogram
	    if(fAccCharge>0) AddCurveToHistogram(-1);

#ifdef FULLLUT
	    fLastIndex=fNEtas-1;
#endif

	    //calculate hough for the new pad
	    Float_t y = CalcY(pad,row);
	    Float_t y2 = y*y;
	    r2 = x2 + y2;
	    Float_t xr2=x/r2;
	    Float_t yr2=y/r2;

	    //AliL3Histogram *hist = fParamSpace[0];
	    //Int_t bb=hist->GetFirstYbin();
	    //if(fNPhi0-1!=hist->GetLastYbin()) cerr << "should not be" << endl;
	    for(Int_t b=0; b<fNPhi0; b++){
	      fLUTKappa[b]=(yr2*fLUT2cosphi0[b]-xr2*fLUT2sinphi0[b]);

	      //Float_t phi=atan2(y,x);
	      //Float_t phi0 = hist->GetBinCenterY(bb);
	      //Float_t R=sqrt(r2);
	      //Float_t kappa = 2*sin(phi - phi0)/R;
	      //cout << fLUTKappa[b] << " " << kappa << endl;
	      //bb++;
	      //Int_t bina=hist->FindBin(fLUTphi0[b],fLUTKappa[b]);
	      //Int_t binb=hist->FindBin(phi0,kappa);
	      //if(bina!=binb) cout << bina << " " << binb << endl;
	    }
	    
	    fLastPad=pad;
	  }

	  //find eta slice
	  Float_t z = CalcZ(time);
	  Float_t z2=z*z;

#ifdef FULLLUT
	  Float_t rz2 = 1 + r2/z2;
	  Int_t eta_index = FindIndex(rz2,fLastIndex);
#else
	  Double_t r = sqrt(r2+z2);
	  Double_t etaval = 0.5 * log((r+z)/(r-z));
	  Int_t eta_index = GetEtaIndex(etaval);
#endif
	  //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;
	  }	  

#ifndef FULLLUT
	  if(fLastIndex==-1) fLastIndex=eta_index;
#endif

	  if(fLastIndex!=eta_index){
	    AddCurveToHistogram(eta_index);
	  }

	  fAccCharge+=charge;
	  //fAccCharge+=1;
	}

      //in case there is charge, store it before moving to the next row
      if(fAccCharge>0) AddCurveToHistogram(-1);
      
      //Move the data pointer to the next padrow:
      AliL3MemHandler::UpdateRowPointer(tempPt);
    }

      if(fAccCharge>0) cerr << "Big error " << endl;
}

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