]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HLT/kalman/AliL3Kalman.cxx
Merged HLT tag v1-2 with ALIROOT tag v3-09-Release.
[u/mrichter/AliRoot.git] / HLT / kalman / AliL3Kalman.cxx
diff --git a/HLT/kalman/AliL3Kalman.cxx b/HLT/kalman/AliL3Kalman.cxx
new file mode 100644 (file)
index 0000000..0e1e146
--- /dev/null
@@ -0,0 +1,363 @@
+#include "AliL3StandardIncludes.h"
+#include <sys/time.h>
+
+#include "AliL3Benchmark.h"
+#ifdef use_aliroot
+#include "AliL3FileHandler.h"
+#else
+#include "AliL3MemHandler.h"
+#endif
+#include "AliL3DataHandler.h"
+#include "AliL3Transform.h"
+#include "AliL3SpacePointData.h"
+#include "AliL3DigitData.h"
+#include "AliL3Kalman.h"
+#include "AliL3Logging.h"
+#include "AliL3TrackArray.h"
+#include "AliL3Track.h"
+#include "AliL3KalmanTrack.h"
+
+#include "TNtuple.h"
+
+/*
+  AliL3Kalman
+*/
+
+ClassImp(AliL3Kalman)
+
+AliL3Kalman::AliL3Kalman(Char_t *datapath, Int_t *slice, Int_t min_clusters = 0)
+{
+  // Constructor
+
+  if (slice)
+    {
+      fMinSlice = slice[0];
+      fMaxSlice = slice[1];
+    }
+  else
+    {
+      fMinSlice = 0;
+      fMaxSlice = 35;
+    }
+  
+  sprintf(fPath,"%s",datapath);
+  fMinPointsOnTrack = 0; //min_clusters;
+  fTracks = 0;
+  fKalmanTracks = 0;
+
+  // NB! fNrow under only for single-patch, must include other possibilities 
+  // later on. ?? Maybe better also to put it in an Init-function
+  fRow[0][0] = 0;
+  fRow[0][1] = AliL3Transform::GetLastRow(-1);
+
+  fBenchmark = 0;
+}  
+
+AliL3Kalman::~AliL3Kalman()
+{
+  // Destructor
+  if (fBenchmark) delete fBenchmark;
+}
+
+void AliL3Kalman::Init()
+{
+  fBenchmark = new AliL3Benchmark();
+}
+
+void AliL3Kalman::LoadTracks(Int_t event, Bool_t sp)
+{
+  // Load tracks from conformal tracker
+  // Must also be possible to take seeds (and clusters) from Hough-transform??
+
+  // Load spacepoints into clusterfile
+
+  Double_t initTime,cpuTime;
+  initTime = GetCpuTime();
+  fBenchmark->Start("Load tracks");
+
+  Char_t fname[1024];
+  AliL3FileHandler *clusterfile[36][6];
+  for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
+    {
+      for(Int_t p=0; p<AliL3Transform::GetNPatches(); p++)
+        {
+          Int_t patch;
+          if(sp==kTRUE)
+            patch=-1;
+          else
+            patch=p;
+         
+          fClusters[s][p] = 0;
+          clusterfile[s][p] = new AliL3FileHandler();
+          if(event<0)
+            sprintf(fname,"%s/points_%d_%d.raw",fPath,s,patch);
+          else
+            sprintf(fname,"%s/points_%d_%d_%d.raw",fPath,event,s,patch);
+          if(!clusterfile[s][p]->SetBinaryInput(fname))
+            {
+              LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
+                <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
+              delete clusterfile[s][p];
+              clusterfile[s][p] = 0;
+              continue;
+            }
+          fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
+          clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
+          clusterfile[s][p]->CloseBinaryInput();
+         if(sp==kTRUE)
+            break;
+         delete clusterfile[s][p];
+        }
+    }
+
+  // Load tracks
+  sprintf(fname,"%s/tracks_%d.raw",fPath,event);
+  AliL3FileHandler *tfile = new AliL3FileHandler();
+  if(!tfile->SetBinaryInput(fname)){
+    LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
+      <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
+    return;
+  }
+  if(fTracks)
+    delete fTracks;
+  fTracks = new AliL3TrackArray();
+  tfile->Binary2TrackArray(fTracks);
+  tfile->CloseBinaryInput();
+  delete tfile;
+  cout << "Number of loaded tracks " << fTracks->GetNTracks() << endl;
+  fBenchmark->Stop("Load tracks");
+  cpuTime = GetCpuTime() - initTime;
+  LOG(AliL3Log::kInformational,"AliL3Kalman::LoadTracks()","Timing")
+    <<"Loading tracks in "<<cpuTime*1000<<" ms"<<ENDLOG;
+
+}
+
+void AliL3Kalman::ProcessTracks()
+{
+
+  Double_t initTime,cpuTime;
+  initTime = GetCpuTime();
+
+  fBenchmark->Start("Process tracks");
+
+  Int_t fEvent = 0;
+  fTracks->QSort();
+
+  fKalmanTracks = new AliL3TrackArray();
+
+  // Make a ntuple to store state vector, covariance matrix and chisquare
+  TNtuple *kalmanTree = new TNtuple("kalmanTree","kalmantracks","x0:x1:x2:x3:x4:c0:c1:c2:c3:c4:c5:c6:c7:c8:c9:c10:c11:c12:c13:c14:chisq");
+  Float_t meas[21];
+
+  // Go through the tracks from conformal or hough tracker
+  for (Int_t iTrack = 0; iTrack < fTracks->GetNTracks(); iTrack++)
+    {
+      AliL3KalmanTrack *kalmantrack = new AliL3KalmanTrack();
+      kalmantrack->Init();
+      AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(iTrack);
+      if (!track) continue;
+      if (track->GetNumberOfPoints() < fMinPointsOnTrack) continue;    
+
+      Bool_t save = kTRUE;
+
+      if (MakeSeed(kalmantrack, track) == 0)
+       {
+         save = kFALSE;
+         continue;
+       }
+
+      if (Propagate(kalmantrack, track) == 0) 
+       {
+         save = kFALSE;
+       }
+
+      if (save) {
+       Float_t x[5]; 
+       kalmantrack->GetStateVector(x);
+       Float_t c[15]; 
+       kalmantrack->GetCovariance(c);
+       Float_t chisq = kalmantrack->GetChisq();
+       meas[0]  = x[0];
+       meas[1]  = x[1];
+       meas[2]  = x[2];
+       meas[3]  = x[3];
+       meas[4]  = x[4];
+       meas[5]  = c[0];
+       meas[6]  = c[1];
+       meas[7]  = c[2];
+       meas[8]  = c[3];
+       meas[9]  = c[4];
+       meas[10] = c[5];
+       meas[11] = c[6];
+       meas[12] = c[7];
+       meas[13] = c[8];
+       meas[14] = c[9];
+       meas[15] = c[10];
+       meas[16] = c[11];
+       meas[17] = c[12];
+       meas[18] = c[13];
+       meas[19] = c[14];
+       meas[20] = chisq;
+       //cout << "Chisq = " << chisq << endl;
+       /*cout << "track " << iTrack << endl; 
+       cout << "x = " << meas[0] << ", " << meas[1] << ", " << meas[2] 
+            << ", " << meas[3] << ", " << meas[4] << ", " << endl;
+       cout << "c = " << endl 
+            << meas[5] << endl
+            << meas[6] << " " << meas[7] << endl
+            << meas[8] << " " << meas[9] << " " << meas[10] << endl 
+            << meas[11] << " " << meas[12] << " " << meas[13] 
+            << " " << meas[14] << endl
+            << meas[15] << " " << meas[16] << " " << meas[17] << " " 
+            << meas[18] << " " << meas[19] << endl;
+            cout << "chisq = " << meas[20] << endl;*/ 
+       //cout << endl;
+       // Add the track to the trackarray      
+       fKalmanTracks->AddLast(kalmantrack);
+       
+       // Fill the ntuple with the state vector, covariance matrix and
+       // chisquare
+       kalmanTree->Fill(meas);
+      }
+
+      delete track;
+      delete kalmantrack;
+    }
+
+  fBenchmark->Stop("Process tracks");
+  cpuTime = GetCpuTime() - initTime;
+  LOG(AliL3Log::kInformational,"AliL3Kalman::ProcessTracks()","Timing")
+    <<"Process tracks in "<<cpuTime*1000<<" ms"<<ENDLOG;
+
+  // Write tracks to binary file
+  if (fWriteOut)
+    {
+      Char_t tname[80];
+      sprintf(tname,"%s/kalmantracks_%d.raw",fWriteOutPath,fEvent);
+      AliL3MemHandler *mem = new AliL3MemHandler();
+      mem->SetBinaryOutput(tname);
+      mem->TrackArray2Binary(fKalmanTracks);
+      mem->CloseBinaryOutput();
+      delete mem;
+    }
+  
+  TFile *out = new TFile("out.root","recreate");      
+  kalmanTree->Write();
+  out->Close();
+
+  delete kalmanTree;
+}
+
+Int_t AliL3Kalman::MakeSeed(AliL3KalmanTrack *kalmantrack, AliL3Track *track)
+{  
+  // Makes a rough state vector and covariance matrix based on the first
+  // space point in the outmost row
+  /*Double_t initTime,cpuTime;
+  initTime = GetCpuTime();
+
+  fBenchmark->Start("Make seed");*/
+
+  UInt_t *hitnum = track->GetHitNumbers();
+  UInt_t id1, id2, id3;
+
+  // Should do something to make sure that points1 2 and 3 are really not empty
+  id1 = hitnum[0];
+  Int_t slice1 = (id1>>25) & 0x7f;
+  Int_t patch1 = (id1>>22) & 0x7;      
+  UInt_t pos1 = id1&0x3fffff;
+  AliL3SpacePointData *points1 = fClusters[slice1][patch1];
+  
+  id2 = hitnum[1];
+  Int_t slice2 = (id2>>25) & 0x7f;
+  Int_t patch2 = (id2>>22) & 0x7;      
+  UInt_t pos2 = id2&0x3fffff;
+  AliL3SpacePointData *points2 = fClusters[slice2][patch2];
+  if (!points2) return 0;
+
+  id3 = hitnum[2];
+  Int_t slice3 = (id3>>25) & 0x7f;
+  Int_t patch3 = (id3>>22) & 0x7;      
+  UInt_t pos3 = id3&0x3fffff;
+  AliL3SpacePointData *points3 = fClusters[slice3][patch3];
+
+  /*fBenchmark->Stop("Make seed");
+  cpuTime = GetCpuTime() - initTime;
+  LOG(AliL3Log::kInformational,"AliL3Kalman::MakeSeed()","Timing")
+  <<"Make seeds in "<<cpuTime*1000<<" ms"<<ENDLOG;*/
+
+  return   kalmantrack->MakeTrackSeed(points1,pos1,points2,pos2,points3,pos3);
+
+}
+
+Int_t AliL3Kalman::Propagate(AliL3KalmanTrack *kalmantrack, AliL3Track *track)
+{
+  // This function should propagte the track to the next layer thera's a 
+  // cluster.
+  // Must do the following steps,
+  // 1. Go to next spacepoint
+  // 2. Must find the layer position (old and new), predict new state vector,
+  //    and covariance matrix (See AliTPCtrack::PropagateTo
+  // 3. Call Update function
+  /*Double_t initTime,cpuTime;
+  initTime = GetCpuTime();
+  fBenchmark->Start("Propagate");*/
+
+  Int_t num_of_clusters = track->GetNumberOfPoints();
+  
+  UInt_t *hitnum = track->GetHitNumbers();
+  UInt_t id;
+
+  for (Int_t icl = 1; icl < num_of_clusters; icl++)
+    {
+      id = hitnum[icl];
+      Int_t slice = (id>>25) & 0x7f;
+      Int_t patch = (id>>22) & 0x7;    
+      UInt_t pos = id&0x3fffff;
+      AliL3SpacePointData *points = fClusters[slice][patch];
+      if (!points) continue;
+      if (kalmantrack->Propagate(points,pos) == 0) return 0;
+    }
+
+  /*fBenchmark->Stop("Propagate");
+  cpuTime = GetCpuTime() - initTime;
+  LOG(AliL3Log::kInformational,"AliL3Kalman::Propagate()","Timing")
+    <<"Propagate in "<<cpuTime*1000<<" ms"<<ENDLOG;*/
+
+  return 1;
+}
+
+void AliL3Kalman::WriteKalmanTrack(AliL3KalmanTrack *kalman)
+{
+  /*  AliL3MemHandler *memory = new AliL3MemHandler();
+  memory->SetBinaryOutput(filename);
+  if(opt=='a'||opt=='i'){  //add intracks
+    for(Int_t i=0;i<merger->GetNIn();i++){
+      AliL3TrackArray *tr=merger->GetInTracks(i);
+      memory->TrackArray2Binary(tr);
+    }
+  }
+
+  if(opt=='o'||opt=='a'){
+    AliL3TrackArray *tr=merger->GetOutTracks();
+    memory->TrackArray2Binary(tr);
+  }
+
+  memory->CloseBinaryOutput();
+
+  return 1;*/
+}
+
+void Display()
+{
+  
+}
+
+Double_t AliL3Kalman::GetCpuTime()
+{
+  //Return the Cputime in seconds.
+ struct timeval tv;
+ gettimeofday( &tv, NULL );
+ return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
+ //return (Double_t)(clock()) / CLOCKS_PER_SEC;
+}