Checking in latest changes.
authorvestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 10 Dec 2001 16:33:59 +0000 (16:33 +0000)
committervestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 10 Dec 2001 16:33:59 +0000 (16:33 +0000)
23 files changed:
HLT/hough/AliL3Defs.h
HLT/hough/AliL3Histogram.cxx
HLT/hough/AliL3Histogram.h
HLT/hough/AliL3Histogram1D.cxx
HLT/hough/AliL3Histogram1D.h
HLT/hough/AliL3Hough.cxx
HLT/hough/AliL3Hough.h
HLT/hough/AliL3HoughDisplay.cxx
HLT/hough/AliL3HoughDisplay.h
HLT/hough/AliL3HoughEval.cxx
HLT/hough/AliL3HoughEval.h
HLT/hough/AliL3HoughLinkDef.h
HLT/hough/AliL3HoughMaxFinder.cxx
HLT/hough/AliL3HoughMaxFinder.h
HLT/hough/AliL3HoughMerger.cxx
HLT/hough/AliL3HoughMerger.h
HLT/hough/AliL3HoughTrack.cxx
HLT/hough/AliL3HoughTrack.h
HLT/hough/AliL3HoughTransformer.cxx
HLT/hough/AliL3HoughTransformer.h
HLT/hough/GetGoodParticles.cxx
HLT/hough/GetGoodParticles.h
HLT/hough/test.C

index e533f2e..0390d12 100644 (file)
@@ -4,13 +4,14 @@
 #include <Rtypes.h>
 
 const Int_t NPatches = 6;
-const Int_t NRowsSlice = 175;
+const Int_t NRowsSlice = 176;
 const Int_t NRows[6][2] = {{0,31},{32,63},{64,91},{92,119},{120,143},{144,175}};
-//const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
+//const Int_t NRows[6][2] = {{ 0, 175},{0,0},{0,0},{0,0},{0,0},{0,0}};
 const Double_t Pi = 3.14159265358979323846;
 const Double_t ToRad = Pi/180.;
 const Int_t MaxNPads = 256;
 const Int_t MaxNTimeBins = 512;
+const Double_t BField = 0.2;
 
 /*fRow[0][0] = 0;     // first row
   fRow[0][1] = 31;
index bbea051..5bb8af5 100644 (file)
@@ -1,10 +1,13 @@
-//Author:        Anders Strand Vestbo
-//Last Modified: 28.6.01
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV
 
 #include "AliL3Logging.h"
 #include "AliL3Histogram.h"
 
-//2D histogram class.
+//_____________________________________________________________
+// AliL3Histogram
+//
+// 2D histogram class
 
 ClassImp(AliL3Histogram)
 
index c43e5a9..77a1eb4 100644 (file)
@@ -69,7 +69,7 @@ class AliL3Histogram {
   Int_t GetNEntries() {return fEntries;}
   
   
-  ClassDef(AliL3Histogram,1)
+  ClassDef(AliL3Histogram,1) //2D histogram class
     
 };
 
index 4911724..24c27a9 100644 (file)
@@ -1,10 +1,13 @@
-//Author:        Anders Strand Vestbo
-//Last Modified: 28.6.01
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV
 
 #include "AliL3Logging.h"
 #include "AliL3Histogram1D.h"
 
-//2D histogram class.
+//_____________________________________________________________
+// AliL3Histogram1D
+//
+// 1D histogram class.
 
 ClassImp(AliL3Histogram1D)
 
@@ -76,6 +79,21 @@ Int_t AliL3Histogram1D::FindBin(Double_t x)
 
 }
 
+Int_t AliL3Histogram1D::GetMaximumBin()
+{
+  Double_t max_value=0;
+  Int_t max_bin=0;
+  for(Int_t i=0; i<fNcells; i++)
+    {
+      if(fContent[i] > max_value)
+       {
+         max_value=fContent[i];
+         max_bin = i;
+       }
+    }
+  return max_bin;
+}
+
 Double_t AliL3Histogram1D::GetBinContent(Int_t bin)
 {
   if(bin >= fNcells)
index f8d3e02..307fc96 100644 (file)
@@ -32,6 +32,7 @@ class AliL3Histogram1D {
   void Reset();
   void Fill(Double_t x,Int_t weight=1);
   void AddBinContent(Int_t bin,Int_t weight);
+  Int_t GetMaximumBin();
   Int_t FindBin(Double_t x);
   Double_t GetBinContent(Int_t bin);
   Double_t GetBinCenter(Int_t bin);
@@ -46,7 +47,7 @@ class AliL3Histogram1D {
   TH1F *GetRootHisto() {return fRootHisto;}
 #endif
   
-  ClassDef(AliL3Histogram1D,1)
+  ClassDef(AliL3Histogram1D,1) //1D histogram class
     
 };
 
index a4d2298..761454f 100644 (file)
@@ -1,10 +1,16 @@
-//Author:        Anders Strand Vestbo
-//Last Modified: 28.6.01
+//$Id$
+
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV 
+
 
 #include <string.h>
 #include <TCanvas.h>
 #include <TFile.h>
 
+#include "AliL3HoughMerger.h"
+#include "AliL3HoughIntMerger.h"
+#include "AliL3HoughGlobalMerger.h"
 #include "AliL3Logging.h"
 #include "AliL3Histogram.h"
 #include "AliL3Hough.h"
 #include "AliL3Defs.h"
 #include "AliL3TrackArray.h"
 #include "AliL3HoughTrack.h"
+#include "AliL3Benchmark.h"
+
+//_____________________________________________________________
+// AliL3Hough
+//
+// Base class for the Hough transform
+//
+
 
 ClassImp(AliL3Hough)
 
 AliL3Hough::AliL3Hough()
 {
+  //Constructor
+  
   fBinary = kFALSE;
   fNEtaSegments = 0;
   fAddHistograms = kFALSE;
-  fRemoveFoundTracks = kFALSE; 
+  fDoIterative = kFALSE; 
   fWriteDigits=kFALSE;
 }
 
@@ -36,7 +52,7 @@ AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments)
   strcpy(fPath,path);
   fNEtaSegments = n_eta_segments;
   fAddHistograms = kFALSE;
-  fRemoveFoundTracks = kFALSE; 
+  fDoIterative = kFALSE; 
   fWriteDigits = kFALSE;
   Init();
 }
@@ -44,92 +60,81 @@ AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments)
 
 AliL3Hough::~AliL3Hough()
 {
-  if(fMemHandler)
-    DeleteMemory();
-  if(fHoughTransformer)
-    DeleteTransformers();
-  if(fEval)
-    DeleteEval();
+  CleanUp();
+  if(fMerger)
+    delete fMerger;
+  if(fInterMerger)
+    delete fInterMerger;
   if(fPeakFinder)
     delete fPeakFinder;
-  if(fTracks)
-    delete fTracks;
-  if(fRootFile)
-    {
-      fRootFile->Close();
-      delete fRootFile;
-    }
-}
-
-void AliL3Hough::DeleteEval()
-{
-  for(Int_t i=0; i<NPatches; i++)
-    {
-      if(!fEval[i]) continue;
-      delete fEval[i];
-    }
-  delete [] fEval;
-}
-
-void AliL3Hough::DeleteTransformers()
-{
-  for(Int_t i=0; i<NPatches; i++)
-    {
-      if(!fHoughTransformer[i]) continue;
-      delete fHoughTransformer[i];
-    }
-  delete [] fHoughTransformer;
 }
 
-void AliL3Hough::DeleteMemory()
+void AliL3Hough::CleanUp()
 {
+  //Cleanup memory
+  
   for(Int_t i=0; i<NPatches; i++)
     {
-      if(!fMemHandler[i]) continue;
-      delete fMemHandler[i];
+      if(fTracks[i]) delete fTracks[i];
+      if(fEval[i]) delete fEval[i];
+      if(fHoughTransformer[i]) delete fHoughTransformer[i];
+      if(fMemHandler[i]) delete fMemHandler[i];
     }
-  delete [] fMemHandler;
+  
+  /*Shitty compiler doesn't allow this:
+    
+  if(fTracks) delete [] fTracks;
+  if(fEval) delete [] fEval;
+  if(fHoughTransformer) delete [] fHoughTransformer;
+  if(fMemHandler) delete [] fMemHandler;
+  */
 }
 
 void AliL3Hough::Init()
 {
   fHoughTransformer = new AliL3HoughTransformer*[NPatches];
   fMemHandler = new AliL3FileHandler*[NPatches];
+  fTracks = new AliL3TrackArray*[NPatches];
+  fEval = new AliL3HoughEval*[NPatches];
   for(Int_t i=0; i<NPatches; i++)
     {
       fHoughTransformer[i] = new AliL3HoughTransformer(1,i,fNEtaSegments);
       fHoughTransformer[i]->CreateHistograms(64,-0.003,0.003,64,-0.26,0.26);
       fHoughTransformer[i]->SetThreshold(3);
+      fEval[i] = new AliL3HoughEval();
+      fTracks[i] = new AliL3TrackArray("AliL3HoughTrack");
       fMemHandler[i] = new AliL3FileHandler();
       if(!fBinary)
        fMemHandler[i]->SetAliInput(fPath);
     }
   fPeakFinder = new AliL3HoughMaxFinder("KappaPhi");
+  fMerger = new AliL3HoughMerger(NPatches);
+  fInterMerger = new AliL3HoughIntMerger();
 }
 
 void AliL3Hough::Process(Int_t minslice,Int_t maxslice)
 {
   //Process all slices [minslice,maxslice].
-
+  
   for(Int_t i=minslice; i<=maxslice; i++)
     {
-      TransformSlice(i);
+      ReadData(i);
+      Transform();
       if(fAddHistograms)
        AddAllHistograms();
       FindTrackCandidates();
-      Evaluate(fRemoveFoundTracks);
+      Evaluate();
       if(fWriteDigits)
        WriteDigits();
     }
 }
 
