]>
Commit | Line | Data |
---|---|---|
31f19f09 | 1 | // $Id$ |
2 | ||
3e87ef69 | 3 | #include "AliL3StandardIncludes.h" |
31f19f09 | 4 | #include "AliL3RootTypes.h" |
3e87ef69 | 5 | #include <sys/time.h> |
31f19f09 | 6 | #include <TNtuple.h> |
7 | #include <TTimer.h> | |
8 | #include "AliL3Track.h" | |
9 | #include "AliL3KalmanTrack.h" | |
3e87ef69 | 10 | #include "AliL3Benchmark.h" |
3e87ef69 | 11 | #include "AliL3MemHandler.h" |
31f19f09 | 12 | #include "AliL3FileHandler.h" |
3e87ef69 | 13 | #include "AliL3DataHandler.h" |
14 | #include "AliL3Transform.h" | |
15 | #include "AliL3SpacePointData.h" | |
16 | #include "AliL3DigitData.h" | |
3e87ef69 | 17 | #include "AliL3Logging.h" |
18 | #include "AliL3TrackArray.h" | |
0a86fbb7 | 19 | #include "AliL3TrackSegmentData.h" |
20 | #include "AliL3InterMerger.h" | |
21 | #include "AliL3TrackMerger.h" | |
31f19f09 | 22 | #include "AliL3Kalman.h" |
de3c3890 | 23 | |
3e87ef69 | 24 | /* |
25 | AliL3Kalman | |
26 | */ | |
27 | ||
28 | ClassImp(AliL3Kalman) | |
29 | ||
0a86fbb7 | 30 | AliL3Kalman::AliL3Kalman(Char_t *datapath, Int_t *slice, Int_t min_clusters = 0){ |
3e87ef69 | 31 | // Constructor |
3e87ef69 | 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] = AliL3Transform::GetLastRow(-1); | |
0a86fbb7 | 52 | fWriteOut = kTRUE; |
3e87ef69 | 53 | fBenchmark = 0; |
b2a02bce | 54 | |
3e87ef69 | 55 | } |
56 | ||
57 | AliL3Kalman::~AliL3Kalman() | |
58 | { | |
59 | // Destructor | |
60 | if (fBenchmark) delete fBenchmark; | |
61 | } | |
62 | ||
63 | void AliL3Kalman::Init() | |
64 | { | |
65 | fBenchmark = new AliL3Benchmark(); | |
66 | } | |
67 | ||
68 | void AliL3Kalman::LoadTracks(Int_t event, Bool_t sp) | |
69 | { | |
0a86fbb7 | 70 | // Load space points and tracks from conformal tracker |
3e87ef69 | 71 | // Must also be possible to take seeds (and clusters) from Hough-transform?? |
72 | ||
3e87ef69 | 73 | Double_t initTime,cpuTime; |
74 | initTime = GetCpuTime(); | |
75 | fBenchmark->Start("Load tracks"); | |
76 | ||
0a86fbb7 | 77 | // Load space points |
3e87ef69 | 78 | Char_t fname[1024]; |
79 | AliL3FileHandler *clusterfile[36][6]; | |
80 | for(Int_t s=fMinSlice; s<=fMaxSlice; s++) | |
81 | { | |
82 | for(Int_t p=0; p<AliL3Transform::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 AliL3FileHandler(); | |
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(AliL3Log::kError,"AliL3Evaluation::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] = (AliL3SpacePointData*)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 | |
de3c3890 | 114 | //sprintf(fname,"%s/kalmantracks_%d.raw",fPath,event); |
3e87ef69 | 115 | sprintf(fname,"%s/tracks_%d.raw",fPath,event); |
116 | AliL3FileHandler *tfile = new AliL3FileHandler(); | |
117 | if(!tfile->SetBinaryInput(fname)){ | |
118 | LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open") | |
119 | <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; | |
120 | return; | |
121 | } | |
122 | if(fTracks) | |
123 | delete fTracks; | |
124 | fTracks = new AliL3TrackArray(); | |
125 | tfile->Binary2TrackArray(fTracks); | |
126 | tfile->CloseBinaryInput(); | |
127 | delete tfile; | |
3e87ef69 | 128 | fBenchmark->Stop("Load tracks"); |
129 | cpuTime = GetCpuTime() - initTime; | |
130 | LOG(AliL3Log::kInformational,"AliL3Kalman::LoadTracks()","Timing") | |
131 | <<"Loading tracks in "<<cpuTime*1000<<" ms"<<ENDLOG; | |
132 | ||
133 | } | |
134 | ||
135 | void AliL3Kalman::ProcessTracks() | |
136 | { | |
0a86fbb7 | 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. | |
b2a02bce | 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; | |
3e87ef69 | 144 | |
145 | Double_t initTime,cpuTime; | |
146 | initTime = GetCpuTime(); | |
147 | ||
148 | fBenchmark->Start("Process tracks"); | |
149 | ||
3e87ef69 | 150 | fTracks->QSort(); |
151 | ||
152 | fKalmanTracks = new AliL3TrackArray(); | |
153 | ||
0bd0c1ef | 154 | // Make a ntuple to store state vector, covariance matrix and chisquare |
0a86fbb7 | 155 | // Will eventually not need a TTree?? |
0bd0c1ef | 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"); |
3e87ef69 | 157 | Float_t meas[21]; |
b2a02bce | 158 | |
3e87ef69 | 159 | // Go through the tracks from conformal or hough tracker |
160 | for (Int_t iTrack = 0; iTrack < fTracks->GetNTracks(); iTrack++) | |
161 | { | |
b2a02bce | 162 | /*Double_t initTime,cpuTime; |
163 | initTime = GetCpuTime(); | |
164 | fBenchmark->Start("Process tracks");*/ | |
0a86fbb7 | 165 | |
3e87ef69 | 166 | AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(iTrack); |
167 | if (!track) continue; | |
168 | if (track->GetNumberOfPoints() < fMinPointsOnTrack) continue; | |
169 | ||
b2a02bce | 170 | AliL3KalmanTrack *kalmantrack = new AliL3KalmanTrack(); |
3e87ef69 | 171 | |
b2a02bce | 172 | Bool_t save = kTRUE; |
0bd0c1ef | 173 | |
de3c3890 | 174 | /*if (InitKalmanTrack(kalmantrack, track) == 0) |
b2a02bce | 175 | { |
0bd0c1ef | 176 | save = kFALSE; |
177 | continue; | |
de3c3890 | 178 | }*/ |
0bd0c1ef | 179 | |
de3c3890 | 180 | if (MakeKalmanSeed(kalmantrack,track) ==0) |
181 | { | |
182 | save = kFALSE; | |
183 | continue; | |
184 | } | |
3e87ef69 | 185 | |
186 | if (Propagate(kalmantrack, track) == 0) | |
187 | { | |
188 | save = kFALSE; | |
189 | } | |
de3c3890 | 190 | |
191 | if (save) {// cout << track->GetPt() << endl; | |
3e87ef69 | 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; | |
0a86fbb7 | 218 | |
3e87ef69 | 219 | // Add the track to the trackarray |
0a86fbb7 | 220 | AliL3Track *outtrack = (AliL3Track*)fKalmanTracks->NextTrack(); |
221 | outtrack->Set(track); | |
de3c3890 | 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*AliL3Transform::GetBField()))); | |
227 | //outtrack->SetPt(1/(TMath::Abs(x[4])/0.0029980*AliL3Transform::GetBField())); | |
228 | outtrack->SetPsi(x[2]); | |
229 | //outtrack->Set(track); | |
0a86fbb7 | 230 | |
3e87ef69 | 231 | // Fill the ntuple with the state vector, covariance matrix and |
0bd0c1ef | 232 | // chisquare |
233 | kalmanTree->Fill(meas); | |
3e87ef69 | 234 | } |
235 | ||
236 | delete track; | |
237 | delete kalmantrack; | |
0a86fbb7 | 238 | |
b2a02bce | 239 | /*fBenchmark->Stop("Process tracks"); |
240 | cpuTime = GetCpuTime() - initTime; | |
241 | LOG(AliL3Log::kInformational,"AliL3Kalman::ProcessTracks()","Timing") | |
242 | <<"Processed track "<<iTrack<<" in "<<cpuTime*1000<<" ms"<<ENDLOG;*/ | |
3e87ef69 | 243 | } |
244 | ||
245 | fBenchmark->Stop("Process tracks"); | |
246 | cpuTime = GetCpuTime() - initTime; | |
247 | LOG(AliL3Log::kInformational,"AliL3Kalman::ProcessTracks()","Timing") | |
248 | <<"Process tracks in "<<cpuTime*1000<<" ms"<<ENDLOG; | |
0a86fbb7 | 249 | |
3e87ef69 | 250 | if (fWriteOut) |
251 | { | |
252 | Char_t tname[80]; | |
253 | sprintf(tname,"%s/kalmantracks_%d.raw",fWriteOutPath,fEvent); | |
254 | AliL3MemHandler *mem = new AliL3MemHandler(); | |
255 | mem->SetBinaryOutput(tname); | |
256 | mem->TrackArray2Binary(fKalmanTracks); | |
257 | mem->CloseBinaryOutput(); | |
258 | delete mem; | |
259 | } | |
b2a02bce | 260 | |
0a86fbb7 | 261 | // This will be removed?? |
262 | TFile *out = new TFile("kalmantracks.root","recreate"); | |
3e87ef69 | 263 | kalmanTree->Write(); |
264 | out->Close(); | |
265 | ||
266 | delete kalmanTree; | |
0a86fbb7 | 267 | |
3e87ef69 | 268 | } |
269 | ||
de3c3890 | 270 | Int_t AliL3Kalman::MakeKalmanSeed(AliL3KalmanTrack *kalmantrack, AliL3Track *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 | AliL3SpacePointData *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 | AliL3SpacePointData *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 | AliL3SpacePointData *points2 = fClusters[slice2][patch2]; | |
294 | ||
295 | return kalmantrack->MakeSeed(track, points0, pos0, slice0, points1, pos1, slice1, points2, pos2, slice2); | |
296 | } | |
297 | ||
b2a02bce | 298 | Int_t AliL3Kalman::InitKalmanTrack(AliL3KalmanTrack *kalmantrack, AliL3Track *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 | AliL3SpacePointData *points = fClusters[slice][patch]; | |
308 | ||
0bd0c1ef | 309 | return kalmantrack->Init(track, points, pos, slice); |
3e87ef69 | 310 | } |
311 | ||
312 | Int_t AliL3Kalman::Propagate(AliL3KalmanTrack *kalmantrack, AliL3Track *track) | |
313 | { | |
0a86fbb7 | 314 | // This function propagtes the kalmantrack to the next cluster of the loaded |
315 | // track | |
3e87ef69 | 316 | Int_t num_of_clusters = track->GetNumberOfPoints(); |
3e87ef69 | 317 | UInt_t *hitnum = track->GetHitNumbers(); |
318 | UInt_t id; | |
0a86fbb7 | 319 | UInt_t badpoint = 0; |
3e87ef69 | 320 | |
321 | for (Int_t icl = 1; icl < num_of_clusters; icl++) | |
322 | { | |
0a86fbb7 | 323 | |
3e87ef69 | 324 | id = hitnum[icl]; |
325 | Int_t slice = (id>>25) & 0x7f; | |
b2a02bce | 326 | Int_t patch = (id>>22) & 0x7; |
3e87ef69 | 327 | UInt_t pos = id&0x3fffff; |
b2a02bce | 328 | |
3e87ef69 | 329 | AliL3SpacePointData *points = fClusters[slice][patch]; |
330 | if (!points) continue; | |
0bd0c1ef | 331 | if (kalmantrack->Propagate(points,pos,slice) == 0) |
0a86fbb7 | 332 | { |
333 | badpoint++; | |
334 | continue; | |
335 | } | |
3e87ef69 | 336 | } |
b2a02bce | 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; | |
3e87ef69 | 341 | } |
342 | ||
343 | Double_t AliL3Kalman::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 | } |