Added support for NEWIO, merged cern-hlt tree, updated to latest track candidate...
[u/mrichter/AliRoot.git] / HLT / comp / AliL3OfflineDataCompressor.cxx
CommitLineData
0a86fbb7 1// @(#) $Id$
2
3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4//*-- Copyright &copy ALICE HLT Group
5
6#include "AliL3StandardIncludes.h"
7
0a86fbb7 8#include <AliTPCParamSR.h>
9#include <AliTPCClustersArray.h>
10#include <AliTPCcluster.h>
11#include <AliTPCClustersRow.h>
12#include <AliTPC.h>
13#include <AliTPCv2.h>
b2a02bce 14#include <AliTPCcluster.h>
15#include <AliTPCtracker.h>
16#include <AliTPCclusterMI.h>
17#include <AliTPCtrackerMI.h>
18#include <AliKalmanTrack.h>
0a86fbb7 19#include <AliRun.h>
20
21#include <TTree.h>
22#include <TFile.h>
23#include <TMath.h>
24#include <TDirectory.h>
25#include <TSystem.h>
26#include <TH2F.h>
27
b2a02bce 28#include "AliL3Transform.h"
29#include "AliL3ModelTrack.h"
30#include "AliL3Compress.h"
31#include "AliL3TrackArray.h"
0a86fbb7 32
33#include "AliL3OfflineDataCompressor.h"
34
35#if GCCVERSION == 3
36using namespace std;
37#endif
38
39//_____________________________________________________________
40//
41// AliL3OfflineDataCompression
42//
43
44
45ClassImp(AliL3OfflineDataCompressor)
46
47AliL3OfflineDataCompressor::AliL3OfflineDataCompressor()
48{
49 fMarian = kFALSE;
50 fTracker=0;
51}
52
53AliL3OfflineDataCompressor::AliL3OfflineDataCompressor(Char_t *path,Bool_t keep,Bool_t writeshape,Bool_t MI)
54 : AliL3DataCompressor(path,keep,writeshape)
55{
56 fMarian = MI;
57 fTracker=0;
58}
59
60AliL3OfflineDataCompressor::~AliL3OfflineDataCompressor()
61{
62 if(fTracker)
63 {
64 fTracker->UnloadClusters();
65 delete fTracker;
66 }
67}
68
69
70void AliL3OfflineDataCompressor::LoadData(Int_t event,Bool_t sp)
71{
72 //Take offline reconstructed tracks as an input.
73 //In this case, no remaining clusters are written.
74
75 fSinglePatch = sp;
76
77 char filename[1024];
78 AliKalmanTrack::SetConvConst(1000/0.299792458/AliL3Transform::GetSolenoidField());
79 if(fMarian==kFALSE)
80 sprintf(filename,"%s/offline/AliTPCclusters.root",fPath);
81 else
82 sprintf(filename,"%s/offline/AliTPCclustersMI.root",fPath);
83
84 TFile *in = TFile::Open(filename);
85 AliTPCParam *param=(AliTPCParam*)in->Get("75x40_100x60_150x60");
86
87 if(fMarian==kFALSE)
88 fTracker = new AliTPCtracker(param);
89 else
90 fTracker = new AliTPCtrackerMI(param);
91 fTracker->SetEventNumber(event);
b2a02bce 92#ifdef asvversion
93 fTracker->LoadClusters();
94#endif
0a86fbb7 95 if(fMarian==kTRUE)
96 {
97 AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
b2a02bce 98#ifdef asvversion
0a86fbb7 99 mitracker->LoadInnerSectors();
100 mitracker->LoadOuterSectors();
b2a02bce 101#endif
0a86fbb7 102 }
103
104 const Int_t MAX=20000;
105 Int_t nentr=0,i=0; TObjArray tarray(MAX);
106 if(fMarian==kFALSE)
107 sprintf(filename,"%s/offline/AliTPCtracks.root",fPath);
108 else
109 sprintf(filename,"%s/offline/AliTPCtracksMI.root",fPath);
110 TFile *tf=TFile::Open(filename);
111
112 char tname[100]; sprintf(tname,"TreeT_TPC_%d",event);
113 TTree *tracktree=(TTree*)tf->Get(tname);
114
115 TBranch *tbranch=tracktree->GetBranch("tracks");
116 nentr=(Int_t)tracktree->GetEntries();
117 AliTPCtrack *iotrack=0;
118
119 for (i=0; i<nentr; i++) {
120 iotrack=new AliTPCtrack;
121 tbranch->SetAddress(&iotrack);
122 tracktree->GetEvent(i);
123 tarray.AddLast(iotrack);
124 }
125 delete tracktree;
126 tf->Close();
127
128 AliL3TrackArray *comptracks = new AliL3TrackArray("AliL3ModelTrack");
129 cout<<"Loaded "<<nentr<<" offline tracks"<<endl;
130 Int_t slice,padrow;
131 Int_t totcounter=0;
132 for(i=0; i<nentr; i++)
133 {
134
135 AliTPCtrack *track=(AliTPCtrack*)tarray.UncheckedAt(i);
136 Int_t nhits = track->GetNumberOfClusters();
137 Int_t idx = track->GetClusterIndex(nhits-1);
138 Int_t sec=(idx&0xff000000)>>24, row=(idx&0x00ff0000)>>16;
139
140 /*
141 TPC sector numbering within the AliTPCtracker class:
142 There are in total 18 inner sectors and 18 outer sectors.
143 This means that one sector includes _both_ sides of the TPC.
144 Example: sec=0 -> sector 0 and 18.
145 sec=18 -> sector 36 and 54
146 */
147
148 if(fMarian==kFALSE)
149 if(sec >= 18)
150 sec += 18;
151
152 AliL3Transform::Sector2Slice(slice,padrow,sec,row);
153 Double_t par[5],xk=AliL3Transform::Row2X(padrow);
154 track->PropagateTo(xk);
155 track->GetExternalParameters(xk,par);
156 Double_t psi = TMath::ASin(par[2]) + track->GetAlpha();
157 if (psi<-TMath::Pi()) psi+=2*TMath::Pi();
158 if (psi>=TMath::Pi()) psi-=2*TMath::Pi();
159 Float_t pt_1=TMath::Abs(par[4]);
160 Int_t charge = 1;
161 if(par[4] > 0)
162 charge=-1;
163
164 Float_t first[3];
165 AliCluster *fcl=0;
166 if(fMarian==kFALSE)
167 fcl= fTracker->GetCluster(idx);
168 else
169 {
170 AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
171 fcl = mitracker->GetClusterMI(idx);
172 }
173 first[0] = xk;
174 first[1] = fcl->GetY();
175 first[2] = fcl->GetZ();
176
177 AliL3Transform::Local2Global(first,slice);
178
179 AliL3ModelTrack *outtrack = (AliL3ModelTrack*)comptracks->NextTrack();
180 outtrack->SetNHits(nhits);
181 outtrack->SetFirstPoint(first[0],first[1],first[2]);
182 outtrack->SetPt(1/pt_1);
183 outtrack->SetPsi(psi);
184 outtrack->SetTgl(par[3]);
185 outtrack->SetCharge(charge);
186 outtrack->CalculateHelix();
187 outtrack->Init(0,-1);
188
189 //cout<<"Loading track with "<<nhits<<" hits"<<endl;
190 for(int j=nhits-1; j>=0; j--)
191 {
192
193 Int_t index = track->GetClusterIndex(j);
194 if(index == 0)
195 continue;
196
197 Float_t xyz[3];
198 Int_t clustercharge =0;
199
200 //AliTPCcluster *cluster = (AliTPCcluster*)tracker->GetCluster(index);
201 AliCluster *cluster=0;
202
203 if(fMarian==kFALSE)
204 cluster = fTracker->GetCluster(index);
205 else
206 {
207 AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
208 cluster = mitracker->GetClusterMI(index);
209 }
210
211 xyz[1] = cluster->GetY();
212 xyz[2] = cluster->GetZ();
213 if(fMarian==kFALSE)
214 {
215 AliTPCcluster *cl = (AliTPCcluster*)cluster;
216 clustercharge = (Int_t)cl->GetQ();
217 }
218 else
219 {
220 AliTPCclusterMI *cl = (AliTPCclusterMI*)cluster;
221 clustercharge = (Int_t)cl->GetQ();
222 }
223
224 cluster->Use();
225
226 sec=(index&0xff000000)>>24; row=(index&0x00ff0000)>>16;
227
228 if(fMarian==kFALSE)
229 {
230 if(sec >= 18)
231 sec += 18;
232
233 if(xyz[2] < 0)
234 sec += 18;
235 }
236
237 //cout<<"sector "<<sec<<" row "<<row<<endl;
238 if(!AliL3Transform::Sector2Slice(slice,padrow,sec,row))
239 exit(5);
240 xyz[0] = AliL3Transform::Row2X(padrow);
241
242 //cout<<"Hit in slice "<<slice<<" padrow "<<padrow<<" index "<<index<<" y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<endl;
243 AliL3Transform::Local2Raw(xyz,sec,row);
244 //cout<<"slice "<<slice<<" padrow "<<padrow<<" pad "<<xyz[1]<<" time "<<xyz[2]<<endl;
245
246 if(xyz[1] < -1 || xyz[1] > AliL3Transform::GetNPads(padrow) ||
247 xyz[2] < -1 || xyz[2] > AliL3Transform::GetNTimeBins())
248 {
249 cerr<<"AliL3DataCompressor::FillOfflineData : Wrong time "<<xyz[2]<<" in slice "
250 <<slice<<" padrow "<<padrow<<endl;
251 cout<<"sector "<<sec<<" row "<<row<<endl;
252 //cout<<"Hit in slice "<<slice<<" padrow "<<padrow<<" y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<endl;
253 cout<<"Track hit "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
254 exit(5);
255 }
256
257 Float_t angle = 0;
258 AliL3Transform::Local2GlobalAngle(&angle,slice);
259 if(!outtrack->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
260 {
261 cerr<<"AliL3DataCompressor::FillOfflineData : Error in crossing point calc on slice "
262 <<slice<<" row "<<padrow<<endl;
263 exit(5);
264 }
265 Float_t xyz_cross[3] = {outtrack->GetPointX(),outtrack->GetPointY(),outtrack->GetPointZ()};
266 AliL3Transform::Global2Raw(xyz_cross,sec,row);
267 /*
268 if(fabs(xyz_cross[1] - xyz[1]) > 10 ||
269 fabs(xyz_cross[2] - xyz[2]) > 10)
270 {
271 cout<<"AliL3DataCompressor::FillOfflineData : Wrong crossing slice "<<slice<<" padrow "
272 <<padrow<<" pad "<<xyz[1]<<" padhit "<<xyz_cross[1]<<" time "<<xyz[2]<<" timehit "<<xyz_cross[2]<<endl;
273 outtrack->Print();
274 exit(5);
275 }
276 */
277 //cout<<" crossing "<<xyz_cross[0]<<" "<<xyz_cross[1]<<" "<<xyz_cross[2]<<endl;
278 outtrack->SetPadHit(padrow,xyz_cross[1]);
279 outtrack->SetTimeHit(padrow,xyz_cross[2]);
280
281 if(fWriteClusterShape)
282 {
283 Float_t angle = outtrack->GetCrossingAngle(padrow,slice);
284 outtrack->SetCrossingAngleLUT(padrow,angle);
285 outtrack->CalculateClusterWidths(padrow,kTRUE);
286 Int_t patch = AliL3Transform::GetPatch(padrow);
287 Float_t sigmaY2 = cluster->GetSigmaY2() / pow(AliL3Transform::GetPadPitchWidth(patch),2);
288 Float_t sigmaZ2 = cluster->GetSigmaZ2() / pow(AliL3Transform::GetZWidth(),2);
289 outtrack->SetCluster(padrow,xyz[1],xyz[2],clustercharge,sigmaY2,sigmaZ2,3);
290 }
291 else
292 outtrack->SetCluster(padrow,xyz[1],xyz[2],clustercharge,0,0,3);
293 totcounter++;
294 outtrack->GetClusterModel(padrow)->fSlice = slice;
b2a02bce 295 fNusedClusters++;
0a86fbb7 296 }
297 }
298
299 cout<<"AliL3DataCompressor::FillOfflineData : Wrote "<<totcounter<<" clusters"<<endl;
300 //Write tracks to file
301 AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
302 comp->WriteFile(comptracks);
303 delete comp;
304 delete comptracks;
305
306}
307
308void AliL3OfflineDataCompressor::WriteRemaining(Bool_t select)
309{
310 //Write remaining clusters (not assigned to any tracks) to file
311
312 if(!fKeepRemaining)
313 return;
314
315 if(select)
316 SelectRemainingClusters();
317
318 Char_t filename[1024];
319
320 if(!fSinglePatch)
321 {
322 cerr<<"AliL3OfflineDataCompressor::WriteRemaining : You have to modify this function when not running singlepatch"<<endl;
323 return;
324 }
b2a02bce 325
326 ofstream idfile;
327 if(fWriteIdsToFile)
328 {
329 sprintf(filename,"%s/comp/remains_ids.txt",fPath);
330 idfile.open(filename);
331 }
332
0a86fbb7 333 cout<<"Writing remaining clusters "<<endl;
334 Int_t nrows = AliL3Transform::GetNRows(),sector,row,sec;
335 AliTPCtracker *tracker = (AliTPCtracker*)fTracker;
336 for(Int_t slice=0; slice<=35; slice++)
337 {
338 sprintf(filename,"%s/comp/remains_%d_%d_%d.raw",fPath,fEvent,slice,-1);
339 FILE *outfile = fopen(filename,"w");
340 if(!outfile)
341 {
342 cerr<<"AliL3OfflineDataCompressor::WriteRemaining : Cannot open file "<<filename<<endl;
343 exit(5);
344 }
345
346 for(Int_t padrow=0; padrow < nrows; padrow++)
347 {
348 AliL3Transform::Slice2Sector(slice,padrow,sector,row);
349 sec=sector;
350
351 if(fMarian == kFALSE)
352 {
353 if(slice >= 18)
354 sec -= 18;
355 if(sec >= 18)
356 sec -= 18;
357 }
358 //cout<<"Getting clusters in sector "<<sec<<" row "<<row<<endl;
b2a02bce 359 Int_t ncl = 0;
360#ifdef asvversion
361 tracker->GetNClusters(sec,row);
362#endif
0a86fbb7 363
364 Int_t counter=0;
365 Int_t j;
366 for(j=0; j<ncl; j++)
367 {
b2a02bce 368 AliTPCcluster *cluster = 0;
369#ifdef asvversion
370 cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
371#endif
0a86fbb7 372 if(cluster->GetZ() < 0 && slice < 18) continue;
373 if(cluster->GetZ() > 0 && slice > 17) continue;
374 if(cluster->IsUsed())
375 continue;
376 counter++;
377 }
378
379 Int_t size = sizeof(AliL3RemainingRow) + counter*sizeof(AliL3RemainingCluster);
380 Byte_t *data = new Byte_t[size];
381 AliL3RemainingRow *tempPt = (AliL3RemainingRow*)data;
382
383 tempPt->fPadRow = padrow;
384 tempPt->fNClusters = counter;
385 //cout<<"Found "<<counter<<" unused out of "<<ncl<<" clusters on slice "<<slice<<" padrow "<<padrow<<endl;
386 Int_t local_counter=0;
387 for(j=0; j<ncl; j++)
388 {
b2a02bce 389 AliTPCcluster *cluster = 0;
390#ifdef asvversion
391 cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
392#endif
0a86fbb7 393 if(cluster->GetZ() < 0 && slice < 18) continue;
394 if(cluster->GetZ() > 0 && slice > 17) continue;
395 if(cluster->IsUsed())
396 continue;
b2a02bce 397
398 if(fWriteIdsToFile)
399 idfile << cluster->GetLabel(0)<<' ';
400
0a86fbb7 401 if(local_counter > counter)
402 {
403 cerr<<"AliL3OfflineDataCompressor::WriterRemaining : array out of range "<<local_counter<<" "<<counter<<endl;
404 return;
405 }
406 Float_t xyz[3] = {AliL3Transform::Row2X(padrow),cluster->GetY(),cluster->GetZ()};
407 AliL3Transform::Local2Raw(xyz,sector,row);
408 //cout<<"y "<<cluster->GetY()<<" z "<<cluster->GetZ()<<" pad "<<xyz[1]<<" time "<<xyz[2]<<endl;
409 tempPt->fClusters[local_counter].fY = cluster->GetY();
410 tempPt->fClusters[local_counter].fZ = cluster->GetZ();
411 tempPt->fClusters[local_counter].fCharge = (UShort_t)cluster->GetQ();
412 tempPt->fClusters[local_counter].fSigmaY2 = cluster->GetSigmaY2();
413 tempPt->fClusters[local_counter].fSigmaZ2 = cluster->GetSigmaZ2();
414 local_counter++;
b2a02bce 415 fNunusedClusters++;
0a86fbb7 416 }
417
418 fwrite(tempPt,size,1,outfile);
419 delete [] data;
420 }
421 fclose(outfile);
422 }
b2a02bce 423 if(fWriteIdsToFile)
424 {
425 idfile << endl;
426 idfile.close();
427 }
0a86fbb7 428}
429
430void AliL3OfflineDataCompressor::SelectRemainingClusters()
431{
432
433 cout<<"Cleaning up clusters"<<endl;
434 Int_t nrows = AliL3Transform::GetNRows();
435 Int_t gap=(Int_t)(0.125*nrows), shift=(Int_t)(0.5*gap);
436
437 Int_t sector,row,sec;
438 AliTPCtracker *tracker = (AliTPCtracker*)fTracker;
439 for(Int_t slice=0; slice<36; slice++)
440 {
441 for(Int_t padrow=0; padrow < nrows; padrow++)
442 {
443 if(padrow >= nrows-1-gap-shift) continue;//save all the clusters in this region
444
445 AliL3Transform::Slice2Sector(slice,padrow,sector,row);
446 sec=sector;
447
448 if(fMarian == kFALSE)
449 {
450 if(slice >= 18)
451 sec -= 18;
452 if(sec >= 18)
453 sec -= 18;
454 }
b2a02bce 455 Int_t ncl = 0;
456#ifdef asvversion
457 ncl=tracker->GetNClusters(sec,row);
458#endif
0a86fbb7 459 for(Int_t j=0; j<ncl; j++)
460 {
b2a02bce 461 AliTPCcluster *cluster = 0; //todo consti (AliTPCcluster*)tracker->GetCluster(sec,row,j);
0a86fbb7 462 if(cluster->GetZ() < 0 && slice < 18) continue;
463 if(cluster->GetZ() > 0 && slice > 17) continue;
464 if(cluster->IsUsed())
465 continue;
466
467 cluster->Use();
468 }
469 }
470
471 }
472}