This commit was generated by cvs2svn to compensate for changes in r3174,
authorvestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Mar 2001 12:51:56 +0000 (12:51 +0000)
committervestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Mar 2001 12:51:56 +0000 (12:51 +0000)
which included commits to RCS files with non-trunk default branches.

23 files changed:
HLT/hough/AliL3Defs.h [new file with mode: 0644]
HLT/hough/AliL3Histogram.cxx [new file with mode: 0644]
HLT/hough/AliL3Histogram.h [new file with mode: 0644]
HLT/hough/AliL3HoughEval.cxx [new file with mode: 0644]
HLT/hough/AliL3HoughEval.h [new file with mode: 0644]
HLT/hough/AliL3HoughLinkDef.h [new file with mode: 0644]
HLT/hough/AliL3HoughMaxFinder.cxx [new file with mode: 0644]
HLT/hough/AliL3HoughMaxFinder.h [new file with mode: 0644]
HLT/hough/AliL3HoughMerge.cxx [new file with mode: 0644]
HLT/hough/AliL3HoughMerge.h [new file with mode: 0644]
HLT/hough/AliL3HoughPixel.cxx [new file with mode: 0644]
HLT/hough/AliL3HoughPixel.h [new file with mode: 0644]
HLT/hough/AliL3HoughTransformer.cxx [new file with mode: 0644]
HLT/hough/AliL3HoughTransformer.h [new file with mode: 0644]
HLT/hough/AliL3Transform.cxx [new file with mode: 0644]
HLT/hough/AliL3Transform.h [new file with mode: 0644]
HLT/hough/Makefile [new file with mode: 0644]
HLT/hough/hough.C [new file with mode: 0644]
HLT/hough/hough_line.C [new file with mode: 0644]
HLT/hough/hough_line_merge.C [new file with mode: 0644]
HLT/hough/hough_merge.C [new file with mode: 0644]
HLT/hough/hough_mergehistos.C [new file with mode: 0644]
HLT/hough/rootlogon.C [new file with mode: 0644]

