45a59de3aec9ef8524864d89fc3dcea89a35a7ba
[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
25 /*
26   AliL3Kalman
27 */
28
29 ClassImp(AliL3Kalman)
30
31 AliL3Kalman::AliL3Kalman(Char_t *datapath, Int_t *slice, Int_t min_clusters = 0){
32   // Constructor
33
34   Int_t fEvent = 0;
35
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 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   //cout << "Number of loaded tracks " << fTracks->GetNTracks() << endl;
131   fBenchmark->Stop("Load tracks");
132   cpuTime = GetCpuTime() - initTime;
133   LOG(AliL3Log::kInformational,"AliL3Kalman::LoadTracks()","Timing")
134     <<"Loading tracks in "<<cpuTime*1000<<" ms"<<ENDLOG;
135
136 }
137
138 void AliL3Kalman::ProcessTracks()
139 {
140   // Run the Kalman filter algorithm on the loaded tracks. 
141   // If the track is OK, the loaded track is saved in file kalmantracks_0.raw
142   // The kalman filter variables (that is the state vector, covariance matrix 
143   // and chi2) is written to root-file kalmantracks.root. 
144
145   Double_t initTime,cpuTime;
146   initTime = GetCpuTime();
147
148   fBenchmark->Start("Process tracks");
149
150   fTracks->QSort();
151
152   fKalmanTracks = new AliL3TrackArray();
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       AliL3KalmanTrack *kalmantrack = new AliL3KalmanTrack();
163       kalmantrack->Init();
164
165       AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(iTrack);
166       if (!track) continue;
167       if (track->GetNumberOfPoints() < fMinPointsOnTrack) continue;    
168
169       Bool_t save = kTRUE;
170
171       if (MakeSeed(kalmantrack, track) == 0)
172         {
173           save = kFALSE;
174           continue;
175         }
176
177       if (Propagate(kalmantrack, track) == 0) 
178         {
179           save = kFALSE;
180         }
181
182       if (save) {
183         Float_t x[5]; 
184         kalmantrack->GetStateVector(x);
185         Float_t c[15]; 
186         kalmantrack->GetCovariance(c);
187         Float_t chisq = kalmantrack->GetChisq();
188         meas[0]  = x[0];
189         meas[1]  = x[1];
190         meas[2]  = x[2];
191         meas[3]  = x[3];
192         meas[4]  = x[4];
193         meas[5]  = c[0];
194         meas[6]  = c[1];
195         meas[7]  = c[2];
196         meas[8]  = c[3];
197         meas[9]  = c[4];
198         meas[10] = c[5];
199         meas[11] = c[6];
200         meas[12] = c[7];
201         meas[13] = c[8];
202         meas[14] = c[9];
203         meas[15] = c[10];
204         meas[16] = c[11];
205         meas[17] = c[12];
206         meas[18] = c[13];
207         meas[19] = c[14];
208         meas[20] = chisq;
209
210         // Add the track to the trackarray      
211         AliL3Track *outtrack = (AliL3Track*)fKalmanTracks->NextTrack();
212         outtrack->Set(track);
213
214         // Fill the ntuple with the state vector, covariance matrix and
215         // chisquare
216         kalmanTree->Fill(meas);
217       }
218
219       delete track;
220       delete kalmantrack;
221
222     }
223
224   //WriteKalmanTracks();
225   fBenchmark->Stop("Process tracks");
226   cpuTime = GetCpuTime() - initTime;
227   LOG(AliL3Log::kInformational,"AliL3Kalman::ProcessTracks()","Timing")
228     <<"Process tracks in "<<cpuTime*1000<<" ms"<<ENDLOG;
229   
230   if (fWriteOut)
231     {
232       Char_t tname[80];
233       sprintf(tname,"%s/kalmantracks_%d.raw",fWriteOutPath,fEvent);
234       AliL3MemHandler *mem = new AliL3MemHandler();
235       mem->SetBinaryOutput(tname);
236       mem->TrackArray2Binary(fKalmanTracks);
237       mem->CloseBinaryOutput();
238       delete mem;
239     }
240
241   // This will be removed??
242   TFile *out = new TFile("kalmantracks.root","recreate");      
243   kalmanTree->Write();
244   out->Close();
245
246   delete kalmanTree;
247   
248 }
249
250 Int_t AliL3Kalman::MakeSeed(AliL3KalmanTrack *kalmantrack, AliL3Track *track)
251 {  
252   // Makes a rough state vector and covariance matrix based on three
253   // space points of the loaded track 
254  
255   UInt_t *hitnum = track->GetHitNumbers();
256   UInt_t id1, id2, id3;
257
258   id1 = hitnum[0];
259   Int_t slice1 = (id1>>25) & 0x7f;
260   Int_t patch1 = (id1>>22) & 0x7;       
261   UInt_t pos1 = id1&0x3fffff;
262   AliL3SpacePointData *points1 = fClusters[slice1][patch1];
263   
264   id2 = hitnum[Int_t(track->GetNHits()/2)];
265   Int_t slice2 = (id2>>25) & 0x7f;
266   Int_t patch2 = (id2>>22) & 0x7;       
267   UInt_t pos2 = id2&0x3fffff;
268   AliL3SpacePointData *points2 = fClusters[slice2][patch2];
269
270   id3 = hitnum[track->GetNHits()-1];
271   Int_t slice3 = (id3>>25) & 0x7f;
272   Int_t patch3 = (id3>>22) & 0x7;       
273   UInt_t pos3 = id3&0x3fffff;
274   AliL3SpacePointData *points3 = fClusters[slice3][patch3];
275
276   return kalmantrack->MakeTrackSeed(points1,pos1,points2,pos2,points3,pos3);
277
278 }
279
280 Int_t AliL3Kalman::Propagate(AliL3KalmanTrack *kalmantrack, AliL3Track *track)
281 {
282   // This function propagtes the kalmantrack to the next cluster of the loaded 
283   // track 
284
285   Int_t num_of_clusters = track->GetNumberOfPoints();
286   
287   UInt_t *hitnum = track->GetHitNumbers();
288   UInt_t id;
289   UInt_t badpoint = 0;
290
291   for (Int_t icl = 1; icl < num_of_clusters; icl++)
292     {
293
294       id = hitnum[icl];
295       Int_t slice = (id>>25) & 0x7f;
296       Int_t patch = (id>>22) & 0x7;     
297       UInt_t pos = id&0x3fffff;
298       AliL3SpacePointData *points = fClusters[slice][patch];
299       if (!points) continue;
300       // Go to next cluster, but don't update track!!??
301       // But if there are several clusters that are no good: return 0;
302       if (kalmantrack->Propagate(points,pos) == 0) 
303         {
304           badpoint++;
305           continue;
306         }
307       //if (kalmantrack->Propagate(points,pos) == 0) return 0;
308     }
309
310   if (badpoint >= UInt_t(num_of_clusters/2)) return 0;
311   else return 1;
312 }
313
314 Double_t AliL3Kalman::GetCpuTime()
315 {
316   //Return the Cputime in seconds.
317  struct timeval tv;
318  gettimeofday( &tv, NULL );
319  return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
320  //return (Double_t)(clock()) / CLOCKS_PER_SEC;
321 }