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