diff --git a/HLT/hough/AliL3Defs.h b/HLT/hough/AliL3Defs.h
new file mode 100644 (file)
index 0000000..c9bfdf0
--- /dev/null
@@ -0,0 +1,12 @@
+#ifndef _ALIL3DEFS_H_
+#define _ALIL3DEFS_H_
+
+#include <Rtypes.h>
+
+const Int_t NRowsSlice = 174;
+const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
+const Double_t Pi = 3.14159265358979323846;
+const Double_t ToRad = Pi/180.;
+
+
+#endif /* _ALIL3DEFS_H_ */
diff --git a/HLT/hough/AliL3Histogram.cxx b/HLT/hough/AliL3Histogram.cxx
new file mode 100644 (file)
index 0000000..ccc5460
--- /dev/null
@@ -0,0 +1,17 @@
+
+#include "AliL3Histogram.h"
+
+ClassImp(AliL3Histogram)
+
+  
+AliL3Histogram::AliL3Histogram()
+{
+
+}
+
+
+AliL3Histogram::~AliL3Histogram()
+{
+  //Destructor
+
+}
diff --git a/HLT/hough/AliL3Histogram.h b/HLT/hough/AliL3Histogram.h
new file mode 100644 (file)
index 0000000..b31e1c4
--- /dev/null
@@ -0,0 +1,20 @@
+#ifndef ALIL3_HOUGHHISTOGRAM
+#define ALIL3_HOUGHHISTOGRAM
+
+#include "AliL3RootTypes.h"
+
+class AliL3HoughTransformer : public TH2F {
+  
+ private:
+
+
+  
+ public:
+  AliL3Histogram();
+  
+
+  ClassDef(AliL3Histogram,1)
+
+};
+
+#endif
diff --git a/HLT/hough/AliL3HoughEval.cxx b/HLT/hough/AliL3HoughEval.cxx
new file mode 100644 (file)
index 0000000..9e0aaca
--- /dev/null
@@ -0,0 +1,205 @@
+#include <TClonesArray.h>
+#include <TH2.h>
+#include <TH1.h>
+#include <TFile.h>
+#include <AliRun.h>
+#include <TParticle.h>
+
+#include "AliL3TrackArray.h"
+#include "AliL3Transform.h"
+#include "AliL3HoughTransformer.h"
+#include "AliL3Defs.h"
+#include "AliL3HoughTrack.h"
+#include "AliL3HoughEval.h"
+
+ClassImp(AliL3HoughEval)
+
+
+AliL3HoughEval::AliL3HoughEval()
+{
+  //Default constructor
+  fTransform = new AliL3Transform();
+  fHoughTransformer = NULL;
+}
+
+
+AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
+{
+  //Constructor
+  fHoughTransformer = transformer;
+  fTransform = new AliL3Transform();
+}
+
+
+AliL3HoughEval::~AliL3HoughEval()
+{
+  //Destructor
+  if(fTransform)
+    delete fTransform;
+}
+
+
+void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist,TH2F *fake)
+{
+  //Look at rawdata along the road specified by the track candidates.
+
+  if(!fHoughTransformer)
+    {
+      printf("AliL3HoughEval: No transformer object\n");
+      return;
+    }
+  
+  AliL3Digits *pixel;
+  Int_t patch = fHoughTransformer->fPatch;
+  Int_t slice = fHoughTransformer->fSlice;
+  Int_t num_of_pads_to_look = 1;
+  Int_t rows_to_miss = 2;
+
+  for(Int_t i=0; i<tracks->GetNTracks(); i++)
+    {
+      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+      if(!track) {printf("No track\n"); return;}
+      
+      Int_t nrow = 0;
+      
+      for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
+       {
+         
+         Float_t xyz[3];
+         //Get the crossing point of the track and current padrow:
+         //if(!track->GetCrossingPoint(slice,padrow,xyz))
+         if(!track->GetCrossingPoint(padrow,xyz))  
+           {
+             printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
+             continue;
+           }
+         
+         if(fake)
+           fake->Fill(xyz[0],xyz[1],1);
+         Int_t npixs = 0;
+         
+         //Get the pixels along the track candidate
+         for(pixel=(AliL3Digits*)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
+           {
+             
+             Int_t sector,row;
+             Float_t xyz_pix[3];
+             fTransform->Slice2Sector(slice,padrow,sector,row);
+             fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
+             
+             
+             //check if we are inside road
+             if(fabs(xyz_pix[1] - xyz[1]) > num_of_pads_to_look*fTransform->GetPadPitchWidthLow()) continue; 
+             npixs++;
+             
+             if(hist)
+               {
+                 //fTransform->Local2Global(xyz_pix,slice);
+                 hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
+               }    
+             
+           }
+         if(npixs > 0)
+           {
+             nrow++;
+           }     
+       }
+      
+      if(nrow < NRows[patch][1]-NRows[patch][0]-rows_to_miss)
+       tracks->Remove(i);
+             
+    }
+  
+  tracks->Compress();
+
+  
+}
+
+void AliL3HoughEval::DisplaySlice(TH2F *hist)
+{
+
+  AliL3Digits *pixel;
+  
+  for(Int_t padrow = 0; padrow < 174; padrow++)
+    {
+      
+      for(pixel=(AliL3Digits*)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
+       {
+         
+         Int_t sector,row;
+         Float_t xyz_pix[3];
+         fTransform->Slice2Sector(fHoughTransformer->fSlice,padrow,sector,row);
+         fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
+         
+         if(hist)
+           {
+             hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
+           }    
+         
+       }
+    }
+  
+}
+
+
+void AliL3HoughEval::CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta)
+{
+  
+  Int_t slice = fSlice;
+  
+  TFile *file = new TFile(rootfile);
+  file->cd();
+
+  AliRun *gAlice = (AliRun*)file->Get("gAlice");
+  gAlice->GetEvent(0);
+  
+  TClonesArray *particles=gAlice->Particles();
+  Int_t n=particles->GetEntriesFast();
+  Float_t torad=TMath::Pi()/180;
+  Float_t phi_min = slice*20 - 10;
+  Float_t phi_max = slice*20 + 10;
+  
+  for (Int_t j=0; j<n; j++) {
+    TParticle *p=(TParticle*)particles->UncheckedAt(j);
+    if (p->GetFirstMother()>=0) continue;  //secondary particle
+    if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
+    Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py();//,pzg=p->Pz();
+    Double_t phi_part = TMath::ATan2(pyg,pxg);
+    if (phi_part < 0) phi_part += 2*TMath::Pi();
+    
+    if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
+    if(ptg<0.100) continue;
+    
+    AliL3HoughTrack *sel_track=0;
+    
+    Double_t min_dist = 10000;
+    
+    for(Int_t t=0; t<merged_tracks->GetNTracks(); t++)
+      {
+       AliL3HoughTrack *track = (AliL3HoughTrack*)merged_tracks->GetCheckedTrack(t);
+       if(!track) {printf("AliL3HoughEval: NO TRACK\n"); continue;}
+       Float_t phi0[2] = {track->GetPhi0(),0};
+       fTransform->Local2GlobalAngle(phi0,slice);
+       Double_t dPt = ptg - track->GetPt();
+       Double_t dPhi0 = phi_part - phi0[0];
+       Double_t distance = pow(dPt,2) + pow(dPhi0,2);
+       if(distance < min_dist)
+         {
+           min_dist = distance;
+           sel_track = track;
+         }      
+       
+      }
+    if(sel_track)
+      sel_track->SetBestMCid(j,min_dist);
+    
+    Double_t dpt = fabs(ptg-sel_track->GetPt());
+    Double_t dphi0 = fabs(phi_part-sel_track->GetPhi0());
+    //printf("Found match, min_dist %f dPt %f dPhi0 %f\n",min_dist,dpt,dphi0);
+  }
+  file->Close();
+  delete file;
+
+  
+}
diff --git a/HLT/hough/AliL3HoughEval.h b/HLT/hough/AliL3HoughEval.h
new file mode 100644 (file)
index 0000000..42ee129
--- /dev/null
@@ -0,0 +1,37 @@
+#ifndef ALIL3_HOUGH_Eval
+#define ALIL3_HOUGH_Eval
+
+#include "AliL3RootTypes.h"
+
+class TClonesArray;
+class TObjArray;
+class TH2F;
+class TH1F;
+class AliL3Transform;
+class AliL3TrackArray;
+class AliL3HoughTransformer;
+
+class AliL3HoughEval : public TObject {
+  
+ private:
+  
+  Int_t fSlice;
+  AliL3HoughTransformer *fHoughTransformer;
+  AliL3Transform *fTransform; //!
+
+ public:
+  AliL3HoughEval(); 
+  AliL3HoughEval(AliL3HoughTransformer *transformer);
+  virtual ~AliL3HoughEval();
+  
+  void LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist=0,TH2F *fake=0);
+  void DisplaySlice(TH2F *hist);
+  void CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta);
+  
+  void SetTransformer(AliL3HoughTransformer *t) {fHoughTransformer=t;}
+
+  ClassDef(AliL3HoughEval,1)
+
+};
+
+#endif
diff --git a/HLT/hough/AliL3HoughLinkDef.h b/HLT/hough/AliL3HoughLinkDef.h
new file mode 100644 (file)
index 0000000..534862b
--- /dev/null
@@ -0,0 +1,15 @@
+#ifdef __CINT__
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliL3HoughTransformer; 
+#pragma link C++ class AliL3HoughPixel; 
+//#pragma link C++ class AliL3HoughTrack; 
+#pragma link C++ class AliL3HoughMaxFinder;
+#pragma link C++ class AliL3HoughEval;
+#pragma link C++ class AliL3HoughMerge;
+
+#endif
+
diff --git a/HLT/hough/AliL3HoughMaxFinder.cxx b/HLT/hough/AliL3HoughMaxFinder.cxx
new file mode 100644 (file)
index 0000000..d36b28c
--- /dev/null
@@ -0,0 +1,113 @@
+#include <string.h>
+
+#include <TH2.h>
+
+#include "AliL3TrackArray.h"
+#include "AliL3HoughTrack.h"
+#include "AliL3HoughMaxFinder.h"
+
+ClassImp(AliL3HoughMaxFinder)
+
+  
+AliL3HoughMaxFinder::AliL3HoughMaxFinder()
+{
+  //Default constructor
+  fThreshold = 0;
+  //fTracks = 0;
+  fHistoType=0;
+}
+
+
+AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype)
+{
+  //Constructor
+
+  //fTracks = new AliL3TrackArray("AliL3HoughTrack");
+  if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
+  if(strcmp(histotype,"DPsi")==0) fHistoType='l';
+  
+}
+
+
+AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
+{
+  //Destructor
+
+}
+
+
+AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_t ref_row)
+{
+  //Locate all the maxima in input histogram.
+  //Maxima is defined as bins with more entries than the
+  //immediately neighbouring bins. 
+  
+  Int_t xmin = hist->GetXaxis()->GetFirst();
+  Int_t xmax = hist->GetXaxis()->GetLast();
+  Int_t ymin = hist->GetYaxis()->GetFirst();
+  Int_t ymax = hist->GetYaxis()->GetLast();
+  Int_t bin[9],track_counter=0;
+  Stat_t value[9];
+  
+  AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
+  AliL3HoughTrack *track;
+
+  for(Int_t xbin=xmin+1; xbin<xmax-1; xbin++)
+    {
+      for(Int_t ybin=ymin+1; ybin<ymax-1; ybin++)
+       {
+         bin[0] = hist->GetBin(xbin-1,ybin-1);
+         bin[1] = hist->GetBin(xbin,ybin-1);
+         bin[2] = hist->GetBin(xbin+1,ybin-1);
+         bin[3] = hist->GetBin(xbin-1,ybin);
+         bin[4] = hist->GetBin(xbin,ybin);
+         bin[5] = hist->GetBin(xbin+1,ybin);
+         bin[6] = hist->GetBin(xbin-1,ybin+1);
+         bin[7] = hist->GetBin(xbin,ybin+1);
+         bin[8] = hist->GetBin(xbin+1,ybin+1);
+         value[0] = hist->GetBinContent(bin[0]);
+         value[1] = hist->GetBinContent(bin[1]);
+         value[2] = hist->GetBinContent(bin[2]);
+         value[3] = hist->GetBinContent(bin[3]);
+         value[4] = hist->GetBinContent(bin[4]);
+         value[5] = hist->GetBinContent(bin[5]);
+         value[6] = hist->GetBinContent(bin[6]);
+         value[7] = hist->GetBinContent(bin[7]);
+         value[8] = hist->GetBinContent(bin[8]);
+         
+         if(value[4] <= fThreshold) continue;//central bin below threshold
+         
+         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 max_x = hist->GetXaxis()->GetBinCenter(xbin);
+             Float_t max_y = hist->GetYaxis()->GetBinCenter(ybin);
+             
+             track = (AliL3HoughTrack*)tracks->NextTrack();
+             
+             
+             switch(fHistoType)
+               {
+               case 'c':
+                 track->SetTrackParameters(max_x,max_y,(Int_t)value[4]);
+                 break;
+               case 'l':
+                 track->SetLineParameters(max_x,max_y,(Int_t)value[4],rowrange,ref_row);
+                 break;
+               default:
+                 printf("AliL3HoughMaxFinder: Error in tracktype\n");
+               }
+             
+             track_counter++;
+             
+
+           }
+         else
+           continue; //not a maxima
+       }
+    }
+  tracks->QSort();
+  return tracks;
+}
diff --git a/HLT/hough/AliL3HoughMaxFinder.h b/HLT/hough/AliL3HoughMaxFinder.h
new file mode 100644 (file)
index 0000000..e937dda
--- /dev/null
@@ -0,0 +1,31 @@
+#ifndef ALIL3_HOUGH_MaxFinder
+#define ALIL3_HOUGH_MaxFinder
+
+#include "AliL3RootTypes.h"
+
+class AliL3TrackArray;
+
+class AliL3HoughMaxFinder : public TObject {
+  
+ private:
+
+  Int_t fThreshold;
+  //AliL3TrackArray *fTracks; //!
+
+  Char_t fHistoType;
+
+ public:
+  AliL3HoughMaxFinder(); 
+  AliL3HoughMaxFinder(Char_t *histotype);
+  virtual ~AliL3HoughMaxFinder();
+
+  AliL3TrackArray *FindMaxima(TH2F *hist,Int_t *rowrange=0,Int_t ref_row=0);
+  void SetThreshold(Int_t f) {fThreshold = f;}
+
+  //AliL3TrackArray *GetTracks() {return fTracks;}
+
+  ClassDef(AliL3HoughMaxFinder,1)
+
+};
+
+#endif
diff --git a/HLT/hough/AliL3HoughMerge.cxx b/HLT/hough/AliL3HoughMerge.cxx
new file mode 100644 (file)
index 0000000..534e329
--- /dev/null
@@ -0,0 +1,71 @@
+#include <TH2.h>
+
+#include "AliL3TrackArray.h"
+#include "AliL3HoughTrack.h"
+#include "AliL3HoughMerge.h"
+#include "AliL3HoughTransformer.h"
+
+ClassImp(AliL3HoughMerge)
+
+  
+AliL3HoughMerge::AliL3HoughMerge()
+{
+  //Default constructor
+
+  fInTracks = NULL;
+  fOutTracks = NULL;
+  fNPatches = 5;
+}
+
+
+AliL3HoughMerge::AliL3HoughMerge(Int_t slice,Int_t row_patches)
+{
+  //Constructor
+
+  //fInTracks = (AliL3TrackArray**)new Byte_t[row_patches*sizeof(AliL3TrackArray*)];
+  fInTracks = new AliL3TrackArray*[row_patches];
+  fNPatches = row_patches;
+  for(Int_t i=0; i<row_patches; i++)
+    fInTracks[i] = new AliL3TrackArray("AliL3HoughTrack");
+
+  fOutTracks = new AliL3TrackArray("AliL3HoughTrack");
+}
+
+
+AliL3HoughMerge::~AliL3HoughMerge()
+{
+  //Destructor
+  if(fInTracks)
+    delete fInTracks;
+  if(fOutTracks)
+    delete fOutTracks;
+  
+}
+
+void AliL3HoughMerge::FillTracks(AliL3TrackArray *tracks,Int_t patch)
+{
+  fInTracks[patch]->AddTracks(tracks); //copies tracks to new trackarray. Does not delete the track objects.
+}
+
+void AliL3HoughMerge::FillHisto(TH2F *merge_hist)
+{
+  
+  for(Int_t pat=0; pat < fNPatches; pat++)
+    {
+      for(Int_t t=0; t<fInTracks[pat]->GetNTracks(); t++)
+       {
+         AliL3HoughTrack *tr = (AliL3HoughTrack*)fInTracks[pat]->GetCheckedTrack(t);
+         if(!tr) {printf("AliL3HoughMerge NO TRACK\n"); continue;}
+         merge_hist->Fill(tr->GetKappa(),tr->GetPhi0(),tr->GetNHits());
+       }
+    }
+}
+
+/*
+void AliL3HoughMerge::MergeLines()
+{
+
+  
+
+}
+*/
diff --git a/HLT/hough/AliL3HoughMerge.h b/HLT/hough/AliL3HoughMerge.h
new file mode 100644 (file)
index 0000000..cf3653b
--- /dev/null
@@ -0,0 +1,31 @@
+#ifndef ALIL3_HOUGH_Merge
+#define ALIL3_HOUGH_Merge
+
+#include "AliL3RootTypes.h"
+
+class TClonesArray;
+class TObjArray;
+class TH2F;
+class AliL3TrackArray;
+
+class AliL3HoughMerge : public TObject {
+  
+ private:
+
+  AliL3TrackArray **fInTracks; //!
+  AliL3TrackArray *fOutTracks; //!
+  
+  Int_t fNPatches;
+ public:
+  AliL3HoughMerge(); 
+  AliL3HoughMerge(Int_t slice,Int_t row_patches=5);
+  virtual ~AliL3HoughMerge();
+
+  void FillTracks(AliL3TrackArray *tracks,Int_t patch);
+  void FillHisto(TH2F *merge_hist);
+  
+  ClassDef(AliL3HoughMerge,1)
+
+};
+
+#endif
diff --git a/HLT/hough/AliL3HoughPixel.cxx b/HLT/hough/AliL3HoughPixel.cxx
new file mode 100644 (file)
index 0000000..f7c89a8
--- /dev/null
@@ -0,0 +1,45 @@
+
+#include "AliL3Transform.h"
+#include "AliL3HoughPixel.h"
+
+ClassImp(AliL3HoughPixel)
+
+AliL3HoughPixel::AliL3HoughPixel()
+{
+
+}
+
+AliL3HoughPixel::AliL3HoughPixel(Int_t pad,Int_t time,Short_t signal)
+  //                            Int_t slice,Int_t padrow)
+{
+  //Constructor
+  fPad = pad;
+  fTime = time;
+  fSignal = signal;
+  /*fSlice = slice;
+  fPadrow = padrow;
+  */
+}
+
+AliL3HoughPixel::~AliL3HoughPixel()
+{
+
+}
+/*
+void AliL3HoughPixel::SetCoordinates()
+{
+
+  AliL3Transform *transform = new AliL3Transform();
+  
+  Float_t xyz[3];
+  Int_t sector,row;
+  transform->Slice2Sector(fSlice,fPadrow,sector,row);
+  transform->Raw2Global(xyz,sector,row,fPad,fTime);
+  fX = xyz[0];
+  fY = xyz[1];
+  fZ = xyz[2];
+  fPhi = transform->GetPhi(xyz);
+  fEta = transform->GetEta(xyz);
+  
+}
+*/
diff --git a/HLT/hough/AliL3HoughPixel.h b/HLT/hough/AliL3HoughPixel.h
new file mode 100644 (file)
index 0000000..734ac77
--- /dev/null
@@ -0,0 +1,45 @@
+#ifndef ALI_PIXEL
+#define ALI_PIXEL
+
+#include <TObject.h>
+
+class AliL3HoughPixel : public TObject {
+
+ private: 
+  Short_t fSignal;
+  Int_t fPad;
+  Int_t fTime;
+  /*Int_t fSlice;
+  Int_t fPadrow;
+
+  Float_t fX;
+  Float_t fY;
+  Float_t fZ;
+  Double_t fPhi;
+  Double_t fEta;
+  */
+ public:
+  
+  AliL3HoughPixel();
+  //AliL3HoughPixel(Int_t pad,Int_t time,Short_t signal,Int_t slice,Int_t padrow);
+  AliL3HoughPixel(Int_t pad,Int_t time,Short_t signal);
+  virtual ~AliL3HoughPixel();
+
+  /*Int_t GetPad() {return fPad;}
+  Int_t GetTime() {return fTime;}
+  Short_t GetSignal() {return fSignal;}
+  Double_t GetPhi() {return fPhi;}
+  Double_t GetEta() {return fEta;}
+  Float_t GetX() {return fX;}
+  Float_t GetY() {return fY;}
+  Float_t GetZ() {return fZ;}
+
+  void SetCoordinates();
+  */
+  AliL3HoughPixel *nextRowPixel;
+
+  ClassDef(AliL3HoughPixel,1)
+};
+
+
+#endif
diff --git a/HLT/hough/AliL3HoughTransformer.cxx b/HLT/hough/AliL3HoughTransformer.cxx
new file mode 100644 (file)
index 0000000..5dfeccd
--- /dev/null
@@ -0,0 +1,580 @@
+#include <math.h>
+
+#include <TFile.h>
+#include <TTree.h>
+#include <TH2.h>
+
+#include "AliTPCParam.h"
+#include "AliSimDigits.h"
+
+#include "AliL3Defs.h"
+#include "AliL3Transform.h"
+#include "AliL3HoughPixel.h"
+#include "AliL3HoughTransformer.h"
+#include "AliL3HoughTrack.h"
+#include "AliL3TrackArray.h"
+
+ClassImp(AliL3HoughTransformer)
+
+AliL3HoughTransformer::AliL3HoughTransformer()
+{
+  //Default constructor
+  fTransform = 0;
+  fEtaMin = 0;
+  fEtaMax = 0;
+  fSlice = 0;
+  fPatch = 0;
+  fRowContainer = 0;
+  fPhiRowContainer = 0;
+  fNumOfPadRows=0;
+  fContainerBounds=0;
+  fNDigits=0;
+  fIndex = 0;
+}
+
+
+AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments)
+{
+  //Constructor
+  
+  fTransform = new AliL3Transform();
+  
+
+  fEtaMin = etarange[0];
+  fEtaMax = etarange[1];
+  fSlice = slice; 
+  fPatch = patch;
+  fNPhiSegments = phi_segments;
+  fNumOfPadRows=NRowsSlice;
+}
+
+
+AliL3HoughTransformer::~AliL3HoughTransformer()
+{
+  //Destructor
+  if(fRowContainer)
+    delete [] fRowContainer;
+  if(fTransform)
+    delete fTransform;
+  if(fPhiRowContainer)
+    delete [] fPhiRowContainer;
+  if(fIndex)
+    {
+      for(Int_t i=0; i<fNDigits; i++)
+       delete [] fIndex[i];
+      delete [] fIndex;
+    }
+}
+
+void AliL3HoughTransformer::InitTemplates(TH2F *hist)
+{
+
+  AliL3Digits *pixel;
+
+  Int_t ymin = hist->GetYaxis()->GetFirst();
+  Int_t ymax = hist->GetYaxis()->GetLast();
+  Int_t nbinsy = hist->GetNbinsY();
+
+  fIndex = new Int_t*[fNDigits];
+  for(Int_t i=0; i<fNDigits; i++)
+    fIndex[i] = new Int_t[nbinsy+1];
+    
+  Int_t sector,row;
+  for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+    {
+      
+      for(pixel=(AliL3Digits*)fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
+       {
+         Float_t xyz[3];
+         fTransform->Slice2Sector(fSlice,padrow,sector,row);
+         fTransform->Raw2Local(xyz,sector,row,pixel->fPad,pixel->fTime);
+         
+         Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+         
+         Double_t phi_pix = fTransform->GetPhi(xyz);
+         Short_t signal = pixel->fCharge;
+         Int_t index = pixel->fIndex;
+         if(index >= fNDigits)
+           printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index);
+         for(Int_t p=ymin; p<=ymax; p++)
+           {
+             
+             Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
+             Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
+             
+             Int_t bin = hist->FindBin(kappa,phi0);
+             if(fIndex[index][p]!=0)
+               printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n");
+             fIndex[index][p] = bin;
+           }
+       }
+    }
+  
+}
+
+
+void AliL3HoughTransformer::CountBins()
+{
+  
+  Int_t middle_row = 87; //middle of the slice
+  
+  Double_t r_in_bundle = fTransform->Row2X(middle_row);
+  //  Double_t phi_min = (fSlice*20 - 10)*ToRad;
+  //Double_t phi_max = (fSlice*20 + 10)*ToRad;
+  Double_t phi_min = -15*ToRad;
+  Double_t phi_max = 15*ToRad;
+
+  Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
+  Double_t min_phi0 = 10000;
+  Double_t max_phi0 = 0;
+  Double_t min_kappa = 100000;
+  Double_t max_kappa = 0;
+  
+  Int_t xbin = 60;
+  Int_t ybin = 60;
+  Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2->
+  Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
+  
+  TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  
+  for(Int_t padrow=NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+    {
+      for(Int_t pad=0; pad < fTransform->GetNPads(padrow); pad++)
+       {
+         for(Int_t time = 0; time < fTransform->GetNTimeBins(); time++)
+           {
+             Float_t xyz[3];
+             Int_t sector,row;
+             fTransform->Slice2Sector(fSlice,padrow,sector,row);
+             fTransform->Raw2Global(xyz,sector,row,pad,time);
+             Double_t eta = fTransform->GetEta(xyz);
+             if(eta < fEtaMin || eta > fEtaMax) continue;
+             fTransform->Global2Local(xyz,sector);
+             Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+             Double_t phi_pix = fTransform->GetPhi(xyz);
+             
+             for(Int_t p=0; p<=fNPhiSegments; p++)
+               {
+                 Double_t phi_in_bundle = phi_min + p*phi_slice;
+                 
+                 Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
+                 
+                 Double_t phi0 = atan(tanPhi0);
+                 //  if(phi0 < 0.55 || phi0 > 0.88) continue;
+                 
+                 //if(phi0 < 0) phi0 = phi0 +2*Pi;
+                 //Double_t kappa = sin(phi_in_bundle - phi0)*2/r_in_bundle;
+                 
+                 Double_t angle = phi_pix - phi0;
+                 Double_t kappa = 2*sin(angle)/r_pix;
+                 histo->Fill(kappa,phi0,1);
+                 if(phi0 < min_phi0)
+                   min_phi0 = phi0;
+                 if(phi0 > max_phi0)
+                   max_phi0 = phi0;
+                 if(kappa < min_kappa)
+                   min_kappa = kappa;
+                 if(kappa > max_kappa)
+                   max_kappa = kappa;
+                                 
+               }
+             
+           }
+         
+       }
+      
+    }
+  Int_t count=0,bi=0;
+    
+  Int_t xmin = histo->GetXaxis()->GetFirst();
+  Int_t xmax = histo->GetXaxis()->GetLast();
+  Int_t ymin = histo->GetYaxis()->GetFirst();
+  Int_t ymax = histo->GetYaxis()->GetLast();
+
+
+  for(Int_t xbin=xmin+1; xbin<xmax; xbin++)
+    {
+      for(Int_t ybin=ymin+1; ybin<ymax; ybin++)
+       {
+         bi++;
+         Int_t bin = histo->GetBin(xbin,ybin);
+         if(histo->GetBinContent(bin)>0)
+           count++;
+       }
+    }
+
+
+  printf("Number of possible tracks in this region %d, bins %d\n",count,bi);
+  printf("Phi, min %f max %f\n",min_phi0,max_phi0);
+  printf("Kappa, min %f max %f\n",min_kappa,max_kappa);
+  histo->Draw("box");
+}
+
+
+void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t middle_row)
+{
+  //Transformation is done with respect to local coordinates in slice.
+
+
+  AliL3Digits *pix1;
+  Int_t sector,row;
+  
+  //Define a common point
+  /*
+    Double_t rowdist1 = fTransform->Row2X(middle_row-1);
+    Double_t rowdist2 = fTransform->Row2X(middle_row);
+    Double_t r_in_bundle = rowdist1 + (rowdist1-rowdist2)/2;
+  */
+  
+  Double_t r_in_bundle = fTransform->Row2X(middle_row);
+  //Make overlap between slices
+  Double_t phi_min = -15*ToRad;
+  Double_t phi_max = 15*ToRad;
+  
+  Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
+  
+  for(Int_t p=0; p <= fNPhiSegments; p++)
+    {
+      Double_t phi_in_bundle = phi_min + p*phi_slice;
+      //printf("phi %f in slice %d patch %d middle row %f\n",phi_in_bundle/ToRad,fSlice,fPatch,r_in_bundle);
+      
+      for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+       //for(Int_t padrow = middle_row; padrow <= 173; padrow++)
+       //for(Int_t padrow = 0; padrow <= middle_row; padrow++)
+       {
+         
+         for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+           {
+             
+             Float_t xyz[3];
+             fTransform->Slice2Sector(fSlice,padrow,sector,row);
+             fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
+             //fTransform->Raw2Global(xyz,sector,row,pix1->fPad,pix1->fTime);
+             
+             Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+             
+             Double_t phi_pix = fTransform->GetPhi(xyz);
+                             
+             //Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
+             Double_t tanPhi0 = (r_in_bundle*sin(phi_pix)-r_pix*sin(phi_in_bundle))/(r_in_bundle*cos(phi_pix)-r_pix*cos(phi_in_bundle));
+             
+             Double_t phi0 = atan(tanPhi0);
+             //if(padrow > middle_row)
+             //phi0 = -phi0;
+             //if(phi0 < 0.55 || phi0 > 0.88) continue;
+             
+             Double_t angle = phi_pix - phi0;
+             Double_t kappa = 2*sin(angle)/r_pix;
+             
+             //Double_t angle = phi_in_bundle - phi0;
+             //Double_t kappa = 2*sin(angle)/r_in_bundle;
+
+             //if(kappa < -0.006 || kappa > 0.006) continue;
+             
+             Short_t signal = pix1->fCharge;
+             
+             hist->Fill(kappa,phi0,signal);
+             
+             
+           }
+         
+       }
+    }
+  
+}
+
+
+void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
+{
+  //Transformation is done with respect to local coordinates in slice.
+  //Transform every pixel into whole phirange, using parametrisation:
+  //kappa = 2*sin(phi-phi0)/R
+
+  printf("Transforming 1 pixel only\n");
+
+  AliL3Digits *pix1;
+  Int_t sector,row;
+  
+  Int_t ymin = hist->GetYaxis()->GetFirst();
+  Int_t ymax = hist->GetYaxis()->GetLast();
+  Int_t nbinsy = hist->GetNbinsY();
+
+  for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+    {
+      
+      for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+       {
+         
+         Short_t signal = pix1->fCharge;
+         Int_t index = pix1->fIndex;
+         
+         for(Int_t p=0; p <= nbinsy; p++)
+           hist->AddBinContent(fIndex[index][p],signal);
+                         
+       }
+      
+    }
+    
+  
+}
+
+/*
+  void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
+  {
+  //Transformation is done with respect to local coordinates in slice.
+  //Transform every pixel into whole phirange, using parametrisation:
+  //kappa = 2*sin(phi-phi0)/R
+  
+  printf("Transforming 1 pixel only\n");
+  
+  AliL3Digits *pix1;
+  Int_t sector,row;
+  
+  Int_t ymin = hist->GetYaxis()->GetFirst();
+  Int_t ymax = hist->GetYaxis()->GetLast();
+  
+  for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+  {
+  
+  for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+  {
+  
+  Float_t xyz[3];
+  fTransform->Slice2Sector(fSlice,padrow,sector,row);
+  fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
+  
+  Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+  
+  Double_t phi_pix = fTransform->GetPhi(xyz);
+  Short_t signal = pix1->fCharge;
+  
+  for(Int_t p=ymin+1; p<=ymax; p++)
+  {
+  Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
+  Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
+  hist->Fill(kappa,phi0,signal);
+  }
+  
+  }
+  
+  }
+  
+  
+  }
+*/
+
+void AliL3HoughTransformer::TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks)
+{
+
+  for(Int_t i=0; i<tracks->GetNTracks(); i++)
+    {
+      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+      if(!track) {printf("AliL3HoughTransformer::TransformLines2Circle : NO TRACK IN ARRAY\n"); continue;}
+     
+      Double_t xmin = fTransform->Row2X(track->GetFirstRow());
+      Double_t xmax = fTransform->Row2X(track->GetLastRow());
+      
+      Double_t a = -1./tan(track->GetPsiLine());
+      Double_t b = track->GetDLine()/sin(track->GetPsiLine());
+      
+      Double_t ymin = a*xmin + b;
+      Double_t ymax = a*xmax + b;
+      
+      Double_t middle_x = xmin + (xmax-xmin)/2;
+      Double_t middle_y = ymin + (ymax-ymin)/2;
+      
+      Double_t r_middle = sqrt(middle_x*middle_x + middle_y*middle_y);
+      Double_t phi = atan2(middle_y,middle_x);
+      Double_t phi0 = 2*phi - track->GetPsiLine();
+      Double_t kappa = 2*sin(phi-phi0)/r_middle;
+      hist->Fill(kappa,phi0,track->GetWeight());
+     
+    }
+
+  
+}
+
+void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw)
+{
+  //Transform every pixel into histogram, using parametrisation:
+  //D = x*cos(psi) + y*sin(psi)
+
+  //  printf("In Transform; rowrange %d %d ref_row %d phirange %f %f\n",rowrange[0],rowrange[1],ref_row,phirange[0],phirange[1]);
+
+  AliL3Digits *pix1;
+  
+  Int_t xmin = hist->GetXaxis()->GetFirst();
+  Int_t xmax = hist->GetXaxis()->GetLast();
+  
+  Double_t x0 = fTransform->Row2X(ref_row);
+  Double_t y0 = 0;
+
+  Int_t sector,row;
+  //for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
+  
+  Double_t phi_min = -10*ToRad;
+  Double_t phi_max = 10*ToRad;
+  Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
+  
+
+  Int_t phi_min_index = (Int_t)((phirange[0]*ToRad-phi_min)/delta_phi);
+  Int_t phi_max_index = (Int_t)((phirange[1]*ToRad-phi_min)/delta_phi);
+    
+  
+  Int_t index;
+  
+  for(Int_t phi=phi_min_index; phi <= phi_max_index; phi++)
+    {
+      for(Int_t padrow = rowrange[0]; padrow <= rowrange[1]; padrow++)
+       {
+         
+         index = phi*fNumOfPadRows + padrow;
+         //printf("Looping index %d\n",index);
+         if(index > fContainerBounds || index < 0)
+           {
+             printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index);
+             return;
+           }
+         //for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+         for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel)
+           {
+             //printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
+             Float_t xyz[3];
+             fTransform->Slice2Sector(fSlice,padrow,sector,row);
+             fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
+                 
+             if(raw)
+               raw->Fill(xyz[0],xyz[1],pix1->fCharge);
+
+             xyz[0] = xyz[0]-x0;
+             xyz[1] = xyz[1]-y0;
+             
+             //printf("calculating...");
+             for(Int_t d=xmin+1; d<=xmax; d++)
+               {
+                 Double_t psi = hist->GetXaxis()->GetBinCenter(d);
+                 Double_t D = xyz[0]*cos(psi) + xyz[1]*sin(psi);
+                 
+                 Short_t signal = pix1->fCharge;
+                 hist->Fill(psi,D,signal);
+                 //printf("Filling histo, psi %f D %f\n",psi,D);
+               }
+             //printf("done\n");
+           }
+         //printf(" done\n");
+       }
+      
+    }
+    
+}
+
+void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
+{
+
+  TFile *file = new TFile(rootfile);
+  file->cd();
+
+  AliTPCParam *param=(AliTPCParam *)file->Get("75x40_100x60");
+  TTree *t=(TTree*)file->Get("TreeD_75x40_100x60");      
+    
+  AliSimDigits da, *digarr=&da;
+  t->GetBranch("Segment")->SetAddress(&digarr);
+  Stat_t num_of_entries=t->GetEntries();
+
+  Int_t digit_counter=0;
+  Float_t xyz[3];
+  Double_t eta,phi;
+  
+  Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
+  printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
+  
+  if(fRowContainer)
+    delete [] fRowContainer;
+  fRowContainer = new AliL3HoughContainer[fNumOfPadRows];
+  memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer));
+
+
+  fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
+  printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
+  if(fPhiRowContainer)
+    delete [] fPhiRowContainer;
+  fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
+  memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
+
+  Double_t phi_min = -10*ToRad;
+  Double_t phi_max = 10*ToRad;
+  Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
+  Int_t index;
+  digit_counter=0;
+
+  for (Int_t i=0; i<num_of_entries; i++) 
+    { 
+      t->GetEvent(i);
+      Int_t sector; 
+      Int_t row;    
+      param->AdjustSectorRow(digarr->GetID(),sector,row);
+      Int_t slice,padrow;
+      fTransform->Sector2Slice(slice,padrow,sector,row);
+      if(slice != fSlice) continue;
+      if(padrow < NRows[fPatch][0]) continue;
+      if(padrow > NRows[fPatch][1]) break;
+      digarr->First();
+      do {
+       Int_t time=digarr->CurrentRow();
+       Int_t pad=digarr->CurrentColumn();
+       Short_t signal=digarr->CurrentDigit();
+       if(time < param->GetMaxTBin()-1 && time > 0)
+         if(digarr->GetDigit(time-1,pad) <= param->GetZeroSup()
+            && digarr->GetDigit(time+1,pad) <= param->GetZeroSup()) 
+           continue;
+       
+       
+       fTransform->Raw2Global(xyz,sector,row,pad,time);
+       eta = fTransform->GetEta(xyz);
+       if(eta < fEtaMin || eta > fEtaMax) continue;
+       fTransform->Global2Local(xyz,sector);
+       
+       
+       phi = fTransform->GetPhi(xyz);
+       if(hist)
+         hist->Fill(xyz[0],xyz[1],signal);
+       
+       AliL3Digits *dig = new AliL3Digits;
+       dig->fIndex = digit_counter;
+       digit_counter++;
+       dig->fCharge = signal;
+       dig->fPad = pad;
+       dig->fTime = time;
+       
+       if(fRowContainer[padrow].first == NULL)
+         fRowContainer[padrow].first = (void*)dig;
+       else
+         ((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
+       fRowContainer[padrow].last = (void*)dig;
+       
+       Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
+       index = phi_index*fNumOfPadRows + padrow;
+       if(phi_index > fContainerBounds || phi_index < 0)
+         {
+           printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
+           continue;
+         }
+               
+       if(fPhiRowContainer[index].first == NULL)
+         fPhiRowContainer[index].first = (void*)dig;
+       else
+         ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
+       fPhiRowContainer[index].last=(void*)dig;
+       
+      }while (digarr->Next());
+      
+    }
+  
+  fNDigits = digit_counter;
+  printf("digitcounter %d\n",digit_counter);
+  printf("Allocated %d bytes to pixels\n",digit_counter*sizeof(AliL3Digits));
+  file->Close();
+  delete file;
+  
+}
diff --git a/HLT/hough/AliL3HoughTransformer.h b/HLT/hough/AliL3HoughTransformer.h
new file mode 100644 (file)
index 0000000..2833e6a
--- /dev/null
@@ -0,0 +1,68 @@
+#ifndef ALIL3_HOUGHTRANSFORMER
+#define ALIL3_HOUGHTRANSFORMER
+
+#include "AliL3RootTypes.h"
+
+struct AliL3Digits
+{
+  UShort_t fCharge;
+  UChar_t fPad;
+  UShort_t fTime;
+  Int_t fIndex;
+  AliL3Digits *nextPhiRowPixel;
+  AliL3Digits *nextRowPixel;
+};
+
+struct AliL3HoughContainer
+{
+  void *first;
+  void *last;
+};
+
+class TH2F;
+class AliL3Transform;
+class AliL3HoughEvaluate;
+class AliL3TrackArray;
+
+class AliL3HoughTransformer : public TObject {
+  
+ private:
+
+  friend class AliL3HoughEval;
+  
+  AliL3Transform *fTransform; //!
+
+  Int_t fNPhiSegments; //Number of patches in phi.
+  Float_t fEtaMin;
+  Float_t fEtaMax;
+  Int_t fNumEtaSegments;
+  Int_t fNumOfPadRows;
+  
+  AliL3HoughContainer *fRowContainer; //!
+  AliL3HoughContainer *fPhiRowContainer; //!
+  Int_t fContainerBounds;
+  Int_t fNDigits;
+  Int_t **fIndex; //!
+  
+  
+  Int_t fSlice;
+  Int_t fPatch;
+
+ public:
+  AliL3HoughTransformer(); 
+  AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments);
+  virtual ~AliL3HoughTransformer();
+  
+  void Transform2Circle(TH2F *hist,Int_t middle_row);
+  void Transform2Circle(TH2F *hist);
+  void Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw=0);
+  void TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks);
+  void GetPixels(Char_t *rootfile,TH2F *hist=0);
+  void InitTemplates(TH2F *hist);
+  void CountBins();
+
+  ClassDef(AliL3HoughTransformer,1)
+
+};
+
+#endif
diff --git a/HLT/hough/AliL3Transform.cxx b/HLT/hough/AliL3Transform.cxx
new file mode 100644 (file)
index 0000000..cdcd6ba
--- /dev/null
@@ -0,0 +1,574 @@
+
+//Author:        Anders Strand Vestbo
+//Author:        Uli Frankenfeld
+//Last Modified: 05.12.2000
+
+#include "AliL3Logging.h"
+#include "AliL3Transform.h"
+//#include <TFile.h>
+#include <math.h>
+//___________________________
+// AliL3Transform
+//
+// Transformation class for ALICE TPC.
+
+ClassImp(AliL3Transform);
+
+
+AliL3Transform::AliL3Transform(){
+  //constructor
+  Init();
+}
+
+
+AliL3Transform::~AliL3Transform(){
+}
+
+void AliL3Transform::Init(){
+  //sector:
+  fNTimeBins = 512;  //?uli
+  fNRowLow = 55;
+  fNRowUp = 119;
+  fNSectorLow = 36;
+  fNSectorUp = 36;
+  fNSector = 72;
+  fPadPitchWidthLow = 0.400000;
+  fPadPitchWidthUp = 0.600000;
+  fZWidth = 0.56599998474121093750;
+  fZSigma = 0.22880849748219134199;
+
+  //slices:
+  fNSlice = 36;
+  fNRow = 174;
+  fPi = 3.14159265358979323846;
+  for(Int_t i=0;i<36;i++){
+    fCos[i] = cos(fPi/9*i);
+    fSin[i] = sin(fPi/9*i);
+  }
+
+  fX[0] = 88.850006103515625;
+  fX[1] = 89.600006103515625;
+  fX[2] = 90.350006103515625;
+  fX[3] = 91.100006103515625;
+  fX[4] = 91.850006103515625;
+  fX[5] = 92.600006103515625;
+  fX[6] = 93.350006103515625;
+  fX[7] = 94.100006103515625;
+  fX[8] = 94.850006103515625;
+  fX[9] = 95.600006103515625;
+  fX[10] = 96.350006103515625;
+  fX[11] = 97.100006103515625;
+  fX[12] = 97.850006103515625;
+  fX[13] = 98.600006103515625;
+  fX[14] = 99.350006103515625;
+  fX[15] = 100.100006103515625;
+  fX[16] = 100.850006103515625;
+  fX[17] = 101.600006103515625;
+  fX[18] = 102.350006103515625;
+  fX[19] = 103.100006103515625;
+  fX[20] = 103.850006103515625;
+  fX[21] = 104.600006103515625;
+  fX[22] = 105.350006103515625;
+  fX[23] = 106.100006103515625;
+  fX[24] = 106.850006103515625;
+  fX[25] = 107.600006103515625;
+  fX[26] = 108.350006103515625;
+  fX[27] = 109.100006103515625;
+  fX[28] = 109.850006103515625;
+  fX[29] = 110.600006103515625;
+  fX[30] = 111.350006103515625;
+  fX[31] = 112.100006103515625;
+  fX[32] = 112.850006103515625;
+  fX[33] = 113.600006103515625;
+  fX[34] = 114.350006103515625;
+  fX[35] = 115.100006103515625;
+  fX[36] = 115.850006103515625;
+  fX[37] = 116.600006103515625;
+  fX[38] = 117.350006103515625;
+  fX[39] = 118.100006103515625;
+  fX[40] = 118.850006103515625;
+  fX[41] = 119.600006103515625;
+  fX[42] = 120.350006103515625;
+  fX[43] = 121.100006103515625;
+  fX[44] = 121.850006103515625;
+  fX[45] = 122.600006103515625;
+  fX[46] = 123.350006103515625;
+  fX[47] = 124.100006103515625;
+  fX[48] = 124.850006103515625;
+  fX[49] = 125.600006103515625;
+  fX[50] = 126.350006103515625;
+  fX[51] = 127.100006103515625;
+  fX[52] = 127.850006103515625;
+  fX[53] = 128.600006103515625;
+  fX[54] = 129.350006103515625;
+  fX[55] = 132.574996948242188;
+  fX[56] = 133.574996948242188;
+  fX[57] = 134.574996948242188;
+  fX[58] = 135.574996948242188;
+  fX[59] = 136.574996948242188;
+  fX[60] = 137.574996948242188;
+  fX[61] = 138.574996948242188;
+  fX[62] = 139.574996948242188;
+  fX[63] = 140.574996948242188;
+  fX[64] = 141.574996948242188;
+  fX[65] = 142.574996948242188;
+  fX[66] = 143.574996948242188;
+  fX[67] = 144.574996948242188;
+  fX[68] = 145.574996948242188;
+  fX[69] = 146.574996948242188;
+  fX[70] = 147.574996948242188;
+  fX[71] = 148.574996948242188;
+  fX[72] = 149.574996948242188;
+  fX[73] = 150.574996948242188;
+  fX[74] = 151.574996948242188;
+  fX[75] = 152.574996948242188;
+  fX[76] = 153.574996948242188;
+  fX[77] = 154.574996948242188;
+  fX[78] = 155.574996948242188;
+  fX[79] = 156.574996948242188;
+  fX[80] = 157.574996948242188;
+  fX[81] = 158.574996948242188;
+  fX[82] = 159.574996948242188;
+  fX[83] = 160.574996948242188;
+  fX[84] = 161.574996948242188;
+  fX[85] = 162.574996948242188;
+  fX[86] = 163.574996948242188;
+  fX[87] = 164.574996948242188;
+  fX[88] = 165.574996948242188;
+  fX[89] = 166.574996948242188;
+  fX[90] = 167.574996948242188;
+  fX[91] = 168.574996948242188;
+  fX[92] = 169.574996948242188;
+  fX[93] = 170.574996948242188;
+  fX[94] = 171.574996948242188;
+  fX[95] = 172.574996948242188;
+  fX[96] = 173.574996948242188;
+  fX[97] = 174.574996948242188;
+  fX[98] = 175.574996948242188;
+  fX[99] = 176.574996948242188;
+  fX[100] = 177.574996948242188;
+  fX[101] = 178.574996948242188;
+  fX[102] = 179.574996948242188;
+  fX[103] = 180.574996948242188;
+  fX[104] = 181.574996948242188;
+  fX[105] = 182.574996948242188;
+  fX[106] = 183.574996948242188;
+  fX[107] = 184.574996948242188;
+  fX[108] = 185.574996948242188;
+  fX[109] = 186.574996948242188;
+  fX[110] = 187.574996948242188;
+  fX[111] = 188.574996948242188;
+  fX[112] = 189.574996948242188;
+  fX[113] = 190.574996948242188;
+  fX[114] = 191.574996948242188;
+  fX[115] = 192.574996948242188;
+  fX[116] = 193.574996948242188;
+  fX[117] = 194.574996948242188;
+  fX[118] = 195.574996948242188;
+  fX[119] = 196.574996948242188;
+  fX[120] = 197.574996948242188;
+  fX[121] = 198.574996948242188;
+  fX[122] = 199.574996948242188;
+  fX[123] = 200.574996948242188;
+  fX[124] = 201.574996948242188;
+  fX[125] = 202.574996948242188;
+  fX[126] = 203.574996948242188;
+  fX[127] = 204.574996948242188;
+  fX[128] = 205.574996948242188;
+  fX[129] = 206.574996948242188;
+  fX[130] = 207.574996948242188;
+  fX[131] = 208.574996948242188;
+  fX[132] = 209.574996948242188;
+  fX[133] = 210.574996948242188;
+  fX[134] = 211.574996948242188;
+  fX[135] = 212.574996948242188;
+  fX[136] = 213.574996948242188;
+  fX[137] = 214.574996948242188;
+  fX[138] = 215.574996948242188;
+  fX[139] = 216.574996948242188;
+  fX[140] = 217.574996948242188;
+  fX[141] = 218.574996948242188;
+  fX[142] = 219.574996948242188;
+  fX[143] = 220.574996948242188;
+  fX[144] = 221.574996948242188;
+  fX[145] = 222.574996948242188;
+  fX[146] = 223.574996948242188;
+  fX[147] = 224.574996948242188;
+  fX[148] = 225.574996948242188;
+  fX[149] = 226.574996948242188;
+  fX[150] = 227.574996948242188;
+  fX[151] = 228.574996948242188;
+  fX[152] = 229.574996948242188;
+  fX[153] = 230.574996948242188;
+  fX[154] = 231.574996948242188;
+  fX[155] = 232.574996948242188;
+  fX[156] = 233.574996948242188;
+  fX[157] = 234.574996948242188;
+  fX[158] = 235.574996948242188;
+  fX[159] = 236.574996948242188;
+  fX[160] = 237.574996948242188;
+  fX[161] = 238.574996948242188;
+  fX[162] = 239.574996948242188;
+  fX[163] = 240.574996948242188;
+  fX[164] = 241.574996948242188;
+  fX[165] = 242.574996948242188;
+  fX[166] = 243.574996948242188;
+  fX[167] = 244.574996948242188;
+  fX[168] = 245.574996948242188;
+  fX[169] = 246.574996948242188;
+  fX[170] = 247.574996948242188;
+  fX[171] = 248.574996948242188;
+  fX[172] = 249.574996948242188;
+  fX[173] = 250.574996948242188;
+  fNPads[0] = 71;
+  fNPads[1] = 71;
+  fNPads[2] = 71;
+  fNPads[3] = 73;
+  fNPads[4] = 73;
+  fNPads[5] = 73;
+  fNPads[6] = 75;
+  fNPads[7] = 75;
+  fNPads[8] = 75;
+  fNPads[9] = 77;
+  fNPads[10] = 77;
+  fNPads[11] = 77;
+  fNPads[12] = 79;
+  fNPads[13] = 79;
+  fNPads[14] = 79;
+  fNPads[15] = 81;
+  fNPads[16] = 81;
+  fNPads[17] = 81;
+  fNPads[18] = 83;
+  fNPads[19] = 83;
+  fNPads[20] = 83;
+  fNPads[21] = 85;
+  fNPads[22] = 85;
+  fNPads[23] = 85;
+  fNPads[24] = 87;
+  fNPads[25] = 87;
+  fNPads[26] = 87;
+  fNPads[27] = 89;
+  fNPads[28] = 89;
+  fNPads[29] = 89;
+  fNPads[30] = 89;
+  fNPads[31] = 91;
+  fNPads[32] = 91;
+  fNPads[33] = 91;
+  fNPads[34] = 93;
+  fNPads[35] = 93;
+  fNPads[36] = 93;
+  fNPads[37] = 95;
+  fNPads[38] = 95;
+  fNPads[39] = 95;
+  fNPads[40] = 97;
+  fNPads[41] = 97;
+  fNPads[42] = 97;
+  fNPads[43] = 99;
+  fNPads[44] = 99;
+  fNPads[45] = 99;
+  fNPads[46] = 101;
+  fNPads[47] = 101;
+  fNPads[48] = 101;
+  fNPads[49] = 103;
+  fNPads[50] = 103;
+  fNPads[51] = 103;
+  fNPads[52] = 105;
+  fNPads[53] = 105;
+  fNPads[54] = 105;
+  fNPads[55] = 73;
+  fNPads[56] = 73;
+  fNPads[57] = 73;
+  fNPads[58] = 75;
+  fNPads[59] = 75;
+  fNPads[60] = 75;
+  fNPads[61] = 75;
+  fNPads[62] = 77;
+  fNPads[63] = 77;
+  fNPads[64] = 77;
+  fNPads[65] = 79;
+  fNPads[66] = 79;
+  fNPads[67] = 79;
+  fNPads[68] = 81;
+  fNPads[69] = 81;
+  fNPads[70] = 81;
+  fNPads[71] = 81;
+  fNPads[72] = 83;
+  fNPads[73] = 83;
+  fNPads[74] = 83;
+  fNPads[75] = 85;
+  fNPads[76] = 85;
+  fNPads[77] = 85;
+  fNPads[78] = 85;
+  fNPads[79] = 87;
+  fNPads[80] = 87;
+  fNPads[81] = 87;
+  fNPads[82] = 89;
+  fNPads[83] = 89;
+  fNPads[84] = 89;
+  fNPads[85] = 91;
+  fNPads[86] = 91;
+  fNPads[87] = 91;
+  fNPads[88] = 91;
+  fNPads[89] = 93;
+  fNPads[90] = 93;
+  fNPads[91] = 93;
+  fNPads[92] = 95;
+  fNPads[93] = 95;
+  fNPads[94] = 95;
+  fNPads[95] = 95;
+  fNPads[96] = 97;
+  fNPads[97] = 97;
+  fNPads[98] = 97;
+  fNPads[99] = 99;
+  fNPads[100] = 99;
+  fNPads[101] = 99;
+  fNPads[102] = 101;
+  fNPads[103] = 101;
+  fNPads[104] = 101;
+  fNPads[105] = 101;
+  fNPads[106] = 103;
+  fNPads[107] = 103;
+  fNPads[108] = 103;
+  fNPads[109] = 105;
+  fNPads[110] = 105;
+  fNPads[111] = 105;
+  fNPads[112] = 105;
+  fNPads[113] = 107;
+  fNPads[114] = 107;
+  fNPads[115] = 107;
+  fNPads[116] = 109;
+  fNPads[117] = 109;
+  fNPads[118] = 109;
+  fNPads[119] = 111;
+  fNPads[120] = 111;
+  fNPads[121] = 111;
+  fNPads[122] = 111;
+  fNPads[123] = 113;
+  fNPads[124] = 113;
+  fNPads[125] = 113;
+  fNPads[126] = 115;
+  fNPads[127] = 115;
+  fNPads[128] = 115;
+  fNPads[129] = 115;
+  fNPads[130] = 117;
+  fNPads[131] = 117;
+  fNPads[132] = 117;
+  fNPads[133] = 119;
+  fNPads[134] = 119;
+  fNPads[135] = 119;
+  fNPads[136] = 121;
+  fNPads[137] = 121;
+  fNPads[138] = 121;
+  fNPads[139] = 121;
+  fNPads[140] = 123;
+  fNPads[141] = 123;
+  fNPads[142] = 123;
+  fNPads[143] = 125;
+  fNPads[144] = 125;
+  fNPads[145] = 125;
+  fNPads[146] = 125;
+  fNPads[147] = 127;
+  fNPads[148] = 127;
+  fNPads[149] = 127;
+  fNPads[150] = 129;
+  fNPads[151] = 129;
+  fNPads[152] = 129;
+  fNPads[153] = 129;
+  fNPads[154] = 131;
+  fNPads[155] = 131;
+  fNPads[156] = 131;
+  fNPads[157] = 133;
+  fNPads[158] = 133;
+  fNPads[159] = 133;
+  fNPads[160] = 135;
+  fNPads[161] = 135;
+  fNPads[162] = 135;
+  fNPads[163] = 135;
+  fNPads[164] = 137;
+  fNPads[165] = 137;
+  fNPads[166] = 137;
+  fNPads[167] = 139;
+  fNPads[168] = 139;
+  fNPads[169] = 139;
+  fNPads[170] = 139;
+  fNPads[171] = 141;
+  fNPads[172] = 141;
+  fNPads[173] = 141;
+}
+
+
+Double_t AliL3Transform::GetEta(Float_t *xyz)
+{
+  Double_t r3 = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]+xyz[2]*xyz[2]);
+  Double_t eta = 0.5 * log((r3+xyz[2])/(r3-xyz[2]));
+  return eta;
+}
+
+Double_t AliL3Transform::GetPhi(Float_t *xyz)
+{
+  
+  Double_t phi = atan2(xyz[1],xyz[0]);
+  //if(phi<0) phi=phi+2*TMath::Pi();
+  return phi;
+}
+
+
+Bool_t AliL3Transform::Slice2Sector(Int_t slice, Int_t slicerow, Int_t & sector, Int_t &row) const{
+  if(slicerow<0&&slicerow>=fNRow) return kFALSE;
+  if(slice<0||slice>=fNSlice) return kFALSE;
+
+  if(slicerow<fNRowLow){
+    sector = slice;
+    row    = slicerow;
+  }
+  else {
+    sector = slice+fNSlice;
+    row    = slicerow-fNRowLow;
+  }
+  return kTRUE;
+}
+
+Bool_t AliL3Transform::Sector2Slice(Int_t & slice, Int_t  sector) const{
+  if(sector<0||sector>=fNSector) return kFALSE;
+  if(sector<fNSectorLow) slice = sector;
+  else          slice = sector - fNSectorLow;
+  return kTRUE;
+}
+
+Bool_t AliL3Transform::Sector2Slice(Int_t & slice, Int_t & slicerow,Int_t  sector, Int_t row) const{
+  if(sector<0||sector>=fNSector||row<0) return kFALSE;
+  if(sector<fNSectorLow){
+    if(row>=fNRowLow) return kFALSE;
+    slice = sector;
+    slicerow = row;
+  }
+  else{
+    if(row>=fNRowUp) return kFALSE;
+    slice = sector - fNSectorLow;
+    slicerow = row + fNRowLow;
+  }
+  return kTRUE;
+}
+
+Double_t AliL3Transform::Row2X(Int_t slicerow){
+  if(slicerow<0||slicerow>=fNRow) return 0;
+  return fX[slicerow];
+}
+
+void AliL3Transform::Local2Global(Float_t *xyz,Int_t slice)
+{
+  //Transformation to global coordinate system
+  Float_t x0 = xyz[0];
+  Float_t y0 = xyz[1];
+  Float_t cs,sn;
+  cs = fCos[slice];
+  sn = fSin[slice];
+  xyz[0]=x0*cs-y0*sn;
+  xyz[1]=x0*sn+y0*cs;
+  xyz[2]=xyz[2];//global z=local z
+}
+
+void AliL3Transform::Local2GlobalAngle(Float_t *angle,Int_t slice){
+  angle[0] = fmod(angle[0]+slice*fPi/9,2*fPi);
+}
+
+void AliL3Transform::Global2LocalAngle(Float_t *angle,Int_t slice){
+  angle[0] = angle[0]-slice*fPi/9;
+  if(angle[0]<0) angle[0]+=2*fPi;
+}
+
+void AliL3Transform::Raw2Local(Float_t *xyz,Int_t sector,Int_t row,Float_t pad,Float_t time)
+{
+  //Transformation from rawdata to local coordinate system
+  
+  Int_t slice,slicerow;
+  Sector2Slice(slice, slicerow, sector, row);  
+
+  xyz[0]=Row2X(slicerow); 
+  Int_t npads= fNPads[slicerow];
+  if(sector<fNSectorLow)
+    xyz[1]=(pad-0.5*(npads-1))*fPadPitchWidthLow;
+  else
+    xyz[1]=(pad-0.5*(npads-1))*fPadPitchWidthUp;
+  xyz[2]=fZWidth*time-3.*fZSigma;
+  Int_t sign=-1;
+  Int_t nis=fNSectorLow;
+  Int_t nos=fNSectorUp;
+  
+  if((sector<nis)/2 || ((sector-nis)<nos/2)) sign=1;
+  xyz[2]=sign*(250.-xyz[2]);
+
+}
+
+void AliL3Transform::Local2Global(Float_t *xyz,Int_t sector,Int_t row)
+{
+  //Transformation to global coordinate system
+  Int_t slice,slicerow;
+  Sector2Slice(slice, slicerow, sector, row);  
+  Float_t r=Row2X(slicerow);
+  Float_t cs = fCos[slice];
+  Float_t sn = fSin[slice];
+
+  xyz[0]=r*cs-xyz[1]*sn;
+  xyz[1]=r*sn+xyz[1]*cs;
+  xyz[2]=xyz[2];//global z=local z
+}
+
+Double_t AliL3Transform::GetMaxY(Int_t slicerow)
+{
+ if(slicerow < fNRowLow)
+     return fPadPitchWidthLow*fNPads[slicerow]/2; 
+ else
+     return fPadPitchWidthUp*fNPads[slicerow]/2;
+}
+
+void AliL3Transform::Global2Local(Float_t *xyz,Int_t sector)
+{
+  Int_t slice;
+  Sector2Slice(slice, sector);  
+  Float_t cs = fCos[slice];
+  Float_t sn = fSin[slice];
+  Float_t x1 = xyz[0]*cs + xyz[1]*sn;
+  Float_t y1 = -xyz[0]*sn + xyz[1]*cs;
+  xyz[0] = x1;
+  xyz[1] = y1;
+}
+
+void AliL3Transform::Raw2Global(Float_t *xyz,Int_t sector,Int_t row,Float_t pad,Float_t time)
+{
+  //Transformation from raw to global coordinates
+  
+  Raw2Local(xyz,sector,row,pad,time);
+  Local2Global(xyz,sector,row);
+}
+
+void AliL3Transform::Local2Raw(Float_t *xyz,Int_t sector,Int_t row)
+{
+  //Transformation from local coordinates to raw
+  
+  Int_t slice,slicerow;
+  Sector2Slice(slice, slicerow, sector, row);  
+   
+  if(sector<fNSectorLow)
+    xyz[1]=xyz[1]/fPadPitchWidthLow+0.5*(fNPads[slicerow]-1);
+  else
+    xyz[1]=xyz[1]/fPadPitchWidthUp+0.5*(fNPads[slicerow]-1);
+  Int_t sign=-1;
+  Int_t nis=fNSectorLow;
+  Int_t nos=fNSectorUp;
+  if ((sector<nis/2) || ((sector-nis)<nos/2)) sign=1; 
+  xyz[2]=250-sign*xyz[2];
+  xyz[2]=(xyz[2]+3.*fZSigma)/fZWidth;
+}
+
+void AliL3Transform::Global2Raw(Float_t *xyz,Int_t sector,Int_t row)
+{
+  //Transformation from global coordinates to raw. 
+
+  Global2Local(xyz,sector);
+  Local2Raw(xyz,sector,row);
+
+}
diff --git a/HLT/hough/AliL3Transform.h b/HLT/hough/AliL3Transform.h
new file mode 100644 (file)
index 0000000..a32bf6c
--- /dev/null
@@ -0,0 +1,65 @@
+#ifndef ALIL3TRANSFORM_H
+#define ALIL3TRANSFORM_H
+
+
+
+//#include "AliTPCParam.h"
+#include "AliL3RootTypes.h"
+class TFile;
+
+class AliL3Transform {
+ private:
+//  AliTPCParam *fParam;
+  Int_t fNTimeBins;
+  Int_t fNRowLow;
+  Int_t fNRowUp;
+  Int_t fNSectorLow;
+  Int_t fNSectorUp;
+  Double_t fPadPitchWidthLow;
+  Double_t fPadPitchWidthUp;
+  Double_t fZWidth;
+  Double_t fZSigma;
+  Int_t fNSector;
+  Int_t fNSlice;
+  Int_t fNRow;
+  Double_t fPi;
+  Double_t fCos[36];
+  Double_t fSin[36];
+  Double_t fX[174];
+  Int_t fNPads[174];
+ public:
+  AliL3Transform();
+  virtual ~AliL3Transform();
+  void Init();
+
+  Double_t GetPadPitchWidthLow() {return fPadPitchWidthLow;}
+  Double_t GetPadPitchWidthUp() {return fPadPitchWidthUp;}
+
+  Bool_t Slice2Sector(Int_t slice, Int_t slicerow, Int_t & sector, Int_t &row) const;
+
+  Bool_t Sector2Slice(Int_t & slice, Int_t  sector) const;
+  Bool_t Sector2Slice(Int_t & slice, Int_t & slicerow,Int_t  sector, Int_t row) const;
+  
+  Double_t Row2X(Int_t slicerow);
+  Int_t GetNPads(Int_t row){return (row<174)?fNPads[row]:0;}
+  Int_t GetNTimeBins(){return fNTimeBins;}
+
+  Double_t GetEta(Float_t *xyz);
+  Double_t GetPhi(Float_t *xyz);
+  Double_t GetMaxY(Int_t slicerow);
+  void Local2Global(Float_t *xyz,Int_t slice);
+  void Local2GlobalAngle(Float_t *angle,Int_t slice);
+  void Global2LocalAngle(Float_t *angle,Int_t slice);
+
+  void Raw2Local(Float_t *xyz,Int_t sector,Int_t row,Float_t pad,Float_t time);
+  void Local2Global(Float_t *xyz,Int_t sector,Int_t row);
+  void Global2Local(Float_t *xyz,Int_t sector);
+  void Raw2Global(Float_t *xyz,Int_t sector,Int_t row,Float_t pad,Float_t time);
+  void Local2Raw(Float_t *xyz,Int_t sector,Int_t row);
+  void Global2Raw(Float_t *xyz,Int_t sector,Int_t row);
+  
+  ClassDef(AliL3Transform,1)
+};
+
+
+#endif
diff --git a/HLT/hough/Makefile b/HLT/hough/Makefile
new file mode 100644 (file)
index 0000000..a4b97f7
--- /dev/null
@@ -0,0 +1,79 @@
+############################### TPC Makefile ##################################
+
+# Include machine specific definitions
+
+include $(ALICE_ROOT)/conf/GeneralDef
+include $(ALICE_ROOT)/conf/MachineDef.$(ALICE_TARGET)
+
+PACKAGE = AliL3Hough
+
+# C++ sources
+
+
+SRCS          = AliL3HoughTransformer.cxx AliL3HoughPixel.cxx  \
+                AliL3HoughMaxFinder.cxx AliL3HoughEval.cxx AliL3HoughMerge.cxx 
+                
+
+# C++ Headers
+
+HDRS          = $(SRCS:.cxx=.h) AliL3Defs.h AliL3HoughLinkDef.h
+# Library dictionary
+
+DICT          = AliL3HoughCint.cxx
+DICTH         = $(DICT:.cxx=.h)
+DICTO         = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(DICT))
+
+# FORTRAN Objectrs
+
+FOBJS         = $(FSRCS:.f=.o)
+
+# C Objects
+
+COBJS         = $(CSRCS:.c=.o)
+
+# C++ Objects
+
+OBJS          = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(SRCS)) $(DICTO)
+
+# C++ compilation flags
+
+CXXFLAGS      = $(CXXOPTS) -g -I$(ROOTSYS)/include -I. -I$(ALICE_ROOT)/include/ -I$(ALICE_ROOT)/TPC \
+               -I$(ALICE_ROOT)/CONTAINERS -I$(ALICE)/level3code/src    
+
+#CXXFLAGS      = $(CXXOPTS) -g -Wall -I$(ROOTSYS)/include -I. -I $(ALICE_ROOT)/TPC -I$(ALICE_ROOT)/include/ -DCOMPILING
+# FORTRAN compilation flags
+
+FFLAGS      = $(FOPT)
+##### TARGETS #####
+# Target
+
+SLIBRARY       = $(ALICE)/mylibs/libAliL3Hough.$(SL)
+
+default:       $(SLIBRARY)
+
+$(ALICE)/mylibs/libAliL3Hough.$(SL):  $(OBJS)
+
+$(DICT):               $(HDRS)
+
+depend:                        $(SRCS)
+
+TOCLEAN                        = $(OBJS) *Cint.h *Cint.cxx
+
+############################### General Macros ################################
+
+include $(ALICE_ROOT)/conf/GeneralMacros
+
+############################ Dependencies #####################################
+
+include tgt_$(ALICE_TARGET)/Make-depend 
+
+###########
+so:
+       rm -fr $(ALICE)/mylibs/libAliL3Hough.so
+clean:
+       rm -fr tgt_Linux/*.o
+       rm -fr $(ALICE)/mylibs/libAliL3Hough.so
+       rm -fr $(DICT) $(DICTH) $(DICTO)
diff --git a/HLT/hough/hough.C b/HLT/hough/hough.C
new file mode 100644 (file)
index 0000000..8b1332d
--- /dev/null
@@ -0,0 +1,143 @@
+void hough(char *rootfile,int patch,Bool_t MC=false)
+{
+  gStyle->SetOptStat(0);
+  
+  Int_t xbin = 60;
+  Int_t ybin = 60;
+  Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2->
+  //Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88
+  Float_t yrange[2] = {-0.26 , 0.26};// -15 - 15
+  //Float_t yrange[2] = {0.55 , 0.88}; //slice 2 0.55->0.88
+  TH2F *hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  hist->GetXaxis()->SetTitle("#kappa [cm^{-1}]");
+  hist->GetYaxis()->SetTitle("#Phi [rad]");
+  SetTH1Options(hist);
+
+  Int_t xr[2] = {0,250};
+  Int_t yr[2] = {-125,125};
+  TH1F *ntracks = new TH1F("ntracks","",100,0,200);
+  TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
+  TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
+  TH2F *fake = new TH2F("fake","",250,0,250,250,-125,125);
+  TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  peaks->SetMarkerStyle(3);
+  peaks->SetMarkerColor(2);
+  
+  int slice = 2,n_phi_segments=30;
+  //float eta[2] = {0.3,0.4};
+  float eta[2] = {0.03,0.04};
+  
+  AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments);
+  a->GetPixels(rootfile,raw);
+  a->InitTemplates(hist);
+  
+  //a->Transform2Circle(hist,78);
+  TStopwatch sw;
+  a->Transform2Circle(hist);
+  sw.Stop();
+  printf("Transformation done in %f ms\n",sw.CpuTime()*1000);
+  
+  AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("KappaPhi");
+  //b->SetThreshold(10000);
+  AliL3TrackArray *tracks = b->FindMaxima(hist);
+    
+  AliL3HoughEval *c = new AliL3HoughEval(a);
+  printf("number of peaks before %d\n",tracks->GetNTracks());
+  
+  /*
+    for(int i=0; i<tracks->GetNTracks(); i++)
+    {
+      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+      if(!track) printf("no track\n");
+      printf("Pt %f Phi %f kappa %f weigth %d charge %d\n",track->GetPt(),track->GetPhi0(),track->GetKappa(),track->GetNHits(),track->GetCharge());
+    }
+  */
+  //c->LookInsideRoad(a,tracks,road,fake,ntracks);
+  c->LookInsideRoad(tracks);
+  printf("Number of trackcandidates %d\n",tracks->GetNTracks());
+
+  AliL3Transform *transform = new AliL3Transform();
+  Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
+  for(Int_t i=0; i<tracks->GetNTracks(); i++)
+    {
+      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+      if(!track) {printf("NO TRACK\n"); continue;}
+      peaks->Fill(track->GetKappa(),track->GetPhi0(),1);
+      //for(int row=0; row<174; row++)
+      for(int row=NRows[patch][0]; row<=NRows[patch][1]; row++)
+       {
+         Float_t xyz[3];
+         track->GetCrossingPoint(row,xyz);
+         //track->GetCrossingPoint(slice,row,xyz);
+           //printf("Track does not cross line at row %d\n",row);
+         //transform->Local2Global(xyz,slice);
+         road->Fill(xyz[0],xyz[1],1);
+       }
+      //float an[2] = {track->GetPhi0(),0};
+      //transform->Local2GlobalAngle(an,slice);
+      printf("pt %f phi0 %f weight %d\n",track->GetPt(),track->GetPhi0(),track->GetWeight());
+    }
+  
+  
+  TCanvas *c1 = new TCanvas("c1","",900,900);
+  SetCanvasOptions(c1);
+  c1->Divide(2,2);
+  c1->cd(1);
+  hist->Draw("lego1");
+  //  peaks->Draw("same");
+
+  //  TCanvas *c2 = new TCanvas("c2","",2);
+  c1->cd(3);
+  raw->Draw("");
+
+  //TCanvas *c3 = new TCanvas("c3","",2);
+  //road->SetMarkerStyle(5);
+  //road->SetMarkerColor(3);
+  c1->cd(4);
+  road->SetMarkerStyle(9);
+  road->Draw();
+  /*
+  TCanvas *c4 = new TCanvas("c4","",2);
+  fake->Draw("cont1");
+
+  TCanvas *c5 = new TCanvas("c5","",2);
+  ntracks->Draw();
+  */
+  delete tracks;
+  delete a;
+  delete b;
+  delete c;
+
+  //----------------------------------------------------------------------------
+  if(MC)
+    {
+      TFile *file = new TFile(rootfile);
+      file->cd();
+      
+      gAlice = (AliRun*)file->Get("gAlice");
+      gAlice->GetEvent(0);
+      
+      TClonesArray *particles=gAlice->Particles();
+      Int_t n=particles->GetEntriesFast();
+      Float_t torad=TMath::Pi()/180;
+      Float_t phi_min = slice*20 - 10;
+      Float_t phi_max = slice*20 + 10;
+      
+      for (Int_t j=0; j<n; j++) {
+       TParticle *p=(TParticle*)particles->UncheckedAt(j);
+       if (p->GetFirstMother()>=0) continue;  //secondary particle
+       if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
+       Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
+       Double_t phi_part = TMath::ATan2(pyg,pxg);
+       if (phi_part < 0) phi_part += 2*TMath::Pi();
+       
+       if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
+       if(ptg<0.100) continue;
+               
+       printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part);
+      }
+      file->Close();
+      delete file;
+    }
+  
+}
diff --git a/HLT/hough/hough_line.C b/HLT/hough/hough_line.C
new file mode 100644 (file)
index 0000000..89e8f41
--- /dev/null
@@ -0,0 +1,126 @@
+void hough_line(char *rootfile,int patch,Bool_t MC=false)
+{
+  gStyle->SetOptStat(0);
+  
+  
+  Double_t pi = TMath::Pi();
+  Int_t xbin = 60;
+  Int_t ybin = 60;
+  Float_t xrange[2] = {-pi/2,pi/2};
+  Float_t yrange[2] = {-40,40};
+  
+  TH2F *hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  hist->GetXaxis()->SetTitle("#Psi");
+  hist->GetYaxis()->SetTitle("D");
+  SetTH1Options(hist);
+
+  Int_t xr[2] = {0,250};
+  Int_t yr[2] = {-125,125};
+  TH1F *ntracks = new TH1F("ntracks","",100,0,200);
+  TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
+  TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
+  TH2F *container = new TH2F("container","",250,0,250,250,yr[0],yr[1]);
+  TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  peaks->SetMarkerStyle(3);
+  peaks->SetMarkerColor(2);
+  
+  int slice = 2,n_phi_segments=20;
+  float eta[2] = {0.3,0.4};
+  //float eta[2] = {0.05,0.06};
+  
+  AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments);
+
+  const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
+  const int ref_row[5]={22,60,93,125,157};
+  int rowrange[2] = {6,12};
+  double phirange[2] = {-10.,10.};
+  a->GetPixels(rootfile,raw);
+  a->Transform2Line(hist,9,rowrange,phirange);
+
+  AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("DPsi");
+  AliL3TrackArray *tracks = b->FindMaxima(hist,rowrange,5);
+
+  TH2F *cross = new TH2F("cross","",250,xr[0],xr[1],250,yr[0],yr[1]);
+  for(int i=0; i<2; i++)
+  //for(int i=0; i<tracks->GetNTracks(); i++)
+    {
+      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+      if(!track) continue;
+      Double_t xy[2];
+      for(int j=0; j<174; j++)
+       {
+         track->GetLineCrossingPoint(j,xy);
+         cross->Fill(xy[0],xy[1],1);
+       }
+      printf("Segment %f %f weight %d\n",track->GetPsiLine(),track->GetDLine(),track->GetWeight());
+    }
+  cross->SetMarkerStyle(6);
+  cross->SetMarkerColor(2);
+  
+  
+  TCanvas *c1 = new TCanvas("c1","",900,900);
+  SetCanvasOptions(c1);
+  c1->Divide(2,2);
+  c1->cd(1);
+  hist->Draw("box");
+  peaks->SetMarkerStyle(3);
+  peaks->SetMarkerColor(2);
+  //peaks->Draw("same");
+
+  c1->cd(2);
+  container->Draw();
+
+  c1->cd(3);
+  raw->Draw("");
+  cross->Draw("same");
+  return;
+  //TCanvas *c3 = new TCanvas("c3","",2);
+  //road->SetMarkerStyle(5);
+  //road->SetMarkerColor(3);
+  c1->cd(4);
+  road->SetMarkerStyle(9);
+  road->Draw();
+  /*
+  TCanvas *c4 = new TCanvas("c4","",2);
+  fake->Draw("cont1");
+
+  TCanvas *c5 = new TCanvas("c5","",2);
+  ntracks->Draw();
+  */
+  delete a;
+  delete b;
+  delete c;
+
+  //----------------------------------------------------------------------------
+  if(MC)
+    {
+      TFile *file = new TFile(rootfile);
+      file->cd();
+      
+      gAlice = (AliRun*)file->Get("gAlice");
+      gAlice->GetEvent(0);
+      
+      TClonesArray *particles=gAlice->Particles();
+      Int_t n=particles->GetEntriesFast();
+      Float_t torad=TMath::Pi()/180;
+      Float_t phi_min = slice*20 - 10;
+      Float_t phi_max = slice*20 + 10;
+      
+      for (Int_t j=0; j<n; j++) {
+       TParticle *p=(TParticle*)particles->UncheckedAt(j);
+       if (p->GetFirstMother()>=0) continue;  //secondary particle
+       if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
+       Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
+       Double_t phi_part = TMath::ATan2(pyg,pxg);
+       if (phi_part < 0) phi_part += 2*TMath::Pi();
+       
+       if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
+       if(ptg<0.100) continue;
+               
+       printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part);
+      }
+      file->Close();
+      delete file;
+    }
+  
+}
diff --git a/HLT/hough/hough_line_merge.C b/HLT/hough/hough_line_merge.C
new file mode 100644 (file)
index 0000000..f24f069
--- /dev/null
@@ -0,0 +1,205 @@
+void hough_line_merge(char *rootfile,int patch,Bool_t MC=false)
+{
+  gStyle->SetOptStat(0);
+  
+  
+  Double_t pi = TMath::Pi();
+  Int_t xbin = 60;
+  Int_t ybin = 60;
+  Float_t xrange[2] = {-pi/2,pi/2};
+  Float_t yrange[2] = {-40,40};
+  
+  TH2F *hist;
+  
+  Int_t xr[2] = {0,250};
+  Int_t yr[2] = {-125,125};
+  TH1F *ntracks = new TH1F("ntracks","",100,0,200);
+  TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
+  TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
+  TH2F *fake = new TH2F("fake","",250,0,250,250,-125,125);
+  TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  peaks->SetMarkerStyle(3);
+  peaks->SetMarkerColor(2);
+  
+  int slice = 2,n_phi_segments=20;
+  float eta[2] = {0.3,0.4};
+  //float eta[2] = {0.05,0.06};
+  AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments);
+  a->GetPixels(rootfile,raw);
+  
+  AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("DPsi");
+  
+  const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
+  //const int ref_row[4]={5,15,25,35};
+  //int rowrange[4][2]={{0,10},{10,20},{20,30},{30,45}};
+  int ref_row[7] = {3,9,15,21,27,33,39};
+  int rowrange[7][2] = {{0,6},{6,12},{12,18},{18,24},{24,30},{30,36},{36,42}};
+  
+  double phirange[7][2] = {{-10.,-5.},{-7.5,-2.5},{-5.,0.},{-2.5,2.5},{0.,5.},{2.5,7.5},{5.,10.}};
+  int npatches = 7;
+  int count=0;
+
+  //AliL3TrackArray **track_pieces = new AliL3TrackArray*[npatches];
+  AliL3TrackArray *track_pieces = new AliL3TrackArray("AliL3HoughTrack");
+
+  double pran[2] = {-10,10};
+  
+  // int pa=3;
+  for(int row_pat=0; row_pat < 7; row_pat++)
+    {
+      
+      for(int pat=0; pat<7; pat++)
+       {
+         
+         TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+         
+         printf("transforming ref_row %d rowrange %d %d phirange %f %f\n",ref_row[row_pat],rowrange[row_pat][0],rowrange[row_pat][1],phirange[pat][0],phirange[pat][1]);
+         
+         a->Transform2Line(histo,ref_row[row_pat],rowrange[row_pat],phirange[pat]);
+         //a->Transform2Line(histo,ref_row[row_pat],rowrange[row_pat],pran);
+
+         AliL3TrackArray *tmp = (AliL3TrackArray*)b->FindMaxima(histo,rowrange[row_pat],ref_row[row_pat]);
+         
+         track_pieces->AddTracks(tmp);
+         
+         //AliL3HoughTrack *tra = (AliL3HoughTrack*)b->GetTracks()->GetCheckedTrack(0);
+         //printf("psi %f rowrange %d %d\n",tra->GetPsiLine(),tra->GetFirstRow(),tra->GetLastRow());
+         
+         delete histo;
+         delete tmp;
+       }
+      
+    }
+  
+  
+  printf("Transformation complete, genereted %d tracks\n",track_pieces->GetNTracks());
+
+    
+  //track_pieces->Sort();
+  TH2F *cross = new TH2F("cross","",250,xr[0],xr[1],250,yr[0],yr[1]);  
+  Double_t xy[2];
+  //for(int i=0; i<npatches; i++)
+  for(int i=0; i<4; i++)
+    {
+      if(track_pieces->GetNTracks()==0) continue;
+      AliL3HoughTrack *track = (AliL3HoughTrack*)track_pieces->GetCheckedTrack(i);
+      if(!track) {printf("NO TRACK\n"); continue;}
+      //for(int j=rowrange[pa][0]; j<=rowrange[pa][1]; j++)
+      /*
+      for(int j=0; j<174; j++)
+       {
+         track->GetLineCrossingPoint(j,xy);
+         cross->Fill(xy[0],xy[1],1);
+       }
+      */
+      //printf("weight %d psi %f D %f rowrange %d %d\n",track->GetWeight(),track->GetPsiLine(),track->GetDLine(),track->GetFirstRow(),track->GetLastRow());
+      
+    }
+  cross->SetMarkerStyle(6);
+  cross->SetMarkerColor(2);
+  //cross->Draw();
+  
+  
+
+  AliL3HoughTransformer *a2 = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments);
+  
+  xbin = 60 ;
+  ybin = 60;
+  xrange[0] = -0.006;
+  xrange[1] = 0.006;
+  yrange[0] = -0.26;
+  yrange[1] = 0.26;
+  hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  
+  TH2F *cross_helix = new TH2F("cross_helix","",250,xr[0],xr[1],250,yr[0],yr[1]);  
+  
+  a2->TransformLines2Circle(hist,(AliL3TrackArray*)track_pieces);
+      
+  AliL3HoughMaxFinder *b2 = new AliL3HoughMaxFinder("KappaPhi");
+  AliL3TrackArray *final_tracks = b2->FindMaxima(hist);
+  
+  printf("Number of track candidates before checking %d\n",final_tracks->GetNTracks());
+  AliL3HoughEval *eval = new AliL3HoughEval(slice);
+  eval->SetTransformer(a);
+  eval->LookInsideRoad(final_tracks);
+  
+  printf("Total number of tracks after checking %d\n",final_tracks->GetNTracks());
+  for(int i=0; i<3; i++)//final_tracks->GetNTracks(); i++)
+    {
+      
+      AliL3HoughTrack *track = (AliL3HoughTrack*)final_tracks->GetCheckedTrack(i);
+      if(!track) continue;
+      
+      for(int r=0; r<174; r++)
+       {
+         Float_t xyz_helix[3];
+         track->GetCrossingPoint(r,xyz_helix);
+         cross_helix->Fill(xyz_helix[0],xyz_helix[1],1);
+       }
+      
+      peaks->Fill(track->GetKappa(),track->GetPhi0(),1);
+      printf("Found track, pt %f phi0 %f\n",track->GetPt(),track->GetPhi0());
+    }
+  cross_helix->SetMarkerStyle(6);
+  cross_helix->SetMarkerColor(2);
+  TCanvas *c1 = new TCanvas("c1","",900,900);
+  SetCanvasOptions(c1);
+  c1->Divide(2,2);
+  c1->cd(1);
+  hist->Draw("box");
+
+  peaks->SetMarkerStyle(3);
+  peaks->SetMarkerColor(2);
+  //peaks->Draw("same");
+
+  c1->cd(2);
+  cross_helix->Draw();
+
+  c1->cd(3);
+  raw->Draw("");
+  cross_helix->DrawCopy("same");
+  //cross->Draw("same");
+
+  c1->cd(4);
+  //hist->Draw("box");
+  raw->DrawCopy();
+  
+  delete track_pieces;
+  delete a;
+  delete b;
+  delete a2;
+  delete eval;
+  //----------------------------------------------------------------------------
+  if(MC)
+    {
+      TFile *file = new TFile(rootfile);
+      file->cd();
+      
+      gAlice = (AliRun*)file->Get("gAlice");
+      gAlice->GetEvent(0);
+      
+      TClonesArray *particles=gAlice->Particles();
+      Int_t n=particles->GetEntriesFast();
+      Float_t torad=TMath::Pi()/180;
+      Float_t phi_min = slice*20 - 10;
+      Float_t phi_max = slice*20 + 10;
+      
+      for (Int_t j=0; j<n; j++) {
+       TParticle *p=(TParticle*)particles->UncheckedAt(j);
+       if (p->GetFirstMother()>=0) continue;  //secondary particle
+       if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
+       Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
+       Double_t phi_part = TMath::ATan2(pyg,pxg);
+       if (phi_part < 0) phi_part += 2*TMath::Pi();
+       
+       if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
+       if(ptg<0.100) continue;
+               
+       printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part);
+      }
+      file->Close();
+      delete file;
+    }
+  
+}
diff --git a/HLT/hough/hough_merge.C b/HLT/hough/hough_merge.C
new file mode 100644 (file)
index 0000000..9eaffcb
--- /dev/null
@@ -0,0 +1,148 @@
+void hough_merge(char *rootfile,Bool_t MC=false)
+{
+  gStyle->SetOptStat(0);
+  
+  Int_t xbin = 60;
+  Int_t ybin = 60;
+  Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.1->4GeV 0.006-0.006
+  //Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88
+  Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
+  //Float_t yrange[2] = {0.55 , 0.88}; //slice 2 0.55->0.88
+
+  Int_t xr[2] = {0,250};
+  Int_t yr[2] = {-125,125};
+  TH1F *ntracks = new TH1F("ntracks","",100,0,200);
+  TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
+  TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
+  TH2F *fake = new TH2F("fake","",300,0,300,250,-125,125);
+  TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  peaks->SetMarkerStyle(3);
+  peaks->SetMarkerColor(2);
+  
+  int slice = 2,patch=0,n_phi_segments=40;
+  float eta[2] = {0.3,0.4};
+  // float eta[2] = {0.04,0.05};
+  
+  AliL3HoughTransformer *a;
+  AliL3HoughMaxFinder *b;
+  AliL3HoughEval *c = new AliL3HoughEval(slice);
+  AliL3HoughMerge *d = new AliL3HoughMerge(slice);
+  TH2F *hist;
+  for(int pat=0; pat<5; pat++)
+    {
+      a = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments);
+      hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+      a->GetPixels(rootfile,raw);
+      a->Transform2Circle(hist,87);
+      
+      b = new AliL3HoughMaxFinder("KappaPsi");
+      AliL3TrackArray *tracks = b->FindMaxima(hist);
+      c->SetTransformer(a);
+      c->LookInsideRoad(tracks);
+      d->FillTracks(tracks,pat);
+      delete a;
+      delete b;
+    }
+  return;
+  
+
+  TH2F *merge_hist = new TH2F("merge_hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  
+  d->FillHisto(merge_hist);
+  
+  
+  b = new AliL3HoughMaxFinder(slice,pat);
+  b->FindMaxima(merge_hist);
+  
+  AliL3TrackArray *mtrs = b->GetTracks();
+  
+  c->CompareMC(rootfile,mtrs,eta);
+  
+  AliL3Transform *transform = new AliL3Transform();
+  for(int tr=0; tr<mtrs->GetNTracks(); tr++)
+    {
+      AliL3HoughTrack *t = (AliL3HoughTrack*)mtrs->GetCheckedTrack(tr);
+      if(!t) {printf("NO TRACK11\n"); continue;}
+      if(t->GetMCid() < 0) {printf("Track %d was fake, weigth %d\n",tr,t->GetNHits()); continue;}
+      printf("Found pt %f weigth %d\n",t->GetPt(),t->GetNHits());
+      //printf("charge %f\n",t->GetCharge());
+      for(Int_t j=0; j<174; j++)
+       {
+         Float_t xyz[3];
+         if(!t->GetCrossingPoint(j,xyz)) 
+           {
+             printf("Track does not cross line\n");
+             continue;
+           }
+         road->Fill(xyz[0],xyz[1],1);
+         
+       }
+       
+      
+    }
+  
+
+  TCanvas *c1 = new TCanvas("c1","",800,800);
+  SetCanvasOptions(c1);
+  c1->Divide(2,2);
+  c1->cd(1);
+  merge_hist->Draw("box");
+  
+  //TCanvas *c2 = new TCanvas("c2","",2);
+  c1->cd(3);
+  raw->Draw();
+
+  //TCanvas *c3 = new TCanvas("c3","",2);
+  c1->cd(4);
+  road->SetMarkerStyle(7);
+  road->Draw();
+
+  delete b;
+  delete d;
+  delete c;
+  //----------------------------------------------------------------------------
+  if(MC)
+    {
+      TFile *file = new TFile(rootfile);
+      file->cd();
+      
+      gAlice = (AliRun*)file->Get("gAlice");
+      gAlice->GetEvent(0);
+      
+      TClonesArray *particles=gAlice->Particles();
+      Int_t n=particles->GetEntriesFast();
+      Float_t torad=TMath::Pi()/180;
+      Float_t phi_min = slice*20 - 10;
+      Float_t phi_max = slice*20 + 10;
+      
+      for (Int_t j=0; j<n; j++) {
+       TParticle *p=(TParticle*)particles->UncheckedAt(j);
+       if (p->GetFirstMother()>=0) continue;  //secondary particle
+       if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
+       Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
+       Double_t phi_part = TMath::ATan2(pyg,pxg);
+       if (phi_part < 0) phi_part += 2*TMath::Pi();
+       
+       if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
+       if(ptg<0.100) continue;
+       int found = 0;
+       for(int m=0; m<mtrs->GetNTracks(); m++)
+         {
+           AliL3HoughTrack *tra = (AliL3HoughTrack*)mtrs->GetCheckedTrack(m);
+           if(tra->GetMCid() != j) continue;
+           found = 1;
+           double dPt = fabs(ptg-tra->GetPt());
+           float  phi0[2] = {tra->GetPhi0(),0};
+           transform->Local2GlobalAngle(phi0,slice);
+           double dPhi0 = fabs(phi_part-phi0[0]);
+           printf("Difference Pt %f Phi0 %f\n",dPt,dPhi0);
+
+         }
+       if(!found) printf("Track %d was not found\n",j);
+       printf("Generated particle %d, with pt %f phi %f\n",p->GetPdgCode(),ptg,phi_part);
+      }
+      file->Close();
+      delete file;
+    }
+  
+}
diff --git a/HLT/hough/hough_mergehistos.C b/HLT/hough/hough_mergehistos.C
new file mode 100644 (file)
index 0000000..854b50b
--- /dev/null
@@ -0,0 +1,63 @@
+void hough_mergehistos(char *rootfile)
+{
+  gStyle->SetOptStat(0);
+  
+  Int_t xbin = 60;
+  Int_t ybin = 60;
+  Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.1->4GeV 0.006-0.006
+  //Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88
+  Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
+  //Float_t yrange[2] = {0.55 , 0.88}; //slice 2 0.55->0.88
+
+  Int_t xr[2] = {0,250};
+  Int_t yr[2] = {-125,125};
+  TH1F *ntracks = new TH1F("ntracks","",100,0,200);
+  TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
+  TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
+  TH2F *fake = new TH2F("fake","",300,0,300,250,-125,125);
+  TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  peaks->SetMarkerStyle(3);
+  peaks->SetMarkerColor(2);
+  
+  int slice = 2,patch=0,n_phi_segments=40;
+  float eta[2] = {0.3,0.4};
+  // float eta[2] = {0.04,0.05};
+  
+  AliL3HoughTransformer *a;
+  AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("KappaPsi");
+  
+  TH2F **hist = new TH2F*[5];
+  
+  for(int pat=0; pat<5; pat++)
+    {
+      a = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments);
+      hist[pat] = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+      a->GetPixels(rootfile,raw);
+      a->InitTemplates(hist[pat]);
+      a->Transform2Circle(hist[pat]);
+      
+      //AliL3TrackArray *tracks = b->FindMaxima(hist[pat]);
+      
+      delete a;
+      
+    }
+  
+  
+  TCanvas *c1 = new TCanvas("c1","",800,800);
+  c1->Divide(2,2);
+  c1->cd(1);
+  hist[0]->DrawCopy("box");
+  c1->cd(2);
+  hist[4]->DrawCopy("box");
+  c1->cd(3);
+  hist[2]->DrawCopy("box");
+  
+  c1->cd(4);
+  hist[4]->Add(hist[0],hist[2]);
+  hist[4]->Draw("box");
+
+  delete b;
+  
+  
+  delete [] hist;
+}
diff --git a/HLT/hough/rootlogon.C b/HLT/hough/rootlogon.C
new file mode 100644 (file)
index 0000000..82a2c4d
--- /dev/null
@@ -0,0 +1,21 @@
+
+{
+  printf("\nWELCOME to the magic world of Level3\n\n"); 
+
+  //load library neccesary for TPC alone 
+  gSystem->Load("$(ROOTSYS)/lib/libPhysics");
+  gSystem->Load("$(ROOTSYS)/lib/libEG");
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libSTEER");
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libCONTAINERS");
+  gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTPC");
+  gSystem->Load("/usr/local/franken/lib/MLUC/lib/linux-i386/libMLUC.so");
+  
+  gSystem->Load("$(ALICE)/mylibs/libAliL3");
+  gSystem->Load("$(ALICE)/mylibs/libAliL3Hough");
+
+  gROOT->LoadMacro("$(ALICE)/dev/XFunct.C");  
+  gStyle->SetStatBorderSize(1);
+  gStyle->SetTitleBorderSize(0);
+  printf("TPC libraries loaded\n");
+  printf("Level3 libraries loaded\n");
+}