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