f9ea20eb49c18c2ab2b63af53de37c9e53293bb9
[u/mrichter/AliRoot.git] / HLT / comp / AliL3ClusterFitter.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
5
6 #include "AliL3StandardIncludes.h"
7
8 #include "AliL3Logging.h"
9 #include "AliL3ClusterFitter.h"
10 #include "AliL3FitUtilities.h"
11 #include "AliL3DigitData.h"
12 #include "AliL3ModelTrack.h"
13 #include "AliL3TrackArray.h"
14 #include "AliL3MemHandler.h"
15 #include "AliL3HoughTrack.h"
16 #include "AliL3SpacePointData.h"
17 #include "AliL3Compress.h"
18
19 #if GCCVERSION == 3
20 using namespace std;
21 #endif
22
23 //_____________________________________________________________
24 //
25 //  AliL3ClusterFitter
26 //
27
28
29 ClassImp(AliL3ClusterFitter)
30
31 Int_t AliL3ClusterFitter::fBadFitError=0;
32 Int_t AliL3ClusterFitter::fFitError=0;
33 Int_t AliL3ClusterFitter::fResultError=0;
34 Int_t AliL3ClusterFitter::fFitRangeError=0;
35
36 AliL3ClusterFitter::AliL3ClusterFitter()
37 {
38   plane=0;
39   fNmaxOverlaps = 3;
40   fChiSqMax=12;
41   fRowMin=-1;
42   fRowMax=-1;
43   fFitted=0;
44   fFailed=0;
45   fYInnerWidthFactor=1;
46   fZInnerWidthFactor=1;
47   fYOuterWidthFactor=1;
48   fZOuterWidthFactor=1;
49   fSeeds=0;
50   fProcessTracks=0;
51   fClusters=0;
52   fNMaxClusters=0;
53   fNClusters=0;
54 }
55
56 AliL3ClusterFitter::AliL3ClusterFitter(Char_t *path)
57 {
58   strcpy(fPath,path);
59   plane=0;
60   fNmaxOverlaps = 3;
61   fChiSqMax=12;
62   fRowMin=-1;
63   fRowMax=-1;
64   fFitted=0;
65   fFailed=0;
66   fYInnerWidthFactor=1;
67   fZInnerWidthFactor=1;
68   fYOuterWidthFactor=1;
69   fZOuterWidthFactor=1;
70   fSeeds=0;
71   fProcessTracks=0;
72   fNMaxClusters=100000;
73   fClusters=0;
74   fNClusters=0;
75 }
76
77 AliL3ClusterFitter::~AliL3ClusterFitter()
78 {
79   if(fSeeds)
80     delete fSeeds;
81   if(fClusters)
82     delete [] fClusters;
83 }
84
85 void AliL3ClusterFitter::Init(Int_t slice,Int_t patch,Int_t *rowrange,AliL3TrackArray *tracks)
86 {
87   //Assuming tracklets found by the line transform
88
89   fSlice=slice;
90   fPatch=patch;
91   
92   if(rowrange[0] > AliL3Transform::GetLastRow(patch) || rowrange[1] < AliL3Transform::GetFirstRow(patch))
93     cerr<<"AliL3ClusterFitter::Init : Wrong rows "<<rowrange[0]<<" "<<rowrange[1]<<endl;
94   fRowMin=rowrange[0];
95   fRowMax=rowrange[1];
96
97   if(fRowMin < 0)
98     fRowMin = 0;
99   if(fRowMax > AliL3Transform::GetLastRow(fPatch))
100     fRowMax = AliL3Transform::GetLastRow(fPatch);
101   
102   fFitted=fFailed=0;
103   
104   Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
105   Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
106   Int_t bounds = ntimes*npads;
107   if(fRow)
108     delete [] fRow;
109   fRow = new Digit[bounds];
110   if(fTracks)
111     delete fTracks;
112   
113   fTracks = new AliL3TrackArray("AliL3ModelTrack");
114   
115   for(Int_t i=0; i<tracks->GetNTracks(); i++)
116     {
117       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
118       if(!track) continue;
119       AliL3ModelTrack *mtrack = (AliL3ModelTrack*)fTracks->NextTrack();
120       mtrack->Init(slice,patch);
121       mtrack->SetTgl(track->GetTgl());
122       mtrack->SetRowRange(rowrange[0],rowrange[1]);
123       for(Int_t j=fRowMin; j<=fRowMax; j++)
124         {
125           Float_t hit[3];
126           track->GetLineCrossingPoint(j,hit);
127           hit[0] += AliL3Transform::Row2X(track->GetFirstRow());
128           Float_t R = sqrt(hit[0]*hit[0] + hit[1]*hit[1]);
129           hit[2] = R*track->GetTgl();
130           Int_t se,ro;
131           AliL3Transform::Slice2Sector(slice,j,se,ro);
132           AliL3Transform::Local2Raw(hit,se,ro);
133           if(hit[1]<0 || hit[1]>=AliL3Transform::GetNPads(j) || hit[2]<0 || hit[2]>=AliL3Transform::GetNTimeBins())
134             {
135               mtrack->SetPadHit(j,-1);
136               mtrack->SetTimeHit(j,-1);
137               continue;
138             }
139           mtrack->SetPadHit(j,hit[1]);
140           mtrack->SetTimeHit(j,hit[2]);
141           mtrack->SetCrossingAngleLUT(j,fabs(track->GetPsiLine() - AliL3Transform::Pi()/2));
142           //if(mtrack->GetCrossingAngleLUT(j) > AliL3Transform::Deg2Rad(20))
143           //  cout<<"Angle "<<mtrack->GetCrossingAngleLUT(j)<<" psiline "<<track->GetPsiLine()*180/3.1415<<endl;
144           mtrack->CalculateClusterWidths(j);
145         }
146     }
147   //  cout<<"Copied "<<fTracks->GetNTracks()<<" tracks "<<endl;
148 }
149
150 void AliL3ClusterFitter::Init(Int_t slice,Int_t patch)
151 {
152   fSlice=slice;
153   fPatch=patch;
154
155   fRowMin=AliL3Transform::GetFirstRow(patch);
156   fRowMax=AliL3Transform::GetLastRow(patch);
157   
158   fFitted=fFailed=0;
159   
160   Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
161   Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
162   Int_t bounds = ntimes*npads;
163   if(fRow)
164     delete [] fRow;
165   fRow = new Digit[bounds];
166   if(fTracks)
167     delete fTracks;
168   fTracks = new AliL3TrackArray("AliL3ModelTrack");  
169
170 }
171
172
173 void AliL3ClusterFitter::LoadSeeds(Int_t *rowrange,Bool_t offline)
174 {
175   
176   cout<<"Loading the seeds"<<endl;
177   Char_t fname[1024];
178   
179   if(offline)
180     sprintf(fname,"%s/offline/tracks_%d.raw",fPath,0);
181   else
182     sprintf(fname,"%s/hough/tracks_%d.raw",fPath,0);
183   
184   cout<<"AliL3ClusterFitter::LoadSeeds : Loading input tracks from "<<fname<<endl;
185   
186   AliL3MemHandler tfile;
187   tfile.SetBinaryInput(fname);
188   
189   if(fSeeds)
190     delete fSeeds;
191   fSeeds = new AliL3TrackArray("AliL3ModelTrack");
192   tfile.Binary2TrackArray(fSeeds);
193   tfile.CloseBinaryInput();
194
195   //if(!offline)
196   //fSeeds->QSort();
197   
198   Int_t clustercount=0;
199   for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
200     {
201       AliL3ModelTrack *track = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
202       if(!track) continue;
203
204       if(!offline)
205         {
206           /*
207           if(track->GetNHits() < 10 || track->GetPt() < 0.08) 
208             {
209               fSeeds->Remove(i);
210               continue;
211             }
212           */
213         }
214       clustercount += track->GetNHits();
215       track->CalculateHelix();
216       
217       Int_t nhits = track->GetNHits();
218       UInt_t *hitids = track->GetHitNumbers();
219
220       Int_t origslice = (hitids[nhits-1]>>25)&0x7f;//Slice of innermost point
221
222       track->Init(origslice,-1);
223       Int_t slice = origslice;
224       
225       //for(Int_t j=rowrange[1]; j>=rowrange[0]; j--)
226       for(Int_t j=rowrange[0]; j<=rowrange[1]; j++)
227         {
228           
229           //Calculate the crossing point between track and padrow
230           Float_t angle = 0; //Perpendicular to padrow in local coordinates
231           AliL3Transform::Local2GlobalAngle(&angle,slice);
232           if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
233             {
234               cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
235               continue;
236               //track->Print();
237               //exit(5);
238             }
239           Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
240           
241           Int_t sector,row;
242           AliL3Transform::Slice2Sector(slice,j,sector,row);
243           AliL3Transform::Global2Raw(xyz_cross,sector,row);
244           //cout<<"Examining slice "<<slice<<" row "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
245           if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //Track leaves the slice
246             {
247             newslice:
248               
249               Int_t tslice=slice;
250               Float_t lastcross=xyz_cross[1];
251               if(xyz_cross[1] > 0)
252                 {
253                   if(slice == 17)
254                     slice=0;
255                   else if(slice == 35)
256                     slice = 18;
257                   else
258                     slice += 1;
259                 }
260               else
261                 {
262                   if(slice == 0)
263                     slice = 17;
264                   else if(slice==18)
265                     slice = 35;
266                   else
267                     slice -= 1;
268                 }
269               if(slice < 0 || slice>35)
270                 {
271                   cerr<<"Wrong slice "<<slice<<" on row "<<j<<endl;
272                   exit(5);
273                 }
274               //cout<<"Track leaving, trying slice "<<slice<<endl;
275               angle=0;
276               AliL3Transform::Local2GlobalAngle(&angle,slice);
277               if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
278                 {
279                   cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
280                   continue;
281                   //track->Print();
282                   //exit(5);
283                 }
284               xyz_cross[0] = track->GetPointX();
285               xyz_cross[1] = track->GetPointY();
286               xyz_cross[2] = track->GetPointZ();
287               Int_t sector,row;
288               AliL3Transform::Slice2Sector(slice,j,sector,row);
289               AliL3Transform::Global2Raw(xyz_cross,sector,row);
290               if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //track is in the borderline
291                 {
292                   if(xyz_cross[1] > 0 && lastcross > 0 || xyz_cross[1] < 0 && lastcross < 0)
293                     goto newslice;
294                   else
295                     {
296                       slice = tslice;//Track is on the border of two slices
297                       continue;
298                     }
299                 }
300             }
301           
302           if(xyz_cross[2] < 0 || xyz_cross[2] >= AliL3Transform::GetNTimeBins())//track goes out of range
303             continue;
304           
305           if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j))
306             {
307               cerr<<"Slice "<<slice<<" padrow "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
308               track->Print();
309               exit(5);
310             }
311           
312           track->SetPadHit(j,xyz_cross[1]);
313           track->SetTimeHit(j,xyz_cross[2]);
314           angle=0;
315           AliL3Transform::Local2GlobalAngle(&angle,slice);
316           Float_t crossingangle = track->GetCrossingAngle(j,slice);
317           track->SetCrossingAngleLUT(j,crossingangle);
318           
319           track->CalculateClusterWidths(j);
320           
321           track->GetClusterModel(j)->fSlice = slice;
322           
323         }
324       memset(track->GetHitNumbers(),0,159*sizeof(UInt_t));
325       track->SetNHits(0);
326     }
327   fSeeds->Compress();
328   
329   AliL3Compress *c = new AliL3Compress(-1,-1,fPath);
330   c->WriteFile(fSeeds,"tracks_before.raw");
331   delete c;
332   
333   cout<<"Loaded "<<fSeeds->GetNTracks()<<" seeds and "<<clustercount<<" clusters"<<endl;
334 }
335
336 void AliL3ClusterFitter::FindClusters()
337 {
338   if(!fTracks)
339     {
340       cerr<<"AliL3ClusterFitter::Process : No tracks"<<endl;
341       return;
342     }
343   if(!fRowData)
344     {
345       cerr<<"AliL3ClusterFitter::Process : No data "<<endl;
346       return;
347     }
348   
349   AliL3DigitRowData *rowPt = fRowData;
350   AliL3DigitData *digPt=0;
351
352   Int_t pad,time;
353   Short_t charge;
354   
355   if(fRowMin < 0)
356     {
357       fRowMin = AliL3Transform::GetFirstRow(fPatch);
358       fRowMax = AliL3Transform::GetLastRow(fPatch);
359     }
360   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
361     {
362       if((Int_t)rowPt->fRow < fRowMin)
363         {
364           AliL3MemHandler::UpdateRowPointer(rowPt);
365           continue;
366         }
367       else if((Int_t)rowPt->fRow > fRowMax)
368         break;
369       else if((Int_t)rowPt->fRow != i)
370         {
371           cerr<<"AliL3ClusterFitter::FindClusters : Mismatching row numbering "<<i<<" "<<rowPt->fRow<<endl;
372           exit(5);
373         }
374       fCurrentPadRow = i;
375       memset((void*)fRow,0,(AliL3Transform::GetNTimeBins()+1)*(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
376       digPt = (AliL3DigitData*)rowPt->fDigitData;
377
378       for(UInt_t j=0; j<rowPt->fNDigit; j++)
379         {
380           pad = digPt[j].fPad;
381           time = digPt[j].fTime;
382           charge = digPt[j].fCharge;
383           if(charge > 1024)
384             charge -= 1024;
385           fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge = charge;
386           fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kFALSE;
387           //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
388         }
389       
390       for(Int_t it=0; it<2; it++)
391         {
392           if(it==0)
393             {
394               fProcessTracks = fSeeds;
395               fSeeding = kTRUE;
396             }
397           else
398             {
399               fProcessTracks = fTracks;
400               fSeeding = kFALSE;
401             }
402           if(!fProcessTracks)
403             continue;
404           
405           for(Int_t k=0; k<fProcessTracks->GetNTracks(); k++)
406             {
407               AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(k);
408               if(!track) continue;
409               
410               if(fSeeding)
411                 if(track->GetClusterModel(i)->fSlice != fSlice) continue;
412               
413               if(track->GetPadHit(i) < 0 || track->GetPadHit(i) > AliL3Transform::GetNPads(i)-1 ||
414                  track->GetTimeHit(i) < 0 || track->GetTimeHit(i) > AliL3Transform::GetNTimeBins()-1)
415                 {
416                   track->SetCluster(i,0,0,0,0,0,0);
417                   continue;
418                 }
419               
420               if(CheckCluster(k) == kFALSE)
421                 fFailed++;
422             }
423         }
424       //FillZeros(rowPt);
425       AliL3MemHandler::UpdateRowPointer(rowPt);
426     }
427   
428   fSeeding = kTRUE;
429   AddClusters();
430   fSeeding = kFALSE;
431   AddClusters();
432     
433   cout<<"Fitted "<<fFitted<<" clusters, failed "<<fFailed<<endl;
434   cout<<"Distribution:"<<endl;
435   cout<<"Bad fit "<<fBadFitError<<endl;
436   cout<<"Fit error "<<fFitError<<endl;
437   cout<<"Result error "<<fResultError<<endl;
438   cout<<"Fit range error "<<fFitRangeError<<endl;
439
440 }
441
442 Bool_t AliL3ClusterFitter::CheckCluster(Int_t trackindex)
443 {
444   //Check if this is a single or overlapping cluster
445   
446   AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(trackindex);
447   
448   Int_t row = fCurrentPadRow;
449   
450   if(track->IsSet(row)) //A fit has already be performed on this one
451     return kTRUE;
452   
453   //Define the cluster region of this hit:
454   Int_t padr[2]={999,-1};
455   Int_t timer[2]={999,-1};
456   
457   if(!SetFitRange(track,padr,timer))
458     {
459       track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
460       
461       if(fDebug)
462         cout<<"Failed to fit cluster at row "<<row<<" pad "<<(Int_t)rint(track->GetPadHit(row))<<" time "
463             <<(Int_t)rint(track->GetTimeHit(row))<<" hitcharge "
464             <<fRow[(AliL3Transform::GetNTimeBins()+1)*(Int_t)rint(track->GetPadHit(row))+(Int_t)rint(track->GetTimeHit(row))].fCharge<<endl;
465       fFitRangeError++;
466       return kFALSE;
467     }
468
469   //Check if any other track contributes to this cluster:
470   //This is done by checking if the tracks are overlapping within
471   //the range defined by the track parameters
472   for(Int_t t=trackindex+1; t<fProcessTracks->GetNTracks(); t++)
473     {
474       AliL3ModelTrack *tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(t);
475       if(!tr) continue;
476       if(fSeeding)
477         if(tr->GetClusterModel(row)->fSlice != fSlice) continue;
478       Int_t xyw = (Int_t)ceil(sqrt(tr->GetParSigmaY2(row))) + 1;
479       Int_t zw = (Int_t)ceil(sqrt(tr->GetParSigmaZ2(row))) + 1;
480       if( 
481          (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
482           tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
483          
484          (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
485           tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
486          
487          (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
488           tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1]) ||
489          
490          (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
491           tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1]) 
492          )
493         {
494           if(SetFitRange(tr,padr,timer)) //Expand the cluster fit range
495             track->SetOverlap(row,t);    //Set overlap
496         }
497     }
498
499   if(fDebug)
500     cout<<"Fitting cluster with "<<track->GetNOverlaps(fCurrentPadRow)<<" overlaps"<<endl;
501   FitClusters(track,padr,timer);
502   return kTRUE;
503 }
504
505 Bool_t AliL3ClusterFitter::SetFitRange(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
506 {
507   Int_t row = fCurrentPadRow;
508   Int_t nt = AliL3Transform::GetNTimeBins()+1;
509   
510   Int_t nsearchbins=0;
511   if(row < 63)
512     nsearchbins=25;
513   else
514     nsearchbins=49;
515   Int_t padloop[49] = {0,0,0,-1,1,-1,1,-1,1,0,0,-1,1,-1,1 ,2,-2,2,-2,2,-2,2,-2,2,-2
516                        ,0,1,2,3,3,3,3,3,3,3
517                        ,2,1,0,-1,-2,-3
518                        ,-3,-3,-3,-3,-3,-3,-1,-1};
519   Int_t timeloop[49] = {0,1,-1,0,0,1,1,-1,-1,2,-2,2,2,-2,-2 ,0,0,1,1,-1,-1,2,2,-2,-2
520                         ,-3,-3,-3,-3,-2,-1,0,1,2,3
521                         ,3,3,3,3,3,3,2,1,0,-1,-2,-3,-3,-3};
522   
523   Int_t padhit = (Int_t)rint(track->GetPadHit(row));
524   Int_t timehit = (Int_t)rint(track->GetTimeHit(row));
525   Int_t padmax=-1;
526   Int_t timemax=-1;
527
528   for(Int_t index=0; index<nsearchbins; index++)
529     {
530       if(IsMaximum(padhit + padloop[index],timehit + timeloop[index])) 
531         {
532           padmax = padhit + padloop[index];
533           timemax = timehit + timeloop[index];
534           break;
535         }
536     }
537
538   /*
539   if(IsMaximum(padhit,timehit))
540     {
541       padmax = padhit;
542       timemax = timehit;
543     }
544   */
545   //Define the cluster region of this hit:
546   //The region we look for, is centered at the local maxima
547   //and expanded around using the parametrized cluster width
548   //according to track parameters.
549   Int_t xyw = (Int_t)ceil(sqrt(track->GetParSigmaY2(row)))+1; 
550   Int_t zw = (Int_t)ceil(sqrt(track->GetParSigmaZ2(row)))+1;  
551   if(padmax>=0 && timemax>=0)
552     {
553       if(fDebug)
554         {
555           cout<<"Expanding cluster range using expected cluster widths: "<<xyw<<" "<<zw
556               <<" and setting local maxima pad "<<padmax<<" time "<<timemax<<endl;
557           if(xyw > 10 || zw > 10)
558             track->Print();
559         }
560       
561       //Set the hit to the local maxima of the cluster.
562       //Store the maxima in the cluster model structure,
563       //-only temporary, it will be overwritten when calling SetCluster.
564       
565       track->GetClusterModel(row)->fDPad = padmax;
566       track->GetClusterModel(row)->fDTime = timemax;
567
568       for(Int_t i=padmax-xyw; i<=padmax+xyw; i++)
569         {
570           for(Int_t j=timemax-zw; j<=timemax+zw; j++)
571             {
572               if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) continue;
573               if(fRow[nt*i+j].fCharge)
574                 {
575                   if(i < padrange[0]) padrange[0]=i;
576                   if(i > padrange[1]) padrange[1]=i;
577                   if(j < timerange[0]) timerange[0]=j;
578                   if(j > timerange[1]) timerange[1]=j;
579                 }
580             }
581         }
582       if(fDebug)
583         cout<<"New padrange "<<padrange[0]<<" "<<padrange[1]<<" "<<" time "<<timerange[0]<<" "<<timerange[1]<<endl;
584       return kTRUE;
585     }
586   return kFALSE;
587 }
588
589 Bool_t AliL3ClusterFitter::IsMaximum(Int_t pad,Int_t time)
590 {
591   if(pad<0 || pad >= AliL3Transform::GetNPads(fCurrentPadRow) ||
592      time<0 || time >= AliL3Transform::GetNTimeBins())
593     return kFALSE;
594   Int_t nt = AliL3Transform::GetNTimeBins()+1;
595   if(fRow[nt*pad+time].fUsed == kTRUE) return kFALSE; //Peak has been assigned before
596   Int_t charge = fRow[nt*pad+time].fCharge;
597   if(charge == 1023 || charge==0) return kFALSE;
598   
599   fRow[nt*pad+time].fUsed = kTRUE;
600   return kTRUE;
601
602   //if(charge < fRow[nt*(pad-1)+(time-1)].fCharge) return kFALSE;
603   if(charge < fRow[nt*(pad)+(time-1)].fCharge) return kFALSE;
604   //if(charge < fRow[nt*(pad+1)+(time-1)].fCharge) return kFALSE;
605   if(charge < fRow[nt*(pad-1)+(time)].fCharge) return kFALSE;
606   if(charge < fRow[nt*(pad+1)+(time)].fCharge) return kFALSE;
607   //if(charge < fRow[nt*(pad-1)+(time+1)].fCharge) return kFALSE;
608   if(charge < fRow[nt*(pad)+(time+1)].fCharge) return kFALSE;
609   //if(charge < fRow[nt*(pad+1)+(time+1)].fCharge) return kFALSE;
610   fRow[nt*pad+time].fUsed = kTRUE;
611   return kTRUE;
612 }
613
614 void AliL3ClusterFitter::FitClusters(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
615 {
616   //Handle single and overlapping clusters
617     
618   //Check whether this cluster has been set before:
619   
620   Int_t size = FIT_PTS;
621   Int_t max_tracks = FIT_MAXPAR/NUM_PARS;
622   if(track->GetNOverlaps(fCurrentPadRow) > max_tracks)
623     {
624       cerr<<"AliL3ClusterFitter::FitOverlappingClusters : Too many overlapping tracks"<<endl;
625       return;
626     }
627   Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
628   
629   //Check if at least one cluster is not already fitted
630   Bool_t all_fitted=kTRUE;
631   
632   Int_t k=-1;
633   while(k < track->GetNOverlaps(fCurrentPadRow))
634     {
635       AliL3ModelTrack *tr=0;
636       if(k==-1)
637         tr = track;
638       else
639         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
640       k++;
641       if(!tr) continue;
642       if(!tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow))//cluster has not been set and is not present
643         {
644           all_fitted = kFALSE;
645           break;
646         }
647     }
648   if(all_fitted)
649     {
650       if(fDebug)
651         cout<<"But all the clusters were already fitted on row "<<fCurrentPadRow<<endl;
652       return;
653     }
654   
655   //Allocate fit parameters array; this is interface to the C code
656   plane = new DPOINT[FIT_PTS];
657   memset(plane,0,FIT_PTS*sizeof(DPOINT));
658
659   Double_t x[FIT_PTS],y[FIT_PTS],s[FIT_PTS];
660   
661   //Fill the fit parameters:
662   Double_t a[FIT_MAXPAR];
663   Int_t lista[FIT_MAXPAR];
664   Double_t dev[FIT_MAXPAR],chisq_f;
665   
666   Int_t fit_pars=0;
667   
668   Int_t n_overlaps=0;
669   k=-1;
670   
671   //Fill the overlapping tracks:
672   while(k < track->GetNOverlaps(fCurrentPadRow))
673     {
674       AliL3ModelTrack *tr=0;
675       if(k==-1)
676         tr = track;
677       else
678         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
679       k++;
680       if(!tr) continue;
681       
682       if(tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow)) continue;//Cluster fit failed before
683       
684       //Use the local maxima as the input to the fitting routine.
685       //The local maxima is temporary stored in the cluster model:
686       Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad);  //rint(tr->GetPadHit(fCurrentPadRow));
687       Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime); //rint(tr->GetTimeHit(fCurrentPadRow));
688       Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fCharge;
689       
690       if(fDebug)
691         cout<<"Fitting track cluster, pad "<<tr->GetPadHit(fCurrentPadRow)<<" time "
692             <<tr->GetTimeHit(fCurrentPadRow)<<" charge "<<charge<<" at local maxima in pad "<<hitpad
693             <<" time "<<hittime<<" xywidth "<<sqrt(tr->GetParSigmaY2(fCurrentPadRow))
694             <<" zwidth "<<sqrt(tr->GetParSigmaZ2(fCurrentPadRow))<<endl;
695       
696       if(charge==0)
697         {
698           cerr<<"Charge still zero!"<<endl;
699           exit(5);
700         }
701             
702       a[n_overlaps*NUM_PARS+2] = hitpad;
703       a[n_overlaps*NUM_PARS+4] = hittime;
704       
705       if(!tr->IsSet(fCurrentPadRow)) //Cluster is not fitted before
706         {
707           a[n_overlaps*NUM_PARS+1] = charge;
708           a[n_overlaps*NUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow)) * GetYWidthFactor();
709           a[n_overlaps*NUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
710           a[n_overlaps*NUM_PARS+6] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
711           lista[n_overlaps*NUM_PARS + 1] = 1;
712           lista[n_overlaps*NUM_PARS + 2] = 1;
713           lista[n_overlaps*NUM_PARS + 3] = 0;
714           lista[n_overlaps*NUM_PARS + 4] = 1;
715           lista[n_overlaps*NUM_PARS + 5] = 0;
716           lista[n_overlaps*NUM_PARS + 6] = 0;
717           fit_pars             += 3;
718         }
719       else  //Cluster was fitted before
720         {
721           if(!tr->IsPresent(fCurrentPadRow))
722             {
723               cerr<<"AliL3ClusterFitter::FindClusters : Cluster not present; there is a bug here"<<endl;
724               exit(5);
725             }
726           Int_t charge;
727           Float_t xywidth,zwidth,pad,time;
728           tr->GetPad(fCurrentPadRow,pad);
729           tr->GetTime(fCurrentPadRow,time);
730           tr->GetClusterCharge(fCurrentPadRow,charge);
731           xywidth = sqrt(tr->GetParSigmaY2(fCurrentPadRow));
732           zwidth = sqrt(tr->GetParSigmaZ2(fCurrentPadRow));
733           if(fDebug)
734             cout<<"Cluster had been fitted before, pad "<<pad<<" time "<<time<<" charge "<<charge<<" width "<<xywidth<<" "<<zwidth<<endl;
735           
736           a[n_overlaps*NUM_PARS+2] = pad;
737           a[n_overlaps*NUM_PARS+4] = time;
738           a[n_overlaps*NUM_PARS+1] = charge;
739           a[n_overlaps*NUM_PARS+3] = sqrt(xywidth) * GetYWidthFactor();
740           a[n_overlaps*NUM_PARS+5] = sqrt(zwidth) * GetZWidthFactor();
741           a[n_overlaps*NUM_PARS+6] = sqrt(zwidth) * GetZWidthFactor();
742
743           lista[n_overlaps*NUM_PARS + 1] = 1;
744           lista[n_overlaps*NUM_PARS + 2] = 0;
745           lista[n_overlaps*NUM_PARS + 3] = 0;
746           lista[n_overlaps*NUM_PARS + 4] = 0;
747           lista[n_overlaps*NUM_PARS + 5] = 0;
748           lista[n_overlaps*NUM_PARS + 6] = 0;
749           fit_pars             += 1;
750         }
751       n_overlaps++;
752     }
753   
754   if(n_overlaps==0) //No clusters here
755     {
756       delete [] plane;
757       return;
758     }
759
760   Int_t pad_num=0;
761   Int_t time_num_max=0;
762   Int_t ndata=0;
763   Int_t tot_charge=0;
764   if(fDebug)
765     cout<<"Padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "<<timerange[0]<<" "<<timerange[1]<<endl;
766   for(Int_t i=padrange[0]; i<=padrange[1]; i++)
767     {
768       Int_t max_charge = 0;
769       Int_t time_num=0;
770       for(Int_t j=timerange[0]; j<=timerange[1]; j++)
771         {
772           Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
773           
774           if(charge <= 0) continue;
775
776           time_num++;
777           if(charge > max_charge)
778             {
779               max_charge = charge;
780               //time_num++;
781             }
782           if(fDebug)
783             cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
784           tot_charge += charge;
785           ndata++;
786           if(ndata >= size)
787             {
788               cerr<<"Too many points; row "<<fCurrentPadRow<<" padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "
789                   <<timerange[0]<<" "<<timerange[1]<<endl;
790               exit(5);
791             }
792
793           plane[ndata].u = (Double_t)i;
794           plane[ndata].v = (Double_t)j;
795           x[ndata]=ndata;
796           y[ndata]=charge;
797           s[ndata]= 1 + sqrt((Double_t)charge);
798         }
799       if(max_charge) //there was charge on this pad
800         pad_num++;
801       if(time_num_max < time_num)
802         time_num_max = time_num;
803     }
804   
805   if(pad_num <= 1 || time_num_max <=1 || n_overlaps > fNmaxOverlaps || ndata <= fit_pars) //too few to do fit
806     {
807       SetClusterfitFalse(track);
808       if(fDebug)
809         cout<<"Too few digits or too many overlaps: "<<pad_num<<" "<<time_num_max<<" "<<n_overlaps<<" ndata "<<ndata<<" fit_pars "<<fit_pars<<endl;
810       delete [] plane;
811       return;
812     }
813
814   
815   Int_t npars = n_overlaps * NUM_PARS;
816   if(fDebug)
817     cout<<"Number of overlapping clusters "<<n_overlaps<<endl;
818   Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisq_f, f2gauss5 );
819   
820   if(ret<0)
821     {
822       SetClusterfitFalse(track);
823       fFailed++;
824       fFitError++;
825       delete [] plane;
826       return;
827       //exit(5);
828     }
829
830   chisq_f /= (ndata-fit_pars);
831   if(fDebug)
832     cout<<"Chisq "<<chisq_f<<endl;
833   
834   k=-1;
835   n_overlaps=0;
836   while(k < track->GetNOverlaps(fCurrentPadRow))
837     {
838       AliL3ModelTrack *tr=0;
839       if(k==-1)
840         tr = track;
841       else
842         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
843       k++;
844       if(!tr) continue;
845       if(!tr->IsPresent(fCurrentPadRow))
846         {
847           if(tr->IsSet(fCurrentPadRow)) continue;//This cluster has been set before
848           
849           if(chisq_f < fChiSqMax)//cluster fit is good enough
850             {
851               tot_charge = (Int_t)(a[n_overlaps*NUM_PARS+1] * a[n_overlaps*NUM_PARS+3] * a[n_overlaps*NUM_PARS+5]);
852               Float_t fpad = a[n_overlaps*NUM_PARS+2];
853               Float_t ftime = a[n_overlaps*NUM_PARS+4];
854               if(tot_charge < 0 || fpad < -1 || fpad > AliL3Transform::GetNPads(fCurrentPadRow) || 
855                  ftime < -1 || ftime > AliL3Transform::GetNTimeBins())
856                 {
857                   if(fDebug)
858                     cout<<"AliL3ClusterFitter::Fatal result(s) in fit; in slice "<<fSlice<<" row "<<fCurrentPadRow
859                         <<"; pad "<<fpad<<" time "<<ftime<<" charge "<<tot_charge<<" xywidth "<<a[n_overlaps*NUM_PARS+3]
860                         <<" zwidth "<<a[n_overlaps*NUM_PARS+5]<<" peakcharge "<<a[n_overlaps*NUM_PARS+1]<<endl;
861                   tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
862                   fFailed++;
863                   fResultError++;
864                   continue;
865                 }
866               
867               tr->SetCluster(fCurrentPadRow,fpad,ftime,tot_charge,0,0,pad_num);
868               if(fDebug)
869                 cout<<"Setting cluster in pad "<<a[n_overlaps*NUM_PARS+2]<<" time "<<a[n_overlaps*NUM_PARS+4]<<" charge "<<tot_charge<<endl;
870               /*
871               //Set the digits to used:
872               for(Int_t i=padrange[0]; i<=padrange[1]; i++)
873               for(Int_t j=timerange[0]; j<=timerange[1]; j++)
874               fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fUsed = kTRUE;
875               */
876               fFitted++;
877             }
878           else //fit was too bad
879             {
880               if(fDebug)
881                 cout<<"Cluster fit was too bad"<<endl;
882               tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
883               fBadFitError++;
884               fFailed++;
885             }
886         }
887       n_overlaps++;
888     }
889   
890   delete [] plane;
891 }
892
893 void AliL3ClusterFitter::SetClusterfitFalse(AliL3ModelTrack *track)
894 {
895   //Cluster fit failed, so set the clusters to all the participating
896   //tracks to zero.
897   
898   Int_t i=-1;
899   Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
900   while(i < track->GetNOverlaps(fCurrentPadRow))
901     {
902       AliL3ModelTrack *tr=0;
903       if(i==-1)
904         tr = track;
905       else
906         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[i]);
907       i++;
908       if(!tr) continue;
909       
910       tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
911     }
912 }
913
914
915 void AliL3ClusterFitter::AddClusters()
916 {
917   if(!fClusters)
918     {
919       fClusters = new AliL3SpacePointData[fNMaxClusters];
920       fNClusters=0;
921     }
922   
923   if(fDebug)
924     cout<<"Writing cluster in slice "<<fSlice<<" patch "<<fPatch<<endl;
925   
926   AliL3TrackArray *tracks=0;
927   if(fSeeding==kTRUE)
928     tracks = fSeeds;
929   else
930     tracks = fTracks;
931   
932   if(!tracks)
933     return;
934   
935   for(Int_t i=0; i<tracks->GetNTracks(); i++)
936     {
937       AliL3ModelTrack *tr = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
938       if(!tr) continue;
939       
940       UInt_t *hitids = tr->GetHitNumbers();
941       Int_t nhits = tr->GetNHits();
942       for(Int_t i=fRowMax; i>=fRowMin; i--)
943         {
944           if(fSeeding)
945             if(tr->GetClusterModel(i)->fSlice != fSlice) continue;
946           if(!tr->IsPresent(i)) continue;
947           fCurrentPadRow = i;
948           Float_t pad,time,xywidth,zwidth;
949           Int_t charge;
950           tr->GetPad(i,pad);
951           tr->GetTime(i,time);
952           tr->GetClusterCharge(i,charge);
953
954           if(pad < -1 || pad >= AliL3Transform::GetNPads(i) || 
955              time < -1 || time >= AliL3Transform::GetNTimeBins())
956             {
957               continue;
958               cout<<"slice "<<fSlice<<" row "<<i<<" pad "<<pad<<" time "<<time<<endl;
959               tr->Print();
960               exit(5);
961             }
962
963           tr->CalculateClusterWidths(i,kTRUE); //Parametrize errors
964           
965           tr->GetXYWidth(i,xywidth);
966           tr->GetZWidth(i,zwidth);
967           Float_t xyz[3];
968           Int_t sector,row;
969           AliL3Transform::Slice2Sector(fSlice,i,sector,row);
970           
971           AliL3Transform::Raw2Global(xyz,sector,row,pad,time);
972           
973           if(fNClusters >= fNMaxClusters)
974             {
975               cerr<<"AliL3ClusterFitter::AddClusters : Too many clusters "<<fNClusters<<endl;
976               exit(5);
977             }
978           fClusters[fNClusters].fX = xyz[0];
979           fClusters[fNClusters].fY = xyz[1];
980           fClusters[fNClusters].fZ = xyz[2];
981           fClusters[fNClusters].fCharge = charge;
982           fClusters[fNClusters].fPadRow = i;
983           Int_t pa = AliL3Transform::GetPatch(i);
984           if(xywidth==0 || zwidth==0)
985             cerr<<"AliL3ClusterFitter::AddClusters : Cluster with zero width"<<endl;
986           if(xywidth>0)
987             fClusters[fNClusters].fSigmaY2 = xywidth*pow(AliL3Transform::GetPadPitchWidth(pa),2);
988           else
989             fClusters[fNClusters].fSigmaY2 = 1;
990           if(zwidth>0)
991             fClusters[fNClusters].fSigmaZ2 = zwidth*pow(AliL3Transform::GetZWidth(),2);
992           else
993             fClusters[fNClusters].fSigmaZ2 = 1;
994           Int_t pat=fPatch;
995           if(fPatch==-1)
996             pat=0;
997           fClusters[fNClusters].fID = fNClusters + ((fSlice&0x7f)<<25)+((pat&0x7)<<22);
998           
999           if(nhits >= 159)
1000             {
1001               cerr<<"AliL3ClusterFitter::AddClusters : Cluster counter of out range "<<nhits<<endl;
1002               exit(5);
1003             }
1004           hitids[nhits++] = fClusters[fNClusters].fID;
1005           
1006 #ifdef do_mc
1007           Int_t trackID[3];
1008           Int_t fpad = (Int_t)rint(pad);
1009           Int_t ftime = (Int_t)rint(time);
1010           if(fpad < 0)
1011             fpad=0;
1012           if(fpad >= AliL3Transform::GetNPads(i))
1013             fpad = AliL3Transform::GetNPads(i)-1;
1014           if(ftime<0)
1015             ftime=0;
1016           if(ftime >= AliL3Transform::GetNTimeBins())
1017             ftime = AliL3Transform::GetNTimeBins()-1;
1018           GetTrackID(fpad,ftime,trackID);
1019           fClusters[fNClusters].fTrackID[0] = trackID[0];
1020           fClusters[fNClusters].fTrackID[1] = trackID[1];
1021           fClusters[fNClusters].fTrackID[2] = trackID[2];
1022 #endif  
1023           //cout<<"Setting id "<<trackID[0]<<" on pad "<<pad<<" time "<<time<<" row "<<i<<endl;
1024           fNClusters++;
1025         }
1026       
1027       //Copy back the number of assigned clusters
1028       tr->SetNHits(nhits);
1029     }
1030 }
1031
1032 void AliL3ClusterFitter::WriteTracks()
1033 {
1034   if(!fSeeds)
1035     return;
1036   
1037   AliL3Compress *c = new AliL3Compress(-1,-1,fPath);
1038   c->WriteFile(fSeeds,"tracks_after.raw");
1039   delete c;
1040   
1041   Int_t clustercount=0;
1042   for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
1043     {
1044       AliL3ModelTrack *tr = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
1045       if(!tr) continue;
1046       if(tr->GetNHits()==0)
1047         fSeeds->Remove(i);
1048       clustercount += tr->GetNHits();
1049       /*
1050         if(tr->GetPt() > 1 && tr->GetNPresentClusters() < 150) 
1051         tr->Print();
1052       */
1053     }
1054   cout<<"Writing "<<clustercount<<" clusters"<<endl;
1055   fSeeds->Compress();
1056   AliL3MemHandler mem;
1057   Char_t filename[1024];
1058   sprintf(filename,"%s/fitter/tracks_0.raw",fPath);
1059   mem.SetBinaryOutput(filename);
1060   mem.TrackArray2Binary(fSeeds);
1061   mem.CloseBinaryOutput();
1062   
1063 }
1064
1065 void AliL3ClusterFitter::WriteClusters()
1066 {
1067   AliL3MemHandler mem;
1068   if(fDebug)
1069     cout<<"Write "<<fNClusters<<" clusters to file"<<endl;
1070   Char_t filename[1024];
1071   sprintf(filename,"%s/fitter/points_0_%d_%d.raw",fPath,fSlice,fPatch);
1072   mem.SetBinaryOutput(filename);
1073   mem.Memory2Binary(fNClusters,fClusters);
1074   mem.CloseBinaryOutput();
1075   mem.Free();
1076   
1077   delete [] fClusters;
1078   fClusters=0;
1079   fNClusters=0;
1080 }