Major bugfix in AliL3Compress::WriteRestoredData(). Also added naming conventions...
[u/mrichter/AliRoot.git] / HLT / comp / AliL3Modeller.cxx
1 //$Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ASV
5
6 #include <stream.h>
7 #include <iostream.h>
8 #include <math.h>
9
10 #include "AliL3Modeller.h"
11 #include "AliL3MemHandler.h"
12 #include "AliL3TrackArray.h"
13 #include "AliL3ModelTrack.h"
14 #include "AliL3DigitData.h"
15 #include "AliL3Transform.h"
16
17 #include "AliL3Defs.h"
18
19 ClassImp(AliL3Modeller)
20
21 AliL3Modeller::AliL3Modeller()
22 {
23   fMemHandler=0;
24   fTracks=0;
25   fTransform=0;
26 }
27
28
29 AliL3Modeller::~AliL3Modeller()
30 {
31   if(fMemHandler)
32     delete fMemHandler;
33   if(fTracks)
34     delete fTracks;
35   if(fTransform)
36     delete fTransform;
37
38 }
39
40 void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *path)
41 {
42   fSlice = slice;
43   fPatch = patch;
44   fPadOverlap=4;
45   fTimeOverlap=4;
46   fTrackThreshold=10;
47   sprintf(fPath,"%s",path);
48   
49   fTransform = new AliL3Transform();
50   fTracks = new AliL3TrackArray("AliL3ModelTrack");
51
52   AliL3MemHandler *file = new AliL3MemHandler();
53   if(!file->SetBinaryInput("tracks.raw"))
54     {
55       cerr<<"AliL3Modeller::Init : Error opening trackfile"<<endl;
56       return;
57     }
58   file->Binary2TrackArray(fTracks);
59   file->CloseBinaryInput();
60   delete file;
61   
62   fTracks->QSort();
63   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
64     {
65       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
66       if(!track) continue;
67       track->Init(fSlice,fPatch);
68       track->Rotate(fSlice,kTRUE);
69       track->CalculateHelix();
70     }    
71   
72   CalculateCrossingPoints();
73   CheckForOverlaps();
74
75   fMemHandler = new AliL3MemHandler();
76   Char_t fname[100];
77   sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
78   if(!fMemHandler->SetBinaryInput(fname))
79     {
80       cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
81       return;
82     }
83   UInt_t ndigits;
84   AliL3DigitRowData *digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
85   
86   SetInputData(digits);
87 }
88
89 void AliL3Modeller::FindClusters()
90 {
91   
92   if(!fTracks)
93     {
94       cerr<<"AliL3Modeller::Process : No tracks"<<endl;
95       return;
96     }
97   if(!fRowData)
98     {
99       cerr<<"AliL3Modeller::Process : No data "<<endl;
100       return;
101     }
102   
103   AliL3DigitRowData *rowPt = fRowData;
104   AliL3DigitData *digPt=0;
105
106   Int_t ntimes = fTransform->GetNTimeBins()+1;
107   Int_t npads = fTransform->GetNPads(NRows[fPatch][1])+1;//Max num of pads.
108   Int_t bounds = ntimes*npads;
109   Digit *row = new Digit[bounds];
110   
111   Int_t seq_charge;
112   Int_t pad,time,index;
113   Short_t charge;
114   Cluster cluster;
115
116   for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
117     {
118       memset((void*)row,0,ntimes*npads*sizeof(Digit));
119       digPt = (AliL3DigitData*)rowPt->fDigitData;
120       for(UInt_t j=0; j<rowPt->fNDigit; j++)
121         {
122           pad = digPt[j].fPad;
123           time = digPt[j].fTime;
124           charge = digPt[j].fCharge;
125           row[ntimes*pad+time].fCharge = charge;
126           row[ntimes*pad+time].fUsed = kFALSE;
127           //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
128         }
129       
130       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
131         {
132           AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
133           if(!track) continue;
134           //if(track->GetOverlap(i)>=0) continue;//Track is overlapping
135           
136           if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetOverlap(i)>=0)
137             {
138               track->SetCluster(i,0,0,0,0,0); //The track has left the patch.
139               continue;
140             }
141           
142           Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
143           Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
144           //cout<<"Checking track on row "<<i<<" with pad "<<hitpad<<" time "<<hittime<<endl;
145           pad = hitpad;
146           time = hittime;
147           Int_t padsign=-1;
148           Int_t timesign=-1;
149           
150           memset(&cluster,0,sizeof(Cluster));
151           
152           while(1)//Process this padrow
153             {
154               if(pad < 0 || pad >= fTransform->GetNPads(i)) 
155                 {
156                   //cout<<"Pad<0 on row "<<i<<endl;
157                   FillCluster(track,&cluster,i);
158                   break;
159                 }
160               seq_charge=0;
161               timesign=-1;
162               time = hittime;
163               while(1) //Process sequence on this pad:
164                 {
165                   if(time < 0) break;
166                   index = ntimes*pad + time;
167                   if(index < 0 || index >= bounds)
168                     {
169                       cerr<<"AliL3Modeller::FindClusters : Index out of range : "<<index
170                         <<" on row "<<i<<" pad "<<pad<<" time "<<time<<endl;
171                       break;
172                     }
173                   charge = row[index].fCharge;
174                   if(charge==0 && timesign==-1)
175                     {
176                       if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
177                         {
178                           time--;
179                           continue;
180                         }
181                       else
182                         {
183                           time = hittime+1;
184                           timesign=1;
185                           continue;
186                         }
187                     }
188                   else if(charge==0 && timesign==1)
189                     {
190                       if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap/2)
191                         {
192                           time++;
193                           continue;
194                         }
195                       else
196                         break;
197                     }
198                   
199                   //cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
200                   
201                   seq_charge += charge;
202                                                   
203                   cluster.fTime += time*charge;
204                   cluster.fPad += pad*charge;
205                   cluster.fCharge += charge;
206                   cluster.fSigmaY2 += pad*pad*charge;
207                   cluster.fSigmaZ2 += time*time*charge;
208                   
209                   row[ntimes*pad+time].fUsed = kTRUE;
210                   time += timesign;
211                 }
212               //cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
213               if(seq_charge)
214                 pad += padsign;
215               else //Nothing more on this pad, goto next pad
216                 {
217                   if(padsign==-1) 
218                     {
219                       if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2)
220                         {
221                           pad--; //In this case, we haven't found anything yet, 
222                         }        //so we will try to expand our search within the natural boundaries.
223                       else
224                         {
225                           pad=hitpad+1; 
226                           padsign=1; 
227                         }
228                       continue;
229                     }
230                   
231                   else if(padsign==1)
232                     {
233                       if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap/2)
234                         {
235                           pad++;     //In this case, we haven't found anything yet, 
236                           continue;  //so we will try to expand our search within the natural boundaries.
237                         }
238                       else //We are out of range, or cluster if finished.
239                         {
240                           //cout<<"Outof range; charge "<<cluster.fCharge<<" paddiff "<<abs(pad-hitpad)<<endl;
241                           FillCluster(track,&cluster,i);
242                           break;
243                         }
244                     }
245                   else //Nothing more in this cluster
246                     {
247                       //cout<<"Filling final cluster"<<endl;
248                       FillCluster(track,&cluster,i);
249                       break;
250                     } 
251                 }
252             }
253           //cout<<"done"<<endl;
254         }
255       
256       FillZeros(rowPt,row);
257       fMemHandler->UpdateRowPointer(rowPt);
258     }
259   delete [] row;
260   cout<<"done processing"<<endl;
261   
262   
263   //Debug:
264   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
265     {
266       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
267       if(!track) continue;
268       if(track->GetNClusters() != NumRows[fPatch])
269         cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<NumRows[fPatch]<<endl<<endl;
270     }
271   
272 }
273
274 void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row)
275 {
276   if(cluster->fCharge==0)
277     {
278       track->SetCluster(row,0,0,0,0,0);
279       return;
280     }
281   Float_t fcharge = (Float_t)cluster->fCharge;
282   Float_t fpad = ((Float_t)cluster->fPad/fcharge);
283   Float_t ftime = ((Float_t)cluster->fTime/fcharge);
284   Float_t sigmaY2,sigmaZ2;
285   CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
286   track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2);
287   
288 }
289
290 void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Digit *row)
291 {
292   //Fill zero where data has been used.
293   
294   Int_t ntimes = fTransform->GetNTimeBins()+1;
295   AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
296   for(UInt_t j=0; j<rowPt->fNDigit; j++)
297     {
298       Int_t pad = digPt[j].fPad;
299       Int_t time = digPt[j].fTime;
300       if(row[ntimes*pad+time].fUsed==kTRUE)
301         digPt[j].fCharge = 0;
302     }
303 }
304
305 void AliL3Modeller::WriteRemaining()
306 {
307   //Write remaining (nonzero) digits to file.
308   
309   AliL3DigitRowData *rowPt;
310   rowPt = (AliL3DigitRowData*)fRowData;
311   Int_t digitcount=0;
312   Int_t ndigits[(NumRows[fPatch])];
313   for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
314     {
315       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
316       ndigits[(i-NRows[fPatch][0])]=0;
317       for(UInt_t j=0; j<rowPt->fNDigit; j++)
318         {
319           if(digPt[j].fCharge==0) continue;
320           digitcount++;
321           ndigits[(i-NRows[fPatch][0])]++;
322         }
323       //cout<<"Difference "<<(int)ndigits[(i-NRows[fPatch][0])]<<" "<<(int)rowPt->fNDigit<<endl;
324       fMemHandler->UpdateRowPointer(rowPt);
325     }
326   
327   Int_t size = digitcount*sizeof(AliL3DigitData) + NumRows[fPatch]*sizeof(AliL3DigitRowData);
328   Byte_t *data = new Byte_t[size];
329   memset(data,0,size);
330   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)data;
331   rowPt = (AliL3DigitRowData*)fRowData;
332   
333   for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
334     {
335       Int_t localcount=0;
336       tempPt->fRow = i;
337       tempPt->fNDigit = ndigits[(i-NRows[fPatch][0])];
338       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
339       for(UInt_t j=0; j<rowPt->fNDigit; j++)
340         {
341           if(digPt[j].fCharge==0) continue;
342           if(localcount >= ndigits[(i-NRows[fPatch][0])])
343             {
344               cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
345               return;
346             }
347           tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
348           tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
349           tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
350           localcount++;
351         }
352       if(ndigits[(i-NRows[fPatch][0])] != localcount)
353         {
354           cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
355           return;
356         }
357       fMemHandler->UpdateRowPointer(rowPt);
358       Byte_t *tmp = (Byte_t*)tempPt;
359       Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-NRows[fPatch][0])]*sizeof(AliL3DigitData);
360       tmp += size;
361       tempPt = (AliL3DigitRowData*)tmp;
362     }
363
364   Char_t fname[100];
365   AliL3MemHandler *mem = new AliL3MemHandler();
366   sprintf(fname,"%s/remains_%d_%d.raw",fPath,fSlice,fPatch);
367   mem->SetBinaryOutput(fname);
368   mem->Memory2CompBinary((UInt_t)NumRows[fPatch],(AliL3DigitRowData*)data);
369   mem->CloseBinaryOutput();
370   delete mem;
371   delete [] data;
372 }
373
374
375 void AliL3Modeller::CalculateCrossingPoints()
376 {
377   cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
378   if(!fTracks)
379     {
380       cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
381       return;
382     }
383   Float_t hit[3];
384   
385   
386   //Remove tracks which are no good:
387   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
388     {
389       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
390       if(!track) continue;
391       if(track->GetFirstPointX() > fTransform->Row2X(NRows[fPatch][1]) || track->GetPt()<0.1)
392         fTracks->Remove(i);
393     }
394   
395   Int_t sector,row;
396   for(Int_t i=NRows[fPatch][1]; i>=NRows[fPatch][0]; i--)
397     {
398       for(Int_t j=0; j<fTracks->GetNTracks(); j++)
399         {
400           AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
401           if(!track) continue;
402                   
403           if(track->GetNHits() < fTrackThreshold)
404             {
405               fTracks->Remove(j);
406               continue;
407             }
408           
409           if(!track->GetCrossingPoint(i,hit)) 
410             {
411               cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
412                 "First point "<<track->GetFirstPointX()<<
413                 " nhits "<<track->GetNHits()<<endl;//" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
414                 //"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
415                 //"--------"<<endl;
416               fTracks->Remove(j);
417               continue;
418             }
419           //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
420           
421           fTransform->Slice2Sector(fSlice,i,sector,row);
422           fTransform->Local2Raw(hit,sector,row);
423           if(hit[1]<0 || hit[1]>fTransform->GetNPads(i) ||
424              hit[2]<0 || hit[2]>fTransform->GetNTimeBins())
425             {//Track is leaving the patch, so flag the track hits (<0)
426               track->SetPadHit(i,-1);
427               track->SetTimeHit(i,-1);
428               continue;
429             }
430             
431           track->SetPadHit(i,hit[1]);
432           track->SetTimeHit(i,hit[2]);
433           
434           //if(hit[1]<0 || hit[2]>445)
435           //if(hit[2]<0 || hit[2]>445)
436           //cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
437           //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
438         }
439     }
440   fTracks->Compress();
441   cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
442 }
443
444 void AliL3Modeller::CheckForOverlaps()
445 {
446   //Flag the tracks that overlap
447   
448   cout<<"Checking for overlaps...";
449   
450   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
451     {
452       AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
453       if(!track1) continue;
454       for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
455         {
456           AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
457           if(!track2) continue;
458           for(Int_t k=NRows[fPatch][0]; k<=NRows[fPatch][1]; k++)
459             {
460               if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
461                  track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
462                 continue;
463               if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) <= fPadOverlap &&
464                  fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) <= fTimeOverlap)
465                 {
466                   track2->SetOverlap(k,i);
467                   track1->SetOverlap(k,j);
468                 }
469             }
470         }
471     }
472   cout<<"done"<<endl;
473 }
474
475
476 void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
477 {
478   
479   Float_t padw,timew;
480   if(fPatch < 3)
481     padw = fTransform->GetPadPitchWidthLow();
482   else
483     padw = fTransform->GetPadPitchWidthUp();
484   Float_t charge = (Float_t)cl->fCharge;
485   Float_t pad = (Float_t)cl->fPad/charge;
486   Float_t time = (Float_t)cl->fTime/charge;
487   Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
488   sigmaY2 = (s2 + 1./12)*padw*padw;
489
490   if(s2 != 0)
491     {
492       sigmaY2 = sigmaY2*0.108;
493       if(fPatch<3)
494         sigmaY2 = sigmaY2*2.07;
495     }
496   
497   s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
498   timew = fTransform->GetZWidth();
499   sigmaZ2 = (s2 +1./12)*timew*timew;
500   if(s2 != 0)
501     {
502       sigmaZ2 = sigmaZ2*0.169;
503       if(fPatch < 3)
504         sigmaZ2 = sigmaZ2*1.77;
505     }
506   
507 }