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