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