new framework
- AliHLTTPCBenchmark.h/.cxx removed 22.04.2008
old benchmark class
+- folder 'tracking' removed 18.03.2010, code has been disabled for quite a while
+++ /dev/null
-// $Id$
-// origin: hough/AliL3Histogram.cxx,v 1.31 Thu Jun 23 17:46:54 2005 UTC by hristov
-
-// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
-//*-- Copyright © ALICE HLT Group
-
-#include "AliHLTStdIncludes.h"
-
-#include "AliHLTTPCLogging.h"
-#include "AliHLTTPCHistogram.h"
-
-#if __GNUC__ >= 3
-using namespace std;
-#endif
-
-/** \class AliHLTTPCHistogram
-<pre>
-//_____________________________________________________________
-// AliHLTTPCHistogram
-//
-// 2D histogram class
-//
-</pre>
-*/
-
-//uncomment if you want overflow checks
-//#define _IFON_
-
-ClassImp(AliHLTTPCHistogram)
-
-AliHLTTPCHistogram::AliHLTTPCHistogram()
-{
- // Default constructor
- fNxbins = 0;
- fNybins = 0;
- fNcells = 0;
- fXmin = 0;
- fYmin = 0;
- fXmax = 0;
- fYmax = 0;
- fBinwidthX = 0;
- fBinwidthY = 0;
- fFirstXbin = 0;
- fLastXbin = 0;
- fFirstYbin = 0;
- fLastYbin = 0;
- fEntries = 0;
- fContent = 0;
- fThreshold = 0;
- fRootHisto = 0;
-}
-
-AliHLTTPCHistogram::AliHLTTPCHistogram(Char_t *name,Char_t */*id*/,
- 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;
- fNybins = nybin;
- fNcells = (nxbin+2)*(nybin+2);
- fXmin = xmin;
- fYmin = ymin;
- fXmax = xmax;
- fYmax = ymax;
- fBinwidthX = (fXmax - fXmin) / fNxbins;
- fBinwidthY = (fYmax - fYmin) / fNybins;
-
- fEntries = 0;
- fFirstXbin = 1;
- fFirstYbin = 1;
- fLastXbin = nxbin;
- fLastYbin = nybin;
- fRootHisto = 0;
- fThreshold = 0;
-
- fContent = new Int_t[fNcells];
- Reset();
-}
-
-AliHLTTPCHistogram::~AliHLTTPCHistogram()
-{
- //Destructor
- if(fContent)
- delete [] fContent;
- if(fRootHisto)
- delete fRootHisto;
-}
-
-void AliHLTTPCHistogram::Reset()
-{
- // Reset histogram contents
- if(fContent)
- for(Int_t i=0; i<fNcells; i++) fContent[i] = 0;
-
- fEntries=0;
-}
-
-void AliHLTTPCHistogram::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)
- return;
-#endif
-
- AddBinContent(bin,weight);
-}
-
-void AliHLTTPCHistogram::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_
- if(bin < 0)
- return;
-#endif
-
- AddBinContent(bin,weight);
-}
-
-void AliHLTTPCHistogram::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_
- if(bin < 0)
- return;
-#endif
-
- AddBinContent(bin,weight);
-}
-
-void AliHLTTPCHistogram::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)
- return;
-#endif
-
- AddBinContent(bin,weight);
-}
-
-Int_t AliHLTTPCHistogram::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_
- if(!xbin || !ybin)
- return -1;
-#endif
-
- return GetBin(xbin,ybin);
-}
-
-Int_t AliHLTTPCHistogram::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_
- if(!xbin || !ybin)
- return -1;
-#endif
-
- return GetLabelBin(xbin,ybin);
-}
-
-Int_t AliHLTTPCHistogram::FindXbin(Double_t x) const
-{
- // Finds the bin which correspond to x
- if(x < fXmin || x > fXmax)
- return 0;
-
- return 1 + (Int_t)(fNxbins*(x-fXmin)/(fXmax-fXmin));
-}
-
-Int_t AliHLTTPCHistogram::FindYbin(Double_t y) const
-{
- // Finds the bin which correspond to y
- if(y < fYmin || y > fYmax)
- return 0;
-
- return 1 + (Int_t)(fNybins*(y-fYmin)/(fYmax-fYmin));
-}
-
-Int_t AliHLTTPCHistogram::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)
- return 0;
-
- return xbin + ybin*(fNxbins+2);
-}
-
-Int_t AliHLTTPCHistogram::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)
- return -1;
-
- return (Int_t)(xbin/2) + ((Int_t)(ybin/2))*((Int_t)((fNxbins+3)/2));
-}
-
-Int_t AliHLTTPCHistogram::GetBinContent(Int_t bin) const
-{
- // Return the bin content
- if(bin >= fNcells)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::GetBinContent","array")<<AliHLTTPCLog::kDec<<
- "bin out of range "<<bin<<ENDLOG;
- return 0;
- }
-
- if(fContent[bin] < fThreshold)
- return 0;
- return fContent[bin];
-}
-
-void AliHLTTPCHistogram::SetBinContent(Int_t xbin,Int_t ybin,Int_t value)
-{
- // Set bin content
- Int_t bin = GetBin(xbin,ybin);
-#ifdef _IFON_
- if(bin == 0)
- return;
-#endif
-
- SetBinContent(bin,value);
-}
-
-void AliHLTTPCHistogram::SetBinContent(Int_t bin,Int_t value)
-{
- // Set bin content
-
- if(bin >= fNcells)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::SetBinContent","array")<<AliHLTTPCLog::kDec<<
- "bin out of range "<<bin<<ENDLOG;
- return;
- }
-
- if(bin == 0)
- return;
- fContent[bin]=value;
-}
-
-void AliHLTTPCHistogram::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)
- return;
-#endif
-
- AddBinContent(bin,weight);
-}
-
-void AliHLTTPCHistogram::AddBinContent(Int_t bin,Int_t weight)
-{
- // Adds weight to bin content
- if(bin < 0 || bin > fNcells)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::AddBinContent","array")<<AliHLTTPCLog::kDec<<
- "bin-value out of range "<<bin<<ENDLOG;
- return;
- }
- if(bin == 0)
- return;
- fEntries++;
- fContent[bin] += weight;
-}
-
-void AliHLTTPCHistogram::Add(AliHLTTPCHistogram *h1,Double_t /*weight*/)
-{
- //Adding two histograms. Should be identical.
-
- if(!h1)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::Add","Pointer")<<
- "Attempting to add a non-existing histogram"<<ENDLOG;
- return;
- }
-
- if(h1->GetNbinsX()!=fNxbins || h1->GetNbinsY()!=fNybins)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::Add","array")<<
- "Mismatch in the number of bins "<<ENDLOG;
- return;
- }
-
- if(h1->GetFirstXbin()!=fFirstXbin || h1->GetLastXbin()!=fLastXbin ||
- h1->GetFirstYbin()!=fFirstYbin || h1->GetLastYbin()!=fLastYbin)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::Add","array")<<
- "Mismatch in the bin numbering "<<ENDLOG;
- return;
- }
-
- for(Int_t bin=0; bin<fNcells; bin++)
- fContent[bin] += h1->GetBinContent(bin);
-
- fEntries += h1->GetNEntries();
-}
-
-Double_t AliHLTTPCHistogram::GetBinCenterX(Int_t xbin) const
-{
- // Returns the position of the center of a bin
- if(xbin < fFirstXbin || xbin > fLastXbin)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::GetBinCenterX","xbin")
- <<"Bin-value out of range "<<xbin<<ENDLOG;
- return -1;
- }
-
- return fXmin + (xbin-0.5) * fBinwidthX;
-}
-
-Double_t AliHLTTPCHistogram::GetBinCenterY(Int_t ybin) const
-{
- // Returns the position of the center of a bin
- if(ybin < fFirstYbin || ybin > fLastYbin)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::GetBinCenterY","ybin")
- <<"Bin-value out of range "<<ybin<<ENDLOG;
- return -1;
- }
-
- return fYmin + (ybin-0.5) * fBinwidthY;
-}
-
-Double_t AliHLTTPCHistogram::GetPreciseBinCenterX(Float_t xbin) const
-{
- // Returns the position of the center of a bin using precise values inside the bin
- if(xbin < (fFirstXbin-1.5) || xbin > (fLastXbin+1.5))
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::GetBinCenterX","xbin")
- <<"Bin-value out of range "<<xbin<<ENDLOG;
- return -1;
- }
- // return fXmin + (xbin-1) * fBinwidthX + 0.5*fBinwidthX;
- return fXmin + (xbin-0.5) * fBinwidthX;
-}
-
-Double_t AliHLTTPCHistogram::GetPreciseBinCenterY(Float_t ybin) const
-{
- // Returns the position of the center of a bin using precise values inside the bin
- if(ybin < (fFirstYbin-1.5) || ybin > (fLastYbin+1.5))
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::GetBinCenterY","ybin")
- <<"Bin-value out of range "<<ybin<<ENDLOG;
- return -1;
- }
- // return fYmin + (ybin-1) * fBinwidthY + 0.5*fBinwidthY;
- return fYmin + (ybin-0.5) * fBinwidthY;
-}
-
-void AliHLTTPCHistogram::Draw(Char_t *option)
-{
- // Fill the contents of the corresponding ROOT histogram and draws it
- if(!fRootHisto)
- CreateRootHisto();
-
- for(Int_t xbin=GetFirstXbin(); xbin<=GetLastXbin(); xbin++)
- {
- for(Int_t ybin=GetFirstYbin(); ybin<=GetLastYbin(); ybin++)
- {
- Int_t bin = GetBin(xbin,ybin);
- fRootHisto->Fill(GetBinCenterX(xbin),GetBinCenterY(ybin),GetBinContent(bin));
- }
- }
-
- //fRootHisto->SetStats(kFALSE);
- fRootHisto->Draw(option);
- return;
-}
-
-void AliHLTTPCHistogram::CreateRootHisto()
-{
- // Create ROOT histogram out of AliHLTTPCHistogram
- fRootHisto = new TH2F(fName,"",fNxbins,fXmin,fXmax,fNybins,fYmin,fYmax);
- return;
-}
-
-ofstream& operator<<(ofstream &o, const AliHLTTPCHistogram &h)
-{
- for(Int_t xbin=h.GetFirstXbin(); xbin<=h.GetLastXbin(); xbin++)
- {
- for(Int_t ybin=h.GetFirstYbin(); ybin<=h.GetLastYbin(); ybin++)
- {
- Int_t bin = h.GetBin(xbin,ybin);
- o << h.GetBinContent(bin) << " ";
- }
- o << endl;
- }
- return o;
-}
+++ /dev/null
-// $Id$
-// origin: hough/AliL3Histogram.h,v 1.28 Wed Sep 29 12:01:17 2004 UTC by cvetan
-
-#ifndef ALIHLTTPCHISTOGRAM_H
-#define ALIHLTTPCHISTOGRAM_H
-
-#include "AliHLTStdIncludes.h"
-#include "AliHLTTPCRootTypes.h"
-
-#include <TH2.h>
-
-class AliHLTTPCHistogram {
- public:
-
- AliHLTTPCHistogram();
- AliHLTTPCHistogram(Char_t *name,Char_t *id,Int_t nxbin,Double_t xmin,Double_t xmax,Int_t nybin,Double_t ymin,Double_t ymax);
- virtual ~AliHLTTPCHistogram();
-
- void Reset();
- virtual void Fill(Double_t x,Double_t y,Int_t weight=1);
- virtual void Fill(Double_t x,Int_t ybin,Int_t weight=1);
- virtual void Fill(Int_t xbin,Double_t y,Int_t weight=1);
- virtual void Fill(Int_t xbin,Int_t ybin,Int_t weight=1);
- virtual Int_t FindBin(Double_t x,Double_t y) const;
- virtual Int_t FindLabelBin(Double_t x,Double_t y) const;
- virtual Int_t FindXbin(Double_t x) const;
- virtual Int_t FindYbin(Double_t y) const;
- Int_t GetBin(Int_t xbin,Int_t ybin) const;
- Int_t GetLabelBin(Int_t xbin,Int_t ybin) const;
- Int_t GetBinContent(Int_t bin) const;
- void SetBinContent(Int_t xbin,Int_t ybin,Int_t value);
- void SetBinContent(Int_t bin,Int_t value);
- void AddBinContent(Int_t xbin,Int_t ybin,Int_t weight);
- void AddBinContent(Int_t bin,Int_t weight);
- void Add(AliHLTTPCHistogram *h1,Double_t weight=1);
- void SetThreshold(Int_t i) {fThreshold = i;}
- void CreateRootHisto();
- virtual void Draw(Char_t *option="hist");
- virtual void Print() const {};
-
- friend ofstream& operator<< (ofstream &o, const AliHLTTPCHistogram &h);
-
- TH2F *GetRootHisto();
-
- Double_t GetXmin() const {return fXmin;}
- Double_t GetXmax() const {return fXmax;}
- Double_t GetYmin() const {return fYmin;}
- Double_t GetYmax() const {return fYmax;}
- virtual Double_t GetBinCenterX(Int_t xbin) const;
- virtual Double_t GetBinCenterY(Int_t ybin) const;
- Double_t GetPreciseBinCenterX(Float_t xbin) const;
- Double_t GetPreciseBinCenterY(Float_t ybin) const;
- Double_t GetBinWidthX() const {return fBinwidthX;}
- Double_t GetBinWidthY() const {return fBinwidthY;}
- Int_t GetFirstXbin() const {return fFirstXbin;}
- Int_t GetLastXbin() const {return fLastXbin;}
- Int_t GetFirstYbin() const {return fFirstYbin;}
- Int_t GetLastYbin() const {return fLastYbin;}
- Int_t GetNbinsX() const {return fNxbins;}
- Int_t GetNbinsY() const {return fNybins;}
- Int_t GetNEntries() const {return fEntries;}
-
- Int_t *fContent; //!
- Int_t *GetContentArray() const {return fContent;}
-
- protected:
- 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
-
- TH2F *fRootHisto; // Corresponding ROOT histogram
-
- private:
- Double_t fBinwidthX; // Bin width of the Hough space
- Double_t fBinwidthY; // Bin width of the Hough space
-
- ClassDef(AliHLTTPCHistogram,1) //2D histogram class
-
-};
-
-inline TH2F *AliHLTTPCHistogram::GetRootHisto()
-{
- if(!fRootHisto)
- {
- STDCERR<<"AliHLTTPCHistogram::GetRootHisto() : You must first Draw histogram before accessing it"<<STDENDL;
- return 0;
- }
- else
- return fRootHisto;
-}
-
-#endif
+++ /dev/null
-// @(#) $Id$
-// origin: hough/AliL3Histogram1D.cxx,v 1.11 Tue Jun 14 10:55:20 2005 UTC by cvetan
-
-// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
-//*-- Copyright © ALICE HLT Group
-
-#include <strings.h>
-#include "AliHLTStdIncludes.h"
-
-#include "AliHLTTPCLogging.h"
-#include "AliHLTTPCHistogram1D.h"
-
-#include <TH1.h>
-
-#if __GNUC__ >= 3
-using namespace std;
-#endif
-
-//_____________________________________________________________
-// AliHLTTPCHistogram1D
-//
-// 1D histogram class.
-
-ClassImp(AliHLTTPCHistogram1D)
-
-AliHLTTPCHistogram1D::AliHLTTPCHistogram1D()
-{
- //default ctor
- fNbins = 0;
- fNcells = 0;
- fEntries = 0;
- fXmin = 0;
- fXmax = 0;
- fRootHisto = 0;
- fThreshold = 0;
- fContent = 0;
-
-}
-
-AliHLTTPCHistogram1D::AliHLTTPCHistogram1D(Char_t *name,Char_t */*id*/,Int_t nxbin,Double_t xmin,Double_t xmax)
-
-{
- //normal ctor
- strcpy(fName,name);
- fNbins = nxbin;
- fNcells = fNbins + 2;
- fEntries = 0;
- fXmin = xmin;
- fXmax = xmax;
- fRootHisto = 0;
- fThreshold = 0;
-
- fContent = new Double_t[fNcells];
- Reset();
-}
-
-AliHLTTPCHistogram1D::~AliHLTTPCHistogram1D()
-{
- //Destructor
- if(fContent)
- delete [] fContent;
- if(fRootHisto)
- delete fRootHisto;
-}
-
-
-void AliHLTTPCHistogram1D::Reset()
-{
- //Reset histogram contents
-#if defined(__DECCXX)
- bzero((char *)fContent,fNcells*sizeof(Double_t));
-#else
- bzero(fContent,fNcells*sizeof(Double_t));
-#endif
- fEntries=0;
-}
-
-void AliHLTTPCHistogram1D::Fill(Double_t x,Int_t weight)
-{
- //Fill a given bin with weight
- Int_t bin = FindBin(x);
- AddBinContent(bin,weight);
-}
-
-
-Int_t AliHLTTPCHistogram1D::FindBin(Double_t x) const
-{
- //Find a given bin
- if(x < fXmin || x > fXmax)
- return 0;
-
- return 1 + (Int_t)(fNbins*(x-fXmin)/(fXmax-fXmin));
-
-}
-
-Int_t AliHLTTPCHistogram1D::GetMaximumBin() const
-{
- //Find the bin with the largest content
- Double_t maxvalue=0;
- Int_t maxbin=0;
- for(Int_t i=0; i<fNcells; i++)
- {
- if(fContent[i] > maxvalue)
- {
- maxvalue=fContent[i];
- maxbin = i;
- }
- }
- return maxbin;
-}
-
-Double_t AliHLTTPCHistogram1D::GetBinContent(Int_t bin) const
-{
- //Get bin content
- if(bin >= fNcells)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::GetBinContent","array")<<AliHLTTPCLog::kDec<<
- "bin out of range "<<bin<<ENDLOG;
- return 0;
- }
-
- if(fContent[bin] < fThreshold)
- return 0;
- return fContent[bin];
-}
-
-
-void AliHLTTPCHistogram1D::SetBinContent(Int_t bin,Int_t value)
-{
- //Set bin content
- if(bin >= fNcells)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::SetBinContent","array")<<AliHLTTPCLog::kDec<<
- "bin out of range "<<bin<<ENDLOG;
- return;
- }
- if(bin == 0)
- return;
- fContent[bin]=value;
-
-}
-
-void AliHLTTPCHistogram1D::AddBinContent(Int_t bin,Int_t weight)
-{
- //Add weight to bin content
- if(bin < 0 || bin > fNcells)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogram::AddBinContent","array")<<AliHLTTPCLog::kDec<<
- "bin-value out of range "<<bin<<ENDLOG;
- return;
- }
- if(bin == 0)
- return;
- fEntries++;
- fContent[bin] += weight;
-}
-
-Double_t AliHLTTPCHistogram1D::GetBinCenter(Int_t bin) const
-{
- //Get bin center
- Double_t binwidth = (fXmax - fXmin) / fNbins;
- return fXmin + (bin-1) * binwidth + 0.5*binwidth;
-
-}
-
-void AliHLTTPCHistogram1D::Draw(Char_t *option)
-{
- //Draw the histogram
- fRootHisto = new TH1F(fName,"",fNbins,fXmin,fXmax);
- for(Int_t bin=0; bin<fNcells; bin++)
- fRootHisto->AddBinContent(bin,GetBinContent(bin));
-
- fRootHisto->Draw(option);
-
-}
+++ /dev/null
-// @(#) $Id$
-// origin: hough/AliL3Histogram1D.h,v 1.5 Thu Jun 17 13:18:42 2004 UTC by cvetan
-
-#ifndef ALIHLTTPCHISTOGRAM1D_H
-#define ALIHLTTPCHISTOGRAM1D_H
-
-#include "AliHLTTPCRootTypes.h"
-
-class TH1F;
-
-class AliHLTTPCHistogram1D {
-
- public:
- AliHLTTPCHistogram1D();
- AliHLTTPCHistogram1D(Char_t *name,Char_t *id,Int_t nxbin,Double_t xmin,Double_t xmax);
- virtual ~AliHLTTPCHistogram1D();
-
- void Reset();
- void Fill(Double_t x,Int_t weight=1);
- void AddBinContent(Int_t bin,Int_t weight);
- Int_t GetMaximumBin() const;
- Int_t FindBin(Double_t x) const;
- Double_t GetBinContent(Int_t bin) const;
- Double_t GetBinCenter(Int_t bin) const;
- Int_t GetNEntries() const {return fEntries;}
-
- void SetBinContent(Int_t bin,Int_t value);
- void SetThreshold(Int_t i) {fThreshold = i;}
-
- void Draw(Char_t *option="hist");
- TH1F *GetRootHisto() {return fRootHisto;}
-
- private:
-
- Double_t *fContent; //!
- Char_t fName[100];//Histogram title
- Int_t fNbins;//Number of bins
- Int_t fNcells;//Number of cells
- Int_t fEntries;//Number of entries
-
- Int_t fThreshold;//Bin content threshold
- Double_t fXmin;//Lower limit in X
- Double_t fXmax;//Upper limit in X
-
-
- TH1F *fRootHisto;//The corresponding ROOT histogram
-
- ClassDef(AliHLTTPCHistogram1D,1) //1D histogram class
-
-};
-
-#endif
+++ /dev/null
-// $Id$
-// origin: hough/AliL3HistogramAdaptive.cxx,v 1.13 Thu Jun 23 17:46:54 2005 UTC by hristov
-
-// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
-//*-- Copyright © ALICE HLT Group
-
-#include "AliHLTStdIncludes.h"
-#include "AliHLTTPCLogging.h"
-#include "AliHLTTPCHistogramAdaptive.h"
-#include "AliHLTTPCTransform.h"
-#include "AliHLTTPCTrack.h"
-
-#if __GNUC__ >= 3
-using namespace std;
-#endif
-
-//_____________________________________________________________
-// AliHLTTPCHistogramAdaptive
-//
-// 2D histogram class adapted for kappa and psi as used in the Circle Hough Transform.
-// The bins in kappa is not linear, but has a width which is specified by argument
-// ptres in the constructor. This gives the relative pt resolution which should
-// be kept throughout the kappa range.
-
-ClassImp(AliHLTTPCHistogramAdaptive)
-
-AliHLTTPCHistogramAdaptive::AliHLTTPCHistogramAdaptive() : AliHLTTPCHistogram()
-{
- //default ctor
- fKappaBins=0;
-}
-
-
-AliHLTTPCHistogramAdaptive::AliHLTTPCHistogramAdaptive(Char_t *name,Double_t minpt,Double_t maxpt,Double_t ptres,
- Int_t nybins,Double_t ymin,Double_t ymax)
-{
- //normal ctor
- strcpy(fName,name);
-
- fPtres = ptres;
- fXmin = -1*AliHLTTPCTransform::GetBFact()*AliHLTTPCTransform::GetBField()/minpt;
- fXmax = AliHLTTPCTransform::GetBFact()*AliHLTTPCTransform::GetBField()/minpt;
-
- fMinPt = minpt;
- fMaxPt = maxpt;
- fNxbins = InitKappaBins();
- fNybins = nybins;
-
- fYmin = ymin;
- fYmax = ymax;
- fFirstXbin=1;
- fFirstYbin=1;
- fLastXbin = fNxbins;
- fLastYbin = fNybins;
- fNcells = (fNxbins+2)*(fNybins+2);
-
- fThreshold=0;
- fContent = new Int_t[fNcells];
- Reset();
-}
-
-AliHLTTPCHistogramAdaptive::~AliHLTTPCHistogramAdaptive()
-{
- //dtor
- if(fKappaBins)
- delete [] fKappaBins;
-}
-
-Int_t AliHLTTPCHistogramAdaptive::InitKappaBins()
-{
- //Here a LUT for the kappa values created. This has to be done since
- //the binwidth in kappa is not constant, but change according to the
- //set relative resolution in pt.
- //Since the kappa values are symmetric about origo, the size of the
- //LUT is half of the total number of bins in kappa direction.
-
- Double_t pt = fMinPt,deltapt,localpt;
- Int_t bin=0;
-
- while(pt < fMaxPt)
- {
- localpt = pt;
- deltapt = fPtres*localpt*localpt;
- pt += 2*deltapt;
- bin++;
- }
- fKappaBins = new Double_t[bin+1];
- pt=fMinPt;
- bin=0;
- fKappaBins[bin] = AliHLTTPCTransform::GetBFact()*AliHLTTPCTransform::GetBField()/fMinPt;
- while(pt < fMaxPt)
- {
- localpt = pt;
- deltapt = fPtres*localpt*localpt;
- pt += 2*deltapt; //*2 because pt +- 1/2*deltapt is one bin
- bin++;
- fKappaBins[bin] = AliHLTTPCTransform::GetBFact()*AliHLTTPCTransform::GetBField()/pt;
- }
- return (bin+1)*2; //Both negative and positive kappa.
-}
-
-
-void AliHLTTPCHistogramAdaptive::Fill(Double_t x,Double_t y,Int_t weight)
-{
- //Fill a given bin in the histogram
- Int_t bin = FindBin(x,y);
- if(bin < 0)
- return;
- AddBinContent(bin,weight);
-
-}
-
-Int_t AliHLTTPCHistogramAdaptive::FindBin(Double_t x,Double_t y) const
-{
- //Find a bin in the histogram
- Int_t xbin = FindXbin(x);
- Int_t ybin = FindYbin(y);
-
- if(!xbin || !ybin)
- return -1;
- return GetBin(xbin,ybin);
-}
-
-Int_t AliHLTTPCHistogramAdaptive::FindXbin(Double_t x) const
-{
- //Find X bin in the histogram
- if(x < fXmin || x > fXmax || fabs(x) < fKappaBins[(fNxbins/2-1)])
- return 0;
-
- //Remember that kappa value is decreasing with bin number!
- //Also, the bin numbering starts at 1 and ends at fNxbins,
- //so the corresponding elements in the LUT is bin - 1.
-
- Int_t bin=0;
- while(bin < fNxbins/2)
- {
- if(fabs(x) <= fKappaBins[bin] && fabs(x) > fKappaBins[bin+1])
- break;
- bin++;
- }
- if(x < 0)
- return bin + 1;
- else
- return fNxbins - bin;
-
-}
-
-Int_t AliHLTTPCHistogramAdaptive::FindYbin(Double_t y) const
-{
- //Find Y bin in the histogram
- if(y < fYmin || y > fYmax)
- return 0;
-
- return 1 + (Int_t)(fNybins*(y-fYmin)/(fYmax-fYmin));
-}
-
-Double_t AliHLTTPCHistogramAdaptive::GetBinCenterX(Int_t xbin) const
-{
- //Returns bin center in X
- if(xbin < fFirstXbin || xbin > fLastXbin)
- {
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHistogramAdaptive::GetBinCenterX","Bin-value")
- <<"XBinvalue out of range "<<xbin<<ENDLOG;
- return 0;
- }
-
- //The bin numbers go from 1 to fNxbins, so the corresponding
- //element in the LUT is xbin - 1. This is the reason why we
- //substract a 1 here:
-
- Int_t bin = xbin;
- bin -= 1;
- if(bin >= fNxbins/2)
- bin = fNxbins - 1 - bin;
-
- //Remember again that the kappa-values are _decreasing_ with bin number.
-
- Double_t binwidth = fKappaBins[bin] - fKappaBins[bin+1];
- Double_t kappa = fKappaBins[bin] - 0.5*binwidth;
- if(xbin < fNxbins/2)
- return -1.*kappa;
- else
- return kappa;
-
-}
-
-Double_t AliHLTTPCHistogramAdaptive::GetBinCenterY(Int_t ybin) const
-{
- //Returns bin center in Y
- if(ybin < fFirstYbin || ybin > fLastYbin)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogramAdaptive::GetBinCenterY","ybin")
- <<"Bin-value out of range "<<ybin<<ENDLOG;
- return -1;
- }
- Double_t binwidth = (fYmax - fYmin) / fNybins;
- return fYmin + (ybin-0.5) * binwidth;
-
-}
-
-
-void AliHLTTPCHistogramAdaptive::Draw(Char_t *option)
-{
- //Draw the histogram
- if(!fRootHisto)
- CreateRootHisto();
-
- Double_t kappa,psi;
- Int_t content,bin;
- for(Int_t i=fFirstXbin; i<=fLastXbin; i++)
- {
- kappa = GetBinCenterX(i);
- for(Int_t j=fFirstYbin; j<=fLastYbin; j++)
- {
- psi = GetBinCenterY(j);
- bin = GetBin(i,j);
- content = GetBinContent(bin);
- fRootHisto->Fill(kappa,psi,content);
- }
- }
- fRootHisto->Draw(option);
- return;
-}
-
-void AliHLTTPCHistogramAdaptive::Print() const
-{
- //Print the contents of the histogram
- cout<<"Printing content of histogram "<<fName<<endl;
- for(Int_t i=0; i<fNcells; i++)
- {
- if(GetBinContent(i)==0) continue;
- cout<<"Bin "<<i<<": "<<GetBinContent(i)<<endl;
- }
-
-}
+++ /dev/null
-// @(#) $Id$
-// origin hough/AliL3HistogramAdaptive.h,v 1.9 Sun Apr 30 16:37:32 2006 UTC by hristov
-
-#ifndef ALIHLTTPCHISTOGRAMADAPTIVE_H
-#define ALIHLTTPCHISTOGRAMADAPTIVE_H
-
-#include "AliHLTTPCRootTypes.h"
-#include "AliHLTTPCHistogram.h"
-
-class AliHLTTPCHistogramAdaptive : public AliHLTTPCHistogram {
-
- public:
- AliHLTTPCHistogramAdaptive();
- AliHLTTPCHistogramAdaptive(Char_t *name,Double_t minpt,Double_t maxpt,Double_t ptres,
- Int_t nybins,Double_t ymin,Double_t ymax);
- ~AliHLTTPCHistogramAdaptive();
-
- void Fill(Double_t x,Double_t y,Int_t weight=1);
- void Fill(Double_t x,Int_t ybin,Int_t weight=1) {
- AliHLTTPCHistogram::Fill(x,ybin,weight);
- }
- void Fill(Int_t xbin,Double_t y,Int_t weight=1) {
- AliHLTTPCHistogram::Fill(xbin,y,weight);
- }
- void Fill(Int_t xbin,Int_t ybin,Int_t weight=1) {
- AliHLTTPCHistogram::Fill(xbin,ybin,weight);
- }
- Int_t FindBin(Double_t x,Double_t y) const;
- Int_t FindXbin(Double_t x) const;
- Int_t FindYbin(Double_t x) const;
- void Draw(Char_t *option = "hist");
- void Print() const;
-
- Double_t GetBinCenterX(Int_t xbin) const;
- Double_t GetBinCenterY(Int_t ybin) const;
-
- private:
- Double_t fPtres;//The desired Pt resolution
- Double_t fMinPt;//Minimum Pt
- Double_t fMaxPt;//Maximum Pt
- Double_t *fKappaBins; //!
-
- Int_t InitKappaBins();
-
- ClassDef(AliHLTTPCHistogramAdaptive,1) //2D histogram class
-
-};
-
-#endif
+++ /dev/null
-// $Id$
-// origin: hough/AliL3Hough.cxx,v 1.50 Tue Mar 28 18:05:12 2006 UTC by alibrary
-
-//**************************************************************************
-//* This file is property of and copyright by the ALICE HLT Project *
-//* ALICE Experiment at CERN, All rights reserved. *
-//* *
-//* Primary Authors: Anders Vestbo, Cvetan Cheshkov *
-//* for The ALICE HLT Project. *
-//* *
-//* Permission to use, copy, modify and distribute this software and its *
-//* documentation strictly for non-commercial purposes is hereby granted *
-//* without fee, provided that the above copyright notice appears in all *
-//* copies and that both the copyright notice and this permission notice *
-//* appear in the supporting documentation. The authors make no claims *
-//* about the suitability of this software for any purpose. It is *
-//* provided "as is" without express or implied warranty. *
-//**************************************************************************
-
-/** @file AliHLTTPCHough.cxx
- @author Anders Vestbo, Cvetan Cheshkov
- @date
- @brief Steering for HLT TPC hough transform tracking algorithms. */
-
-
-#include "AliHLTStdIncludes.h"
-#include "AliLog.h"
-#include <sys/time.h>
-
-#include "AliHLTTPCLogging.h"
-
-#ifdef HAVE_ALIHLTHOUGHMERGER
-#include "AliHLTHoughMerger.h"
-#endif //HAVE_ALIHLTHOUGHMERGER
-
-#ifdef HAVE_ALIHLTHOUGHINTMERGER
-#include "AliHLTHoughIntMerger.h"
-#endif //HAVE_ALIHLTHOUGHINTMERGER
-
-#ifdef HAVE_ALIHLTHOUGHGLOBALMERGER
-#include "AliHLTHoughGlobalMerger.h"
-#endif //HAVE_ALIHLTHOUGHGLOBALMERGER
-
-#include "AliHLTTPCHistogram.h"
-#include "AliHLTTPCHough.h"
-
-#ifdef HAVE_ALIHLTHOUGHTRANSFORMERDEFAULT
-// the original AliHLTHoughBaseTransformer has been renamed to
-// AliHLTTPCHoughTransformer and AliHLTHoughTransformer to
-// AliHLTHoughTransformerDefault, but the latter is not yet
-// migrated
-#include "AliHLTHoughTransformer.h"
-#endif // HAVE_ALIHLTHOUGHTRANSFORMERDEFAULT
-
-#ifdef HAVE_ALIHLTHOUGHCLUSTERTRANSFORMER
-#include "AliHLTHoughClusterTransformer.h"
-#endif // HAVE_ALIHLTHOUGHCLUSTERTRANSFORMER
-
-#ifdef HAVE_ALIHLTHOUGHTRANSFORMERLUT
-#include "AliHLTHoughTransformerLUT.h"
-#endif // HAVE_ALIHLTHOUGHTRANSFORMERLUT
-
-#ifdef HAVE_ALIHLTHOUGHTRANSFORMERVHDL
-#include "AliHLTHoughTransformerVhdl.h"
-#endif // HAVE_ALIHLTHOUGHTRANSFORMERVHDL
-
-#include "AliHLTTPCHoughTransformerRow.h"
-
-#include "AliHLTTPCHoughMaxFinder.h"
-#include "AliHLTTPCBenchmark.h"
-#include "AliHLTTPCFileHandler.h"
-
-#ifdef HAVE_ALIHLTDATAHANDLER
-#include "AliHLTDataHandler.h"
-#endif // HAVE_ALIHLTDATAHANDLER
-
-//#include "AliHLTDigitData.h"
-#include "AliHLTTPCHoughEval.h"
-#include "AliHLTTPCTransform.h"
-#include "AliHLTTPCTrackArray.h"
-#include "AliHLTTPCHoughTrack.h"
-
-#ifdef HAVE_ALIHLTDDLDATAFILEHANDLER
-#include "AliHLTDDLDataFileHandler.h"
-#endif // HAVE_ALIHLTDDLDATAFILEHANDLER
-
-#include "AliHLTTPCHoughKalmanTrack.h"
-
-#ifdef HAVE_THREAD
-#include "TThread.h"
-#endif // HAVE_THREAD
-#include <AliRunLoader.h>
-#include <AliRawEvent.h>
-#include <AliESDEvent.h>
-#include <AliESDtrack.h>
-#include <AliESDHLTtrack.h>
-
-#if __GNUC__ >= 3
-using namespace std;
-#endif
-
-/** ROOT macro for the implementation of ROOT specific class methods */
-ClassImp(AliHLTTPCHough);
-
-AliHLTTPCHough::AliHLTTPCHough()
-{
- //Constructor
-
- fBinary = kFALSE;
- fAddHistograms = kFALSE;
- fDoIterative = kFALSE;
- fWriteDigits = kFALSE;
- fUse8bits = kFALSE;
-
- fMemHandler = 0;
- fHoughTransformer = 0;
- fEval = 0;
- fPeakFinder = 0;
- fTracks = 0;
- fGlobalTracks = 0;
- fMerger = 0;
- fInterMerger = 0;
- fGlobalMerger = 0;
- fBenchmark = 0;
-
- fNEtaSegments = 0;
- fNPatches = 0;
- fLastPatch =-1;
- fVersion = 0;
- fCurrentSlice = 0;
- fEvent = 0;
-
- fKappaSpread = 6;
- fPeakRatio = 0.5;
- fInputFile = 0;
- fInputPtr = 0;
- fRawEvent = 0;
-
- SetTransformerParams();
- SetThreshold();
- SetNSaveIterations();
- SetPeakThreshold();
- //just be sure that index is empty for new event
- AliHLTTPCFileHandler::CleanStaticIndex();
- fRunLoader = 0;
-#ifdef HAVE_THREAD
- fThread = 0;
-#endif // HAVE_THREAD
-}
-
-AliHLTTPCHough::AliHLTTPCHough(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 = netasegments;
- fAddHistograms = kFALSE;
- fDoIterative = kFALSE;
- fWriteDigits = kFALSE;
- fUse8bits = bit8;
- fVersion = tv;
- fKappaSpread=6;
- fPeakRatio=0.5;
- if(!fBinary) {
- if(infile) {
- fInputFile = infile;
- fInputPtr = 0;
- }
- else {
- fInputFile = 0;
- fInputPtr = ptr;
- }
- }
- else {
- fInputFile = 0;
- fInputPtr = 0;
- }
- //just be sure that index is empty for new event
- AliHLTTPCFileHandler::CleanStaticIndex();
- fRunLoader = 0;
-#ifdef HAVE_THREAD
- fThread = 0;
-#endif // HAVE_THREAD
-}
-
-AliHLTTPCHough::~AliHLTTPCHough()
-{
- //dtor
-
- CleanUp();
-
-#ifdef HAVE_ALIHLTHOUGHMERGER
- if(fMerger)
- delete fMerger;
-#endif //HAVE_ALIHLTHOUGHMERGER
- //cout << "Cleaned class merger " << endl;
-#ifdef HAVE_ALIHLTHOUGHINTMERGER
- if(fInterMerger)
- delete fInterMerger;
-#endif //HAVE_ALIHLTHOUGHINTMERGER
- //cout << "Cleaned class inter " << endl;
- if(fPeakFinder)
- delete fPeakFinder;
- //cout << "Cleaned class peak " << endl;
-#ifdef HAVE_ALIHLTHOUGHGLOBALMERGER
- if(fGlobalMerger)
- delete fGlobalMerger;
-#endif //HAVE_ALIHLTHOUGHGLOBALMERGER
- //cout << "Cleaned class global " << endl;
- if(fBenchmark)
- delete fBenchmark;
- //cout << "Cleaned class bench " << endl;
- if(fGlobalTracks)
- delete fGlobalTracks;
- //cout << "Cleaned class globaltracks " << endl;
-#ifdef HAVE_THREAD
- if(fThread) {
- // fThread->Delete();
- delete fThread;
- fThread = 0;
- }
-#endif // HAVE_THREAD
-}
-
-void AliHLTTPCHough::CleanUp()
-{
- //Cleanup memory
-
- for(Int_t i=0; i<fNPatches; i++)
- {
- if(fTracks[i]) delete fTracks[i];
- //cout << "Cleaned tracks " << i << endl;
- if(fEval[i]) delete fEval[i];
- //cout << "Cleaned eval " << i << endl;
- if(fHoughTransformer[i]) delete fHoughTransformer[i];
- //cout << "Cleaned traf " << i << endl;
- if(fMemHandler[i]) delete fMemHandler[i];
- //cout << "Cleaned mem " << i << endl;
- }
-
- if(fTracks) delete [] fTracks;
- //cout << "Cleaned class tracks " << endl;
- if(fEval) delete [] fEval;
- //cout << "Cleaned class eval " << endl;
- if(fHoughTransformer) delete [] fHoughTransformer;
- //cout << "Cleaned cleass trafo " << endl;
- if(fMemHandler) delete [] fMemHandler;
- //cout << "Cleaned class mem " << endl;
-}
-
-void AliHLTTPCHough::Init(Int_t netasegments,Int_t tv,AliRawEvent *rawevent,Float_t zvertex)
-{
- //Normal constructor
- fNEtaSegments = netasegments;
- fVersion = tv;
- fRawEvent = rawevent;
- fZVertex = zvertex;
-
- Init();
-}
-
-void AliHLTTPCHough::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 AliHLTTPCHough
- fBinary = binary;
- strcpy(fPath,path);
- fNEtaSegments = netasegments;
- fWriteDigits = kFALSE;
- fUse8bits = bit8;
- fVersion = tv;
- if(!fBinary) {
- if(infile) {
- fInputFile = infile;
- fInputPtr = 0;
- }
- else {
- fInputFile = 0;
- fInputPtr = ptr;
- }
- }
- else {
- fInputFile = 0;
- fInputPtr = 0;
- }
- fZVertex = zvertex;
-
- Init(); //do the rest
-}
-
-void AliHLTTPCHough::Init(Bool_t doit, Bool_t addhists)
-{
- // Init
- fDoIterative = doit;
- fAddHistograms = addhists;
-
- fNPatches = AliHLTTPCTransform::GetNPatches();
- fHoughTransformer = new AliHLTTPCHoughTransformer*[fNPatches];
- fMemHandler = new AliHLTTPCMemHandler*[fNPatches];
-
- fTracks = new AliHLTTPCTrackArray*[fNPatches];
- fEval = new AliHLTTPCHoughEval*[fNPatches];
-
- fGlobalTracks = new AliHLTTPCTrackArray("AliHLTTPCHoughTrack");
-
- AliHLTTPCHoughTransformer *lasttransformer = 0;
-
- for(Int_t i=0; i<fNPatches; i++)
- {
- switch (fVersion){ //choose Transformer
- case 1:
-#ifdef HAVE_ALIHLTHOUGHTRANSFORMERLUT
- fHoughTransformer[i] = new AliHLTHoughTransformerLUT(0,i,fNEtaSegments);
-#else
- AliErrorClassStream() << "AliHLTHoughTransformerLUT not compiled" << endl;
-#endif // HAVE_ALIHLTHOUGHTRANSFORMERLUT
- break;
- case 2:
-#ifdef HAVE_ALIHLTHOUGHCLUSTERTRANSFORMER
- fHoughTransformer[i] = new AliHLTHoughClusterTransformer(0,i,fNEtaSegments);
-#else
- AliErrorClassStream() << "AliHLTHoughClusterTransformer not compiled" << endl;
-#endif // HAVE_ALIHLTHOUGHCLUSTERTRANSFORMER
- break;
- case 3:
-#ifdef HAVE_ALIHLTHOUGHTRANSFORMERVHDL
- fHoughTransformer[i] = new AliHLTHoughTransformerVhdl(0,i,fNEtaSegments,fNSaveIterations);
-#else
- AliErrorClassStream() << "AliHLTHoughTransformerVhdl not compiled" << endl;
-#endif // HAVE_ALIHLTHOUGHTRANSFORMERVHDL
- break;
- case 4:
- fHoughTransformer[i] = new AliHLTTPCHoughTransformerRow(0,i,fNEtaSegments,kFALSE,fZVertex);
- break;
- default:
-#ifdef HAVE_ALIHLTHOUGHTRANSFORMERDEFAULT
- fHoughTransformer[i] = new AliHLTHoughTransformerDefault(0,i,fNEtaSegments,kFALSE,kFALSE);
-#else
- AliErrorClassStream() << "AliHLTHoughTransformerDefault not compiled" << endl;
-#endif // HAVE_ALIHLTHOUGHTRANSFORMERDEFAULT
- }
-
- fHoughTransformer[i]->SetLastTransformer(lasttransformer);
- lasttransformer = fHoughTransformer[i];
- // fHoughTransformer[i]->CreateHistograms(fNBinX[i],fLowPt[i],fNBinY[i],-fPhi[i],fPhi[i]);
- fHoughTransformer[i]->CreateHistograms(fNBinX[i],-fLowPt[i],fLowPt[i],fNBinY[i],-fPhi[i],fPhi[i]);
- //fHoughTransformer[i]->CreateHistograms(fLowPt[i],fUpperPt[i],fPtRes[i],fNBinY[i],fPhi[i]);
-
- fHoughTransformer[i]->SetLowerThreshold(fThreshold[i]);
- fHoughTransformer[i]->SetUpperThreshold(100);
-
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHough::Init","Version")
- <<"Initializing Hough transformer version "<<fVersion<<ENDLOG;
-
- fEval[i] = new AliHLTTPCHoughEval();
- fTracks[i] = new AliHLTTPCTrackArray("AliHLTTPCHoughTrack");
- if(fUse8bits) {
-#ifdef HAVE_ALIHLTDATAHANDLER
- fMemHandler[i] = new AliHLTDataHandler();
-#else //!HAVE_ALIHLTDATAHANDLER
- AliErrorClassStream() << "AliHLTDataHandler not compiled" << endl;
-#endif // HAVE_ALIHLTDATAHANDLER
- } else
- {
- if(!fRawEvent) {
- if(!fInputFile) {
- if(!fInputPtr) {
- /* In case of reading digits file */
- fMemHandler[i] = new AliHLTTPCFileHandler(kTRUE); //use static index
- if(!fBinary) {
- if(!fRunLoader) {
- Char_t filename[1024];
- sprintf(filename,"%s/digitfile.root",fPath);
- fMemHandler[i]->SetAliInput(filename);
- }
- else {
- fMemHandler[i]->SetAliInput(fRunLoader);
- }
- }
- }
- else {
- /* In case of reading from DATE */
-#ifdef HAVE_ALIHLTDDLDATAFILEHANDLER
- fMemHandler[i] = new AliHLTDDLDataFileHandler();
- fMemHandler[i]->SetReaderInput(fInputPtr,-1);
-#else //!HAVE_ALIHLTDDLDATAFILEHANDLER
- AliErrorClassStream() << "AliHLTDDLDataFileHandler not compiled" << endl;
-#endif //HAVE_ALIHLTDDLDATAFILEHANDLER
- }
- }
- else {
- /* In case of reading rawdata from ROOT file */
-#ifdef HAVE_ALIHLTDDLDATAFILEHANDLER
- fMemHandler[i] = new AliHLTDDLDataFileHandler();
- fMemHandler[i]->SetReaderInput(fInputFile);
-#else //!HAVE_ALIHLTDDLDATAFILEHANDLER
- AliErrorClassStream() << "AliHLTDDLDataFileHandler not compiled" << endl;
-#endif //HAVE_ALIHLTDDLDATAFILEHANDLER
- }
- }
- else {
- /* In case of reading rawdata using AliRawEvent */
-#ifdef HAVE_ALIHLTDDLDATAFILEHANDLER
- fMemHandler[i] = new AliHLTDDLDataFileHandler();
- fMemHandler[i]->SetReaderInput(fRawEvent);
-#else //!HAVE_ALIHLTDDLDATAFILEHANDLER
- AliErrorClassStream() << "AliHLTDDLDataFileHandler not compiled" << endl;
-#endif //HAVE_ALIHLTDDLDATAFILEHANDLER
- }
- }
- }
-
- fPeakFinder = new AliHLTTPCHoughMaxFinder("KappaPhi",50000);
- if(fVersion!=4) {
-#ifdef HAVE_ALIHLTHOUGHMERGER
- fMerger = new AliHLTHoughMerger(fNPatches);
-#else
- AliErrorClassStream() << "AliHLTHoughMerger not compiled" << endl;
-#endif //HAVE_ALIHLTHOUGHMERGER
-#ifdef HAVE_ALIHLTHOUGHINTMERGER
- fInterMerger = new AliHLTHoughIntMerger();
-#else
- AliErrorClassStream() << "AliHLTHoughIntMerger not compiled" << endl;
-#endif //HAVE_ALIHLTHOUGHINTMERGER
- }
- else {
- fMerger = 0;
- fInterMerger = 0;
- }
- fGlobalMerger = 0;
- fBenchmark = new AliHLTTPCBenchmark();
-}
-
-void AliHLTTPCHough::SetTransformerParams(Float_t ptres,Float_t ptmin,Float_t ptmax,Int_t ny,Int_t patch)
-{
- // Setup the parameters for the Hough Transformer
- // This includes the bin size and limits for
- // the parameter space histograms
-
- Int_t mrow;
- Float_t psi=0;
- if(patch==-1)
- mrow = 80;
- else
- mrow = AliHLTTPCTransform::GetLastRow(patch);
- if(ptmin)
- {
- Double_t lineradius = sqrt(pow(AliHLTTPCTransform::Row2X(mrow),2) + pow(AliHLTTPCTransform::GetMaxY(mrow),2));
- Double_t kappa = -1*AliHLTTPCTransform::GetBField()*AliHLTTPCTransform::GetBFact()/ptmin;
- psi = AliHLTTPCTransform::Deg2Rad(10) - asin(lineradius*kappa/2);
- cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
- }
-
- if(patch==-1)
- {
- Int_t i=0;
- while(i < 6)
- {
- fPtRes[i] = ptres;
- fLowPt[i] = ptmin;
- fUpperPt[i] = ptmax;
- fNBinY[i] = ny;
- fPhi[i] = psi;
- fNBinX[i]=0;
- i++;
- }
- return;
- }
-
- fPtRes[patch] = ptres;
- fLowPt[patch] = ptmin;
- fUpperPt[patch] = ptmax;
- fNBinY[patch] = ny;
- fPhi[patch] = psi;
-}
-/*
-void AliHLTTPCHough::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(AliHLTTPCTransform::Row2X(mrow),2) + pow(AliHLTTPCTransform::GetMaxY(mrow),2));
- Double_t kappa = -1*AliHLTTPCTransform::GetBField()*AliHLTTPCTransform::GetBFact()/ptmin;
- Double_t psi = AliHLTTPCTransform::Deg2Rad(10) - asin(lineradius*kappa/2);
- cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
-
- Int_t i=0;
- while(i < 6)
- {
- fLowPt[i] = ptmin;
- fNBinY[i] = ny;
- fNBinX[i] = nx;
- fPhi[i] = psi;
- i++;
- }
-}
-*/
-void AliHLTTPCHough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t /*patch*/)
-{
- // Setup the parameters for the Hough Transformer
-
- Double_t lineradius = 1.0/(AliHLTTPCHoughTransformerRow::GetBeta1()*sqrt(1.0+tan(AliHLTTPCTransform::Pi()*10/180)*tan(AliHLTTPCTransform::Pi()*10/180)));
- Double_t alpha1 = AliHLTTPCHoughTransformerRow::GetBeta1()*tan(AliHLTTPCTransform::Pi()*10/180);
- Double_t kappa = 1*AliHLTTPCTransform::GetBField()*AliHLTTPCTransform::GetBFact()/(ptmin*0.9);
- Double_t psi = AliHLTTPCTransform::Deg2Rad(10) - asin(lineradius*kappa/2);
- // cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
- Double_t alpha2 = alpha1 - (AliHLTTPCHoughTransformerRow::GetBeta1()-AliHLTTPCHoughTransformerRow::GetBeta2())*tan(psi);
- // cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<" in patch "<<patch<<endl;
-
- Int_t i=0;
- while(i < 6)
- {
- fLowPt[i] = 1.1*alpha1;
- fNBinY[i] = ny;
- fNBinX[i] = nx;
- fPhi[i] = alpha2;
- i++;
- }
-}
-
-void AliHLTTPCHough::CalcTransformerParams(Float_t ptmin)
-{
- // Setup the parameters for the Row Hough Transformer
- // Automatically adjusts the number of bins in X and Y in a way
- // that the size of the hough bin is 2x (in X) and 2.5 (in Y) the
- // size of the tpc pads
-
- Double_t lineradius = 1.0/(AliHLTTPCHoughTransformerRow::GetBeta1()*sqrt(1.0+tan(AliHLTTPCTransform::Pi()*10/180)*tan(AliHLTTPCTransform::Pi()*10/180)));
- Double_t alpha1 = AliHLTTPCHoughTransformerRow::GetBeta1()*tan(AliHLTTPCTransform::Pi()*10/180);
- Double_t kappa = 1*AliHLTTPCTransform::GetBField()*AliHLTTPCTransform::GetBFact()/(ptmin*0.9);
- Double_t psi = AliHLTTPCTransform::Deg2Rad(10) - asin(lineradius*kappa/2);
- // cout<<"Calculated psi range "<<psi<<endl;
- Double_t alpha2 = alpha1 - (AliHLTTPCHoughTransformerRow::GetBeta1()-AliHLTTPCHoughTransformerRow::GetBeta2())*tan(psi);
- alpha1 *= 1.1;
- // cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<endl;
-
- Double_t sizex = 2.0*AliHLTTPCTransform::GetPadPitchWidthLow()*AliHLTTPCHoughTransformerRow::GetBeta1()*AliHLTTPCHoughTransformerRow::GetBeta1();
- Double_t sizey = 2.5*AliHLTTPCTransform::GetPadPitchWidthUp()*AliHLTTPCHoughTransformerRow::GetBeta2()*AliHLTTPCHoughTransformerRow::GetBeta2();
-
- Int_t nx = 2*(Int_t)(alpha1/sizex)+1;
- Int_t ny = 2*(Int_t)(alpha2/sizey)+1;
- // cout<<"Calculated number of bins "<<nx<<" "<<ny<<endl;
-
- Int_t i=0;
- while(i < 6)
- {
- fLowPt[i] = alpha1;
- fNBinY[i] = ny;
- fNBinX[i] = nx;
- fPhi[i] = alpha2;
- i++;
- }
-}
-
-void AliHLTTPCHough::SetTransformerParams(Int_t nx,Int_t ny,Float_t lpt,Float_t phi)
-{
- // SetTransformerParams
-
- Int_t i=0;
- while(i < 6)
- {
- fLowPt[i] = lpt;
- fNBinY[i] = ny;
- fNBinX[i] = nx;
- fPhi[i] = phi;
- i++;
- }
-}
-
-void AliHLTTPCHough::SetThreshold(Int_t t3,Int_t patch)
-{
- // Set digits threshold
- if(patch==-1)
- {
- Int_t i=0;
- while(i < 6)
- fThreshold[i++]=t3;
- return;
- }
- fThreshold[patch]=t3;
-}
-
-void AliHLTTPCHough::SetPeakThreshold(Int_t threshold,Int_t patch)
-{
- // Set Peak Finder threshold
- if(patch==-1)
- {
- Int_t i=0;
- while(i < 6)
- fPeakThreshold[i++]=threshold;
- return;
- }
- fPeakThreshold[patch]=threshold;
-}
-
-void AliHLTTPCHough::DoBench(Char_t *name)
-{
- fBenchmark->Analyze(name);
-}
-
-void AliHLTTPCHough::Process(Int_t minslice,Int_t maxslice)
-{
- //Process all slices [minslice,maxslice].
-#ifdef HAVE_ALIHLTHOUGHGLOBALMERGER
- fGlobalMerger = new AliHLTHoughGlobalMerger(minslice,maxslice);
-#else
- return;
-#endif //HAVE_ALIHLTHOUGHGLOBALMERGER
-
- for(Int_t i=minslice; i<=maxslice; i++)
- {
- ReadData(i);
- Transform();
- if(fAddHistograms) {
- if(fVersion != 4)
- AddAllHistograms();
- else
- AddAllHistogramsRows();
- }
- FindTrackCandidates();
- //Evaluate();
- //fGlobalMerger->FillTracks(fTracks[0],i);
- }
-}
-
-void AliHLTTPCHough::ReadData(Int_t slice,Int_t eventnr)
-{
- //Read data from files, binary or root.
-
- if(fEvent!=eventnr) //just be sure that index is empty for new event
- AliHLTTPCFileHandler::CleanStaticIndex();
- fCurrentSlice = slice;
-
- for(Int_t i=0; i<fNPatches; i++)
- {
- fMemHandler[i]->Free();
- UInt_t ndigits=0;
- AliHLTTPCDigitRowData *digits =0;
- Char_t name[256];
- fMemHandler[i]->Init(slice,i);
- if(fBinary)//take input data from binary files
- {
- if(fUse8bits)
- sprintf(name,"%s/binaries/digits_c8_%d_%d_%d.raw",fPath,eventnr,slice,i);
- else
- sprintf(name,"%s/binaries/digits_%d_%d_%d.raw",fPath,eventnr,slice,i);
-
- fMemHandler[i]->SetBinaryInput(name);
- digits = (AliHLTTPCDigitRowData *)fMemHandler[i]->CompBinary2Memory(ndigits);
- fMemHandler[i]->CloseBinaryInput();
- }
- else //read data from root file
- {
- if(fEvent!=eventnr)
- fMemHandler[i]->FreeDigitsTree();//or else the new event is not loaded
- digits=(AliHLTTPCDigitRowData *)fMemHandler[i]->AliAltroDigits2Memory(ndigits,eventnr);
- }
-
- //Set the pointer to the TPCRawStream in case of fast raw data reading
- fHoughTransformer[i]->SetTPCRawStream(fMemHandler[i]->GetTPCRawStream());
-
- //set input data and init transformer
- fHoughTransformer[i]->SetInputData(ndigits,digits);
- fHoughTransformer[i]->Init(slice,i,fNEtaSegments);
- }
-
- fEvent=eventnr;
-}
-
-void AliHLTTPCHough::Transform(Int_t *rowrange)
-{
- //Transform all data given to the transformer within the given slice
- //(after ReadData(slice))
-
- Double_t initTime,cpuTime;
- initTime = GetCpuTime();
- Int_t patchorder[6] = {5,2,0,1,3,4}; //The order in which patches are processed
- // Int_t patchorder[6] = {0,1,2,3,4,5}; //The order in which patches are processed
- // Int_t patchorder[6] = {5,4,3,2,1,0}; //The order in which patches are processed
- // Int_t patchorder[6] = {5,2,4,3,1,0}; //The order in which patches are processed
- fLastPatch=-1;
- for(Int_t i=0; i<fNPatches; i++)
- {
- // In case of Row transformer reset the arrays only once
- if((fVersion != 4) || (i == 0)) {
- fBenchmark->Start("Hough Reset");
- fHoughTransformer[0]->Reset();//Reset the histograms
- fBenchmark->Stop("Hough Reset");
- }
- fBenchmark->Start("Hough Transform");
- PrepareForNextPatch(patchorder[i]);
- if(!rowrange) {
- char buf[256];
- sprintf(buf,"Patch %d",patchorder[i]);
- fBenchmark->Start(buf);
- fHoughTransformer[patchorder[i]]->SetLastPatch(fLastPatch);
- fHoughTransformer[patchorder[i]]->TransformCircle();
- fBenchmark->Stop(buf);
- }
- else
- fHoughTransformer[i]->TransformCircleC(rowrange,1);
- fBenchmark->Stop("Hough Transform");
- fLastPatch=patchorder[i];
- }
- cpuTime = GetCpuTime() - initTime;
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHough::Transform()","Timing")
- <<"Transform done in average per patch of "<<cpuTime*1000/fNPatches<<" ms"<<ENDLOG;
-}
-
-void AliHLTTPCHough::MergePatches()
-{
- // Merge patches if they are not summed
- if(fAddHistograms) //Nothing to merge here
- return;
- if (fMerger==NULL) return;
-#ifdef HAVE_ALIHLTHOUGHMERGER
- fMerger->MergePatches(kTRUE);
-#endif // HAVE_ALIHLTHOUGHMERGER
-}
-
-void AliHLTTPCHough::MergeInternally()
-{
- // Merge patches internally
- if (fMerger==NULL) return;
-#ifdef HAVE_ALIHLTHOUGHINTMERGER
- if(fAddHistograms)
- fInterMerger->FillTracks(fTracks[0]);
- else {
-#ifdef HAVE_ALIHLTHOUGHMERGER
- fInterMerger->FillTracks(fMerger->GetOutTracks());
-#endif // HAVE_ALIHLTHOUGHMERGER
- }
-
- fInterMerger->MMerge();
-#endif // HAVE_ALIHLTHOUGHINTMERGER
-}
-
-void AliHLTTPCHough::ProcessSliceIter()
-{
- //Process current slice (after ReadData(slice)) iteratively.
-
- if(!fAddHistograms)
- {
- if (fMerger==NULL) return;
- for(Int_t i=0; i<fNPatches; i++)
- {
- ProcessPatchIter(i);
-#ifdef HAVE_ALIHLTHOUGHMERGER
- fMerger->FillTracks(fTracks[i],i); //Copy tracks to merger
-#endif // HAVE_ALIHLTHOUGHMERGER
- }
- }
- else
- {
- for(Int_t i=0; i<10; i++)
- {
- Transform();
- AddAllHistograms();
- InitEvaluate();
- AliHLTTPCHoughTransformer *tr = fHoughTransformer[0];
- for(Int_t j=0; j<fNEtaSegments; j++)
- {
- AliHLTTPCHistogram *hist = tr->GetHistogram(j);
- if(hist->GetNEntries()==0) continue;
- fPeakFinder->Reset();
- fPeakFinder->SetHistogram(hist);
- fPeakFinder->FindAbsMaxima();
- AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)fTracks[0]->NextTrack();
- track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
- track->SetEtaIndex(j);
- track->SetEta(tr->GetEta(j,fCurrentSlice));
- for(Int_t k=0; k<fNPatches; k++)
- {
- fEval[i]->SetNumOfPadsToLook(2);
- fEval[i]->SetNumOfRowsToMiss(2);
- fEval[i]->RemoveFoundTracks();
- /*
- Int_t nrows=0;
- if(!fEval[i]->LookInsideRoad(track,nrows))
- {
- fTracks[0]->Remove(fTracks[0]->GetNTracks()-1);
- fTracks[0]->Compress();
- }
- */
- }
- }
-
- }
-
- }
-}
-
-void AliHLTTPCHough::ProcessPatchIter(Int_t patch)
-{
- //Process patch in a iterative way.
- //transform + peakfinding + evaluation + transform +...
-
- Int_t numoftries = 5;
- AliHLTTPCHoughTransformer *tr = fHoughTransformer[patch];
- AliHLTTPCTrackArray *tracks = fTracks[patch];
- tracks->Reset();
- AliHLTTPCHoughEval *ev = fEval[patch];
- ev->InitTransformer(tr);
- //ev->RemoveFoundTracks();
- ev->SetNumOfRowsToMiss(3);
- ev->SetNumOfPadsToLook(2);
- AliHLTTPCHistogram *hist;
- for(Int_t t=0; t<numoftries; t++)
- {
- tr->Reset();
- tr->TransformCircle();
- for(Int_t i=0; i<fNEtaSegments; i++)
- {
- hist = tr->GetHistogram(i);
- if(hist->GetNEntries()==0) continue;
- fPeakFinder->Reset();
- fPeakFinder->SetHistogram(hist);
- fPeakFinder->FindAbsMaxima();
- //fPeakFinder->FindPeak1();
- AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)tracks->NextTrack();
- track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
- track->SetEtaIndex(i);
- track->SetEta(tr->GetEta(i,fCurrentSlice));
- /*
- Int_t nrows=0;
- if(!ev->LookInsideRoad(track,nrows))
- {
- tracks->Remove(tracks->GetNTracks()-1);
- tracks->Compress();
- }
- */
- }
- }
- fTracks[0]->QSort();
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHough::ProcessPatch","NTracks")
- <<AliHLTTPCLog::kDec<<"Found "<<tracks->GetNTracks()<<" tracks in patch "<<patch<<ENDLOG;
-}
-
-void AliHLTTPCHough::AddAllHistograms()
-{
- //Add the histograms within one etaslice.
- //Resulting histogram are in patch=0.
-
- Double_t initTime,cpuTime;
- initTime = GetCpuTime();
- fBenchmark->Start("Add Histograms");
- for(Int_t i=0; i<fNEtaSegments; i++)
- {
- AliHLTTPCHistogram *hist0 = fHoughTransformer[0]->GetHistogram(i);
- for(Int_t j=1; j<fNPatches; j++)
- {
- AliHLTTPCHistogram *hist = fHoughTransformer[j]->GetHistogram(i);
- hist0->Add(hist);
- }
- }
- fBenchmark->Stop("Add Histograms");
- fAddHistograms = kTRUE;
- cpuTime = GetCpuTime() - initTime;
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHough::AddAllHistograms()","Timing")
- <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
-}
-
-void AliHLTTPCHough::AddAllHistogramsRows()
-{
- //Add the histograms within one etaslice.
- //Resulting histogram are in patch=0.
-
- Double_t initTime,cpuTime;
- initTime = GetCpuTime();
- fBenchmark->Start("Add HistogramsRows");
-
- UChar_t lastpatchlastrow = AliHLTTPCTransform::GetLastRowOnDDL(fLastPatch)+1;
-
- UChar_t *tracklastrow = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetTrackLastRow();
-
- for(Int_t i=0; i<fNEtaSegments; i++)
- {
- UChar_t *gapcount = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetGapCount(i);
- UChar_t *currentrowcount = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetCurrentRowCount(i);
-
- AliHLTTPCHistogram *hist = fHoughTransformer[0]->GetHistogram(i);
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
- Int_t nxbins = hist->GetNbinsX()+2;
-
- for(Int_t ybin=ymin; ybin<=ymax; ybin++)
- {
- for(Int_t xbin=xmin; xbin<=xmax; xbin++)
- {
- Int_t bin = xbin + ybin*nxbins; //Int_t bin = hist->GetBin(xbin,ybin);
- if(gapcount[bin] < MAX_N_GAPS) {
- if(tracklastrow[bin] > lastpatchlastrow) {
- if(lastpatchlastrow > currentrowcount[bin])
- gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
- }
- else {
- if(tracklastrow[bin] > currentrowcount[bin])
- gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
- }
- if(gapcount[bin] < MAX_N_GAPS)
- hist->AddBinContent(bin,(159-gapcount[bin]));
- }
- }
- }
- }
-
- fBenchmark->Stop("Add HistogramsRows");
- fAddHistograms = kTRUE;
- cpuTime = GetCpuTime() - initTime;
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHough::AddAllHistogramsRows()","Timing")
- <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
-}
-
-void AliHLTTPCHough::PrepareForNextPatch(Int_t nextpatch)
-{
- // Prepare the parameter space for the processing of
- // the next read patch. According to the already
- // accumulated number of gaps in parameter space
- // bins, the routine updates the dynamic
- // pointers used in order to jump rapidly during the
- // filling of the parameter space.
-
- char buf[256];
- sprintf(buf,"Prepare For Patch %d",nextpatch);
- fBenchmark->Start(buf);
-
- UChar_t lastpatchlastrow;
- if(fLastPatch == -1)
- lastpatchlastrow = 0;
- else
- lastpatchlastrow = AliHLTTPCTransform::GetLastRowOnDDL(fLastPatch)+1;
- UChar_t nextpatchfirstrow;
- if(nextpatch==0)
- nextpatchfirstrow = 0;
- else
- nextpatchfirstrow = AliHLTTPCTransform::GetFirstRowOnDDL(nextpatch)-1;
-
- UChar_t *trackfirstrow = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetTrackFirstRow();
- UChar_t *tracklastrow = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetTrackLastRow();
-
- for(Int_t i=0; i<fNEtaSegments; i++)
- {
- UChar_t *gapcount = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetGapCount(i);
- UChar_t *currentrowcount = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetCurrentRowCount(i);
- UChar_t *prevbin = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetPrevBin(i);
- UChar_t *nextbin = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetNextBin(i);
- UChar_t *nextrow = ((AliHLTTPCHoughTransformerRow *)fHoughTransformer[0])->GetNextRow(i);
-
- AliHLTTPCHistogram *hist = fHoughTransformer[0]->GetHistogram(i);
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
- Int_t nxbins = hist->GetNbinsX()+2;
-
- if(fLastPatch != -1) {
- UChar_t lastyvalue = 0;
- Int_t endybin = ymin - 1;
- for(Int_t ybin=nextrow[ymin]; ybin<=ymax; ybin = nextrow[++ybin])
- {
- UChar_t lastxvalue = 0;
- UChar_t maxvalue = 0;
- Int_t endxbin = xmin - 1;
- for(Int_t xbin=xmin; xbin<=xmax; xbin++)
- {
- Int_t bin = xbin + ybin*nxbins;
- UChar_t value = 0;
- if(gapcount[bin] < MAX_N_GAPS) {
- if(tracklastrow[bin] > lastpatchlastrow) {
- if(lastpatchlastrow > currentrowcount[bin])
- gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
- }
- else {
- if(tracklastrow[bin] > currentrowcount[bin])
- gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
- }
- if(gapcount[bin] < MAX_N_GAPS) {
- value = 1;
- maxvalue = 1;
- if(trackfirstrow[bin] < nextpatchfirstrow)
- currentrowcount[bin] = nextpatchfirstrow;
- else
- currentrowcount[bin] = trackfirstrow[bin];
- }
- }
- if(value > 0)
- {
- nextbin[xbin + ybin*nxbins] = (UChar_t)xbin;
- prevbin[xbin + ybin*nxbins] = (UChar_t)xbin;
- if(value > lastxvalue)
- {
- UChar_t *tempnextbin = nextbin + endxbin + 1 + ybin*nxbins;
- memset(tempnextbin,(UChar_t)xbin,xbin-endxbin-1);
- }
- endxbin = xbin;
- }
- else
- {
- prevbin[xbin + ybin*nxbins] = (UChar_t)endxbin;
- }
- lastxvalue = value;
- }
- UChar_t *tempnextbin = nextbin + endxbin + 1 + ybin*nxbins;
- memset(tempnextbin,(UChar_t)(xmax+1),xmax-endxbin);
- if(maxvalue > 0)
- {
- nextrow[ybin] = (UChar_t)ybin;
- if(maxvalue > lastyvalue)
- {
- UChar_t *tempnextrow = nextrow + endybin + 1;
- memset(tempnextrow,(UChar_t)ybin,ybin-endybin-1);
- }
- endybin = ybin;
- }
- lastyvalue = maxvalue;
- }
- UChar_t *tempnextrow = nextrow + endybin + 1;
- memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin+1);
- }
- else {
- UChar_t lastyvalue = 0;
- Int_t endybin = ymin - 1;
- for(Int_t ybin=ymin; ybin<=ymax; ybin++)
- {
- UChar_t maxvalue = 0;
- for(Int_t xbin=xmin; xbin<=xmax; xbin++)
- {
- Int_t bin = xbin + ybin*nxbins;
- if(gapcount[bin] < MAX_N_GAPS) {
- maxvalue = 1;
- if(trackfirstrow[bin] < nextpatchfirstrow)
- currentrowcount[bin] = nextpatchfirstrow;
- else
- currentrowcount[bin] = trackfirstrow[bin];
- }
- }
- if(maxvalue > 0)
- {
- nextrow[ybin] = (UChar_t)ybin;
- if(maxvalue > lastyvalue)
- {
- UChar_t *tempnextrow = nextrow + endybin + 1;
- memset(tempnextrow,(UChar_t)ybin,ybin-endybin-1);
- }
- endybin = ybin;
- }
- lastyvalue = maxvalue;
- }
- UChar_t *tempnextrow = nextrow + endybin + 1;
- memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin+1);
- }
- }
-
- fBenchmark->Stop(buf);
-}
-
-void AliHLTTPCHough::AddTracks()
-{
- // Add current slice slice tracks to the global list of found tracks
- if(!fTracks[0])
- {
- cerr<<"AliHLTTPCHough::AddTracks : No tracks"<<endl;
- return;
- }
- AliHLTTPCTrackArray *tracks = fTracks[0];
- for(Int_t i=0; i<tracks->GetNTracks(); i++)
- {
- AliHLTTPCTrack *track = tracks->GetCheckedTrack(i);
- if(!track) continue;
- if(track->GetNHits()!=1) cerr<<"NHITS "<<track->GetNHits()<<endl;
- UInt_t *ids = track->GetHitNumbers();
- ids[0] = (fCurrentSlice&0x7f)<<25;
- }
-
- fGlobalTracks->AddTracks(fTracks[0],0,fCurrentSlice);
-}
-
-void AliHLTTPCHough::FindTrackCandidatesRow()
-{
- // Find AliHLTTPCHoughTransformerRow track candidates
- if(fVersion != 4) {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHough::FindTrackCandidatesRow()","")
- <<"Incompatible Peak Finder version!"<<ENDLOG;
- return;
- }
-
- //Look for peaks in histograms, and find the track candidates
- Int_t npatches;
- if(fAddHistograms)
- npatches = 1; //Histograms have been added.
- else
- npatches = fNPatches;
-
- Double_t initTime,cpuTime;
- initTime = GetCpuTime();
- fBenchmark->Start("Find Maxima");
- for(Int_t i=0; i<npatches; i++)
- {
- AliHLTTPCHoughTransformer *tr = fHoughTransformer[i];
- AliHLTTPCHistogram *h = tr->GetHistogram(0);
- Float_t deltax = h->GetBinWidthX()*AliHLTTPCHoughTransformerRow::GetDAlpha();
- Float_t deltay = h->GetBinWidthY()*AliHLTTPCHoughTransformerRow::GetDAlpha();
- Float_t deltaeta = (tr->GetEtaMax()-tr->GetEtaMin())/tr->GetNEtaSegments()*AliHLTTPCHoughTransformerRow::GetDEta();
- Float_t zvertex = tr->GetZVertex();
- fTracks[i]->Reset();
- fPeakFinder->Reset();
-
- for(Int_t j=0; j<fNEtaSegments; j++)
- {
- AliHLTTPCHistogram *hist = tr->GetHistogram(j);
- if(hist->GetNEntries()==0) continue;
- fPeakFinder->SetHistogram(hist);
- fPeakFinder->SetEtaSlice(j);
- fPeakFinder->SetTrackLUTs(((AliHLTTPCHoughTransformerRow *)tr)->GetTrackNRows(),((AliHLTTPCHoughTransformerRow *)tr)->GetTrackFirstRow(),((AliHLTTPCHoughTransformerRow *)tr)->GetTrackLastRow(),((AliHLTTPCHoughTransformerRow *)tr)->GetNextRow(j));
-#ifdef do_mc
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHough::FindTrackCandidates()","")
- <<"Starting "<<j<<" etaslice"<<ENDLOG;
-#endif
- fPeakFinder->SetThreshold(fPeakThreshold[i]);
- fPeakFinder->FindAdaptedRowPeaks(1,0,0);//Maxima finder for HoughTransformerRow
- }
-
- for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
- {
- // if(fPeakFinder->GetWeight(k) < 0) continue;
- AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)fTracks[i]->NextTrack();
- Double_t starteta = tr->GetEta(fPeakFinder->GetStartEta(k),fCurrentSlice);
- Double_t endeta = tr->GetEta(fPeakFinder->GetEndEta(k),fCurrentSlice);
- Double_t eta = (starteta+endeta)/2.0;
- track->SetTrackParametersRow(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),eta,fPeakFinder->GetWeight(k));
- track->SetPterr(deltax); track->SetPsierr(deltay); track->SetTglerr(deltaeta);
- track->SetBinXY(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetXPeakSize(k),fPeakFinder->GetYPeakSize(k));
- track->SetZ0(zvertex);
- Int_t etaindex = (fPeakFinder->GetStartEta(k)+fPeakFinder->GetEndEta(k))/2;
- track->SetEtaIndex(etaindex);
- Int_t rows[2];
- ((AliHLTTPCHoughTransformerRow *)tr)->GetTrackLength(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),rows);
- track->SetRowRange(rows[0],rows[1]);
- track->SetSector(fCurrentSlice);
- track->SetSlice(fCurrentSlice);
-#ifdef do_mc
- Int_t label = tr->GetTrackID(etaindex,fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k));
- track->SetMCid(label);
-#endif
- }
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHough::FindTrackCandidates()","")
- <<"Found "<<fTracks[i]->GetNTracks()<<" tracks in slice "<<fCurrentSlice<<ENDLOG;
- fTracks[i]->QSort();
- }
- fBenchmark->Stop("Find Maxima");
- cpuTime = GetCpuTime() - initTime;
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHough::FindTrackCandidates()","Timing")
- <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
-}
-
-void AliHLTTPCHough::FindTrackCandidates()
-{
- // Find AliHLTTPCHoughTransformer track candidates
- if(fVersion == 4) {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHough::FindTrackCandidatesRow()","")
- <<"Incompatible Peak Finder version!"<<ENDLOG;
- return;
- }
-
- Int_t npatches;
- if(fAddHistograms)
- npatches = 1; //Histograms have been added.
- else
- npatches = fNPatches;
-
- Double_t initTime,cpuTime;
- initTime = GetCpuTime();
- fBenchmark->Start("Find Maxima");
- for(Int_t i=0; i<npatches; i++)
- {
- AliHLTTPCHoughTransformer *tr = fHoughTransformer[i];
- fTracks[i]->Reset();
-
- for(Int_t j=0; j<fNEtaSegments; j++)
- {
- AliHLTTPCHistogram *hist = tr->GetHistogram(j);
- if(hist->GetNEntries()==0) continue;
- fPeakFinder->Reset();
- fPeakFinder->SetHistogram(hist);
-#ifdef do_mc
- cout<<"Starting "<<j<<" etaslice"<<endl;
-#endif
- fPeakFinder->SetThreshold(fPeakThreshold[i]);
- fPeakFinder->FindAdaptedPeaks(fKappaSpread,fPeakRatio);
-
- for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
- {
- AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)fTracks[i]->NextTrack();
- track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
- track->SetEtaIndex(j);
- track->SetEta(tr->GetEta(j,fCurrentSlice));
- track->SetRowRange(AliHLTTPCTransform::GetFirstRow(0),AliHLTTPCTransform::GetLastRow(5));
- }
- }
- cout<<"Found "<<fTracks[i]->GetNTracks()<<" tracks in patch "<<i<<endl;
- fTracks[i]->QSort();
- }
- fBenchmark->Stop("Find Maxima");
- cpuTime = GetCpuTime() - initTime;
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHough::FindTrackCandidates()","Timing")
- <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
-}
-
-void AliHLTTPCHough::InitEvaluate()
-{
- //Pass the transformer objects to the AliHLTTPCHoughEval objects:
- //This will provide the evaluation objects with all the necessary
- //data and parameters it needs.
-
- for(Int_t i=0; i<fNPatches; i++)
- fEval[i]->InitTransformer(fHoughTransformer[i]);
-}
-
-Int_t AliHLTTPCHough::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.
- //If histograms were not added, the check is done locally in patch,
- //meaning that nrowstomiss is the number of padrows the road can miss with respect
- //to the number of rows in the patch.
- //If the histograms were added, the comparison is done globally in the _slice_,
- //meaing that nrowstomiss is the number of padrows the road can miss with
- //respect to the total number of padrows in the slice.
- //
- //Return value = number of tracks which were removed (only in case of fAddHistograms)
-
- if(!fTracks[0])
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHough::Evaluate","Track Array")
- <<"No tracks to work with..."<<ENDLOG;
- return 0;
- }
-
- Int_t removedtracks=0;
- AliHLTTPCTrackArray *tracks=0;
-
- if(fAddHistograms)
- {
- tracks = fTracks[0];
- for(Int_t i=0; i<tracks->GetNTracks(); i++)
- {
- AliHLTTPCTrack *track = tracks->GetCheckedTrack(i);
- if(!track) continue;
- track->SetNHits(0);
- }
- }
-
- for(Int_t i=0; i<fNPatches; i++)
- EvaluatePatch(i,roadwidth,nrowstomiss);
-
- //Here we check the tracks globally;
- //how many good rows (padrows with signal)
- //did it cross in the slice
- if(fAddHistograms)
- {
- for(Int_t j=0; j<tracks->GetNTracks(); j++)
- {
- AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)tracks->GetCheckedTrack(j);
-
- if(track->GetNHits() < AliHLTTPCTransform::GetNRows() - nrowstomiss)
- {
- tracks->Remove(j);
- removedtracks++;
- }
- }
- tracks->Compress();
- tracks->QSort();
- }
-
- return removedtracks;
-}
-
-void AliHLTTPCHough::EvaluatePatch(Int_t i,Int_t roadwidth,Int_t nrowstomiss)
-{
- //Evaluate patch i.
-
- fEval[i]->InitTransformer(fHoughTransformer[i]);
- fEval[i]->SetNumOfPadsToLook(roadwidth);
- fEval[i]->SetNumOfRowsToMiss(nrowstomiss);
- //fEval[i]->RemoveFoundTracks();
-
- AliHLTTPCTrackArray *tracks=0;
-
- if(!fAddHistograms)
- tracks = fTracks[i];
- else
- tracks = fTracks[0];
-
- Int_t nrows=0;
- for(Int_t j=0; j<tracks->GetNTracks(); j++)
- {
- AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)tracks->GetCheckedTrack(j);
- if(!track)
- {
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHough::EvaluatePatch","Track array")
- <<"Track object missing!"<<ENDLOG;
- continue;
- }
- nrows=0;
- Int_t rowrange[2] = {AliHLTTPCTransform::GetFirstRow(i),AliHLTTPCTransform::GetLastRow(i)};
- Bool_t result = fEval[i]->LookInsideRoad(track,nrows,rowrange);
- if(fAddHistograms)
- {
- Int_t pre=track->GetNHits();
- track->SetNHits(pre+nrows);
- }
- else//the track crossed too few good padrows (padrows with signal) in the patch, so remove it
- {
- if(result == kFALSE)
- tracks->Remove(j);
- }
- }
-
- tracks->Compress();
-
-}
-
-void AliHLTTPCHough::MergeEtaSlices()
-{
- //Merge tracks found in neighbouring eta slices.
- //Removes the track with the lower weight.
-
- fBenchmark->Start("Merge Eta-slices");
- AliHLTTPCTrackArray *tracks = fTracks[0];
- if(!tracks)
- {
- cerr<<"AliHLTTPCHough::MergeEtaSlices : No tracks "<<endl;
- return;
- }
- for(Int_t j=0; j<tracks->GetNTracks(); j++)
- {
- AliHLTTPCHoughTrack *track1 = (AliHLTTPCHoughTrack*)tracks->GetCheckedTrack(j);
- if(!track1) continue;
- for(Int_t k=j+1; k<tracks->GetNTracks(); k++)
- {
- AliHLTTPCHoughTrack *track2 = (AliHLTTPCHoughTrack*)tracks->GetCheckedTrack(k);
- if(!track2) continue;
- if(abs(track1->GetEtaIndex() - track2->GetEtaIndex()) != 1) continue;
- if(fabs(track1->GetKappa()-track2->GetKappa()) < 0.006 &&
- fabs(track1->GetPsi()- track2->GetPsi()) < 0.1)
- {
- //cout<<"Merging track in slices "<<track1->GetEtaIndex()<<" "<<track2->GetEtaIndex()<<endl;
- if(track1->GetWeight() > track2->GetWeight())
- tracks->Remove(k);
- else
- tracks->Remove(j);
- }
- }
- }
- fBenchmark->Stop("Merge Eta-slices");
- tracks->Compress();
-}
-
-void AliHLTTPCHough::WriteTracks(Char_t *path)
-{
- // Write found tracks into file
- //cout<<"AliHLTTPCHough::WriteTracks : Sorting the tracsk"<<endl;
- //fGlobalTracks->QSort();
-
- Char_t filename[1024];
- sprintf(filename,"%s/tracks_%d.raw",path,fEvent);
- AliHLTTPCMemHandler mem;
- mem.SetBinaryOutput(filename);
- mem.TrackArray2Binary(fGlobalTracks);
- mem.CloseBinaryOutput();
- fGlobalTracks->Reset();
-}
-
-void AliHLTTPCHough::WriteTracks(Int_t slice,Char_t *path)
-{
- // Write found tracks slice by slice into file
-
- AliHLTTPCMemHandler mem;
- Char_t fname[100];
- if(fAddHistograms)
- {
- sprintf(fname,"%s/tracks_ho_%d_%d.raw",path,fEvent,slice);
- mem.SetBinaryOutput(fname);
- mem.TrackArray2Binary(fTracks[0]);
- mem.CloseBinaryOutput();
- }
- else
- {
- for(Int_t i=0; i<fNPatches; i++)
- {
- sprintf(fname,"%s/tracks_ho_%d_%d_%d.raw",path,fEvent,slice,i);
- mem.SetBinaryOutput(fname);
- mem.TrackArray2Binary(fTracks[i]);
- mem.CloseBinaryOutput();
- }
- }
-}
-
-Int_t AliHLTTPCHough::FillESD(AliESDEvent *esd)
-{
- // Fill the found hough transform tracks
- // into the ESD. The tracks are stored as
- // AliESDHLTtrack objects.
-
- // No TPC PID so far,assuming pions
- Double_t prob[AliPID::kSPECIES];
- for(Int_t i=0;i<AliPID::kSPECIES;i++) {
- if(i==AliPID::kPion) prob[i]=1.0;
- else prob[i]=0.1;
- }
-
- if(!fGlobalTracks) return 0;
- Int_t nglobaltracks = 0;
- for(Int_t i=0; i<fGlobalTracks->GetNTracks(); i++)
- {
- AliHLTTPCHoughTrack *tpt = (AliHLTTPCHoughTrack *)fGlobalTracks->GetCheckedTrack(i);
- if(!tpt) continue;
-
- if(tpt->GetWeight()<0) continue;
- AliHLTTPCHoughKalmanTrack *tpctrack = new AliHLTTPCHoughKalmanTrack(*tpt);
- if(!tpctrack) continue;
- AliESDtrack *esdtrack2 = new AliESDtrack() ;
- esdtrack2->UpdateTrackParams(tpctrack,AliESDtrack::kTPCin);
- esdtrack2->SetESDpid(prob);
- esd->AddTrack(esdtrack2);
- nglobaltracks++;
- delete esdtrack2;
- delete tpctrack;
- }
- return nglobaltracks;
-}
-
-void AliHLTTPCHough::WriteDigits(Char_t *outfile)
-{
- //Write the current data to a new rootfile.
- for(Int_t i=0; i<fNPatches; i++)
- {
- AliHLTTPCDigitRowData *tempPt = (AliHLTTPCDigitRowData*)fHoughTransformer[i]->GetDataPointer();
- fMemHandler[i]->AliDigits2RootFile(tempPt,outfile);
- }
-}
-
-Double_t AliHLTTPCHough::GetCpuTime()
-{
- //Return the Cputime in seconds.
- struct timeval tv;
- gettimeofday( &tv, NULL );
- return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
-}
-
-void *AliHLTTPCHough::ProcessInThread(void *args)
-{
- // Called in case Hough transform tracking
- // is executed in a thread
-
- AliHLTTPCHough *instance = (AliHLTTPCHough *)args;
- Int_t minslice = instance->GetMinSlice();
- Int_t maxslice = instance->GetMaxSlice();
- for(Int_t i=minslice; i<=maxslice; i++)
- {
- instance->ReadData(i,0);
- instance->Transform();
- instance->AddAllHistogramsRows();
- instance->FindTrackCandidatesRow();
- instance->AddTracks();
- }
- return (void *)0;
-}
-
-void AliHLTTPCHough::StartProcessInThread(Int_t minslice,Int_t maxslice)
-{
- // Starts the Hough transform tracking as a
- // separate thread. Takes as parameters the
- // range of TPC slices (sectors) to be reconstructed
-
-#ifdef HAVE_THREAD
- if(!fThread) {
- char buf[255];
- sprintf(buf,"houghtrans_%d_%d",minslice,maxslice);
- SetMinMaxSlices(minslice,maxslice);
- // fThread = new TThread(buf,(void (*) (void *))&ProcessInThread,(void *)this);
- fThread = new TThread(buf,&ProcessInThread,(void *)this);
- fThread->Run();
- }
-#else // HAVE_THREAD
- AliErrorClassStream() << "thread support not compiled" << endl;
-#endif // HAVE_THREAD
- return;
-}
-
-Int_t AliHLTTPCHough::WaitForThreadFinish()
-{
- // Routine is used in case we run the
- // Hough transform tracking in several
- // threads and want to sync them before
- // writing the results to the ESD
-
-#ifdef HAVE_THREAD
-#if ROOT_VERSION_CODE < 262403
- return TThread::Join(fThread->GetId());
-#else
- return fThread->Join(fThread->GetId());
-#endif
-#endif // HAVE_THREAD
- AliErrorClassStream() << "thread support not compiled" << endl;
- return 0;
-}
+++ /dev/null
-//-*- Mode: C++ -*-
-// $Id$
-// origin hough/AliL3Hough.h,v 1.31 Fri Feb 25 07:32:13 2005 UTC by cvetan
-
-#ifndef ALIHLTTPCHOUGH_H
-#define ALIHLTTPCHOUGH_H
-
-//* This file is property of and copyright by the ALICE HLT Project *
-//* ALICE Experiment at CERN, All rights reserved. *
-//* See cxx source for full Copyright notice *
-
-/** @file AliHLTTPCHough.h
- @author Anders Vestbo, Cvetan Cheshkov
- @date
- @brief Steering for HLT TPC hough transform tracking algorithms. */
-
-#include "AliHLTTPCRootTypes.h"
-
-class AliHLTTPCHoughMaxFinder;
-class AliHLTTPCHoughTransformer;
-class AliHLTTPCHistogram;
-class AliHLTTPCMemHandler;
-class AliHLTTPCFileHandler;
-class AliHLTTPCHoughEval;
-class AliHLTTPCTrackArray;
-class AliHLTTPCHoughMerger;
-class AliHLTTPCHoughIntMerger;
-class AliHLTTPCHoughGlobalMerger;
-class AliHLTTPCBenchmark;
-
-#ifdef HAVE_THREAD
-class TThread;
-#endif // HAVE_THREAD
-class AliRunLoader;
-class AliRawEvent;
-class AliESDEvent;
-class AliESDHLTtrack;
-
-/**
- * @class AliHLTTPCHough
- * Interface class for the HLT TPC Hough transform tracking algorithms
- *
- * Example how to use:
- *<pre>
- * AliHLTTPCHough *hough = new AliHLTTPCHough(path,kTRUE,NumberOfEtaSegments);
- * hough->ReadData(slice);
- * hough->Transform();
- * hough->FindTrackCandidates();
- *
- * AliHLTTPCTrackArray *tracks = hough->GetTracks(patch);
- *
- *</pre>
-*/
-class AliHLTTPCHough {
- public:
-
- AliHLTTPCHough();
- AliHLTTPCHough(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 ~AliHLTTPCHough();
-
- void SetRunLoader(AliRunLoader *runloader) {fRunLoader = runloader;}
-
- void Init(Int_t netasegments,Int_t tv,AliRawEvent *rawevent,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 *rowrange = 0);
- void ProcessSliceIter();
- void ProcessPatchIter(Int_t patch);
- void MergePatches();
- void MergeInternally();
- void MergeEtaSlices();
-
- void FindTrackCandidates();
- void FindTrackCandidatesRow();
- void AddAllHistograms();
- void AddAllHistogramsRows();
- void PrepareForNextPatch(Int_t nextpatch);
- 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);
- Int_t FillESD(AliESDEvent *esd);
- void WriteDigits(Char_t *outfile="output_digits.root");
- void InitEvaluate();
- void DoBench(Char_t *filename);
- void AddTracks();
-
- //Setters
- void SetNEtaSegments(Int_t i) {fNEtaSegments = i;}
- void SetAddHistograms() {fAddHistograms = kTRUE;}
- void DoIterative() {fDoIterative = kTRUE;}
- void SetWriteDigits() {fWriteDigits = kTRUE;}
- void SetTransformerParams(Float_t ptres=0,Float_t ptmin=0,Float_t ptmax=0,Int_t ny=0,Int_t patch=-1);
- //{fPtRes=ptres;fNBinY=ny;fLowPt=ptmin;fUpperPt=ptmax;fPhi=psi;}
- void SetTransformerParams(Int_t nx,Int_t ny,Float_t lpt,Int_t patch);
- void CalcTransformerParams(Float_t lpt);
- void SetTransformerParams(Int_t nx,Int_t ny,Float_t lpt,Float_t phi);
- //{fNBinX=nx;fNBinY=ny;fLowPt=lpt;fPhi=phi;}
- void SetThreshold(Int_t t=3,Int_t patch=-1);
- void SetNSaveIterations(Int_t t=10) {fNSaveIterations=t;}
- void SetPeakThreshold(Int_t threshold=0,Int_t patch=-1);
-
- void SetPeakParameters(Int_t kspread,Float_t pratio) {fKappaSpread=kspread; fPeakRatio=pratio;}
-
- //Getters
- AliHLTTPCHoughTransformer *GetTransformer(Int_t i) const {if(!fHoughTransformer[i]) return 0; return fHoughTransformer[i];}
- AliHLTTPCTrackArray *GetTracks(Int_t i) const {if(!fTracks[i]) return 0; return fTracks[i];}
- AliHLTTPCHoughEval *GetEval(Int_t i) const {if(!fEval[i]) return 0; return fEval[i];}
- AliHLTTPCHoughMerger *GetMerger() const{if(!fMerger) return 0; return fMerger;}
- AliHLTTPCHoughIntMerger *GetInterMerger() const {if(!fInterMerger) return 0; return fInterMerger;}
- AliHLTTPCMemHandler *GetMemHandler(Int_t i) const {if(!fMemHandler[i]) return 0; return fMemHandler[i];}
- AliHLTTPCHoughMaxFinder *GetMaxFinder() const {return fPeakFinder;}
-
- //Special methods for executing Hough Transform as a thread
- static void *ProcessInThread(void *args);
- void StartProcessInThread(Int_t minslice,Int_t maxslice);
- Int_t WaitForThreadFinish();
- void SetMinMaxSlices(Int_t minslice,Int_t maxslice) {fMinSlice = minslice; fMaxSlice = maxslice;}
- Int_t GetMinSlice() const {return fMinSlice;}
- Int_t GetMaxSlice() const {return fMaxSlice;}
-
- private:
- /// copy constructor not permitted
- AliHLTTPCHough(const AliHLTTPCHough);
- /// assignment operator not permitted
- AliHLTTPCHough& operator=(const AliHLTTPCHough);
-
- Char_t *fInputFile;//!
- Char_t *fInputPtr;//!
- AliRawEvent *fRawEvent;//!
- 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 fLastPatch; //The index of the last processed patch
- Int_t fVersion; //which HoughTransformer to use
- Int_t fCurrentSlice; // Current TPC slice (sector)
- 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
-
- Int_t fMinSlice; // First TPC slice (sector) to process while running in a thread
- Int_t fMaxSlice; // Last TPC slice (sector) to process while running in a thread
-
- AliHLTTPCMemHandler **fMemHandler; //!
- AliHLTTPCHoughTransformer **fHoughTransformer; //!
- AliHLTTPCHoughEval **fEval; //!
- AliHLTTPCHoughMaxFinder *fPeakFinder; //!
- AliHLTTPCTrackArray **fTracks; //!
- AliHLTTPCTrackArray *fGlobalTracks; //!
- AliHLTTPCHoughMerger *fMerger; //!
- AliHLTTPCHoughIntMerger *fInterMerger; //!
- AliHLTTPCHoughGlobalMerger *fGlobalMerger; //!
- AliHLTTPCBenchmark *fBenchmark; //!
-
- AliRunLoader *fRunLoader; // Run Loader
-
- void CleanUp();
- Double_t GetCpuTime();
-
-#ifdef HAVE_THREAD
- TThread *fThread; //! Pointer to the TThread object in case of running in a thread
-#endif // HAVE_THREAD
-
- ClassDef(AliHLTTPCHough,1) //Hough transform base class
-};
-
-#endif
+++ /dev/null
-// @(#) $Id$
-// origin: hough/AliL3HoughEval.cxx,v 1.28 Thu Jun 17 10:36:14 2004 UTC by cvetan
-
-// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
-//*-- Copyright © ALICE HLT Group
-
-#include "AliHLTStdIncludes.h"
-
-#include <TH1.h>
-#include <TFile.h>
-
-#include "AliHLTTPCLogging.h"
-#include "AliHLTTPCHoughEval.h"
-#include "AliHLTTPCMemHandler.h"
-#include "AliHLTTPCTrackArray.h"
-#include "AliHLTTPCHoughTransformer.h"
-#include "AliHLTTPCDigitData.h"
-#include "AliHLTTPCHoughTrack.h"
-#include "AliHLTTPCTransform.h"
-#include "AliHLTTPCHistogram.h"
-#include "AliHLTTPCHistogram1D.h"
-
-#if __GNUC__ == 3
-using namespace std;
-#endif
-
-/** /class AliHLTTPCHoughEval
-//<pre>
-//_____________________________________________________________
-// AliHLTTPCHoughEval
-//
-// Evaluation class for tracklets produced by the Hough transform.
-//
-</pre>
-*/
-
-ClassImp(AliHLTTPCHoughEval)
-
-AliHLTTPCHoughEval::AliHLTTPCHoughEval()
-{
- //default ctor
- fRemoveFoundTracks = kFALSE;
- fNumOfPadsToLook = 1;
- fNumOfRowsToMiss = 1;
- fEtaHistos=0;
- fRowPointers = 0;
-}
-
-
-AliHLTTPCHoughEval::~AliHLTTPCHoughEval()
-{
- //dtor
- fHoughTransformer = 0;
- if(fRowPointers)
- {
- for(Int_t i=0; i<fNrows; i++)
- fRowPointers[i] = 0;
- delete [] fRowPointers;
- }
-}
-
-void AliHLTTPCHoughEval::InitTransformer(AliHLTTPCHoughTransformer *transformer)
-{
- //Init hough transformer
- fHoughTransformer = transformer;
- fSlice = fHoughTransformer->GetSlice();
- fPatch = fHoughTransformer->GetPatch();
- fNrows = AliHLTTPCTransform::GetLastRow(fPatch) - AliHLTTPCTransform::GetFirstRow(fPatch) + 1;
- fNEtaSegments = fHoughTransformer->GetNEtaSegments();
- fEtaMin = fHoughTransformer->GetEtaMin();
- fEtaMax = fHoughTransformer->GetEtaMax();
- fZVertex = fHoughTransformer->GetZVertex();
- GenerateLUT();
-}
-
-void AliHLTTPCHoughEval::GenerateLUT()
-{
- //Generate a Look-up table, to limit the access to raw data
-
- if(!fRowPointers)
- fRowPointers = new AliHLTTPCDigitRowData*[fNrows];
-
- AliHLTTPCDigitRowData *tempPt = (AliHLTTPCDigitRowData*)fHoughTransformer->GetDataPointer();
- if(!tempPt)
- printf("\nAliHLTTPCHoughEval::GenerateLUT : Zero data pointer\n");
-
- for(Int_t i=AliHLTTPCTransform::GetFirstRow(fPatch); i<=AliHLTTPCTransform::GetLastRow(fPatch); i++)
- {
- Int_t prow = i - AliHLTTPCTransform::GetFirstRow(fPatch);
- fRowPointers[prow] = tempPt;
- AliHLTTPCMemHandler::UpdateRowPointer(tempPt);
- }
-
-}
-
-Bool_t AliHLTTPCHoughEval::LookInsideRoad(AliHLTTPCHoughTrack *track,Int_t &nrowscrossed,Int_t *rowrange,Bool_t remove)
-{
- //Look at rawdata along the road specified by the track candidates.
- //If track is good, return true, if not return false.
-
- Int_t sector,row;
-
- Int_t nrow=0,npixs=0;//,rows_crossed=0;
- Float_t xyz[3];
-
- Int_t totalcharge=0;//total charge along the road
-
- //for(Int_t padrow = AliHLTTPCTransform::GetFirstRow(fPatch); padrow <= AliHLTTPCTransform::GetLastRow(fPatch); padrow++)
- for(Int_t padrow = rowrange[0]; padrow<=rowrange[1]; padrow++)
- {
- Int_t prow = padrow - AliHLTTPCTransform::GetFirstRow(fPatch);
- if(track->IsHelix())
- {
- if(!track->GetCrossingPoint(padrow,xyz))
- {
- continue;
- }
- }
- else
- {
- track->GetLineCrossingPoint(padrow,xyz);
- xyz[0] += AliHLTTPCTransform::Row2X(track->GetFirstRow());
- Float_t r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
- xyz[2] = r*track->GetTgl();
- }
-
- AliHLTTPCTransform::Slice2Sector(fSlice,padrow,sector,row);
- AliHLTTPCTransform::Local2Raw(xyz,sector,row);
-
- npixs=0;
-
- //Get the timebins for this pad
- AliHLTTPCDigitRowData *tempPt = fRowPointers[prow];
- if(!tempPt)
- {
- printf("AliHLTTPCHoughEval::LookInsideRoad : Zero data pointer\n");
- continue;
- }
-
- //Look at both sides of the pad:
- for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
- {
- AliHLTTPCDigitData *digPt = tempPt->fDigitData;
- for(UInt_t j=0; j<tempPt->fNDigit; j++)
- {
- Int_t pad = digPt[j].fPad;
- Int_t charge = digPt[j].fCharge;
- if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
- if(pad < p) continue;
- if(pad > p) break;
- UShort_t time = digPt[j].fTime;
- Double_t eta = AliHLTTPCTransform::GetEta(fSlice,padrow,pad,time);
- Int_t pixelindex = fHoughTransformer->GetEtaIndex(eta);
- if(pixelindex != track->GetEtaIndex()) continue;
- totalcharge += digPt[j].fCharge;
- if(remove)
- digPt[j].fCharge = 0; //Erease the track from image
- npixs++;
- }
- }
-
- if(npixs > 1)//At least 2 digits on this padrow
- {
- nrow++;
- }
- }
- if(remove)
- return kTRUE;
-
- nrowscrossed += nrow; //Update the number of rows crossed.
-
- if(nrow >= rowrange[1]-rowrange[0]+1 - fNumOfRowsToMiss)//this was a good track
- {
- if(fRemoveFoundTracks)
- {
- Int_t dummy=0;
- LookInsideRoad(track,dummy,rowrange,kTRUE);
- }
- return kTRUE;
- }
- else
- return kFALSE;
-}
-
-void AliHLTTPCHoughEval::FindEta(AliHLTTPCTrackArray *tracks)
-{
- //Find the corresponding eta slice hough space
- Int_t sector,row;
- Float_t xyz[3];
-
- Int_t ntracks = tracks->GetNTracks();
- fEtaHistos = new AliHLTTPCHistogram1D*[ntracks];
-
- Char_t hname[100];
- for(Int_t i=0; i<ntracks; i++)
- {
- sprintf(hname,"etahist_%d",i);
- fEtaHistos[i] = new AliHLTTPCHistogram1D(hname,hname,100,0,1);
- }
- Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
-
- for(Int_t ntr=0; ntr<ntracks; ntr++)
- {
- AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)tracks->GetCheckedTrack(ntr);
- if(!track) continue;
- for(Int_t padrow = AliHLTTPCTransform::GetFirstRow(fPatch); padrow <= AliHLTTPCTransform::GetLastRow(fPatch); padrow++)
- {
- Int_t prow = padrow - AliHLTTPCTransform::GetFirstRow(fPatch);
-
- if(!track->GetCrossingPoint(padrow,xyz))
- {
- printf("AliHLTTPCHoughEval::LookInsideRoad : Track does not cross line!!\n");
- continue;
- }
-
- AliHLTTPCTransform::Slice2Sector(fSlice,padrow,sector,row);
- AliHLTTPCTransform::Local2Raw(xyz,sector,row);
-
- //Get the timebins for this pad
- AliHLTTPCDigitRowData *tempPt = fRowPointers[prow];
- if(!tempPt)
- {
- printf("AliHLTTPCHoughEval::LookInsideRoad : Zero data pointer\n");
- continue;
- }
-
- //Look at both sides of the pad:
- for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
- {
- AliHLTTPCDigitData *digPt = tempPt->fDigitData;
- for(UInt_t j=0; j<tempPt->fNDigit; j++)
- {
- UChar_t pad = digPt[j].fPad;
- Int_t charge = digPt[j].fCharge;
- if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
- if(pad < p) continue;
- if(pad > p) break;
- UShort_t time = digPt[j].fTime;
- Double_t eta = AliHLTTPCTransform::GetEta(fSlice,padrow,pad,time);
- Int_t pixelindex = (Int_t)(eta/etaslice);
- if(pixelindex > track->GetEtaIndex()+1) continue;
- if(pixelindex < track->GetEtaIndex()-1) break;
- fEtaHistos[ntr]->Fill(eta,digPt[j].fCharge);
- }
- }
- }
- }
-
- for(Int_t i=0; i<ntracks; i++)
- {
- AliHLTTPCHistogram1D *hist = fEtaHistos[i];
- Int_t maxbin = hist->GetMaximumBin();
- Double_t maxvalue = hist->GetBinContent(maxbin);
- AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)tracks->GetCheckedTrack(i);
- if(!track) continue;
- if(hist->GetBinContent(maxbin-1)<maxvalue && hist->GetBinContent(maxbin+1)<maxvalue)
- {
- track->SetWeight((Int_t)maxvalue,kTRUE);
- track->SetEta(hist->GetBinCenter(maxbin));
- track->SetNHits(track->GetWeight());
- }
- else
- {
- track->SetWeight(0);
- tracks->Remove(i); //remove this track, because it was not a peak
- }
- }
- tracks->Compress();
-
- //for(Int_t i=0; i<ntracks; i++)
- //delete fEtaHistos[i];
- //delete []Â fEtaHistos;
-}
-
-void AliHLTTPCHoughEval::DisplayEtaSlice(Int_t etaindex,AliHLTTPCHistogram *hist)
-{
- //Display the current raw data inside the (slice,patch)
-
- if(!hist)
- {
- printf("AliHLTTPCHoughEval::DisplayEtaSlice : No input histogram!\n");
- return;
- }
-
- for(Int_t padrow = AliHLTTPCTransform::GetFirstRow(fPatch); padrow <= AliHLTTPCTransform::GetLastRow(fPatch); padrow++)
- {
- Int_t prow = padrow - AliHLTTPCTransform::GetFirstRow(fPatch);
-
- AliHLTTPCDigitRowData *tempPt = fRowPointers[prow];
- if(!tempPt)
- {
- printf("AliHLTTPCHoughEval::DisplayEtaSlice : Zero data pointer\n");
- continue;
- }
-
- AliHLTTPCDigitData *digPt = tempPt->fDigitData;
- if((Int_t)tempPt->fRow != padrow)
- {
- printf("\nAliHLTTPCHoughEval::DisplayEtaSlice : Mismatching padrows!!!\n");
- return;
- }
- for(UInt_t j=0; j<tempPt->fNDigit; j++)
- {
- UChar_t pad = digPt[j].fPad;
- UChar_t charge = digPt[j].fCharge;
- UShort_t time = digPt[j].fTime;
- if((Int_t)charge <= fHoughTransformer->GetLowerThreshold() || (Int_t)charge >= fHoughTransformer->GetUpperThreshold()) continue;
- Float_t xyz[3];
- Int_t sector,row;
- AliHLTTPCTransform::Slice2Sector(fSlice,padrow,sector,row);
- AliHLTTPCTransform::Raw2Local(xyz,sector,row,pad,time);
- xyz[2] -= fZVertex;
- Double_t eta = AliHLTTPCTransform::GetEta(xyz);
- Int_t pixelindex = fHoughTransformer->GetEtaIndex(eta);//(Int_t)(eta/etaslice);
- if(pixelindex != etaindex) continue;
- hist->Fill(xyz[0],xyz[1],charge);
- }
- }
-
-}
-
-void AliHLTTPCHoughEval::CompareMC(AliHLTTPCTrackArray */*tracks*/,Char_t */*trackfile*/,Int_t /*threshold*/)
-{
- /*
- struct GoodTrack goodtracks[15000];
- Int_t nt=0;
- ifstream in(trackfile);
- if(in)
- {
- printf("Reading good tracks from file %s\n",trackfile);
- while (in>>goodtracks[nt].label>>goodtracks[nt].code>>
- goodtracks[nt].px>>goodtracks[nt].py>>goodtracks[nt].pz>>
- goodtracks[nt].pt>>goodtracks[nt].eta>>goodtracks[nt].nhits)
- {
- nt++;
- if (nt==15000)
- {
- cerr<<"Too many good tracks"<<endl;
- break;
- }
- }
- if (!in.eof())
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughEval::CompareMC","Input file")
- <<"Error in file reading"<<ENDLOG;
- return;
- }
- }
- else
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughEval::CompareMC","Input")
- <<"No input trackfile "<<trackfile<<ENDLOG;
- }
-
- Int_t *particles = new Int_t[fNEtaSegments];
- Int_t *ftracks = new Int_t[fNEtaSegments];
- for(Int_t i=0; i<fNEtaSegments; i++)
- {
- particles[i]=0;
- ftracks[i]=0;
- }
-
- TH1F *ptgood = new TH1F("ptgood","ptgood",5,0,2);
- TH1F *ptfound = new TH1F("ptfound","ptgood",5,0,2);
- TH1F *pteff = new TH1F("pteff","pteff",5,0,2);
- TH1F *etafound = new TH1F("etafound","etafound",5,0,1);
- TH1F *etagood = new TH1F("etagood","etagood",5,0,1);
- TH1F *etaeff = new TH1F("etaeff","etaeff",5,0,1);
-
- Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
- for(Int_t i=0; i<tracks->GetNTracks(); i++)
- {
- AliHLTTPCHoughTrack *tr = (AliHLTTPCHoughTrack*)tracks->GetCheckedTrack(i);
- if(!tr) continue;
- if(tr->GetWeight()<threshold) continue;
- Int_t trackindex = tr->GetEtaIndex();
- if(trackindex <0 || trackindex >= fNEtaSegments) continue;
- ftracks[trackindex]++;
- ptfound->Fill(tr->GetPt());
- etafound->Fill(tr->GetEta());
- }
- for(Int_t i=0; i<nt; i++)
- {
- if(goodtracks[i].nhits < 174) continue;
- if(goodtracks[i].pt < 0.2) continue;
- Int_t particleindex = (Int_t)(goodtracks[i].eta/etaslice);
- if(particleindex < 0 || particleindex >= fNEtaSegments) continue;
- particles[particleindex]++;
- ptgood->Fill(goodtracks[i].pt);
- etagood->Fill(goodtracks[i].eta);
- }
-
- Double_t found=0;
- Double_t good =0;
- for(Int_t i=0; i<fNEtaSegments; i++)
- {
- //printf("Slice %d : Found tracks %d, good tracks %d\n",i,ftracks[i],particles[i]);
- found += ftracks[i];
- good += particles[i];
- }
- printf("And the total efficiency was: %f\n",found/good);
-
- ptgood->Sumw2(); ptfound->Sumw2();
- etagood->Sumw2(); etafound->Sumw2();
- pteff->Divide(ptfound,ptgood,1,1,"b");
- etaeff->Divide(etafound,etagood,1,1,"b");
- TFile *file = TFile::Open("eff.root","RECREATE");
- ptgood->Write();
- ptfound->Write();
- pteff->Write();
- etafound->Write();
- etagood->Write();
- etaeff->Write();
- file->Close();
-
- delete [] particles;
- delete [] ftracks;
- */
-}
+++ /dev/null
-// @(#) $Id$
-// origin: hough/AliL3HoughEval.h,v 1.17 Thu Jun 17 10:36:14 2004 UTC by cvetan
-
-#ifndef ALIHLTTPCHOUGHEVAL_H
-#define ALIHLTTPCHOUGHEVAL_H
-
-#include "AliHLTTPCRootTypes.h"
-
-class AliHLTTPCTrackArray;
-class AliHLTTPCHoughTransformer;
-class AliHLTTPCHoughTrack;
-class AliHLTTPCDigitRowData;
-class AliHLTTPCHistogram;
-class AliHLTTPCHistogram1D;
-
-class AliHLTTPCHoughEval {
-
- public:
- AliHLTTPCHoughEval();
- virtual ~AliHLTTPCHoughEval();
-
- void InitTransformer(AliHLTTPCHoughTransformer *transformer);
- void GenerateLUT();
- void DisplayEtaSlice(Int_t etaindex,AliHLTTPCHistogram *hist);
- Bool_t LookInsideRoad(AliHLTTPCHoughTrack *track,Int_t &nrowscrossed,Int_t *rowrange,Bool_t remove=kFALSE);
- void CompareMC(AliHLTTPCTrackArray *tracks,Char_t *goodtracks="good_tracks",Int_t treshold=0);
- void FindEta(AliHLTTPCTrackArray *tracks);
-
- //Getters
- AliHLTTPCHistogram1D *GetEtaHisto(Int_t i) {if(!fEtaHistos) return 0; if(!fEtaHistos[i]) return 0; return fEtaHistos[i];}
-
- //Setters:
- void RemoveFoundTracks() {fRemoveFoundTracks = kTRUE;}
- void SetNumOfRowsToMiss(Int_t i) {fNumOfRowsToMiss = i;}
- void SetNumOfPadsToLook(Int_t i) {fNumOfPadsToLook = i;}
- void SetSlice(Int_t i) {fSlice=i;}
- void SetZVertex(Float_t zvertex) {fZVertex=zvertex;}
-
- private:
-
- Int_t fSlice;//Index of the slice being processed
- Int_t fPatch;//Index of the patch being processed
- Int_t fNrows;//Number of padrows inside the patch
- Int_t fNEtaSegments;//Number of eta slices
- Double_t fEtaMin;//Minimum allowed eta
- Double_t fEtaMax;//Maximum allowed eta
- Int_t fNumOfPadsToLook;//Padrow search window
- Int_t fNumOfRowsToMiss;//Maximum numbers of padrow which could be missed
- AliHLTTPCHistogram1D **fEtaHistos; //!
- Float_t fZVertex;//Z position of the primary vertex
-
- //Flags
- Bool_t fRemoveFoundTracks;//Remove the found tracks or not?
-
- AliHLTTPCHoughTransformer *fHoughTransformer; //!
- AliHLTTPCDigitRowData **fRowPointers; //!
-
- ClassDef(AliHLTTPCHoughEval,1) //Hough transform verfication class
-
-};
-
-#endif
+++ /dev/null
-// @(#) $Id$
-// origin: hough/AliL3HoughKalmanTrack.cxx,v 1.3 Tue Sep 5 08:45:27 2006 UTC by hristov
-
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-//-------------------------------------------------------------------------
-// Implementation of the HLT TPC Hough Kalman track class
-//
-// Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
-//-------------------------------------------------------------------------
-
-#include "AliHLTTPCHoughKalmanTrack.h"
-
-#include "AliHLTStdIncludes.h"
-#include "AliHLTTPCHoughTrack.h"
-#include "AliHLTTPCHoughTransformer.h"
-#include "AliHLTTPCHoughTransformerRow.h"
-#include "AliHLTTPCHistogram.h"
-
-Int_t CalcExternalParams(const AliHLTTPCHoughTrack& t, Double_t deltax, Double_t deltay, Double_t deltaeta, const Double_t zvertex, const Double_t xhit, Double_t xx[5]);
-
-ClassImp(AliHLTTPCHoughKalmanTrack)
-
-//____________________________________________________________________________
-AliHLTTPCHoughKalmanTrack::AliHLTTPCHoughKalmanTrack(const AliHLTTPCHoughTrack& t) throw (const Char_t *)
- : AliTPCtrack()
-{
- // The method constructs an AliHLTTPCHoughKalmanTrack object
- // from an HLT TPC Hough track
-
- SetChi2(0.);
- SetNumberOfClusters(t.GetLastRow()-t.GetFirstRow());
- SetLabel(t.GetMCid());
- SetFakeRatio(0.);
- SetMass(0.13957);
-
- fdEdx=0;
- Double_t alpha = fmod((t.GetSector()+0.5)*(2*TMath::Pi()/18),2*TMath::Pi());
- if (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
- else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
-
- const Double_t xhit = 82.97;
- const Double_t zvertex = t.GetFirstPointZ();
- Double_t par[5];
- Double_t deltax = t.GetPterr();
- Double_t deltay = t.GetPsierr();
- Double_t deltaeta = t.GetTglerr();
- if(CalcExternalParams(t,0,0,0,zvertex,xhit,par)==0) throw "AliHLTTPCHoughKalmanTrack: conversion failed !\n";
-
- Double_t cnv=1./(GetBz()*kB2C);
-
- //and covariance matrix
- //For the moment estimate the covariance matrix numerically
- Double_t xx1[5];
- if(CalcExternalParams(t,deltax,0,0,zvertex,xhit,xx1)==0) throw "AliHLTTPCHoughKalmanTrack: conversion failed !\n";
- Double_t xx2[5];
- if(CalcExternalParams(t,0,deltay,0,zvertex,xhit,xx2)==0) throw "AliHLTTPCHoughKalmanTrack: conversion failed !\n";
- Double_t xx3[5];
- if(CalcExternalParams(t,0,0,deltaeta,zvertex,xhit,xx3)==0) throw "AliHLTTPCHoughKalmanTrack: conversion failed !\n";
-
- Double_t dx1[5],dx2[5],dx3[5];
- for(Int_t i=0;i<5;i++) {
- dx1[i]=xx1[i]-par[i];
- dx2[i]=xx2[i]-par[i];
- dx3[i]=xx3[i]-par[i];
- }
-
- Double_t cov[15]={
- dx1[0]*dx1[0]+dx2[0]*dx2[0],
- 0., dx3[1]*dx3[1],
- 0., 0., dx1[2]*dx1[2]+dx2[2]*dx2[2],
- 0., dx3[3]*dx3[1], 0., dx3[3]*dx3[3],
- 0., 0., 0., 0., dx1[4]*dx1[4]+dx2[4]*dx2[4]
- };
- /*
- fC20=dx1[2]*dx1[0]+dx2[2]*dx2[0];
- fC40=dx1[4]*dx1[0]+dx2[4]*dx2[0];
- fC42=dx1[4]*dx1[2]+dx2[4]*dx2[2];
- fC33=dx3[3]*dx3[3];
- fC11=dx3[1]*dx3[1];
- fC31=dx3[3]*dx3[1];
- fC10=fC30=fC21=fC41=fC32=fC43=0;
- fC20=fC42=fC40=0;
- */
- cov[10]*=cnv; cov[11]*=cnv; cov[12]*=cnv; cov[13]*=cnv; cov[14]*=(cnv*cnv);
- par[4]*=cnv;
-
- Set(xhit,alpha,par,cov);
-
-}
-
-//____________________________________________________________________________
-Int_t CalcExternalParams(const AliHLTTPCHoughTrack& t, Double_t deltax, Double_t deltay, Double_t deltaeta, const Double_t zvertex, const Double_t xhit, Double_t xx[5])
-{
- // Translate the parameters of the Hough tracks into
- // AliKalmanTrack paramters
-
- //First get the emiision angle and track curvature
- Double_t binx = t.GetBinX()+deltax;
- Double_t biny = t.GetBinY()+deltay;
- Double_t psi = atan((binx-biny)/(AliHLTTPCHoughTransformerRow::GetBeta1()-AliHLTTPCHoughTransformerRow::GetBeta2()));
- Double_t kappa = 2.0*(binx*cos(psi)-AliHLTTPCHoughTransformerRow::GetBeta1()*sin(psi));
- Double_t radius = 1./kappa;
-
- //Local y coordinate
- Double_t centerx = -1.*radius*sin(psi);
- Double_t centery = radius*cos(psi);
- Double_t aa = (xhit - centerx)*(xhit - centerx);
- Double_t r2 = radius*radius;
- if(aa > r2) return 0;
- Double_t aa2 = sqrt(r2 - aa);
- Double_t y1 = centery + aa2;
- Double_t y2 = centery - aa2;
- Double_t yhit = y1;
- if(fabs(y2) < fabs(y1)) yhit = y2;
-
- //Local z coordinate
- Double_t stot = sqrt(xhit*xhit+yhit*yhit);
- Double_t zhit;
-
- //Lambda
- Double_t eta=t.GetPseudoRapidity()+deltaeta;
- Double_t theta = 2*atan(exp(-1.*eta));
- Double_t tanl = 1./tan(theta);
- zhit = zvertex + stot*tanl;
-
- xx[0] = yhit;
- xx[1] = zhit;
- xx[2] = (xhit-centerx)/radius;
- xx[3] = tanl;
- xx[4] = kappa;
- return 1;
-}
+++ /dev/null
-// @(#) $Id$
-// origin: hough/AliL3HoughKalmanTrack.h,v 1.1 Thu Mar 31 04:48:58 2005 UTC by cvetan
-
-#ifndef ALIHLTTPCHOUGHKALMANTRACK_H
-#define ALIHLTTPCHOUGHKALMANTRACK_H
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice */
-
-//-------------------------------------------------------------------------
-// High Level Trigger TPC Hough Kalman Track Class
-//
-// Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
-//-------------------------------------------------------------------------
-
-
-/*****************************************************************************
- * October 11, 2004 *
- * The class inherits from the off-line AliTPCtrack class. *
- * It is used to transform AliHLTHoughTrack into AliTPCTrack, which is *
- * then stored as AliESDtrack object in the ESD *
- *****************************************************************************/
-
-#include <AliTPCtrack.h>
-
-class AliHLTTPCHoughTrack;
-class AliHLTTPCHoughTransformer;
-
-class AliHLTTPCHoughKalmanTrack : public AliTPCtrack {
-public:
- AliHLTTPCHoughKalmanTrack(const AliHLTTPCHoughTrack& t) throw (const Char_t *);
-
- ClassDef(AliHLTTPCHoughKalmanTrack,1) //TPC TPC Hough track
-};
-
-#endif
+++ /dev/null
-// @(#) $Id$
-// origin: hough/AliL3HoughMaxFinder.cxx,v 1.13 Tue Mar 28 18:05:12 2006 UTC by alibrary
-
-// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
-//*-- Copyright © ALICE HLT Group
-
-#include <strings.h>
-
-#include <TNtuple.h>
-#include <TFile.h>
-
-#include "AliHLTTPCLogging.h"
-#include "AliHLTTPCHoughMaxFinder.h"
-#include "AliHLTTPCHistogram.h"
-#include "AliHLTTPCTrackArray.h"
-#include "AliHLTTPCHoughTrack.h"
-#include "AliHLTTPCTransform.h"
-#include "AliHLTTPCHoughTransformerRow.h"
-
-#if __GNUC__ >= 3
-using namespace std;
-#endif
-
-/** \class AliHLTTPCHoughMaxFinder
-<pre>
-//_____________________________________________________________
-// AliHLTTPCHoughMaxFinder
-//
-// Maximum finder
-//
-</pre>
-*/
-
-ClassImp(AliHLTTPCHoughMaxFinder)
-
-
-AliHLTTPCHoughMaxFinder::AliHLTTPCHoughMaxFinder()
-{
- //Default constructor
- fThreshold = 0;
- fCurrentEtaSlice = -1;
- fHistoType=0;
- fXPeaks=0;
- fYPeaks=0;
- fWeight=0;
- fNPeaks=0;
- fN1PeaksPrevEtaSlice=0;
- fN2PeaksPrevEtaSlice=0;
- fSTARTXPeaks=0;
- fSTARTYPeaks=0;
- fENDXPeaks=0;
- fENDYPeaks=0;
- fSTARTETAPeaks=0;
- fENDETAPeaks=0;
- fNMax=0;
- fGradX=1;
- fGradY=1;
- fNtuppel = 0;
-}
-
-AliHLTTPCHoughMaxFinder::AliHLTTPCHoughMaxFinder(Char_t *histotype,Int_t nmax,AliHLTTPCHistogram *hist)
-{
- //Constructor
-
- //fTracks = new AliHLTTPCTrackArray("AliHLTTPCHoughTrack");
- if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
- if(strcmp(histotype,"DPsi")==0) fHistoType='l';
-
- fCurrentEtaSlice = -1;
-
- if(hist)
- fCurrentHisto = hist;
-
- fGradX=1;
- fGradY=1;
- fNMax=nmax;
- fXPeaks = new Float_t[fNMax];
- fYPeaks = new Float_t[fNMax];
- fWeight = new Int_t[fNMax];
- fSTARTXPeaks = new Int_t[fNMax];
- fSTARTYPeaks = new Int_t[fNMax];
- fENDXPeaks = new Int_t[fNMax];
- fENDYPeaks = new Int_t[fNMax];
- fSTARTETAPeaks = new Int_t[fNMax];
- fENDETAPeaks = new Int_t[fNMax];
- fNtuppel = 0;
- fThreshold=0;
-}
-
-AliHLTTPCHoughMaxFinder::~AliHLTTPCHoughMaxFinder()
-{
- //Destructor
- if(fXPeaks)
- delete [] fXPeaks;
- if(fYPeaks)
- delete [] fYPeaks;
- if(fWeight)
- delete [] fWeight;
- if(fSTARTXPeaks)
- delete [] fSTARTXPeaks;
- if(fSTARTYPeaks)
- delete [] fSTARTYPeaks;
- if(fENDXPeaks)
- delete [] fENDXPeaks;
- if(fENDYPeaks)
- delete [] fENDYPeaks;
- if(fSTARTETAPeaks)
- delete [] fSTARTETAPeaks;
- if(fENDETAPeaks)
- delete [] fENDETAPeaks;
- if(fNtuppel)
- delete fNtuppel;
-}
-
-void AliHLTTPCHoughMaxFinder::Reset()
-{
- // Method to reinit the Peak Finder
- for(Int_t i=0; i<fNMax; i++)
- {
- fXPeaks[i]=0;
- fYPeaks[i]=0;
- fWeight[i]=0;
- fSTARTXPeaks[i]=0;
- fSTARTYPeaks[i]=0;
- fENDXPeaks[i]=0;
- fENDYPeaks[i]=0;
- fSTARTETAPeaks[i]=0;
- fENDETAPeaks[i]=0;
- }
- fNPeaks=0;
- fN1PeaksPrevEtaSlice=0;
- fN2PeaksPrevEtaSlice=0;
-}
-
-void AliHLTTPCHoughMaxFinder::CreateNtuppel()
-{
- // Fill a NTuple with the peak parameters
- //content#; neighbouring bins of the peak.
- fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
- fNtuppel->SetDirectory(0);
-}
-
-void AliHLTTPCHoughMaxFinder::WriteNtuppel(Char_t *filename)
-{
- // Write the NTuple with the peak parameters
- TFile *file = TFile::Open(filename,"RECREATE");
- if(!file)
- {
- cerr<<"AliHLTTPCHoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl;
- return;
- }
- fNtuppel->Write();
- file->Close();
-}
-
-void AliHLTTPCHoughMaxFinder::FindAbsMaxima()
-{
- // Simple Peak Finder in the Hough space
- if(!fCurrentHisto)
- {
- cerr<<"AliHLTTPCHoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
- return;
- }
- AliHLTTPCHistogram *hist = fCurrentHisto;
-
- if(hist->GetNEntries() == 0)
- return;
-
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
- Int_t bin;
- Double_t value,maxvalue=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>maxvalue)
- {
- maxvalue = value;
- maxxbin = xbin;
- maxybin = ybin;
- }
- }
- }
-
- if(maxvalue == 0)
- return;
-
- if(fNPeaks > fNMax)
- {
- cerr<<"AliHLTTPCHoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl;
- return;
- }
-
- Double_t maxx = hist->GetBinCenterX(maxxbin);
- Double_t maxy = hist->GetBinCenterY(maxybin);
- fXPeaks[fNPeaks] = maxx;
- fYPeaks[fNPeaks] = maxy;
- fWeight[fNPeaks] = (Int_t)maxvalue;
-
- fNPeaks++;
- if(fNtuppel)
- {
- 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(maxx,maxy,maxvalue,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
- }
-}
-
-void AliHLTTPCHoughMaxFinder::FindBigMaxima()
-{
- // Another Peak finder
- AliHLTTPCHistogram *hist = fCurrentHisto;
-
- if(hist->GetNEntries() == 0)
- return;
-
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
- 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++)
- {
- binindex=0;
- for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
- {
- for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
- {
- 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==binindex) break;
- b++;
- //printf("b %d\n",b);
- }
- if(b == binindex)
- {
- //Found maxima
- if(fNPeaks > fNMax)
- {
- cerr<<"AliHLTTPCHoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
- return;
- }
-
- Double_t maxx = hist->GetBinCenterX(xbin);
- Double_t maxy = hist->GetBinCenterY(ybin);
- fXPeaks[fNPeaks] = maxx;
- fYPeaks[fNPeaks] = maxy;
- fNPeaks++;
- }
- }
- }
-}
-
-void AliHLTTPCHoughMaxFinder::FindMaxima(Int_t threshold)
-{
- //Locate all the maxima in input histogram.
- //Maxima is defined as bins with more entries than the
- //immediately neighbouring bins.
-
- if(fCurrentHisto->GetNEntries() == 0)
- return;
-
- Int_t xmin = fCurrentHisto->GetFirstXbin();
- Int_t xmax = fCurrentHisto->GetLastXbin();
- Int_t ymin = fCurrentHisto->GetFirstYbin();
- Int_t ymax = fCurrentHisto->GetLastYbin();
- Int_t bin[9];
- Double_t value[9];
-
- //Float_t max_kappa = 0.001;
- //Float_t max_phi0 = 0.08;
-
- for(Int_t xbin=xmin+1; xbin<=xmax-1; xbin++)
- {
- for(Int_t ybin=ymin+1; ybin<=ymax-1; ybin++)
- {
- bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1);
- bin[1] = fCurrentHisto->GetBin(xbin,ybin-1);
- bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1);
- bin[3] = fCurrentHisto->GetBin(xbin-1,ybin);
- bin[4] = fCurrentHisto->GetBin(xbin,ybin);
- bin[5] = fCurrentHisto->GetBin(xbin+1,ybin);
- bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1);
- bin[7] = fCurrentHisto->GetBin(xbin,ybin+1);
- bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1);
- value[0] = fCurrentHisto->GetBinContent(bin[0]);
- value[1] = fCurrentHisto->GetBinContent(bin[1]);
- value[2] = fCurrentHisto->GetBinContent(bin[2]);
- value[3] = fCurrentHisto->GetBinContent(bin[3]);
- value[4] = fCurrentHisto->GetBinContent(bin[4]);
- value[5] = fCurrentHisto->GetBinContent(bin[5]);
- value[6] = fCurrentHisto->GetBinContent(bin[6]);
- value[7] = fCurrentHisto->GetBinContent(bin[7]);
- value[8] = fCurrentHisto->GetBinContent(bin[8]);
-
-
- if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2]
- && value[4]>value[3] && value[4]>value[5] && value[4]>value[6]
- && value[4]>value[7] && value[4]>value[8])
- {
- //Found a local maxima
- 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)
- {
- cout<<"AliHLTTPCHoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
- return;
- }
-
- //Check the gradient:
- if(value[3]/value[4] > fGradX && value[5]/value[4] > fGradX)
- continue;
-
- if(value[1]/value[4] > fGradY && value[7]/value[4] > fGradY)
- continue;
-
- fXPeaks[fNPeaks] = maxx;
- fYPeaks[fNPeaks] = maxy;
- fWeight[fNPeaks] = (Int_t)value[4];
- fNPeaks++;
-
- /*
- //Check if the peak is overlapping with a previous:
- Bool_t bigger = kFALSE;
- for(Int_t p=0; p<fNPeaks; p++)
- {
- 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] = maxx;
- fYPeaks[p] = maxy;
- fWeight[p] = (Int_t)value[4];
- }
- else
- continue; //previous peak is bigger.
- }
- }
- if(!bigger) //there were no overlapping peaks.
- {
- fXPeaks[fNPeaks] = maxx;
- fYPeaks[fNPeaks] = maxy;
- fWeight[fNPeaks] = (Int_t)value[4];
- fNPeaks++;
- }
- */
- }
- }
- }
-
-}
-
-struct AliHLTTPCWindow
-{
- Int_t fStart; // Start
- Int_t fSum; // Sum
-};
-
-void AliHLTTPCHoughMaxFinder::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
- //in windows (size controlled by kappawindow) summing over each psi-bin.
- //These pre-preaks are then matched between neighbouring kappa-bins to
- //look for real 2D peaks exhbiting the typical cross-shape in the Hough circle transform.
- //The maximum bin within this region is marked as the peak itself, and
- //a few checks is performed to avoid the clear fake peaks (asymmetry check etc.)
-
-
- AliHLTTPCHistogram *hist = fCurrentHisto;
-
- if(!hist)
- {
- cerr<<"AliHLTTPCHoughMaxFinder : No histogram!"<<endl;
- return;
- }
-
- if(hist->GetNEntries() == 0)
- return;
-
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
-
-
- //Start by looking for pre-peaks:
-
- AliHLTTPCWindow **localmaxima = new AliHLTTPCWindow*[hist->GetNbinsY()];
-
- Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
- Int_t n,lastsum,sum;
- Bool_t sumwasrising;
- for(Int_t ybin=ymin; ybin<=ymax; ybin++)
- {
- localmaxima[ybin-ymin] = new AliHLTTPCWindow[hist->GetNbinsX()];
- nmaxs[ybin-ymin] = 0;
- sumwasrising=0;
- lastsum=0;
- n=0;
- for(Int_t xbin=xmin; xbin<=xmax-kappawindow; xbin++)
- {
- sum=0;
- for(Int_t lbin=xbin; lbin<xbin+kappawindow; lbin++)
- sum += hist->GetBinContent(hist->GetBin(lbin,ybin));
-
- if(sum < lastsum)
- {
- if(sum > fThreshold)
- if(sumwasrising)//Previous sum was a local maxima
- {
- localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fStart = xbin-1;
- localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fSum = lastsum;
- nmaxs[ybin-ymin]++;
- }
-
- sumwasrising=0;
- }
- else if(sum > 0)
- sumwasrising=1;
- lastsum=sum;
- }
- }
-
- Int_t match=0;
- Int_t *starts = new Int_t[hist->GetNbinsY()+1];
- Int_t *maxs = new Int_t[hist->GetNbinsY()+1];
-
- for(Int_t ybin=ymax; ybin >= ymin+1; ybin--)
- {
- for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
- {
- 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=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));
- maxybin = ybin;
- maxxbin = k;
- }
-
- //start expanding in the psi-direction:
-
- Int_t lb = localmaxima[ybin-ymin][i].fStart;
- //Int_t ystart=ybin;
- 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 "<<localmaxima[ybin-ymin][i].sum<<endl;
- while(yl >= ymin)
- {
- Bool_t found=0;
- for(Int_t j=0; j<nmaxs[yl-ymin]; j++)
- {
- 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 "<<localmaxima[yl-ymin][j].start<<" sum "<<localmaxima[yl-ymin][j].sum<<endl;
-
- Int_t lmaxvalue=0,lmaxxbin=0;
- 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)
- {
- maxvalue = hist->GetBinContent(hist->GetBin(k,yl));
- maxxbin = k;
- maxybin = yl;
- maxwindow = j;
- }
- if(hist->GetBinContent(hist->GetBin(k,yl)) > lmaxvalue)//local maxima value
- {
- lmaxvalue=hist->GetBinContent(hist->GetBin(k,yl));
- lmaxxbin=k;
- }
- }
- nybins++;
- starts[yl] = localmaxima[yl-ymin][j].fStart;
- maxs[yl] = lmaxxbin;
- localmaxima[yl-ymin][j].fSum=-1; //Mark as used
- found=1;
- 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.
- }
- }
- if(!found || yl == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
- {
- if(nybins > 4)
- {
- //cout<<"ystart "<<ystart<<" and nybins "<<nybins<<endl;
-
- Bool_t truepeak=kTRUE;
-
- //cout<<"Maxima found at xbin "<<maxxbin<<" ybin "<<maxybin<<" value "<<maxvalue<<endl;
- //cout<<"Starting to sum at xbin "<<starts[maxybin-ymin]<<endl;
-
-
- //Look in a window on both sides to probe the asymmetry
- Float_t right=0,left=0;
- for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
- {
- for(Int_t r=maxybin+1; r<=maxybin+3; r++)
- {
- right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
- }
- }
-
- for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
- {
- for(Int_t r=maxybin+1; r<=maxybin+3; r++)
- {
- left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
- }
- }
-
- //cout<<"ratio "<<right/left<<endl;
-
- Float_t upperratio=1,lowerratio=1;
- if(left)
- upperratio = right/left;
-
- right=left=0;
- for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
- {
- for(Int_t r=maxybin-1; r>=maxybin-3; r--)
- {
- right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
- }
- }
-
- for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
- {
- for(Int_t r=maxybin-1; r>=maxybin-3; r--)
- {
- left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
- }
- }
-
- //cout<<"ratio "<<left/right<<endl;
-
- if(right)
- lowerratio = left/right;
-
- if(upperratio > cutratio || lowerratio > cutratio)
- truepeak=kFALSE;
-
- if(truepeak)
- {
-
- fXPeaks[fNPeaks] = hist->GetBinCenterX(maxxbin);
- fYPeaks[fNPeaks] = hist->GetBinCenterY(maxybin);
- fWeight[fNPeaks] = maxvalue;
- fNPeaks++;
-
- /*
- //Calculate the peak using weigthed means:
- Float_t sum=0;
- fYPeaks[fNPeaks]=0;
- for(Int_t k=maxybin-1; k<=maxybin+1; k++)
- {
- Float_t lsum = 0;
- for(Int_t l=starts[k]; l<starts[k]+kappawindow; l++)
- {
- lsum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
- sum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
- }
- fYPeaks[fNPeaks] += lsum*hist->GetBinCenterY(k);
- }
- fYPeaks[fNPeaks] /= sum;
- Int_t ybin1,ybin2;
- if(fYPeaks[fNPeaks] < hist->GetBinCenterY(hist->FindYbin(fYPeaks[fNPeaks])))
- {
- ybin1 = hist->FindYbin(fYPeaks[fNPeaks])-1;
- ybin2 = ybin1+1;
- }
- else
- {
- ybin1 = hist->FindYbin(fYPeaks[fNPeaks]);
- ybin2 = ybin1+1;
- }
-
- Float_t kappa1=0,kappa2=0;
- sum=0;
- for(Int_t k=starts[ybin1]; k<starts[ybin1] + kappawindow; k++)
- {
- kappa1 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin1));
- sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin1));
- }
- kappa1 /= sum;
- sum=0;
- for(Int_t k=starts[ybin2]; k<starts[ybin2] + kappawindow; k++)
- {
- kappa2 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin2));
- sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin2));
- }
- kappa2 /= sum;
-
- fXPeaks[fNPeaks] = ( kappa1*( hist->GetBinCenterY(ybin2) - fYPeaks[fNPeaks] ) +
- kappa2*( fYPeaks[fNPeaks] - hist->GetBinCenterY(ybin1) ) ) /
- (hist->GetBinCenterY(ybin2) - hist->GetBinCenterY(ybin1));
-
- fNPeaks++;
- */
- }
- }
- break;
- }
- else
- yl--;//Search continues...
- }
- }
- }
-
- for(Int_t i=0; i<hist->GetNbinsY(); i++)
- delete localmaxima[i];
-
- delete [] localmaxima;
- delete [] nmaxs;
- delete [] starts;
- delete [] maxs;
-}
-
-struct AliHLTTPCPreYPeak
-{
- 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
-};
-
-void AliHLTTPCHoughMaxFinder::FindAdaptedRowPeaks(Int_t kappawindow,Int_t xsize,Int_t ysize)
-{
- // Peak finder which is working over the Hough Space provided by the AliHLTTPCHoughTransformerRow class
- AliHLTTPCHistogram *hist = fCurrentHisto;
-
- if(!hist)
- {
- cerr<<"AliHLTTPCHoughMaxFinder : No histogram!"<<endl;
- return;
- }
-
- if(hist->GetNEntries() == 0)
- return;
-
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
- Int_t nxbins = hist->GetNbinsX()+2;
- Int_t *content = hist->GetContentArray();
-
- //Start by looking for pre-peaks:
-
- AliHLTTPCPreYPeak **localmaxima = new AliHLTTPCPreYPeak*[hist->GetNbinsY()];
-
- Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
- memset(nmaxs,0,hist->GetNbinsY()*sizeof(Short_t));
- Int_t lastvalue=0,value=0;
- for(Int_t ybin=fNextRow[ymin]; ybin<=ymax; ybin = fNextRow[ybin+1])
- {
- localmaxima[ybin-ymin] = new AliHLTTPCPreYPeak[nxbins-2];
- lastvalue = 0;
- Bool_t found = 0;
- for(Int_t xbin=xmin; xbin<=xmax; xbin++)
- {
- value = content[xbin + nxbins*ybin]; //value = hist->GetBinContent(xbin + nxbins*ybin); //value = hist->GetBinContent(hist->GetBin(xbin,ybin));
- if(value > 0)
- {
- if((value - lastvalue) > 1)
- {
- 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 - lastvalue) <= 1)
- {
- if(found) {
- 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;
- }
- }
- if((value - lastvalue) < -1)
- {
- if(found) {
- localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value;
- nmaxs[ybin-ymin]++;
- found = 0;
- }
- }
- }
- else
- {
- if(found) {
- localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value;
- nmaxs[ybin-ymin]++;
- found = 0;
- }
- }
- lastvalue = value;
-
- }
- if(found) {
- localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = 0;
- nmaxs[ybin-ymin]++;
- }
- }
-
- AliHLTTPCPre2DPeak 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 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 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 localy=ybin-1,nybins=1;
-
- while(localy >= ymin)
- {
- Bool_t found=0;
- for(Int_t j=0; j<nmaxs[localy-ymin]; j++)
- {
- if( (localmaxima[localy-ymin][j].fStartPosition <= (tempxend + kappawindow)) && (localmaxima[localy-ymin][j].fEndPosition >= (tempxstart - kappawindow)))
- {
- if((localmaxima[localy-ymin][j].fMinValue <= localmaxvalue) && (localmaxima[localy-ymin][j].fMaxValue >= localminvalue))
- {
- 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(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 || localy == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
- {
- if((nybins > ysize) && ((localxend-localxstart+1) > xsize) && (localmaxvalue > localprevvalue) && (localmaxvalue > localnextvalue) && (localmaxvalue > localleftvalue) && (localmaxvalue > localrightvalue))
- // if((nybins > ysize) && ((localxend-localxstart+1) > xsize))
- {
- 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)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
- localy--;//Search continues...
- }
- }
- }
-
- //remove fake tracks
-
- for(Int_t i = 0; i < (nmaxima - 1); i++)
- {
- // if(maxima[i].fWeight < 0) continue;
- for(Int_t j = i + 1; j < nmaxima; j++)
- {
- // if(maxima[j].fWeight < 0) continue;
- MergeRowPeaks(&maxima[i],&maxima[j],5.0);
- }
- }
-
- //merge tracks in neighbour eta slices
- Int_t currentnpeaks = fNPeaks;
- for(Int_t i = 0; i < nmaxima; i++) {
- if(maxima[i].fWeight < 0) continue;
- Bool_t merged = kFALSE;
- for(Int_t j = fN1PeaksPrevEtaSlice; j < fN2PeaksPrevEtaSlice; j++) {
- // if(fWeight[j] < 0) continue;
- // Merge only peaks with limited size in eta
- if((fENDETAPeaks[j]-fSTARTETAPeaks[j]) >= 2) continue;
- if((maxima[i].fStartX <= fENDXPeaks[j]+1) && (maxima[i].fEndX >= fSTARTXPeaks[j]-1) && (maxima[i].fStartY <= fENDYPeaks[j]+1) && (maxima[i].fEndY >= fSTARTYPeaks[j]-1)){
- //merge
- merged = kTRUE;
- if(fWeight[j] > 0) {
- 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);
- fSTARTXPeaks[fNPeaks] = maxima[i].fStartX;
- fSTARTYPeaks[fNPeaks] = maxima[i].fStartY;
- fENDXPeaks[fNPeaks] = maxima[i].fEndX;
- fENDYPeaks[fNPeaks] = maxima[i].fEndY;
-
- fWeight[fNPeaks] = abs((Int_t)maxima[i].fWeight) + abs(fWeight[j]);
- fSTARTETAPeaks[fNPeaks] = fSTARTETAPeaks[j];
- fENDETAPeaks[fNPeaks] = fCurrentEtaSlice;
- fNPeaks++;
- }
- fWeight[j] = -abs(fWeight[j]);
- }
- }
- fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].fX);
- fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].fY);
- if(!merged)
- fWeight[fNPeaks] = abs((Int_t)maxima[i].fWeight);
- else
- fWeight[fNPeaks] = -abs((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++;
- }
-
- fN1PeaksPrevEtaSlice = currentnpeaks;
- fN2PeaksPrevEtaSlice = fNPeaks;
-
- for(Int_t ybin=fNextRow[ymin]; ybin<=ymax; ybin = fNextRow[ybin+1])
- delete [] localmaxima[ybin-ymin];
-
- delete [] localmaxima;
- delete [] nmaxs;
-}
-
-void AliHLTTPCHoughMaxFinder::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 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 xbinsides.
-
- if(!fCurrentHisto)
- {
- printf("AliHLTTPCHoughMaxFinder::FindPeak1 : No input histogram\n");
- return;
- }
- if(fCurrentHisto->GetNEntries()==0)
- return;
-
- //Int_t ywindow=2;
- //Int_t xbinsides=1;
-
- //Float_t max_kappa = 0.001;
- //Float_t max_phi0 = 0.08;
-
- Int_t maxsum=0;
-
- Int_t xmin = fCurrentHisto->GetFirstXbin();
- Int_t xmax = fCurrentHisto->GetLastXbin();
- Int_t ymin = fCurrentHisto->GetFirstYbin();
- Int_t ymax = fCurrentHisto->GetLastYbin();
- Int_t nbinsx = fCurrentHisto->GetNbinsX()+1;
-
- AliHLTTPCAxisWindow **windowPt = new AliHLTTPCAxisWindow*[nbinsx];
- AliHLTTPCAxisWindow **anotherPt = new AliHLTTPCAxisWindow*[nbinsx];
-
- for(Int_t i=0; i<nbinsx; i++)
- {
- windowPt[i] = new AliHLTTPCAxisWindow;
-#if defined(__DECCXX)
- bzero((char *)windowPt[i],sizeof(AliHLTTPCAxisWindow));
-#else
- bzero((void*)windowPt[i],sizeof(AliHLTTPCAxisWindow));
-#endif
- anotherPt[i] = windowPt[i];
- }
-
- for(Int_t xbin=xmin; xbin<=xmax; xbin++)
- {
- maxsum = 0;
- for(Int_t ybin=ymin; ybin<=ymax-ywindow; ybin++)
- {
- Int_t suminwindow=0;
- for(Int_t b=ybin; b<ybin+ywindow; b++)
- {
- //inside window
- Int_t bin = fCurrentHisto->GetBin(xbin,b);
- suminwindow += (Int_t)fCurrentHisto->GetBinContent(bin);
- }
-
- if(suminwindow > maxsum)
- {
- maxsum = suminwindow;
- windowPt[xbin]->fYmin = ybin;
- windowPt[xbin]->fYmax = ybin + ywindow;
- windowPt[xbin]->fWeight = suminwindow;
- windowPt[xbin]->fXbin = xbin;
- }
- }
- }
-
- //Sort the windows according to the weight
- SortPeaks(windowPt,0,nbinsx);
-
- Float_t top,butt;
- for(Int_t i=0; i<nbinsx; i++)
- {
- top=butt=0;
- 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]->fWeight > anotherPt[xbin]->fWeight ||
- anotherPt[xbin+1]->fWeight > anotherPt[xbin]->fWeight)
- continue;
-
- 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]->fXbin,j);
- top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
- butt += fCurrentHisto->GetBinContent(bin);
- }
-
- if(butt < fThreshold)
- continue;
-
- 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;
- fNPeaks++;
- if(fNPeaks==fNMax)
- {
- cerr<<"AliHLTTPCHoughMaxFinder::FindPeak1 : Peak array out of range!!!"<<endl;
- break;
- }
- }
-
-
- //Improve the peaks by including the region around in x.
- Float_t ytop,ybutt;
- Int_t prev;
- Int_t w;
- for(Int_t i=0; i<fNPeaks; i++)
- {
- Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]);
- if(xbin - xbinsides < xmin || xbin + xbinsides > xmax) continue;
- top=butt=0;
- ytop=0,ybutt=0;
- w=0;
- prev = xbin - xbinsides+1;
- for(Int_t j=xbin-xbinsides; j<=xbin+xbinsides; j++)
- {
- /*
- //Check if the windows are overlapping:
- if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
- if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
- prev = j;
- */
-
- top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->fWeight;
- butt += anotherPt[j]->fWeight;
-
- 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));
- ybutt += fCurrentHisto->GetBinContent(bin);
- w+=(Int_t)fCurrentHisto->GetBinContent(bin);
- }
- }
-
- fXPeaks[i] = top/butt;
- fYPeaks[i] = ytop/ybutt;
- fWeight[i] = w;
- //cout<<"Setting weight "<<w<<" kappa "<<fXPeaks[i]<<" phi0 "<<fYPeaks[i]<<endl;
-
- /*
- //Check if this peak is overlapping with a previous:
- for(Int_t p=0; p<i; p++)
- {
- //cout<<fabs(fXPeaks[p] - fXPeaks[i])<<" "<<fabs(fYPeaks[p] - fYPeaks[i])<<endl;
- if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa &&
- fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0)
- {
- fWeight[i]=0;
- //break;
- }
- }
- */
- }
-
- for(Int_t i=0; i<nbinsx; i++)
- delete windowPt[i];
- delete [] windowPt;
- delete [] anotherPt;
-}
-
-void AliHLTTPCHoughMaxFinder::SortPeaks(struct AliHLTTPCAxisWindow **a,Int_t first,Int_t last)
-{
- //General sorting routine
- //Sort according to PeakCompare()
-
- static struct AliHLTTPCAxisWindow *tmp;
- static int i; // "static" to save stack space
- int j;
-
- while (last - first > 1) {
- i = first;
- j = last;
- for (;;) {
- while (++i < last && PeakCompare(a[i], a[first]) < 0)
- ;
- while (--j > first && PeakCompare(a[j], a[first]) > 0)
- ;
- if (i >= j)
- break;
-
- tmp = a[i];
- a[i] = a[j];
- a[j] = tmp;
- }
- if (j == first) {
- ++first;
- continue;
- }
- tmp = a[first];
- a[first] = a[j];
- a[j] = tmp;
- if (j - first < last - (j + 1)) {
- SortPeaks(a, first, j);
- first = j + 1; // QSort(j + 1, last);
- } else {
- SortPeaks(a, j + 1, last);
- last = j; // QSort(first, j);
- }
- }
-
-}
-
-Int_t AliHLTTPCHoughMaxFinder::PeakCompare(struct AliHLTTPCAxisWindow *a,struct AliHLTTPCAxisWindow *b) const
-{
- // Peak comparison based on peaks weight
- if(a->fWeight < b->fWeight) return 1;
- if(a->fWeight > b->fWeight) return -1;
- return 0;
-}
-
-void AliHLTTPCHoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
-{
- //Attempt of a more sophisticated peak finder.
- //Finds the best peak in the histogram, and returns the corresponding
- //track object.
-
- if(!fCurrentHisto)
- {
- printf("AliHLTTPCHoughMaxFinder::FindPeak : No histogram!!\n");
- return;
- }
- AliHLTTPCHistogram *hist = fCurrentHisto;
- if(hist->GetNEntries()==0)
- return;
-
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
- Int_t nbinsx = hist->GetNbinsX()+1;
-
- Int_t *m = 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;
- mlow[i]=0;
- mup[i]=0;
- }
-
- Int_t maxx=0,sum=0,maxxbin=0,bin=0;
-
- for(Int_t xbin=xmin; xbin<=xmax; xbin++)
- {
- for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
- {
- sum = 0;
- for(Int_t y=ybin; y <= ybin+t1; y++)
- {
- if(y>ymax) break;
- //Inside window
- bin = hist->GetBin(xbin,y);
- sum += (Int_t)hist->GetBinContent(bin);
-
- }
- if(sum > m[xbin]) //Max value locally in this xbin
- {
- m[xbin]=sum;
- mlow[xbin]=ybin;
- mup[xbin]=ybin + t1;
- }
-
- }
-
- if(m[xbin] > maxx) //Max value globally in x-direction
- {
- maxxbin = xbin;
- maxx = m[xbin];//sum;
- }
- }
- //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 xlow=0,xup=0;
-
- for(Int_t xbin=maxxbin-1; xbin >= xmin; xbin--)
- {
- if(m[xbin] < maxx*t2)
- {
- xlow = xbin+1;
- break;
- }
- }
- for(Int_t xbin = maxxbin+1; xbin <=xmax; xbin++)
- {
- if(m[xbin] < maxx*t2)
- {
- xup = xbin-1;
- break;
- }
- }
-
- Double_t top=0,butt=0,value,xpeak;
- if(xup - xlow + 1 > t3)
- {
- t1 -= 1;
- printf("\nxrange out if limit xup %d xlow %d t1 %d\n\n",xlow,xup,t1);
- if(t1 > 1)
- goto recompute;
- else
- {
- xpeak = hist->GetBinCenterX(maxxbin);
- goto moveon;
- }
- }
-
- //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=xlow; xbin <= xup; xbin++)
- {
- value = hist->GetBinCenterX(xbin);
- top += value*m[xbin];
- butt += m[xbin];
- }
- xpeak = top/butt;
-
- moveon:
-
- //Find the peak in y direction:
- Int_t xl = hist->FindXbin(xpeak);
- if(hist->GetBinCenterX(xl) > xpeak)
- xl--;
-
- Int_t xu = xl + 1;
-
- if(hist->GetBinCenterX(xl) > xpeak || hist->GetBinCenterX(xu) <= xpeak)
- printf("\nAliHLTTPCHoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(xl),xpeak,hist->GetBinCenterX(xu));
-
- //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(xl),hist->GetBinCenterX(xu));
-
- value=top=butt=0;
-
- //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=mlow[xl]; ybin <= mup[xl]; ybin++)
- {
- value = hist->GetBinCenterY(ybin);
- bin = hist->GetBin(xl,ybin);
- top += value*hist->GetBinContent(bin);
- butt += hist->GetBinContent(bin);
- }
- Double_t ypeaklow = top/butt;
-
- //printf("ypeaklow %f\n",ypeaklow);
-
- value=top=butt=0;
- for(Int_t ybin=mlow[xu]; ybin <= mup[xu]; ybin++)
- {
- value = hist->GetBinCenterY(ybin);
- bin = hist->GetBin(xu,ybin);
- top += value*hist->GetBinContent(bin);
- butt += hist->GetBinContent(bin);
- }
- Double_t ypeakup = top/butt;
-
- //printf("ypeakup %f\n",ypeakup);
-
- Double_t xvalueup = hist->GetBinCenterX(xu);
- Double_t xvaluelow = hist->GetBinCenterX(xl);
-
- Double_t ypeak = (ypeaklow*(xvalueup - xpeak) + ypeakup*(xpeak - xvaluelow))/(xvalueup - xvaluelow);
-
-
- //Find the weight:
- //bin = hist->FindBin(xpeak,ypeak);
- //Int_t weight = (Int_t)hist->GetBinContent(bin);
-
- //AliHLTTPCHoughTrack *track = new AliHLTTPCHoughTrack();
- //track->SetTrackParameters(xpeak,ypeak,weight);
- fXPeaks[fNPeaks]=xpeak;
- fYPeaks[fNPeaks]=ypeak;
- fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin);
- fNPeaks++;
-
- delete [] m;
- delete [] mlow;
- delete [] mup;
-
- //return track;
-}
-
-Float_t AliHLTTPCHoughMaxFinder::GetXPeakSize(Int_t i) const
-{
- // Get X size of a peak
- if(i<0 || i>fNMax)
- {
- STDCERR<<"AliHLTTPCHoughMaxFinder::GetXPeakSize : Invalid index "<<i<<STDENDL;
- return 0;
- }
- Float_t binwidth = fCurrentHisto->GetBinWidthX();
- return binwidth*(fENDXPeaks[i]-fSTARTXPeaks[i]+1);
-}
-
-Float_t AliHLTTPCHoughMaxFinder::GetYPeakSize(Int_t i) const
-{
- // Get Y size of a peak
- if(i<0 || i>fNMax)
- {
- STDCERR<<"AliHLTTPCHoughMaxFinder::GetYPeak : Invalid index "<<i<<STDENDL;
- return 0;
- }
- Float_t binwidth = fCurrentHisto->GetBinWidthY();
- return binwidth*(fENDYPeaks[i]-fSTARTYPeaks[i]+1);
-}
-
-Bool_t AliHLTTPCHoughMaxFinder::MergeRowPeaks(AliHLTTPCPre2DPeak *maxima1, AliHLTTPCPre2DPeak *maxima2, Float_t distance)
-{
- // Check the distance between tracks corresponding to given Hough space peaks and if the
- // distance is smaller than some threshold value marks the smaller peak as fake
- AliHLTTPCHistogram *hist = fCurrentHisto;
- Int_t nxbins = hist->GetNbinsX()+2;
-
- Int_t xtrack1=0,xtrack2=0,ytrack1=0,ytrack2=0;
- Int_t deltax = 9999;
- for(Int_t ix1 = maxima1->fStartX; ix1 <= maxima1->fEndX; ix1++) {
- for(Int_t ix2 = maxima2->fStartX; ix2 <= maxima2->fEndX; ix2++) {
- if(abs(ix1 - ix2) < deltax) {
- deltax = abs(ix1 - ix2);
- xtrack1 = ix1;
- xtrack2 = ix2;
- }
- }
- }
- Int_t deltay = 9999;
- for(Int_t iy1 = maxima1->fStartY; iy1 <= maxima1->fEndY; iy1++) {
- for(Int_t iy2 = maxima2->fStartY; iy2 <= maxima2->fEndY; iy2++) {
- if(abs(iy1 - iy2) < deltay) {
- deltay = abs(iy1 - iy2);
- ytrack1 = iy1;
- ytrack2 = iy2;
- }
- }
- }
- Int_t firstrow1 = fTrackFirstRow[xtrack1 + nxbins*ytrack1];
- Int_t lastrow1 = fTrackLastRow[xtrack1 + nxbins*ytrack1];
- Int_t firstrow2 = fTrackFirstRow[xtrack1 + nxbins*ytrack1];
- Int_t lastrow2 = fTrackLastRow[xtrack1 + nxbins*ytrack1];
- Int_t firstrow,lastrow;
- if(firstrow1 < firstrow2)
- firstrow = firstrow2;
- else
- firstrow = firstrow1;
-
- if(lastrow1 > lastrow2)
- lastrow = lastrow2;
- else
- lastrow = lastrow1;
-
- AliHLTTPCHoughTrack track1;
- Float_t x1 = hist->GetPreciseBinCenterX(xtrack1);
- Float_t y1 = hist->GetPreciseBinCenterY(ytrack1);
- Float_t psi1 = atan((x1-y1)/(AliHLTTPCHoughTransformerRow::GetBeta1()-AliHLTTPCHoughTransformerRow::GetBeta2()));
- Float_t kappa1 = 2.0*(x1*cos(psi1)-AliHLTTPCHoughTransformerRow::GetBeta1()*sin(psi1));
- track1.SetTrackParameters(kappa1,psi1,1);
- Float_t firsthit1[3];
- if(!track1.GetCrossingPoint(firstrow,firsthit1)) return kFALSE;
- Float_t lasthit1[3];
- if(!track1.GetCrossingPoint(lastrow,lasthit1)) return kFALSE;
-
- AliHLTTPCHoughTrack track2;
- Float_t x2 = hist->GetPreciseBinCenterX(xtrack2);
- Float_t y2 = hist->GetPreciseBinCenterY(ytrack2);
- Float_t psi2 = atan((x2-y2)/(AliHLTTPCHoughTransformerRow::GetBeta1()-AliHLTTPCHoughTransformerRow::GetBeta2()));
- Float_t kappa2 = 2.0*(x2*cos(psi2)-AliHLTTPCHoughTransformerRow::GetBeta1()*sin(psi2));
- track2.SetTrackParameters(kappa2,psi2,1);
- Float_t firsthit2[3];
- if(!track2.GetCrossingPoint(firstrow,firsthit2)) return kFALSE;
- Float_t lasthit2[3];
- if(!track2.GetCrossingPoint(lastrow,lasthit2)) return kFALSE;
-
- Float_t padpitchlow = AliHLTTPCTransform::GetPadPitchWidth(AliHLTTPCTransform::GetPatch(firstrow));
- Float_t padpitchup = AliHLTTPCTransform::GetPadPitchWidth(AliHLTTPCTransform::GetPatch(lastrow));
- // check the distance between tracks at the edges
- // cout<<"Check "<<firsthit1[1]<<" "<<firsthit2[1]<<" "<<padpitchlow<<" "<<lasthit1[1]<<" "<<lasthit2[1]<<" "<<padpitchup<<" "<<xtrack1<<" "<<ytrack1<<" "<<xtrack2<<" "<<ytrack2<<endl;
- if((fabs(firsthit1[1]-firsthit2[1])/padpitchlow + fabs(lasthit1[1]-lasthit2[1])/padpitchup) < distance) {
- if(maxima1->fSizeX*maxima1->fSizeY > maxima2->fSizeX*maxima2->fSizeY)
- maxima2->fWeight = -fabs(maxima2->fWeight);
- if(maxima1->fSizeX*maxima1->fSizeY < maxima2->fSizeX*maxima2->fSizeY)
- maxima1->fWeight = -fabs(maxima1->fWeight);
- if(maxima1->fSizeX*maxima1->fSizeY == maxima2->fSizeX*maxima2->fSizeY) {
- if(maxima1->fStartX > maxima2->fStartX)
- maxima1->fStartX = maxima2->fStartX;
- if(maxima1->fStartY > maxima2->fStartY)
- maxima1->fStartY = maxima2->fStartY;
- if(maxima1->fEndX < maxima2->fEndX)
- maxima1->fEndX = maxima2->fEndX;
- if(maxima1->fEndY < maxima2->fEndY)
- maxima1->fEndY = maxima2->fEndY;
- maxima1->fX = ((Float_t)maxima1->fStartX + (Float_t)maxima1->fEndX)/2.0;
- maxima1->fY = ((Float_t)maxima1->fStartY + (Float_t)maxima1->fEndY)/2.0;
- maxima1->fSizeX = (maxima1->fEndX - maxima1->fStartX + 1);
- maxima1->fSizeY = (maxima1->fEndY - maxima1->fStartY + 1);
- maxima2->fWeight = -fabs(maxima2->fWeight);
- }
- return kTRUE;
- }
- return kFALSE;
-}
+++ /dev/null
-// @(#) $Id$
-// origin: hough/AliL3HoughMaxFinder.h,v 1.24 Mon Nov 8 11:31:13 2004 UTC by cvetan
-
-#ifndef ALIHLTTPCHOUGHMAXFINDER_H
-#define ALIHLTTPCHOUGHMAXFINDER_H
-
-#include "AliHLTTPCRootTypes.h"
-#include "AliHLTStdIncludes.h"
-
-class AliHLTTPCHistogram;
-class AliHLTTPCTrackArray;
-class AliHLTTPCHoughTrack;
-class TNtuple;
-
-struct AliHLTTPCAxisWindow
-{
- Int_t fYmin; // min Y
- Int_t fYmax; // max Y
- Int_t fXbin; // X bin
- Int_t fWeight; // weight
-};
-
-struct AliHLTTPCPre2DPeak
-{
- 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
-};
-
-class AliHLTTPCHoughMaxFinder {
-
- public:
- AliHLTTPCHoughMaxFinder();
- AliHLTTPCHoughMaxFinder(Char_t *histotype,Int_t nmax,AliHLTTPCHistogram *hist=0);
- virtual ~AliHLTTPCHoughMaxFinder();
- void Reset();
-
- void CreateNtuppel();
- void WriteNtuppel(Char_t *filename);
-
- //Simple maxima finders:
- void FindAbsMaxima();
- void FindBigMaxima();
- void FindMaxima(Int_t threshold=0);
- 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 ywindow=2,Int_t xbinsides=1);
- void SortPeaks(struct AliHLTTPCAxisWindow **a,Int_t first,Int_t last);
- Int_t PeakCompare(struct AliHLTTPCAxisWindow *a,struct AliHLTTPCAxisWindow *b) const;
-
- //Setters:
- void SetGradient(Float_t x,Float_t y) {fGradX=x; fGradY=y;}
- void SetThreshold(Int_t f) {fThreshold = f;}
- void SetHistogram(AliHLTTPCHistogram *hist) {fCurrentHisto = hist;}
- void SetTrackLUTs(UChar_t *tracknrows, UChar_t *trackfirstrow, UChar_t *tracklastrow, UChar_t *nextrow) {fTrackNRows = tracknrows; fTrackFirstRow = trackfirstrow; fTrackLastRow = tracklastrow; fNextRow = nextrow;}
- void SetEtaSlice(Int_t etaslice) {fCurrentEtaSlice = etaslice;}
-
- //Getters:
- 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;}
-
- //Method for merging of peaks produced by AliHLTTPCHoughTransfromerRow
- Bool_t MergeRowPeaks(AliHLTTPCPre2DPeak *maxima1, AliHLTTPCPre2DPeak *maxima2,Float_t distance);
-
- private:
-
- Int_t fThreshold; // Threshold for Peak Finder
- Int_t fCurrentEtaSlice; // Current eta slice being processed
- AliHLTTPCHistogram *fCurrentHisto; //!
-
- UChar_t *fTrackNRows; //!
- UChar_t *fTrackFirstRow; //!
- UChar_t *fTrackLastRow; //!
- UChar_t *fNextRow; //!
-
- 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
-
- TNtuple *fNtuppel; //!
-
- ClassDef(AliHLTTPCHoughMaxFinder,1) //Maximum finder class
-
-};
-
-inline Float_t AliHLTTPCHoughMaxFinder::GetXPeak(Int_t i) const
-{
- if(i<0 || i>fNMax)
- {
- STDCERR<<"AliHLTTPCHoughMaxFinder::GetXPeak : Invalid index "<<i<<STDENDL;
- return 0;
- }
- return fXPeaks[i];
-}
-
-inline Float_t AliHLTTPCHoughMaxFinder::GetYPeak(Int_t i) const
-{
- if(i<0 || i>fNMax)
- {
- STDCERR<<"AliHLTTPCHoughMaxFinder::GetYPeak : Invalid index "<<i<<STDENDL;
- return 0;
- }
- return fYPeaks[i];
-
-}
-
-inline Int_t AliHLTTPCHoughMaxFinder::GetWeight(Int_t i) const
-{
- if(i<0 || i>fNMax)
- {
- STDCERR<<"AliHLTTPCHoughMaxFinder::GetWeight : Invalid index "<<i<<STDENDL;
- return 0;
- }
- return fWeight[i];
-}
-
-inline Int_t AliHLTTPCHoughMaxFinder::GetStartEta(Int_t i) const
-{
- if(i<0 || i>fNMax)
- {
- STDCERR<<"AliHLTTPCHoughMaxFinder::GetStartEta : Invalid index "<<i<<STDENDL;
- return 0;
- }
- return fSTARTETAPeaks[i];
-}
-
-inline Int_t AliHLTTPCHoughMaxFinder::GetEndEta(Int_t i) const
-{
- if(i<0 || i>fNMax)
- {
- STDCERR<<"AliHLTTPCHoughMaxFinder::GetStartEta : Invalid index "<<i<<STDENDL;
- return 0;
- }
- return fENDETAPeaks[i];
-}
-
-#endif
-
+++ /dev/null
-// @(#) $Id$
-// origin: hough/AliL3HoughTrack.cxx,v 1.23 ue Mar 28 18:05:12 2006 UTC by alibrary
-
-// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
-//*-- Copyright © ALICE HLT Group
-
-#include "AliHLTStdIncludes.h"
-
-#include "AliHLTTPCLogging.h"
-#include "AliHLTTPCHoughTrack.h"
-#include "AliHLTTPCTransform.h"
-#include "AliHLTTPCHoughTransformerRow.h"
-
-#if __GNUC__ >= 3
-using namespace std;
-#endif
-
-/** \class AliHLTTPCHoughTrack
-<pre>
-//_____________________________________________________________
-// AliHLTTPCHoughTrack
-//
-// Track class for Hough tracklets
-//
-</pre>
-*/
-
-ClassImp(AliHLTTPCHoughTrack);
-
-
-AliHLTTPCHoughTrack::AliHLTTPCHoughTrack() : AliHLTTPCTrack()
-{
- //Constructor
-
- fWeight = 0;
- fMinDist=0;
- fDLine = 0;
- fPsiLine = 0;
- fIsHelix = true;
- fEtaIndex = -1;
- fEta = 0;
- ComesFromMainVertex(kTRUE);
-}
-
-AliHLTTPCHoughTrack::~AliHLTTPCHoughTrack()
-{
- //dtor
-}
-
-void AliHLTTPCHoughTrack::Set(AliHLTTPCTrack *track)
-{
- //Basically copy constructor
- AliHLTTPCHoughTrack *tpt = (AliHLTTPCHoughTrack*)track;
- SetTrackParameters(tpt->GetKappa(),tpt->GetPsi(),tpt->GetWeight());
- SetEtaIndex(tpt->GetEtaIndex());
- SetEta(tpt->GetEta());
- SetTgl(tpt->GetTgl());
- SetPsi(tpt->GetPsi());
- SetPterr(tpt->GetPterr());
- SetTglerr(tpt->GetTglerr());
- SetPsierr(tpt->GetPsierr());
- SetCenterX(tpt->GetCenterX());
- SetCenterY(tpt->GetCenterY());
- SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
- SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
- SetCharge(tpt->GetCharge());
- SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
- SetSlice(tpt->GetSlice());
- SetHits(tpt->GetNHits(),(UInt_t *)tpt->GetHitNumbers());
- SetMCid(tpt->GetMCid());
- SetBinXY(tpt->GetBinX(),tpt->GetBinY(),tpt->GetSizeX(),tpt->GetSizeY());
- SetSector(tpt->GetSector());
- return;
-
-// fWeight = tpt->GetWeight();
-// fDLine = tpt->GetDLine();
-// fPsiLine = tpt->GetPsiLine();
-// SetNHits(tpt->GetWeight());
-// SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
-// fIsHelix = false;
-}
-
-Int_t AliHLTTPCHoughTrack::Compare(const AliHLTTPCTrack *tpt) const
-{
- //Compare 2 hough tracks according to their weight
- AliHLTTPCHoughTrack *track = (AliHLTTPCHoughTrack*)tpt;
- if(track->GetWeight() < GetWeight()) return 1;
- if(track->GetWeight() > GetWeight()) return -1;
- return 0;
-}
-
-void AliHLTTPCHoughTrack::SetEta(Double_t f)
-{
- //Set eta, and calculate fTanl, which is the tan of dipangle
-
- fEta = f;
- Double_t theta = 2*atan(exp(-1.*fEta));
- Double_t dipangle = AliHLTTPCTransform::PiHalf() - theta;
- Double_t tgl = tan(dipangle);
- SetTgl(tgl);
-}
-
-void AliHLTTPCHoughTrack::UpdateToFirstRow()
-{
- //Update the track parameters to the point where track cross
- //its first padrow.`
-
- //Get the crossing point with the first padrow:
- Float_t xyz[3];
- if(!GetCrossingPoint(GetFirstRow(),xyz))
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTrack::UpdateToFirstRow()","Track parameters")
- <<AliHLTTPCLog::kDec<<"Track does not cross padrow "<<GetFirstRow()<<" centerx "
- <<GetCenterX()<<" centery "<<GetCenterY()<<" Radius "<<GetRadius()<<" tgl "<<GetTgl()<<ENDLOG;
-
- //printf("Track with eta %f tgl %f crosses at x %f y %f z %f on padrow %d\n",GetEta(),GetTgl(),xyz[0],xyz[1],xyz[2],GetFirstRow());
- //printf("Before: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
-
- Double_t radius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
-
- //Get the track parameters
-
- /*
- Double_t x0 = GetR0() * cos(GetPhi0()) ;
- Double_t y0 = GetR0() * sin(GetPhi0()) ;
- */
- Double_t rc = GetRadius();//fabs(GetPt()) / AliHLTTPCTransform::GetBFieldValue();
- Double_t tPhi0 = GetPsi() + GetCharge() * AliHLTTPCTransform::PiHalf() / abs(GetCharge()) ;
- Double_t xc = GetCenterX();//x0 - rc * cos(tPhi0) ;
- Double_t yc = GetCenterY();//y0 - rc * sin(tPhi0) ;
-
- //Check helix and cylinder intersect
- Double_t fac1 = xc*xc + yc*yc ;
- Double_t sfac = sqrt( fac1 ) ;
-
- if ( fabs(sfac-rc) > radius || fabs(sfac+rc) < radius ) {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughTrack::UpdateToFirstRow","Tracks")<<AliHLTTPCLog::kDec<<
- "Track does not intersect"<<ENDLOG;
- return;
- }
-
- //Find intersection
- Double_t fac2 = (radius*radius + fac1 - rc*rc) / (2.00 * radius * sfac ) ;
- Double_t phi = atan2(yc,xc) + GetCharge()*acos(fac2) ;
- Double_t td = atan2(radius*sin(phi) - yc,radius*cos(phi) - xc) ;
-
- //Intersection in z
- if ( td < 0 ) td = td + AliHLTTPCTransform::TwoPi();
- Double_t deltat = fmod((-GetCharge()*td + GetCharge()*tPhi0),AliHLTTPCTransform::TwoPi());
- if ( deltat < 0. ) deltat += AliHLTTPCTransform::TwoPi();
- else if ( deltat > AliHLTTPCTransform::TwoPi() ) deltat -= AliHLTTPCTransform::TwoPi();
- Double_t z = GetZ0() + rc * GetTgl() * deltat ;
-
- Double_t xExtra = radius * cos(phi) ;
- Double_t yExtra = radius * sin(phi) ;
-
- Double_t tPhi = atan2(yExtra-yc,xExtra-xc);
-
- //if ( tPhi < 0 ) tPhi += 2. * M_PI ;
- Double_t tPsi = tPhi - GetCharge() * AliHLTTPCTransform::PiHalf() / abs(GetCharge()) ;
- if ( tPsi > AliHLTTPCTransform::TwoPi() ) tPsi -= AliHLTTPCTransform::TwoPi() ;
- else if ( tPsi < 0. ) tPsi += AliHLTTPCTransform::TwoPi();
-
- //And finally, update the track parameters
- SetR0(radius);
- SetPhi0(phi);
- SetZ0(z);
- SetPsi(tPsi);
- SetFirstPoint(xyz[0],xyz[1],z);
- //printf("After: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
-
- //printf("First point set %f %f %f\n",xyz[0],xyz[1],z);
-
- //Also, set the coordinates of the point where track crosses last padrow:
- GetCrossingPoint(GetLastRow(),xyz);
- SetLastPoint(xyz[0],xyz[1],xyz[2]);
- //printf("last point %f %f %f\n",xyz[0],xyz[1],xyz[2]);
-}
-
-void AliHLTTPCHoughTrack::SetTrackParameters(Double_t kappa,Double_t eangle,Int_t weight)
-{
- //Set track parameters - sort of ctor
- fWeight = weight;
- fMinDist = 100000;
- SetKappa(kappa);
- Double_t pt = fabs(AliHLTTPCTransform::GetBFieldValue()/kappa);
- SetPt(pt);
- Double_t radius = 1/fabs(kappa);
- SetRadius(radius);
- SetFirstPoint(0,0,0);
- SetPsi(eangle); //Psi = emission angle when first point is vertex
- SetPhi0(0); //not defined for vertex reference point
- SetR0(0);
- Double_t charge = -1.*kappa;
- SetCharge((Int_t)copysign(1.,charge));
- Double_t trackPhi0 = GetPsi() + charge*0.5*AliHLTTPCTransform::Pi()/fabs(charge);
- Double_t xc = GetFirstPointX() - GetRadius() * cos(trackPhi0) ;
- Double_t yc = GetFirstPointY() - GetRadius() * sin(trackPhi0) ;
- SetCenterX(xc);
- SetCenterY(yc);
- SetNHits(1); //just for the trackarray IO
- fIsHelix = true;
-}
-
-void AliHLTTPCHoughTrack::SetTrackParametersRow(Double_t alpha1,Double_t alpha2,Double_t eta,Int_t weight)
-{
- //Set track parameters for HoughTransformerRow
- //This includes curvature,emission angle and eta
- Double_t psi = atan((alpha1-alpha2)/(AliHLTTPCHoughTransformerRow::GetBeta1()-AliHLTTPCHoughTransformerRow::GetBeta2()));
- Double_t kappa = 2.0*(alpha1*cos(psi)-AliHLTTPCHoughTransformerRow::GetBeta1()*sin(psi));
- SetTrackParameters(kappa,psi,weight);
-
- Double_t zovr;
- Double_t etaparam1 = AliHLTTPCHoughTransformerRow::GetEtaCalcParam1();
- Double_t etaparam2 = AliHLTTPCHoughTransformerRow::GetEtaCalcParam2();
- if(eta>0)
- zovr = (etaparam1 - sqrt(etaparam1*etaparam1 - 4.*etaparam2*eta))/(2.*etaparam2);
- else
- zovr = -1.*(etaparam1 - sqrt(etaparam1*etaparam1 + 4.*etaparam2*eta))/(2.*etaparam2);
- Double_t r = sqrt(1.+zovr*zovr);
- Double_t exacteta = 0.5*log((1+zovr/r)/(1-zovr/r));
- SetEta(exacteta);
-}
-
-void AliHLTTPCHoughTrack::SetLineParameters(Double_t psi,Double_t D,Int_t weight,Int_t *rowrange,Int_t /*ref_row*/)
-{
- //Initialize a track piece, not yet a track
- //Used in case of straight line transformation
-
- //Transform line parameters to coordinate system of slice:
-
- // D = D + fTransform->Row2X(ref_row)*cos(psi);
-
- fDLine = D;
- fPsiLine = psi;
- fWeight = weight;
- SetNHits(1);
- SetRowRange(rowrange[0],rowrange[1]);
- fIsHelix = false;
-}
-
-void AliHLTTPCHoughTrack::SetBestMCid(Int_t mcid,Double_t mindist)
-{
- //Finds and set the closest mc label
- if(mindist < fMinDist)
- {
- fMinDist = mindist;
- SetMCid(mcid);
- }
-}
-
-void AliHLTTPCHoughTrack::GetLineCrossingPoint(Int_t padrow,Float_t *xy)
-{
- //Returns the crossing point of the track with a given padrow
- if(fIsHelix)
- {
- printf("AliHLTTPCHoughTrack::GetLineCrossingPoint : Track is not a line\n");
- return;
- }
-
- Float_t xhit = AliHLTTPCTransform::Row2X(padrow) - AliHLTTPCTransform::Row2X(GetFirstRow());
- Float_t a = -1/tan(fPsiLine);
- Float_t b = fDLine/sin(fPsiLine);
- Float_t yhit = a*xhit + b;
- xy[0] = xhit;
- xy[1] = yhit;
-}
+++ /dev/null
-// @(#) $Id$
-// origin hough/AliL3HoughTrack.h,v 1.8 Thu Mar 31 04:48:58 2005 UTC by cvetan
-
-#ifndef ALIHLTTPCHOUGHTRACK_H
-#define ALIHLTTPCHOUGHTRACK_H
-
-#include "AliHLTTPCTrack.h"
-
-class AliHLTTPCHoughTrack : public AliHLTTPCTrack {
-
- public:
- AliHLTTPCHoughTrack();
- virtual ~AliHLTTPCHoughTrack();
-
- virtual void Set(AliHLTTPCTrack *track);
- virtual Int_t Compare(const AliHLTTPCTrack *track) const;
-
- Bool_t IsHelix() const {return fIsHelix;}
- void UpdateToFirstRow();
- void SetTrackParameters(Double_t kappa,Double_t eangle,Int_t weight);
- void SetTrackParametersRow(Double_t alpha1,Double_t alpha2,Double_t eta,Int_t weight);
- void SetLineParameters(Double_t psi,Double_t D,Int_t weight,Int_t *rowrange,Int_t refrow);
-
- Int_t GetWeight() const {return fWeight;}
- Double_t GetPsiLine() const {return fPsiLine;}
- Double_t GetDLine() const {return fDLine;}
-
- Int_t GetEtaIndex() const {return fEtaIndex;}
- Double_t GetEta() const {return fEta;}
- Int_t GetSlice() const {return fSlice;}
- void GetLineCrossingPoint(Int_t padrow,Float_t *xy);
-
- Float_t GetBinX() const {return fBinX;}
- Float_t GetBinY() const {return fBinY;}
- Float_t GetSizeX() const {return fSizeX;}
- Float_t GetSizeY() const {return fSizeY;}
-
- void SetHelixTrue() {fIsHelix=kTRUE;}
- void SetSlice(Int_t slice) {fSlice=slice;}
- void SetEta(Double_t f);
- void SetWeight(Int_t i,Bool_t update=kFALSE) {if(update) fWeight+= i; else fWeight = i;}
- void SetEtaIndex(Int_t f) {fEtaIndex = f;}
- void SetBestMCid(Int_t f,Double_t mindist);
- void SetDLine(Double_t f) {fDLine=f;}
- void SetPsiLine(Double_t f) {fPsiLine=f;}
-
- void SetBinXY(Float_t binx,Float_t biny,Float_t sizex,Float_t sizey) {fBinX = binx; fBinY = biny; fSizeX = sizex; fSizeY = sizey;}
-
- private:
-
- Double_t fMinDist;//Minimum distance to a generated track while associating mc label
- Int_t fWeight;//Track weight
- Int_t fEtaIndex;//Eta slice index
- Double_t fEta;//Track Eta
- Int_t fSlice; //The slice where this track was found
-
- Double_t fDLine;//??
- Double_t fPsiLine;//??
-
- Bool_t fIsHelix;//Is the track helix or not?
-
- Float_t fBinX,fBinY,fSizeX,fSizeY;//Size and position of the hough space peak
-
- ClassDef(AliHLTTPCHoughTrack,1) //Track class for Hough tracklets
-
-};
-
-#endif
+++ /dev/null
-// $Id$
-// origin: src/AliL3TPCtracker.cxx 1.2 Fri Jun 10 04:25:00 2005 UTC by cvetan
-
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-/** @file AliHLTTPCHoughTracker.cxx
- @author Cvetan Cheshkov
- @date
- @brief Implementation of the HLT TPC hough transform tracker. */
-
-//-------------------------------------------------------------------------
-// Implementation of the HLT TPC hough transform tracker class
-//
-// It reads directly TPC digits using runloader and runs the HT
-// algorithm over them.
-// It stores the reconstructed hough tracks in the HLT ESD using the
-// the off-line AliESDtrack format.
-//
-// Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
-//-------------------------------------------------------------------------
-
-#include "AliESDEvent.h"
-#include "AliRunLoader.h"
-#include "AliHLTTPCHoughTracker.h"
-#include "AliHLTTPCHough.h"
-#include "AliLog.h"
-#include "AliHLTTPCTransform.h"
-#include "AliHLTTPCLog.h"
-
-/** ROOT macro for the implementation of ROOT specific class methods */
-ClassImp(AliHLTTPCHoughTracker);
-
-/**
- * creator
- * used for instantiation when the library is loaded dynamically
- */
-extern "C" AliTracker* CreateHLTTPCHoughTrackerInstance(AliRunLoader* runLoader)
-{
- return new AliHLTTPCHoughTracker(runLoader);
-}
-
-AliHLTTPCHoughTracker::AliHLTTPCHoughTracker(AliRunLoader *runLoader):AliTracker()
-{
- //--------------------------------------------------------------
- // Constructor
- //--------------------------------------------------------------
-
- if(AliHLTTPCTransform::GetVersion() == AliHLTTPCTransform::kVdefault) {
- Bool_t isinit=AliHLTTPCTransform::Init(runLoader);
- if(!isinit) AliWarning("Could not init AliHLTTPCTransform settings, using defaults!");
- }
-
- fRunLoader = runLoader;
-}
-
-Int_t AliHLTTPCHoughTracker::Clusters2Tracks(AliESDEvent *event)
-{
- //--------------------------------------------------------------------
- // This method reconstructs HLT TPC Hough tracks
- //--------------------------------------------------------------------
-
- if (!fRunLoader) {
- AliError("Missing runloader!");
- return kTRUE;
- }
- Int_t iEvent = fRunLoader->GetEventNumber();
-
- Float_t ptmin = 0.1*AliHLTTPCTransform::GetSolenoidField();
-
- Float_t zvertex = GetZ();
-
- AliInfo(Form("Hough Transform will run with ptmin=%f and zvertex=%f",ptmin,zvertex));
-
- // increase logging level temporarily to avoid bunches of info messages
- AliHLTTPCLog::TLogLevel loglevelbk=AliHLTTPCLog::fgLevel;
- AliHLTTPCLog::fgLevel=AliHLTTPCLog::kWarning;
-
- AliHLTTPCHough *hough = new AliHLTTPCHough();
-
- hough->SetThreshold(4);
- hough->CalcTransformerParams(ptmin);
- hough->SetPeakThreshold(70,-1);
- hough->SetRunLoader(fRunLoader);
- hough->Init("./", kFALSE, 100, kFALSE,4,0,0,zvertex);
- hough->SetAddHistograms();
-
- for(Int_t slice=0; slice<=35; slice++)
- {
- hough->ReadData(slice,iEvent);
- hough->Transform();
- hough->AddAllHistogramsRows();
- hough->FindTrackCandidatesRow();
- hough->AddTracks();
- }
-
- Int_t ntrk = hough->FillESD(event);
-
- Info("Clusters2Tracks","Number of found tracks: %d\n",ntrk);
-
- delete hough;
- AliHLTTPCLog::fgLevel=loglevelbk;
-
- return 0;
-}
+++ /dev/null
-// $Id$
-// origin: src/AliL3TPCtracker.h 1.1 Thu Mar 31 04:48:58 2005 UTC by cvetan
-
-#ifndef ALIL3TPCHOUGHTRACKER_H
-#define ALIL3TPCHOUGHTRACKER_H
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice */
-
-/** @file AliHLTTPCHoughTracker.h
- @author Cvetan Cheshkov
- @date
- @brief Implementation of the HLT TPC hough transform tracker. */
-
-//-------------------------------------------------------------------------
-// High Level Trigger TPC tracker
-// This class encapsulates the Hough transform HLT tracking
-// algorithm. It is used to call the algorithm inside the off-line
-// reconstruction chain. So far the tracker uses AliRunLoader to
-// to get the TPC digits. In the future all the references to
-// runloaders will be removed and the tracker will take as an input
-// the digits tree.
-//
-// Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
-//-------------------------------------------------------------------------
-
-#include "AliTracker.h"
-
-class AliRunLoader;
-class AliESDEvent;
-
-//-------------------------------------------------------------------------
-class AliHLTTPCHoughTracker : public AliTracker {
-public:
- AliHLTTPCHoughTracker(AliRunLoader *runLoader);
-
- Int_t Clusters2Tracks(AliESDEvent *event);
-
- Int_t PropagateBack(AliESDEvent */*event*/) {return 0;}
- Int_t RefitInward(AliESDEvent */*event*/) {return 0;}
- Int_t LoadClusters(TTree */*cf*/) {return 0;}
- void UnloadClusters() {return;}
-
- AliCluster *GetCluster(Int_t /*index*/) const {return NULL;}
-
-private:
- AliRunLoader *fRunLoader; // Pointer to the runloader
-
- ClassDef(AliHLTTPCHoughTracker,1) //HLT TPC Hough tracker
-};
-
-#endif
+++ /dev/null
-// @(#) $Id$
-// origin hough/AliL3HoughBaseTransformer.cxx,v 1.16 Tue Mar 22 13:11:58 2005 UTC by cvetan
-
-// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
-//*-- Copyright © ALICE HLT Group
-//-------------------------------------------------------------------------
-// Implementation of the AliHLTTPCHoughTransformer class
-// that is the base class for AliHLTHoughTransformer,
-// AliHLTHoughTransformerVhdl, AliHLTHoughTransformerGlobal,
-// AliHLTHoughTransformerRow
-//-------------------------------------------------------------------------
-
-#include "AliHLTTPCHoughTransformer.h"
-#include "AliTPCRawStream.h"
-
-/** \class AliHLTTPCHoughTransformer
-<pre>
-//_____________________________________________________________
-// AliHLTTPCHoughTransformer
-//
-// The base class for implementations of Hough Transform on ALICE TPC data.
-//
-// This is an abstract class, and is only meant to provide the interface
-// to the different implementations.
-//
-</pre>
-*/
-
-ClassImp(AliHLTTPCHoughTransformer)
-
-AliHLTTPCHoughTransformer::AliHLTTPCHoughTransformer()
- :
- fLastTransformer(NULL),
- fSlice(0),
- fPatch(0),
- fLastPatch(-1),
- fNEtaSegments(0),
- fEtaMin(0),
- fEtaMax(0),
- fLowerThreshold(0),
- fUpperThreshold(1023),
- fDigitRowData(NULL),
- fZVertex(0.0)
-{
- //Default constructor
-}
-
-AliHLTTPCHoughTransformer::AliHLTTPCHoughTransformer(Int_t slice,Int_t patch,Int_t netasegments,Float_t zvertex)
- :
- fLastTransformer(NULL),
- fSlice(0),
- fPatch(0),
- fLastPatch(-1),
- fNEtaSegments(0),
- fEtaMin(0),
- fEtaMax(0),
- fLowerThreshold(3),
- fUpperThreshold(1023),
- fDigitRowData(NULL),
- fZVertex(zvertex)
-{
- //normal ctor
- Init(slice,patch,netasegments);
-}
-
-AliHLTTPCHoughTransformer::~AliHLTTPCHoughTransformer()
-{
- //dtor
-}
-
-void AliHLTTPCHoughTransformer::Init(Int_t slice,Int_t patch,Int_t netasegments,Int_t /*n_seqs*/)
-{
- //Transformer init
- fSlice = slice;
- fPatch = patch;
- fLastPatch = -1;
- fNEtaSegments = netasegments;
- fEtaMin = 0;
- fEtaMax = fSlice < 18 ? 1. : -1.;
-}
+++ /dev/null
-// @(#) $Id$
-// origin hough/AliL3HoughBaseTransformer.h,v 1.22 Tue Mar 22 13:11:58 2005 UTC by cvetan
-
-#ifndef ALIHLTTPCHOUGHBASETRANSFORMER_H
-#define ALIHLTTPCHOUGHBASETRANSFORMER_H
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice */
-
-/** @file AliHLTTPCHoughTransformer.h
- @author Anders Vestbo
- @date
- @brief Base class for TPC hough tracking algorithms
- */
-
-//-------------------------------------------------------------------------
-// Class AliHLTTPCHoughTransformer
-// This is the base class for all the Hough Transformer tracking
-// algorithms for HLT.
-//-------------------------------------------------------------------------
-
-#include "AliHLTStdIncludes.h"
-#include "AliHLTTPCRootTypes.h"
-
-#ifdef do_mc
-const UInt_t MaxTrack=120;
-struct AliHLTTrackIndex {
- Int_t fLabel[MaxTrack];//MC label
- UChar_t fNHits[MaxTrack];//Number of different mc labels
- UChar_t fCurrentRow[MaxTrack];//Index of the current row while filling Hough space
-};
-typedef struct AliHLTTrackIndex AliHLTTrackIndex;
-#endif
-
-class AliHLTTPCDigitRowData;
-class AliHLTTPCHistogram;
-class AliTPCRawStream;
-
-/**
- * @class AliHLTTPCHoughTransformer
- * Base class for TPC hough tracking transfomers.
- */
-class AliHLTTPCHoughTransformer {
-
- public:
-
- /** standard constructor */
- AliHLTTPCHoughTransformer();
- /** constructor */
- AliHLTTPCHoughTransformer(Int_t slice,Int_t patch,Int_t netasegments,Float_t zvertex=0.0);
- /** standard destructor */
- virtual ~AliHLTTPCHoughTransformer();
-
- void SetInputData(UInt_t /*ndigits*/,AliHLTTPCDigitRowData *ptr) {fDigitRowData = ptr;}
-
- //this is for adaptave histograms
- virtual void CreateHistograms(Float_t /*ptmin*/,Float_t /*ptmax*/,Float_t /*pres*/,Int_t /*nybin*/,Float_t /*psi*/)
- {STDCERR<<"Adaptive histograms are not supported for this Transformer!"<<STDENDL;}
-
- virtual void CreateHistograms(Int_t nxbin,Float_t ptmin,Int_t nybin,Float_t phimin,Float_t phimax) = 0;
- virtual void CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,Int_t nybin,Float_t ymin,Float_t ymax) = 0;
-
- virtual void Reset() = 0;
- virtual void TransformCircle()
- {STDCERR<<"TransformCircle is not defined for this transformer!"<<STDENDL;}
- virtual void TransformCircle(Int_t */*rowRange*/,Int_t /*every*/)
- {STDCERR<<"TransformCircle is not defined for this transformer!"<<STDENDL;}
- virtual void TransformCircleC(Int_t */*rowRange*/,Int_t /*every*/)
- {STDCERR<<"TransformCircleC is not defined for this transformer!"<<STDENDL;}
- virtual void TransformLine(Int_t */*rowrange*/=0,Float_t */*phirange*/=0)
- {STDCERR<<"TransformLine is not defined for this Transformer!"<<STDENDL;}
- virtual void TransformLineC(Int_t */*rowrange*/,Float_t */*phirange*/)
- {STDCERR<<"TransformLineC is not defined for this Transformer!"<<STDENDL;}
-
- //Getters
- Int_t GetSlice() const {return fSlice;}
- Int_t GetPatch() const {return fPatch;}
- Int_t GetLastPatch() const {return fLastPatch;}
- AliHLTTPCHoughTransformer *GetLastTransfromer() const {return fLastTransformer;}
- Int_t GetNEtaSegments() const {return fNEtaSegments;}
- Int_t GetLowerThreshold() const {return fLowerThreshold;}
- Int_t GetUpperThreshold() const {return fUpperThreshold;}
- Double_t GetEtaMin() const {return fEtaMin;}
- Double_t GetEtaMax() const {return fEtaMax;}
- Float_t GetZVertex() const {return fZVertex;}
-
- AliHLTTPCDigitRowData *GetDataPointer() {return fDigitRowData;}
-
- virtual Int_t GetEtaIndex(Double_t eta) const = 0;
- virtual void GetEtaIndexes(Double_t /*eta*/,Int_t */*indexes*/) const
- {STDCERR<<"GetEtaIndexes not implemented for this Transformer class"<<STDENDL;}
- virtual AliHLTTPCHistogram *GetHistogram(Int_t etaindex) = 0;
- virtual Double_t GetEta(Int_t etaindex,Int_t slice) const = 0;
-
- virtual Int_t GetTrackID(Int_t /*etaindex*/,Double_t /*kappa*/,Double_t /*psi*/) const {
- STDCERR<<"GetTrackID not implemented for this Transformer class"<<STDENDL;
- return -1;
- }
-
- //setters
- virtual void Init(Int_t slice=0,Int_t patch=0,Int_t netasegments=100,Int_t nsegs=-1);
- void SetLowerThreshold(Int_t i) {fLowerThreshold = i;}
- void SetUpperThreshold(Int_t i) {fUpperThreshold = i;}
- void SetLastPatch(Int_t i) {fLastPatch = i;}
- void SetLastTransformer(AliHLTTPCHoughTransformer *transformer) {fLastTransformer = transformer;}
-
- virtual void SetTPCRawStream(AliTPCRawStream */*rawstream*/){};
-
- virtual void Print(){};
-
- protected:
-
- AliHLTTPCHoughTransformer *fLastTransformer;//Pointer to the previous hough transformer
-
- private:
- /** copy constructor prohibited */
- AliHLTTPCHoughTransformer(const AliHLTTPCHoughTransformer&);
- /** assignment operator prohibited */
- AliHLTTPCHoughTransformer& operator=(const AliHLTTPCHoughTransformer&);
-
- Int_t fSlice;//Index of the current slice being processed
- Int_t fPatch;//Index of the current patch being processed
- Int_t fLastPatch;//Index of the last processed patch
- Int_t fNEtaSegments;//Number of eta slices
- Double_t fEtaMin;//Minimum allowed eta
- Double_t fEtaMax;//Maximum allowed eta
- Int_t fLowerThreshold;//Lower threshold for digits amplitude
- Int_t fUpperThreshold;//Upper threshold for digits amplitude
-
- AliHLTTPCDigitRowData *fDigitRowData; //!
-
- Float_t fZVertex;//Z position of the primary vertex
-
- ClassDef(AliHLTTPCHoughTransformer,1) //Hough transformation base class
-
-};
-
-#endif
-
+++ /dev/null
-// @(#) $Id$
-// origin hough/AliL3HoughTransformerRow.cxx,v 1.21 Tue Mar 28 18:05:12 2006 UTC by alibrary
-
-/** @file AliHLTTPCHoughTransformerRow.cxx
- @author Cvetan Cheshkov <mailto:cvetan.cheshkov@cern.ch>
- @date
- @brief Implementation of fast HLT TPC hough transform tracking. */
-
-//_____________________________________________________________
-// AliHLTTPCHoughTransformerRow
-//
-// Impelementation of the so called "TPC rows Hough transformation" class
-//
-// Transforms the TPC data into the hough space and counts the missed TPC
-// rows corresponding to each track cnadidate - hough space bin
-
-#include "AliHLTStdIncludes.h"
-
-#include "AliHLTTPCLogging.h"
-#include "AliHLTTPCMemHandler.h"
-#include "AliHLTTPCTransform.h"
-#include "AliHLTTPCDigitData.h"
-#include "AliHLTTPCHistogramAdaptive.h"
-#include "AliHLTTPCHoughTrack.h"
-#include "AliHLTTPCHoughTransformerRow.h"
-#include "AliTPCRawStream.h"
-
-#if __GNUC__ >= 3
-using namespace std;
-#endif
-
-/** ROOT macro for the implementation of ROOT specific class methods */
-ClassImp(AliHLTTPCHoughTransformerRow);
-
-Float_t AliHLTTPCHoughTransformerRow::fgBeta1 = 1.0/AliHLTTPCTransform::Row2X(84);
-Float_t AliHLTTPCHoughTransformerRow::fgBeta2 = 1.0/(AliHLTTPCTransform::Row2X(158)*(1.0+tan(AliHLTTPCTransform::Pi()*10/180)*tan(AliHLTTPCTransform::Pi()*10/180)));
-Float_t AliHLTTPCHoughTransformerRow::fgDAlpha = 0.22;
-Float_t AliHLTTPCHoughTransformerRow::fgDEta = 0.40;
-Double_t AliHLTTPCHoughTransformerRow::fgEtaCalcParam1 = 1.0289;
-Double_t AliHLTTPCHoughTransformerRow::fgEtaCalcParam2 = 0.15192;
-Double_t AliHLTTPCHoughTransformerRow::fgEtaCalcParam3 = 1./(32.*600.*600.);
-
-AliHLTTPCHoughTransformerRow::AliHLTTPCHoughTransformerRow()
-{
- //Default constructor
- fParamSpace = 0;
-
- fGapCount = 0;
- fCurrentRowCount = 0;
-#ifdef do_mc
- fTrackID = 0;
-#endif
- fTrackNRows = 0;
- fTrackFirstRow = 0;
- fTrackLastRow = 0;
- fInitialGapCount = 0;
-
- fPrevBin = 0;
- fNextBin = 0;
- fNextRow = 0;
-
- fStartPadParams = 0;
- fEndPadParams = 0;
- fLUTr = 0;
- fLUTforwardZ = 0;
- fLUTbackwardZ = 0;
-
-}
-
-AliHLTTPCHoughTransformerRow::AliHLTTPCHoughTransformerRow(Int_t slice,Int_t patch,Int_t netasegments,Bool_t /*DoMC*/,Float_t zvertex) : AliHLTTPCHoughTransformer(slice,patch,netasegments,zvertex)
-{
- //Normal constructor
- fParamSpace = 0;
-
- fGapCount = 0;
- fCurrentRowCount = 0;
-#ifdef do_mc
- fTrackID = 0;
-#endif
-
- fTrackNRows = 0;
- fTrackFirstRow = 0;
- fTrackLastRow = 0;
- fInitialGapCount = 0;
-
- fPrevBin = 0;
- fNextBin = 0;
- fNextRow = 0;
-
- fStartPadParams = 0;
- fEndPadParams = 0;
- fLUTr = 0;
- fLUTforwardZ = 0;
- fLUTbackwardZ = 0;
-
-}
-
-AliHLTTPCHoughTransformerRow::~AliHLTTPCHoughTransformerRow()
-{
- //Destructor
- if(fLastTransformer) return;
- DeleteHistograms();
-#ifdef do_mc
- if(fTrackID)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fTrackID[i]) continue;
- delete fTrackID[i];
- }
- delete [] fTrackID;
- fTrackID = 0;
- }
-#endif
-
- if(fGapCount)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fGapCount[i]) continue;
- delete [] fGapCount[i];
- }
- delete [] fGapCount;
- fGapCount = 0;
- }
- if(fCurrentRowCount)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(fCurrentRowCount[i])
- delete [] fCurrentRowCount[i];
- }
- delete [] fCurrentRowCount;
- fCurrentRowCount = 0;
- }
- if(fTrackNRows)
- {
- delete [] fTrackNRows;
- fTrackNRows = 0;
- }
- if(fTrackFirstRow)
- {
- delete [] fTrackFirstRow;
- fTrackFirstRow = 0;
- }
- if(fTrackLastRow)
- {
- delete [] fTrackLastRow;
- fTrackLastRow = 0;
- }
- if(fInitialGapCount)
- {
- delete [] fInitialGapCount;
- fInitialGapCount = 0;
- }
- if(fPrevBin)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fPrevBin[i]) continue;
- delete [] fPrevBin[i];
- }
- delete [] fPrevBin;
- fPrevBin = 0;
- }
- if(fNextBin)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fNextBin[i]) continue;
- delete [] fNextBin[i];
- }
- delete [] fNextBin;
- fNextBin = 0;
- }
- if(fNextRow)
- {
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fNextRow[i]) continue;
- delete [] fNextRow[i];
- }
- delete [] fNextRow;
- fNextRow = 0;
- }
- if(fStartPadParams)
- {
- for(Int_t i = AliHLTTPCTransform::GetFirstRow(0); i<=AliHLTTPCTransform::GetLastRow(5); i++)
- {
- if(!fStartPadParams[i]) continue;
- delete [] fStartPadParams[i];
- }
- delete [] fStartPadParams;
- fStartPadParams = 0;
- }
- if(fEndPadParams)
- {
- for(Int_t i = AliHLTTPCTransform::GetFirstRow(0); i<=AliHLTTPCTransform::GetLastRow(5); i++)
- {
- if(!fEndPadParams[i]) continue;
- delete [] fEndPadParams[i];
- }
- delete [] fEndPadParams;
- fEndPadParams = 0;
- }
- if(fLUTr)
- {
- for(Int_t i = AliHLTTPCTransform::GetFirstRow(0); i<=AliHLTTPCTransform::GetLastRow(5); i++)
- {
- if(!fLUTr[i]) continue;
- delete [] fLUTr[i];
- }
- delete [] fLUTr;
- fLUTr = 0;
- }
- if(fLUTforwardZ)
- {
- delete[] fLUTforwardZ;
- fLUTforwardZ=0;
- }
- if(fLUTbackwardZ)
- {
- delete[] fLUTbackwardZ;
- fLUTbackwardZ=0;
- }
-}
-
-void AliHLTTPCHoughTransformerRow::DeleteHistograms()
-{
- // Clean up
- if(!fParamSpace)
- return;
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fParamSpace[i]) continue;
- delete fParamSpace[i];
- }
- delete [] fParamSpace;
-}
-
-void AliHLTTPCHoughTransformerRow::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
- if(fLastTransformer) {
- SetTransformerArrays((AliHLTTPCHoughTransformerRow *)fLastTransformer);
- return;
- }
- fParamSpace = new AliHLTTPCHistogram*[GetNEtaSegments()];
-
- Char_t histname[256];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- sprintf(histname,"paramspace_%d",i);
- fParamSpace[i] = new AliHLTTPCHistogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
- }
-#ifdef do_mc
- {
- AliHLTTPCHistogram *hist = fParamSpace[0];
- Int_t ncellsx = (hist->GetNbinsX()+3)/2;
- Int_t ncellsy = (hist->GetNbinsY()+3)/2;
- Int_t ncells = ncellsx*ncellsy;
- if(!fTrackID)
- {
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(AliHLTTrackIndex)<<" bytes to fTrackID"<<ENDLOG;
- fTrackID = new AliHLTTrackIndex*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fTrackID[i] = new AliHLTTrackIndex[ncells];
- }
- }
-#endif
- AliHLTTPCHistogram *hist = fParamSpace[0];
- Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
- if(!fGapCount)
- {
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fGapCount"<<ENDLOG;
- fGapCount = new UChar_t*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fGapCount[i] = new UChar_t[ncells];
- }
- if(!fCurrentRowCount)
- {
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fCurrentRowCount"<<ENDLOG;
- fCurrentRowCount = new UChar_t*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fCurrentRowCount[i] = new UChar_t[ncells];
- }
- if(!fPrevBin)
- {
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fPrevBin"<<ENDLOG;
- fPrevBin = new UChar_t*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fPrevBin[i] = new UChar_t[ncells];
- }
- if(!fNextBin)
- {
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fNextBin"<<ENDLOG;
- fNextBin = new UChar_t*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fNextBin[i] = new UChar_t[ncells];
- }
- Int_t ncellsy = hist->GetNbinsY()+2;
- if(!fNextRow)
- {
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<GetNEtaSegments()*ncellsy*sizeof(UChar_t)<<" bytes to fNextRow"<<ENDLOG;
- fNextRow = new UChar_t*[GetNEtaSegments()];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fNextRow[i] = new UChar_t[ncellsy];
- }
-
- if(!fTrackNRows)
- {
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackNRows"<<ENDLOG;
- fTrackNRows = new UChar_t[ncells];
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackFirstRow"<<ENDLOG;
- fTrackFirstRow = new UChar_t[ncells];
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackLastRow"<<ENDLOG;
- fTrackLastRow = new UChar_t[ncells];
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fInitialGapCount"<<ENDLOG;
- fInitialGapCount = new UChar_t[ncells];
-
-
- AliHLTTPCHoughTrack track;
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t xmiddle = (hist->GetNbinsX()+1)/2;
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
- Int_t nxbins = hist->GetNbinsX()+2;
- Int_t nxgrid = (hist->GetNbinsX()+3)/2+1;
- Int_t nygrid = hist->GetNbinsY()+3;
-
- AliHLTTrackLength *tracklength = new AliHLTTrackLength[nxgrid*nygrid];
- memset(tracklength,0,nxgrid*nygrid*sizeof(AliHLTTrackLength));
-
- for(Int_t ybin=ymin-1; ybin<=(ymax+1); ybin++)
- {
- for(Int_t xbin=xmin-1; xbin<=xmiddle; xbin++)
- {
- fTrackNRows[xbin + ybin*nxbins] = 255;
- for(Int_t deltay = 0; deltay <= 1; deltay++) {
- for(Int_t deltax = 0; deltax <= 1; deltax++) {
-
- AliHLTTrackLength *curtracklength = &tracklength[(xbin + deltax) + (ybin + deltay)*nxgrid];
- UInt_t maxfirstrow = 0;
- UInt_t maxlastrow = 0;
- Float_t maxtrackpt = 0;
- if(curtracklength->fIsFilled) {
- maxfirstrow = curtracklength->fFirstRow;
- maxlastrow = curtracklength->fLastRow;
- maxtrackpt = curtracklength->fTrackPt;
- }
- else {
- Float_t xtrack = hist->GetPreciseBinCenterX((Float_t)xbin+0.5*(Float_t)(2*deltax-1));
- Float_t ytrack = hist->GetPreciseBinCenterY((Float_t)ybin+0.5*(Float_t)(2*deltay-1));
-
- Float_t psi = atan((xtrack-ytrack)/(fgBeta1-fgBeta2));
- Float_t kappa = 2.0*(xtrack*cos(psi)-fgBeta1*sin(psi));
- track.SetTrackParameters(kappa,psi,1);
- maxtrackpt = track.GetPt();
- if(maxtrackpt < 0.9*0.1*AliHLTTPCTransform::GetSolenoidField())
- {
- maxfirstrow = maxlastrow = 0;
- curtracklength->fIsFilled = kTRUE;
- curtracklength->fFirstRow = maxfirstrow;
- curtracklength->fLastRow = maxlastrow;
- curtracklength->fTrackPt = maxtrackpt;
- }
- else
- {
- Bool_t firstrow = kFALSE;
- UInt_t curfirstrow = 0;
- UInt_t curlastrow = 0;
-
- Double_t centerx = track.GetCenterX();
- Double_t centery = track.GetCenterY();
- Double_t radius = track.GetRadius();
-
- for(Int_t j=AliHLTTPCTransform::GetFirstRow(0); j<=AliHLTTPCTransform::GetLastRow(5); j++)
- {
- Float_t hit[3];
- // if(!track.GetCrossingPoint(j,hit)) continue;
- hit[0] = AliHLTTPCTransform::Row2X(j);
- Double_t aa = (hit[0] - centerx)*(hit[0] - centerx);
- Double_t r2 = radius*radius;
- if(aa > r2)
- continue;
-
- Double_t aa2 = sqrt(r2 - aa);
- Double_t y1 = centery + aa2;
- Double_t y2 = centery - aa2;
- hit[1] = y1;
- if(fabs(y2) < fabs(y1)) hit[1] = y2;
-
- hit[2] = 0;
-
- AliHLTTPCTransform::LocHLT2Raw(hit,0,j);
- hit[1] += 0.5;
- if(hit[1]>=0 && hit[1]<AliHLTTPCTransform::GetNPads(j))
- {
- if(!firstrow) {
- curfirstrow = j;
- firstrow = kTRUE;
- }
- curlastrow = j;
- }
- else {
- if(firstrow) {
- firstrow = kFALSE;
- if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
- maxfirstrow = curfirstrow;
- maxlastrow = curlastrow;
- }
- }
- }
- }
- if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
- maxfirstrow = curfirstrow;
- maxlastrow = curlastrow;
- }
-
- curtracklength->fIsFilled = kTRUE;
- curtracklength->fFirstRow = maxfirstrow;
- curtracklength->fLastRow = maxlastrow;
- curtracklength->fTrackPt = maxtrackpt;
- }
- }
- if((maxlastrow-maxfirstrow) < fTrackNRows[xbin + ybin*nxbins]) {
- fTrackNRows[xbin + ybin*nxbins] = maxlastrow-maxfirstrow;
- fInitialGapCount[xbin + ybin*nxbins] = 1;
- if((maxlastrow-maxfirstrow+1)<=MIN_TRACK_LENGTH)
- fInitialGapCount[xbin + ybin*nxbins] = MAX_N_GAPS+1;
- if(maxtrackpt < 0.9*0.1*AliHLTTPCTransform::GetSolenoidField())
- fInitialGapCount[xbin + ybin*nxbins] = MAX_N_GAPS;
- fTrackFirstRow[xbin + ybin*nxbins] = maxfirstrow;
- fTrackLastRow[xbin + ybin*nxbins] = maxlastrow;
-
- Int_t xmirror = xmax - xbin + 1;
- Int_t ymirror = ymax - ybin + 1;
- fTrackNRows[xmirror + ymirror*nxbins] = fTrackNRows[xbin + ybin*nxbins];
- fInitialGapCount[xmirror + ymirror*nxbins] = fInitialGapCount[xbin + ybin*nxbins];
- fTrackFirstRow[xmirror + ymirror*nxbins] = fTrackFirstRow[xbin + ybin*nxbins];
- fTrackLastRow[xmirror + ymirror*nxbins] = fTrackLastRow[xbin + ybin*nxbins];
- }
- }
- }
- // cout<<" fTrackNRows "<<(Int_t)fInitialGapCount[xbin + ybin*nxbins]<<" "<<xbin<<" "<<ybin<<" "<<(Int_t)fTrackNRows[xbin + ybin*nxbins]<<" "<<(Int_t)fTrackFirstRow[xbin + ybin*nxbins]<<" "<<(Int_t)fTrackLastRow[xbin + ybin*nxbins]<<" "<<endl;
- }
- }
- delete [] tracklength;
- }
-
- if(!fStartPadParams)
- {
- Int_t nrows = AliHLTTPCTransform::GetLastRow(5) - AliHLTTPCTransform::GetFirstRow(0) + 1;
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating about "<<nrows*100*sizeof(AliHLTPadHoughParams)<<" bytes to fStartPadParams"<<ENDLOG;
- fStartPadParams = new AliHLTPadHoughParams*[nrows];
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating about "<<nrows*100*sizeof(AliHLTPadHoughParams)<<" bytes to fEndPadParams"<<ENDLOG;
- fEndPadParams = new AliHLTPadHoughParams*[nrows];
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating about "<<nrows*100*sizeof(Float_t)<<" bytes to fLUTr"<<ENDLOG;
- fLUTr = new Float_t*[nrows];
-
- Float_t beta1 = fgBeta1;
- Float_t beta2 = fgBeta2;
- Float_t beta1minusbeta2 = fgBeta1 - fgBeta2;
- Float_t ymin = hist->GetYmin();
- Float_t histbin = hist->GetBinWidthY();
- Float_t xmin = hist->GetXmin();
- Float_t xmax = hist->GetXmax();
- Float_t xbin = (xmax-xmin)/hist->GetNbinsX();
- Int_t firstbinx = hist->GetFirstXbin();
- Int_t lastbinx = hist->GetLastXbin();
- Int_t nbinx = hist->GetNbinsX()+2;
- Int_t firstbin = hist->GetFirstYbin();
- Int_t lastbin = hist->GetLastYbin();
- for(Int_t i=AliHLTTPCTransform::GetFirstRow(0); i<=AliHLTTPCTransform::GetLastRow(5); i++)
- {
- Int_t npads = AliHLTTPCTransform::GetNPads(i);
- Int_t ipatch = AliHLTTPCTransform::GetPatch(i);
- Double_t padpitch = AliHLTTPCTransform::GetPadPitchWidth(ipatch);
- Float_t x = AliHLTTPCTransform::Row2X(i);
- Float_t x2 = x*x;
-
- fStartPadParams[i] = new AliHLTPadHoughParams[npads];
- fEndPadParams[i] = new AliHLTPadHoughParams[npads];
- fLUTr[i] = new Float_t[npads];
- for(Int_t pad=0; pad<npads; pad++)
- {
- Float_t y = (pad-0.5*(npads-1))*padpitch;
- fLUTr[i][pad] = sqrt(x2 + y*y);
- Float_t starty = (pad-0.5*npads)*padpitch;
- Float_t r1 = x2 + starty*starty;
- Float_t xoverr1 = x/r1;
- Float_t startyoverr1 = starty/r1;
- Float_t endy = (pad-0.5*(npads-2))*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 alpha1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
- Float_t deltaalpha1 = b1*histbin/xbin;
- if(b1<0)
- alpha1 += deltaalpha1;
- Float_t alpha2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
- Float_t deltaalpha2 = b2*histbin/xbin;
- if(b2>=0)
- alpha2 += deltaalpha2;
-
- fStartPadParams[i][pad].fAlpha = alpha1;
- fStartPadParams[i][pad].fDeltaAlpha = deltaalpha1;
- fEndPadParams[i][pad].fAlpha = alpha2;
- fEndPadParams[i][pad].fDeltaAlpha = deltaalpha2;
-
- //Find the first and last bin rows to be filled
- Bool_t binfound1 = kFALSE;
- Bool_t binfound2 = kFALSE;
- Int_t firstbin1 = lastbin;
- Int_t firstbin2 = lastbin;
- Int_t lastbin1 = firstbin;
- Int_t lastbin2 = firstbin;
- for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
- {
- Int_t binx1 = 1 + (Int_t)alpha1;
- if(binx1<=lastbinx) {
- UChar_t initialgapcount;
- if(binx1>=firstbinx)
- initialgapcount = fInitialGapCount[binx1 + b*nbinx];
- else
- initialgapcount = fInitialGapCount[firstbinx + b*nbinx];
- if(initialgapcount != MAX_N_GAPS) {
- if(!binfound1) {
- firstbin1 = b;
- binfound1 = kTRUE;
- }
- lastbin1 = b;
- }
- }
- Int_t binx2 = 1 + (Int_t)alpha2;
- if(binx2>=firstbin) {
- UChar_t initialgapcount;
- if(binx2<=lastbinx)
- initialgapcount = fInitialGapCount[binx2 + b*nbinx];
- else
- initialgapcount = fInitialGapCount[lastbinx + b*nbinx];
- if(initialgapcount != MAX_N_GAPS) {
- if(!binfound2) {
- firstbin2 = b;
- binfound2 = kTRUE;
- }
- lastbin2 = b;
- }
- }
- }
- fStartPadParams[i][pad].fFirstBin=firstbin1;
- fStartPadParams[i][pad].fLastBin=lastbin1;
- fEndPadParams[i][pad].fFirstBin=firstbin2;
- fEndPadParams[i][pad].fLastBin=lastbin2;
- }
- }
- }
-
- //create lookup table for z of the digits
- if(!fLUTforwardZ)
- {
- Int_t ntimebins = AliHLTTPCTransform::GetNTimeBins();
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<ntimebins*sizeof(Float_t)<<" bytes to fLUTforwardZ"<<ENDLOG;
- fLUTforwardZ = new Float_t[ntimebins];
- LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
- <<"Transformer: Allocating "<<ntimebins*sizeof(Float_t)<<" bytes to fLUTbackwardZ"<<ENDLOG;
- fLUTbackwardZ = new Float_t[ntimebins];
- for(Int_t i=0; i<ntimebins; i++){
- Float_t z;
- z=AliHLTTPCTransform::GetZFast(0,i,GetZVertex());
- fLUTforwardZ[i]=z;
- z=AliHLTTPCTransform::GetZFast(18,i,GetZVertex());
- fLUTbackwardZ[i]=z;
- }
- }
-}
-
-void AliHLTTPCHoughTransformerRow::Reset()
-{
- //Reset all the histograms. Should be done when processing new slice
- if(fLastTransformer) return;
- if(!fParamSpace)
- {
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformer::Reset","Histograms")
- <<"No histograms to reset"<<ENDLOG;
- return;
- }
-
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fParamSpace[i]->Reset();
-
-#ifdef do_mc
- {
- AliHLTTPCHistogram *hist = fParamSpace[0];
- Int_t ncellsx = (hist->GetNbinsX()+3)/2;
- 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(AliHLTTrackIndex));
- }
-#endif
- AliHLTTPCHistogram *hist = fParamSpace[0];
- Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- memcpy(fGapCount[i],fInitialGapCount,ncells*sizeof(UChar_t));
- memcpy(fCurrentRowCount[i],fTrackFirstRow,ncells*sizeof(UChar_t));
- }
-}
-
-Int_t AliHLTTPCHoughTransformerRow::GetEtaIndex(Double_t eta) const
-{
- //Return the histogram index of the corresponding eta.
-
- Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
- Double_t index = (eta-GetEtaMin())/etaslice;
- return (Int_t)index;
-}
-
-inline AliHLTTPCHistogram *AliHLTTPCHoughTransformerRow::GetHistogram(Int_t etaindex)
-{
- // Return a pointer to the histogram which contains etaindex eta slice
- if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
- return 0;
- if(!fParamSpace[etaindex])
- return 0;
- return fParamSpace[etaindex];
-}
-
-Double_t AliHLTTPCHoughTransformerRow::GetEta(Int_t etaindex,Int_t /*slice*/) const
-{
- // Return eta calculated in the middle of the eta slice
- Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
- Double_t eta=0;
- eta=(Double_t)((etaindex+0.5)*etaslice);
- return eta;
-}
-
-void AliHLTTPCHoughTransformerRow::TransformCircle()
-{
- // This method contains the hough transformation
- // Depending on the option selected, it reads as an input
- // either the preloaded array with digits or directly raw
- // data stream (option fast_raw)
- if(GetDataPointer())
- TransformCircleFromDigitArray();
- else if(fTPCRawStream)
- TransformCircleFromRawStream();
-}
-void AliHLTTPCHoughTransformerRow::TransformCircleFromDigitArray()
-{
- //Do the Hough Transform
-
- //Load the parameters used by the fast calculation of eta
- Double_t etaparam1 = GetEtaCalcParam1();
- Double_t etaparam2 = GetEtaCalcParam2();
- Double_t etaparam3 = GetEtaCalcParam3();
-
- Int_t netasegments = GetNEtaSegments();
- Double_t etamin = GetEtaMin();
- Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
-
- Int_t lowerthreshold = GetLowerThreshold();
-
- //Assumes that all the etaslice histos are the same!
- AliHLTTPCHistogram *h = fParamSpace[0];
- Int_t firstbiny = h->GetFirstYbin();
- Int_t firstbinx = h->GetFirstXbin();
- Int_t lastbinx = h->GetLastXbin();
- Int_t nbinx = h->GetNbinsX()+2;
-
- UChar_t lastpad;
- Int_t lastetaindex=-1;
- AliHLTEtaRow *etaclust = new AliHLTEtaRow[netasegments];
-
- AliHLTTPCDigitRowData *tempPt = GetDataPointer();
- if(!tempPt)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughTransformer::TransformCircle","Data")
- <<"No input data "<<ENDLOG;
- return;
- }
-
- Int_t ipatch = GetPatch();
- Int_t ilastpatch = GetLastPatch();
- Int_t islice = GetSlice();
- Float_t *lutz;
- if(islice < 18)
- lutz = fLUTforwardZ;
- else
- lutz = fLUTbackwardZ;
-
- //Loop over the padrows:
- for(UChar_t i=AliHLTTPCTransform::GetFirstRow(ipatch); i<=AliHLTTPCTransform::GetLastRow(ipatch); i++)
- {
- lastpad = 255;
- //Flush eta clusters array
- memset(etaclust,0,netasegments*sizeof(AliHLTEtaRow));
-
- Float_t radius=0;
-
- //Get the data on this padrow:
- AliHLTTPCDigitData *digPt = tempPt->fDigitData;
- if((Int_t)i != (Int_t)tempPt->fRow)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughTransformerRow::TransformCircle","Data")
- <<"AliHLTTPCHoughTransform::TransformCircle : Mismatching padrow numbering "<<(Int_t)i<<" "<<(Int_t)tempPt->fRow<<ENDLOG;
- continue;
- }
- // cout<<" Starting row "<<i<<endl;
- //Loop over the data on this padrow:
- for(UInt_t j=0; j<tempPt->fNDigit; j++)
- {
- UShort_t charge = digPt[j].fCharge;
- if((Int_t)charge <= lowerthreshold)
- continue;
- UChar_t pad = digPt[j].fPad;
- UShort_t time = digPt[j].fTime;
-
- if(pad != lastpad)
- {
- radius = fLUTr[(Int_t)i][(Int_t)pad];
- lastetaindex = -1;
- }
-
- Float_t z = lutz[(Int_t)time];
- Double_t radiuscorr = radius*(1.+etaparam3*radius*radius);
- Double_t zovr = z/radiuscorr;
- Double_t eta = (etaparam1-etaparam2*fabs(zovr))*zovr;
- //Get the corresponding index, which determines which histogram to fill:
- Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
-
-#ifndef do_mc
- if(etaindex == lastetaindex) continue;
-#endif
- // cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
-
- if(etaindex < 0 || etaindex >= netasegments)
- continue;
-
- if(!etaclust[etaindex].fIsFound)
- {
- etaclust[etaindex].fStartPad = pad;
- etaclust[etaindex].fEndPad = pad;
- etaclust[etaindex].fIsFound = 1;
-#ifdef do_mc
- FillClusterMCLabels(digPt[j],&etaclust[etaindex]);
-#endif
- continue;
- }
- else
- {
- if(pad <= (etaclust[etaindex].fEndPad + 1))
- {
- etaclust[etaindex].fEndPad = pad;
-#ifdef do_mc
- FillClusterMCLabels(digPt[j],&etaclust[etaindex]);
-#endif
- }
- else
- {
- FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
-
- etaclust[etaindex].fStartPad = pad;
- etaclust[etaindex].fEndPad = pad;
-
-#ifdef do_mc
- memset(etaclust[etaindex].fMcLabels,0,MaxTrack);
- FillClusterMCLabels(digPt[j],&etaclust[etaindex]);
-#endif
- }
- }
- lastpad = pad;
- lastetaindex = etaindex;
- }
- //Write remaining clusters
- for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
- {
- //Check for empty row
- if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
-
- FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
- }
-
- //Move the data pointer to the next padrow:
- AliHLTTPCMemHandler::UpdateRowPointer(tempPt);
- }
-
- delete [] etaclust;
-}
-
-void AliHLTTPCHoughTransformerRow::TransformCircleFromRawStream()
-{
- //Do the Hough Transform
-
- //Load the parameters used by the fast calculation of eta
- Double_t etaparam1 = GetEtaCalcParam1();
- Double_t etaparam2 = GetEtaCalcParam2();
- Double_t etaparam3 = GetEtaCalcParam3();
-
- Int_t netasegments = GetNEtaSegments();
- Double_t etamin = GetEtaMin();
- Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
-
- Int_t lowerthreshold = GetLowerThreshold();
-
- //Assumes that all the etaslice histos are the same!
- AliHLTTPCHistogram *h = fParamSpace[0];
- Int_t firstbiny = h->GetFirstYbin();
- Int_t firstbinx = h->GetFirstXbin();
- Int_t lastbinx = h->GetLastXbin();
- Int_t nbinx = h->GetNbinsX()+2;
-
- Int_t lastetaindex = -1;
- AliHLTEtaRow *etaclust = new AliHLTEtaRow[netasegments];
-
- if(!fTPCRawStream)
- {
- LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughTransformer::TransformCircle","Data")
- <<"No input data "<<ENDLOG;
- return;
- }
-
- Int_t ipatch = GetPatch();
- UChar_t rowmin = AliHLTTPCTransform::GetFirstRowOnDDL(ipatch);
- UChar_t rowmax = AliHLTTPCTransform::GetLastRowOnDDL(ipatch);
- // Int_t ntimebins = AliHLTTPCTransform::GetNTimeBins();
- Int_t ilastpatch = GetLastPatch();
- Int_t islice = GetSlice();
- Float_t *lutz;
- if(islice < 18)
- lutz = fLUTforwardZ;
- else
- lutz = fLUTbackwardZ;
-
- //Flush eta clusters array
- memset(etaclust,0,netasegments*sizeof(AliHLTEtaRow));
-
- UChar_t i=0;
- Int_t npads=0;
- Float_t radius=0;
- UChar_t pad=0;
-
- //Loop over the rawdata:
- while (fTPCRawStream->Next()) {
-
- if(fTPCRawStream->IsNewSector() || fTPCRawStream->IsNewRow()) {
-
- //Write remaining clusters
- for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
- {
- //Check for empty row
- if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
-
- FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
- }
-
- Int_t sector=fTPCRawStream->GetSector();
- Int_t row=fTPCRawStream->GetRow();
- Int_t slice,srow;
- AliHLTTPCTransform::Sector2Slice(slice,srow,sector,row);
- if(slice!=islice){
- LOG(AliHLTTPCLog::kError,"AliHLTDDLDataFileHandler::DDLDigits2Memory","Slice")
- <<AliHLTTPCLog::kDec<<"Found slice "<<slice<<", expected "<<islice<<ENDLOG;
- continue;
- }
-
- i=(UChar_t)srow;
- npads = AliHLTTPCTransform::GetNPads(srow)-1;
-
- //Flush eta clusters array
- memset(etaclust,0,netasegments*sizeof(AliHLTEtaRow));
-
- radius=0;
-
- }
-
- if((i<rowmin)||(i>rowmax))continue;
-
- // cout<<" Starting row "<<(UInt_t)i<<endl;
- //Loop over the data on this padrow:
- if(fTPCRawStream->IsNewRow() || fTPCRawStream->IsNewPad()) {
- pad=fTPCRawStream->GetPad();
- /*
- if((pad<0)||(pad>=(npads+1))){
- LOG(AliHLTTPCLog::kError,"AliHLTDDLDataFileHandler::DDLDigits2Memory","Pad")
- <<AliHLTTPCLog::kDec<<"Pad value out of bounds "<<pad<<" "
- <<npads+1<<ENDLOG;
- continue;
- }
- */
- radius = fLUTr[(Int_t)i][(Int_t)pad];
- lastetaindex = -1;
- }
-
- UShort_t time=fTPCRawStream->GetTime();
- /*
- if((time<0)||(time>=ntimebins)){
- LOG(AliHLTTPCLog::kError,"AliHLTDDLDataFileHandler::DDLDigits2Memory","Time")
- <<AliHLTTPCLog::kDec<<"Time out of bounds "<<time<<" "
- <<AliHLTTPCTransform::GetNTimeBins()<<ENDLOG;
- continue;
- }
- */
-
- if(fTPCRawStream->GetSignal() <= lowerthreshold)
- continue;
-
- Float_t z = lutz[(Int_t)time];
- Double_t radiuscorr = radius*(1.+etaparam3*radius*radius);
- Double_t zovr = z/radiuscorr;
- Double_t eta = (etaparam1-etaparam2*fabs(zovr))*zovr;
- //Get the corresponding index, which determines which histogram to fill:
- Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
-
-#ifndef do_mc
- if(etaindex == lastetaindex) continue;
-#endif
- // cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
-
- if(etaindex < 0 || etaindex >= netasegments)
- continue;
-
- if(!etaclust[etaindex].fIsFound)
- {
- etaclust[etaindex].fStartPad = pad;
- etaclust[etaindex].fEndPad = pad;
- etaclust[etaindex].fIsFound = 1;
- continue;
- }
- else
- {
- if(pad <= (etaclust[etaindex].fEndPad + 1))
- {
- etaclust[etaindex].fEndPad = pad;
- }
- else
- {
- FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
-
- etaclust[etaindex].fStartPad = pad;
- etaclust[etaindex].fEndPad = pad;
-
- }
- }
- lastetaindex = etaindex;
- }
-
- //Write remaining clusters
- for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
- {
- //Check for empty row
- if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
-
- FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
- }
-
- delete [] etaclust;
-}
-
-#ifndef do_mc
-Int_t AliHLTTPCHoughTransformerRow::GetTrackID(Int_t /*etaindex*/,Double_t /*alpha1*/,Double_t /*alpha2*/) const
-{
- // Does nothing if do_mc undefined
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::GetTrackID","Data")
- <<"Flag switched off"<<ENDLOG;
- return -1;
-#else
-Int_t AliHLTTPCHoughTransformerRow::GetTrackID(Int_t etaindex,Double_t alpha1,Double_t alpha2) const
-{
- // Returns the MC label for a given peak found in the Hough space
- if(etaindex < 0 || etaindex > GetNEtaSegments())
- {
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::GetTrackID","Data")
- <<"Wrong etaindex "<<etaindex<<ENDLOG;
- return -1;
- }
- AliHLTTPCHistogram *hist = fParamSpace[etaindex];
- Int_t bin = hist->FindLabelBin(alpha1,alpha2);
- if(bin==-1) {
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::GetTrackID()","")
- <<"Track candidate outside Hough space boundaries: "<<alpha1<<" "<<alpha2<<ENDLOG;
- return -1;
- }
- Int_t label=-1;
- Int_t max=0;
- for(UInt_t i=0; i<(MaxTrack-1); i++)
- {
- Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
- if(nhits == 0) break;
- if(nhits > max)
- {
- max = nhits;
- label = fTrackID[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[etaindex][bin].fNHits[i];
- if(nhits == 0) break;
- if(nhits > max2)
- {
- if(fTrackID[etaindex][bin].fLabel[i]!=label) {
- max2 = nhits;
- label2 = fTrackID[etaindex][bin].fLabel[i];
- }
- }
- }
- if(max2 !=0 ) {
- LOG(AliHLTTPCLog::kDebug,"AliHLTTPCHoughTransformerRow::GetTrackID()","")
- <<" TrackID"<<" label "<<label<<" max "<<max<<" label2 "<<label2<<" max2 "<<max2<<" "<<(Float_t)max2/(Float_t)max<<" "<<fTrackID[etaindex][bin].fLabel[MaxTrack-1]<<" "<<(Int_t)fTrackID[etaindex][bin].fNHits[MaxTrack-1]<<ENDLOG;
- }
- return label;
-#endif
-}
-
-Int_t AliHLTTPCHoughTransformerRow::GetTrackLength(Double_t alpha1,Double_t alpha2,Int_t *rows) const
-{
- // Returns the track length for a given peak found in the Hough space
-
- AliHLTTPCHistogram *hist = fParamSpace[0];
- Int_t bin = hist->FindBin(alpha1,alpha2);
- if(bin==-1) {
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::GetTrackLength()","")
- <<"Track candidate outside Hough space boundaries: "<<alpha1<<" "<<alpha2<<ENDLOG;
- return -1;
- }
- rows[0] = fTrackFirstRow[bin];
- rows[1] = fTrackLastRow[bin];
-
- return 0;
-}
-
-inline void AliHLTTPCHoughTransformerRow::FillClusterRow(UChar_t i,Int_t binx1,Int_t binx2,UChar_t *ngaps2,UChar_t *currentrow2,UChar_t *lastrow2
-#ifdef do_mc
- ,AliHLTEtaRow etaclust,AliHLTTrackIndex *trackid
-#endif
- )
-{
- // The method is a part of the fast hough transform.
- // It fills one row of the hough space.
- // It is called by FillCluster() method inside the
- // loop over alpha2 bins
-
- for(Int_t bin=binx1;bin<=binx2;bin++)
- {
- if(ngaps2[bin] < MAX_N_GAPS) {
- if(i < lastrow2[bin] && i > currentrow2[bin]) {
- ngaps2[bin] += (i-currentrow2[bin]-1);
- currentrow2[bin]=i;
- }
-#ifdef do_mc
- if(i < lastrow2[bin] && i >= currentrow2[bin]) {
- for(UInt_t t=0;t<(MaxTrack-1); t++)
- {
- Int_t label = etaclust.fMcLabels[t];
- if(label == 0) break;
- UInt_t c;
- Int_t tempbin2 = (Int_t)(bin/2);
- for(c=0; c<(MaxTrack-1); c++)
- if(trackid[tempbin2].fLabel[c] == label || trackid[tempbin2].fNHits[c] == 0)
- break;
- trackid[tempbin2].fLabel[c] = label;
- if(trackid[tempbin2].fCurrentRow[c] != i) {
- trackid[tempbin2].fNHits[c]++;
- trackid[tempbin2].fCurrentRow[c] = i;
- }
- }
- }
-#endif
- }
- }
-
-}
-
-inline void AliHLTTPCHoughTransformerRow::FillCluster(UChar_t i,Int_t etaindex,AliHLTEtaRow *etaclust,Int_t ilastpatch,Int_t firstbinx,Int_t lastbinx,Int_t nbinx,Int_t firstbiny)
-{
- // The method is a part of the fast hough transform.
- // It fills a TPC cluster into the hough space.
-
- UChar_t *ngaps = fGapCount[etaindex];
- UChar_t *currentrow = fCurrentRowCount[etaindex];
- UChar_t *lastrow = fTrackLastRow;
- UChar_t *prevbin = fPrevBin[etaindex];
- UChar_t *nextbin = fNextBin[etaindex];
- UChar_t *nextrow = fNextRow[etaindex];
-#ifdef do_mc
- AliHLTTrackIndex *trackid = fTrackID[etaindex];
-#endif
-
- //Do the transformation:
- AliHLTPadHoughParams *startparams = &fStartPadParams[(Int_t)i][etaclust[etaindex].fStartPad];
- AliHLTPadHoughParams *endparams = &fEndPadParams[(Int_t)i][etaclust[etaindex].fEndPad];
-
- Float_t alpha1 = startparams->fAlpha;
- Float_t deltaalpha1 = startparams->fDeltaAlpha;
- Float_t alpha2 = endparams->fAlpha;
- Float_t deltaalpha2 = endparams->fDeltaAlpha;
-
- Int_t firstbin1 = startparams->fFirstBin;
- Int_t firstbin2 = endparams->fFirstBin;
- Int_t firstbin = firstbin1;
- if(firstbin>firstbin2) firstbin = firstbin2;
-
- Int_t lastbin1 = startparams->fLastBin;
- Int_t lastbin2 = endparams->fLastBin;
- Int_t lastbin = lastbin1;
- if(lastbin<lastbin2) lastbin = lastbin2;
-
- alpha1 += (firstbin-firstbiny)*deltaalpha1;
- alpha2 += (firstbin-firstbiny)*deltaalpha2;
-
- //Fill the histogram along the alpha2 range
- if(ilastpatch == -1) {
- for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
- {
- Int_t binx1 = 1 + (Int_t)alpha1;
- if(binx1>lastbinx) continue;
- if(binx1<firstbinx) binx1 = firstbinx;
- Int_t binx2 = 1 + (Int_t)alpha2;
- if(binx2<firstbinx) continue;
- if(binx2>lastbinx) binx2 = lastbinx;
-#ifdef do_mc
- if(binx2<binx1) {
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::TransformCircle()","")
- <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<etaclust[etaindex].fStartPad<<" "<<etaclust[etaindex].fEndPad<<ENDLOG;
- }
-#endif
- Int_t tempbin = b*nbinx;
- UChar_t *ngaps2 = ngaps + tempbin;
- UChar_t *currentrow2 = currentrow + tempbin;
- UChar_t *lastrow2 = lastrow + tempbin;
-#ifdef do_mc
- Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2));
- AliHLTTrackIndex *trackid2 = trackid + tempbin2;
-#endif
- FillClusterRow(i,binx1,binx2,ngaps2,currentrow2,lastrow2
-#ifdef do_mc
- ,etaclust[etaindex],trackid2
-#endif
- );
- }
- }
- else {
- for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
- {
- Int_t binx1 = 1 + (Int_t)alpha1;
- if(binx1>lastbinx) continue;
- if(binx1<firstbinx) binx1 = firstbinx;
- Int_t binx2 = 1 + (Int_t)alpha2;
- if(binx2<firstbinx) continue;
- if(binx2>lastbinx) binx2 = lastbinx;
-#ifdef do_mc
- if(binx2<binx1) {
- LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::TransformCircle()","")
- <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<etaclust[etaindex].fStartPad<<" "<<etaclust[etaindex].fEndPad<<ENDLOG;
- }
-#endif
- if(nextrow[b] > b) {
- Int_t deltab = (nextrow[b] - b - 1);
- b += deltab;
- alpha1 += deltaalpha1*deltab;
- alpha2 += deltaalpha2*deltab;
- continue;
- }
- Int_t tempbin = b*nbinx;
- binx1 = (UInt_t)nextbin[binx1 + tempbin];
- binx2 = (UInt_t)prevbin[binx2 + tempbin];
- if(binx2<binx1) continue;
- UChar_t *ngaps2 = ngaps + tempbin;
- UChar_t *currentrow2 = currentrow + tempbin;
- UChar_t *lastrow2 = lastrow + tempbin;
-#ifdef do_mc
- Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2));
- AliHLTTrackIndex *trackid2 = trackid + tempbin2;
-#endif
- FillClusterRow(i,binx1,binx2,ngaps2,currentrow2,lastrow2
-#ifdef do_mc
- ,etaclust[etaindex],trackid2
-#endif
- );
- }
- }
-}
-
-#ifdef do_mc
-inline void AliHLTTPCHoughTransformerRow::FillClusterMCLabels(AliHLTTPCDigitData digpt,AliHLTEtaRow *etaclust)
-{
- // The method is a part of the fast hough transform.
- // It fills the MC labels of a TPC cluster into a
- // special hough space array.
- for(Int_t t=0; t<3; t++)
- {
- Int_t label = digpt.fTrackID[t];
- if(label < 0) break;
- UInt_t c;
- for(c=0; c<(MaxTrack-1); c++)
- if(etaclust->fMcLabels[c] == label || etaclust->fMcLabels[c] == 0)
- break;
-
- etaclust->fMcLabels[c] = label;
- }
-}
-#endif
-
-void AliHLTTPCHoughTransformerRow::SetTransformerArrays(AliHLTTPCHoughTransformerRow *tr)
-{
- // In case of sequential filling of the hough space, the method is used to
- // transmit the pointers to the hough arrays from one transformer to the
- // following one.
-
- fGapCount = tr->fGapCount;
- fCurrentRowCount = tr->fCurrentRowCount;
-#ifdef do_mc
- fTrackID = tr->fTrackID;
-#endif
- fTrackNRows = tr->fTrackNRows;
- fTrackFirstRow = tr->fTrackFirstRow;
- fTrackLastRow = tr->fTrackLastRow;
- fInitialGapCount = tr->fInitialGapCount;
-
- fPrevBin = tr->fPrevBin;
- fNextBin = tr->fNextBin;
- fNextRow = tr->fNextRow;
-
- fStartPadParams = tr->fStartPadParams;
- fEndPadParams = tr->fEndPadParams;
- fLUTr = tr->fLUTr;
- fLUTforwardZ = tr->fLUTforwardZ;
- fLUTbackwardZ = tr->fLUTbackwardZ;
-
- fParamSpace = tr->fParamSpace;
-
- return;
-}
+++ /dev/null
-// @(#) $Id$
-// origin hough/AliL3HoughTransformerRow.h,v 1.15 Sun Apr 30 16:37:32 2006 UTC by hristov
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice */
-
-/** @file AliHLTTPCHoughTransformerRow.h
- @author Cvetan Cheshkov
- @date
- @brief Implementation of fast HLT TPC hough transform tracking. */
-
-#ifndef ALIHLTTPCHOUGHTRANSFORMERROW_H
-#define ALIHLTTPCHOUGHTRANSFORMERROW_H
-
-#include "AliHLTTPCRootTypes.h"
-#include "AliHLTTPCHoughTransformer.h"
-
-#define MAX_N_GAPS 5
-#define MIN_TRACK_LENGTH 70
-
-class AliHLTTPCDigitData;
-class AliHLTTPCHistogram;
-
-class AliHLTTPCHoughTransformerRow : public AliHLTTPCHoughTransformer {
-
- public:
- /** standard constructor */
- AliHLTTPCHoughTransformerRow();
- /** constructor */
- AliHLTTPCHoughTransformerRow(Int_t slice,Int_t patch,Int_t netasegments,Bool_t DoMC=kFALSE,Float_t zvertex=0.0);
- /** standard destructor */
- virtual ~AliHLTTPCHoughTransformerRow();
-
- struct AliHLTEtaRow {
- UChar_t fStartPad; //First pad in the cluster
- UChar_t fEndPad; //Last pad in the cluster
- Bool_t fIsFound; //Is the cluster already found
-#ifdef do_mc
- Int_t fMcLabels[MaxTrack]; //Array to store mc labels inside cluster
-#endif
- };
-
- struct AliHLTPadHoughParams {
- // Parameters which represent given pad in the hough space
- // Used in order to avoid as much as possible floating
- // point operations during the hough transform
- Float_t fAlpha; // Starting value for the hough parameter alpha1
- Float_t fDeltaAlpha; // Slope of alpha1
- Int_t fFirstBin; // First alpha2 bin to be filled
- Int_t fLastBin; // Last alpha2 bin to be filled
- };
-
- struct AliHLTTrackLength {
- // Structure is used for temporarely storage of the LUT
- // which contains the track lengths associated to each hough
- // space bin
- Bool_t fIsFilled; // Is bin already filled?
- UInt_t fFirstRow; // First TPC row crossed by the track
- UInt_t fLastRow; // Last TPC row crossed by the track
- Float_t fTrackPt; // Pt of the track
- };
-
-
-
- void CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t pres,Int_t nybin,Float_t psi) {
- AliHLTTPCHoughTransformer::CreateHistograms(ptmin,ptmax,pres,nybin,psi);
- }
- void CreateHistograms(Int_t /*nxbin*/,Float_t /*ptmin*/,Int_t /*nybin*/,Float_t /*phimin*/,Float_t /*phimax*/)
- {STDCERR<<"This method for creation of parameter space histograms is not supported for this Transformer!"<<STDENDL;}
- void CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
- Int_t nybin,Float_t ymin,Float_t ymax);
- void Reset();
- void TransformCircle();
- void TransformCircle(Int_t *rowRange,Int_t every) {
- AliHLTTPCHoughTransformer::TransformCircle(rowRange,every);
- }
-
- Int_t GetEtaIndex(Double_t eta) const;
- AliHLTTPCHistogram *GetHistogram(Int_t etaindex);
- Double_t GetEta(Int_t etaindex,Int_t slice) const;
- Int_t GetTrackID(Int_t etaindex,Double_t alpha1,Double_t alpha2) const;
- Int_t GetTrackLength(Double_t alpha1,Double_t alpha2,Int_t *rows) const;
- UChar_t *GetGapCount(Int_t etaindex) const { return fGapCount[etaindex]; }
- UChar_t *GetCurrentRowCount(Int_t etaindex) const { return fCurrentRowCount[etaindex]; }
- UChar_t *GetPrevBin(Int_t etaindex) const { return fPrevBin[etaindex]; }
- UChar_t *GetNextBin(Int_t etaindex) const { return fNextBin[etaindex]; }
- UChar_t *GetNextRow(Int_t etaindex) const { return fNextRow[etaindex]; }
- UChar_t *GetTrackNRows() const { return fTrackNRows; }
- UChar_t *GetTrackFirstRow() const { return fTrackFirstRow; }
- UChar_t *GetTrackLastRow() const { return fTrackLastRow; }
- static Float_t GetBeta1() {return fgBeta1;}
- static Float_t GetBeta2() {return fgBeta2;}
- static Float_t GetDAlpha() {return fgDAlpha;}
- static Float_t GetDEta() {return fgDEta;}
- static Double_t GetEtaCalcParam1() {return fgEtaCalcParam1;}
- static Double_t GetEtaCalcParam2() {return fgEtaCalcParam2;}
- static Double_t GetEtaCalcParam3() {return fgEtaCalcParam3;}
-
- void SetTPCRawStream(AliTPCRawStream *rawstream) {fTPCRawStream=rawstream;}
-
- private:
- /** copy constructor prohibited */
- AliHLTTPCHoughTransformerRow(const AliHLTTPCHoughTransformerRow&);
- /** assignment operator prohibited */
- AliHLTTPCHoughTransformerRow& operator=(const AliHLTTPCHoughTransformerRow&);
-
- UChar_t **fGapCount; //!
- UChar_t **fCurrentRowCount; //!
-#ifdef do_mc
- AliHLTTrackIndex **fTrackID; //!
-#endif
-
- UChar_t *fTrackNRows; //!
- UChar_t *fTrackFirstRow; //!
- UChar_t *fTrackLastRow; //!
- UChar_t *fInitialGapCount; //!
-
- UChar_t **fPrevBin; //!
- UChar_t **fNextBin; //!
- UChar_t **fNextRow; //!
-
- AliHLTPadHoughParams **fStartPadParams; //!
- AliHLTPadHoughParams **fEndPadParams; //!
- Float_t **fLUTr; //!
-
- Float_t *fLUTforwardZ; //!
- Float_t *fLUTbackwardZ; //!
-
- AliHLTTPCHistogram **fParamSpace; //!
-
- void TransformCircleFromDigitArray();
- void TransformCircleFromRawStream();
-
- void DeleteHistograms(); //Method to clean up the histograms containing Hough space
-
- inline void FillClusterRow(UChar_t i,Int_t binx1,Int_t binx2,UChar_t *ngaps2,UChar_t *currentrow2,UChar_t *lastrow2
-#ifdef do_mc
- ,AliHLTEtaRow etaclust,AliHLTTrackIndex *trackid
-#endif
- );
- inline void FillCluster(UChar_t i,Int_t etaindex,AliHLTEtaRow *etaclust,Int_t ilastpatch,Int_t firstbinx,Int_t lastbinx,Int_t nbinx,Int_t firstbiny);
-#ifdef do_mc
- inline void FillClusterMCLabels(AliHLTDigitData digpt,AliHLTEtaRow *etaclust);
-#endif
-
- void SetTransformerArrays(AliHLTTPCHoughTransformerRow *tr);
-
- static Float_t fgBeta1,fgBeta2; // Two curves which define the Hough space
- static Float_t fgDAlpha, fgDEta; // Correlation factor between Hough space bin size and resolution
- static Double_t fgEtaCalcParam1, fgEtaCalcParam2; // Parameters used for fast calculation of eta during the binning of Hough space
- static Double_t fgEtaCalcParam3; // Parameter used during the eta binning of the Hough Space in order to account for finite track radii
-
- AliTPCRawStream *fTPCRawStream; // Pointer to the raw stream in case of fast reading of the raw data (fast_raw flag)
-
- ClassDef(AliHLTTPCHoughTransformerRow,1) //TPC Rows Hough transformation class
-
-};
-
-#endif
-
-
-
-
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff