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