Added support for NEWIO, merged cern-hlt tree, updated to latest track candidate...
[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]=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]=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   tfile.Binary2TrackArray(fSeeds);
242   tfile.CloseBinaryInput();
243
244   //if(!offline)
245   //fSeeds->QSort();
246   
247   Int_t clustercount=0;
248   for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
249     {
250       AliL3ModelTrack *track = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
251       if(!track) continue;
252
253       if(!offline)
254         {
255           if(i==0) cerr<<"AliL3ClusterFitter::LoadSeeds : Cutting on pt of 4GeV!!"<<endl;
256           if(track->GetPt() > 4.) 
257             {
258               fSeeds->Remove(i);
259               continue;
260             }
261         }
262       clustercount += track->GetNHits();
263       track->CalculateHelix();
264       
265       Int_t nhits = track->GetNHits();
266       UInt_t *hitids = track->GetHitNumbers();
267
268       Int_t origslice = (hitids[nhits-1]>>25)&0x7f;//Slice of innermost point
269
270       track->Init(origslice,-1);
271       Int_t slice = origslice;
272       
273       //for(Int_t j=rowrange[1]; j>=rowrange[0]; j--)
274       for(Int_t j=rowrange[0]; j<=rowrange[1]; j++)
275         {
276           
277           //Calculate the crossing point between track and padrow
278           Float_t angle = 0; //Perpendicular to padrow in local coordinates
279           AliL3Transform::Local2GlobalAngle(&angle,slice);
280           if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
281             {
282               //cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
283               continue;
284               //track->Print();
285               //exit(5);
286             }
287           Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
288           
289           Int_t sector,row;
290           AliL3Transform::Slice2Sector(slice,j,sector,row);
291           AliL3Transform::Global2Raw(xyz_cross,sector,row);
292           //cout<<"Examining slice "<<slice<<" row "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
293           if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //Track leaves the slice
294             {
295             newslice:
296               
297               Int_t tslice=slice;
298               Float_t lastcross=xyz_cross[1];
299               if(xyz_cross[1] > 0)
300                 {
301                   if(slice == 17)
302                     slice=0;
303                   else if(slice == 35)
304                     slice = 18;
305                   else
306                     slice += 1;
307                 }
308               else
309                 {
310                   if(slice == 0)
311                     slice = 17;
312                   else if(slice==18)
313                     slice = 35;
314                   else
315                     slice -= 1;
316                 }
317               if(slice < 0 || slice>35)
318                 {
319                   cerr<<"Wrong slice "<<slice<<" on row "<<j<<endl;
320                   exit(5);
321                 }
322               //cout<<"Track leaving, trying slice "<<slice<<endl;
323               angle=0;
324               AliL3Transform::Local2GlobalAngle(&angle,slice);
325               if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
326                 {
327                   cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
328                   continue;
329                   //track->Print();
330                   //exit(5);
331                 }
332               xyz_cross[0] = track->GetPointX();
333               xyz_cross[1] = track->GetPointY();
334               xyz_cross[2] = track->GetPointZ();
335               Int_t sector,row;
336               AliL3Transform::Slice2Sector(slice,j,sector,row);
337               AliL3Transform::Global2Raw(xyz_cross,sector,row);
338               if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //track is in the borderline
339                 {
340                   if(xyz_cross[1] > 0 && lastcross > 0 || xyz_cross[1] < 0 && lastcross < 0)
341                     goto newslice;
342                   else
343                     {
344                       slice = tslice;//Track is on the border of two slices
345                       continue;
346                     }
347                 }
348             }
349           
350           if(xyz_cross[2] < 0 || xyz_cross[2] >= AliL3Transform::GetNTimeBins())//track goes out of range
351             continue;
352           
353           if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j))
354             {
355               cerr<<"Slice "<<slice<<" padrow "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
356               track->Print();
357               exit(5);
358             }
359           
360           track->SetPadHit(j,xyz_cross[1]);
361           track->SetTimeHit(j,xyz_cross[2]);
362           angle=0;
363           AliL3Transform::Local2GlobalAngle(&angle,slice);
364           Float_t crossingangle = track->GetCrossingAngle(j,slice);
365           track->SetCrossingAngleLUT(j,crossingangle);
366           
367           track->CalculateClusterWidths(j);
368           
369           track->GetClusterModel(j)->fSlice = slice;
370           
371         }
372       memset(track->GetHitNumbers(),0,159*sizeof(UInt_t));//Reset the hitnumbers
373       track->SetNHits(0);
374     }
375   fSeeds->Compress();
376   
377   cout<<"Loaded "<<fSeeds->GetNTracks()<<" seeds and "<<clustercount<<" clusters"<<endl;
378 }
379
380 void AliL3ClusterFitter::FindClusters()
381 {
382   if(!fTracks)
383     {
384       cerr<<"AliL3ClusterFitter::Process : No tracks"<<endl;
385       return;
386     }
387   if(!fRowData)
388     {
389       cerr<<"AliL3ClusterFitter::Process : No data "<<endl;
390       return;
391     }
392   
393   AliL3DigitRowData *rowPt = fRowData;
394   AliL3DigitData *digPt=0;
395
396   Int_t pad,time;
397   Short_t charge;
398   
399   if(fRowMin < 0)
400     {
401       fRowMin = AliL3Transform::GetFirstRow(fPatch);
402       fRowMax = AliL3Transform::GetLastRow(fPatch);
403     }
404   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
405     {
406       if((Int_t)rowPt->fRow < fRowMin)
407         {
408           AliL3MemHandler::UpdateRowPointer(rowPt);
409           continue;
410         }
411       else if((Int_t)rowPt->fRow > fRowMax)
412         break;
413       else if((Int_t)rowPt->fRow != i)
414         {
415           cerr<<"AliL3ClusterFitter::FindClusters : Mismatching row numbering "<<i<<" "<<rowPt->fRow<<endl;
416           exit(5);
417         }
418       fCurrentPadRow = i;
419       memset((void*)fRow,0,(AliL3Transform::GetNTimeBins()+1)*(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
420       digPt = (AliL3DigitData*)rowPt->fDigitData;
421
422       for(UInt_t j=0; j<rowPt->fNDigit; j++)
423         {
424           pad = digPt[j].fPad;
425           time = digPt[j].fTime;
426           charge = digPt[j].fCharge;
427           fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge = charge;
428           fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kFALSE;
429           //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
430         }
431       
432       for(Int_t it=0; it<2; it++)
433         {
434           if(it==0)
435             {
436               fProcessTracks = fSeeds;
437               fSeeding = kTRUE;
438             }
439           else
440             {
441               fProcessTracks = fTracks;
442               fSeeding = kFALSE;
443             }
444           if(!fProcessTracks)
445             continue;
446           
447           for(Int_t k=0; k<fProcessTracks->GetNTracks(); k++)
448             {
449               AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(k);
450               if(!track) continue;
451               
452               if(fSeeding)
453                 if(track->GetClusterModel(i)->fSlice != fSlice) continue;
454               
455               if(track->GetPadHit(i) < 0 || track->GetPadHit(i) > AliL3Transform::GetNPads(i)-1 ||
456                  track->GetTimeHit(i) < 0 || track->GetTimeHit(i) > AliL3Transform::GetNTimeBins()-1)
457                 {
458                   track->SetCluster(i,0,0,0,0,0,0);
459                   continue;
460                 }
461               
462               if(CheckCluster(k) == kFALSE)
463                 fFailed++;
464             }
465         }
466       AliL3MemHandler::UpdateRowPointer(rowPt);
467     }
468   
469   fSeeding = kTRUE;
470   AddClusters();
471   fSeeding = kFALSE;
472   AddClusters();
473     
474   cout<<"Fitted "<<fFitted<<" clusters, failed "<<fFailed<<endl;
475   cout<<"Distribution:"<<endl;
476   cout<<"Bad fit "<<fBadFitError<<endl;
477   cout<<"Fit error "<<fFitError<<endl;
478   cout<<"Result error "<<fResultError<<endl;
479   cout<<"Fit range error "<<fFitRangeError<<endl;
480
481 }
482
483 Bool_t AliL3ClusterFitter::CheckCluster(Int_t trackindex)
484 {
485   //Check if this is a single or overlapping cluster
486   
487   AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(trackindex);
488   
489   Int_t row = fCurrentPadRow;
490   
491   if(track->IsSet(row)) //A fit has already be performed on this one
492     return kTRUE;
493   
494   //Define the cluster region of this hit:
495   Int_t padr[2]={999,-1};
496   Int_t timer[2]={999,-1};
497   
498   if(!SetFitRange(track,padr,timer))
499     {
500       track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
501       
502       if(fDebug)
503         cout<<"Failed to fit cluster at row "<<row<<" pad "<<(Int_t)rint(track->GetPadHit(row))<<" time "
504             <<(Int_t)rint(track->GetTimeHit(row))<<" hitcharge "
505             <<fRow[(AliL3Transform::GetNTimeBins()+1)*(Int_t)rint(track->GetPadHit(row))+(Int_t)rint(track->GetTimeHit(row))].fCharge<<endl;
506       fFitRangeError++;
507       return kFALSE;
508     }
509
510   //Check if any other track contributes to this cluster:
511   //This is done by checking if the tracks are overlapping within
512   //the range defined by the track parameters
513   
514   for(Int_t t=trackindex+1; t<fProcessTracks->GetNTracks(); t++)
515     {
516       AliL3ModelTrack *tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(t);
517       if(!tr) continue;
518       if(fSeeding)
519         if(tr->GetClusterModel(row)->fSlice != fSlice) continue;
520
521       Int_t xyw = (Int_t)ceil(sqrt(tr->GetParSigmaY2(row))*GetYWidthFactor()); 
522       Int_t zw = (Int_t)ceil(sqrt(tr->GetParSigmaZ2(row))*GetZWidthFactor()); 
523       
524       if( 
525          (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
526           tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
527          
528          (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
529           tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
530          
531          (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
532           tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1]) ||
533          
534          (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
535           tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1]) 
536          )
537         {
538           if(SetFitRange(tr,padr,timer)) //Expand the cluster fit range
539             track->SetOverlap(row,t);    //Set overlap
540         }
541     }
542
543   if(fDebug)
544     cout<<"Fitting cluster with "<<track->GetNOverlaps(fCurrentPadRow)<<" overlaps"<<endl;
545   FitClusters(track,padr,timer);
546   return kTRUE;
547 }
548
549 Bool_t AliL3ClusterFitter::SetFitRange(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
550 {
551   Int_t row = fCurrentPadRow;
552   Int_t nt = AliL3Transform::GetNTimeBins()+1;
553   
554   Int_t nsearchbins=0;
555   if(row < AliL3Transform::GetNRowLow())
556     nsearchbins=25;
557   else if(row < AliL3Transform::GetNRowLow() + AliL3Transform::GetNRowUp1())
558     nsearchbins=49;
559   else
560     nsearchbins=49;
561   
562   /*
563   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
564                        ,0,1,2,3,3,3,3,3,3,3
565                        ,2,1,0,-1,-2,-3
566                        ,-3,-3,-3,-3,-3,-3,-1,-1};
567   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
568                         ,-3,-3,-3,-3,-2,-1,0,1,2,3
569                         ,3,3,3,3,3,3,2,1,0,-1,-2,-3,-3,-3};
570   */
571   
572   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
573                        ,-3,3,-3,3,-3,3,-3,3,-3,3
574                        ,0,0,-1,-1,1,1,-2,-2,2,2,-3,-3,3,3};
575   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
576                         ,0,0,-1,-1,1,1,-2,-2,2,2
577                         ,-3,3,-3,3,-3,3,-3,3,-3,3,-3,3,-3,3};
578   
579   Int_t padhit = (Int_t)rint(track->GetPadHit(row));
580   Int_t timehit = (Int_t)rint(track->GetTimeHit(row));
581   Int_t padmax=-1;
582   Int_t timemax=-1;
583
584   for(Int_t index=0; index<nsearchbins; index++)
585     {
586       if(IsMaximum(padhit + padloop[index],timehit + timeloop[index])) 
587         {
588           padmax = padhit + padloop[index];
589           timemax = timehit + timeloop[index];
590           break;
591         }
592     }
593
594
595   //Define the cluster region of this hit:
596   //The region we look for, is centered at the local maxima
597   //and expanded around using the parametrized cluster width
598   //according to track parameters.
599   
600   Int_t xyw = (Int_t)ceil(sqrt(track->GetParSigmaY2(row))*GetYWidthFactor());
601   Int_t zw = (Int_t)ceil(sqrt(track->GetParSigmaZ2(row))*GetZWidthFactor());
602   
603   if(padmax>=0 && timemax>=0)
604     {
605       if(fDebug)
606         {
607           cout<<"Expanding cluster range using expected cluster widths: "<<xyw<<" "<<zw
608               <<" and setting local maxima pad "<<padmax<<" time "<<timemax<<endl;
609           if(xyw > 10 || zw > 10)
610             track->Print();
611         }
612       
613       //Set the hit to the local maxima of the cluster.
614       //Store the maxima in the cluster model structure,
615       //-only temporary, it will be overwritten when calling SetCluster.
616       
617       track->GetClusterModel(row)->fDPad = padmax;
618       track->GetClusterModel(row)->fDTime = timemax;
619
620       for(Int_t i=padmax-xyw; i<=padmax+xyw; i++)
621         {
622           for(Int_t j=timemax-zw; j<=timemax+zw; j++)
623             {
624               if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) continue;
625               if(fRow[nt*i+j].fCharge)
626                 {
627                   if(i < padrange[0]) padrange[0]=i;
628                   if(i > padrange[1]) padrange[1]=i;
629                   if(j < timerange[0]) timerange[0]=j;
630                   if(j > timerange[1]) timerange[1]=j;
631                 }
632             }
633         }
634       if(fDebug)
635         cout<<"New padrange "<<padrange[0]<<" "<<padrange[1]<<" "<<" time "<<timerange[0]<<" "<<timerange[1]<<endl;
636       return kTRUE;
637     }
638   return kFALSE;
639 }
640
641 Bool_t AliL3ClusterFitter::IsMaximum(Int_t pad,Int_t time)
642 {
643   if(pad<0 || pad >= AliL3Transform::GetNPads(fCurrentPadRow) ||
644      time<0 || time >= AliL3Transform::GetNTimeBins())
645     return kFALSE;
646   Int_t nt = AliL3Transform::GetNTimeBins()+1;
647   if(fRow[nt*pad+time].fUsed == kTRUE) return kFALSE; //Peak has been assigned before
648   Int_t charge = fRow[nt*pad+time].fCharge;
649   if(charge == 1023 || charge==0) return kFALSE;
650   
651   //fRow[nt*pad+time].fUsed = kTRUE;
652   //return kTRUE;
653
654   if(charge < fRow[nt*(pad-1)+(time-1)].fCharge) return kFALSE;
655   if(charge < fRow[nt*(pad)+(time-1)].fCharge) return kFALSE;
656   if(charge < fRow[nt*(pad+1)+(time-1)].fCharge) return kFALSE;
657   if(charge < fRow[nt*(pad-1)+(time)].fCharge) return kFALSE;
658   if(charge < fRow[nt*(pad+1)+(time)].fCharge) return kFALSE;
659   if(charge < fRow[nt*(pad-1)+(time+1)].fCharge) return kFALSE;
660   if(charge < fRow[nt*(pad)+(time+1)].fCharge) return kFALSE;
661   if(charge < fRow[nt*(pad+1)+(time+1)].fCharge) return kFALSE;
662   fRow[nt*pad+time].fUsed = kTRUE;
663   return kTRUE;
664 }
665
666 void AliL3ClusterFitter::FitClusters(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
667 {
668   //Handle single and overlapping clusters
669     
670   //Check whether this cluster has been set before:
671   
672   Int_t size = FIT_PTS;
673   Int_t max_tracks = FIT_MAXPAR/NUM_PARS;
674   if(track->GetNOverlaps(fCurrentPadRow) > max_tracks)
675     {
676       cerr<<"AliL3ClusterFitter::FitOverlappingClusters : Too many overlapping tracks"<<endl;
677       return;
678     }
679   Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
680   
681   //Check if at least one cluster is not already fitted
682   Bool_t all_fitted=kTRUE;
683   
684   Int_t k=-1;
685   while(k < track->GetNOverlaps(fCurrentPadRow))
686     {
687       AliL3ModelTrack *tr=0;
688       if(k==-1)
689         tr = track;
690       else
691         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
692       k++;
693       if(!tr) continue;
694       if(!tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow))//cluster has not been set and is not present
695         {
696           all_fitted = kFALSE;
697           break;
698         }
699     }
700   if(all_fitted)
701     {
702       if(fDebug)
703         cout<<"But all the clusters were already fitted on row "<<fCurrentPadRow<<endl;
704       return;
705     }
706   
707   //Allocate fit parameters array; this is interface to the C code
708   plane = new DPOINT[FIT_PTS];
709   memset(plane,0,FIT_PTS*sizeof(DPOINT));
710
711   Double_t x[FIT_PTS],y[FIT_PTS],s[FIT_PTS];
712   
713   //Fill the fit parameters:
714   Double_t a[FIT_MAXPAR];
715   Int_t lista[FIT_MAXPAR];
716   Double_t dev[FIT_MAXPAR],chisq_f;
717   
718   Int_t fit_pars=0;
719   
720   Int_t n_overlaps=0;
721   k=-1;
722   
723   //Fill the overlapping tracks:
724   while(k < track->GetNOverlaps(fCurrentPadRow))
725     {
726       AliL3ModelTrack *tr=0;
727       if(k==-1)
728         tr = track;
729       else
730         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
731       k++;
732       if(!tr) continue;
733       
734       if(tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow)) continue;//Cluster fit failed before
735       
736       //Use the local maxima as the input to the fitting routine.
737       //The local maxima is temporary stored in the cluster model:
738       Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad);  
739       Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime);
740       Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fCharge;
741       
742       if(fDebug)
743         cout<<"Fitting track cluster, pad "<<tr->GetPadHit(fCurrentPadRow)<<" time "
744             <<tr->GetTimeHit(fCurrentPadRow)<<" charge "<<charge<<" at local maxima in pad "<<hitpad
745             <<" time "<<hittime<<" xywidth "<<sqrt(tr->GetParSigmaY2(fCurrentPadRow))
746             <<" zwidth "<<sqrt(tr->GetParSigmaZ2(fCurrentPadRow))<<endl;
747       
748       if(charge==0)
749         {
750           cerr<<"Charge still zero!"<<endl;
751           exit(5);
752         }
753             
754       a[n_overlaps*NUM_PARS+2] = hitpad;
755       a[n_overlaps*NUM_PARS+4] = hittime;
756       
757       if(!tr->IsSet(fCurrentPadRow)) //Cluster is not fitted before
758         {
759           a[n_overlaps*NUM_PARS+1] = charge;
760           a[n_overlaps*NUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow)) * GetYWidthFactor();
761           a[n_overlaps*NUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
762           a[n_overlaps*NUM_PARS+6] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
763           lista[n_overlaps*NUM_PARS + 1] = 1;
764           lista[n_overlaps*NUM_PARS + 2] = 1;
765           lista[n_overlaps*NUM_PARS + 3] = 0;
766           lista[n_overlaps*NUM_PARS + 4] = 1;
767           lista[n_overlaps*NUM_PARS + 5] = 0;
768           lista[n_overlaps*NUM_PARS + 6] = 0;
769           fit_pars             += 3;
770         }
771       else  //Cluster was fitted before
772         {
773           if(!tr->IsPresent(fCurrentPadRow))
774             {
775               cerr<<"AliL3ClusterFitter::FindClusters : Cluster not present; there is a bug here"<<endl;
776               exit(5);
777             }
778           Int_t charge;
779           Float_t xywidth,zwidth,pad,time;
780           tr->GetPad(fCurrentPadRow,pad);
781           tr->GetTime(fCurrentPadRow,time);
782           tr->GetClusterCharge(fCurrentPadRow,charge);
783           xywidth = sqrt(tr->GetParSigmaY2(fCurrentPadRow));
784           zwidth = sqrt(tr->GetParSigmaZ2(fCurrentPadRow));
785           if(fDebug)
786             cout<<"Cluster had been fitted before, pad "<<pad<<" time "<<time<<" charge "<<charge<<" width "<<xywidth<<" "<<zwidth<<endl;
787           
788           a[n_overlaps*NUM_PARS+2] = pad;
789           a[n_overlaps*NUM_PARS+4] = time;
790           a[n_overlaps*NUM_PARS+1] = charge;
791           a[n_overlaps*NUM_PARS+3] = sqrt(xywidth) * GetYWidthFactor();
792           a[n_overlaps*NUM_PARS+5] = sqrt(zwidth) * GetZWidthFactor();
793           a[n_overlaps*NUM_PARS+6] = sqrt(zwidth) * GetZWidthFactor();
794
795           lista[n_overlaps*NUM_PARS + 1] = 1;
796           lista[n_overlaps*NUM_PARS + 2] = 0;
797           lista[n_overlaps*NUM_PARS + 3] = 0;
798           lista[n_overlaps*NUM_PARS + 4] = 0;
799           lista[n_overlaps*NUM_PARS + 5] = 0;
800           lista[n_overlaps*NUM_PARS + 6] = 0;
801           fit_pars             += 1;
802         }
803       n_overlaps++;
804     }
805   
806   if(n_overlaps==0) //No clusters here
807     {
808       delete [] plane;
809       return;
810     }
811
812   Int_t pad_num=0;
813   Int_t time_num_max=0;
814   Int_t ndata=0;
815   Int_t tot_charge=0;
816   if(fDebug)
817     cout<<"Padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "<<timerange[0]<<" "<<timerange[1]<<endl;
818   for(Int_t i=padrange[0]; i<=padrange[1]; i++)
819     {
820       Int_t max_charge = 0;
821       Int_t time_num=0;
822       for(Int_t j=timerange[0]; j<=timerange[1]; j++)
823         {
824           Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
825           
826           if(charge <= 0) continue;
827
828           time_num++;
829           if(charge > max_charge)
830             {
831               max_charge = charge;
832               //time_num++;
833             }
834           if(fDebug)
835             cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
836           tot_charge += charge;
837           ndata++;
838           if(ndata >= size)
839             {
840               cerr<<"Too many points; row "<<fCurrentPadRow<<" padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "
841                   <<timerange[0]<<" "<<timerange[1]<<endl;
842               exit(5);
843             }
844
845           plane[ndata].u = (Double_t)i;
846           plane[ndata].v = (Double_t)j;
847           x[ndata]=ndata;
848           y[ndata]=charge;
849           s[ndata]= 1 + sqrt((Double_t)charge);
850         }
851       if(max_charge) //there was charge on this pad
852         pad_num++;
853       if(time_num_max < time_num)
854         time_num_max = time_num;
855     }
856   
857   if(pad_num <= 1 || time_num_max <=1 || n_overlaps > fNmaxOverlaps || ndata <= fit_pars) //too few to do fit
858     {
859       SetClusterfitFalse(track);
860       if(fDebug)
861         cout<<"Too few digits or too many overlaps: "<<pad_num<<" "<<time_num_max<<" "<<n_overlaps<<" ndata "<<ndata<<" fit_pars "<<fit_pars<<endl;
862       delete [] plane;
863       return;
864     }
865
866   
867   Int_t npars = n_overlaps * NUM_PARS;
868   if(fDebug)
869     cout<<"Number of overlapping clusters "<<n_overlaps<<endl;
870   Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisq_f, f2gauss5 );
871   
872   if(ret<0)
873     {
874       SetClusterfitFalse(track);
875       fFailed++;
876       fFitError++;
877       delete [] plane;
878       return;
879       //exit(5);
880     }
881
882   chisq_f /= (ndata-fit_pars);
883   if(fDebug)
884     cout<<"Chisq "<<chisq_f<<endl;
885   
886   Bool_t overlapping=kFALSE;
887   if(track->GetNOverlaps(fCurrentPadRow) > 0)//There is a overlap
888     overlapping=kTRUE;
889   
890   k=-1;
891   n_overlaps=0;
892   while(k < track->GetNOverlaps(fCurrentPadRow))
893     {
894       AliL3ModelTrack *tr=0;
895       if(k==-1)
896         tr = track;
897       else
898         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
899       k++;
900       if(!tr) continue;
901       if(!tr->IsPresent(fCurrentPadRow))
902         {
903           if(tr->IsSet(fCurrentPadRow)) continue;//This cluster has been set before
904           
905           if(chisq_f < fChiSqMax[(Int_t)overlapping])//cluster fit is good enough
906             {
907               tot_charge = (Int_t)(a[n_overlaps*NUM_PARS+1] * a[n_overlaps*NUM_PARS+3] * a[n_overlaps*NUM_PARS+5]);
908               Float_t fpad = a[n_overlaps*NUM_PARS+2];
909               Float_t ftime = a[n_overlaps*NUM_PARS+4];
910               if(tot_charge < 0 || fpad < -1 || fpad > AliL3Transform::GetNPads(fCurrentPadRow) || 
911                  ftime < -1 || ftime > AliL3Transform::GetNTimeBins())
912                 {
913                   if(fDebug)
914                     cout<<"AliL3ClusterFitter::Fatal result(s) in fit; in slice "<<fSlice<<" row "<<fCurrentPadRow
915                         <<"; pad "<<fpad<<" time "<<ftime<<" charge "<<tot_charge<<" xywidth "<<a[n_overlaps*NUM_PARS+3]
916                         <<" zwidth "<<a[n_overlaps*NUM_PARS+5]<<" peakcharge "<<a[n_overlaps*NUM_PARS+1]<<endl;
917                   tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
918                   fFailed++;
919                   fResultError++;
920                   continue;
921                 }
922               
923               tr->SetCluster(fCurrentPadRow,fpad,ftime,tot_charge,0,0,pad_num);
924               if(fDebug)
925                 cout<<"Setting cluster in pad "<<a[n_overlaps*NUM_PARS+2]<<" time "<<a[n_overlaps*NUM_PARS+4]<<" charge "<<tot_charge<<endl;
926               /*
927               //Set the digits to used:
928               for(Int_t i=padrange[0]; i<=padrange[1]; i++)
929               for(Int_t j=timerange[0]; j<=timerange[1]; j++)
930               fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fUsed = kTRUE;
931               */
932               fFitted++;
933             }
934           else //fit was too bad
935             {
936               if(fDebug)
937                 cout<<"Cluster fit was too bad"<<endl;
938               tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
939               fBadFitError++;
940               fFailed++;
941             }
942         }
943       n_overlaps++;
944     }
945   
946   delete [] plane;
947 }
948
949 void AliL3ClusterFitter::SetClusterfitFalse(AliL3ModelTrack *track)
950 {
951   //Cluster fit failed, so set the clusters to all the participating
952   //tracks to zero.
953   
954   Int_t i=-1;
955   Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
956   while(i < track->GetNOverlaps(fCurrentPadRow))
957     {
958       AliL3ModelTrack *tr=0;
959       if(i==-1)
960         tr = track;
961       else
962         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[i]);
963       i++;
964       if(!tr) continue;
965       
966       //Set the digit data to unused, so it can be fitted to another bastard:
967       Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad);  
968       Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime);
969       fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fUsed = kFALSE;
970       
971       tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
972     }
973
974 }
975
976
977 void AliL3ClusterFitter::AddClusters()
978 {
979   if(!fClusters)
980     {
981       fClusters = new AliL3SpacePointData[fNMaxClusters];
982       fNClusters=0;
983     }
984   
985   if(fDebug)
986     cout<<"Writing cluster in slice "<<fSlice<<" patch "<<fPatch<<endl;
987   
988   AliL3TrackArray *tracks=0;
989   if(fSeeding==kTRUE)
990     tracks = fSeeds;
991   else
992     tracks = fTracks;
993   
994   if(!tracks)
995     return;
996   
997   for(Int_t i=0; i<tracks->GetNTracks(); i++)
998     {
999       AliL3ModelTrack *tr = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
1000       if(!tr) continue;
1001       
1002       UInt_t *hitids = tr->GetHitNumbers();
1003       Int_t nhits = tr->GetNHits();
1004       for(Int_t i=fRowMax; i>=fRowMin; i--)
1005         {
1006           if(fSeeding)
1007             if(tr->GetClusterModel(i)->fSlice != fSlice) continue;
1008           if(!tr->IsPresent(i)) continue;
1009           fCurrentPadRow = i;
1010           Float_t pad,time,xywidth,zwidth;
1011           Int_t charge;
1012           tr->GetPad(i,pad);
1013           tr->GetTime(i,time);
1014           tr->GetClusterCharge(i,charge);
1015
1016           if(pad < -1 || pad >= AliL3Transform::GetNPads(i) || 
1017              time < -1 || time >= AliL3Transform::GetNTimeBins())
1018             {
1019               continue;
1020               cout<<"slice "<<fSlice<<" row "<<i<<" pad "<<pad<<" time "<<time<<endl;
1021               tr->Print();
1022               exit(5);
1023             }
1024
1025           tr->CalculateClusterWidths(i,kTRUE); //Parametrize errors
1026           
1027           tr->GetXYWidth(i,xywidth);
1028           tr->GetZWidth(i,zwidth);
1029           Float_t xyz[3];
1030           Int_t sector,row;
1031           AliL3Transform::Slice2Sector(fSlice,i,sector,row);
1032           
1033           AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
1034           
1035           if(fNClusters >= fNMaxClusters)
1036             {
1037               cerr<<"AliL3ClusterFitter::AddClusters : Too many clusters "<<fNClusters<<endl;
1038               exit(5);
1039             }
1040           
1041           fClusters[fNClusters].fX = xyz[0];
1042           fClusters[fNClusters].fY = xyz[1];
1043           fClusters[fNClusters].fZ = xyz[2];
1044           fClusters[fNClusters].fCharge = charge;
1045           fClusters[fNClusters].fPadRow = i;
1046           Int_t pa = AliL3Transform::GetPatch(i);
1047           if(xywidth==0 || zwidth==0)
1048             cerr<<"AliL3ClusterFitter::AddClusters : Cluster with zero width"<<endl;
1049           if(xywidth>0)
1050             fClusters[fNClusters].fSigmaY2 = xywidth*pow(AliL3Transform::GetPadPitchWidth(pa),2);
1051           else
1052             fClusters[fNClusters].fSigmaY2 = 1;
1053           if(zwidth>0)
1054             fClusters[fNClusters].fSigmaZ2 = zwidth*pow(AliL3Transform::GetZWidth(),2);
1055           else
1056             fClusters[fNClusters].fSigmaZ2 = 1;
1057           Int_t pat=fPatch;
1058           if(fPatch==-1)
1059             pat=0;
1060           fClusters[fNClusters].fID = fNClusters + ((fSlice&0x7f)<<25)+((pat&0x7)<<22);
1061           
1062           if(nhits >= AliL3Transform::GetNRows())
1063             {
1064               cerr<<"AliL3ClusterFitter::AddClusters : Cluster counter of out range "<<nhits<<endl;
1065               exit(5);
1066             }
1067           hitids[nhits++] = fClusters[fNClusters].fID;
1068           
1069 #ifdef do_mc
1070           Int_t trackID[3];
1071           Int_t fpad = (Int_t)rint(pad);
1072           Int_t ftime = (Int_t)rint(time);
1073           if(fpad < 0)
1074             fpad=0;
1075           if(fpad >= AliL3Transform::GetNPads(i))
1076             fpad = AliL3Transform::GetNPads(i)-1;
1077           if(ftime<0)
1078             ftime=0;
1079           if(ftime >= AliL3Transform::GetNTimeBins())
1080             ftime = AliL3Transform::GetNTimeBins()-1;
1081           GetTrackID(fpad,ftime,trackID);
1082           fClusters[fNClusters].fTrackID[0] = trackID[0];
1083           fClusters[fNClusters].fTrackID[1] = trackID[1];
1084           fClusters[fNClusters].fTrackID[2] = trackID[2];
1085 #endif  
1086           //cout<<"Setting id "<<trackID[0]<<" on pad "<<pad<<" time "<<time<<" row "<<i<<endl;
1087           fNClusters++;
1088         }
1089       
1090       //Copy back the number of assigned clusters
1091       tr->SetNHits(nhits);
1092       
1093     }
1094 }
1095
1096 void AliL3ClusterFitter::WriteTracks(Int_t min_hits)
1097 {
1098   if(!fSeeds)
1099     return;
1100   
1101   AliL3TrackArray *fakes = new AliL3TrackArray();
1102   
1103   Int_t clustercount=0;
1104   for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
1105     {
1106       AliL3ModelTrack *tr = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
1107       if(!tr) continue;
1108       if(tr->GetNHits() < min_hits)
1109         {
1110           fakes->AddLast(tr);
1111           fSeeds->Remove(i);
1112         }
1113       clustercount += tr->GetNHits();
1114     }
1115   
1116   cout<<"Writing "<<clustercount<<" clusters"<<endl;
1117   fSeeds->Compress();
1118   AliL3MemHandler mem;
1119   Char_t filename[1024];
1120   sprintf(filename,"%s/fitter/tracks_%d.raw",fPath,fEvent);
1121   mem.SetBinaryOutput(filename);
1122   mem.TrackArray2Binary(fSeeds);
1123   mem.CloseBinaryOutput();
1124   
1125   //Write the fake tracks to its own file
1126   mem.Free();
1127   sprintf(filename,"%s/fitter/tracks_fakes_%d.raw",fPath,fEvent);
1128   mem.SetBinaryOutput(filename);
1129   mem.TrackArray2Binary(fakes);
1130   mem.CloseBinaryOutput();
1131   delete fakes;
1132 }
1133
1134 void AliL3ClusterFitter::WriteClusters(Bool_t global)
1135 {
1136   AliL3MemHandler mem;
1137   if(fDebug)
1138     cout<<"Write "<<fNClusters<<" clusters to file"<<endl;
1139   Char_t filename[1024];
1140   sprintf(filename,"%s/fitter/points_%d_%d_%d.raw",fPath,fEvent,fSlice,fPatch);
1141   mem.SetBinaryOutput(filename);
1142   if(global == kTRUE)
1143     mem.Transform(fNClusters,fClusters,fSlice);
1144   mem.Memory2Binary(fNClusters,fClusters);
1145   mem.CloseBinaryOutput();
1146   mem.Free();
1147   
1148   delete [] fClusters;
1149   fClusters=0;
1150   fNClusters=0;
1151 }
1152