HLT/kalman/AliHLTKalman.cxx

   1 // $Id$
   2
   3 #include "AliHLTStandardIncludes.h"
   4 #include "AliHLTRootTypes.h"
   5 #include <sys/time.h>
   6 #include <TNtuple.h>
   7 #include <TTimer.h>
   8 #include "AliHLTTrack.h"
   9 #include "AliHLTKalmanTrack.h"
  10 #include "AliHLTBenchmark.h"
  11 #include "AliHLTMemHandler.h"
  12 #include "AliHLTFileHandler.h"
  13 #include "AliHLTDataHandler.h"
  14 #include "AliHLTTransform.h"
  15 #include "AliHLTSpacePointData.h"
  16 #include "AliHLTDigitData.h"
  17 #include "AliHLTLogging.h"
  18 #include "AliHLTTrackArray.h"
  19 #include "AliHLTTrackSegmentData.h"
  20 #include "AliHLTInterMerger.h"
  21 #include "AliHLTTrackMerger.h"
  22 #include "AliHLTKalman.h"
  23
  24 /*
  25   AliHLTKalman
  26 */
  27
  28 ClassImp(AliHLTKalman)
  29
  30 AliHLTKalman::AliHLTKalman(Char_t *datapath, Int_t *slice, Int_t min_clusters = 0){
  31   // Constructor
  32   if (slice)
  33     {
  34       fMinSlice = slice[0];
  35       fMaxSlice = slice[1];
  36     }
  37   else
  38     {
  39       fMinSlice = 0;
  40       fMaxSlice = 35;
  41     }
  42
  43   sprintf(fPath,"%s",datapath);
  44   fMinPointsOnTrack = 0; //min_clusters;
  45   fTracks = 0;
  46   fKalmanTracks = 0;
  47
  48   // NB! fNrow under only for single-patch, must include other possibilities
  49   // later on. ?? Maybe better also to put it in an Init-function
  50   fRow[0][0] = 0;
  51   fRow[0][1] = AliHLTTransform::GetLastRow(-1);
  52   fWriteOut = kTRUE;
  53   fBenchmark = 0;
  54
  55 }
  56
  57 AliHLTKalman::~AliHLTKalman()
  58 {
  59   // Destructor
  60   if (fBenchmark) delete fBenchmark;
  61 }
  62
  63 void AliHLTKalman::Init()
  64 {
  65   fBenchmark = new AliHLTBenchmark();
  66 }
  67
  68 void AliHLTKalman::LoadTracks(Int_t event, Bool_t sp)
  69 {
  70   // Load space points and tracks from conformal tracker
  71   // Must also be possible to take seeds (and clusters) from Hough-transform??
  72
  73   Double_t initTime,cpuTime;
  74   initTime = GetCpuTime();
  75   fBenchmark->Start("Load tracks");
  76
  77   // Load space points
  78   Char_t fname[1024];
  79   AliHLTFileHandler *clusterfile[36][6];
  80   for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
  81     {
  82       for(Int_t p=0; p<AliHLTTransform::GetNPatches(); p++)
  83         {
  84           Int_t patch;
  85           if(sp==kTRUE)
  86             patch=-1;
  87           else
  88             patch=p;
  89
  90           fClusters[s][p] = 0;
  91           clusterfile[s][p] = new AliHLTFileHandler();
  92           if(event<0)
  93             sprintf(fname,"%s/points_%d_%d.raw",fPath,s,patch);
  94           else
  95             sprintf(fname,"%s/points_%d_%d_%d.raw",fPath,event,s,patch);
  96           if(!clusterfile[s][p]->SetBinaryInput(fname))
  97             {
  98               LOG(AliHLTLog::kError,"AliHLTEvaluation::Setup","File Open")
  99                 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 100               delete clusterfile[s][p];
 101               clusterfile[s][p] = 0;
 102               continue;
 103             }
 104           fClusters[s][p] = (AliHLTSpacePointData*)clusterfile[s][p]->Allocate();
 105           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
 106           clusterfile[s][p]->CloseBinaryInput();
 107           if(sp==kTRUE)
 108             break;
 109           delete clusterfile[s][p];
 110         }
 111     }
 112
 113   // Load tracks
 114   //sprintf(fname,"%s/kalmantracks_%d.raw",fPath,event);
 115   sprintf(fname,"%s/tracks_%d.raw",fPath,event);
 116   AliHLTFileHandler *tfile = new AliHLTFileHandler();
 117   if(!tfile->SetBinaryInput(fname)){
 118     LOG(AliHLTLog::kError,"AliHLTEvaluation::Setup","File Open")
 119       <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
 120     return;
 121   }
 122   if(fTracks)
 123     delete fTracks;
 124   fTracks = new AliHLTTrackArray();
 125   tfile->Binary2TrackArray(fTracks);
 126   tfile->CloseBinaryInput();
 127   delete tfile;
 128   fBenchmark->Stop("Load tracks");
 129   cpuTime = GetCpuTime() - initTime;
 130   LOG(AliHLTLog::kInformational,"AliHLTKalman::LoadTracks()","Timing")
 131     <<"Loading tracks in "<<cpuTime*1000<<" ms"<<ENDLOG;
 132
 133 }
 134
 135 void AliHLTKalman::ProcessTracks()
 136 {
 137   // Run the Kalman filter algorithm on the loaded tracks.
 138   // If the track is OK, the loaded track is saved in file kalmantracks_0.raw
 139   // The kalman filter variables (that is the state vector, covariance matrix
 140   // and chi2) is written to root-file kalmantracks.root.
 141   // Should be extended to correct the track parameters of the loaded track
 142   // or save the kalmantrack instead of the loaded track.??
 143   UInt_t fEvent = 0;
 144
 145   Double_t initTime,cpuTime;
 146   initTime = GetCpuTime();
 147
 148   fBenchmark->Start("Process tracks");
 149
 150   fTracks->QSort();
 151
 152   fKalmanTracks = new AliHLTTrackArray();
 153
 154   // Make a ntuple to store state vector, covariance matrix and chisquare
 155   // Will eventually not need a TTree??
 156   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");
 157   Float_t meas[21];
 158
 159   // Go through the tracks from conformal or hough tracker
 160   for (Int_t iTrack = 0; iTrack < fTracks->GetNTracks(); iTrack++)
 161     {
 162       /*Double_t initTime,cpuTime;
 163       initTime = GetCpuTime();
 164       fBenchmark->Start("Process tracks");*/
 165
 166       AliHLTTrack *track = (AliHLTTrack*)fTracks->GetCheckedTrack(iTrack);
 167       if (!track) continue;
 168       if (track->GetNumberOfPoints() < fMinPointsOnTrack) continue;
 169
 170       AliHLTKalmanTrack *kalmantrack = new AliHLTKalmanTrack();
 171
 172       Bool_t save = kTRUE;
 173
 174       /*if (InitKalmanTrack(kalmantrack, track) == 0)
 175         {
 176           save = kFALSE;
 177           continue;
 178           }*/
 179
 180       if (MakeKalmanSeed(kalmantrack,track) ==0)
 181         {
 182           save = kFALSE;
 183           continue;
 184         }
 185
 186       if (Propagate(kalmantrack, track) == 0)
 187         {
 188           save = kFALSE;
 189         }
 190
 191       if (save) {// cout << track->GetPt() << endl;
 192         Float_t x[5];
 193         kalmantrack->GetStateVector(x);
 194         Float_t c[15];
 195         kalmantrack->GetCovariance(c);
 196         Float_t chisq = kalmantrack->GetChisq();
 197         meas[0]  = x[0];
 198         meas[1]  = x[1];
 199         meas[2]  = x[2];
 200         meas[3]  = x[3];
 201         meas[4]  = x[4];
 202         meas[5]  = c[0];
 203         meas[6]  = c[1];
 204         meas[7]  = c[2];
 205         meas[8]  = c[3];
 206         meas[9]  = c[4];
 207         meas[10] = c[5];
 208         meas[11] = c[6];
 209         meas[12] = c[7];
 210         meas[13] = c[8];
 211         meas[14] = c[9];
 212         meas[15] = c[10];
 213         meas[16] = c[11];
 214         meas[17] = c[12];
 215         meas[18] = c[13];
 216         meas[19] = c[14];
 217         meas[20] = chisq;
 218
 219         // Add the track to the trackarray
 220         AliHLTTrack *outtrack = (AliHLTTrack*)fKalmanTracks->NextTrack();
 221         outtrack->Set(track);
 222         // SET THE PARAMETERS ACCORDING TO KALMAN FILTER
 223         outtrack->SetTgl(x[3]);
 224         // The factor 2 in the expression for Pt is not included in the similar offline expression. However
 225         // it should be like this if I use a factor 1/2 in the calculation of par4??
 226         //outtrack->SetPt(1/(2*TMath::Abs(1e-9*TMath::Abs(x[4])/x[4] + x[4])/(0.0029980*AliHLTTransform::GetBField())));
 227         //outtrack->SetPt(1/(TMath::Abs(x[4])/0.0029980*AliHLTTransform::GetBField()));
 228         outtrack->SetPsi(x[2]);
 229         //outtrack->Set(track);
 230
 231         // Fill the ntuple with the state vector, covariance matrix and
 232         // chisquare
 233         kalmanTree->Fill(meas);
 234       }
 235
 236       delete track;
 237       delete kalmantrack;
 238
 239       /*fBenchmark->Stop("Process tracks");
 240       cpuTime = GetCpuTime() - initTime;
 241       LOG(AliHLTLog::kInformational,"AliHLTKalman::ProcessTracks()","Timing")
 242       <<"Processed track "<<iTrack<<" in "<<cpuTime*1000<<" ms"<<ENDLOG;*/
 243     }
 244
 245   fBenchmark->Stop("Process tracks");
 246   cpuTime = GetCpuTime() - initTime;
 247   LOG(AliHLTLog::kInformational,"AliHLTKalman::ProcessTracks()","Timing")
 248     <<"Process tracks in "<<cpuTime*1000<<" ms"<<ENDLOG;
 249
 250   if (fWriteOut)
 251     {
 252       Char_t tname[80];
 253       sprintf(tname,"%s/kalmantracks_%d.raw",fWriteOutPath,fEvent);
 254       AliHLTMemHandler *mem = new AliHLTMemHandler();
 255       mem->SetBinaryOutput(tname);
 256       mem->TrackArray2Binary(fKalmanTracks);
 257       mem->CloseBinaryOutput();
 258       delete mem;
 259     }
 260
 261   // This will be removed??
 262   TFile *out = new TFile("kalmantracks.root","recreate");
 263   kalmanTree->Write();
 264   out->Close();
 265
 266   delete kalmanTree;
 267
 268 }
 269
 270 Int_t AliHLTKalman::MakeKalmanSeed(AliHLTKalmanTrack *kalmantrack, AliHLTTrack *track)
 271 {
 272   Int_t num_of_clusters = track->GetNumberOfPoints();
 273
 274   UInt_t *hitnum = track->GetHitNumbers();
 275   UInt_t id;
 276
 277   id = hitnum[0];
 278   Int_t slice0 = (id>>25) & 0x7f;
 279   Int_t patch0 = (id>>22) & 0x7;
 280   UInt_t pos0 = id&0x3fffff;
 281   AliHLTSpacePointData *points0 = fClusters[slice0][patch0];
 282
 283   id = hitnum[Int_t(num_of_clusters/2)];
 284   Int_t slice1 = (id>>25) & 0x7f;
 285   Int_t patch1 = (id>>22) & 0x7;
 286   UInt_t pos1 = id&0x3fffff;
 287   AliHLTSpacePointData *points1 = fClusters[slice1][patch1];
 288
 289   id = hitnum[num_of_clusters-1];
 290   Int_t slice2 = (id>>25) & 0x7f;
 291   Int_t patch2 = (id>>22) & 0x7;
 292   UInt_t pos2 = id&0x3fffff;
 293   AliHLTSpacePointData *points2 = fClusters[slice2][patch2];
 294
 295   return kalmantrack->MakeSeed(track, points0, pos0, slice0, points1, pos1, slice1, points2, pos2, slice2);
 296 }
 297
 298 Int_t AliHLTKalman::InitKalmanTrack(AliHLTKalmanTrack *kalmantrack, AliHLTTrack *track)
 299 {
 300   UInt_t *hitnum = track->GetHitNumbers();
 301   UInt_t id;
 302
 303   id = hitnum[0];
 304   Int_t slice = (id>>25) & 0x7f;
 305   Int_t patch = (id>>22) & 0x7;
 306   UInt_t pos = id&0x3fffff;
 307   AliHLTSpacePointData *points = fClusters[slice][patch];
 308
 309   return kalmantrack->Init(track, points, pos, slice);
 310 }
 311
 312 Int_t AliHLTKalman::Propagate(AliHLTKalmanTrack *kalmantrack, AliHLTTrack *track)
 313 {
 314   // This function propagtes the kalmantrack to the next cluster of the loaded
 315   // track
 316   Int_t num_of_clusters = track->GetNumberOfPoints();
 317   UInt_t *hitnum = track->GetHitNumbers();
 318   UInt_t id;
 319   UInt_t badpoint = 0;
 320
 321   for (Int_t icl = 1; icl < num_of_clusters; icl++)
 322     {
 323
 324       id = hitnum[icl];
 325       Int_t slice = (id>>25) & 0x7f;
 326       Int_t patch = (id>>22) & 0x7;
 327       UInt_t pos = id&0x3fffff;
 328
 329       AliHLTSpacePointData *points = fClusters[slice][patch];
 330       if (!points) continue;
 331       if (kalmantrack->Propagate(points,pos,slice) == 0)
 332         {
 333           badpoint++;
 334           continue;
 335         }
 336     }
 337
 338   // If too many clusters are missing, the track is no good
 339   if (badpoint >= UInt_t(num_of_clusters*0.8)) return 0;
 340   return 1;
 341 }
 342
 343 Double_t AliHLTKalman::GetCpuTime()
 344 {
 345   //Return the Cputime in seconds.
 346  struct timeval tv;
 347  gettimeofday( &tv, NULL );
 348  return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
 349  //return (Double_t)(clock()) / CLOCKS_PER_SEC;
 350 }