New version of SPD raw-data reconstruction. The format now correponds to the actual...
[u/mrichter/AliRoot.git] / 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 }