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