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