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