-void AliL3Hough::TransformSlice(Int_t slice)
+void AliL3Hough::ReadData(Int_t slice)
 {
-  
+  //Read data from files, binary or root.
+
   for(Int_t i=0; i<NPatches; i++)
     {
-      //Reset memories
-      fHoughTransformer[i]->Reset();      
       fMemHandler[i]->Free();
       UInt_t ndigits=0;
       AliL3DigitRowData *digits =0;
@@ -147,11 +152,97 @@ void AliL3Hough::TransformSlice(Int_t slice)
          digits=(AliL3DigitRowData *)fMemHandler[i]->AliDigits2Memory(ndigits); 
        }
       fHoughTransformer[i]->SetInputData(ndigits,digits);
+    }
+}
+
+void AliL3Hough::Transform()
+{
+  //Transform all data given to the transformer within the given slice
+  //(after ReadData(slice))
+
+  Double_t initTime,finalTime;
+  for(Int_t i=0; i<NPatches; i++)
+    {
+      fHoughTransformer[i]->Reset();//Reset the histograms
+      initTime = AliL3Benchmark::GetCpuTime();
       fHoughTransformer[i]->TransformCircle();
+      finalTime = AliL3Benchmark::GetCpuTime();
+      LOG(AliL3Log::kInformational,"AliL3Hough::Transform","Timing")
+       <<AliL3Log::kDec<<"Transform finished in "<<(finalTime-initTime)*1000<<"ms"<<ENDLOG;
+    }
+}
+
+void AliL3Hough::MergePatches()
+{
+  if(fAddHistograms) //Nothing to merge here
+    return;
+  AliL3Transform *tr = new AliL3Transform();
+  fMerger->SetTransformer(tr);
+  fMerger->MergePatches(kTRUE);
+  delete tr;
+}
+
+void AliL3Hough::MergeInternally()
+{
+  if(fAddHistograms)
+    fInterMerger->FillTracks(fTracks[0]);
+  else
+    fInterMerger->FillTracks(fMerger->GetOutTracks());
+  
+  fInterMerger->MMerge();
+}
+
+void AliL3Hough::ProcessSliceIter()
+{
+  //Process current slice (after ReadData(slice)) iteratively.
+  
+  for(Int_t i=0; i<NPatches; i++)
+    {
+      ProcessPatchIter(i);
+      fMerger->FillTracks(fTracks[i],i); //Copy tracks to merger
     }
   
 }
 
+void AliL3Hough::ProcessPatchIter(Int_t patch)
+{
+  //Process patch in a iterative way. 
+  //transform + peakfinding + evaluation + transform +...
+
+  Int_t num_of_tries = 10;
+  AliL3HoughTransformer *tr = fHoughTransformer[patch];
+  AliL3TrackArray *tracks = fTracks[patch];
+  tracks->Reset();
+  AliL3HoughEval *ev = fEval[patch];
+  ev->InitTransformer(tr);
+  ev->RemoveFoundTracks();
+  ev->SetNumOfRowsToMiss(2);
+  AliL3Histogram *hist;
+  for(Int_t t=0; t<num_of_tries; t++)
+    {
+      tr->Reset();
+      tr->TransformCircle();
+      for(Int_t i=0; i<fNEtaSegments; i++)
+       {
+         hist = tr->GetHistogram(i);
+         if(hist->GetNEntries()==0) continue;
+         fPeakFinder->SetHistogram(hist);
+         Int_t n=1;
+         Int_t x[n],y[n];
+         fPeakFinder->FindAbsMaxima(*x,*y);
+         AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
+         track->SetTrackParameters(hist->GetBinCenterX(*x),hist->GetBinCenterY(*y),1);
+         if(!ev->LookInsideRoad(track,i))
+           {   
+             tracks->Remove(tracks->GetNTracks()-1);
+             tracks->Compress();
+           }
+       }
+    }
+  LOG(AliL3Log::kInformational,"AliL3Hough::ProcessPatch","NTracks")
+    <<AliL3Log::kDec<<"Found "<<tracks->GetNTracks()<<" tracks in patch "<<patch<<ENDLOG;
+}
+
 AliL3Histogram *AliL3Hough::AddHistograms(Int_t eta_index)
 {
 
@@ -169,9 +260,7 @@ void AliL3Hough::AddAllHistograms()
 {
   //Add the histograms within one etaslice.
   //Resulting histogram are in patch=0.
-  
-  LOG(AliL3Log::kDebug,"AliL3Hough::AddAllHistograms","Progress")
-    <<"Adding all histograms"<<ENDLOG;
+
   for(Int_t i=0; i<fNEtaSegments; i++)
     {
       AliL3Histogram *hist0 = fHoughTransformer[0]->GetHistogram(i);
@@ -181,90 +270,111 @@ void AliL3Hough::AddAllHistograms()
          hist0->Add(hist);
        }
     }
+  fAddHistograms = kTRUE;
 }
 
 void AliL3Hough::FindTrackCandidates()
 {
   //Look for peaks in histograms, and find the track candidates
   
-  if(fTracks)
-  {
-    LOG(AliL3Log::kDebug,"AliL3Hough::FindTrackCandidates","Track array")
-      <<"Deleting old track array"<<ENDLOG;
-    delete fTracks;
-  }
-  fTracks = new AliL3TrackArray("AliL3HoughTrack");
-  
   Int_t n_patches;
   if(fAddHistograms)
     n_patches = 1; //Histograms has been added.
   else
     n_patches = NPatches;
-
+  
   for(Int_t i=0; i<n_patches; i++)
     {
       AliL3HoughTransformer *tr = fHoughTransformer[i];
+      fTracks[i]->Reset();
       for(Int_t j=0; j<fNEtaSegments; j++)
        {
          AliL3Histogram *hist = tr->GetHistogram(j);
+         if(hist->GetNEntries()==0) continue;
          fPeakFinder->SetHistogram(hist);
          Int_t n=10;
          Float_t x[10];
          Float_t y[10];
-         fPeakFinder->FindPeak1(x,y,n);
+         Int_t weight[10];
+         fPeakFinder->FindPeak1(x,y,weight,n);
          for(Int_t k=0; k<n; k++)
            {
-             AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks->NextTrack();
-             track->SetTrackParameters(x[k],y[k],1);
+             if(weight[k] == 0) continue;
+             AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
+             track->SetTrackParameters(x[k],y[k],weight[k]);
              track->SetEtaIndex(j);
+             track->SetEta((Double_t)(j*tr->GetEtaSlice()));
+             track->SetRowRange(NRows[0][0],NRows[5][1]);
            }
        }
+      fTracks[i]->QSort();
     }
-  
 }
 
-void AliL3Hough::Evaluate(Bool_t remove)
+void AliL3Hough::Evaluate(Int_t road_width)
 {
   //Evaluate the tracks, by looking along the road in the raw data.
-  //You may choose to remove the found tracks from the image.
   
-  if(!fTracks)
+
+  if(!fTracks[0])
     {
-      LOG(AliL3Log::kError,"AliL3Hough::Evaluate","Track array")
-       <<AliL3Log::kHex<<"No tracks to work on "<<(Int_t)fTracks<<ENDLOG;
+      LOG(AliL3Log::kError,"AliL3Hough::Evaluate","Track Array")
+       <<"No tracks to work with..."<<ENDLOG;
       return;
     }
   
-  if(fEval)
-    {
-      LOG(AliL3Log::kDebug,"AliL3Hough::Evaluate","Evaluate object")
-       <<"Deleting old AliL3HoughEval objects"<<ENDLOG;
-      DeleteEval();
-    }
-  
-  fEval = new AliL3HoughEval*[NPatches];
+  printf("Number of tracks before evaluation %d\n",fTracks[0]->GetNTracks());
+  AliL3TrackArray *tracks;
   for(Int_t i=0; i<NPatches; i++)
     {
-      fEval[i] = new AliL3HoughEval(fHoughTransformer[i]);
-      if(remove)
-       fEval[i]->RemoveFoundTracks();
-      for(Int_t j=0; j<fTracks->GetNTracks(); j++)
+      fEval[i]->InitTransformer(fHoughTransformer[i]);
+      fEval[i]->SetNumOfRowsToMiss(1);
+      fEval[i]->SetNumOfPadsToLook(road_width);
+      if(fAddHistograms)
+       tracks = fTracks[0];
+      else
+       tracks = fTracks[i];
+      for(Int_t j=0; j<tracks->GetNTracks(); j++)
        {
-         AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks->GetCheckedTrack(j);
+         AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
          if(!track)
            {
-             printf("AliL3Hough::Evaluate : Missing track object...\n");
+             LOG(AliL3Log::kWarning,"AliL3Hough::Evaluate","Track array")
+               <<"Track object missing!"<<ENDLOG;
              continue;
            }
-         //fEval[i]->LookInsideRoad(track,track->GetEtaIndex());
          if(!fEval[i]->LookInsideRoad(track,track->GetEtaIndex()))
-           fTracks->Remove(j);
+           tracks->Remove(j);
+         if(fAddHistograms)
+           track->SetRowRange(NRows[0][0],NRows[5][1]);//All rows included
        }
-      fTracks->Compress();
+      tracks->Compress();
+      tracks->QSort(); //Sort the tracks according to weight
+      
+      if(!fAddHistograms)
+       fMerger->FillTracks(tracks,i); //Copy tracks to the track merger
     }
   
 }
 
+void AliL3Hough::EvaluateWithEta()
+{
+  if(!fTracks[0])
+    {
+      printf("AliL3Hough::EvaluateWithEta: NO TRACKS\n");
+      return;
+    }
+  printf("Number of tracks before evaluation %d\n",fTracks[0]->GetNTracks());
+  for(Int_t i=0; i<NPatches; i++)
+    {
+      fEval[i]->InitTransformer(fHoughTransformer[i]);
+      fEval[i]->FindEta(fTracks[0]);
+    }
+  fMerger->FillTracks(fTracks[0],0);
+  fMerger->MergeEtaSlices(0);
+}
+
 void AliL3Hough::WriteDigits(Char_t *outfile)
 {
   //Write the current data to a new rootfile.
index 3f0208d..9b5d7f6 100644 (file)
@@ -10,27 +10,27 @@ class AliL3FileHandler;
 class AliL3HoughEval;
 class AliL3Transform;
 class AliL3TrackArray;
-class TFile;
+class AliL3HoughMerger;
+class AliL3HoughIntMerger;
 
-class AliL3Hough : public TObject {
+class AliL3Hough {
   
  private:
   Char_t fPath[256];
   Bool_t fBinary;
   Bool_t fAddHistograms;
-  Bool_t fRemoveFoundTracks;
+  Bool_t fDoIterative;
   Bool_t fWriteDigits;
   Int_t fNEtaSegments;
   AliL3FileHandler **fMemHandler; //!
   AliL3HoughTransformer **fHoughTransformer; //!
   AliL3HoughEval **fEval; //!
   AliL3HoughMaxFinder *fPeakFinder; //!
-  AliL3TrackArray *fTracks; //!
-  TFile *fRootFile; //!
+  AliL3TrackArray **fTracks; //!
+  AliL3HoughMerger *fMerger; //!
+  AliL3HoughIntMerger *fInterMerger; //!
   
-  void DeleteEval();
-  void DeleteTransformers();
-  void DeleteMemory();
+  void CleanUp();
   void Init();
   
  public:
@@ -40,24 +40,35 @@ class AliL3Hough : public TObject {
   virtual ~AliL3Hough();
   
   void Process(Int_t minslice,Int_t maxslice);
-  void TransformSlice(Int_t slice);
+  void ReadData(Int_t slice);
+  void Transform();
+  void ProcessSliceIter();
+  void ProcessPatchIter(Int_t patch);
+  void MergePatches();
+  void MergeInternally();
+
   void FindTrackCandidates();
   AliL3Histogram *AddHistograms(Int_t eta_index);
   void AddAllHistograms();
-  void Evaluate(Bool_t remove=kFALSE);
+  void Evaluate(Int_t road_width=1);
+  void EvaluateWithEta();
   void WriteDigits(Char_t *outfile="output_digits.root");
   
   //Setters
   void SetNEtaSegments(Int_t i) {fNEtaSegments = i;}
   void SetAddHistograms() {fAddHistograms = kTRUE;}
-  void SetRemoveFoundTracks() {fRemoveFoundTracks = kTRUE;}
+  void DoIterative() {fDoIterative = kTRUE;}
   void SetWriteDigits() {fWriteDigits = kTRUE;}
   
   //Getters
   AliL3HoughTransformer *GetTransformer(Int_t i) {if(!fHoughTransformer[i]) return 0; return fHoughTransformer[i];}
-  AliL3TrackArray *GetTracks() {if(!fTracks) return 0; return fTracks;}
+  AliL3TrackArray *GetTracks(Int_t i) {if(!fTracks[i]) return 0; return fTracks[i];}
+  AliL3HoughEval *GetEval(Int_t i) {if(!fEval[i]) return 0; return fEval[i];}
+  AliL3HoughMerger *GetMerger() {if(!fMerger) return 0; return fMerger;}
+  AliL3HoughIntMerger *GetInterMerger() {if(!fInterMerger) return 0; return fInterMerger;}
+  AliL3FileHandler *GetMemHandler(Int_t i) {if(!fMemHandler[i]) return 0; return fMemHandler[i];}
 
-  ClassDef(AliL3Hough,1)
+  ClassDef(AliL3Hough,1) //Hough transform base class
 
 };
 
index 14d0aab..a239115 100644 (file)
@@ -1,6 +1,7 @@
-//Author:        Anders Strand Vestbo
-//Last Modified: 12.01.2001
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV 
 
+#include <iostream.h>
 #include <TCanvas.h>
 #include <TView.h>
 #include <TPolyMarker3D.h>
 #include "AliL3Transform.h"
 #include "AliL3HoughTrack.h"
 #include "AliL3TrackArray.h"
-#include "AliL3Logging.h"
 #include "AliL3HoughDisplay.h"
 #include "AliL3Defs.h"
+#include "AliL3MemHandler.h"
+
+//_____________________________________________________________
+// Display class for Hough transform code
 
 ClassImp(AliL3HoughDisplay)
 
 
 AliL3HoughDisplay::AliL3HoughDisplay()
 {
-  //Ctor. Specify which slices you want to look at.
 
   fTracks = 0;
-
+  fDigitRowData = 0;
+  fNDigitRowData = 0;
+  fShowSlice = -1;
+  fPatch = -1;
   Init();
 }
 
@@ -39,8 +45,7 @@ void AliL3HoughDisplay::Init()
 {
   TFile *file = TFile::Open("/prog/alice/data/GEO/alice.geom");
   if(!file->IsOpen())
-    LOG(AliL3Log::kError,"AliL3HoughDisplay::AliL3HoughDisplay","File Open")
-      <<"Geometry file alice.geom does not exist"<<ENDLOG;
+    cerr<<"AliL3HoughDisplay::AliL3HoughDisplay : Geometry file alice.geom does not exist"<<endl;
   fGeom = (TGeometry*)file->Get("AliceGeom");
   file->Close();
   fTransform = new AliL3Transform();
@@ -67,7 +72,47 @@ void AliL3HoughDisplay::GenerateHits(AliL3HoughTrack *track,Float_t *x,Float_t *
     }
 }
 
-void AliL3HoughDisplay::DisplayTracks()
+TPolyMarker3D *AliL3HoughDisplay::LoadDigits()
+{
+  
+  AliL3DigitRowData *tempPt = fDigitRowData;
+  if(!tempPt)
+    {
+      cerr<<"AliL3HoughDisplay::LoadDigits : No data"<<endl;
+      return 0;
+    }
+  
+  Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
+  Int_t count=0;
+  for(UInt_t i=0; i<nrows; i++)
+    {
+      count += tempPt->fNDigit;
+      AliL3MemHandler::UpdateRowPointer(tempPt);
+    }
+  tempPt = fDigitRowData;
+  TPolyMarker3D *pm = new TPolyMarker3D(count);
+  Float_t xyz[3];
+  Int_t sector,row;
+  count=0;
+  for(UInt_t i=0; i<nrows; i++)
+    {
+      AliL3DigitData *digPt = tempPt->fDigitData;
+      Int_t padrow = (Int_t)tempPt->fRow;
+      for(UInt_t j=0; j<tempPt->fNDigit; j++)
+       {
+         fTransform->Slice2Sector(fShowSlice,padrow,sector,row);
+         fTransform->Raw2Global(xyz,sector,row,(Int_t)digPt->fPad,(Int_t)digPt->fTime);
+         pm->SetPoint(count,xyz[0],xyz[1],xyz[2]);
+         count++;
+       }
+      AliL3MemHandler::UpdateRowPointer(tempPt);
+    }
+
+  cout<<"Displaying "<<count<<" digits"<<endl;
+  return pm;
+}
+
+void AliL3HoughDisplay::DisplayEvent()
 {
   //Display the found tracks.
   
@@ -82,7 +127,7 @@ void AliL3HoughDisplay::DisplayTracks()
   
   TView *v = new TView(1);
   v->SetRange(-430,-560,-430,430,560,1710);
-  
+
   c1->Clear();
   c1->SetFillColor(1);
   c1->SetTheta(90.);
@@ -102,7 +147,7 @@ void AliL3HoughDisplay::DisplayTracks()
       for(Int_t h=0; h<n; h++)
        pm->SetPoint(h,x[h],y[h],z[h]);
       
-      pm->SetMarkerColor(2);
+      pm->SetMarkerColor(1);
       pm->Draw();
       TPolyLine3D *current_line = &(line[j]);
       current_line = new TPolyLine3D(n,x,y,z,"");
@@ -111,6 +156,13 @@ void AliL3HoughDisplay::DisplayTracks()
       
     }
   
+  if(fShowSlice>=0)
+    {
+      TPolyMarker3D *pm = LoadDigits();
+      pm->SetMarkerColor(2);
+      pm->Draw("same");
+    }
+  
   fGeom->Draw("same");
   c1->x3d();
 }
index b2e3d2d..21f6ed1 100644 (file)
@@ -6,6 +6,9 @@
 class TGeometry;
 class AliL3TrackArray;
 class AliL3Transform;
+class AliL3DigitRowData;
+class TPolyMarker3D;
+class AliL3HoughTrack;
 
 class AliL3HoughDisplay {
 
@@ -14,18 +17,33 @@ class AliL3HoughDisplay {
   TGeometry *fGeom; //!
   AliL3TrackArray *fTracks; //!
   AliL3Transform *fTransform; //!
-
+  AliL3DigitRowData *fDigitRowData;  //!
+  UInt_t fNDigitRowData; //!
+  Int_t fShowSlice; 
+  Int_t fPatch;
+  
   void GenerateHits(AliL3HoughTrack *track,Float_t *x,Float_t *y,Float_t *z,Int_t &n);
   void Init();
-  
+  TPolyMarker3D *LoadDigits();
+
  public:
   AliL3HoughDisplay();
   virtual ~AliL3HoughDisplay();
 
-  void DisplayTracks();
+  void DisplayEvent();
   void SetTracks(AliL3TrackArray *tracks) {fTracks = tracks;}
+  void ShowData(AliL3DigitRowData *data,UInt_t size,Int_t slice,Int_t patch);
 
   ClassDef(AliL3HoughDisplay,1) 
 };
 
+inline void AliL3HoughDisplay::ShowData(AliL3DigitRowData *data,UInt_t size,Int_t slice,Int_t patch)
+{
+  fShowSlice = slice;
+  fPatch = patch;
+  fDigitRowData = data;
+  fNDigitRowData = size;
+}
+
+
 #endif
index 7e145d8..33623b4 100644 (file)
@@ -1,5 +1,5 @@
-//Author:        Anders Strand Vestbo
-//Last Modified: 28.6.01
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV 
 
 #include <math.h>
 #include <TH1.h>
 #include "AliL3Histogram1D.h"
 #include "AliL3Defs.h"
 
+//_____________________________________________________________
+// AliL3HoughEval
+//
+// Evaluation class for tracklets produced by the Hough transform.
+
 ClassImp(AliL3HoughEval)
 
 AliL3HoughEval::AliL3HoughEval()
@@ -28,7 +33,7 @@ AliL3HoughEval::AliL3HoughEval()
   fNumOfPadsToLook = 1;
   fNumOfRowsToMiss = 1;
   fEtaHistos=0;
-  
+  fRowPointers = 0;
 }
 
 
@@ -38,7 +43,11 @@ AliL3HoughEval::~AliL3HoughEval()
   if(fTransform)
     delete fTransform;
   if(fRowPointers)
-    delete [] fRowPointers;
+    {
+      for(Int_t i=0; i<fNrows; i++)
+       fRowPointers[i] = 0;
+      delete [] fRowPointers;
+    }
 }
 
 void AliL3HoughEval::InitTransformer(AliL3HoughTransformer *transformer)
@@ -57,7 +66,8 @@ void AliL3HoughEval::GenerateLUT()
 {
   //Generate a Look-up table, to limit the access to raw data
   
-  fRowPointers = new AliL3DigitRowData*[fNrows];
+  if(!fRowPointers)
+    fRowPointers = new AliL3DigitRowData*[fNrows];
 
   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer->GetDataPointer();
   if(!tempPt)
@@ -90,7 +100,7 @@ Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Boo
       
       if(!track->GetCrossingPoint(padrow,xyz))  
        {
-         printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
+         printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!; pt %f phi0 %f\n",track->GetPt(),track->GetPhi0());
          continue;
        }
       
@@ -112,7 +122,7 @@ Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Boo
          AliL3DigitData *digPt = tempPt->fDigitData;
          for(UInt_t j=0; j<tempPt->fNDigit; j++)
            {
-             if(digPt->fCharge <= fHoughTransformer->GetThreshold()) continue;
+             //if(digPt->fCharge <= fHoughTransformer->GetThreshold()) continue;
              UChar_t pad = digPt[j].fPad;
              if(pad < p) continue;
              if(pad > p) break;
@@ -142,6 +152,8 @@ Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Boo
       track->SetEtaIndex(eta_index);
       track->SetWeight(total_charge,kTRUE);
       track->SetEta(eta_track);
+      track->SetRowRange(NRows[fPatch][0],NRows[fPatch][1]);
+      track->SetSlice(fSlice);
       if(fRemoveFoundTracks)
        LookInsideRoad(track,eta_index,kTRUE);
       return kTRUE;
@@ -286,7 +298,7 @@ void AliL3HoughEval::DisplayEtaSlice(Int_t eta_index,AliL3Histogram *hist)
   
 }
 
-void AliL3HoughEval::CompareMC(AliL3TrackArray *tracks,Char_t *trackfile)
+void AliL3HoughEval::CompareMC(AliL3TrackArray *tracks,Char_t *trackfile,Int_t threshold)
 {
   
   struct GoodTrack goodtracks[15000];
@@ -339,7 +351,7 @@ void AliL3HoughEval::CompareMC(AliL3TrackArray *tracks,Char_t *trackfile)
     {
       AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
       if(!tr) continue;
-      //if(tr->GetWeight()<14000) continue;
+      if(tr->GetWeight()<threshold) continue;
       Int_t trackindex = tr->GetEtaIndex();
       if(trackindex <0 || trackindex >= fNEtaSegments) continue;
       ftracks[trackindex]++;
index f0ce488..e437981 100644 (file)
@@ -40,7 +40,7 @@ class AliL3HoughEval : public TObject {
   void GenerateLUT();
   void DisplayEtaSlice(Int_t eta_index,AliL3Histogram *hist);
   Bool_t LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove=kFALSE);
-  void CompareMC(AliL3TrackArray *tracks,Char_t *goodtracks="good_tracks");
+  void CompareMC(AliL3TrackArray *tracks,Char_t *goodtracks="good_tracks",Int_t treshold=0);
   void FindEta(AliL3TrackArray *tracks);
   
   //Getters
@@ -51,7 +51,7 @@ class AliL3HoughEval : public TObject {
   void SetNumOfRowsToMiss(Int_t i) {fNumOfRowsToMiss = i;}
   void SetNumOfPadsToLook(Int_t i) {fNumOfPadsToLook = i;}
 
-  ClassDef(AliL3HoughEval,1)
+  ClassDef(AliL3HoughEval,1) //Hough transform verfication class
 
 };
 
index 30688ec..79cd6f6 100644 (file)
@@ -13,6 +13,9 @@
 #pragma link C++ class AliL3Histogram;
 #pragma link C++ class AliL3Histogram1D;
 #pragma link C++ class AliL3ClusterFinder;
+#pragma link C++ class AliL3HoughIntMerger;
+#pragma link C++ class AliL3HoughGlobalMerger;
+#pragma link C++ class AliL3HoughDisplay;
 #pragma link C++ function GetGoodParticles(Int_t,Int_t,char *,char *);
 
 #endif
index c30bdd3..1b6472b 100644 (file)
@@ -1,5 +1,5 @@
-//Author:        Anders Strand Vestbo
-//Last Modified: 28.6.01
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV 
 
 #include <string.h>
 #include <math.h>
 #include "AliL3HoughTrack.h"
 #include "AliL3HoughMaxFinder.h"
 
+//_____________________________________________________________
+// AliL3HoughMaxFinder
+//
+// Maximum finder
+
 ClassImp(AliL3HoughMaxFinder)
 
   
@@ -496,10 +501,16 @@ AliL3HoughTrack *AliL3HoughMaxFinder::FindPeakLine(Double_t rho,Double_t theta)
 
 
 
-void AliL3HoughMaxFinder::FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t &n)
+void AliL3HoughMaxFinder::FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t *weight,Int_t &n,
+                                   Int_t y_window,Int_t x_bin_sides)
 {
-  //testing mutliple peakfinding
-  
+  //Testing mutliple peakfinding.
+  //The algorithm searches the histogram for prepreaks by looking in windows
+  //for each bin on the xaxis. The size of these windows is controlled by y_window.
+  //Then the prepreaks are sorted according to their weight (sum inside window),
+  //and the peak positions are calculated by taking the weighted mean in both
+  //x and y direction. The size of the peak in x-direction is controlled by x_bin_sides.
+
   Int_t max_entries = n;
   n=0;
   if(!fCurrentHisto)
@@ -507,8 +518,12 @@ void AliL3HoughMaxFinder::FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t &n)
       printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
       return;
     }  
-  Int_t y_window=2;
-  Int_t x_bin_sides=1;
+  //Int_t y_window=2;
+  //Int_t x_bin_sides=1;
+  
+  Float_t max_kappa = 0.001;
+  Float_t max_phi0 = 0.05;
+  
   Int_t max_sum=0;
   
   Int_t xmin = fCurrentHisto->GetFirstXbin();
@@ -565,9 +580,9 @@ void AliL3HoughMaxFinder::FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t &n)
       
       if(xbin<xmin || xbin > xmax-1) continue;
       
-      //Check if this is really a peak
-      if(anotherPt[xbin-1]->weight > windowPt[i]->weight ||
-        anotherPt[xbin+1]->weight > windowPt[i]->weight)
+      //Check if this is really a local maxima
+      if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight ||
+        anotherPt[xbin+1]->weight > anotherPt[xbin]->weight)
        continue;
 
       for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++)
@@ -586,27 +601,50 @@ void AliL3HoughMaxFinder::FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t &n)
   
   //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<n; i++)
     {
       Int_t xbin = fCurrentHisto->FindXbin(xpeaks[i]);
       if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue;
       top=butt=0;
       ytop=0,ybutt=0;    
+      w=0;
+      prev = xbin - x_bin_sides+1;
       for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++)
        {
-         if(anotherPt[j]->ymin > anotherPt[xbin]->ymax) continue;
-         if(anotherPt[j]->ymax < anotherPt[xbin]->ymin) continue;
+         //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]->weight;
          butt += anotherPt[j]->weight;
+         
          for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++)
            {
              Int_t bin = fCurrentHisto->GetBin(j,k);
              ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
              ybutt += fCurrentHisto->GetBinContent(bin);
+             w+=fCurrentHisto->GetBinContent(bin);
            }
        }
+      
       xpeaks[i] = top/butt;
       ypeaks[i] = ytop/ybutt;
+      weight[i] = w;
+      
+      //Check if this peak is overlapping with a previous:
+      for(Int_t p=0; p<i-1; p++)
+       {
+         if(fabs(xpeaks[p] - xpeaks[i]) < max_kappa ||
+            fabs(ypeaks[p] - ypeaks[i]) < max_phi0)
+           {
+             weight[i]=0;
+             break;
+           }
+       }
+      
     }
   
   for(Int_t i=0; i<nbinsx; i++)
@@ -666,7 +704,7 @@ Int_t AliL3HoughMaxFinder::PeakCompare(struct AxisWindow *a,struct AxisWindow *b
 
 }
 
-AliL3HoughTrack *AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
+void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3,Float_t &kappa,Float_t &phi0)
 {
   //Attempt of a more sophisticated peak finder.
   //Finds the best peak in the histogram, and returns the corresponding
@@ -675,7 +713,7 @@ AliL3HoughTrack *AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
   if(!fCurrentHisto)
     {
       printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
-      return 0;
+      return;
     }
   AliL3Histogram *hist = fCurrentHisto;
 
@@ -830,15 +868,16 @@ AliL3HoughTrack *AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
   bin = hist->FindBin(x_peak,y_peak);
   Int_t weight = (Int_t)hist->GetBinContent(bin);
 
-  AliL3HoughTrack *track = new AliL3HoughTrack();
-  track->SetTrackParameters(x_peak,y_peak,weight);
-  
+  //AliL3HoughTrack *track = new AliL3HoughTrack();
+  //track->SetTrackParameters(x_peak,y_peak,weight);
+  kappa = x_peak;
+  phi0 = y_peak;
   
   delete [] m;
   delete [] m_low;
   delete [] m_up;
   
-  return track;
+  //return track;
 
     
 }
index 229f918..15a0813 100644 (file)
@@ -34,10 +34,10 @@ class AliL3HoughMaxFinder : public TObject {
   AliL3TrackArray *FindMaxima(AliL3Histogram *hist,Int_t *rowrange=0,Int_t ref_row=0);
   AliL3TrackArray *LookForPeaks(AliL3Histogram *hist,Int_t nbins);
   
-  AliL3HoughTrack *FindPeak(Int_t t1,Double_t t2,Int_t t3);
+  void FindPeak(Int_t t1,Double_t t2,Int_t t3,Float_t &kappa,Float_t &phi0);
   AliL3HoughTrack *FindPeakLine(Double_t rho,Double_t theta);
   AliL3HoughTrack *CalculatePeakInWindow(Int_t *maxbin,Int_t t0,Int_t t1,Double_t t2,Int_t t3);
-  void FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t &n);
+  void FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t *weight,Int_t &n,Int_t y_window=2,Int_t x_bin_sides=1);
   void SortPeaks(struct AxisWindow **a,Int_t first,Int_t last);
   Int_t PeakCompare(struct AxisWindow *a,struct AxisWindow *b);
   
@@ -45,7 +45,7 @@ class AliL3HoughMaxFinder : public TObject {
   
   void SetHistogram(AliL3Histogram *hist) {fCurrentHisto = hist;}
   
-  ClassDef(AliL3HoughMaxFinder,1)
+  ClassDef(AliL3HoughMaxFinder,1) //Maximum finder class
 
 };
 
index 8cf7c42..dbabfc8 100644 (file)
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV 
 
 #include "AliL3Logging.h"
 #include "AliL3Defs.h"
+#include "AliL3Transform.h"
 #include "AliL3TrackArray.h"
 #include "AliL3HoughTrack.h"
 #include "AliL3HoughMerger.h"
 #include "AliL3HoughTransformer.h"
 
