-#include <math.h>
-
-#include <TFile.h>
-#include <TTree.h>
-#include <TH2.h>
+// @(#) $Id$
+
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright © ALICE HLT Group
+
+/** \class AliL3HoughTransformer
+<pre>
+//_____________________________________________________________
+// AliL3HoughTransformer
+//
+// Hough transformation class
+//
+</pre>
+*/
-#include "AliTPCParam.h"
-#include "AliSimDigits.h"
+#include "AliL3StandardIncludes.h"
-#include "AliL3Defs.h"
-#include "AliL3Transform.h"
-#include "AliL3HoughPixel.h"
+#include "AliL3Logging.h"
#include "AliL3HoughTransformer.h"
-#include "AliL3HoughTrack.h"
-#include "AliL3TrackArray.h"
+#include "AliL3MemHandler.h"
+#include "AliL3Transform.h"
+#include "AliL3DigitData.h"
+#include "AliL3HistogramAdaptive.h"
+
+#if __GNUC__ >= 3
+using namespace std;
+#endif
ClassImp(AliL3HoughTransformer)
AliL3HoughTransformer::AliL3HoughTransformer()
{
//Default constructor
- fTransform = 0;
- fEtaMin = 0;
- fEtaMax = 0;
- fSlice = 0;
- fPatch = 0;
- fRowContainer = 0;
- fPhiRowContainer = 0;
- fNumOfPadRows=0;
- fContainerBounds=0;
- fNDigits=0;
- fIndex = 0;
+ fParamSpace = 0;
+ fDoMC = kFALSE;;
+ fEtaOverlap=kFALSE;
+#ifdef do_mc
+ fTrackID = 0;
+#endif
+}
+
+AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t netasegments,Bool_t DoEtaOverlap,Bool_t /*DoMC*/) : AliL3HoughBaseTransformer(slice,patch,netasegments)
+{
+ //Normal constructor
+ fParamSpace = 0;
+ fDoMC = kFALSE;
+ fEtaOverlap = DoEtaOverlap;
+ fDoMC=kFALSE;
+#ifdef do_mc
+ fTrackID = 0;
+#endif
}
+AliL3HoughTransformer::~AliL3HoughTransformer()
+{
+ // Dtor
+ DeleteHistograms();
+#ifdef do_mc
+ if(fTrackID)
+ {
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ if(!fTrackID[i]) continue;
+ delete fTrackID[i];
+ }
+ delete [] fTrackID;
+ }
+#endif
+}
-AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments)
+void AliL3HoughTransformer::DeleteHistograms()
{
- //Constructor
+ // Clean up
+ if(!fParamSpace)
+ return;
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ if(!fParamSpace[i]) continue;
+ delete fParamSpace[i];
+ }
+ delete [] fParamSpace;
+ fParamSpace = 0;
+}
+
+void AliL3HoughTransformer::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!!
- fTransform = new AliL3Transform();
+ if(ptmin > ptmax)
+ {
+ cerr<<"AliL3HoughTransformer::CreateHistograms: Error in ptrange "<<ptmin<<" "<<ptmax<<endl;
+ return;
+ }
+ if(psi < 0)
+ {
+ cerr<<"AliL3HoughTransformer::CreateHistograms: Wrong psi-angle "<<psi<<endl;
+ return;
+ }
+ fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
+ Char_t histname[256];
+ Int_t i;
+ for(i=0; i<GetNEtaSegments(); i++)
+ {
+ sprintf(histname,"paramspace_%d",i);
+ fParamSpace[i] = new AliL3HistogramAdaptive(histname,ptmin,ptmax,ptres,nybin,-psi,psi);
+ }
+}
- fEtaMin = etarange[0];
- fEtaMax = etarange[1];
- fSlice = slice;
- fPatch = patch;
- fNPhiSegments = phi_segments;
- fNumOfPadRows=NRowsSlice;
+void AliL3HoughTransformer::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
+ //ptmin = mimium Pt of track (corresponding to maximum kappa)
+ //phimin = mimimum phi0
+ //phimax = maximum phi0
+
+ 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 AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
+ Int_t nybin,Float_t ymin,Float_t ymax)
+{
+ //Create the histograms (parameter space).
+ //nxbin = #bins in X
+ //nybin = #bins in Y
+ //xmin xmax ymin ymax = histogram limits in X and Y
+
+ fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
+
+ Char_t histname[256];
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ {
+ sprintf(histname,"paramspace_%d",i);
+ //fParamSpace[i] = new AliL3HistogramAdaptive(histname,0.5,1.5,0.05,nybin,ymin,ymax);
+ 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<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(AliL3TrackIndex)<<" bytes to fTrackID"<<endl;
+ fTrackID = new AliL3TrackIndex*[GetNEtaSegments()];
+ for(Int_t i=0; i<GetNEtaSegments(); i++)
+ fTrackID[i] = new AliL3TrackIndex[ncells];
+ }
+#endif
+}
-AliL3HoughTransformer::~AliL3HoughTransformer()
+void AliL3HoughTransformer::Reset()
{
- //Destructor
- if(fRowContainer)
- delete [] fRowContainer;
- if(fTransform)
- delete fTransform;
- if(fPhiRowContainer)
- delete [] fPhiRowContainer;
- if(fIndex)
+ //Reset all the histograms
+
+ if(!fParamSpace)
{
- for(Int_t i=0; i<fNDigits; i++)
- delete [] fIndex[i];
- delete [] fIndex;
+ LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
+ <<"No histograms to reset"<<ENDLOG;
+ return;
}
+
+ for(Int_t i=0; i<GetNEtaSegments(); 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<GetNEtaSegments(); i++)
+ memset(fTrackID[i],0,ncells*sizeof(AliL3TrackIndex));
+ }
+#endif
}
-void AliL3HoughTransformer::InitTemplates(TH2F *hist)
+Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta) const
{
+ //Return the histogram index of the corresponding eta.
+
+ Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
+ Double_t index = (eta-GetEtaMin())/etaslice;
+ return (Int_t)index;
+}
- AliL3Digits *pixel;
+void AliL3HoughTransformer::GetEtaIndexes(Double_t eta,Int_t *indexes) const
+{
+ //Return histogram indexes in case of overlapping etaslices.
+
+ Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
+ Int_t index = (Int_t)((eta-GetEtaMin())/etaslice);
+ if(index%2 == 0)
+ {
+ indexes[0] = index;
+ indexes[1] = index - 1;
+ }
+ else
+ {
+ indexes[0] = index - 1;
+ indexes[1] = index;
+ }
+}
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
- Int_t nbinsy = hist->GetNbinsY();
+AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t etaindex)
+{
+ // Return a pointer to the histogram which contains etaindex eta slice
+ if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
+ return 0;
+ if(!fParamSpace[etaindex])
+ return 0;
+ return fParamSpace[etaindex];
+}
- fIndex = new Int_t*[fNDigits];
- for(Int_t i=0; i<fNDigits; i++)
- fIndex[i] = new Int_t[nbinsy+1];
-
- Int_t sector,row;
- for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+Double_t AliL3HoughTransformer::GetEta(Int_t etaindex,Int_t /*slice*/) const
+{
+ // Return eta calculated in the middle of the eta slice
+ Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
+ Double_t eta=0;
+ if(fEtaOverlap)
+ {
+ Int_t index = etaindex + 1;
+ eta=(Double_t)((index)*etaslice);
+ }
+ else
+ eta=(Double_t)((etaindex+0.5)*etaslice);
+ return eta;
+}
+
+void AliL3HoughTransformer::TransformCircle()
+{
+ //Transform the input data with a circle HT.
+ //The function loops over all the data, and transforms each pixel with the equations:
+ //
+ //kappa = 2/r*sin(phi - phi0)
+ //
+ //where r = sqrt(x*x +y*y), and phi = arctan(y/x)
+ //
+ //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 calculated
+ //and the proper histogram index is found by GetEtaIndex(eta).
+
+
+ AliL3DigitRowData *tempPt = GetDataPointer();
+ if(!tempPt)
+ {
+ LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
+ <<"No input data "<<ENDLOG;
+ return;
+ }
+
+ //Loop over the padrows:
+ for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
{
+ //Get the data on this padrow:
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ if(i != (Int_t)tempPt->fRow)
+ {
+ cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
+ continue;
+ }
- for(pixel=(AliL3Digits*)fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
+ //Loop over the data on this padrow:
+ 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((Int_t)charge <= GetLowerThreshold())
+ continue;
+
+ if((Int_t)charge > GetUpperThreshold())
+ charge = GetUpperThreshold();
+
+ Int_t sector,row;
Float_t xyz[3];
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Local(xyz,sector,row,pixel->fPad,pixel->fTime);
+
+ //Transform data to local cartesian coordinates:
+ AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
+ AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
+
+ //Calculate the eta:
+ Double_t eta = AliL3Transform::GetEta(xyz);
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+ //Get the corresponding index, which determines which histogram to fill:
+ Int_t etaindex = GetEtaIndex(eta);
+
+ if(etaindex < 0 || etaindex >= GetNEtaSegments())
+ continue;
- Double_t phi_pix = fTransform->GetPhi(xyz);
- Short_t signal = pixel->fCharge;
- Int_t index = pixel->fIndex;
- if(index >= fNDigits)
- printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index);
- for(Int_t p=ymin; p<=ymax; p++)
+ //Get the correct histogrampointer:
+ AliL3Histogram *hist = fParamSpace[etaindex];
+ if(!hist)
{
-
- Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
- Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
-
- Int_t bin = hist->FindBin(kappa,phi0);
- if(fIndex[index][p]!=0)
- printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n");
- fIndex[index][p] = bin;
+ cerr<<"AliL3HoughTransformer::TransformCircle : Error getting histogram in index "<<etaindex<<endl;
+ continue;
+ }
+
+ //Do the transformation:
+ Float_t r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
+ Float_t phi = AliL3Transform::GetPhi(xyz);
+
+
+ //Fill the histogram along the phirange
+ for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
+ {
+ Float_t phi0 = hist->GetBinCenterY(b);
+ Float_t kappa = 2*sin(phi - phi0)/r;
+ //hist->Fill(kappa,phi0,(int)rint(log((Float_t)charge)));
+ hist->Fill(kappa,phi0,charge);
+ //hist->Fill(kappa,phi0,1);
+#ifdef do_mc
+ 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[etaindex][bin].fLabel[c] == label || fTrackID[etaindex][bin].fNHits[c] == 0)
+ break;
+ if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
+ fTrackID[etaindex][bin].fLabel[c] = label;
+ fTrackID[etaindex][bin].fNHits[c]++;
+ }
+ }
+#endif
}
}
+
+ //Move the data pointer to the next padrow:
+ AliL3MemHandler::UpdateRowPointer(tempPt);
}
-
}
+struct AliL3Digit {
+ Int_t fRow; // Digit padrow
+ Double_t fR; // Digit radius in local coordinate system
+ Double_t fPhi; // Digit Phi angle in local coordinate system
+ Int_t fCharge; // Digit charge
+ AliL3Digit *fNext; // Next digit
+};
-void AliL3HoughTransformer::CountBins()
+struct AliL3EtaContainer {
+ AliL3Digit *fFirst; //First digit
+ AliL3Digit *fLast; //Last digit
+};
+
+void AliL3HoughTransformer::TransformCircleC(Int_t *rowrange,Int_t every)
{
+ //Circle transform, using combinations of every 2 points lying
+ //on different padrows and within the same etaslice.
- Int_t middle_row = 87; //middle of the slice
+ AliL3DigitRowData *tempPt = GetDataPointer();
+ if(!tempPt)
+ LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
+ <<"No input data "<<ENDLOG;
- Double_t r_in_bundle = fTransform->Row2X(middle_row);
- // Double_t phi_min = (fSlice*20 - 10)*ToRad;
- //Double_t phi_max = (fSlice*20 + 10)*ToRad;
- Double_t phi_min = -15*ToRad;
- Double_t phi_max = 15*ToRad;
-
- Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
- Double_t min_phi0 = 10000;
- Double_t max_phi0 = 0;
- Double_t min_kappa = 100000;
- Double_t max_kappa = 0;
+ Int_t minrow = AliL3Transform::GetFirstRow(GetPatch());
+ Int_t maxrow = AliL3Transform::GetLastRow(GetPatch());
+ if(rowrange)
+ {
+ minrow = rowrange[0];
+ maxrow = rowrange[1];
+ if(minrow < AliL3Transform::GetFirstRow(GetPatch()) || minrow >= AliL3Transform::GetLastRow(GetPatch()))
+ minrow = AliL3Transform::GetFirstRow(GetPatch());
+ if(maxrow < AliL3Transform::GetFirstRow(GetPatch()) || maxrow >= AliL3Transform::GetLastRow(GetPatch()))
+ maxrow = AliL3Transform::GetLastRow(GetPatch());
+ if(minrow > maxrow || maxrow==minrow)
+ {
+ cerr<<"AliL3HoughTransformer::TransformCircleC : Bad row range "<<minrow<<" "<<maxrow<<endl;
+ return;
+ }
+ }
+ else
+ {
+ minrow = AliL3Transform::GetFirstRow(GetPatch());
+ maxrow = AliL3Transform::GetLastRow(GetPatch());
+ }
+
+ 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));
- Int_t xbin = 60;
- Int_t ybin = 60;
- Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2->
- Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
+ AliL3Digit *digits = new AliL3Digit[counter];
+ cout<<"Allocating "<<counter*sizeof(AliL3Digit)<<" bytes to digitsarray"<<endl;
+ memset(digits,0,counter*sizeof(AliL3Digit));
+
+ Int_t sector,row,totcharge,pad,time,charge;
+ Double_t r1,r2,phi1,phi2,eta,kappa,phi0;
+ Float_t xyz[3];
- TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+ counter=0;
+ tempPt = GetDataPointer();
- for(Int_t padrow=NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+ cout<<"Calculating digits in patch "<<GetPatch()<<endl;
+ for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
{
- for(Int_t pad=0; pad < fTransform->GetNPads(padrow); pad++)
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ for(UInt_t di=0; di<tempPt->fNDigit; di++)
{
- for(Int_t time = 0; time < fTransform->GetNTimeBins(); time++)
+ charge = digPt[di].fCharge;
+ pad = digPt[di].fPad;
+ time = digPt[di].fTime;
+ AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
+ AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
+ eta = AliL3Transform::GetEta(xyz);
+
+ digits[counter].fRow = i;
+ digits[counter].fR = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
+ digits[counter].fPhi = atan2(xyz[1],xyz[0]);
+ digits[counter].fCharge = charge;
+
+ if(!fEtaOverlap)
{
- Float_t xyz[3];
- Int_t sector,row;
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Global(xyz,sector,row,pad,time);
- Double_t eta = fTransform->GetEta(xyz);
- if(eta < fEtaMin || eta > fEtaMax) continue;
- fTransform->Global2Local(xyz,sector);
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
- Double_t phi_pix = fTransform->GetPhi(xyz);
+ Int_t etaindex = GetEtaIndex(eta);
- for(Int_t p=0; p<=fNPhiSegments; p++)
+ Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
+
+ if(index > 0 && index < bound)
{
- Double_t phi_in_bundle = phi_min + p*phi_slice;
-
- Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
-
- Double_t phi0 = atan(tanPhi0);
- // if(phi0 < 0.55 || phi0 > 0.88) continue;
-
- //if(phi0 < 0) phi0 = phi0 +2*Pi;
- //Double_t kappa = sin(phi_in_bundle - phi0)*2/r_in_bundle;
-
- Double_t angle = phi_pix - phi0;
- Double_t kappa = 2*sin(angle)/r_pix;
- histo->Fill(kappa,phi0,1);
- if(phi0 < min_phi0)
- min_phi0 = phi0;
- if(phi0 > max_phi0)
- max_phi0 = phi0;
- if(kappa < min_kappa)
- min_kappa = kappa;
- if(kappa > max_kappa)
- max_kappa = kappa;
-
+ if(etaPt[index].fFirst == 0)
+ etaPt[index].fFirst = &digits[counter];
+ else
+ (etaPt[index].fLast)->fNext = &digits[counter];
+ etaPt[index].fLast = &digits[counter];
+ }
+ }
+ else
+ {
+ Int_t etaindex[2];
+ GetEtaIndexes(eta,etaindex);
+ Int_t index[2];
+ index[0] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex[0];
+ index[1] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex[1];
+ if(index[0] == index[1])
+ {
+ cerr<<"Same etaindexes "<<index[0]<<" "<<index[1]<<endl;
+ exit(5);
+ }
+
+ Int_t ind = index[0];
+ if(ind > 0 && ind < bound)
+ {
+ if(etaPt[ind].fFirst == 0)
+ etaPt[ind].fFirst = &digits[counter];
+ else
+ (etaPt[ind].fLast)->fNext = &digits[counter];
+ etaPt[ind].fLast = &digits[counter];
}
+ ind = index[1];
+ if(ind > 0 && ind < bound)
+ {
+ if(etaPt[ind].fFirst == 0)
+ etaPt[ind].fFirst = &digits[counter];
+ else
+ (etaPt[ind].fLast)->fNext = &digits[counter];
+ etaPt[ind].fLast = &digits[counter];
+ }
}
-
- }
-
- }
- Int_t count=0,bi=0;
-
- Int_t xmin = histo->GetXaxis()->GetFirst();
- Int_t xmax = histo->GetXaxis()->GetLast();
- Int_t ymin = histo->GetYaxis()->GetFirst();
- Int_t ymax = histo->GetYaxis()->GetLast();
-
- for(Int_t xbin=xmin+1; xbin<xmax; xbin++)
- {
- for(Int_t ybin=ymin+1; ybin<ymax; ybin++)
- {
- bi++;
- Int_t bin = histo->GetBin(xbin,ybin);
- if(histo->GetBinContent(bin)>0)
- count++;
+ counter++;
}
+ AliL3MemHandler::UpdateRowPointer(tempPt);
}
-
-
- printf("Number of possible tracks in this region %d, bins %d\n",count,bi);
- printf("Phi, min %f max %f\n",min_phi0,max_phi0);
- printf("Kappa, min %f max %f\n",min_kappa,max_kappa);
- histo->Draw("box");
-}
-
-
-void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t middle_row)
-{
- //Transformation is done with respect to local coordinates in slice.
-
-
- AliL3Digits *pix1;
- Int_t sector,row;
-
- //Define a common point
-
- /*
- Double_t rowdist1 = fTransform->Row2X(middle_row-1);
- Double_t rowdist2 = fTransform->Row2X(middle_row);
- Double_t r_in_bundle = rowdist1 + (rowdist1-rowdist2)/2;
- */
- Double_t r_in_bundle = fTransform->Row2X(middle_row);
- //Make overlap between slices
- Double_t phi_min = -15*ToRad;
- Double_t phi_max = 15*ToRad;
+ cout<<"Doing the combinatorics"<<endl;
- Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
+ AliL3Digit *dPt1,*dPt2;
- for(Int_t p=0; p <= fNPhiSegments; p++)
+ for(Int_t e=0; e<GetNEtaSegments(); e++)
{
- Double_t phi_in_bundle = phi_min + p*phi_slice;
- //printf("phi %f in slice %d patch %d middle row %f\n",phi_in_bundle/ToRad,fSlice,fPatch,r_in_bundle);
-
- for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
- //for(Int_t padrow = middle_row; padrow <= 173; padrow++)
- //for(Int_t padrow = 0; padrow <= middle_row; padrow++)
+ for(Int_t i=minrow; i<=maxrow; i+=every)
{
+ Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
- for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+ for(dPt1 = (AliL3Digit*)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3Digit*)dPt1->fNext)
{
-
- Float_t xyz[3];
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
- //fTransform->Raw2Global(xyz,sector,row,pix1->fPad,pix1->fTime);
-
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
-
- Double_t phi_pix = fTransform->GetPhi(xyz);
-
- //Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
- Double_t tanPhi0 = (r_in_bundle*sin(phi_pix)-r_pix*sin(phi_in_bundle))/(r_in_bundle*cos(phi_pix)-r_pix*cos(phi_in_bundle));
-
- Double_t phi0 = atan(tanPhi0);
- //if(padrow > middle_row)
- //phi0 = -phi0;
- //if(phi0 < 0.55 || phi0 > 0.88) continue;
-
- Double_t angle = phi_pix - phi0;
- Double_t kappa = 2*sin(angle)/r_pix;
-
- //Double_t angle = phi_in_bundle - phi0;
- //Double_t kappa = 2*sin(angle)/r_in_bundle;
-
- //if(kappa < -0.006 || kappa > 0.006) continue;
-
- Short_t signal = pix1->fCharge;
-
- hist->Fill(kappa,phi0,signal);
-
-
+ for(Int_t j=i+every; j<=maxrow; j+=every)
+ {
+ 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);
+ }
+
+ //Get the correct histogrampointer:
+ AliL3Histogram *hist = fParamSpace[e];
+ if(!hist)
+ {
+ printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",i);
+ continue;
+ }
+
+ //Do the transform:
+ r1 = dPt1->fR;
+ phi1 = dPt1->fPhi;
+ r2 = dPt2->fR;
+ phi2 = dPt2->fPhi;
+ phi0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
+ kappa = 2*sin(phi2-phi0)/r2;
+ totcharge = dPt1->fCharge + dPt2->fCharge;
+ hist->Fill(kappa,phi0,totcharge);
+
+ }
+ }
}
-
}
}
-
-}
+ cout<<"done"<<endl;
+ delete [] etaPt;
+ delete [] digits;
-void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
+}
+
+void AliL3HoughTransformer::TransformLine(Int_t *rowrange,Float_t *phirange)
{
- //Transformation is done with respect to local coordinates in slice.
- //Transform every pixel into whole phirange, using parametrisation:
- //kappa = 2*sin(phi-phi0)/R
+ //Do a line transform on the data.
- printf("Transforming 1 pixel only\n");
- AliL3Digits *pix1;
- Int_t sector,row;
+ AliL3DigitRowData *tempPt = GetDataPointer();
+ if(!tempPt)
+ {
+ LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
+ <<"No input data "<<ENDLOG;
+ return;
+ }
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
- Int_t nbinsy = hist->GetNbinsY();
+ Int_t minrow = AliL3Transform::GetFirstRow(GetPatch());
+ Int_t maxrow = AliL3Transform::GetLastRow(GetPatch());
+ if(rowrange)
+ {
+ minrow = rowrange[0];
+ maxrow = rowrange[1];
+ if(minrow < AliL3Transform::GetFirstRow(GetPatch()) || minrow >= AliL3Transform::GetLastRow(GetPatch()))
+ minrow = AliL3Transform::GetFirstRow(GetPatch());
+ if(maxrow < AliL3Transform::GetFirstRow(GetPatch()) || maxrow >= AliL3Transform::GetLastRow(GetPatch()))
+ maxrow = AliL3Transform::GetLastRow(GetPatch());
+ if(minrow > maxrow || maxrow==minrow)
+ {
+ cerr<<"AliL3HoughTransformer::TransformCircleC : Bad row range "<<minrow<<" "<<maxrow<<endl;
+ return;
+ }
+ }
- for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+ for(Int_t i=minrow; i<=maxrow; i++)
{
-
- for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+ AliL3DigitData *digPt = tempPt->fDigitData;
+ 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())
+ 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);
+
+ if(phirange)
+ {
+ Float_t phi = AliL3Transform::GetPhi(xyz);
+ if(phi < phirange[0] || phi > phirange[1])
+ continue;
+ }
+ Float_t eta = AliL3Transform::GetEta(xyz);
+ Int_t etaindex = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
+ if(etaindex < 0 || etaindex >= GetNEtaSegments())
+ continue;
- Short_t signal = pix1->fCharge;
- Int_t index = pix1->fIndex;
+ xyz[0] = xyz[0] - AliL3Transform::Row2X(minrow);
- for(Int_t p=0; p <= nbinsy; p++)
- hist->AddBinContent(fIndex[index][p],signal);
-
+ //Get the correct histogram:
+ AliL3Histogram *hist = fParamSpace[etaindex];
+ if(!hist)
+ {
+ printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",etaindex);
+ continue;
+ }
+ for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
+ {
+ Double_t theta = hist->GetBinCenterX(xbin);
+ Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
+ hist->Fill(theta,rho,charge);
+ }
}
-
+ AliL3MemHandler::UpdateRowPointer(tempPt);
}
-
}
-/*
- void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
- {
- //Transformation is done with respect to local coordinates in slice.
- //Transform every pixel into whole phirange, using parametrisation:
- //kappa = 2*sin(phi-phi0)/R
-
- printf("Transforming 1 pixel only\n");
-
- AliL3Digits *pix1;
- Int_t sector,row;
-
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
-
- for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
- {
-
- for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
- {
-
- Float_t xyz[3];
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
-
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
-
- Double_t phi_pix = fTransform->GetPhi(xyz);
- Short_t signal = pix1->fCharge;
-
- for(Int_t p=ymin+1; p<=ymax; p++)
- {
- Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
- Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
- hist->Fill(kappa,phi0,signal);
- }
-
- }
-
- }
-
-
- }
-*/
-
-void AliL3HoughTransformer::TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks)
+struct AliL3LDigit {
+ Int_t fRow; // Digit rowpad
+ Int_t fCharge; // Digit charge
+ Float_t fY; // Y position of the digit in the local coor system
+ AliL3LDigit *fNext; // Next digit
+};
+struct AliL3LEtaContainer {
+ AliL3LDigit *fFirst; //First digit
+ AliL3LDigit *fLast; //Last digit
+};
+void AliL3HoughTransformer::TransformLineC(Int_t *rowrange,Float_t *phirange)
{
+ //Circle transform ??
+ AliL3DigitRowData *tempPt = GetDataPointer();
+ if(!tempPt)
+ LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
+ <<"No input data "<<ENDLOG;
- for(Int_t i=0; i<tracks->GetNTracks(); i++)
- {
- AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
- if(!track) {printf("AliL3HoughTransformer::TransformLines2Circle : NO TRACK IN ARRAY\n"); continue;}
-
- Double_t xmin = fTransform->Row2X(track->GetFirstRow());
- Double_t xmax = fTransform->Row2X(track->GetLastRow());
-
- Double_t a = -1./tan(track->GetPsiLine());
- Double_t b = track->GetDLine()/sin(track->GetPsiLine());
-
- Double_t ymin = a*xmin + b;
- Double_t ymax = a*xmax + b;
-
- Double_t middle_x = xmin + (xmax-xmin)/2;
- Double_t middle_y = ymin + (ymax-ymin)/2;
- Double_t r_middle = sqrt(middle_x*middle_x + middle_y*middle_y);
- Double_t phi = atan2(middle_y,middle_x);
- Double_t phi0 = 2*phi - track->GetPsiLine();
- Double_t kappa = 2*sin(phi-phi0)/r_middle;
- hist->Fill(kappa,phi0,track->GetWeight());
-
+ Int_t counter=0;
+ for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
+ {
+ counter += tempPt->fNDigit;
+ AliL3MemHandler::UpdateRowPointer(tempPt);
}
-
-}
-
-void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw)
-{
- //Transform every pixel into histogram, using parametrisation:
- //D = x*cos(psi) + y*sin(psi)
-
- // printf("In Transform; rowrange %d %d ref_row %d phirange %f %f\n",rowrange[0],rowrange[1],ref_row,phirange[0],phirange[1]);
-
- AliL3Digits *pix1;
-
- Int_t xmin = hist->GetXaxis()->GetFirst();
- Int_t xmax = hist->GetXaxis()->GetLast();
+ Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
+ AliL3LEtaContainer *etaPt = new AliL3LEtaContainer[bound];
+ memset(etaPt,0,bound*sizeof(AliL3LEtaContainer));
- Double_t x0 = fTransform->Row2X(ref_row);
- Double_t y0 = 0;
+ AliL3LDigit *digits = new AliL3LDigit[counter];
+ cout<<"Allocating "<<counter*sizeof(AliL3LDigit)<<" bytes to digitsarray"<<endl;
+ memset(digits,0,counter*sizeof(AliL3LDigit));
Int_t sector,row;
- //for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
-
- Double_t phi_min = -10*ToRad;
- Double_t phi_max = 10*ToRad;
- Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
-
-
- Int_t phi_min_index = (Int_t)((phirange[0]*ToRad-phi_min)/delta_phi);
- Int_t phi_max_index = (Int_t)((phirange[1]*ToRad-phi_min)/delta_phi);
-
+ Float_t xyz[3];
- Int_t index;
+ counter=0;
+ tempPt = GetDataPointer();
- for(Int_t phi=phi_min_index; phi <= phi_max_index; phi++)
+ cout<<"Calculating digits in patch "<<GetPatch()<<endl;
+ for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
{
- for(Int_t padrow = rowrange[0]; padrow <= rowrange[1]; padrow++)
+ 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].fCharge = charge;
+
+ Int_t etaindex = GetEtaIndex(eta);
+ Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
- index = phi*fNumOfPadRows + padrow;
- //printf("Looping index %d\n",index);
- if(index > fContainerBounds || index < 0)
+ if(index > 0 && index < bound)
{
- printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index);
- return;
+ if(etaPt[index].fFirst == 0)
+ etaPt[index].fFirst = &digits[counter];
+ else
+ (etaPt[index].fLast)->fNext = &digits[counter];
+ etaPt[index].fLast = &digits[counter];
}
- //for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
- for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel)
+ counter++;
+ }
+ AliL3MemHandler::UpdateRowPointer(tempPt);
+ }
+
+ cout<<"Doing the combinatorics"<<endl;
+
+ AliL3LDigit *dPt1,*dPt2;
+
+ 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 = (AliL3LDigit*)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3LDigit*)dPt1->fNext)
{
- //printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
- Float_t xyz[3];
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
-
- if(raw)
- raw->Fill(xyz[0],xyz[1],pix1->fCharge);
-
- xyz[0] = xyz[0]-x0;
- xyz[1] = xyz[1]-y0;
-
- //printf("calculating...");
- for(Int_t d=xmin+1; d<=xmax; d++)
+ for(Int_t j=i+1; j<=rowrange[1]; j++)
{
- Double_t psi = hist->GetXaxis()->GetBinCenter(d);
- Double_t D = xyz[0]*cos(psi) + xyz[1]*sin(psi);
+ Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
- Short_t signal = pix1->fCharge;
- hist->Fill(psi,D,signal);
- //printf("Filling histo, psi %f D %f\n",psi,D);
+ for(dPt2 = (AliL3LDigit*)etaPt[index2].fFirst; dPt2 != 0; dPt2 = (AliL3LDigit*)dPt2->fNext)
+ {
+ if(dPt1->fRow == dPt2->fRow)
+ {
+ cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
+ exit(5);
+ }
+
+ //Get the correct histogrampointer:
+ AliL3Histogram *hist = fParamSpace[e];
+ if(!hist)
+ {
+ printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",i);
+ continue;
+ }
+
+ //Do the transform:
+ float x1 = AliL3Transform::Row2X(dPt1->fRow) - AliL3Transform::Row2X(rowrange[0]);
+ float x2 = AliL3Transform::Row2X(dPt2->fRow) - AliL3Transform::Row2X(rowrange[0]);
+ float y1 = dPt1->fY;
+ float y2 = dPt2->fY;
+ float theta = atan2(x2-x1,y1-y2);
+ float rho = x1*cos(theta)+y1*sin(theta);
+ hist->Fill(theta,rho,1);//dPt1->charge+dPt2->charge);
+ }
}
- //printf("done\n");
}
- //printf(" done\n");
}
-
}
-
+
+ cout<<"done"<<endl;
+ delete [] etaPt;
+ delete [] digits;
}
-void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
+#ifdef do_mc
+Int_t AliL3HoughTransformer::GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi) const
{
-
- TFile *file = new TFile(rootfile);
- file->cd();
-
- AliTPCParam *param=(AliTPCParam *)file->Get("75x40_100x60");
- TTree *t=(TTree*)file->Get("TreeD_75x40_100x60");
-
- AliSimDigits da, *digarr=&da;
- t->GetBranch("Segment")->SetAddress(&digarr);
- Stat_t num_of_entries=t->GetEntries();
-
- Int_t digit_counter=0;
- Float_t xyz[3];
- Double_t eta,phi;
-
- Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
- printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
-
- if(fRowContainer)
- delete [] fRowContainer;
- fRowContainer = new AliL3HoughContainer[fNumOfPadRows];
- memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer));
-
-
- fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
- printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
- if(fPhiRowContainer)
- delete [] fPhiRowContainer;
- fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
- memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
-
- Double_t phi_min = -10*ToRad;
- Double_t phi_max = 10*ToRad;
- Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
- Int_t index;
- digit_counter=0;
-
- for (Int_t i=0; i<num_of_entries; i++)
- {
- t->GetEvent(i);
- Int_t sector;
- Int_t row;
- param->AdjustSectorRow(digarr->GetID(),sector,row);
- Int_t slice,padrow;
- fTransform->Sector2Slice(slice,padrow,sector,row);
- if(slice != fSlice) continue;
- if(padrow < NRows[fPatch][0]) continue;
- if(padrow > NRows[fPatch][1]) break;
- digarr->First();
- do {
- Int_t time=digarr->CurrentRow();
- Int_t pad=digarr->CurrentColumn();
- Short_t signal=digarr->CurrentDigit();
- if(time < param->GetMaxTBin()-1 && time > 0)
- if(digarr->GetDigit(time-1,pad) <= param->GetZeroSup()
- && digarr->GetDigit(time+1,pad) <= param->GetZeroSup())
- continue;
-
-
- fTransform->Raw2Global(xyz,sector,row,pad,time);
- eta = fTransform->GetEta(xyz);
- if(eta < fEtaMin || eta > fEtaMax) continue;
- fTransform->Global2Local(xyz,sector);
-
-
- phi = fTransform->GetPhi(xyz);
- if(hist)
- hist->Fill(xyz[0],xyz[1],signal);
-
- AliL3Digits *dig = new AliL3Digits;
- dig->fIndex = digit_counter;
- digit_counter++;
- dig->fCharge = signal;
- dig->fPad = pad;
- dig->fTime = time;
-
- if(fRowContainer[padrow].first == NULL)
- fRowContainer[padrow].first = (void*)dig;
- else
- ((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
- fRowContainer[padrow].last = (void*)dig;
-
- Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
- index = phi_index*fNumOfPadRows + padrow;
- if(phi_index > fContainerBounds || phi_index < 0)
- {
- printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
- continue;
- }
-
- if(fPhiRowContainer[index].first == NULL)
- fPhiRowContainer[index].first = (void*)dig;
- else
- ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
- fPhiRowContainer[index].last=(void*)dig;
-
- }while (digarr->Next());
-
+ // Returns the MC label for a given peak found in the Hough space
+ if(!fDoMC)
+ {
+ cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
+ return -1;
}
- fNDigits = digit_counter;
- printf("digitcounter %d\n",digit_counter);
- printf("Allocated %d bytes to pixels\n",digit_counter*sizeof(AliL3Digits));
- file->Close();
- delete file;
+ if(etaindex < 0 || etaindex > GetNEtaSegments())
+ {
+ cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<etaindex<<endl;
+ return -1;
+ }
+ AliL3Histogram *hist = fParamSpace[etaindex];
+ 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[etaindex][bin].fNHits[i];
+ if(nhits == 0) break;
+ if(nhits > max)
+ {
+ max = nhits;
+ label = fTrackID[etaindex][bin].fLabel[i];
+ }
+ }
+ //nhits = max;
+ return label;
+#else
+ Int_t AliL3HoughTransformer::GetTrackID(Int_t /*etaindex*/,Double_t /*kappa*/,Double_t /*psi*/) const
+{
+ // Returns the MC label for a given peak found in the Hough space
+ if(!fDoMC)
+ {
+ cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
+ return -1;
+ }
+ cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
+ return -1;
+#endif
}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff