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