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