+//_____________________________________________________________
+// AliL3HoughMerger
+//
+// Patch merging class for Hough tracklets
+
 ClassImp(AliL3HoughMerger)
 
   
 AliL3HoughMerger::AliL3HoughMerger()
 {
   //Default constructor
-
 }
 
 
-AliL3HoughMerger::AliL3HoughMerger(Int_t nsubsectors)
+AliL3HoughMerger::AliL3HoughMerger(Int_t nsubsectors) : AliL3Merger(nsubsectors,"AliL3HoughTrack")
 {
   //Constructor
-  fNIn = nsubsectors;
-  SetArray();
+  Is2Global(kFALSE);
+  SetParameters(0.001,0.1,0.05);
 }
 
 
 AliL3HoughMerger::~AliL3HoughMerger()
 {
-  DeleteArray();
-  if(fOutTracks)
-    delete fOutTracks;
 }
 
-void AliL3HoughMerger::DeleteArray()
+void AliL3HoughMerger::FillTracks(AliL3TrackArray *tracks,Int_t patch)
 {
-  if(!fInTracks)
-    return;
-  for(Int_t i=0; i<fNIn; i++)
-    {
-      if(!fInTracks[i]) continue;
-      delete fInTracks;
-    }
-  //delete [] fInTracks;
+  if(tracks->GetNTracks()==0)
+    LOG(AliL3Log::kWarning,"AliL3HoughMerger::FillTracks","Track Array")
+      <<"Adding empty track array"<<ENDLOG;
+  
+  GetInTracks(patch)->AddTracks(tracks,kFALSE);//Copy tracks
+  printf("Filling %d tracks to merger\n",tracks->GetNTracks());
 }
 
-void AliL3HoughMerger::SetArray()
+void AliL3HoughMerger::SetParameters(Double_t maxkappa,Double_t maxpsi,Double_t maxphi0)
 {
-  fInTracks = new AliL3TrackArray*[fNIn];
-  for(Int_t i=0; i<fNIn; i++)
-    fInTracks[i] = new AliL3TrackArray("AliL3HoughTrack");
-  fOutTracks = new AliL3TrackArray("AliL3HoughTrack");
+  fMaxKappa = maxkappa;
+  fMaxPsi = maxpsi;
+  fMaxPhi0 = maxphi0;
 }
 
-void AliL3HoughMerger::FillTracks(AliL3TrackArray *tracks,Int_t patch)
+Bool_t AliL3HoughMerger::IsTrack(AliL3Track *innertrack,AliL3Track *outertrack)
 {
-  if(tracks->GetNTracks()==0)
-    LOG(AliL3Log::kWarning,"AliL3HoughMerger::FillTracks","Track Array")
-      <<"Adding empty track array"<<ENDLOG;
-  fInTracks[patch]->AddTracks(tracks,kFALSE); //Copy tracks
+  //Check if the tracks can be merged, called by the track merger
+  
+  AliL3HoughTrack *tr1 = (AliL3HoughTrack*)innertrack;
+  AliL3HoughTrack *tr2 = (AliL3HoughTrack*)outertrack;
+  
+  if( (!tr1->IsPoint()) || (!tr2->IsPoint()) )  return kFALSE; 
+  if(abs(tr1->GetEtaIndex() - tr2->GetEtaIndex()) > 1) return kFALSE;
+  if(tr1->GetCharge() != tr2->GetCharge()) return kFALSE;
+  if(fabs(tr1->GetPhi0() - tr2->GetPhi0()) > fMaxPhi0) return kFALSE;
+  if(fabs(tr1->GetKappa() - tr2->GetKappa()) > fMaxKappa) return kFALSE;
+  
+  /*
+    if( (!tr1->IsPoint()) || (!tr2->IsPoint()) )  return kFALSE; 
+    if(fabs(innertrack->GetPointY()-outertrack->GetPointY()) >fMaxY) return kFALSE;
+    if(fabs(innertrack->GetPointZ()-outertrack->GetPointZ()) >fMaxZ) return kFALSE;
+    if(fabs(innertrack->GetKappa()-outertrack->GetKappa())   >fMaxKappa) return kFALSE;
+    if(GetAngle(innertrack->GetPointPsi(),outertrack->GetPointPsi()) >fMaxPsi) return kFALSE;
+    if(fabs(innertrack->GetTgl()-outertrack->GetTgl()) >fMaxTgl) return kFALSE;
+  */
+  
+  return kTRUE;//Tracks could be merged
 }
 
-void AliL3HoughMerger::MergePatches()
+void AliL3HoughMerger::AddTrack(AliL3TrackArray *mergedtrack,AliL3Track *track)
 {
+  AliL3Track *t[1];
+  t[0] = track;
+  MultiMerge(mergedtrack,t,1);
+}
+
+AliL3Track *AliL3HoughMerger::MultiMerge(AliL3TrackArray *mergedtrack,AliL3Track **tracks, Int_t ntrack)
+{
+  //Called by the track merger
+
+  AliL3HoughTrack *newtrack = (AliL3HoughTrack*)mergedtrack->NextTrack();
+  AliL3HoughTrack **trs = (AliL3HoughTrack**)tracks;
+  Int_t weight=0;
+
+  //Sum up the total weight:
+  for(Int_t i=ntrack-1; i>=0; i--)
+    weight += trs[i]->GetWeight();
   
+  AliL3HoughTrack *tpt=trs[0];//This is the innermost track
+  AliL3HoughTrack *tpl=trs[ntrack-1];
+  newtrack->SetTrackParameters(tpt->GetKappa(),tpt->GetPhi0(),weight);
+  newtrack->SetEtaIndex(tpt->GetEtaIndex());
+  newtrack->SetEta(tpt->GetEta());
+  newtrack->SetPsi(tpt->GetPsi());
+  newtrack->SetCenterX(tpt->GetCenterX());
+  newtrack->SetCenterY(tpt->GetCenterY());
+  newtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
+  newtrack->SetLastPoint(tpl->GetLastPointX(),tpl->GetLastPointY(),tpl->GetLastPointZ());
+  newtrack->SetCharge(tpt->GetCharge());
+  newtrack->SetRowRange(tpt->GetFirstRow(),tpl->GetLastRow());
   
+  return (AliL3Track*)newtrack;
+
 }
 
