AliL3Histogram::AliL3Histogram()
{
+ // Default constructor
fNxbins = 0;
fNybins = 0;
fNcells = 0;
Int_t nxbin,Double_t xmin,Double_t xmax,
Int_t nybin,Double_t ymin,Double_t ymax)
{
+ // Normal constructor
strcpy(fName,name);
fNxbins = nxbin;
void AliL3Histogram::Reset()
{
+ // Reset histogram contents
if(fContent)
for(Int_t i=0; i<fNcells; i++) fContent[i] = 0;
void AliL3Histogram::Fill(Double_t x,Double_t y,Int_t weight)
{
+ // Fill the weight into a bin which correspond to x and y
Int_t bin = FindBin(x,y);
#ifdef _IFON_
if(bin < 0)
void AliL3Histogram::Fill(Double_t x,Int_t ybin,Int_t weight)
{
+ // Fill the weight into a bin which correspond to x and ybin
Int_t xbin = FindXbin(x);
Int_t bin = GetBin(xbin,ybin);
#ifdef _IFON_
void AliL3Histogram::Fill(Int_t xbin,Double_t y,Int_t weight)
{
+ // Fill the weight into a bin which correspond to xbin and y
Int_t ybin = FindYbin(y);
Int_t bin = GetBin(xbin,ybin);
#ifdef _IFON_
void AliL3Histogram::Fill(Int_t xbin,Int_t ybin,Int_t weight)
{
+ // Fill the weight into a bin which correspond to xbin and ybin
Int_t bin = GetBin(xbin,ybin);
#ifdef _IFON_
if(bin < 0)
Int_t AliL3Histogram::FindBin(Double_t x,Double_t y) const
{
+ // Finds the bin which correspond to x and y
Int_t xbin = FindXbin(x);
Int_t ybin = FindYbin(y);
#ifdef _IFON_
Int_t AliL3Histogram::FindLabelBin(Double_t x,Double_t y) const
{
+ // Returns the corresponding bin with the mc labels
Int_t xbin = FindXbin(x);
Int_t ybin = FindYbin(y);
#ifdef _IFON_
Int_t AliL3Histogram::FindXbin(Double_t x) const
{
+ // Finds the bin which correspond to x
if(x < fXmin || x > fXmax)
return 0;
Int_t AliL3Histogram::FindYbin(Double_t y) const
{
+ // Finds the bin which correspond to y
if(y < fYmin || y > fYmax)
return 0;
Int_t AliL3Histogram::GetBin(Int_t xbin,Int_t ybin) const
{
+ // Returns the bin which correspond to xbin and ybin
if(xbin < fFirstXbin || xbin > fLastXbin)
return 0;
if(ybin < fFirstYbin || ybin > fLastYbin)
Int_t AliL3Histogram::GetLabelBin(Int_t xbin,Int_t ybin) const
{
+ // Returns the corresponding bin with the mc labels
if(xbin < fFirstXbin || xbin > fLastXbin)
return -1;
if(ybin < fFirstYbin || ybin > fLastYbin)
Int_t AliL3Histogram::GetBinContent(Int_t bin) const
{
+ // Return the bin content
if(bin >= fNcells)
{
LOG(AliL3Log::kError,"AliL3Histogram::GetBinContent","array")<<AliL3Log::kDec<<
void AliL3Histogram::SetBinContent(Int_t xbin,Int_t ybin,Int_t value)
{
+ // Set bin content
Int_t bin = GetBin(xbin,ybin);
#ifdef _IFON_
if(bin == 0)
void AliL3Histogram::SetBinContent(Int_t bin,Int_t value)
{
+ // Set bin content
if(bin >= fNcells)
{
void AliL3Histogram::AddBinContent(Int_t xbin,Int_t ybin,Int_t weight)
{
+ // Adds weight to bin content
Int_t bin = GetBin(xbin,ybin);
#ifdef _IFON_
if(bin == 0)
void AliL3Histogram::AddBinContent(Int_t bin,Int_t weight)
{
+ // Adds weight to bin content
if(bin < 0 || bin > fNcells)
{
LOG(AliL3Log::kError,"AliL3Histogram::AddBinContent","array")<<AliL3Log::kDec<<
Double_t AliL3Histogram::GetBinCenterX(Int_t xbin) const
{
+ // Returns the position of the center of a bin
if(xbin < fFirstXbin || xbin > fLastXbin)
{
LOG(AliL3Log::kError,"AliL3Histogram::GetBinCenterX","xbin")
Double_t AliL3Histogram::GetBinCenterY(Int_t ybin) const
{
+ // Returns the position of the center of a bin
if(ybin < fFirstYbin || ybin > fLastYbin)
{
LOG(AliL3Log::kError,"AliL3Histogram::GetBinCenterY","ybin")
Double_t AliL3Histogram::GetPreciseBinCenterX(Float_t xbin) const
{
+ // Returns the position of the center of a bin using precise values inside the bin
if(xbin < (fFirstXbin-0.5) || xbin > (fLastXbin+0.5))
{
LOG(AliL3Log::kError,"AliL3Histogram::GetBinCenterX","xbin")
Double_t AliL3Histogram::GetPreciseBinCenterY(Float_t ybin) const
{
+ // Returns the position of the center of a bin using precise values inside the bin
if(ybin < (fFirstYbin-0.5) || ybin > (fLastYbin+0.5))
{
LOG(AliL3Log::kError,"AliL3Histogram::GetBinCenterY","ybin")
void AliL3Histogram::Draw(Char_t *option)
{
+ // Fill the contents of the corresponding ROOT histogram and draws it
#ifdef use_root
if(!fRootHisto)
CreateRootHisto();
void AliL3Histogram::CreateRootHisto()
{
+ // Create ROOT histogram out of AliL3Histogram
#ifdef use_root
fRootHisto = new TH2F(fName,"",fNxbins,fXmin,fXmax,fNybins,fYmin,fYmax);
return;
// @(#) $Id$
-#ifndef ALIL3_HISTOGRAM
-#define ALIL3_HISTOGRAM
+#ifndef ALIL3HISTOGRAM_H
+#define ALIL3HISTOGRAM_H
#include "AliL3StandardIncludes.h"
#include "AliL3RootTypes.h"
#endif
class AliL3Histogram {
-
- private:
- Double_t fBinwidthX;
- Double_t fBinwidthY;
-
- protected:
- Int_t *fContent; //!
- Char_t fName[100];
- Int_t fNxbins;
- Int_t fNybins;
- Int_t fNcells;
- Int_t fEntries;
- Int_t fFirstXbin;
- Int_t fFirstYbin;
- Int_t fLastXbin;
- Int_t fLastYbin;
- Int_t fThreshold;
-
- Double_t fXmin;
- Double_t fYmin;
- Double_t fXmax;
- Double_t fYmax;
-
-#ifdef use_root
- TH2F *fRootHisto;
-#endif
-
public:
AliL3Histogram();
Int_t GetNbinsX() const {return fNxbins;}
Int_t GetNbinsY() const {return fNybins;}
Int_t GetNEntries() const {return fEntries;}
+
+ protected:
+ Int_t *fContent; //!
+ Char_t fName[100]; // Name of the histogram
+ Int_t fNxbins; // Number of bins in the histogram
+ Int_t fNybins; // Number of bins in the histogram
+ Int_t fNcells; // Overall number of bins in the histogram
+ Int_t fEntries; // Number of entries in the histogram
+ Int_t fFirstXbin; // First active bin
+ Int_t fFirstYbin; // First active bin
+ Int_t fLastXbin; // Last active bin
+ Int_t fLastYbin; // Last active bin
+ Int_t fThreshold; // Bin content threshold
+
+ Double_t fXmin; // Lower limit in X
+ Double_t fYmin; // Lower limit in Y
+ Double_t fXmax; // Upper limit in X
+ Double_t fYmax; // Upper limit in Y
+
+#ifdef use_root
+ TH2F *fRootHisto; // Corresponding ROOT histogram
+#endif
+
+ private:
+ Double_t fBinwidthX; // Bin width of the Hough space
+ Double_t fBinwidthY; // Bin width of the Hough space
ClassDef(AliL3Histogram,1) //2D histogram class
#endif
}
-AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr)
+AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t netasegments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr)
{
+ //Normal constructor
fBinary = binary;
strcpy(fPath,path);
- fNEtaSegments = n_eta_segments;
+ fNEtaSegments = netasegments;
fAddHistograms = kFALSE;
fDoIterative = kFALSE;
fWriteDigits = kFALSE;
//cout << "Cleaned class mem " << endl;
}
-void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr,Float_t zvertex)
+void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t netasegments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr,Float_t zvertex)
{
+ //Normal init of the AliL3Hough
fBinary = binary;
strcpy(fPath,path);
- fNEtaSegments = n_eta_segments;
+ fNEtaSegments = netasegments;
fWriteDigits = kFALSE;
fUse8bits = bit8;
fVersion = tv;
void AliL3Hough::Init(Bool_t doit, Bool_t addhists)
{
+ // Init
fDoIterative = doit;
fAddHistograms = addhists;
void AliL3Hough::SetTransformerParams(Float_t ptres,Float_t ptmin,Float_t ptmax,Int_t ny,Int_t patch)
{
-
+ // Setup the parameters for the Hough Transformer
Int_t mrow;
Float_t psi=0;
if(patch==-1)
/*
void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t patch)
{
+ // Setup the parameters for the Hough Transformer
Int_t mrow=80;
Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
*/
void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t /*patch*/)
{
+ // Setup the parameters for the Hough Transformer
Int_t mrow=79;
void AliL3Hough::SetThreshold(Int_t t3,Int_t patch)
{
+ // Set digits threshold
if(patch==-1)
{
Int_t i=0;
void AliL3Hough::SetPeakThreshold(Int_t threshold,Int_t patch)
{
+ // Set Peak Finder threshold
if(patch==-1)
{
Int_t i=0;
fEvent=eventnr;
}
-void AliL3Hough::Transform(Int_t *row_range)
+void AliL3Hough::Transform(Int_t *rowrange)
{
//Transform all data given to the transformer within the given slice
//(after ReadData(slice))
if((fVersion != 4) || (i == 0))
fHoughTransformer[i]->Reset();//Reset the histograms
fBenchmark->Start("Hough Transform");
- if(!row_range)
+ if(!rowrange)
fHoughTransformer[i]->TransformCircle();
else
- fHoughTransformer[i]->TransformCircleC(row_range,1);
+ fHoughTransformer[i]->TransformCircleC(rowrange,1);
fBenchmark->Stop("Hough Transform");
}
cpuTime = GetCpuTime() - initTime;
void AliL3Hough::MergePatches()
{
+ // Merge patches if they are not summed
if(fAddHistograms) //Nothing to merge here
return;
fMerger->MergePatches(kTRUE);
void AliL3Hough::MergeInternally()
{
+ // Merge patches internally
if(fAddHistograms)
fInterMerger->FillTracks(fTracks[0]);
else
//Process patch in a iterative way.
//transform + peakfinding + evaluation + transform +...
- Int_t num_of_tries = 5;
+ Int_t numoftries = 5;
AliL3HoughBaseTransformer *tr = fHoughTransformer[patch];
AliL3TrackArray *tracks = fTracks[patch];
tracks->Reset();
ev->SetNumOfRowsToMiss(3);
ev->SetNumOfPadsToLook(2);
AliL3Histogram *hist;
- for(Int_t t=0; t<num_of_tries; t++)
+ for(Int_t t=0; t<numoftries; t++)
{
tr->Reset();
tr->TransformCircle();
void AliL3Hough::AddTracks()
{
+ // Add current slice slice tracks to the global list of found tracks
if(!fTracks[0])
{
cerr<<"AliL3Hough::AddTracks : No tracks"<<endl;
void AliL3Hough::FindTrackCandidatesRow()
{
+ // Find AliL3HoughTransformerRow track candidates
if(fVersion != 4) {
LOG(AliL3Log::kError,"AliL3Hough::FindTrackCandidatesRow()","")
<<"Incompatible Peak Finder version!"<<ENDLOG;
}
//Look for peaks in histograms, and find the track candidates
- Int_t n_patches;
+ Int_t npatches;
if(fAddHistograms)
- n_patches = 1; //Histograms have been added.
+ npatches = 1; //Histograms have been added.
else
- n_patches = fNPatches;
+ npatches = fNPatches;
Double_t initTime,cpuTime;
initTime = GetCpuTime();
fBenchmark->Start("Find Maxima");
- for(Int_t i=0; i<n_patches; i++)
+ for(Int_t i=0; i<npatches; i++)
{
AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
fTracks[i]->Reset();
void AliL3Hough::FindTrackCandidates()
{
+ // Find AliL3HoughTransformer track candidates
if(fVersion == 4) {
LOG(AliL3Log::kError,"AliL3Hough::FindTrackCandidatesRow()","")
<<"Incompatible Peak Finder version!"<<ENDLOG;
return;
}
- Int_t n_patches;
+ Int_t npatches;
if(fAddHistograms)
- n_patches = 1; //Histograms have been added.
+ npatches = 1; //Histograms have been added.
else
- n_patches = fNPatches;
+ npatches = fNPatches;
Double_t initTime,cpuTime;
initTime = GetCpuTime();
fBenchmark->Start("Find Maxima");
- for(Int_t i=0; i<n_patches; i++)
+ for(Int_t i=0; i<npatches; i++)
{
AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
fTracks[i]->Reset();
fEval[i]->InitTransformer(fHoughTransformer[i]);
}
-Int_t AliL3Hough::Evaluate(Int_t road_width,Int_t nrowstomiss)
+Int_t AliL3Hough::Evaluate(Int_t roadwidth,Int_t nrowstomiss)
{
//Evaluate the tracks, by looking along the road in the raw data.
//If track does not cross all padrows - rows2miss, it is removed from the arrray.
return 0;
}
- Int_t removed_tracks=0;
+ Int_t removedtracks=0;
AliL3TrackArray *tracks=0;
if(fAddHistograms)
}
for(Int_t i=0; i<fNPatches; i++)
- EvaluatePatch(i,road_width,nrowstomiss);
+ EvaluatePatch(i,roadwidth,nrowstomiss);
//Here we check the tracks globally;
//how many good rows (padrows with signal)
if(track->GetNHits() < AliL3Transform::GetNRows() - nrowstomiss)
{
tracks->Remove(j);
- removed_tracks++;
+ removedtracks++;
}
}
tracks->Compress();
tracks->QSort();
}
- return removed_tracks;
+ return removedtracks;
}
-void AliL3Hough::EvaluatePatch(Int_t i,Int_t road_width,Int_t nrowstomiss)
+void AliL3Hough::EvaluatePatch(Int_t i,Int_t roadwidth,Int_t nrowstomiss)
{
//Evaluate patch i.
fEval[i]->InitTransformer(fHoughTransformer[i]);
- fEval[i]->SetNumOfPadsToLook(road_width);
+ fEval[i]->SetNumOfPadsToLook(roadwidth);
fEval[i]->SetNumOfRowsToMiss(nrowstomiss);
//fEval[i]->RemoveFoundTracks();
void AliL3Hough::WriteTracks(Char_t *path)
{
+ // Write found tracks into file
//cout<<"AliL3Hough::WriteTracks : Sorting the tracsk"<<endl;
//fGlobalTracks->QSort();
void AliL3Hough::WriteTracks(Int_t slice,Char_t *path)
{
+ // Write found tracks slice by slice into file
AliL3MemHandler mem;
Char_t fname[100];
void AliL3Hough::WriteDigits(Char_t *outfile)
{
-#ifdef use_aliroot
//Write the current data to a new rootfile.
+#ifdef use_aliroot
for(Int_t i=0; i<fNPatches; i++)
{
// @(#) $Id$
-#ifndef ALIL3HOUGH
-#define ALIL3HOUGH
+#ifndef ALIL3HOUGH_H
+#define ALIL3HOUGH_H
#include "AliL3RootTypes.h"
#endif
class AliL3Hough {
-
- private:
- Char_t *fInputFile;//!
- Char_t *fInputPtr;//!
- Char_t fPath[1024];
- Bool_t fBinary;
- Bool_t fAddHistograms;
- Bool_t fDoIterative;
- Bool_t fWriteDigits;
- Bool_t fUse8bits;
- Int_t fNEtaSegments;
- Int_t fNPatches;
- Int_t fVersion; //which HoughTransformer to use
- Int_t fCurrentSlice;
- Int_t fEvent;
-
- Int_t fPeakThreshold[6];
- Float_t fLowPt[6];
- Float_t fUpperPt[6];
- Float_t fPtRes[6];
- Float_t fPhi[6];
- Int_t fNBinX[6];
- Int_t fNBinY[6];
- Int_t fThreshold[6];
- Int_t fNSaveIterations; //for HoughtransformerVhdl
-
- //parameters for the peak finder:
- Int_t fKappaSpread;
- Float_t fPeakRatio;
-
- Float_t fZVertex;
-
- AliL3MemHandler **fMemHandler; //!
- AliL3HoughBaseTransformer **fHoughTransformer; //!
- AliL3HoughEval **fEval; //!
- AliL3HoughMaxFinder *fPeakFinder; //!
- AliL3TrackArray **fTracks; //!
- AliL3TrackArray *fGlobalTracks; //!
- AliL3HoughMerger *fMerger; //!
- AliL3HoughIntMerger *fInterMerger; //!
- AliL3HoughGlobalMerger *fGlobalMerger; //!
- AliL3Benchmark *fBenchmark; //!
-
-#ifdef use_newio
- AliRunLoader *fRunLoader;
-#endif
-
- void CleanUp();
- Double_t GetCpuTime();
-
public:
AliL3Hough();
- AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments=100,Bool_t bit8=kFALSE,Int_t tv=0,Char_t *infile=0,Char_t *ptr=0);
+ AliL3Hough(Char_t *path,Bool_t binary,Int_t netasegments=100,Bool_t bit8=kFALSE,Int_t tv=0,Char_t *infile=0,Char_t *ptr=0);
virtual ~AliL3Hough();
#ifdef use_newio
void SetRunLoader(AliRunLoader *runloader) {fRunLoader = runloader;}
#endif
- void Init(Char_t *path,Bool_t binary,Int_t n_eta_segments=100,Bool_t bit8=kFALSE,Int_t tv=0,Char_t *infile=0,Char_t *ptr=0,Float_t zvertex=0.0);
+ void Init(Char_t *path,Bool_t binary,Int_t netasegments=100,Bool_t bit8=kFALSE,Int_t tv=0,Char_t *infile=0,Char_t *ptr=0,Float_t zvertex=0.0);
void Init(Bool_t doit=kFALSE, Bool_t addhists=kFALSE);
void Process(Int_t minslice,Int_t maxslice);
void ReadData(Int_t slice,Int_t eventnr=0);
- void Transform(Int_t *row_range = 0);
+ void Transform(Int_t *rowrange = 0);
void ProcessSliceIter();
void ProcessPatchIter(Int_t patch);
void MergePatches();
void FindTrackCandidatesRow();
void AddAllHistograms();
void AddAllHistogramsRows();
- Int_t Evaluate(Int_t road_width=1,Int_t nrowstomiss=1);
- void EvaluatePatch(Int_t i,Int_t road_width,Int_t nrowstomiss);
+ Int_t Evaluate(Int_t roadwidth=1,Int_t nrowstomiss=1);
+ void EvaluatePatch(Int_t i,Int_t roadwidth,Int_t nrowstomiss);
void WriteTracks(Int_t slice,Char_t *path="./");
void WriteTracks(Char_t *path);
void WriteDigits(Char_t *outfile="output_digits.root");
AliL3HoughIntMerger *GetInterMerger() {if(!fInterMerger) return 0; return fInterMerger;}
AliL3MemHandler *GetMemHandler(Int_t i) {if(!fMemHandler[i]) return 0; return fMemHandler[i];}
AliL3HoughMaxFinder *GetMaxFinder() {return fPeakFinder;}
+
+ private:
+ Char_t *fInputFile;//!
+ Char_t *fInputPtr;//!
+ Char_t fPath[1024]; // Path to the files
+ Bool_t fBinary; // Is input binary
+ Bool_t fAddHistograms; // Add all patch histograms at the end or not
+ Bool_t fDoIterative; // Iterative or not
+ Bool_t fWriteDigits; // Write Digits or not
+ Bool_t fUse8bits; // Use 8 bits or not
+ Int_t fNEtaSegments; // Number of eta slices
+ Int_t fNPatches; // Number of patches
+ Int_t fVersion; //which HoughTransformer to use
+ Int_t fCurrentSlice; // Current eta slice
+ Int_t fEvent; // Current event number
+
+ Int_t fPeakThreshold[6]; // Threshold for the peak finder
+ Float_t fLowPt[6]; // Lower limit on Pt
+ Float_t fUpperPt[6]; // Upper limit on Pt
+ Float_t fPtRes[6]; // Desired Pt resolution
+ Float_t fPhi[6]; // Limit on the emission angle
+ Int_t fNBinX[6]; // Number of bins in the Hough space
+ Int_t fNBinY[6]; // Number of bins in the Hough space
+ Int_t fThreshold[6]; // Threshold for digits
+ Int_t fNSaveIterations; //for HoughtransformerVhdl
+
+ //parameters for the peak finder:
+ Int_t fKappaSpread; // Kappa spread
+ Float_t fPeakRatio; // Peak ratio
+
+ Float_t fZVertex; // Z position of the primary vertex
+
+ AliL3MemHandler **fMemHandler; //!
+ AliL3HoughBaseTransformer **fHoughTransformer; //!
+ AliL3HoughEval **fEval; //!
+ AliL3HoughMaxFinder *fPeakFinder; //!
+ AliL3TrackArray **fTracks; //!
+ AliL3TrackArray *fGlobalTracks; //!
+ AliL3HoughMerger *fMerger; //!
+ AliL3HoughIntMerger *fInterMerger; //!
+ AliL3HoughGlobalMerger *fGlobalMerger; //!
+ AliL3Benchmark *fBenchmark; //!
+
+#ifdef use_newio
+ AliRunLoader *fRunLoader; // Run Loader
+#endif
+
+ void CleanUp();
+ Double_t GetCpuTime();
ClassDef(AliL3Hough,1) //Hough transform base class
};
void AliL3HoughMaxFinder::Reset()
{
+ // Method to reinit the Peak Finder
for(Int_t i=0; i<fNMax; i++)
{
fXPeaks[i]=0;
void AliL3HoughMaxFinder::CreateNtuppel()
{
+ // Fill a NTuple with the peak parameters
#ifndef no_root
//content#; neighbouring bins of the peak.
fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
{
+ // Write the NTuple with the peak parameters
#ifndef no_root
TFile *file = TFile::Open(filename,"RECREATE");
if(!file)
void AliL3HoughMaxFinder::FindAbsMaxima()
{
-
+ // Simple Peak Finder in the Hough space
if(!fCurrentHisto)
{
cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
Int_t ymin = hist->GetFirstYbin();
Int_t ymax = hist->GetLastYbin();
Int_t bin;
- Double_t value,max_value=0;
+ Double_t value,maxvalue=0;
- Int_t max_xbin=0,max_ybin=0;
+ Int_t maxxbin=0,maxybin=0;
for(Int_t xbin=xmin; xbin<=xmax; xbin++)
{
for(Int_t ybin=ymin; ybin<=ymax; ybin++)
{
bin = hist->GetBin(xbin,ybin);
value = hist->GetBinContent(bin);
- if(value>max_value)
+ if(value>maxvalue)
{
- max_value = value;
- max_xbin = xbin;
- max_ybin = ybin;
+ maxvalue = value;
+ maxxbin = xbin;
+ maxybin = ybin;
}
}
}
- if(max_value == 0)
+ if(maxvalue == 0)
return;
if(fNPeaks > fNMax)
return;
}
- Double_t max_x = hist->GetBinCenterX(max_xbin);
- Double_t max_y = hist->GetBinCenterY(max_ybin);
- fXPeaks[fNPeaks] = max_x;
- fYPeaks[fNPeaks] = max_y;
- fWeight[fNPeaks] = (Int_t)max_value;
+ Double_t maxx = hist->GetBinCenterX(maxxbin);
+ Double_t maxy = hist->GetBinCenterY(maxybin);
+ fXPeaks[fNPeaks] = maxx;
+ fYPeaks[fNPeaks] = maxy;
+ fWeight[fNPeaks] = (Int_t)maxvalue;
fNPeaks++;
#ifndef no_root
if(fNtuppel)
{
- Int_t bin3 = hist->GetBin(max_xbin-1,max_ybin);
- Int_t bin5 = hist->GetBin(max_xbin+1,max_ybin);
- Int_t bin1 = hist->GetBin(max_xbin,max_ybin-1);
- Int_t bin7 = hist->GetBin(max_xbin,max_ybin+1);
+ Int_t bin3 = hist->GetBin(maxxbin-1,maxybin);
+ Int_t bin5 = hist->GetBin(maxxbin+1,maxybin);
+ Int_t bin1 = hist->GetBin(maxxbin,maxybin-1);
+ Int_t bin7 = hist->GetBin(maxxbin,maxybin+1);
- fNtuppel->Fill(max_x,max_y,max_value,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
+ fNtuppel->Fill(maxx,maxy,maxvalue,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
}
#endif
}
void AliL3HoughMaxFinder::FindBigMaxima()
{
-
+ // Another Peak finder
AliL3Histogram *hist = fCurrentHisto;
if(hist->GetNEntries() == 0)
Int_t xmax = hist->GetLastXbin();
Int_t ymin = hist->GetFirstYbin();
Int_t ymax = hist->GetLastYbin();
- Int_t bin[25],bin_index;
+ Int_t bin[25],binindex;
Double_t value[25];
for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
{
for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
{
- bin_index=0;
+ binindex=0;
for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
{
for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
{
- bin[bin_index]=hist->GetBin(xb,yb);
- value[bin_index]=hist->GetBinContent(bin[bin_index]);
- bin_index++;
+ bin[binindex]=hist->GetBin(xb,yb);
+ value[binindex]=hist->GetBinContent(bin[binindex]);
+ binindex++;
}
}
if(value[12]==0) continue;
Int_t b=0;
while(1)
{
- if(value[b] > value[12] || b==bin_index) break;
+ if(value[b] > value[12] || b==binindex) break;
b++;
//printf("b %d\n",b);
}
- if(b == bin_index)
+ if(b == binindex)
{
//Found maxima
if(fNPeaks > fNMax)
return;
}
- Double_t max_x = hist->GetBinCenterX(xbin);
- Double_t max_y = hist->GetBinCenterY(ybin);
- fXPeaks[fNPeaks] = max_x;
- fYPeaks[fNPeaks] = max_y;
+ Double_t maxx = hist->GetBinCenterX(xbin);
+ Double_t maxy = hist->GetBinCenterY(ybin);
+ fXPeaks[fNPeaks] = maxx;
+ fYPeaks[fNPeaks] = maxy;
fNPeaks++;
}
}
&& value[4]>value[7] && value[4]>value[8])
{
//Found a local maxima
- Float_t max_x = fCurrentHisto->GetBinCenterX(xbin);
- Float_t max_y = fCurrentHisto->GetBinCenterY(ybin);
+ Float_t maxx = fCurrentHisto->GetBinCenterX(xbin);
+ Float_t maxy = fCurrentHisto->GetBinCenterY(ybin);
if((Int_t)value[4] <= threshold) continue;//central bin below threshold
if(fNPeaks >= fNMax)
if(value[1]/value[4] > fGradY && value[7]/value[4] > fGradY)
continue;
- fXPeaks[fNPeaks] = max_x;
- fYPeaks[fNPeaks] = max_y;
+ fXPeaks[fNPeaks] = maxx;
+ fYPeaks[fNPeaks] = maxy;
fWeight[fNPeaks] = (Int_t)value[4];
fNPeaks++;
Bool_t bigger = kFALSE;
for(Int_t p=0; p<fNPeaks; p++)
{
- if(fabs(max_x - fXPeaks[p]) < max_kappa && fabs(max_y - fYPeaks[p]) < max_phi0)
+ if(fabs(maxx - fXPeaks[p]) < max_kappa && fabs(maxy - fYPeaks[p]) < max_phi0)
{
bigger = kTRUE;
if(value[4] > fWeight[p]) //this peak is bigger.
{
- fXPeaks[p] = max_x;
- fYPeaks[p] = max_y;
+ fXPeaks[p] = maxx;
+ fYPeaks[p] = maxy;
fWeight[p] = (Int_t)value[4];
}
else
}
if(!bigger) //there were no overlapping peaks.
{
- fXPeaks[fNPeaks] = max_x;
- fYPeaks[fNPeaks] = max_y;
+ fXPeaks[fNPeaks] = maxx;
+ fYPeaks[fNPeaks] = maxy;
fWeight[fNPeaks] = (Int_t)value[4];
fNPeaks++;
}
}
-struct Window
+struct AliL3Window
{
- Int_t start;
- Int_t sum;
+ Int_t fStart; // Start
+ Int_t fSum; // Sum
};
-void AliL3HoughMaxFinder::FindAdaptedPeaks(Int_t kappawindow,Float_t cut_ratio)
+void AliL3HoughMaxFinder::FindAdaptedPeaks(Int_t kappawindow,Float_t cutratio)
{
//Peak finder which looks for peaks with a certain shape.
//The first step involves a pre-peak finder, which looks for peaks
//Start by looking for pre-peaks:
- Window **local_maxima = new Window*[hist->GetNbinsY()];
+ AliL3Window **localmaxima = new AliL3Window*[hist->GetNbinsY()];
Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
- Int_t n,last_sum,sum;
- Bool_t sum_was_rising;
+ Int_t n,lastsum,sum;
+ Bool_t sumwasrising;
for(Int_t ybin=ymin; ybin<=ymax; ybin++)
{
- local_maxima[ybin-ymin] = new Window[hist->GetNbinsX()];
+ localmaxima[ybin-ymin] = new AliL3Window[hist->GetNbinsX()];
nmaxs[ybin-ymin] = 0;
- sum_was_rising=0;
- last_sum=0;
+ sumwasrising=0;
+ lastsum=0;
n=0;
for(Int_t xbin=xmin; xbin<=xmax-kappawindow; xbin++)
{
for(Int_t lbin=xbin; lbin<xbin+kappawindow; lbin++)
sum += hist->GetBinContent(hist->GetBin(lbin,ybin));
- if(sum < last_sum)
+ if(sum < lastsum)
{
if(sum > fThreshold)
- if(sum_was_rising)//Previous sum was a local maxima
+ if(sumwasrising)//Previous sum was a local maxima
{
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].start = xbin-1;
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].sum = last_sum;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fStart = xbin-1;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fSum = lastsum;
nmaxs[ybin-ymin]++;
}
- sum_was_rising=0;
+ sumwasrising=0;
}
else if(sum > 0)
- sum_was_rising=1;
- last_sum=sum;
+ sumwasrising=1;
+ lastsum=sum;
}
}
{
for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
{
- Int_t lw = local_maxima[ybin-ymin][i].sum;
+ Int_t lw = localmaxima[ybin-ymin][i].fSum;
if(lw<0)
continue; //already used
Int_t maxvalue=0,maxybin=0,maxxbin=0,maxwindow=0;
- for(Int_t k=local_maxima[ybin-ymin][i].start; k<local_maxima[ybin-ymin][i].start + kappawindow; k++)
+ for(Int_t k=localmaxima[ybin-ymin][i].fStart; k<localmaxima[ybin-ymin][i].fStart + kappawindow; k++)
if(hist->GetBinContent(hist->GetBin(k,ybin)) > maxvalue)
{
maxvalue = hist->GetBinContent(hist->GetBin(k,ybin));
//start expanding in the psi-direction:
- Int_t lb = local_maxima[ybin-ymin][i].start;
+ Int_t lb = localmaxima[ybin-ymin][i].fStart;
//Int_t ystart=ybin;
- starts[ybin] = local_maxima[ybin-ymin][i].start;
+ starts[ybin] = localmaxima[ybin-ymin][i].fStart;
maxs[ybin] = maxxbin;
Int_t yl=ybin-1,nybins=1;
- //cout<<"Starting search at ybin "<<ybin<<" start "<<lb<<" with sum "<<local_maxima[ybin-ymin][i].sum<<endl;
+ //cout<<"Starting search at ybin "<<ybin<<" start "<<lb<<" with sum "<<localmaxima[ybin-ymin][i].sum<<endl;
while(yl >= ymin)
{
Bool_t found=0;
for(Int_t j=0; j<nmaxs[yl-ymin]; j++)
{
- if( local_maxima[yl-ymin][j].start - lb < 0) continue;
- if( local_maxima[yl-ymin][j].start < lb + kappawindow + match &&
- local_maxima[yl-ymin][j].start >= lb && local_maxima[yl-ymin][j].sum > 0)
+ if( localmaxima[yl-ymin][j].fStart - lb < 0) continue;
+ if( localmaxima[yl-ymin][j].fStart < lb + kappawindow + match &&
+ localmaxima[yl-ymin][j].fStart >= lb && localmaxima[yl-ymin][j].fSum > 0)
{
- //cout<<"match at ybin "<<yl<<" yvalue "<<hist->GetBinCenterY(yl)<<" start "<<local_maxima[yl-ymin][j].start<<" sum "<<local_maxima[yl-ymin][j].sum<<endl;
+ //cout<<"match at ybin "<<yl<<" yvalue "<<hist->GetBinCenterY(yl)<<" start "<<localmaxima[yl-ymin][j].start<<" sum "<<localmaxima[yl-ymin][j].sum<<endl;
Int_t lmaxvalue=0,lmaxxbin=0;
- for(Int_t k=local_maxima[yl-ymin][j].start; k<local_maxima[yl-ymin][j].start + kappawindow; k++)
+ for(Int_t k=localmaxima[yl-ymin][j].fStart; k<localmaxima[yl-ymin][j].fStart + kappawindow; k++)
{
if(hist->GetBinContent(hist->GetBin(k,yl)) > maxvalue)
{
}
}
nybins++;
- starts[yl] = local_maxima[yl-ymin][j].start;
+ starts[yl] = localmaxima[yl-ymin][j].fStart;
maxs[yl] = lmaxxbin;
- local_maxima[yl-ymin][j].sum=-1; //Mark as used
+ localmaxima[yl-ymin][j].fSum=-1; //Mark as used
found=1;
- lb = local_maxima[yl-ymin][j].start;
+ lb = localmaxima[yl-ymin][j].fStart;
break;//Since we found a match in this bin, we dont have to search it anymore, goto next bin.
}
}
//cout<<"ratio "<<right/left<<endl;
- Float_t upper_ratio=1,lower_ratio=1;
+ Float_t upperratio=1,lowerratio=1;
if(left)
- upper_ratio = right/left;
+ upperratio = right/left;
right=left=0;
for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
//cout<<"ratio "<<left/right<<endl;
if(right)
- lower_ratio = left/right;
+ lowerratio = left/right;
- if(upper_ratio > cut_ratio || lower_ratio > cut_ratio)
+ if(upperratio > cutratio || lowerratio > cutratio)
truepeak=kFALSE;
if(truepeak)
}
for(Int_t i=0; i<hist->GetNbinsY(); i++)
- delete local_maxima[i];
+ delete localmaxima[i];
- delete [] local_maxima;
+ delete [] localmaxima;
delete [] nmaxs;
delete [] starts;
delete [] maxs;
}
-struct PreYPeak
+struct AliL3PreYPeak
{
- Int_t start_position;
- Int_t end_position;
- Int_t min_value;
- Int_t max_value;
- Int_t prev_value;
- Int_t left_value;
- Int_t right_value;
+ Int_t fStartPosition; // Start position in X
+ Int_t fEndPosition; // End position in X
+ Int_t fMinValue; // Minimum value inside the prepeak
+ Int_t fMaxValue; // Maximum value inside the prepeak
+ Int_t fPrevValue; // Neighbour values
+ Int_t fLeftValue; // Neighbour values
+ Int_t fRightValue; // Neighbour values
};
-struct Pre2DPeak
+struct AliL3Pre2DPeak
{
- Float_t x;
- Float_t y;
- Float_t size_x;
- Float_t size_y;
- Int_t start_x;
- Int_t start_y;
- Int_t end_x;
- Int_t end_y;
- Float_t weight;
+ Float_t fX; // X coordinate of the preak
+ Float_t fY; // Y coordinate of the preak
+ Float_t fSizeX; // Size of the peak
+ Float_t fSizeY; // Size of the peak
+ Int_t fStartX; // Start position of the peak
+ Int_t fStartY; // Start position of the peak
+ Int_t fEndX; // End position of the peak
+ Int_t fEndY; // End position of the peak
+ Float_t fWeight; // Weight assigned to the peak
};
void AliL3HoughMaxFinder::FindAdaptedRowPeaks(Int_t kappawindow,Int_t xsize,Int_t ysize)
{
-
+ // Peak finder which is working over the Hough Space provided by the AliL3HoughTransformerRow class
AliL3Histogram *hist = fCurrentHisto;
if(!hist)
//Start by looking for pre-peaks:
- PreYPeak **local_maxima = new PreYPeak*[hist->GetNbinsY()];
+ AliL3PreYPeak **localmaxima = new AliL3PreYPeak*[hist->GetNbinsY()];
Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
- Int_t last_value=0,value=0;
+ Int_t lastvalue=0,value=0;
for(Int_t ybin=ymin; ybin<=ymax; ybin++)
{
- local_maxima[ybin-ymin] = new PreYPeak[hist->GetNbinsX()];
+ localmaxima[ybin-ymin] = new AliL3PreYPeak[hist->GetNbinsX()];
nmaxs[ybin-ymin] = 0;
- last_value = 0;
+ lastvalue = 0;
Bool_t found = 0;
for(Int_t xbin=xmin; xbin<=xmax; xbin++)
{
value = hist->GetBinContent(hist->GetBin(xbin,ybin));
if(value > 0)
{
- if(abs(value - last_value) > 1)
+ if(abs(value - lastvalue) > 1)
{
if(found) {
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].right_value = value;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value;
nmaxs[ybin-ymin]++;
}
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].start_position = xbin;
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].end_position = xbin;
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].min_value = value;
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].max_value = value;
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].prev_value = 0;
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].left_value = last_value;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fStartPosition = xbin;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fEndPosition = xbin;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue = value;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue = value;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fPrevValue = 0;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fLeftValue = lastvalue;
found = 1;
}
- if(abs(value - last_value) <= 1)
+ if(abs(value - lastvalue) <= 1)
{
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].end_position = xbin;
- if(value>local_maxima[ybin-ymin][nmaxs[ybin-ymin]].max_value)
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].max_value = value;
- if(value<local_maxima[ybin-ymin][nmaxs[ybin-ymin]].min_value)
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].min_value = value;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fEndPosition = xbin;
+ if(value>localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue)
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue = value;
+ if(value<localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue)
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue = value;
}
}
- last_value = value;
+ lastvalue = value;
}
if(found) {
- local_maxima[ybin-ymin][nmaxs[ybin-ymin]].right_value = value;
+ localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value;
nmaxs[ybin-ymin]++;
}
}
- Pre2DPeak maxima[500];
+ AliL3Pre2DPeak maxima[500];
Int_t nmaxima = 0;
for(Int_t ybin=ymax; ybin >= ymin; ybin--)
{
for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
{
- Int_t local_min_value = local_maxima[ybin-ymin][i].min_value;
- Int_t local_max_value = local_maxima[ybin-ymin][i].max_value;
- Int_t local_prev_value = local_maxima[ybin-ymin][i].prev_value;
- Int_t local_next_value = 0;
- Int_t local_left_value = local_maxima[ybin-ymin][i].left_value;
- Int_t local_right_value = local_maxima[ybin-ymin][i].right_value;
-
- if(local_min_value<0)
+ Int_t localminvalue = localmaxima[ybin-ymin][i].fMinValue;
+ Int_t localmaxvalue = localmaxima[ybin-ymin][i].fMaxValue;
+ Int_t localprevvalue = localmaxima[ybin-ymin][i].fPrevValue;
+ Int_t localnextvalue = 0;
+ Int_t localleftvalue = localmaxima[ybin-ymin][i].fLeftValue;
+ Int_t localrightvalue = localmaxima[ybin-ymin][i].fRightValue;
+
+ if(localminvalue<0)
continue; //already used
//start expanding in the psi-direction:
- Int_t local_x_start = local_maxima[ybin-ymin][i].start_position;
- Int_t local_x_end = local_maxima[ybin-ymin][i].end_position;
- Int_t temp_x_start = local_maxima[ybin-ymin][i].start_position;
- Int_t temp_x_end = local_maxima[ybin-ymin][i].end_position;
+ Int_t localxstart = localmaxima[ybin-ymin][i].fStartPosition;
+ Int_t localxend = localmaxima[ybin-ymin][i].fEndPosition;
+ Int_t tempxstart = localmaxima[ybin-ymin][i].fStartPosition;
+ Int_t tempxend = localmaxima[ybin-ymin][i].fEndPosition;
- Int_t local_y=ybin-1,nybins=1;
+ Int_t localy=ybin-1,nybins=1;
- while(local_y >= ymin)
+ while(localy >= ymin)
{
Bool_t found=0;
- for(Int_t j=0; j<nmaxs[local_y-ymin]; j++)
+ for(Int_t j=0; j<nmaxs[localy-ymin]; j++)
{
- if( (local_maxima[local_y-ymin][j].start_position <= (temp_x_end + kappawindow)) && (local_maxima[local_y-ymin][j].end_position >= (temp_x_start - kappawindow)))
+ if( (localmaxima[localy-ymin][j].fStartPosition <= (tempxend + kappawindow)) && (localmaxima[localy-ymin][j].fEndPosition >= (tempxstart - kappawindow)))
{
- if(((local_maxima[local_y-ymin][j].min_value <= local_max_value) && (local_maxima[local_y-ymin][j].min_value >= local_min_value)) ||
- ((local_maxima[local_y-ymin][j].max_value >= local_min_value) && (local_maxima[local_y-ymin][j].max_value <= local_max_value)))
+ if(((localmaxima[localy-ymin][j].fMinValue <= localmaxvalue) && (localmaxima[localy-ymin][j].fMinValue >= localminvalue)) ||
+ ((localmaxima[localy-ymin][j].fMaxValue >= localminvalue) && (localmaxima[localy-ymin][j].fMaxValue <= localmaxvalue)))
{
- if(local_maxima[local_y-ymin][j].end_position > local_x_end)
- local_x_end = local_maxima[local_y-ymin][j].end_position;
- if(local_maxima[local_y-ymin][j].start_position < local_x_start)
- local_x_start = local_maxima[local_y-ymin][j].start_position;
- temp_x_start = local_maxima[local_y-ymin][j].start_position;
- temp_x_end = local_maxima[local_y-ymin][j].end_position;
- if(local_maxima[local_y-ymin][j].min_value < local_min_value)
- local_min_value = local_maxima[local_y-ymin][j].min_value;
- if(local_maxima[local_y-ymin][j].max_value > local_max_value)
- local_max_value = local_maxima[local_y-ymin][j].max_value;
- if(local_maxima[local_y-ymin][j].right_value > local_right_value)
- local_right_value = local_maxima[local_y-ymin][j].right_value;
- if(local_maxima[local_y-ymin][j].left_value > local_left_value)
- local_left_value = local_maxima[local_y-ymin][j].left_value;
- local_maxima[local_y-ymin][j].min_value = -1;
+ if(localmaxima[localy-ymin][j].fEndPosition > localxend)
+ localxend = localmaxima[localy-ymin][j].fEndPosition;
+ if(localmaxima[localy-ymin][j].fStartPosition < localxstart)
+ localxstart = localmaxima[localy-ymin][j].fStartPosition;
+ tempxstart = localmaxima[localy-ymin][j].fStartPosition;
+ tempxend = localmaxima[localy-ymin][j].fEndPosition;
+ if(localmaxima[localy-ymin][j].fMinValue < localminvalue)
+ localminvalue = localmaxima[localy-ymin][j].fMinValue;
+ if(localmaxima[localy-ymin][j].fMaxValue > localmaxvalue)
+ localmaxvalue = localmaxima[localy-ymin][j].fMaxValue;
+ if(localmaxima[localy-ymin][j].fRightValue > localrightvalue)
+ localrightvalue = localmaxima[localy-ymin][j].fRightValue;
+ if(localmaxima[localy-ymin][j].fLeftValue > localleftvalue)
+ localleftvalue = localmaxima[localy-ymin][j].fLeftValue;
+ localmaxima[localy-ymin][j].fMinValue = -1;
found = 1;
nybins++;
break;
}
else
{
- if(local_max_value > local_maxima[local_y-ymin][j].prev_value)
- local_maxima[local_y-ymin][j].prev_value = local_max_value;
- if(local_maxima[local_y-ymin][j].max_value > local_next_value)
- local_next_value = local_maxima[local_y-ymin][j].max_value;
+ if(localmaxvalue > localmaxima[localy-ymin][j].fPrevValue)
+ localmaxima[localy-ymin][j].fPrevValue = localmaxvalue;
+ if(localmaxima[localy-ymin][j].fMaxValue > localnextvalue)
+ localnextvalue = localmaxima[localy-ymin][j].fMaxValue;
}
}
}
- if(!found || local_y == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
+ if(!found || localy == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
{
- if((nybins > ysize) && ((local_x_end-local_x_start+1) > xsize) && (local_max_value > local_prev_value) && (local_max_value > local_next_value) && (local_max_value > local_left_value) && (local_max_value > local_right_value))
- // if((nybins > ysize) && ((local_x_end-local_x_start+1) > xsize))
+ if((nybins > ysize) && ((localxend-localxstart+1) > xsize) && (localmaxvalue > localprevvalue) && (localmaxvalue > localnextvalue) && (localmaxvalue > localleftvalue) && (localmaxvalue > localrightvalue))
+ // if((nybins > ysize) && ((localxend-localxstart+1) > xsize))
{
- maxima[nmaxima].x = ((Float_t)local_x_start+(Float_t)local_x_end)/2.0;
- maxima[nmaxima].y = ((Float_t)ybin+(Float_t)(local_y+1))/2.0;
- maxima[nmaxima].size_x = local_x_end-local_x_start+1;
- maxima[nmaxima].size_y = nybins;
- maxima[nmaxima].weight = (local_min_value+local_max_value)/2;
- maxima[nmaxima].start_x = local_x_start;
- maxima[nmaxima].end_x = local_x_end;
- maxima[nmaxima].start_y = local_y +1;
- maxima[nmaxima].end_y = ybin;
+ maxima[nmaxima].fX = ((Float_t)localxstart+(Float_t)localxend)/2.0;
+ maxima[nmaxima].fY = ((Float_t)ybin+(Float_t)(localy+1))/2.0;
+ maxima[nmaxima].fSizeX = localxend-localxstart+1;
+ maxima[nmaxima].fSizeY = nybins;
+ maxima[nmaxima].fWeight = (localminvalue+localmaxvalue)/2;
+ maxima[nmaxima].fStartX = localxstart;
+ maxima[nmaxima].fEndX = localxend;
+ maxima[nmaxima].fStartY = localy +1;
+ maxima[nmaxima].fEndY = ybin;
#ifdef do_mc
- // cout<<"Peak found at: "<<((Float_t)local_x_start+(Float_t)local_x_end)/2.0<<" "<<((Float_t)ybin+(Float_t)(local_y+1))/2.0<<" "<<local_min_value<<" "<<local_max_value<<" "<<" with weight "<<(local_min_value+local_max_value)/2<<" and size "<<local_x_end-local_x_start+1<<" by "<<nybins<<endl;
+ // cout<<"Peak found at: "<<((Float_t)localxstart+(Float_t)localxend)/2.0<<" "<<((Float_t)ybin+(Float_t)(localy+1))/2.0<<" "<<localminvalue<<" "<<localmaxvalue<<" "<<" with weight "<<(localminvalue+localmaxvalue)/2<<" and size "<<localxend-localxstart+1<<" by "<<nybins<<endl;
#endif
nmaxima++;
}
break;
}
else
- local_y--;//Search continues...
+ localy--;//Search continues...
}
}
}
for(Int_t i = 0; i < (nmaxima - 1); i++)
{
- if(maxima[i].weight < 0) continue;
+ if(maxima[i].fWeight < 0) continue;
for(Int_t j = i + 1; j < nmaxima; j++)
{
- if(maxima[j].weight < 0) continue;
+ if(maxima[j].fWeight < 0) continue;
Int_t xtrack1=0,xtrack2=0,ytrack1=0,ytrack2=0;
Int_t deltax = 9999;
- for(Int_t ix1 = maxima[i].start_x; ix1 <= maxima[i].end_x; ix1++) {
- for(Int_t ix2 = maxima[j].start_x; ix2 <= maxima[j].end_x; ix2++) {
+ for(Int_t ix1 = maxima[i].fStartX; ix1 <= maxima[i].fEndX; ix1++) {
+ for(Int_t ix2 = maxima[j].fStartX; ix2 <= maxima[j].fEndX; ix2++) {
if(abs(ix1 - ix2) < deltax) {
deltax = abs(ix1 - ix2);
xtrack1 = ix1;
}
}
Int_t deltay = 9999;
- for(Int_t iy1 = maxima[i].start_y; iy1 <= maxima[i].end_y; iy1++) {
- for(Int_t iy2 = maxima[j].start_y; iy2 <= maxima[j].end_y; iy2++) {
+ for(Int_t iy1 = maxima[i].fStartY; iy1 <= maxima[i].fEndY; iy1++) {
+ for(Int_t iy2 = maxima[j].fStartY; iy2 <= maxima[j].fEndY; iy2++) {
if(abs(iy1 - iy2) < deltay) {
deltay = abs(iy1 - iy2);
ytrack1 = iy1;
#endif
//cvetan please check!!! I added a cast to Int_t
if((fabs(firsthit1[1]-firsthit2[1]) < 3.0*padpitchlow) && (fabs(lasthit1[1]-lasthit2[1]) < 3.0*padpitchup)) {
- if(maxima[i].size_x*maxima[i].size_y > maxima[j].size_x*maxima[j].size_y)
- maxima[j].weight = -maxima[j].weight;
- if(maxima[i].size_x*maxima[i].size_y < maxima[j].size_x*maxima[j].size_y)
- maxima[i].weight = -maxima[i].weight;
+ if(maxima[i].fSizeX*maxima[i].fSizeY > maxima[j].fSizeX*maxima[j].fSizeY)
+ maxima[j].fWeight = -maxima[j].fWeight;
+ if(maxima[i].fSizeX*maxima[i].fSizeY < maxima[j].fSizeX*maxima[j].fSizeY)
+ maxima[i].fWeight = -maxima[i].fWeight;
#ifdef do_mc
- // cout<<"Merge peaks "<<i<<" "<<j<<" "<<maxima[i].weight<<" "<<maxima[j].weight<<endl;
+ // cout<<"Merge peaks "<<i<<" "<<j<<" "<<maxima[i].fWeight<<" "<<maxima[j].fWeight<<endl;
#endif
}
}
//merge tracks in neighbour eta slices
/*
for(Int_t i = 0; i < nmaxima; i++) {
- if(maxima[i].weight > 0) {
- fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].x);
- fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].y);
- fWeight[fNPeaks] = (Int_t)maxima[i].weight;
+ if(maxima[i].fWeight > 0) {
+ fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].fX);
+ fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].fY);
+ fWeight[fNPeaks] = (Int_t)maxima[i].fWeight;
#ifdef do_mc
- cout<<"Final Peak found at: "<<maxima[i].x<<" "<<maxima[i].y<<" "<<" "<<fXPeaks[fNPeaks]<<" "<<fYPeaks[fNPeaks]<<" with weight "<<fWeight[fNPeaks]<<" and size "<<maxima[i].size_x<<" by "<<maxima[i].size_y<<endl;
+ cout<<"Final Peak found at: "<<maxima[i].fX<<" "<<maxima[i].fY<<" "<<" "<<fXPeaks[fNPeaks]<<" "<<fYPeaks[fNPeaks]<<" with weight "<<fWeight[fNPeaks]<<" and size "<<maxima[i].fSizeX<<" by "<<maxima[i].fSizeY<<endl;
#endif
fNPeaks++;
}
Int_t currentnpeaks = fNPeaks;
for(Int_t i = 0; i < nmaxima; i++) {
- if(maxima[i].weight < 0) continue;
+ if(maxima[i].fWeight < 0) continue;
Bool_t merged = kFALSE;
for(Int_t j = fN1PeaksPrevEtaSlice; j < fN2PeaksPrevEtaSlice; j++) {
if(fWeight[j] < 0) continue;
if((fENDETAPeaks[j]-fSTARTETAPeaks[j]) >= 1) continue;
- if((maxima[i].start_x <= fENDXPeaks[j]+1) && (maxima[i].end_x >= fSTARTXPeaks[j]-1)) {
- if((maxima[i].start_y <= fENDYPeaks[j]+1) && (maxima[i].end_y >= fSTARTYPeaks[j]-1)) {
+ if((maxima[i].fStartX <= fENDXPeaks[j]+1) && (maxima[i].fEndX >= fSTARTXPeaks[j]-1)) {
+ if((maxima[i].fStartY <= fENDYPeaks[j]+1) && (maxima[i].fEndY >= fSTARTYPeaks[j]-1)) {
//merge
merged = kTRUE;
- fXPeaks[fNPeaks] = (hist->GetPreciseBinCenterX(maxima[i].x)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fXPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
- fYPeaks[fNPeaks] = (hist->GetPreciseBinCenterY(maxima[i].y)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fYPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
- fWeight[fNPeaks] = (Int_t)maxima[i].weight + fWeight[j];
- fSTARTXPeaks[fNPeaks] = maxima[i].start_x;
- fSTARTYPeaks[fNPeaks] = maxima[i].start_y;
- fENDXPeaks[fNPeaks] = maxima[i].end_x;
- fENDYPeaks[fNPeaks] = maxima[i].end_y;
+ fXPeaks[fNPeaks] = (hist->GetPreciseBinCenterX(maxima[i].fX)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fXPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
+ fYPeaks[fNPeaks] = (hist->GetPreciseBinCenterY(maxima[i].fY)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fYPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
+ fWeight[fNPeaks] = (Int_t)maxima[i].fWeight + fWeight[j];
+ fSTARTXPeaks[fNPeaks] = maxima[i].fStartX;
+ fSTARTYPeaks[fNPeaks] = maxima[i].fStartY;
+ fENDXPeaks[fNPeaks] = maxima[i].fEndX;
+ fENDYPeaks[fNPeaks] = maxima[i].fEndY;
fSTARTETAPeaks[fNPeaks] = fSTARTETAPeaks[j];
fENDETAPeaks[fNPeaks] = fCurrentEtaSlice;
fNPeaks++;
}
}
}
- fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].x);
- fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].y);
+ fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].fX);
+ fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].fY);
if(!merged)
- fWeight[fNPeaks] = (Int_t)maxima[i].weight;
+ fWeight[fNPeaks] = (Int_t)maxima[i].fWeight;
else
- fWeight[fNPeaks] = -(Int_t)maxima[i].weight;
- fSTARTXPeaks[fNPeaks] = maxima[i].start_x;
- fSTARTYPeaks[fNPeaks] = maxima[i].start_y;
- fENDXPeaks[fNPeaks] = maxima[i].end_x;
- fENDYPeaks[fNPeaks] = maxima[i].end_y;
+ fWeight[fNPeaks] = -(Int_t)maxima[i].fWeight;
+ fSTARTXPeaks[fNPeaks] = maxima[i].fStartX;
+ fSTARTYPeaks[fNPeaks] = maxima[i].fStartY;
+ fENDXPeaks[fNPeaks] = maxima[i].fEndX;
+ fENDYPeaks[fNPeaks] = maxima[i].fEndY;
fSTARTETAPeaks[fNPeaks] = fCurrentEtaSlice;
fENDETAPeaks[fNPeaks] = fCurrentEtaSlice;
fNPeaks++;
fN2PeaksPrevEtaSlice = fNPeaks;
for(Int_t i=0; i<hist->GetNbinsY(); i++)
- delete local_maxima[i];
+ delete localmaxima[i];
- delete [] local_maxima;
+ delete [] localmaxima;
delete [] nmaxs;
}
-void AliL3HoughMaxFinder::FindPeak1(Int_t y_window,Int_t x_bin_sides)
+void AliL3HoughMaxFinder::FindPeak1(Int_t ywindow,Int_t xbinsides)
{
//Testing mutliple peakfinding.
//The algorithm searches the histogram for prepreaks by looking in windows
- //for each bin on the xaxis. The size of these windows is controlled by y_window.
+ //for each bin on the xaxis. The size of these windows is controlled by ywindow.
//Then the prepreaks are sorted according to their weight (sum inside window),
//and the peak positions are calculated by taking the weighted mean in both
- //x and y direction. The size of the peak in x-direction is controlled by x_bin_sides.
+ //x and y direction. The size of the peak in x-direction is controlled by xbinsides.
if(!fCurrentHisto)
{
if(fCurrentHisto->GetNEntries()==0)
return;
- //Int_t y_window=2;
- //Int_t x_bin_sides=1;
+ //Int_t ywindow=2;
+ //Int_t xbinsides=1;
//Float_t max_kappa = 0.001;
//Float_t max_phi0 = 0.08;
- Int_t max_sum=0;
+ Int_t maxsum=0;
Int_t xmin = fCurrentHisto->GetFirstXbin();
Int_t xmax = fCurrentHisto->GetLastXbin();
Int_t ymax = fCurrentHisto->GetLastYbin();
Int_t nbinsx = fCurrentHisto->GetNbinsX()+1;
- AxisWindow **windowPt = new AxisWindow*[nbinsx];
- AxisWindow **anotherPt = new AxisWindow*[nbinsx];
+ AliL3AxisWindow **windowPt = new AliL3AxisWindow*[nbinsx];
+ AliL3AxisWindow **anotherPt = new AliL3AxisWindow*[nbinsx];
for(Int_t i=0; i<nbinsx; i++)
{
- windowPt[i] = new AxisWindow;
+ windowPt[i] = new AliL3AxisWindow;
#if defined(__DECCXX)
- bzero((char *)windowPt[i],sizeof(AxisWindow));
+ bzero((char *)windowPt[i],sizeof(AliL3AxisWindow));
#else
- bzero((void*)windowPt[i],sizeof(AxisWindow));
+ bzero((void*)windowPt[i],sizeof(AliL3AxisWindow));
#endif
anotherPt[i] = windowPt[i];
}
for(Int_t xbin=xmin; xbin<=xmax; xbin++)
{
- max_sum = 0;
- for(Int_t ybin=ymin; ybin<=ymax-y_window; ybin++)
+ maxsum = 0;
+ for(Int_t ybin=ymin; ybin<=ymax-ywindow; ybin++)
{
- Int_t sum_in_window=0;
- for(Int_t b=ybin; b<ybin+y_window; b++)
+ Int_t suminwindow=0;
+ for(Int_t b=ybin; b<ybin+ywindow; b++)
{
//inside window
Int_t bin = fCurrentHisto->GetBin(xbin,b);
- sum_in_window += (Int_t)fCurrentHisto->GetBinContent(bin);
+ suminwindow += (Int_t)fCurrentHisto->GetBinContent(bin);
}
- if(sum_in_window > max_sum)
+ if(suminwindow > maxsum)
{
- max_sum = sum_in_window;
- windowPt[xbin]->ymin = ybin;
- windowPt[xbin]->ymax = ybin + y_window;
- windowPt[xbin]->weight = sum_in_window;
- windowPt[xbin]->xbin = xbin;
+ maxsum = suminwindow;
+ windowPt[xbin]->fYmin = ybin;
+ windowPt[xbin]->fYmax = ybin + ywindow;
+ windowPt[xbin]->fWeight = suminwindow;
+ windowPt[xbin]->fXbin = xbin;
}
}
}
for(Int_t i=0; i<nbinsx; i++)
{
top=butt=0;
- Int_t xbin = windowPt[i]->xbin;
+ Int_t xbin = windowPt[i]->fXbin;
if(xbin<xmin || xbin > xmax-1) continue;
//Check if this is really a local maxima
- if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight ||
- anotherPt[xbin+1]->weight > anotherPt[xbin]->weight)
+ if(anotherPt[xbin-1]->fWeight > anotherPt[xbin]->fWeight ||
+ anotherPt[xbin+1]->fWeight > anotherPt[xbin]->fWeight)
continue;
- for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++)
+ for(Int_t j=windowPt[i]->fYmin; j<windowPt[i]->fYmax; j++)
{
//Calculate the mean in y direction:
- Int_t bin = fCurrentHisto->GetBin(windowPt[i]->xbin,j);
+ Int_t bin = fCurrentHisto->GetBin(windowPt[i]->fXbin,j);
top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
butt += fCurrentHisto->GetBinContent(bin);
}
if(butt < fThreshold)
continue;
- fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->xbin);
+ fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->fXbin);
fYPeaks[fNPeaks] = top/butt;
fWeight[fNPeaks] = (Int_t)butt;
//cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
for(Int_t i=0; i<fNPeaks; i++)
{
Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]);
- if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue;
+ if(xbin - xbinsides < xmin || xbin + xbinsides > xmax) continue;
top=butt=0;
ytop=0,ybutt=0;
w=0;
- prev = xbin - x_bin_sides+1;
- for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++)
+ prev = xbin - xbinsides+1;
+ for(Int_t j=xbin-xbinsides; j<=xbin+xbinsides; j++)
{
/*
//Check if the windows are overlapping:
prev = j;
*/
- top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->weight;
- butt += anotherPt[j]->weight;
+ top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->fWeight;
+ butt += anotherPt[j]->fWeight;
- for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++)
+ for(Int_t k=anotherPt[j]->fYmin; k<anotherPt[j]->fYmax; k++)
{
Int_t bin = fCurrentHisto->GetBin(j,k);
ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
delete [] anotherPt;
}
-void AliL3HoughMaxFinder::SortPeaks(struct AxisWindow **a,Int_t first,Int_t last)
+void AliL3HoughMaxFinder::SortPeaks(struct AliL3AxisWindow **a,Int_t first,Int_t last)
{
//General sorting routine
//Sort according to PeakCompare()
- static struct AxisWindow *tmp;
+ static struct AliL3AxisWindow *tmp;
static int i; // "static" to save stack space
int j;
}
-Int_t AliL3HoughMaxFinder::PeakCompare(struct AxisWindow *a,struct AxisWindow *b)
+Int_t AliL3HoughMaxFinder::PeakCompare(struct AliL3AxisWindow *a,struct AliL3AxisWindow *b) const
{
- if(a->weight < b->weight) return 1;
- if(a->weight > b->weight) return -1;
+ // Peak comparison based on peaks weight
+ if(a->fWeight < b->fWeight) return 1;
+ if(a->fWeight > b->fWeight) return -1;
return 0;
}
Int_t nbinsx = hist->GetNbinsX()+1;
Int_t *m = new Int_t[nbinsx];
- Int_t *m_low = new Int_t[nbinsx];
- Int_t *m_up = new Int_t[nbinsx];
+ Int_t *mlow = new Int_t[nbinsx];
+ Int_t *mup = new Int_t[nbinsx];
recompute: //this is a goto.
for(Int_t i=0; i<nbinsx; i++)
{
m[i]=0;
- m_low[i]=0;
- m_up[i]=0;
+ mlow[i]=0;
+ mup[i]=0;
}
- Int_t max_x=0,sum=0,max_xbin=0,bin=0;
+ Int_t maxx=0,sum=0,maxxbin=0,bin=0;
for(Int_t xbin=xmin; xbin<=xmax; xbin++)
{
if(sum > m[xbin]) //Max value locally in this xbin
{
m[xbin]=sum;
- m_low[xbin]=ybin;
- m_up[xbin]=ybin + t1;
+ mlow[xbin]=ybin;
+ mup[xbin]=ybin + t1;
}
}
- if(m[xbin] > max_x) //Max value globally in x-direction
+ if(m[xbin] > maxx) //Max value globally in x-direction
{
- max_xbin = xbin;
- max_x = m[xbin];//sum;
+ maxxbin = xbin;
+ maxx = m[xbin];//sum;
}
}
- //printf("max_xbin %d max_x %d m_low %d m_up %d\n",max_xbin,max_x,m_low[max_xbin],m_up[max_xbin]);
- //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
+ //printf("maxxbin %d maxx %d mlow %d mup %d\n",maxxbin,maxx,mlow[maxxbin],mup[maxxbin]);
+ //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[maxxbin]),hist->GetBinCenterY(mup[maxxbin]));
//Determine a width in the x-direction
- Int_t x_low=0,x_up=0;
+ Int_t xlow=0,xup=0;
- for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
+ for(Int_t xbin=maxxbin-1; xbin >= xmin; xbin--)
{
- if(m[xbin] < max_x*t2)
+ if(m[xbin] < maxx*t2)
{
- x_low = xbin+1;
+ xlow = xbin+1;
break;
}
}
- for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
+ for(Int_t xbin = maxxbin+1; xbin <=xmax; xbin++)
{
- if(m[xbin] < max_x*t2)
+ if(m[xbin] < maxx*t2)
{
- x_up = xbin-1;
+ xup = xbin-1;
break;
}
}
- Double_t top=0,butt=0,value,x_peak;
- if(x_up - x_low + 1 > t3)
+ Double_t top=0,butt=0,value,xpeak;
+ if(xup - xlow + 1 > t3)
{
t1 -= 1;
- printf("\nxrange out if limit x_up %d x_low %d t1 %d\n\n",x_low,x_up,t1);
+ printf("\nxrange out if limit xup %d xlow %d t1 %d\n\n",xlow,xup,t1);
if(t1 > 1)
goto recompute;
else
{
- x_peak = hist->GetBinCenterX(max_xbin);
+ xpeak = hist->GetBinCenterX(maxxbin);
goto moveon;
}
}
- //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
- //printf("Spread in x %d\n",x_up-x_low +1);
+ //printf("xlow %f xup %f\n",hist->GetBinCenterX(xlow),hist->GetBinCenterX(xup));
+ //printf("Spread in x %d\n",xup-xlow +1);
//Now, calculate the center of mass in x-direction
- for(Int_t xbin=x_low; xbin <= x_up; xbin++)
+ for(Int_t xbin=xlow; xbin <= xup; xbin++)
{
value = hist->GetBinCenterX(xbin);
top += value*m[xbin];
butt += m[xbin];
}
- x_peak = top/butt;
+ xpeak = top/butt;
moveon:
//Find the peak in y direction:
- Int_t x_l = hist->FindXbin(x_peak);
- if(hist->GetBinCenterX(x_l) > x_peak)
- x_l--;
+ Int_t xl = hist->FindXbin(xpeak);
+ if(hist->GetBinCenterX(xl) > xpeak)
+ xl--;
- Int_t x_u = x_l + 1;
+ Int_t xu = xl + 1;
- if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
- printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
+ if(hist->GetBinCenterX(xl) > xpeak || hist->GetBinCenterX(xu) <= xpeak)
+ printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(xl),xpeak,hist->GetBinCenterX(xu));
- //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
+ //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(xl),hist->GetBinCenterX(xu));
value=top=butt=0;
- //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
- //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
+ //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[xl]),hist->GetBinCenterY(mup[xl]));
+ //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[xu]),hist->GetBinCenterY(mup[xu]));
- for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
+ for(Int_t ybin=mlow[xl]; ybin <= mup[xl]; ybin++)
{
value = hist->GetBinCenterY(ybin);
- bin = hist->GetBin(x_l,ybin);
+ bin = hist->GetBin(xl,ybin);
top += value*hist->GetBinContent(bin);
butt += hist->GetBinContent(bin);
}
- Double_t y_peak_low = top/butt;
+ Double_t ypeaklow = top/butt;
- //printf("y_peak_low %f\n",y_peak_low);
+ //printf("ypeaklow %f\n",ypeaklow);
value=top=butt=0;
- for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
+ for(Int_t ybin=mlow[xu]; ybin <= mup[xu]; ybin++)
{
value = hist->GetBinCenterY(ybin);
- bin = hist->GetBin(x_u,ybin);
+ bin = hist->GetBin(xu,ybin);
top += value*hist->GetBinContent(bin);
butt += hist->GetBinContent(bin);
}
- Double_t y_peak_up = top/butt;
+ Double_t ypeakup = top/butt;
- //printf("y_peak_up %f\n",y_peak_up);
+ //printf("ypeakup %f\n",ypeakup);
- Double_t x_value_up = hist->GetBinCenterX(x_u);
- Double_t x_value_low = hist->GetBinCenterX(x_l);
+ Double_t xvalueup = hist->GetBinCenterX(xu);
+ Double_t xvaluelow = hist->GetBinCenterX(xl);
- Double_t y_peak = (y_peak_low*(x_value_up - x_peak) + y_peak_up*(x_peak - x_value_low))/(x_value_up - x_value_low);
+ Double_t ypeak = (ypeaklow*(xvalueup - xpeak) + ypeakup*(xpeak - xvaluelow))/(xvalueup - xvaluelow);
//Find the weight:
- //bin = hist->FindBin(x_peak,y_peak);
+ //bin = hist->FindBin(xpeak,ypeak);
//Int_t weight = (Int_t)hist->GetBinContent(bin);
//AliL3HoughTrack *track = new AliL3HoughTrack();
- //track->SetTrackParameters(x_peak,y_peak,weight);
- fXPeaks[fNPeaks]=x_peak;
- fYPeaks[fNPeaks]=y_peak;
+ //track->SetTrackParameters(xpeak,ypeak,weight);
+ fXPeaks[fNPeaks]=xpeak;
+ fYPeaks[fNPeaks]=ypeak;
fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin);
fNPeaks++;
delete [] m;
- delete [] m_low;
- delete [] m_up;
+ delete [] mlow;
+ delete [] mup;
//return track;
}
-Float_t AliL3HoughMaxFinder::GetXPeakSize(Int_t i)
+Float_t AliL3HoughMaxFinder::GetXPeakSize(Int_t i) const
{
+ // Get X size of a peak
if(i<0 || i>fNMax)
{
STDCERR<<"AliL3HoughMaxFinder::GetXPeakSize : Invalid index "<<i<<STDENDL;
return 0;
}
- Float_t BinWidth = fCurrentHisto->GetBinWidthX();
- return BinWidth*(fENDXPeaks[i]-fSTARTXPeaks[i]+1);
+ Float_t binwidth = fCurrentHisto->GetBinWidthX();
+ return binwidth*(fENDXPeaks[i]-fSTARTXPeaks[i]+1);
}
-Float_t AliL3HoughMaxFinder::GetYPeakSize(Int_t i)
+Float_t AliL3HoughMaxFinder::GetYPeakSize(Int_t i) const
{
+ // Get Y size of a peak
if(i<0 || i>fNMax)
{
STDCERR<<"AliL3HoughMaxFinder::GetYPeak : Invalid index "<<i<<STDENDL;
return 0;
}
- Float_t BinWidth = fCurrentHisto->GetBinWidthY();
- return BinWidth*(fENDYPeaks[i]-fSTARTYPeaks[i]+1);
+ Float_t binwidth = fCurrentHisto->GetBinWidthY();
+ return binwidth*(fENDYPeaks[i]-fSTARTYPeaks[i]+1);
}
// @(#) $Id$
-#ifndef ALIL3_HOUGH_MaxFinder
-#define ALIL3_HOUGH_MaxFinder
+#ifndef ALIL3HOUGHMAXFINDER_H
+#define ALIL3HOUGHMAXFINDER_H
#include "AliL3RootTypes.h"
#include "AliL3StandardIncludes.h"
class AliL3HoughTrack;
class TNtuple;
-struct AxisWindow
+struct AliL3AxisWindow
{
- Int_t ymin;
- Int_t ymax;
- Int_t xbin;
- Int_t weight;
+ Int_t fYmin; // min Y
+ Int_t fYmax; // max Y
+ Int_t fXbin; // X bin
+ Int_t fWeight; // weight
};
class AliL3HoughMaxFinder {
-
- private:
-
- Int_t fThreshold;
- Int_t fCurrentEtaSlice;
- AliL3Histogram *fCurrentHisto; //!
-
- Float_t fGradX;
- Float_t fGradY;
- Float_t *fXPeaks; //!
- Float_t *fYPeaks; //!
- Int_t *fSTARTXPeaks; //!
- Int_t *fSTARTYPeaks; //!
- Int_t *fENDXPeaks; //!
- Int_t *fENDYPeaks; //!
- Int_t *fSTARTETAPeaks; //!
- Int_t *fENDETAPeaks; //!
- Int_t *fWeight; //!
- Int_t fN1PeaksPrevEtaSlice;
- Int_t fN2PeaksPrevEtaSlice;
- Int_t fNPeaks;
- Int_t fNMax;
-
- Char_t fHistoType;
-
-#ifndef no_root
- TNtuple *fNtuppel; //!
-#endif
public:
AliL3HoughMaxFinder();
void FindAbsMaxima();
void FindBigMaxima();
void FindMaxima(Int_t threshold=0);
- void FindAdaptedPeaks(Int_t nkappawindow,Float_t cut_ratio);
+ void FindAdaptedPeaks(Int_t nkappawindow,Float_t cutratio);
//Peak finder for HoughTransformerRow
void FindAdaptedRowPeaks(Int_t kappawindow,Int_t xsize,Int_t ysize);
//More sophisticated peak finders:
void FindPeak(Int_t t1,Double_t t2,Int_t t3);
- void FindPeak1(Int_t y_window=2,Int_t x_bin_sides=1);
- void SortPeaks(struct AxisWindow **a,Int_t first,Int_t last);
- Int_t PeakCompare(struct AxisWindow *a,struct AxisWindow *b);
+ void FindPeak1(Int_t ywindow=2,Int_t xbinsides=1);
+ void SortPeaks(struct AliL3AxisWindow **a,Int_t first,Int_t last);
+ Int_t PeakCompare(struct AliL3AxisWindow *a,struct AliL3AxisWindow *b) const;
//Setters:
void SetGradient(Float_t x,Float_t y) {fGradX=x; fGradY=y;}
void SetEtaSlice(Int_t etaslice) {fCurrentEtaSlice = etaslice;}
//Getters:
- Float_t GetXPeak(Int_t i);
- Float_t GetYPeak(Int_t i);
- Float_t GetXPeakSize(Int_t i);
- Float_t GetYPeakSize(Int_t i);
- Int_t GetWeight(Int_t i);
- Int_t GetStartEta(Int_t i);
- Int_t GetEndEta(Int_t i);
- Int_t GetEntries() {return fNPeaks;}
+ Float_t GetXPeak(Int_t i) const;
+ Float_t GetYPeak(Int_t i) const;
+ Float_t GetXPeakSize(Int_t i) const;
+ Float_t GetYPeakSize(Int_t i) const;
+ Int_t GetWeight(Int_t i) const;
+ Int_t GetStartEta(Int_t i) const;
+ Int_t GetEndEta(Int_t i) const;
+ Int_t GetEntries() const {return fNPeaks;}
+
+ private:
+
+ Int_t fThreshold; // Threshold for Peak Finder
+ Int_t fCurrentEtaSlice; // Current eta slice being processed
+ AliL3Histogram *fCurrentHisto; //!
+
+ Float_t fGradX; // Gradient threshold inside Peak Finder
+ Float_t fGradY; // Gradient threshold inside Peak Finder
+ Float_t *fXPeaks; //!
+ Float_t *fYPeaks; //!
+ Int_t *fSTARTXPeaks; //!
+ Int_t *fSTARTYPeaks; //!
+ Int_t *fENDXPeaks; //!
+ Int_t *fENDYPeaks; //!
+ Int_t *fSTARTETAPeaks; //!
+ Int_t *fENDETAPeaks; //!
+ Int_t *fWeight; //!
+ Int_t fN1PeaksPrevEtaSlice; // Index of the first peak in the previous eta slice
+ Int_t fN2PeaksPrevEtaSlice; // Index of the last peak in the previous eta slice
+ Int_t fNPeaks; // Index of the last accumulated peak
+ Int_t fNMax; // Maximum allowed number of peaks
+
+ Char_t fHistoType; // Histogram type
+
+#ifndef no_root
+ TNtuple *fNtuppel; //!
+#endif
ClassDef(AliL3HoughMaxFinder,1) //Maximum finder class
};
-inline Float_t AliL3HoughMaxFinder::GetXPeak(Int_t i)
+inline Float_t AliL3HoughMaxFinder::GetXPeak(Int_t i) const
{
if(i<0 || i>fNMax)
{
return fXPeaks[i];
}
-inline Float_t AliL3HoughMaxFinder::GetYPeak(Int_t i)
+inline Float_t AliL3HoughMaxFinder::GetYPeak(Int_t i) const
{
if(i<0 || i>fNMax)
{
}
-inline Int_t AliL3HoughMaxFinder::GetWeight(Int_t i)
+inline Int_t AliL3HoughMaxFinder::GetWeight(Int_t i) const
{
if(i<0 || i>fNMax)
{
return fWeight[i];
}
-inline Int_t AliL3HoughMaxFinder::GetStartEta(Int_t i)
+inline Int_t AliL3HoughMaxFinder::GetStartEta(Int_t i) const
{
if(i<0 || i>fNMax)
{
return fSTARTETAPeaks[i];
}
-inline Int_t AliL3HoughMaxFinder::GetEndEta(Int_t i)
+inline Int_t AliL3HoughMaxFinder::GetEndEta(Int_t i) const
{
if(i<0 || i>fNMax)
{
#endif
}
-AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoEtaOverlap,Bool_t /*DoMC*/) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
+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;
AliL3HoughTransformer::~AliL3HoughTransformer()
{
+ // Dtor
DeleteHistograms();
#ifdef do_mc
if(fTrackID)
void AliL3HoughTransformer::DeleteHistograms()
{
+ // Clean up
if(!fParamSpace)
return;
for(Int_t i=0; i<GetNEtaSegments(); i++)
}
}
-void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Float_t pt_min,
+void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Float_t ptmin,
Int_t nybin,Float_t phimin,Float_t phimax)
{
//Create the histograms (parameter space).
//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
- //phi_max = maximum phi0
+ //ptmin = mimium Pt of track (corresponding to maximum kappa)
+ //phimin = mimimum phi0
+ //phimax = maximum phi0
- Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
+ 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()];
}
}
-inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t eta_index)
+inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t etaindex)
{
- if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
+ // Return a pointer to the histogram which contains etaindex eta slice
+ if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
return 0;
- if(!fParamSpace[eta_index])
+ if(!fParamSpace[etaindex])
return 0;
- return fParamSpace[eta_index];
+ return fParamSpace[etaindex];
}
-Double_t AliL3HoughTransformer::GetEta(Int_t eta_index,Int_t /*slice*/)
+Double_t AliL3HoughTransformer::GetEta(Int_t etaindex,Int_t /*slice*/)
{
- Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
+ // 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 = eta_index + 1;
- eta=(Double_t)((index)*eta_slice);
+ Int_t index = etaindex + 1;
+ eta=(Double_t)((index)*etaslice);
}
else
- eta=(Double_t)((eta_index+0.5)*eta_slice);
+ eta=(Double_t)((etaindex+0.5)*etaslice);
return eta;
}
//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)
+ //kappa = 2/r*sin(phi - phi0)
//
- //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
+ //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
Double_t eta = AliL3Transform::GetEta(xyz);
//Get the corresponding index, which determines which histogram to fill:
- Int_t eta_index = GetEtaIndex(eta);
+ Int_t etaindex = GetEtaIndex(eta);
- if(eta_index < 0 || eta_index >= GetNEtaSegments())
+ if(etaindex < 0 || etaindex >= GetNEtaSegments())
continue;
//Get the correct histogrampointer:
- AliL3Histogram *hist = fParamSpace[eta_index];
+ AliL3Histogram *hist = fParamSpace[etaindex];
if(!hist)
{
- cerr<<"AliL3HoughTransformer::TransformCircle : Error getting histogram in index "<<eta_index<<endl;
+ 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 r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
Float_t phi = AliL3Transform::GetPhi(xyz);
for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
{
Float_t phi0 = hist->GetBinCenterY(b);
- Float_t kappa = 2*sin(phi - phi0)/R;
+ 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);
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)
+ 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[eta_index][bin].fLabel[c] = label;
- fTrackID[eta_index][bin].fNHits[c]++;
+ fTrackID[etaindex][bin].fLabel[c] = label;
+ fTrackID[etaindex][bin].fNHits[c]++;
}
}
#endif
}
}
-struct Digit {
- Int_t row;
- Double_t r;
- Double_t phi;
- Int_t charge;
- Digit *next;
+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
};
-struct EtaContainer {
- Digit *first;
- Digit *last;
+struct AliL3EtaContainer {
+ AliL3Digit *fFirst; //First digit
+ AliL3Digit *fLast; //Last digit
};
-void AliL3HoughTransformer::TransformCircleC(Int_t *row_range,Int_t every)
+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 minrow = AliL3Transform::GetFirstRow(GetPatch());
Int_t maxrow = AliL3Transform::GetLastRow(GetPatch());
- if(row_range)
+ if(rowrange)
{
- minrow = row_range[0];
- maxrow = row_range[1];
+ 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()))
}
Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
- EtaContainer *etaPt = new EtaContainer[bound];
- memset(etaPt,0,bound*sizeof(EtaContainer));
+ AliL3EtaContainer *etaPt = new AliL3EtaContainer[bound];
+ memset(etaPt,0,bound*sizeof(AliL3EtaContainer));
- Digit *digits = new Digit[counter];
- cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
- memset(digits,0,counter*sizeof(Digit));
+ AliL3Digit *digits = new AliL3Digit[counter];
+ cout<<"Allocating "<<counter*sizeof(AliL3Digit)<<" bytes to digitsarray"<<endl;
+ memset(digits,0,counter*sizeof(AliL3Digit));
- Int_t sector,row,tot_charge,pad,time,charge;
- Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
+ Int_t sector,row,totcharge,pad,time,charge;
+ Double_t r1,r2,phi1,phi2,eta,kappa,phi0;
Float_t xyz[3];
counter=0;
AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
eta = AliL3Transform::GetEta(xyz);
- digits[counter].row = i;
- digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
- digits[counter].phi = atan2(xyz[1],xyz[0]);
- digits[counter].charge = charge;
+ 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)
{
- Int_t eta_index = GetEtaIndex(eta);
+ Int_t etaindex = GetEtaIndex(eta);
- Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index;
+ Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
if(index > 0 && index < bound)
{
- if(etaPt[index].first == 0)
- etaPt[index].first = &digits[counter];
+ if(etaPt[index].fFirst == 0)
+ etaPt[index].fFirst = &digits[counter];
else
- (etaPt[index].last)->next = &digits[counter];
- etaPt[index].last = &digits[counter];
+ (etaPt[index].fLast)->fNext = &digits[counter];
+ etaPt[index].fLast = &digits[counter];
}
}
else
{
- Int_t eta_index[2];
- GetEtaIndexes(eta,eta_index);
+ Int_t etaindex[2];
+ GetEtaIndexes(eta,etaindex);
Int_t index[2];
- index[0] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index[0];
- index[1] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index[1];
+ 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;
Int_t ind = index[0];
if(ind > 0 && ind < bound)
{
- if(etaPt[ind].first == 0)
- etaPt[ind].first = &digits[counter];
+ if(etaPt[ind].fFirst == 0)
+ etaPt[ind].fFirst = &digits[counter];
else
- (etaPt[ind].last)->next = &digits[counter];
- etaPt[ind].last = &digits[counter];
+ (etaPt[ind].fLast)->fNext = &digits[counter];
+ etaPt[ind].fLast = &digits[counter];
}
ind = index[1];
if(ind > 0 && ind < bound)
{
- if(etaPt[ind].first == 0)
- etaPt[ind].first = &digits[counter];
+ if(etaPt[ind].fFirst == 0)
+ etaPt[ind].fFirst = &digits[counter];
else
- (etaPt[ind].last)->next = &digits[counter];
- etaPt[ind].last = &digits[counter];
+ (etaPt[ind].fLast)->fNext = &digits[counter];
+ etaPt[ind].fLast = &digits[counter];
}
}
cout<<"Doing the combinatorics"<<endl;
- Digit *dPt1,*dPt2;
+ AliL3Digit *dPt1,*dPt2;
for(Int_t e=0; e<GetNEtaSegments(); e++)
{
{
Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
- for(dPt1 = (Digit*)etaPt[index1].first; dPt1 != 0; dPt1 = (Digit*)dPt1->next)
+ for(dPt1 = (AliL3Digit*)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3Digit*)dPt1->fNext)
{
for(Int_t j=i+every; j<=maxrow; j+=every)
{
Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
- for(dPt2 = (Digit*)etaPt[index2].first; dPt2 != 0; dPt2 = (Digit*)dPt2->next)
+ for(dPt2 = (AliL3Digit*)etaPt[index2].fFirst; dPt2 != 0; dPt2 = (AliL3Digit*)dPt2->fNext)
{
- if(dPt1->row == dPt2->row)
+ if(dPt1->fRow == dPt2->fRow)
{
cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
exit(5);
}
//Do the transform:
- r1 = dPt1->r;
- phi1 = dPt1->phi;
- r2 = dPt2->r;
- phi2 = dPt2->phi;
- phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
- kappa = 2*sin(phi2-phi_0)/r2;
- tot_charge = dPt1->charge + dPt2->charge;
- hist->Fill(kappa,phi_0,tot_charge);
+ 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);
}
}
}
-void AliL3HoughTransformer::TransformLine(Int_t *row_range,Float_t *phirange)
+void AliL3HoughTransformer::TransformLine(Int_t *rowrange,Float_t *phirange)
{
//Do a line transform on the data.
Int_t minrow = AliL3Transform::GetFirstRow(GetPatch());
Int_t maxrow = AliL3Transform::GetLastRow(GetPatch());
- if(row_range)
+ if(rowrange)
{
- minrow = row_range[0];
- maxrow = row_range[1];
+ 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()))
continue;
}
Float_t eta = AliL3Transform::GetEta(xyz);
- Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
- if(eta_index < 0 || eta_index >= GetNEtaSegments())
+ Int_t etaindex = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
+ if(etaindex < 0 || etaindex >= GetNEtaSegments())
continue;
xyz[0] = xyz[0] - AliL3Transform::Row2X(minrow);
//Get the correct histogram:
- AliL3Histogram *hist = fParamSpace[eta_index];
+ AliL3Histogram *hist = fParamSpace[etaindex];
if(!hist)
{
- printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
+ printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",etaindex);
continue;
}
for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
}
-struct LDigit {
- Int_t row;
- Int_t charge;
- Float_t y;
- LDigit *next;
+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 LEtaContainer {
- LDigit *first;
- LDigit *last;
+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")
}
Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
- LEtaContainer *etaPt = new LEtaContainer[bound];
- memset(etaPt,0,bound*sizeof(LEtaContainer));
+ AliL3LEtaContainer *etaPt = new AliL3LEtaContainer[bound];
+ memset(etaPt,0,bound*sizeof(AliL3LEtaContainer));
- LDigit *digits = new LDigit[counter];
- cout<<"Allocating "<<counter*sizeof(LDigit)<<" bytes to digitsarray"<<endl;
- memset(digits,0,counter*sizeof(LDigit));
+ AliL3LDigit *digits = new AliL3LDigit[counter];
+ cout<<"Allocating "<<counter*sizeof(AliL3LDigit)<<" bytes to digitsarray"<<endl;
+ memset(digits,0,counter*sizeof(AliL3LDigit));
Int_t sector,row;
Float_t xyz[3];
Float_t phi = atan2(xyz[1],xyz[0]);
if(phi < phirange[0] || phi > phirange[1]) continue;
- digits[counter].row = i;
- digits[counter].y = xyz[1];
- digits[counter].charge = charge;
+ digits[counter].fRow = i;
+ digits[counter].fY = xyz[1];
+ digits[counter].fCharge = charge;
- Int_t eta_index = GetEtaIndex(eta);
- Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index;
+ Int_t etaindex = GetEtaIndex(eta);
+ Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
if(index > 0 && index < bound)
{
- if(etaPt[index].first == 0)
- etaPt[index].first = &digits[counter];
+ if(etaPt[index].fFirst == 0)
+ etaPt[index].fFirst = &digits[counter];
else
- (etaPt[index].last)->next = &digits[counter];
- etaPt[index].last = &digits[counter];
+ (etaPt[index].fLast)->fNext = &digits[counter];
+ etaPt[index].fLast = &digits[counter];
}
counter++;
}
cout<<"Doing the combinatorics"<<endl;
- LDigit *dPt1,*dPt2;
+ AliL3LDigit *dPt1,*dPt2;
for(Int_t e=0; e<GetNEtaSegments(); e++)
{
{
Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
- for(dPt1 = (LDigit*)etaPt[index1].first; dPt1 != 0; dPt1 = (LDigit*)dPt1->next)
+ for(dPt1 = (AliL3LDigit*)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3LDigit*)dPt1->fNext)
{
for(Int_t j=i+1; j<=rowrange[1]; j++)
{
Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
- for(dPt2 = (LDigit*)etaPt[index2].first; dPt2 != 0; dPt2 = (LDigit*)dPt2->next)
+ for(dPt2 = (AliL3LDigit*)etaPt[index2].fFirst; dPt2 != 0; dPt2 = (AliL3LDigit*)dPt2->fNext)
{
- if(dPt1->row == dPt2->row)
+ if(dPt1->fRow == dPt2->fRow)
{
cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
exit(5);
}
//Do the transform:
- float x1 = AliL3Transform::Row2X(dPt1->row) - AliL3Transform::Row2X(rowrange[0]);
- float x2 = AliL3Transform::Row2X(dPt2->row) - AliL3Transform::Row2X(rowrange[0]);
- float y1 = dPt1->y;
- float y2 = dPt2->y;
+ 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);
delete [] digits;
}
-Int_t AliL3HoughTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
+Int_t AliL3HoughTransformer::GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi)
{
+ // Returns the MC label for a given peak found in the Hough space
if(!fDoMC)
{
cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
}
#ifdef do_mc
- if(eta_index < 0 || eta_index > GetNEtaSegments())
+ if(etaindex < 0 || etaindex > GetNEtaSegments())
{
- cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
+ cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<etaindex<<endl;
return -1;
}
- AliL3Histogram *hist = fParamSpace[eta_index];
+ 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[eta_index][bin].fNHits[i];
+ Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
if(nhits == 0) break;
if(nhits > max)
{
max = nhits;
- label = fTrackID[eta_index][bin].fLabel[i];
+ label = fTrackID[etaindex][bin].fLabel[i];
}
}
//nhits = max;
// @(#) $Id$
-#ifndef ALIL3_HOUGHTRANSFORMER
-#define ALIL3_HOUGHTRANSFORMER
+#ifndef ALIL3HOUGHTRANSFORMER_H
+#define ALIL3HOUGHTRANSFORMER_H
#include "AliL3RootTypes.h"
#include "AliL3HoughBaseTransformer.h"
class AliL3Histogram;
class AliL3HoughTransformer : public AliL3HoughBaseTransformer {
-
- private:
-
- AliL3Histogram **fParamSpace; //!
-#ifdef do_mc
- TrackIndex **fTrackID; //!
-#endif
- Bool_t fDoMC;
- Bool_t fEtaOverlap;
-
- void DeleteHistograms();
public:
AliL3HoughTransformer();
- AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoEtaOverlap=kFALSE,Bool_t DoMC=kFALSE);
+ AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t netasegments,Bool_t DoEtaOverlap=kFALSE,Bool_t DoMC=kFALSE);
virtual ~AliL3HoughTransformer();
void CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t ptres,Int_t nybin,Float_t psi);
Int_t nybin,Float_t ymin,Float_t ymax);
void Reset();
void TransformCircle();
- void TransformCircleC(Int_t *row_range,Int_t every=1);
+ void TransformCircleC(Int_t *rowrange,Int_t every=1);
void TransformLine(Int_t *rowrange=0,Float_t *phirange=0);
void TransformLineC(Int_t *rowrange,Float_t *phirange);
Int_t GetEtaIndex(Double_t eta);
void GetEtaIndexes(Double_t eta,Int_t *indexes);
- AliL3Histogram *GetHistogram(Int_t eta_index);
- Double_t GetEta(Int_t eta_index,Int_t slice);
- Int_t GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi);
+ AliL3Histogram *GetHistogram(Int_t etaindex);
+ Double_t GetEta(Int_t etaindex,Int_t slice);
+ Int_t GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi);
+ private:
+
+ AliL3Histogram **fParamSpace; //!
+#ifdef do_mc
+ TrackIndex **fTrackID; //!
+#endif
+ Bool_t fDoMC; // Calculate mc labels or not
+ Bool_t fEtaOverlap; // Allow overlapping of eta slice or not
+
+ void DeleteHistograms();
+
ClassDef(AliL3HoughTransformer,1) //Normal Hough transformation class
};
ClassImp(AliL3HoughTransformerRow)
-UChar_t **AliL3HoughTransformerRow::fRowCount = 0;
-UChar_t **AliL3HoughTransformerRow::fGapCount = 0;
-UChar_t **AliL3HoughTransformerRow::fCurrentRowCount = 0;
+UChar_t **AliL3HoughTransformerRow::fgRowCount = 0;
+UChar_t **AliL3HoughTransformerRow::fgGapCount = 0;
+UChar_t **AliL3HoughTransformerRow::fgCurrentRowCount = 0;
#ifdef do_mc
-TrackIndex **AliL3HoughTransformerRow::fTrackID = 0;
+TrackIndex **AliL3HoughTransformerRow::fgTrackID = 0;
#endif
-UChar_t *AliL3HoughTransformerRow::fTrackNRows = 0;
-UChar_t *AliL3HoughTransformerRow::fTrackFirstRow = 0;
-UChar_t *AliL3HoughTransformerRow::fTrackLastRow = 0;
+UChar_t *AliL3HoughTransformerRow::fgTrackNRows = 0;
+UChar_t *AliL3HoughTransformerRow::fgTrackFirstRow = 0;
+UChar_t *AliL3HoughTransformerRow::fgTrackLastRow = 0;
-Float_t AliL3HoughTransformerRow::fBeta1 = AliL3Transform::Row2X(79)/pow(AliL3Transform::Row2X(79),2);
-Float_t AliL3HoughTransformerRow::fBeta2 = (AliL3Transform::Row2X(158)+6.0)/pow((AliL3Transform::Row2X(158)+6.0),2);
+Float_t AliL3HoughTransformerRow::fgBeta1 = AliL3Transform::Row2X(79)/pow(AliL3Transform::Row2X(79),2);
+Float_t AliL3HoughTransformerRow::fgBeta2 = (AliL3Transform::Row2X(158)+6.0)/pow((AliL3Transform::Row2X(158)+6.0),2);
AliL3HoughTransformerRow::AliL3HoughTransformerRow()
{
fLUTbackwardZ2=0;
}
-AliL3HoughTransformerRow::AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t /*DoMC*/,Float_t zvertex) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments,zvertex)
+AliL3HoughTransformerRow::AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t netasegments,Bool_t /*DoMC*/,Float_t zvertex) : AliL3HoughBaseTransformer(slice,patch,netasegments,zvertex)
{
//Normal constructor
fParamSpace = 0;
AliL3HoughTransformerRow::~AliL3HoughTransformerRow()
{
+ //Destructor
DeleteHistograms();
#ifdef do_mc
- if(fTrackID)
+ if(fgTrackID)
{
for(Int_t i=0; i<GetNEtaSegments(); i++)
{
- if(!fTrackID[i]) continue;
- delete fTrackID[i];
+ if(!fgTrackID[i]) continue;
+ delete fgTrackID[i];
}
- delete [] fTrackID;
- fTrackID = 0;
+ delete [] fgTrackID;
+ fgTrackID = 0;
}
#endif
- if(fRowCount)
+ if(fgRowCount)
{
for(Int_t i=0; i<GetNEtaSegments(); i++)
{
- if(!fRowCount[i]) continue;
- delete [] fRowCount[i];
+ if(!fgRowCount[i]) continue;
+ delete [] fgRowCount[i];
}
- delete [] fRowCount;
- fRowCount = 0;
+ delete [] fgRowCount;
+ fgRowCount = 0;
}
- if(fGapCount)
+ if(fgGapCount)
{
for(Int_t i=0; i<GetNEtaSegments(); i++)
{
- if(!fGapCount[i]) continue;
- delete [] fGapCount[i];
+ if(!fgGapCount[i]) continue;
+ delete [] fgGapCount[i];
}
- delete [] fGapCount;
- fGapCount = 0;
+ delete [] fgGapCount;
+ fgGapCount = 0;
}
- if(fCurrentRowCount)
+ if(fgCurrentRowCount)
{
for(Int_t i=0; i<GetNEtaSegments(); i++)
{
- if(fCurrentRowCount[i])
- delete [] fCurrentRowCount[i];
+ if(fgCurrentRowCount[i])
+ delete [] fgCurrentRowCount[i];
}
- delete [] fCurrentRowCount;
- fCurrentRowCount = 0;
+ delete [] fgCurrentRowCount;
+ fgCurrentRowCount = 0;
}
- if(fTrackNRows)
+ if(fgTrackNRows)
{
- delete [] fTrackNRows;
- fTrackNRows = 0;
+ delete [] fgTrackNRows;
+ fgTrackNRows = 0;
}
- if(fTrackFirstRow)
+ if(fgTrackFirstRow)
{
- delete [] fTrackFirstRow;
- fTrackFirstRow = 0;
+ delete [] fgTrackFirstRow;
+ fgTrackFirstRow = 0;
}
- if(fTrackLastRow)
+ if(fgTrackLastRow)
{
- delete [] fTrackLastRow;
- fTrackLastRow = 0;
+ delete [] fgTrackLastRow;
+ fgTrackLastRow = 0;
}
}
void AliL3HoughTransformerRow::DeleteHistograms()
{
+ // Clean up
delete[] fLUTforwardZ;
delete[] fLUTforwardZ2;
delete[] fLUTbackwardZ;
delete [] fParamSpace;
}
-void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t pt_min,
+void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t ptmin,
Int_t nybin,Float_t phimin,Float_t phimax)
{
//Create the histograms (parameter space).
//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
//phi_max = maximum phi0
- Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
+ 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 AliL3HoughTransformerRow::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];
Int_t ncellsx = (hist->GetNbinsX()+3)/2;
Int_t ncellsy = (hist->GetNbinsY()+3)/2;
Int_t ncells = ncellsx*ncellsy;
- if(!fTrackID)
+ if(!fgTrackID)
{
LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<ENDLOG;
- fTrackID = new TrackIndex*[GetNEtaSegments()];
+ <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fgTrackID"<<ENDLOG;
+ fgTrackID = new TrackIndex*[GetNEtaSegments()];
for(Int_t i=0; i<GetNEtaSegments(); i++)
- fTrackID[i] = new TrackIndex[ncells];
+ fgTrackID[i] = new TrackIndex[ncells];
}
}
#endif
AliL3Histogram *hist = fParamSpace[0];
Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
- if(!fRowCount)
+ if(!fgRowCount)
{
LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fRowCount"<<ENDLOG;
- fRowCount = new UChar_t*[GetNEtaSegments()];
+ <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgRowCount"<<ENDLOG;
+ fgRowCount = new UChar_t*[GetNEtaSegments()];
for(Int_t i=0; i<GetNEtaSegments(); i++)
- fRowCount[i] = new UChar_t[ncells];
+ fgRowCount[i] = new UChar_t[ncells];
}
- if(!fGapCount)
+ if(!fgGapCount)
{
LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fGapCount"<<ENDLOG;
- fGapCount = new UChar_t*[GetNEtaSegments()];
+ <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgGapCount"<<ENDLOG;
+ fgGapCount = new UChar_t*[GetNEtaSegments()];
for(Int_t i=0; i<GetNEtaSegments(); i++)
- fGapCount[i] = new UChar_t[ncells];
+ fgGapCount[i] = new UChar_t[ncells];
}
- if(!fCurrentRowCount)
+ if(!fgCurrentRowCount)
{
LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fCurrentRowCount"<<ENDLOG;
- fCurrentRowCount = new UChar_t*[GetNEtaSegments()];
+ <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgCurrentRowCount"<<ENDLOG;
+ fgCurrentRowCount = new UChar_t*[GetNEtaSegments()];
for(Int_t i=0; i<GetNEtaSegments(); i++)
- fCurrentRowCount[i] = new UChar_t[ncells];
+ fgCurrentRowCount[i] = new UChar_t[ncells];
}
- if(!fTrackNRows)
+ if(!fgTrackNRows)
{
LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackNRows"<<ENDLOG;
- fTrackNRows = new UChar_t[ncells];
+ <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackNRows"<<ENDLOG;
+ fgTrackNRows = new UChar_t[ncells];
LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackFirstRow"<<ENDLOG;
- fTrackFirstRow = new UChar_t[ncells];
+ <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackFirstRow"<<ENDLOG;
+ fgTrackFirstRow = new UChar_t[ncells];
LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackLastRow"<<ENDLOG;
- fTrackLastRow = new UChar_t[ncells];
+ <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackLastRow"<<ENDLOG;
+ fgTrackLastRow = new UChar_t[ncells];
AliL3HoughTrack track;
Int_t xmin = hist->GetFirstXbin();
{
//cvetan: we get strange warning on gcc-2.95
//warning: large integer implicitly truncated to unsigned type
- fTrackNRows[xbin + ybin*nxbins] = 99999;
+ fgTrackNRows[xbin + ybin*nxbins] = 99999;
for(Int_t deltay = -1; deltay <= 1; deltay += 2) {
for(Int_t deltax = -1; deltax <= 1; deltax += 2) {
Float_t xtrack = hist->GetPreciseBinCenterX((Float_t)xbin+0.5*(Float_t)deltax);
Float_t ytrack = hist->GetPreciseBinCenterY((Float_t)ybin+0.5*(Float_t)deltay);
- Float_t psi = atan((xtrack-ytrack)/(fBeta1-fBeta2));
- Float_t kappa = 2.0*(xtrack*cos(psi)-fBeta1*sin(psi));
+ Float_t psi = atan((xtrack-ytrack)/(fgBeta1-fgBeta2));
+ Float_t kappa = 2.0*(xtrack*cos(psi)-fgBeta1*sin(psi));
track.SetTrackParameters(kappa,psi,1);
- Bool_t first_row = kFALSE;
+ Bool_t firstrow = kFALSE;
UInt_t maxfirstrow = 0;
UInt_t maxlastrow = 0;
UInt_t curfirstrow = 0;
AliL3Transform::LocHLT2Raw(hit,0,j);
if(hit[1]>=0 && hit[1]<AliL3Transform::GetNPads(j))
{
- if(!first_row) {
+ if(!firstrow) {
curfirstrow = j;
- first_row = kTRUE;
+ firstrow = kTRUE;
}
curlastrow = j;
}
else {
- if(first_row) {
- first_row = kFALSE;
+ if(firstrow) {
+ firstrow = kFALSE;
if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
maxfirstrow = curfirstrow;
maxlastrow = curlastrow;
maxfirstrow = curfirstrow;
maxlastrow = curlastrow;
}
- if((maxlastrow-maxfirstrow) < fTrackNRows[xbin + ybin*nxbins]) {
- fTrackNRows[xbin + ybin*nxbins] = maxlastrow-maxfirstrow;
- fTrackFirstRow[xbin + ybin*nxbins] = maxfirstrow;
- fTrackLastRow[xbin + ybin*nxbins] = maxlastrow;
+ if((maxlastrow-maxfirstrow) < fgTrackNRows[xbin + ybin*nxbins]) {
+ fgTrackNRows[xbin + ybin*nxbins] = maxlastrow-maxfirstrow;
+ fgTrackFirstRow[xbin + ybin*nxbins] = maxfirstrow;
+ fgTrackLastRow[xbin + ybin*nxbins] = maxlastrow;
}
}
}
- // cout<<" fTrackNRows "<<xbin<<" "<<ybin<<" "<<(Int_t)fTrackNRows[xbin + ybin*nxbins]<<" "<<endl;
+ // cout<<" fgTrackNRows "<<xbin<<" "<<ybin<<" "<<(Int_t)fgTrackNRows[xbin + ybin*nxbins]<<" "<<endl;
}
}
}
void AliL3HoughTransformerRow::Reset()
{
- //Reset all the histograms
+ //Reset all the histograms. Should be done when processing new slice
if(!fParamSpace)
{
Int_t ncellsy = (hist->GetNbinsY()+3)/2;
Int_t ncells = ncellsx*ncellsy;
for(Int_t i=0; i<GetNEtaSegments(); i++)
- memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
+ memset(fgTrackID[i],0,ncells*sizeof(TrackIndex));
}
#endif
AliL3Histogram *hist = fParamSpace[0];
Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
for(Int_t i=0; i<GetNEtaSegments(); i++)
{
- memset(fRowCount[i],0,ncells*sizeof(UChar_t));
- memset(fGapCount[i],1,ncells*sizeof(UChar_t));
- // memset(fCurrentRowCount[i],160,ncells*sizeof(UChar_t));
- memcpy(fCurrentRowCount[i],fTrackFirstRow,ncells*sizeof(UChar_t));
+ memset(fgRowCount[i],0,ncells*sizeof(UChar_t));
+ memset(fgGapCount[i],1,ncells*sizeof(UChar_t));
+ // memset(fgCurrentRowCount[i],160,ncells*sizeof(UChar_t));
+ memcpy(fgCurrentRowCount[i],fgTrackFirstRow,ncells*sizeof(UChar_t));
}
}
return (Int_t)index;
}
-inline AliL3Histogram *AliL3HoughTransformerRow::GetHistogram(Int_t eta_index)
+inline AliL3Histogram *AliL3HoughTransformerRow::GetHistogram(Int_t etaindex)
{
- if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
+ // Return a pointer to the histogram which contains etaindex eta slice
+ if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
return 0;
- if(!fParamSpace[eta_index])
+ if(!fParamSpace[etaindex])
return 0;
- return fParamSpace[eta_index];
+ return fParamSpace[etaindex];
}
-Double_t AliL3HoughTransformerRow::GetEta(Int_t eta_index,Int_t /*slice*/)
+Double_t AliL3HoughTransformerRow::GetEta(Int_t etaindex,Int_t /*slice*/)
{
- Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
+ // Return eta calculated in the middle of the eta slice
+ Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
Double_t eta=0;
- eta=(Double_t)((eta_index+0.5)*eta_slice);
+ eta=(Double_t)((etaindex+0.5)*etaslice);
return eta;
}
-struct EtaRow {
- UChar_t start_pad;
- UChar_t end_pad;
- Bool_t found;
- Float_t start_y;
+struct AliL3EtaRow {
+ UChar_t fStartPad; //First pad in the cluster
+ UChar_t fEndPad; //Last pad in the cluster
+ Bool_t fIsFound; //Is the cluster already found
+ Float_t fStartY; //Y position of the first pad in the cluster
#ifdef do_mc
- Int_t mc_labels[MaxTrack];
+ Int_t fMcLabels[MaxTrack]; //Array to store mc labels inside cluster
#endif
};
void AliL3HoughTransformerRow::TransformCircle()
{
- Float_t beta1 = fBeta1;
- Float_t beta2 = fBeta2;
- Float_t beta1_minus_beta2 = fBeta1 - fBeta2;
+ //Do the Hough Transform
+ Float_t beta1 = fgBeta1;
+ Float_t beta2 = fgBeta2;
+ Float_t beta1minusbeta2 = fgBeta1 - fgBeta2;
- Int_t n_eta_segments = GetNEtaSegments();
- Double_t eta_min = GetEtaMin();
- Double_t eta_slice = (GetEtaMax() - eta_min)/n_eta_segments;
+ Int_t netasegments = GetNEtaSegments();
+ Double_t etamin = GetEtaMin();
+ Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
- Int_t lower_threshold = GetLowerThreshold();
+ Int_t lowerthreshold = GetLowerThreshold();
//Assumes that all the etaslice histos are the same!
AliL3Histogram *h = fParamSpace[0];
- Float_t y_min = h->GetYmin();
+ Float_t ymin = h->GetYmin();
//Float_t y_max = h->GetYmax();
- Float_t hist_bin = h->GetBinWidthY();
- Int_t first_bin = h->GetFirstYbin();
- Int_t last_bin = h->GetLastYbin();
- Float_t x_min = h->GetXmin();
- Float_t x_max = h->GetXmax();
- Float_t x_bin = (x_max-x_min)/h->GetNbinsX();
- Int_t first_binx = h->GetFirstXbin()-1;
- Int_t last_binx = h->GetLastXbin()+1;
+ Float_t histbin = h->GetBinWidthY();
+ Int_t firstbin = h->GetFirstYbin();
+ Int_t lastbin = h->GetLastYbin();
+ Float_t xmin = h->GetXmin();
+ Float_t xmax = h->GetXmax();
+ Float_t xbin = (xmax-xmin)/h->GetNbinsX();
+ Int_t firstbinx = h->GetFirstXbin()-1;
+ Int_t lastbinx = h->GetLastXbin()+1;
Int_t nbinx = h->GetNbinsX()+2;
- UChar_t last_pad;
- Int_t last_eta_index;
- EtaRow *eta_clust = new EtaRow[n_eta_segments];
+ UChar_t lastpad;
+ Int_t lastetaindex;
+ AliL3EtaRow *etaclust = new AliL3EtaRow[netasegments];
AliL3DigitRowData *tempPt = GetDataPointer();
if(!tempPt)
}
Int_t ipatch = GetPatch();
- Double_t pad_pitch = AliL3Transform::GetPadPitchWidth(ipatch);
+ Double_t padpitch = AliL3Transform::GetPadPitchWidth(ipatch);
Int_t islice = GetSlice();
Float_t *fLUTz;
Float_t *fLUTz2;
{
Int_t npads = AliL3Transform::GetNPads((Int_t)i)-1;
- last_pad = 255;
+ lastpad = 255;
//Flush eta clusters array
- memset(eta_clust,0,n_eta_segments*sizeof(EtaRow));
+ memset(etaclust,0,netasegments*sizeof(AliL3EtaRow));
Float_t x = AliL3Transform::Row2X((Int_t)i);
Float_t x2 = x*x;
for(UInt_t j=0; j<tempPt->fNDigit; j++)
{
UShort_t charge = digPt[j].fCharge;
- if((Int_t)charge <= lower_threshold)
+ if((Int_t)charge <= lowerthreshold)
continue;
UChar_t pad = digPt[j].fPad;
UShort_t time = digPt[j].fTime;
- if(pad != last_pad)
+ if(pad != lastpad)
{
- y = (pad-0.5*npads)*pad_pitch;
+ y = (pad-0.5*npads)*padpitch;
Float_t y2 = y*y;
r2 = x2 + y2;
- last_eta_index = -1;
+ lastetaindex = -1;
}
//Transform data to local cartesian coordinates:
Double_t r = sqrt(r2+z2);
Double_t eta = 0.5 * log((r+z)/(r-z));
//Get the corresponding index, which determines which histogram to fill:
- Int_t eta_index = (Int_t)((eta-eta_min)/eta_slice);
+ Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
#ifndef do_mc
- if(eta_index == last_eta_index) continue;
+ if(etaindex == lastetaindex) continue;
#endif
- // cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<eta_index<<endl;
+ // cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
- if(eta_index < 0 || eta_index >= n_eta_segments)
+ if(etaindex < 0 || etaindex >= netasegments)
continue;
- if(!eta_clust[eta_index].found)
+ if(!etaclust[etaindex].fIsFound)
{
- eta_clust[eta_index].start_pad = pad;
- eta_clust[eta_index].end_pad = pad;
- eta_clust[eta_index].start_y = y - pad_pitch/2.0;
- eta_clust[eta_index].found = 1;
+ etaclust[etaindex].fStartPad = pad;
+ etaclust[etaindex].fEndPad = pad;
+ etaclust[etaindex].fStartY = y - padpitch/2.0;
+ etaclust[etaindex].fIsFound = 1;
#ifdef do_mc
if(fDoMC)
{
if(label < 0) break;
UInt_t c;
for(c=0; c<(MaxTrack-1); c++)
- if(eta_clust[eta_index].mc_labels[c] == label || eta_clust[eta_index].mc_labels[c] == 0)
+ if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
break;
// if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
- eta_clust[eta_index].mc_labels[c] = label;
+ etaclust[etaindex].fMcLabels[c] = label;
}
}
#endif
}
else
{
- if(pad <= (eta_clust[eta_index].end_pad + 1))
+ if(pad <= (etaclust[etaindex].fEndPad + 1))
{
- eta_clust[eta_index].end_pad = pad;
+ etaclust[etaindex].fEndPad = pad;
#ifdef do_mc
if(fDoMC)
{
if(label < 0) break;
UInt_t c;
for(c=0; c<(MaxTrack-1); c++)
- if(eta_clust[eta_index].mc_labels[c] == label || eta_clust[eta_index].mc_labels[c] == 0)
+ if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
break;
// if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
- eta_clust[eta_index].mc_labels[c] = label;
+ etaclust[etaindex].fMcLabels[c] = label;
}
}
#endif
}
else
{
- Bool_t fill_cluster = kTRUE;
- if(fill_cluster) {
+ Bool_t fillcluster = kTRUE;
+ if(fillcluster) {
- UChar_t *nrows = fRowCount[eta_index];
- UChar_t *ngaps = fGapCount[eta_index];
- UChar_t *currentrow = fCurrentRowCount[eta_index];
- UChar_t *lastrow = fTrackLastRow;
+ UChar_t *nrows = fgRowCount[etaindex];
+ UChar_t *ngaps = fgGapCount[etaindex];
+ UChar_t *currentrow = fgCurrentRowCount[etaindex];
+ UChar_t *lastrow = fgTrackLastRow;
//Do the transformation:
- Float_t start_y = eta_clust[eta_index].start_y;
- if(eta_clust[eta_index].start_pad == 0)
- start_y -= 0.0*pad_pitch;
- Float_t R1 = x2 + start_y*start_y;
- Float_t x_over_R1 = x/R1;
- Float_t start_y_over_R1 = start_y/R1;
- Float_t end_y = (eta_clust[eta_index].end_pad-0.5*(npads-1))*pad_pitch;
- if(eta_clust[eta_index].end_pad == npads)
- end_y += 0.0*pad_pitch;
- Float_t R2 = x2 + end_y*end_y;
- Float_t x_over_R2 = x/R2;
- Float_t end_y_over_R2 = end_y/R2;
- Float_t A1 = beta1_minus_beta2/(x_over_R1-beta2);
- Float_t B1 = (x_over_R1-beta1)/(x_over_R1-beta2);
- Float_t A2 = beta1_minus_beta2/(x_over_R2-beta2);
- Float_t B2 = (x_over_R2-beta1)/(x_over_R2-beta2);
-
- Float_t kappa1 = (A1*start_y_over_R1+B1*y_min-x_min)/x_bin;
- Float_t delta_kappa1 = B1*hist_bin/x_bin;
- if(B1<0)
- kappa1 += delta_kappa1;
- Float_t kappa2 = (A2*end_y_over_R2+B2*y_min-x_min)/x_bin;
- Float_t delta_kappa2 = B2*hist_bin/x_bin;
- if(B2>=0)
- kappa2 += delta_kappa2;
+ Float_t starty = etaclust[etaindex].fStartY;
+ if(etaclust[etaindex].fStartPad == 0)
+ starty -= 0.0*padpitch;
+ Float_t r1 = x2 + starty*starty;
+ Float_t xoverr1 = x/r1;
+ Float_t startyoverr1 = starty/r1;
+ Float_t endy = (etaclust[etaindex].fEndPad-0.5*(npads-1))*padpitch;
+ if(etaclust[etaindex].fEndPad == npads)
+ endy += 0.0*padpitch;
+ Float_t r2 = x2 + endy*endy;
+ Float_t xoverr2 = x/r2;
+ Float_t endyoverr2 = endy/r2;
+ Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
+ Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
+ Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
+ Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
+
+ Float_t kappa1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
+ Float_t deltaalpha1 = b1*histbin/xbin;
+ if(b1<0)
+ kappa1 += deltaalpha1;
+ Float_t kappa2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
+ Float_t deltaalpha2 = b2*histbin/xbin;
+ if(b2>=0)
+ kappa2 += deltaalpha2;
//Fill the histogram along the phirange
- for(Int_t b=first_bin; b<=last_bin; b++, kappa1 += delta_kappa1, kappa2 += delta_kappa2)
+ for(Int_t b=firstbin; b<=lastbin; b++, kappa1 += deltaalpha1, kappa2 += deltaalpha2)
{
Int_t binx1 = 1 + (Int_t)kappa1;
- if(binx1>last_binx) continue;
- if(binx1<first_binx) binx1 = first_binx;
+ if(binx1>lastbinx) continue;
+ if(binx1<firstbinx) binx1 = firstbinx;
Int_t binx2 = 1 + (Int_t)kappa2;
- if(binx2<first_binx) continue;
- if(binx2>last_binx) binx2 = last_binx;
+ if(binx2<firstbinx) continue;
+ if(binx2>lastbinx) binx2 = lastbinx;
#ifdef do_mc
if(binx2<binx1) {
LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
- <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<start_y<<" "<<end_y<<ENDLOG;
+ <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<starty<<" "<<endy<<ENDLOG;
}
#endif
- Int_t temp_bin = b*nbinx;
- UChar_t *nrows2 = nrows + temp_bin;
- UChar_t *ngaps2 = ngaps + temp_bin;
- UChar_t *currentrow2 = currentrow + temp_bin;
- UChar_t *lastrow2 = lastrow + temp_bin;
+ Int_t tempbin = b*nbinx;
+ UChar_t *nrows2 = nrows + tempbin;
+ UChar_t *ngaps2 = ngaps + tempbin;
+ UChar_t *currentrow2 = currentrow + tempbin;
+ UChar_t *lastrow2 = lastrow + tempbin;
for(Int_t bin=binx1;bin<=binx2;bin++)
{
if(ngaps2[bin] < MAX_N_GAPS) {
{
for(UInt_t t=0;t<(MaxTrack-1); t++)
{
- Int_t label = eta_clust[eta_index].mc_labels[t];
+ Int_t label = etaclust[etaindex].fMcLabels[t];
if(label == 0) break;
UInt_t c;
- Int_t temp_bin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
+ Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
for(c=0; c<(MaxTrack-1); c++)
- if(fTrackID[eta_index][temp_bin2].fLabel[c] == label || fTrackID[eta_index][temp_bin2].fNHits[c] == 0)
+ if(fgTrackID[etaindex][tempbin2].fLabel[c] == label || fgTrackID[etaindex][tempbin2].fNHits[c] == 0)
break;
// if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
- fTrackID[eta_index][temp_bin2].fLabel[c] = label;
- if(fTrackID[eta_index][temp_bin2].fCurrentRow[c] != i) {
- fTrackID[eta_index][temp_bin2].fNHits[c]++;
- fTrackID[eta_index][temp_bin2].fCurrentRow[c] = i;
+ fgTrackID[etaindex][tempbin2].fLabel[c] = label;
+ if(fgTrackID[etaindex][tempbin2].fCurrentRow[c] != i) {
+ fgTrackID[etaindex][tempbin2].fNHits[c]++;
+ fgTrackID[etaindex][tempbin2].fCurrentRow[c] = i;
}
}
}
}
- eta_clust[eta_index].start_pad = pad;
- eta_clust[eta_index].end_pad = pad;
- eta_clust[eta_index].start_y = y - pad_pitch/2.0;
+ etaclust[etaindex].fStartPad = pad;
+ etaclust[etaindex].fEndPad = pad;
+ etaclust[etaindex].fStartY = y - padpitch/2.0;
#ifdef do_mc
if(fDoMC)
{
- memset(eta_clust[eta_index].mc_labels,0,MaxTrack);
+ memset(etaclust[etaindex].fMcLabels,0,MaxTrack);
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-1); c++)
- if(eta_clust[eta_index].mc_labels[c] == label || eta_clust[eta_index].mc_labels[c] == 0)
+ if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
break;
// if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
- eta_clust[eta_index].mc_labels[c] = label;
+ etaclust[etaindex].fMcLabels[c] = label;
}
}
#endif
}
}
- last_pad = pad;
- last_eta_index = eta_index;
+ lastpad = pad;
+ lastetaindex = etaindex;
}
//Write remaining clusters
- for(Int_t eta_index = 0;eta_index < n_eta_segments;eta_index++)
+ for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
{
//Check for empty row
- if((eta_clust[eta_index].start_pad == 0) && (eta_clust[eta_index].end_pad == 0)) continue;
+ if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
- UChar_t *nrows = fRowCount[eta_index];
- UChar_t *ngaps = fGapCount[eta_index];
- UChar_t *currentrow = fCurrentRowCount[eta_index];
- UChar_t *lastrow = fTrackLastRow;
+ UChar_t *nrows = fgRowCount[etaindex];
+ UChar_t *ngaps = fgGapCount[etaindex];
+ UChar_t *currentrow = fgCurrentRowCount[etaindex];
+ UChar_t *lastrow = fgTrackLastRow;
//Do the transformation:
- Float_t start_y = eta_clust[eta_index].start_y;
- if(eta_clust[eta_index].start_pad == 0)
- start_y -= 0.0*pad_pitch;
- Float_t R1 = x2 + start_y*start_y;
- Float_t x_over_R1 = x/R1;
- Float_t start_y_over_R1 = start_y/R1;
- Float_t end_y = (eta_clust[eta_index].end_pad-0.5*(npads-1))*pad_pitch;
- if(eta_clust[eta_index].end_pad == npads)
- end_y += 0.0*pad_pitch;
- Float_t R2 = x2 + end_y*end_y;
- Float_t x_over_R2 = x/R2;
- Float_t end_y_over_R2 = end_y/R2;
- Float_t A1 = beta1_minus_beta2/(x_over_R1-beta2);
- Float_t B1 = (x_over_R1-beta1)/(x_over_R1-beta2);
- Float_t A2 = beta1_minus_beta2/(x_over_R2-beta2);
- Float_t B2 = (x_over_R2-beta1)/(x_over_R2-beta2);
-
- Float_t kappa1 = (A1*start_y_over_R1+B1*y_min-x_min)/x_bin;
- Float_t delta_kappa1 = B1*hist_bin/x_bin;
- if(B1<0)
- kappa1 += delta_kappa1;
- Float_t kappa2 = (A2*end_y_over_R2+B2*y_min-x_min)/x_bin;
- Float_t delta_kappa2 = B2*hist_bin/x_bin;
- if(B2>=0)
- kappa2 += delta_kappa2;
+ Float_t starty = etaclust[etaindex].fStartY;
+ if(etaclust[etaindex].fStartPad == 0)
+ starty -= 0.0*padpitch;
+ Float_t r1 = x2 + starty*starty;
+ Float_t xoverr1 = x/r1;
+ Float_t startyoverr1 = starty/r1;
+ Float_t endy = (etaclust[etaindex].fEndPad-0.5*(npads-1))*padpitch;
+ if(etaclust[etaindex].fEndPad == npads)
+ endy += 0.0*padpitch;
+ Float_t r2 = x2 + endy*endy;
+ Float_t xoverr2 = x/r2;
+ Float_t endyoverr2 = endy/r2;
+ Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
+ Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
+ Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
+ Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
+
+ Float_t kappa1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
+ Float_t deltaalpha1 = b1*histbin/xbin;
+ if(b1<0)
+ kappa1 += deltaalpha1;
+ Float_t kappa2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
+ Float_t deltaalpha2 = b2*histbin/xbin;
+ if(b2>=0)
+ kappa2 += deltaalpha2;
//Fill the histogram along the phirange
- for(Int_t b=first_bin; b<=last_bin; b++, kappa1 += delta_kappa1, kappa2 += delta_kappa2)
+ for(Int_t b=firstbin; b<=lastbin; b++, kappa1 += deltaalpha1, kappa2 += deltaalpha2)
{
Int_t binx1 = 1 + (Int_t)kappa1;
- if(binx1>last_binx) continue;
- if(binx1<first_binx) binx1 = first_binx;
+ if(binx1>lastbinx) continue;
+ if(binx1<firstbinx) binx1 = firstbinx;
Int_t binx2 = 1 + (Int_t)kappa2;
- if(binx2<first_binx) continue;
- if(binx2>last_binx) binx2 = last_binx;
+ if(binx2<firstbinx) continue;
+ if(binx2>lastbinx) binx2 = lastbinx;
#ifdef do_mc
if(binx2<binx1) {
LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
- <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<start_y<<" "<<end_y<<ENDLOG;
+ <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<starty<<" "<<endy<<ENDLOG;
}
#endif
- Int_t temp_bin = b*nbinx;
- UChar_t *nrows2 = nrows + temp_bin;
- UChar_t *ngaps2 = ngaps + temp_bin;
- UChar_t *currentrow2 = currentrow + temp_bin;
- UChar_t *lastrow2 = lastrow + temp_bin;
+ Int_t tempbin = b*nbinx;
+ UChar_t *nrows2 = nrows + tempbin;
+ UChar_t *ngaps2 = ngaps + tempbin;
+ UChar_t *currentrow2 = currentrow + tempbin;
+ UChar_t *lastrow2 = lastrow + tempbin;
for(Int_t bin=binx1;bin<=binx2;bin++)
{
if(ngaps2[bin] < MAX_N_GAPS) {
{
for(UInt_t t=0;t<(MaxTrack-1); t++)
{
- Int_t label = eta_clust[eta_index].mc_labels[t];
+ Int_t label = etaclust[etaindex].fMcLabels[t];
if(label == 0) break;
UInt_t c;
- Int_t temp_bin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
+ Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
for(c=0; c<(MaxTrack-1); c++)
- if(fTrackID[eta_index][temp_bin2].fLabel[c] == label || fTrackID[eta_index][temp_bin2].fNHits[c] == 0)
+ if(fgTrackID[etaindex][tempbin2].fLabel[c] == label || fgTrackID[etaindex][tempbin2].fNHits[c] == 0)
break;
// if(c == MaxTrack-1) cout<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
- fTrackID[eta_index][temp_bin2].fLabel[c] = label;
- if(fTrackID[eta_index][temp_bin2].fCurrentRow[c] != i) {
- fTrackID[eta_index][temp_bin2].fNHits[c]++;
- fTrackID[eta_index][temp_bin2].fCurrentRow[c] = i;
+ fgTrackID[etaindex][tempbin2].fLabel[c] = label;
+ if(fgTrackID[etaindex][tempbin2].fCurrentRow[c] != i) {
+ fgTrackID[etaindex][tempbin2].fNHits[c]++;
+ fgTrackID[etaindex][tempbin2].fCurrentRow[c] = i;
}
}
}
AliL3MemHandler::UpdateRowPointer(tempPt);
}
- delete [] eta_clust;
+ delete [] etaclust;
}
-Int_t AliL3HoughTransformerRow::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
+Int_t AliL3HoughTransformerRow::GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi)
{
+ // Returns the MC label for a given peak found in the Hough space
if(!fDoMC)
{
LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
}
#ifdef do_mc
- if(eta_index < 0 || eta_index > GetNEtaSegments())
+ if(etaindex < 0 || etaindex > GetNEtaSegments())
{
LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
- <<"Wrong etaindex "<<eta_index<<ENDLOG;
+ <<"Wrong etaindex "<<etaindex<<ENDLOG;
return -1;
}
- AliL3Histogram *hist = fParamSpace[eta_index];
+ AliL3Histogram *hist = fParamSpace[etaindex];
Int_t bin = hist->FindLabelBin(kappa,psi);
if(bin==-1) {
LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
Int_t max=0;
for(UInt_t i=0; i<(MaxTrack-1); i++)
{
- Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
+ Int_t nhits=fgTrackID[etaindex][bin].fNHits[i];
if(nhits == 0) break;
if(nhits > max)
{
max = nhits;
- label = fTrackID[eta_index][bin].fLabel[i];
+ label = fgTrackID[etaindex][bin].fLabel[i];
}
}
Int_t label2=-1;
Int_t max2=0;
for(UInt_t i=0; i<(MaxTrack-1); i++)
{
- Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
+ Int_t nhits=fgTrackID[etaindex][bin].fNHits[i];
if(nhits == 0) break;
if(nhits > max2)
{
- if(fTrackID[eta_index][bin].fLabel[i]!=label) {
+ if(fgTrackID[etaindex][bin].fLabel[i]!=label) {
max2 = nhits;
- label2 = fTrackID[eta_index][bin].fLabel[i];
+ label2 = fgTrackID[etaindex][bin].fLabel[i];
}
}
}
- LOG(AliL3Log::kDebug,"AliL3HoughTransformerRow::GetTrackID()","")
- <<" TrackID"<<" label "<<label<<" max "<<max<<" label2 "<<label2<<" max2 "<<max2<<" "<<(Float_t)max2/(Float_t)max<<" "<<fTrackID[eta_index][bin].fLabel[MaxTrack-1]<<" "<<(Int_t)fTrackID[eta_index][bin].fNHits[MaxTrack-1]<<ENDLOG;
+ if(max2 !=0 ) {
+ LOG(AliL3Log::kDebug,"AliL3HoughTransformerRow::GetTrackID()","")
+ <<" TrackID"<<" label "<<label<<" max "<<max<<" label2 "<<label2<<" max2 "<<max2<<" "<<(Float_t)max2/(Float_t)max<<" "<<fgTrackID[etaindex][bin].fLabel[MaxTrack-1]<<" "<<(Int_t)fgTrackID[etaindex][bin].fNHits[MaxTrack-1]<<ENDLOG;
+ }
return label;
#endif
LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
return -1;
}
-UChar_t *AliL3HoughTransformerRow::GetRowCount(Int_t eta_index)
+UChar_t *AliL3HoughTransformerRow::GetRowCount(Int_t etaindex)
{
- return fRowCount[eta_index];
+ return fgRowCount[etaindex];
}
-UChar_t *AliL3HoughTransformerRow::GetGapCount(Int_t eta_index)
+UChar_t *AliL3HoughTransformerRow::GetGapCount(Int_t etaindex)
{
- return fGapCount[eta_index];
+ return fgGapCount[etaindex];
}
-UChar_t *AliL3HoughTransformerRow::GetCurrentRowCount(Int_t eta_index)
+UChar_t *AliL3HoughTransformerRow::GetCurrentRowCount(Int_t etaindex)
{
- return fCurrentRowCount[eta_index];
+ return fgCurrentRowCount[etaindex];
}
UChar_t *AliL3HoughTransformerRow::GetTrackNRows()
{
- return fTrackNRows;
+ return fgTrackNRows;
}
UChar_t *AliL3HoughTransformerRow::GetTrackFirstRow()
{
- return fTrackFirstRow;
+ return fgTrackFirstRow;
}
UChar_t *AliL3HoughTransformerRow::GetTrackLastRow()
{
- return fTrackLastRow;
+ return fgTrackLastRow;
}
// @(#) $Id$
-#ifndef ALIL3_HOUGHTRANSFORMERROW
-#define ALIL3_HOUGHTRANSFORMERROW
+#ifndef ALIL3HOUGHTRANSFORMERROW_H
+#define ALIL3HOUGHTRANSFORMERROW_H
#include "AliL3RootTypes.h"
#include "AliL3HoughBaseTransformer.h"
class AliL3Histogram;
class AliL3HoughTransformerRow : public AliL3HoughBaseTransformer {
-
- private:
-
- AliL3Histogram **fParamSpace; //!
-
-#ifdef do_mc
- static TrackIndex **fTrackID; //!
-#endif
- Bool_t fDoMC;
-
- void DeleteHistograms();
-
- static UChar_t **fRowCount; //!
- static UChar_t **fGapCount; //!
- static UChar_t **fCurrentRowCount; //!
-
- static UChar_t *fTrackNRows; //!
- static UChar_t *fTrackFirstRow; //!
- static UChar_t *fTrackLastRow; //!
-
- Float_t *fLUTforwardZ; //!
- Float_t *fLUTforwardZ2; //!
- Float_t *fLUTbackwardZ; //!
- Float_t *fLUTbackwardZ2; //!
-
- static Float_t fBeta1,fBeta2;
public:
AliL3HoughTransformerRow();
- AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoMC=kFALSE,Float_t zvertex=0.0);
+ AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t netasegments,Bool_t DoMC=kFALSE,Float_t zvertex=0.0);
virtual ~AliL3HoughTransformerRow();
//void CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t ptres,Int_t nybin,Float_t psi);
void TransformCircle();
Int_t GetEtaIndex(Double_t eta);
- AliL3Histogram *GetHistogram(Int_t eta_index);
- Double_t GetEta(Int_t eta_index,Int_t slice);
- Int_t GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi);
- UChar_t *GetRowCount(Int_t eta_index);
- UChar_t *GetGapCount(Int_t eta_index);
- UChar_t *GetCurrentRowCount(Int_t eta_index);
+ AliL3Histogram *GetHistogram(Int_t etaindex);
+ Double_t GetEta(Int_t etaindex,Int_t slice);
+ Int_t GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi);
+ UChar_t *GetRowCount(Int_t etaindex);
+ UChar_t *GetGapCount(Int_t etaindex);
+ UChar_t *GetCurrentRowCount(Int_t etaindex);
static UChar_t *GetTrackNRows();
static UChar_t *GetTrackFirstRow();
static UChar_t *GetTrackLastRow();
- static Float_t GetBeta1() {return fBeta1;}
- static Float_t GetBeta2() {return fBeta2;}
+ static Float_t GetBeta1() {return fgBeta1;}
+ static Float_t GetBeta2() {return fgBeta2;}
+
+ private:
+
+ AliL3Histogram **fParamSpace; //!
+
+#ifdef do_mc
+ static TrackIndex **fgTrackID; //!
+#endif
+ Bool_t fDoMC; // Do MC labels or not
+
+ void DeleteHistograms(); //Method to clean up the histograms containing Hough space
+
+ static UChar_t **fgRowCount; //!
+ static UChar_t **fgGapCount; //!
+ static UChar_t **fgCurrentRowCount; //!
+
+ static UChar_t *fgTrackNRows; //!
+ static UChar_t *fgTrackFirstRow; //!
+ static UChar_t *fgTrackLastRow; //!
+
+ Float_t *fLUTforwardZ; //!
+ Float_t *fLUTforwardZ2; //!
+ Float_t *fLUTbackwardZ; //!
+ Float_t *fLUTbackwardZ2; //!
+
+ static Float_t fgBeta1,fgBeta2; // Two curves which define the Hough space
ClassDef(AliL3HoughTransformerRow,1) //TPC Rows Hough transformation class
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff