Added support for NEWIO, merged cern-hlt tree, updated to latest track candidate...
[u/mrichter/AliRoot.git] / HLT / comp / AliL3OfflineDataCompressor.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
5
6 #include "AliL3StandardIncludes.h"
7
8 #include <AliTPCParamSR.h>
9 #include <AliTPCClustersArray.h>
10 #include <AliTPCcluster.h>
11 #include <AliTPCClustersRow.h>
12 #include <AliTPC.h>
13 #include <AliTPCv2.h>
14 #include <AliTPCcluster.h>
15 #include <AliTPCtracker.h>
16 #include <AliTPCclusterMI.h>
17 #include <AliTPCtrackerMI.h>
18 #include <AliKalmanTrack.h>
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
28 #include "AliL3Transform.h"
29 #include "AliL3ModelTrack.h"
30 #include "AliL3Compress.h"
31 #include "AliL3TrackArray.h"
32
33 #include "AliL3OfflineDataCompressor.h"
34
35 #if GCCVERSION == 3
36 using namespace std;
37 #endif
38
39 //_____________________________________________________________
40 //
41 //  AliL3OfflineDataCompression
42 //
43
44
45 ClassImp(AliL3OfflineDataCompressor)
46
47 AliL3OfflineDataCompressor::AliL3OfflineDataCompressor()
48 {
49   fMarian = kFALSE;
50   fTracker=0;
51 }
52
53 AliL3OfflineDataCompressor::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
60 AliL3OfflineDataCompressor::~AliL3OfflineDataCompressor()
61 {
62   if(fTracker)
63     {
64       fTracker->UnloadClusters();
65       delete fTracker;
66     }
67 }
68
69
70 void 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);
92 #ifdef asvversion
93   fTracker->LoadClusters();
94 #endif  
95   if(fMarian==kTRUE)
96     {
97       AliTPCtrackerMI *mitracker = (AliTPCtrackerMI*)fTracker;
98 #ifdef asvversion
99       mitracker->LoadInnerSectors();
100       mitracker->LoadOuterSectors();
101 #endif
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;
295           fNusedClusters++;
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
308 void 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     }
325   
326   ofstream idfile;
327   if(fWriteIdsToFile)
328     {
329       sprintf(filename,"%s/comp/remains_ids.txt",fPath);
330       idfile.open(filename);
331     }
332   
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;
359           Int_t ncl = 0;
360 #ifdef asvversion
361           tracker->GetNClusters(sec,row);
362 #endif
363                   
364           Int_t counter=0;
365           Int_t j;
366           for(j=0; j<ncl; j++)
367             {
368               AliTPCcluster *cluster = 0;
369 #ifdef asvversion
370               cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
371 #endif
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             {
389               AliTPCcluster *cluster = 0;
390 #ifdef asvversion
391               cluster=(AliTPCcluster*)tracker->GetCluster(sec,row,j);
392 #endif
393               if(cluster->GetZ() < 0 && slice < 18) continue;
394               if(cluster->GetZ() > 0 && slice > 17) continue;
395               if(cluster->IsUsed())
396                 continue;
397               
398               if(fWriteIdsToFile)
399                 idfile << cluster->GetLabel(0)<<' ';
400               
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++;
415               fNunusedClusters++;
416             }
417           
418           fwrite(tempPt,size,1,outfile);
419           delete [] data;
420         }
421       fclose(outfile);
422     }
423   if(fWriteIdsToFile)
424     {
425       idfile << endl;
426       idfile.close();
427     }
428 }
429
430 void 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             }
455           Int_t ncl = 0;
456 #ifdef asvversion
457           ncl=tracker->GetNClusters(sec,row);
458 #endif
459           for(Int_t j=0; j<ncl; j++)
460             {
461               AliTPCcluster *cluster = 0; //todo consti (AliTPCcluster*)tracker->GetCluster(sec,row,j);
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 }