Full Hough->deconvolution->fit chain implemented (C.Cheshkov)

diff --git a/MONITOR/AliMonitorHLTHough.cxx b/MONITOR/AliMonitorHLTHough.cxx
index 540f906..da310f5 100644 (file)
--- a/MONITOR/AliMonitorHLTHough.cxx
+++ b/MONITOR/AliMonitorHLTHough.cxx
@@ -30,7 +30,9 @@
 #include <stdlib.h>
 #include <AliL3MemHandler.h>
 #include <AliL3TrackArray.h>
+#include <AliL3SpacePointData.h>
 #include <AliL3HoughTrack.h>
+#include <AliL3Transform.h>
 #endif
 
 //_____________________________________________________________________________
@@ -64,6 +66,25 @@ void AliMonitorHLTHough::CreateHistos(TFolder* folder)
 
   fFolder = folder->AddFolder("HLTHOUGH", "HLTHOUGH");
 
+  fClustersCharge = CreateHisto1("ClustersCharge", 
+                                "charge distribution of clusters", 
+                                100, 0, 1000, "charge", "#Delta N/N",
+                                AliMonitorHisto::kNormEvents);
+
+  Int_t nRows = fParam->GetNRowLow() + fParam->GetNRowUp();
+  fNClustersVsRow = CreateHisto1("NClustersVsRow", 
+                                "mean number of clusters per pad row", 
+                                nRows, -0.5, nRows-0.5,
+                                "pad row", "<N_{clusters}>",
+                                AliMonitorHisto::kNormEvents);
+
+  Int_t nSector = fParam->GetNInnerSector();
+  fNClustersVsSector = CreateHisto1("NClustersVsSector", 
+                                   "mean number of clusters per sector", 
+                                   nSector, -0.5, nSector-0.5, 
+                                   "sector", "<N_{clusters}>",
+                                   AliMonitorHisto::kNormEvents);
+
   fNTracks = CreateTrend("NTracks", "number of tracks per event", 
                         "N_{tracks}");
 
@@ -79,6 +100,20 @@ void AliMonitorHLTHough::CreateHistos(TFolder* folder)
                           120, 0, 360, "#phi [#circ]", "#Delta N/N",
                           AliMonitorHisto::kNormEntries);
 
+  fTrackNHits = CreateHisto1("TrackNHits", "Number of hits per track", 
+                          200, 0, 200, "N_{hits}", "#Delta N/N",
+                          AliMonitorHisto::kNormNone);
+
+  fTrackDEdxVsP = CreateHisto2("TrackDEdxVsP", "dE/dx of tracks", 
+                              100, 0, 3, 100, 0, 1000, 
+                              "p [GeV/c]", "dE/dx", "#Delta N/N",
+                              AliMonitorHisto::kNormEntries);
+
+  fTrackDEdx = CreateHisto1("TrackDEdx", "dE/dx for tracks with 0.4<p<1.0 GeV/c", 
+                              50, 0, 300, 
+                              "dE/dx", "#Delta N/N",
+                              AliMonitorHisto::kNormEntries);
+
   fTrackEtaVsPhi = CreateHisto2("TrackEtaVsPhi", "#phi vs #eta", 
                               20, -1, 1, 25, 0, 360, 
                               "#eta", "#phi", "#Delta N/N",
@@ -102,38 +137,111 @@ void AliMonitorHLTHough::FillHistos(AliRunLoader* /*runLoader*/,
   Warning("FillHistos", "the code was compiled without HLT support");
 
 #else
-  AliL3MemHandler memHandler[36];
-  Int_t nHoughTracks = 0;
+  AliL3MemHandler clusterHandler[36][6];
+  AliL3SpacePointData *clusters[36][6];
   for (Int_t iSector = 0; iSector < fParam->GetNInnerSector(); iSector++) {
-    char fileName[256];
-    sprintf(fileName, "hlt/tracks_ho_%d.raw", iSector);
-    if (!memHandler[iSector].SetBinaryInput(fileName)) {
-      Warning("FillHistos", "could not open file hlt/tracks.raw");
-      return;
-    }
-    AliL3TrackArray* tracks = new AliL3TrackArray;
-    memHandler[iSector].Binary2TrackArray(tracks);
-
-    nHoughTracks += tracks->GetNTracks();
-    for (Int_t iTrack = 0; iTrack < tracks->GetNTracks(); iTrack++) {
-      AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(iTrack);
-      if(!track) continue;
-
-      track->CalculateHelix();
-      track->Rotate(iSector);
-
-      fTrackPt->Fill(track->GetPt());
-      fTrackEta->Fill(track->GetPseudoRapidity());
-      fTrackPhi->Fill(track->GetPsi() * TMath::RadToDeg());
-      if(track->GetPt()>3.) {
-       fTrackEtaVsPhi->Fill(track->GetPseudoRapidity(),track->GetPsi() * TMath::RadToDeg());
-       fPtEtaVsPhi->Fill(track->GetPseudoRapidity(),track->GetPsi() * TMath::RadToDeg(),track->GetPt());
+    for (Int_t iPatch = 0; iPatch < 6; iPatch++) {
+      char fileName[256];
+      sprintf(fileName, "hlt/fitter/points_%d_%d.raw", iSector,iPatch);
+      if (!clusterHandler[iSector][iPatch].SetBinaryInput(fileName)) {
+       Warning("FillHistos", "could not open file %s", fileName);
+       continue;
       }
+      clusters[iSector][iPatch] = (AliL3SpacePointData*) clusterHandler[iSector][iPatch].Allocate();
+      UInt_t nClusters = 0;
+      clusterHandler[iSector][iPatch].Binary2Memory(nClusters, clusters[iSector][iPatch]);
+
+      for (UInt_t iCluster = 0; iCluster < nClusters; iCluster++) {
+       AliL3SpacePointData& cluster = clusters[iSector][iPatch][iCluster];
+       fClustersCharge->Fill(cluster.fCharge);
+       fNClustersVsRow->Fill(cluster.fPadRow);
+       fNClustersVsSector->Fill(iSector);
+      }
+
+      clusterHandler[iSector][iPatch].CloseBinaryInput();
+    }
+  }
+
+  fNClustersVsSector->ScaleErrorBy(10.);
+
+
+  AliL3MemHandler memHandler;
+  Int_t nHoughTracks = 0;
+
+  char fileName[256];
+  sprintf(fileName, "hlt/fitter/tracks.raw");
+  if (!memHandler.SetBinaryInput(fileName)) {
+    Warning("FillHistos", "could not open file hlt/fitter/tracks.raw");
+    return;
+  }
+  AliL3TrackArray* tracks = new AliL3TrackArray;
+  memHandler.Binary2TrackArray(tracks);
+
+  nHoughTracks += tracks->GetNTracks();
+  for (Int_t iTrack = 0; iTrack < tracks->GetNTracks(); iTrack++) {
+    AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(iTrack);
+    if(!track) continue;
+
+    track->CalculateHelix();
+    //  track->Rotate(iSector);
+
+    fTrackPt->Fill(track->GetPt());
+    fTrackEta->Fill(track->GetPseudoRapidity());
+    fTrackPhi->Fill(track->GetPsi() * TMath::RadToDeg());
+    if(track->GetPt()>3.) {
+      fTrackEtaVsPhi->Fill(track->GetPseudoRapidity(),track->GetPsi() * TMath::RadToDeg());
+      fPtEtaVsPhi->Fill(track->GetPseudoRapidity(),track->GetPsi() * TMath::RadToDeg(),track->GetPt());
+    }
+    fTrackNHits->Fill(track->GetNHits());
+
+    // Track dEdx
+    Int_t nc=track->GetNHits();
+    UInt_t *hits = track->GetHitNumbers();
+    Float_t sampleDEdx[159];
+    for (Int_t iHit = 0; iHit < nc; iHit++) {
+      UInt_t hitID = hits[iHit];
+      Int_t iSector = (hitID>>25) & 0x7f;
+      Int_t iPatch = (hitID>>22) & 0x7;
+      UInt_t position = hitID&0x3fffff;
+      UChar_t padrow = clusters[iSector][iPatch][position].fPadRow;
+      Float_t pWidth = AliL3Transform::GetPadPitchWidthLow();
+      if (padrow>63)
+       pWidth = AliL3Transform::GetPadPitchWidthUp(); 
+      Float_t corr=1.; if (padrow>63) corr=0.67;
+      sampleDEdx[iHit] = clusters[iSector][iPatch][position].fCharge/pWidth*corr;
+      Double_t crossingangle = track->GetCrossingAngle(padrow,iSector);
+      Double_t s = sin(crossingangle);
+      Double_t t = track->GetTgl();
+      sampleDEdx[iHit] *= TMath::Sqrt((1-s*s)/(1+t*t));
     }
-    fNTracks->Fill(nHoughTracks);
 
-    delete tracks;
-    memHandler[iSector].CloseBinaryInput();
+    /* Cook dEdx */
+    Int_t i;
+    Int_t swap;//stupid sorting
+    do {
+      swap=0;
+      for (i=0; i<nc-1; i++) {
+       if (sampleDEdx[i]<=sampleDEdx[i+1]) continue;
+       Float_t tmp=sampleDEdx[i];
+       sampleDEdx[i]=sampleDEdx[i+1]; sampleDEdx[i+1]=tmp;
+       swap++;
+      }
+    } while (swap);
+
+    Double_t low=0.05; Double_t up=0.7;
+    Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
+    Float_t trackDEdx=0;
+    for (i=nl; i<=nu; i++) trackDEdx += sampleDEdx[i];
+    trackDEdx /= (nu-nl+1);
+
+    fTrackDEdxVsP->Fill(track->GetP(),trackDEdx);
+    if(track->GetP()>0.4 && track->GetP()<1.0)
+      fTrackDEdx->Fill(trackDEdx);
   }
+  fNTracks->Fill(nHoughTracks);
+
+  delete tracks;
+  memHandler.CloseBinaryInput();
+
 #endif
 }

diff --git a/MONITOR/AliMonitorHLTHough.h b/MONITOR/AliMonitorHLTHough.h
index 9d4df4e..2ad266b 100644 (file)
--- a/MONITOR/AliMonitorHLTHough.h
+++ b/MONITOR/AliMonitorHLTHough.h
@@ -22,10 +22,16 @@ public:
 private:
   AliTPCParam*     fParam;              // TPC parameters
 
+  AliMonitorHisto* fClustersCharge;     // charge distribution of HLT clusters
+  AliMonitorHisto* fNClustersVsRow;     // mean number of HLT clusters per pad row
+  AliMonitorHisto* fNClustersVsSector;  // mean number of HLT clusters per sector
   AliMonitorTrend* fNTracks;            // number of HLT tracks per event
   AliMonitorHisto* fTrackPt;            // pt distribution of HLT tracks
   AliMonitorHisto* fTrackEta;           // eta distribution of HLT tracks
   AliMonitorHisto* fTrackPhi;           // phi distribution of HLT tracks
+  AliMonitorHisto* fTrackNHits;         // number of hits per HLT track
+  AliMonitorHisto* fTrackDEdxVsP;       // dedx distribution of HLT tracks
+  AliMonitorHisto* fTrackDEdx;          // dedx distribution of HLT tracks for a given momentum region
   AliMonitorHisto* fTrackEtaVsPhi;      // phi vs eta for HLT tracks
   AliMonitorHisto* fPtEtaVsPhi;         // phi vs eta for HLT tracks
 

diff --git a/MONITOR/AliMonitorProcess.cxx b/MONITOR/AliMonitorProcess.cxx
index 7144ddb..b663163 100644 (file)
--- a/MONITOR/AliMonitorProcess.cxx
+++ b/MONITOR/AliMonitorProcess.cxx
@@ -52,9 +52,15 @@
 #ifdef ALI_HLT
 #include <AliLevel3.h>
 #include <AliL3Transform.h>
+#include <AliL3Track.h>
+#include <AliL3TrackArray.h>
 #include <AliL3StandardIncludes.h>
 #include <AliL3HoughMaxFinder.h>
+#include <AliL3HoughBaseTransformer.h>
 #include <AliL3Hough.h>
+#include <AliL3ClusterFitter.h>
+#include <AliL3Vertex.h>
+#include <AliL3Fitter.h>
 #endif
 
 ClassImp(AliMonitorProcess) 
@@ -188,6 +194,7 @@ AliMonitorProcess::~AliMonitorProcess()
 
 #ifdef ALI_HLT
   delete fHLT;
+  delete fHLTHough;
 #endif
 }
 
@@ -619,10 +626,12 @@ void AliMonitorProcess::CreateHLTHough(const char* fileName)
   strcpy(name, fileName);
 
   fHLTHough = new AliL3Hough();
-  fHLTHough->SetTransformerParams(64,64,0.1,30);
-  fHLTHough->Init("./", kFALSE, 20, kFALSE,0,name);
+  fHLTHough->SetThreshold(4);
+  fHLTHough->SetTransformerParams(140,150,0.2,-1);
+  fHLTHough->SetPeakThreshold(6000,-1);
+  fHLTHough->Init("./", kFALSE, 50, kFALSE,0,name);
   fHLTHough->SetAddHistograms();
-  fHLTHough->GetMaxFinder()->SetThreshold(55000); // or 14000 ?
+  //  fHLTHough->GetMaxFinder()->SetThreshold(14000);
 
 }
 #endif
@@ -683,6 +692,10 @@ Bool_t AliMonitorProcess::ReconstructHLTHough(
   return kTRUE;
 
 #else
+  gSystem->Exec("rm -rf hlt/hough");
+  gSystem->MakeDirectory("hlt/hough");
+  gSystem->Exec("rm -rf hlt/fitter");
+  gSystem->MakeDirectory("hlt/fitter");
   if (!fHLTHough) return kFALSE;
 
   //  fHLTHough->Process(0, 35);
@@ -694,12 +707,77 @@ Bool_t AliMonitorProcess::ReconstructHLTHough(
       //      if(fHLTHough->fAddHistograms)
       fHLTHough->AddAllHistograms();
       fHLTHough->FindTrackCandidates();
-      fHLTHough->WriteTracks(i,"./hlt");
+      fHLTHough->AddTracks();
     }
+  fHLTHough->WriteTracks("./hlt/hough");
+
+  // Run cluster fitter
+  AliL3ClusterFitter *fitter = new AliL3ClusterFitter("./hlt");
+
+  // Set debug flag for the cluster fitter
+  //  fitter->Debug();
+
+  // Setting fitter parameters
+  fitter->SetInnerWidthFactor(1,1.5);
+  fitter->SetOuterWidthFactor(1,1.5);
+  fitter->SetNmaxOverlaps(5);
+  
+  //fitter->SetChiSqMax(5,kFALSE); //isolated clusters
+  fitter->SetChiSqMax(5,kTRUE);  //overlapping clusters
+
+  Int_t rowrange[2] = {0,AliL3Transform::GetNRows()-1};
+
+  // Takes input from global hough tracks produced by HT
+  fitter->LoadSeeds(rowrange,kFALSE,iEvent);
+
+  for(Int_t islice = 0; islice <= 35; islice++)
+    {
+      for(Int_t ipatch = 0; ipatch < AliL3Transform::GetNPatches(); ipatch++)
+       {
+         fitter->Init(islice,ipatch);
+
+         // Read digits
+         AliL3DigitRowData *digits = fHLTHough->GetTransformer(ipatch)->GetDataPointer();
+         fitter->SetInputData(digits);
+         fitter->FindClusters();
+         fitter->WriteClusters();
+       }
+    }
+
+  // Refit of the clusters
+  AliL3Vertex vertex;
+  //The seeds are the input tracks from circle HT
+  AliL3TrackArray *tracks = fitter->GetSeeds();
+  AliL3Fitter *ft = new AliL3Fitter(&vertex,1);
+
+  // Temporary solution !
+  //  sprintf(command, "rename _5.raw _-1.raw hlt/fitter/points*.raw", iEvent);
+  //  gSystem->Exec(command);
+
+  ft->LoadClusters("./hlt/fitter/",iEvent,kFALSE);
+  for(Int_t i=0; i<tracks->GetNTracks(); i++)
+    {
+      AliL3Track *track = tracks->GetCheckedTrack(i);
+      if(!track) continue;
+      if(track->GetNHits() < 40) continue;
+      ft->SortTrackClusters(track);
+      ft->FitHelix(track);
+      ft->UpdateTrack(track);
+    }
+  delete ft;
+        
+  //Write the final tracks
+  fitter->WriteTracks(5);
+
+  delete fitter;
 
   // remove the event number from the file names
   char command[256];
-  sprintf(command, "rename tracks_ho_%d tracks_ho hlt/*.raw", iEvent);
+  sprintf(command, "rename tracks_%d tracks hlt/hough/*.raw", iEvent);
+  gSystem->Exec(command);
+  sprintf(command, "rename tracks_%d tracks hlt/fitter/*.raw", iEvent);
+  gSystem->Exec(command);
+  sprintf(command, "rename points_%d points hlt/fitter/*.raw", iEvent);
   gSystem->Exec(command);
   return kTRUE;
 #endif