-void AliL3HoughMerger::MergeEtaSlices(Int_t patch)
+void AliL3HoughMerger::MergePatches(Bool_t slow)
 {
-  AliL3TrackArray *tracks = fInTracks[patch];
-  Int_t ntracks = tracks->GetNTracks();
-  printf("Number of tracks to merging %d\n",ntracks);
-  AliL3HoughTrack *tr1,*tr2;
-  Double_t ptdiff=0.01;
-  Double_t phi0diff=0.01;
-  for(Int_t i=0; i<ntracks; i++)
+  //Merge tracks from across the patches.
+  
+  fSlow = slow;
+  AliL3TrackArray *tracks;
+  AliL3HoughTrack *track;
+  for(Int_t i=0; i<GetNIn(); i++)
     {
-      tr1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
-      if(!tr1) continue;
-      for(Int_t j=0; j<ntracks; j++)
+      tracks = GetInTracks(i);
+      for(Int_t j=0; j<tracks->GetNTracks(); j++)
        {
-         tr2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
-         if(!tr2) continue;
-         if(tr1==tr2) continue;
-         if(tr1->GetEtaIndex() == tr2->GetEtaIndex()) continue;
-         if(tr1->GetEtaIndex() == tr2->GetEtaIndex()-1 || tr1->GetEtaIndex() == tr2->GetEtaIndex()+1)
-           if(fabs(tr1->GetPt()-tr2->GetPt())<ptdiff && fabs(tr1->GetPhi0()-tr2->GetPhi0())<phi0diff) 
-             MergeTracks(tracks,i,j);
+         track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
+         if(!track) continue;
+         track->UpdateToFirstRow();
        }
     }
-  tracks->Compress();
+  Merge();
   
-  printf("Number of tracks that was not merged: %d\n",tracks->GetNTracks());
-  //Add rest of the tracks, which was not merged
-  for(Int_t i=0; i<tracks->GetNTracks(); i++)
-    {
-      tr1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
-      if(!tr1) continue;
-      tr2 = (AliL3HoughTrack*)fOutTracks->NextTrack();
-      tr2->Set(tr1);
+}
+
+void AliL3HoughMerger::Merge()
+{
+  Double_t edge0 = PI/18.;
+  Double_t edge1 = 2*PI - edge0;
+  AliL3TrackArray *ttt = GetOutTracks();
+  
+  Int_t subsec = GetNIn() - 2; 
+  for(Int_t i=subsec;i>=0;i--){
+    AliL3TrackArray *tout = GetOutTracks();
+    if(i==subsec) tout = GetInTracks(subsec+1);
+    AliL3TrackArray *tin = GetInTracks(i);
+    Double_t xval = fTransformer->Row2X(NRows[i][1]);
+    Double_t xmax = fTransformer->Row2X(NRows[i+1][1]);
+    Double_t ymax = xval*tan(edge0);
+    for(Int_t out=0;out<tout->GetNTracks();out++){
+      AliL3Track *outtrack=tout->GetCheckedTrack(out);
+      if(!outtrack) continue;
+      //outtrack->CalculateHelix();
+      outtrack->CalculatePoint(xval);
+      if(outtrack->IsPoint()&&fabs(outtrack->GetPointY())>ymax){
+       tout->Remove(out);
+      }
     }
-  printf("Total number of tracks after merging: %d\n",fOutTracks->GetNTracks());
+    //    tout->Compress();
+    for(Int_t in=0;in<tin->GetNTracks();in++){
+      AliL3Track *intrack=(AliL3Track*)tin->GetTrack(in);
+      //intrack->CalculateHelix();
+      intrack->CalculatePoint(xval);
+    }
+    tin->QSort();
+    tout->QSort();
+
+    if(fSlow) SlowMerge(ttt,tin,tout,xval);
+    else Merge(ttt,tin,tout);
+
+    /*
+    //Add the tracks that cross the sector boundary:
+    for(Int_t in=0;in<tin->GetNTracks();in++){
+      AliL3Track *intrack=(AliL3Track*)tin->GetCheckedTrack(in);
+      if(!intrack) continue;
+      if(intrack->CalculateEdgePoint(edge0)){
+        if(intrack->GetPointX()<xmax ){
+          AddTrack(ttt,intrack);
+          tin->Remove(in);
+        }
+      } 
+      else if(intrack->CalculateEdgePoint(edge1)){
+        if(intrack->GetPointX()<xmax ){
+          AddTrack(ttt,intrack);
+          tin->Remove(in);
+        }
+      }
+    }
+    */
+  } // end subsector loop
+  LOG(AliL3Log::kInformational,"AliL3HoughMerger::Merge","Result")
+    <<AliL3Log::kDec<<"Total Merged Tracks: "<<GetOutTracks()->GetNPresent()
+    <<ENDLOG;
 }
 
-void AliL3HoughMerger::MergeTracks(AliL3TrackArray *intracks,Int_t i,Int_t j)
+Int_t AliL3HoughMerger::Merge(AliL3TrackArray* mergedtrack,AliL3TrackArray *tracksin,AliL3TrackArray *tracksout)
 {
-  AliL3HoughTrack *newtrack = (AliL3HoughTrack*)fOutTracks->NextTrack();
-  AliL3HoughTrack *mergetrack = (AliL3HoughTrack*)intracks->GetCheckedTrack(i);
-  AliL3HoughTrack *mergetrack2 = (AliL3HoughTrack*)intracks->GetCheckedTrack(j);
-  if(!mergetrack || !mergetrack2)
-    {
-      printf("\nALiL3HoughMerger::MergeTracks : NO TRACK!!!\n");
-      return;
+  
+  AliL3Track *tracks[2];
+  const Int_t  kNOut=tracksout->GetNTracks();
+  const Int_t  kNIn =tracksin->GetNTracks();
+  const Int_t  kNMerged =mergedtrack->GetNTracks();
+
+  Bool_t *ismatchedin  = new Bool_t[kNIn];
+  for(Int_t in =0;in<kNIn;in++)
+    ismatchedin[in]=kFALSE;
+  Bool_t *ismatchedout = new Bool_t[kNOut];
+  for(Int_t out =0;out<kNOut;out++)
+    ismatchedout[out] = kFALSE;
+  for(Int_t out =0;out<kNOut;out++){
+    AliL3Track *outertrack=(AliL3Track*)tracksout->GetCheckedTrack(out);
+    if(!outertrack) continue;
+    for(Int_t in =0;in<kNIn;in++){
+      if(ismatchedin[in]) continue;
+      AliL3Track *innertrack=(AliL3Track*)tracksin->GetCheckedTrack(in);
+      if(!innertrack) continue;
+      if(outertrack==innertrack) continue;
+      
+      if(IsTrack(innertrack,outertrack)) //They can be merged
+       {
+         tracks[0]=innertrack; tracks[1]=outertrack; 
+         SortTracks(tracks,2); //Sort the tracks according to minimum x-point
+         //if(tracks[0]->GetLastPointX()<tracks[1]->GetFirstPointX()){
+         MultiMerge(mergedtrack,tracks,2);
+         tracksout->Remove(out);
+         tracksin->Remove(in);
+         ismatchedin[in]=kTRUE;
+         ismatchedout[out]=kTRUE;
+         break;
+         // }
+       }
     }
-  Int_t w = mergetrack->GetWeight() + mergetrack2->GetWeight();
-  newtrack->SetTrackParameters(mergetrack->GetKappa(),mergetrack->GetPhi0(),w);
-  newtrack->SetEtaIndex(mergetrack->GetEtaIndex());
-  newtrack->SetEta(mergetrack->GetEta());
+  }
   
-  intracks->Remove(i);
-  intracks->Remove(j);
+  Int_t nmerged = mergedtrack->GetNTracks()-kNMerged;
+  LOG(AliL3Log::kInformational,"AliL3HoughMerger::Merge","Result")
+    <<AliL3Log::kDec<<"Merged Tracks: "<<nmerged<<ENDLOG;
+  delete[] ismatchedin;
+  delete[] ismatchedout;
+  return nmerged;
+}
+
+void AliL3HoughMerger::SlowMerge(AliL3TrackArray *mergedtrack,AliL3TrackArray *tracksin,AliL3TrackArray *tracksout,Double_t xval)
+{
+  void *ntuple=GetNtuple();
+  const Int_t  kNOut=tracksout->GetNTracks();
+  const Int_t  kNIn =tracksin->GetNTracks();
+  const Int_t  kNMerged =mergedtrack->GetNTracks();
+  AliL3Track *tracks[2];
+  Bool_t merge = kTRUE;
+  while(merge){
+    Int_t inmin=-1,outmin=-1;
+    Double_t min=10;
+    for(Int_t out=0;out<kNOut;out++){
+    AliL3Track *outertrack=tracksout->GetCheckedTrack(out);
+    if(!outertrack) continue;
+      for(Int_t in=0;in<kNIn;in++){
+        AliL3Track *innertrack=tracksin->GetCheckedTrack(in);
+        if(!innertrack) continue;
+        Double_t diff = TrackDiff(innertrack,outertrack);
+        if(diff>=0&&diff<min){
+          min=diff;
+          inmin=in;
+          outmin=out; 
+        }
+      } 
+    }
+    if(inmin>=0&&outmin>=0){
+      AliL3Track *outertrack=tracksout->GetTrack(outmin);
+      AliL3Track *innertrack=tracksin->GetTrack(inmin);
+      tracks[0]=innertrack;
+      tracks[1]=outertrack;
+      SortTracks(tracks,2);
+      Print(tracks);
+      MultiMerge(mergedtrack,tracks,2);
+      outertrack->CalculatePoint(xval);
+      innertrack->CalculatePoint(xval);
+      FillNtuple(ntuple,innertrack,outertrack);
+      tracksout->Remove(outmin);
+      tracksin->Remove(inmin);
+      //      tracksout->Compress();
+      //      tracksin->Compress(); 
+    }
+    else merge = kFALSE;
+  }
+  LOG(AliL3Log::kInformational,"AliL3HoughMerger::SlowMerge","Result")
+    <<AliL3Log::kDec<<"Merged Tracks: "
+    <<mergedtrack->GetNTracks()-kNMerged<<ENDLOG;
+  char name[256] = "ntuple_t.root";
+  for(Int_t i=0;i<GetNIn();i++)
+    if(tracksin==GetInTracks(i))
+      sprintf(name,"ntuple_t_%d.root",i);
+  WriteNtuple(name,ntuple);
+}
+
+void AliL3HoughMerger::Print(AliL3Track **tracks)
+{
+  AliL3HoughTrack *tr1 = (AliL3HoughTrack*)tracks[0];
+  AliL3HoughTrack *tr2 = (AliL3HoughTrack*)tracks[1];
+  Double_t kappadiff = fabs(tr1->GetKappa()-tr2->GetKappa());
+  Double_t phi0diff = fabs(tr1->GetPhi0()-tr2->GetPhi0());
+  cout << "---------Difference in merged tracks---------"<<endl;
+  cout << "Kappa: "<<kappadiff<<" Phi0 : "<<phi0diff<<endl;
   
 }
index 8b0ba95..015a5da 100644 (file)
@@ -2,36 +2,43 @@
 #define ALIL3_HOUGH_Merger
 
 #include "AliL3RootTypes.h"
+#include "AliL3Merger.h"
 
 class AliL3TrackArray;
+class AliL3Track;
 
-class AliL3HoughMerger {
+class AliL3HoughMerger : public AliL3Merger {
   
  private:
+  Double_t fMaxY;
+  Double_t fMaxZ;
+  Double_t fMaxKappa;
+  Double_t fMaxPsi;
+  Double_t fMaxTgl;
+  Double_t fMaxPhi0;
+  Bool_t fSlow;
+
+  void Merge();
+  Int_t Merge(AliL3TrackArray* mergedtrack,AliL3TrackArray *tracksin,AliL3TrackArray *tracksout);
+  void SlowMerge(AliL3TrackArray *mergedtrack,AliL3TrackArray *tracksin,AliL3TrackArray *tracksout,Double_t xval);
 
-  Int_t fNIn;
-  AliL3TrackArray **fInTracks; //!
-  AliL3TrackArray *fOutTracks; //!
-
-  void SetArray();
-  void DeleteArray();
-  
  public:
   AliL3HoughMerger(); 
   AliL3HoughMerger(Int_t nsubsectors);
   virtual ~AliL3HoughMerger();
   
+  virtual Bool_t IsTrack(AliL3Track *innertrack,AliL3Track *outertrack);
+  virtual AliL3Track *MultiMerge(AliL3TrackArray *mergedtrack,AliL3Track **tracks, Int_t ntrack);
+  virtual void AddTrack(AliL3TrackArray *mergedtrack,AliL3Track *track);
   void FillTracks(AliL3TrackArray *tracks,Int_t patch);
-  void MergePatches();
+  
+  void MergePatches(Bool_t slow=kTRUE);
   void MergeEtaSlices(Int_t patch);
   void MergeTracks(AliL3TrackArray *intracks,Int_t i,Int_t j);
-
+  void Print(AliL3Track **tracks);
+  void SetParameters(Double_t maxkappa=0.001, Double_t maxpsi=0.05,Double_t maxphi0=0.1);
   
-  //Getters
-  AliL3TrackArray *GetOutTracks() {return fOutTracks;}
-  AliL3TrackArray *GetInTracks(Int_t i) {if(!fInTracks) return 0; if(!fInTracks[i]) return 0; return fInTracks[i];}
-
-  ClassDef(AliL3HoughMerger,1)
+  ClassDef(AliL3HoughMerger,1) //Patch merger for Hough tracklets
 
 };
 
index ff006c3..ba45579 100644 (file)
@@ -1,3 +1,7 @@
+//$Id$
+
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV 
 
 #include <math.h>
 
@@ -6,6 +10,11 @@
 #include "AliL3Defs.h"
 #include "AliL3HoughTrack.h"
 
+//_____________________________________________________________
+// AliL3HoughTrack
+//
+// Track class for Hough tracklets
+
 ClassImp(AliL3HoughTrack)
 
 
@@ -44,6 +53,7 @@ void AliL3HoughTrack::Set(AliL3Track *track)
   SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
   SetCharge(tpt->GetCharge());
   SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
+  SetSlice(tpt->GetSlice());
   return;
 
   fWeight = tpt->GetWeight();
@@ -66,7 +76,8 @@ Int_t AliL3HoughTrack::Compare(const AliL3Track *tpt) const
 
 void AliL3HoughTrack::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 = Pi/2 - theta;
@@ -97,17 +108,20 @@ void AliL3HoughTrack::UpdateToFirstRow()
   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()) / ( BFACT * BField )  ;
   Double_t tPhi0 = GetPsi() + GetCharge() * 0.5 * pi / fabs(GetCharge()) ;
-  Double_t x0    = GetR0() * cos(GetPhi0()) ;
-  Double_t y0    = GetR0() * sin(GetPhi0()) ;
-  Double_t rc    = fabs(GetPt()) / ( BFACT * bField )  ;
-  Double_t xc    = x0 - rc * cos(tPhi0) ;
-  Double_t yc    = y0 - rc * sin(tPhi0) ;
+  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(AliL3Log::kError,"AliL3HoughTrack::UpdateToFirstRow","Tracks")<<AliL3Log::kDec<<
       "Track does not intersect"<<ENDLOG;
@@ -161,7 +175,7 @@ void AliL3HoughTrack::SetTrackParameters(Double_t kappa,Double_t phi,Int_t weigh
   fMinDist = 100000;
   SetKappa(kappa);
   SetPhi0(phi);
-  Double_t pt = fabs(BFACT*bField/kappa);
+  Double_t pt = fabs(BFACT*BField/kappa);
   SetPt(pt);
   Double_t radius = 1/fabs(kappa);
   SetRadius(radius);
index 93bb0c2..9261fb2 100644 (file)
@@ -14,6 +14,7 @@ class AliL3HoughTrack : public AliL3Track {
   Int_t fWeight;
   Int_t fEtaIndex;
   Double_t fEta;
+  Int_t fSlice; //The slice where this track was found
 
   Double_t fDLine;
   Double_t fPsiLine;
@@ -31,14 +32,17 @@ class AliL3HoughTrack : public AliL3Track {
   void UpdateToFirstRow();
   void SetTrackParameters(Double_t kappa,Double_t phi,Int_t weight);  
   void SetLineParameters(Double_t psi,Double_t D,Int_t weight,Int_t *rowrange,Int_t ref_row);
+
   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,Double_t *xy);
-  
+
+  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;}
@@ -46,7 +50,7 @@ class AliL3HoughTrack : public AliL3Track {
   void SetDLine(Double_t f) {fDLine=f;}
   void SetPsiLine(Double_t f) {fPsiLine=f;}
 
-  ClassDef(AliL3HoughTrack,1)
+  ClassDef(AliL3HoughTrack,1) //Track class for Hough tracklets
 
 };
 
index f9a81f7..7723a69 100644 (file)
@@ -1,5 +1,5 @@
-//Author:        Anders Strand Vestbo
-//Last Modified: 14.09.01
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV 
 
 #include "AliL3MemHandler.h"
 #include "AliL3Logging.h"
@@ -9,6 +9,11 @@
 #include "AliL3DigitData.h"
 #include "AliL3Histogram.h"
 
+//_____________________________________________________________
+// AliL3HoughTransformer
+//
+// Hough transformation class
+
 ClassImp(AliL3HoughTransformer)
 
 AliL3HoughTransformer::AliL3HoughTransformer()
@@ -158,6 +163,95 @@ void AliL3HoughTransformer::TransformCircle()
     }
 }
 
+void AliL3HoughTransformer::TransformCircleC()
+{
+  //Circle transform, using combinations of every 2 points lying
+  //on different padrows and within the same etaslice.
+  
+
+  Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
+  AliL3DigitRowData **rowPt = new AliL3DigitRowData*[nrows];
+  
+  AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
+  if(!tempPt)
+    printf("\nAliL3HoughTransformer::TransformCircleC() : Zero data pointer\n");
+  
+  Int_t prow;
+  for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
+    {
+      prow = i - NRows[fPatch][0];
+      rowPt[prow] = tempPt;
+      AliL3MemHandler::UpdateRowPointer(tempPt);
+    }
+  Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
+
+  AliL3DigitData *digPt;
+  Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
+  UShort_t charge1,charge2,time;
+  UChar_t pad;
+  Float_t xyz[3];
+  Int_t sector,row,eta_index1,eta_index2,tot_charge;
+  for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
+    {
+      prow = i - NRows[fPatch][0];
+      digPt = rowPt[prow]->fDigitData;
+      for(UInt_t di=0; di<rowPt[prow]->fNDigit; di++)
+       {
+         charge1 = digPt[di].fCharge;
+         pad = digPt[di].fPad;
+         time = digPt[di].fTime;
+         if(charge1 <= fThreshold)
+           continue;
+         fTransform->Slice2Sector(fSlice,i,sector,row);
+         fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
+         eta = fTransform->GetEta(xyz);
+         eta_index1 = (Int_t)((eta-fEtaMin)/etaslice);
+         if(eta_index1 < 0 || eta_index1 >= fNEtaSegments)
+           continue;
+         r1 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
+         phi1 = atan2(xyz[1],xyz[0]);
+         
+         //Get the correct histogrampointer:
+         AliL3Histogram *hist = fParamSpace[eta_index1];
+         if(!hist)
+           {
+             printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",eta_index1);
+             continue;
+           }
+
+         for(Int_t j=i+1; j<NRows[fPatch][1]; j++)
+           {
+             prow = j - NRows[fPatch][0];
+             digPt = rowPt[prow]->fDigitData;
+             for(UInt_t ni=0; ni<rowPt[prow]->fNDigit; ni++)
+               {
+                 charge2 = digPt[ni].fCharge;
+                 pad = digPt[ni].fPad;
+                 time = digPt[ni].fTime;
+                 if(charge2 <= fThreshold)
+                   continue;
+                 fTransform->Slice2Sector(fSlice,j,sector,row);
+                 fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
+                 eta = fTransform->GetEta(xyz);
+                 eta_index2 = (Int_t)((eta-fEtaMin)/etaslice);
+                 if(eta_index2 != eta_index1)
+                   continue;
+                 r2 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
+                 phi2 = atan2(xyz[1],xyz[0]);
+                 
+                 phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
+                 kappa = 2*sin(phi2-phi_0)/r2;
+                 tot_charge = charge1+charge2;
+                 //printf("Filling kappa %f phi %f charge %d\n",kappa,phi_0,tot_charge);
+                 hist->Fill(kappa,phi_0,tot_charge);
+               }
+           }
+       }
+    }
+
+  delete [] rowPt;
+}
+
 void AliL3HoughTransformer::TransformLine()
 {
   //Do a line transform on the data.
index f418e12..588df41 100644 (file)
@@ -39,6 +39,7 @@ class AliL3HoughTransformer : public TObject {
                        Int_t nybin=64,Double_t ymin=-0.26,Double_t ymax=0.26);
   void Reset();
   void TransformCircle();
+  void TransformCircleC();
   void TransformLine();
 
   //Getters
@@ -48,19 +49,20 @@ class AliL3HoughTransformer : public TObject {
   Int_t GetThreshold() {return fThreshold;}
   Double_t GetEtaMin() {return fEtaMin;}
   Double_t GetEtaMax() {return fEtaMax;}
+  Double_t GetEtaSlice() {return (fEtaMax - fEtaMin)/fNEtaSegments;}
   void *GetDataPointer() {return (void*)fDigitRowData;}
   AliL3Histogram *GetHistogram(Int_t eta_index);
   
   //setters
   void SetThreshold(Int_t i) {fThreshold = i;}
 
-  ClassDef(AliL3HoughTransformer,1)
+  ClassDef(AliL3HoughTransformer,1) //Hough transformation class
 
 };
 
 inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t eta_index)
 {
-  if(!fParamSpace || eta_index >= fNEtaSegments)
+  if(!fParamSpace || eta_index >= fNEtaSegments || eta_index < 0)
     return 0;
   if(!fParamSpace[eta_index])
     return 0;
index fad654e..18119d4 100644 (file)
@@ -5,10 +5,10 @@
 #endif
 
 
-void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitfile)
+void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitfile,Int_t event)
 {
 
-  Int_t good_number = 70;
+  Int_t good_number = 10; //Minimum number of points on a good track
   
   struct GoodTrack goodtracks[15000];
   Int_t nt=0;
@@ -28,7 +28,7 @@ void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitf
       return;
     }
   gAlice = (AliRun*)evfile->Get("gAlice");
-  gAlice->GetEvent(0);
+  gAlice->GetEvent(event);
 
   AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
   AliTPCParam *param = (AliTPCParam*)evfile->Get("75x40_100x60");
@@ -52,7 +52,9 @@ void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitf
   AliL3Transform *transform = new AliL3Transform();
 
   TFile *digfile = TFile::Open(digitfile);
-  TTree *tree = (TTree*)digfile->Get("TreeD_75x40_100x60_0");
+  Char_t dname[100];
+  sprintf(dname,"TreeD_75x40_100x60_%d",event);
+  TTree *tree=(TTree*)digfile->Get(dname);
   AliSimDigits da, *digits=&da;
   tree->GetBranch("Segment")->SetAddress(&digits);
   
@@ -79,7 +81,7 @@ void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitf
        if (idx0>=0 && dig>=zero) count[idx0]+=1;
        if (idx1>=0 && dig>=zero) count[idx1]+=1;
        if (idx2>=0 && dig>=zero) count[idx2]+=1;
-         } while (digits->Next());
+      } while (digits->Next());
       
       for (Int_t j=0; j<np; j++) 
        {
@@ -97,7 +99,7 @@ void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitf
   
   TTree *TH=gAlice->TreeH();
   Int_t npart=(Int_t)TH->GetEntries();
-    
+  cout<<"Found "<<npart<<" particles in this event"<<endl;
 
 
   while (npart--) {
@@ -120,13 +122,12 @@ void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitf
     if (hit1->fQ != 0.) continue;
     Int_t i=hit->Track();
     TParticle *p = (TParticle*)gAlice->Particle(i);
-    
+
     if (p->GetFirstMother()>=0) continue;  //secondary particle
     if (good[i] < good_number) continue;
     if (p->Pt() < 0.1) continue;
     if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
-    printf("checking paricle %d with %d hits \n",i,good[i]);
-    
+    cout<<"Checking particle with "<<good[i]<<" hits, and pt: "<<p->Pt()<<endl;
     goodtracks[nt].eta = p->Eta();
     goodtracks[nt].label=i;
     goodtracks[nt].code=p->GetPdgCode();
@@ -137,7 +138,7 @@ void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitf
     
     if (hit0) delete hit0;
   }
-  
+  cout<<"Found "<<nt<<" good tracks in this event"<<endl;
   ofstream out(fname);
   if(out) 
     {
index 12e93e9..25310dc 100644 (file)
@@ -13,6 +13,6 @@ struct GoodTrack
 };
 
 
-void GetGoodParticles(Int_t,Int_t,Char_t *,Char_t *);
+void GetGoodParticles(Int_t minslice,Int_t maxslice,Char_t *eventfile,Char_t *digifile,Int_t event=0);
 
 #endif
index d86f566..b8be7e8 100644 (file)
@@ -1,7 +1,13 @@
-void test(char *file="/prog/alice/data/Rawdata/6_patch/hg_42105_s1-3/",bool bin=true)
+void test(char *file="/prog/alice/data/Rawdata/1_patch/pp/event_0/",bool bin=true)
 {
+  AliL3Logger l;
+//  l.UnSet(AliL3Logger::kDebug);
+//  l.UnSet(AliL3Logger::kAll);
+  l.Set(AliL3Logger::kAll);
+  //l.UseStdout();
+  l.UseStream();
   
-  Int_t fNEtaSegments = 1 ;
+  Int_t fNEtaSegments = 1;
   
   Char_t histname[50];
   Int_t i;
@@ -14,16 +20,15 @@ void test(char *file="/prog/alice/data/Rawdata/6_patch/hg_42105_s1-3/",bool bin=
   
   fTransform = new AliL3Transform();
   
-  fHoughTransformer->CreateHistograms();//64,0,3.1415,128,-120,120);
+  fHoughTransformer->CreateHistograms();//64,-0.012,0.012,64,-0.26,0.26);
   fHoughTransformer->SetThreshold(10);
   fTracks = new AliL3TrackArray("AliL3HoughTrack");
-  
-  
+    
   TH2F *road = new TH2F("road","",250,0,250,250,-125,125);
   TH2F *peaks = new TH2F("peaks","",64,-0.006,0.006,64,-0.26,0.26);
   peaks->SetMarkerStyle(3);
   peaks->SetMarkerColor(2);
-  road->SetMarkerStyle(6);
+  road->SetMarkerStyle(5);
   road->SetMarkerColor(2);
   
   Char_t name[256];
@@ -44,19 +49,21 @@ void test(char *file="/prog/alice/data/Rawdata/6_patch/hg_42105_s1-3/",bool bin=
   else //read data from root file
     {
       const Int_t NRows[6][2] = {{0,31},{32,63},{64,91},{92,119},{120,143},{144,175}};
+      const int nr[2] = {0,175};
       fMemHandler->Free();
-      TFile *rootfile = new TFile(file);
-      fMemHandler->SetAliInput(rootfile);
+      fMemHandler->SetAliInput(file);
       fMemHandler->Init(slice,patch,NRows[patch]);
       fMemHandler->Init(fTransform);
       digits=(AliL3DigitRowData *)fMemHandler->AliDigits2Memory(ndigits); 
-      rootfile->Close();
-      delete rootfile;
+      fMemHandler->CloseAliInput();
     }
+    
   fHoughTransformer->SetInputData(ndigits,digits);
-  fEval = new AliL3HoughEval(fHoughTransformer);
-  fEval->SetNumOfRowsToMiss(2);
+  fEval = new AliL3HoughEval();
+  fEval->InitTransformer(fHoughTransformer);
+  //fEval->SetNumOfRowsToMiss(0);
   //fEval->SetNumOfPadsToLook(2);
+  //fEval->RemoveFoundTracks();
   AliL3HoughTrack *track;
   AliL3Histogram *histo;
   Int_t good_count;
@@ -68,9 +75,7 @@ void test(char *file="/prog/alice/data/Rawdata/6_patch/hg_42105_s1-3/",bool bin=
       fHoughTransformer->TransformCircle();
       double final = AliL3Benchmark::GetCpuTime();
       printf("done in %f ms\n",(final-init)*1000);
-      
       good_count=0;
-      /*
       for(Int_t e=0; e<fNEtaSegments; e++)
        {
          histo = (AliL3Histogram*)fHoughTransformer->GetHistogram(e);
@@ -78,28 +83,37 @@ void test(char *file="/prog/alice/data/Rawdata/6_patch/hg_42105_s1-3/",bool bin=
          
          fMaxFinder->SetHistogram(histo);
          Int_t n=10;
+         Int_t weight[10];
          Float_t x[10];
          Float_t y[10];
-         fMaxFinder->FindPeak1(x,y,n);
+         fMaxFinder->FindPeak1(x,y,weight,n);
          track = new AliL3HoughTrack();
          track->SetTrackParameters(x[0],y[0],1);
          
-         //if(!fEval->LookInsideRoad(track,e))
-         //continue;
-         for(int t=0; t<176; t++)
+         //      if(!fEval->LookInsideRoad(track,e))
+         // continue;
+         
+         
+         if(e==eind)
            {
-             float xyz_tr[3];
-             track->GetCrossingPoint(t,xyz_tr);
-             road->Fill(xyz_tr[0],xyz_tr[1],1);
+             for(int t=0; t<176; t++)
+               {
+                 if(t%10) continue;
+                 float xyz_tr[3];
+                 track->GetCrossingPoint(t,xyz_tr);
+                 road->Fill(xyz_tr[0],xyz_tr[1],1);
+               }
            }
+         delete track;
          
        }
-      */
+         
       break;
       if(good_count==0)
        break;
     }
   
+
   image = new AliL3Histogram("image","",250,0,250,250,-125,125);
   fEval->DisplayEtaSlice(eind,image);
   
@@ -109,55 +123,113 @@ void test(char *file="/prog/alice/data/Rawdata/6_patch/hg_42105_s1-3/",bool bin=
   histo = (AliL3Histogram*)fHoughTransformer->GetHistogram(eind);
   if(!histo)
     {printf("No histogram\n"); return;}
-  fHoughTransformer->GetHistogram(eind)->Draw("box");
-  //  peaks->Draw("same");
+  histo->Draw("box");
+
+  peaks->Draw("same");
   c1->cd(2);
   image->Draw("");
-  //road->Draw("same");
+  road->Draw("same");
 } 
 
-void process(char *path="/prog/alice/data/Rawdata/6_patch/hg_42105_s1-3/",bool bin=true)
+void process(char *digitfile="/prog/alice/data/Rawdata/6_patch/hg_42105_s1-3/",bool bin=true,char *trackfile="good_tracks")
 {
-  double torad = 3.1415/180;
-  a = new AliL3Hough(path,bin,1);
-  a->TransformSlice(1);
+  AliL3Logger l;
+  //  l.UnSet(AliL3Logger::kDebug);
+  //  l.UnSet(AliL3Logger::kAll);
+  //l.Set(AliL3Logger::kError);
+  //l.UseStdout();
+  l.UseStream();
   
-  hist = (AliL3Histogram*)a->AddHistograms();
-  
-  //hist->SetThreshold(10000);
+  double torad = 3.1415/180;
+  int slice=1;
+  a = new AliL3Hough(digitfile,bin,100);
+  a->ReadData(slice);
+  a->Transform();
+  a->AddAllHistograms();
+  a->FindTrackCandidates();
+  a->Evaluate(3);
+  a->MergePatches();
+  
+  //a->MergeInternally();
+  
+  //tracks = (AliL3TrackArray*)(a->GetInterMerger())->GetOutTracks();
+  //tracks = (AliL3TrackArray*)(a->GetMerger())->GetOutTracks();
+  tracks = (AliL3TrackArray*)a->GetTracks(0);
+  //a->GetEval(0)->CompareMC(tracks,"good_tracks_hg4000_s1");
+  
+  
+  d = new AliL3HoughDisplay();
+  d->SetTracks(tracks);
+  UInt_t size;
+  int patch = 0;
+  data = (AliL3DigitRowData*)a->GetMemHandler(patch)->GetDataPointer(size);
+  //d->ShowData(data,size,slice,patch);
+  //return;
+  d->DisplayEvent();
+  //return;
+  
+  for(int i=0; i<tracks->GetNTracks(); i++)
+    {
+      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+      if(!track) continue;
+      printf("found pt %f phi0 %f eta %f weight %d rowange %d %d\n",track->GetPt(),track->GetPhi0(),track->GetEta(),track->GetWeight(),track->GetFirstRow(),track->GetLastRow());
+    }
+  return;
+
+}
+
+void HistogramParticles(char *trackfile)
+{
+  struct GoodTrack 
+  {
+    
+    Int_t label;
+    Double_t eta;
+    Int_t code;
+    Double_t px,py,pz;
+    Double_t pt;
+    Int_t nhits;
+  };
+  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;
+           }
+       }
+    }
   
-  b = new AliL3HoughMaxFinder("KappaPhi");
-  b->SetHistogram(hist);
+  int fNEtaSegments = 150;
+  Double_t etaslice = (0.9 - 0)/fNEtaSegments;
+  Int_t *particles = new Int_t[fNEtaSegments];
+  for(Int_t i=0; i<fNEtaSegments; i++)
+    particles[i]=0;
   
-  Int_t xbin,ybin;
+  for(Int_t i=0; i<nt; i++)
+    {
+      //if(goodtracks[i].nhits < 150) continue;
+      if(goodtracks[i].pt < 0.5) continue;
+      Int_t particleindex = (Int_t)(goodtracks[i].eta/etaslice);
+      if(particleindex < 0 || particleindex >= fNEtaSegments) continue;
+      particles[particleindex]++;
+    }
   
-  Int_t n=10;
-  Float_t x[10];
-  Float_t y[10];
-  b->FindPeak1(x,y,n);
-  printf("peak at pt %f phi0 %f\n",0.2*0.003/x[0],y[0]/torad);
+  hist = new TH1F("hist","",21,0,20);
+  for(int i=0; i<fNEtaSegments; i++)
+    hist->Fill(particles[i]);
   
-  track = new AliL3HoughTrack();
-  track->SetTrackParameters(x[0],y[0],1);
   
-  image = new AliL3Histogram("image","",250,0,250,250,-125,125);
-  a->Evaluate(image);
-  TH2F *road = new TH2F("road","",250,0,250,250,-125,125);
-  road->SetMarkerStyle(5);
-  road->SetMarkerColor(2);
+  c1 = new TCanvas("c1","",2);
+  hist->Draw();
   
-  float xyz[3];
-  for(int i=0; i<176; i++)
-    {
-      if(i%10) continue;
-      track->GetCrossingPoint(i,xyz);
-      road->Fill(xyz[0],xyz[1],1);
-    }
-  c1 = new TCanvas("c1","",1000,500);
-  c1->Divide(2);
-  c1->cd(1);
-  hist->Draw("box");
-  c1->cd(2);
-  image->Draw();
-  road->Draw("same");
 }