-//$Id$
+// @(#) $Id$
// Author: Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
-//*-- Copyright & copy CL
+//*-- Copyright © ALICE HLT Group
+/** /class AliL3HoughTransformerLUT
+//<pre>
+//_____________________________________________________________
+// AliL3HoughTransformerLUT
+//
+// Hough transformation class using Look-UP-Tables
+//
+//</pre>
+*/
#include "AliL3StandardIncludes.h"
#include "AliL3MemHandler.h"
#include "AliL3DigitData.h"
#include "AliL3Histogram.h"
+#include "AliL3HistogramAdaptive.h"
-#if GCCVERSION == 3
+#if __GNUC__ >= 3
using namespace std;
#endif
-//_____________________________________________________________
-// AliL3HoughTransformerLUT
-//
-// Hough transformation class using Look-UP-Tables
-//
-
ClassImp(AliL3HoughTransformerLUT)
AliL3HoughTransformerLUT::AliL3HoughTransformerLUT() : AliL3HoughBaseTransformer()
{
+ // default contructor
+ fParamSpace=0;
fMinRow=0;
fMaxRow=0;
fLUTX=0;
fLUTY=0;
fLUTEta=0;
+ fLUTEtaReal=0;
fLUTphi0=0;
fLUT2sinphi0=0;
fLUT2cosphi0=0;
+ fLUTKappa=0;
+
+ fLastPad=0;
+ fLastIndex=0;
+ fAccCharge=0;
+ fX=fY=0.;
}
-AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
+AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t nEtaSegments) : AliL3HoughBaseTransformer(slice,patch,nEtaSegments)
{
+ // constructor
+ fParamSpace=0;
fMinRow=0;
fMaxRow=0;
fLUTX=0;
fLUTY=0;
fLUTEta=0;
+ fLUTEtaReal=0;
fLUTphi0=0;
fLUT2sinphi0=0;
fLUT2cosphi0=0;
+ fLUTKappa=0;
+
+ fLastPad=0;
+ fLastIndex=0;
+ fAccCharge=0;
+ fX=fY=0.;
- Init(slice,patch,n_eta_segments);
+ Init(slice,patch,nEtaSegments);
}
AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
{
+ // destructor
DeleteHistograms();
-#ifdef do_mc
- if(fTrackID)
- {
- for(Int_t i=0; i<fNEtas; i++)
- {
- if(!fTrackID[i]) continue;
- delete fTrackID[i];
- }
- delete [] fTrackID;
- }
-#endif
-
if(fNRows){
delete[] fLUTX;
delete[] fLUTY;
if(fNEtas){
delete[] fLUTEta;
+ delete[] fLUTEtaReal;
fNEtas=0;
}
}
void AliL3HoughTransformerLUT::DeleteHistograms()
{
+ // deletes all histograms
if(fNPhi0){
delete[] fLUT2sinphi0;
delete[] fLUT2cosphi0;
delete[] fLUTphi0;
+ delete[] fLUTKappa;
fNPhi0=0;
fLUTphi0=0;
fLUT2sinphi0=0;
fLUT2cosphi0=0;
+ fLUTKappa=0;
}
if(!fParamSpace){
}
}
-void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t n_eta_segments)
+void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t nEtaSegments,Int_t /*nSeqs*/)
{
- AliL3HoughBaseTransformer::Init(slice,patch,n_eta_segments);
+ // Initialization
+ AliL3HoughBaseTransformer::Init(slice,patch,nEtaSegments);
//delete old LUT tables
if(fNRows){
}
if(fNEtas){
delete[] fLUTEta;
+ delete[] fLUTEtaReal;
fNEtas=0;
}
fMinRow=AliL3Transform::GetFirstRow(patch);;
fMaxRow=AliL3Transform::GetLastRow(patch);;
fNRows=AliL3Transform::GetNRows(patch);;
- fNEtas=n_eta_segments;
+ fNEtas=nEtaSegments;
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;
+ fZSign = slice < 18 ? 1:-1; //see AliL3Transformer
fSlice = slice;
fZLengthPlusOff=AliL3Transform::GetZLength()+AliL3Transform::GetZOffset();
fTimeWidth=AliL3Transform::GetZWidth();
else
fPadPitch=AliL3Transform::GetPadPitchWidthUp();
- Float_t etamax_=GetEtaMax();
- Float_t etamin_=GetEtaMin();
- fEtaSlice=(etamax_-etamin_)/n_eta_segments;
+ Float_t etaMax=GetEtaMax();
+ Float_t etaMin=GetEtaMin();
+ fEtaSlice=(etaMax-etaMin)/nEtaSegments;
//lookup tables for X and Y
fLUTX=new Float_t[fNRows];
}
//lookup tables for rz2s <=> etas
+ //will only be used ifdef FULLLUT
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;
+ Float_t eta=etaMin+(rr+1)*fEtaSlice;
fLUTEta[rr]=CalcRoverZ2(eta);
- //cout << rr << ": " << eta << " " << fLUTEta[rr] << " " << GetEta(rr,fSlice) << endl;
+ 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)
+void AliL3HoughTransformerLUT::CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t ptres,
+ Int_t nybin,Float_t psi)
+{
+ //Create histograms.
+ //_Only_ to be used in case of the adaptive histograms!
+ //phimax is given in radians!!
+
+ if(ptmin > ptmax)
+ {
+ cerr<<"AliL3HoughTransformerLUT::CreateHistograms: Error in ptrange "<<ptmin<<" "<<ptmax<<endl;
+ return;
+ }
+ if(psi < 0)
+ {
+ cerr<<"AliL3HoughTransformerLUT::CreateHistograms: Wrong psi-angle "<<psi<<endl;
+ return;
+ }
+
+ fParamSpace = new AliL3Histogram*[fNEtas];
+ Char_t histname[256];
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ sprintf(histname,"paramspace_%d",i);
+ fParamSpace[i] = new AliL3HistogramAdaptive(histname,ptmin,ptmax,ptres,nybin,-psi,psi);
+ }
+
+ //create lookup table for sin and cos
+ fNPhi0=nybin;
+ fLUTphi0=new Float_t[fNPhi0];
+ fLUT2sinphi0=new Float_t[fNPhi0];
+ fLUT2cosphi0=new Float_t[fNPhi0];
+ fLUTKappa=new Float_t[fNPhi0];
+ AliL3Histogram *hist=fParamSpace[0];
+ Int_t i=0;
+ for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
+ {
+ Float_t phi0 = hist->GetBinCenterY(b);
+ fLUTphi0[i]=phi0;
+ fLUT2sinphi0[i]=2.*sin(phi0);
+ fLUT2cosphi0[i]=2.*cos(phi0);
+ fLUTKappa[i]=0.;
+ i++;
+ //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
+ }
+}
+
+void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Float_t ptMin,Int_t nybin,Float_t phimin,Float_t phimax)
{
//Create the histograms (parameter space).
//These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
//The arguments give the range and binning;
//nxbin = #bins in kappa
//nybin = #bins in phi0
- //pt_min = mimium Pt of track (corresponding to maximum kappa)
+ //ptMin = 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);
+ Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/ptMin;
+ //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)
+void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,Int_t nybin,Float_t ymin,Float_t ymax)
{
fParamSpace = new AliL3Histogram*[fNEtas];
fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
}
-#ifdef do_mc
- if(fDoMC)
- {
- 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
-
//create lookup table for sin and cos
- fNPhi0=nybin+1;
-
+ fNPhi0=nybin;
fLUTphi0=new Float_t[fNPhi0];
fLUT2sinphi0=new Float_t[fNPhi0];
fLUT2cosphi0=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);
- //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
- }
+ fLUTKappa=new Float_t[fNPhi0];
+
+ AliL3Histogram *hist=fParamSpace[0];
+ Int_t i=0;
+ for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
+ {
+ Float_t phi0 = hist->GetBinCenterY(b);
+ 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;
+ i++;
+ }
}
void AliL3HoughTransformerLUT::Reset()
for(Int_t i=0; i<fNEtas; i++)
fParamSpace[i]->Reset();
-#ifdef do_mc
- if(fDoMC)
- {
- 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
}
-Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
+Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta) const
{
//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 index=0; //could improve search through devide and conquere strategy
- while((index<fNEtas)&&(rz2<=fLUTEta[index])){
- index++;
- //cout << index << ": " << rz2 << " " << fLUTEta[index] << endl;
- }
- return index;
-}
-
-inline AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t eta_index)
+AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t etaIndex)
{
- if(!fParamSpace || eta_index >= fNEtas || eta_index < 0)
+ // gets hitogram
+ if(!fParamSpace || etaIndex >= fNEtas || etaIndex < 0)
return 0;
- if(!fParamSpace[eta_index])
+ if(!fParamSpace[etaIndex])
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];
+ return fParamSpace[etaIndex];
}
-Float_t AliL3HoughTransformerLUT::CalcY(Int_t pad,Int_t row)
+Double_t AliL3HoughTransformerLUT::GetEta(Int_t etaIndex,Int_t slice) const
{
- 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;
+ // gets eta
+ if(etaIndex >= fNEtas || etaIndex < 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;
+ 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[etaIndex]);
}
void AliL3HoughTransformerLUT::TransformCircle()
//
//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).
+ //After a day of testing it is proven that this h
AliL3DigitRowData *tempPt = GetDataPointer();
if(!tempPt)
<<"No input data "<<ENDLOG;
return;
}
-
+
+ Int_t lowch=GetLowerThreshold();
+ Int_t highch=GetUpperThreshold();
+
//Loop over the padrows:
for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
{
continue;
}
+
+ //start a new row
+ fLastPad=-1;
+ //store x for this row
+ fX = CalcX(row);
+ //accumulate charge per etaslice
+ fAccCharge=0;
+ //last histogram
+ 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())
+ if((Int_t)charge <= lowch)
continue;
+ if ((Int_t)charge > highch)
+ charge=highch;
UChar_t pad = digPt[j].fPad;
UShort_t time = digPt[j].fTime;
- Float_t x = CalcX(row);
- Float_t y = CalcY(pad,row);
- Float_t z = CalcZ(time);
+ if(fLastPad!=pad){ //only update if necessary
+
+ //if there is accumalated charge,
+ //put it in the old histogram
+ //using the old X and Y values
+ if(fAccCharge>0) AddCurveToHistogram(-1); //fLastIndex will be set to -1
+#ifdef FULLLUT
+ fLastIndex=fNEtas-1;
+#endif
+
+ //calculate Y for the new pad
+ fY = CalcY(pad,row);
+ //remember new pad
+ fLastPad=pad;
+ }
+
+#ifdef FULLLUT
//find eta slice
- Float_t rz2 = 1 + (x*x+y*y)/(z*z);
- Int_t eta_index = FindIndex(rz2);
+ Float_t z = CalcZ(time);
- if(eta_index < 0 || eta_index >= fNEtas){
- //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
+ Float_t z2=z*z;
+ Float_t rz2 = 1 + r2/z2;
+ Int_t etaIndex = FindIndex(rz2,fLastIndex);
+#else
+ Float_t z = CalcZ(time);
+ Double_t r = sqrt(fX*fX+fY*fY+z*z);
+ Double_t eta = 0.5 * log((r+z)/(r-z));
+ Int_t etaIndex = GetEtaIndex(eta);
+#endif
+ if(etaIndex < 0 || etaIndex >= fNEtas){
+ LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<etaIndex<<"."<<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;
- }
-
- //Fill the histogram along the phirange
- Float_t R2=1/(x*x+y*y);
- for(Int_t b=0; b<fNPhi0; b++){
- Float_t kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
-
- hist->Fill(kappa,fLUTphi0[b],charge);
- //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
-
-#ifdef do_mcc
- if(fDoMC)
- {
- 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]++;
- }
- }
+#ifndef FULLLUT
+ if(fLastIndex==-1) fLastIndex=etaIndex;
#endif
+
+ if(fLastIndex!=etaIndex){ //enter old values first
+ AddCurveToHistogram(etaIndex);
}
+
+ 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);
}
}
-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
- 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
- cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
- return -1;
-}
-
void AliL3HoughTransformerLUT::Print()
{
+ // debug printout
cout << "fSlice: " << GetSlice() << endl;
cout << "fPatch: " << GetPatch() << endl;
cout << "fSector: " << fSector << 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;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff