Corrected casting
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //-------------------------------------------------------
19 //          Implementation of the TPC clusterer
20 //
21 //   Origin: Marian Ivanov 
22 //-------------------------------------------------------
23
24 #include "AliTPCReconstructor.h"
25 #include "AliTPCclustererMI.h"
26 #include "AliTPCclusterMI.h"
27 #include <TObjArray.h>
28 #include <TFile.h>
29 #include "AliTPCClustersArray.h"
30 #include "AliTPCClustersRow.h"
31 #include "AliTPCRawStream.h"
32 #include "AliDigits.h"
33 #include "AliSimDigits.h"
34 #include "AliTPCParam.h"
35 #include "AliRawReader.h"
36 #include "AliTPCRawStream.h"
37 #include "AliRunLoader.h"
38 #include "AliLoader.h"
39 #include "Riostream.h"
40 #include <TTree.h>
41
42 ClassImp(AliTPCclustererMI)
43
44
45
46 AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par)
47 {
48   fInput =0;
49   fOutput=0;
50   fParam = par;
51 }
52 void AliTPCclustererMI::SetInput(TTree * tree)
53 {
54   //
55   // set input tree with digits
56   //
57   fInput = tree;  
58   if  (!fInput->GetBranch("Segment")){
59     cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
60     fInput=0;
61     return;
62   }
63 }
64
65 void AliTPCclustererMI::SetOutput(TTree * tree) 
66 {
67   //
68   //
69   fOutput= tree;  
70   AliTPCClustersRow clrow;
71   AliTPCClustersRow *pclrow=&clrow;  
72   clrow.SetClass("AliTPCclusterMI");
73   clrow.SetArray(1); // to make Clones array
74   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);    
75 }
76
77
78 Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
79   // sigma y2 = in digits  - we don't know the angle
80   Float_t z = iz*fParam->GetZWidth();
81   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
82     (fPadWidth*fPadWidth);
83   Float_t sres = 0.25;
84   Float_t res = sd2+sres;
85   return res;
86 }
87
88
89 Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
90   //sigma z2 = in digits - angle estimated supposing vertex constraint
91   Float_t z = iz*fZWidth;
92   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
93   Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth);
94   angular*=angular;
95   angular/=12.;
96   Float_t sres = fParam->GetZSigma()/fZWidth;
97   sres *=sres;
98   Float_t res = angular +sd2+sres;
99   return res;
100 }
101
102 void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Int_t *bins, UInt_t /*m*/,
103 AliTPCclusterMI &c) 
104 {
105   Int_t i0=k/max;  //central pad
106   Int_t j0=k%max;  //central time bin
107
108   // set pointers to data
109   //Int_t dummy[5] ={0,0,0,0,0};
110   Int_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
111   Int_t * resmatrix[5];
112   for (Int_t di=-2;di<=2;di++){
113     matrix[di+2]  =  &bins[k+di*max];
114     resmatrix[di+2]  =  &fResBins[k+di*max];
115   }
116   //build matrix with virtual charge
117   Float_t sigmay2= GetSigmaY2(j0);
118   Float_t sigmaz2= GetSigmaZ2(j0);
119
120   Float_t vmatrix[5][5];
121   vmatrix[2][2] = matrix[2][0];
122   c.SetType(0);
123   c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
124   for (Int_t di =-1;di <=1;di++)
125     for (Int_t dj =-1;dj <=1;dj++){
126       Float_t amp = matrix[di+2][dj];
127       if ( (amp<2) && (fLoop<2)){
128         // if under threshold  - calculate virtual charge
129         Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
130         amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
131         if (amp>2) amp = 2;
132         vmatrix[2+di][2+dj]=amp;
133         vmatrix[2+2*di][2+2*dj]=0;
134         if ( (di*dj)!=0){       
135           //DIAGONAL ELEMENTS
136           vmatrix[2+2*di][2+dj] =0;
137           vmatrix[2+di][2+2*dj] =0;
138         }
139         continue;
140       }
141       if (amp<4){
142         //if small  amplitude - below  2 x threshold  - don't consider other one        
143         vmatrix[2+di][2+dj]=amp;
144         vmatrix[2+2*di][2+2*dj]=0;  // don't take to the account next bin
145         if ( (di*dj)!=0){       
146           //DIAGONAL ELEMENTS
147           vmatrix[2+2*di][2+dj] =0;
148           vmatrix[2+di][2+2*dj] =0;
149         }
150         continue;
151       }
152       //if bigger then take everything
153       vmatrix[2+di][2+dj]=amp;
154       vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;      
155       if ( (di*dj)!=0){       
156           //DIAGONAL ELEMENTS
157           vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
158           vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
159         }      
160     }
161
162
163   
164   Float_t sumw=0;
165   Float_t sumiw=0;
166   Float_t sumi2w=0;
167   Float_t sumjw=0;
168   Float_t sumj2w=0;
169   //
170   for (Int_t i=-2;i<=2;i++)
171     for (Int_t j=-2;j<=2;j++){
172       Float_t amp = vmatrix[i+2][j+2];
173
174       sumw    += amp;
175       sumiw   += i*amp;
176       sumi2w  += i*i*amp;
177       sumjw   += j*amp;
178       sumj2w  += j*j*amp;
179     }    
180   //   
181   Float_t meani = sumiw/sumw;
182   Float_t mi2   = sumi2w/sumw-meani*meani;
183   Float_t meanj = sumjw/sumw;
184   Float_t mj2   = sumj2w/sumw-meanj*meanj;
185   //
186   Float_t ry = mi2/sigmay2;
187   Float_t rz = mj2/sigmaz2;
188   
189   //
190   if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
191   if ( (ry <1.2) && (rz<1.2) ) {
192     //if cluster looks like expected 
193     //+1.2 deviation from expected sigma accepted
194     //    c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
195
196     meani +=i0;
197     meanj +=j0;
198     //set cluster parameters
199     c.SetQ(sumw);
200     c.SetY(meani*fPadWidth); 
201     c.SetZ(meanj*fZWidth); 
202     c.SetSigmaY2(mi2);
203     c.SetSigmaZ2(mj2);
204     AddCluster(c);
205     //remove cluster data from data
206     for (Int_t di=-2;di<=2;di++)
207       for (Int_t dj=-2;dj<=2;dj++){
208         resmatrix[di+2][dj] -= Int_t(vmatrix[di+2][dj+2]);
209         if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0;
210       }
211     resmatrix[2][0] =0;
212     return;     
213   }
214   //
215   //unfolding when neccessary  
216   //
217   
218   Int_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
219   Int_t dummy[7]={0,0,0,0,0,0};
220   for (Int_t di=-3;di<=3;di++){
221     matrix2[di+3] =  &bins[k+di*max];
222     if ((k+di*max)<3)  matrix2[di+3] = &dummy[3];
223     if ((k+di*max)>fMaxBin-3)  matrix2[di+3] = &dummy[3];
224   }
225   Float_t vmatrix2[5][5];
226   Float_t sumu;
227   Float_t overlap;
228   UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
229   //
230   //  c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
231   meani +=i0;
232   meanj +=j0;
233   //set cluster parameters
234   c.SetQ(sumu);
235   c.SetY(meani*fPadWidth); 
236   c.SetZ(meanj*fZWidth); 
237   c.SetSigmaY2(mi2);
238   c.SetSigmaZ2(mj2);
239   c.SetType(Char_t(overlap)+1);
240   AddCluster(c);
241
242   //unfolding 2
243   meani-=i0;
244   meanj-=j0;
245   if (gDebug>4)
246     printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
247 }
248
249
250
251 void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj, 
252                                       Float_t & sumu, Float_t & overlap )
253 {
254   //
255   //unfold cluster from input matrix
256   //data corresponding to cluster writen in recmatrix
257   //output meani and meanj
258
259   //take separatelly y and z
260
261   Float_t sum3i[7] = {0,0,0,0,0,0,0};
262   Float_t sum3j[7] = {0,0,0,0,0,0,0};
263
264   for (Int_t k =0;k<7;k++)
265     for (Int_t l = -1; l<=1;l++){
266       sum3i[k]+=matrix2[k][l];
267       sum3j[k]+=matrix2[l+3][k-3];
268     }
269   Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
270   //
271   //unfold  y 
272   Float_t sum3wi    = 0;  //charge minus overlap
273   Float_t sum3wio   = 0;  //full charge
274   Float_t sum3iw    = 0;  //sum for mean value
275   for (Int_t dk=-1;dk<=1;dk++){
276     sum3wio+=sum3i[dk+3];
277     if (dk==0){
278       sum3wi+=sum3i[dk+3];     
279     }
280     else{
281       Float_t ratio =1;
282       if (  ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
283            sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
284         Float_t xm2 = sum3i[-dk+3];
285         Float_t xm1 = sum3i[+3];
286         Float_t x1  = sum3i[2*dk+3];
287         Float_t x2  = sum3i[3*dk+3];    
288         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);      
289         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
290         ratio = w11/(w11+w12);   
291         for (Int_t dl=-1;dl<=1;dl++)
292           mratio[dk+1][dl+1] *= ratio;
293       }
294       Float_t amp = sum3i[dk+3]*ratio;
295       sum3wi+=amp;
296       sum3iw+= dk*amp;      
297     }
298   }
299   meani = sum3iw/sum3wi;
300   Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
301
302
303
304   //unfold  z 
305   Float_t sum3wj    = 0;  //charge minus overlap
306   Float_t sum3wjo   = 0;  //full charge
307   Float_t sum3jw    = 0;  //sum for mean value
308   for (Int_t dk=-1;dk<=1;dk++){
309     sum3wjo+=sum3j[dk+3];
310     if (dk==0){
311       sum3wj+=sum3j[dk+3];     
312     }
313     else{
314       Float_t ratio =1;
315       if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
316            (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
317         Float_t xm2 = sum3j[-dk+3];
318         Float_t xm1 = sum3j[+3];
319         Float_t x1  = sum3j[2*dk+3];
320         Float_t x2  = sum3j[3*dk+3];    
321         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);      
322         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
323         ratio = w11/(w11+w12);   
324         for (Int_t dl=-1;dl<=1;dl++)
325           mratio[dl+1][dk+1] *= ratio;
326       }
327       Float_t amp = sum3j[dk+3]*ratio;
328       sum3wj+=amp;
329       sum3jw+= dk*amp;      
330     }
331   }
332   meanj = sum3jw/sum3wj;
333   Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;  
334   overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);  
335   sumu = (sum3wj+sum3wi)/2.;
336   
337   if (overlap ==3) {
338     //if not overlap detected remove everything
339     for (Int_t di =-2; di<=2;di++)
340       for (Int_t dj =-2; dj<=2;dj++){
341         recmatrix[di+2][dj+2] = matrix2[3+di][dj];
342       }
343   }
344   else{
345     for (Int_t di =-1; di<=1;di++)
346       for (Int_t dj =-1; dj<=1;dj++){
347         Float_t ratio =1;
348         if (mratio[di+1][dj+1]==1){
349           recmatrix[di+2][dj+2]     = matrix2[3+di][dj];
350           if (TMath::Abs(di)+TMath::Abs(dj)>1){
351             recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
352             recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
353           }       
354           recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
355         }
356         else
357           {
358             //if we have overlap in direction
359             recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];    
360             if (TMath::Abs(di)+TMath::Abs(dj)>1){
361               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
362               recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
363               //
364               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
365               recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
366             }
367             else{
368               ratio =  recmatrix[di+2][dj+2]/matrix2[3][0];
369               recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
370             }
371           }
372       }
373   }
374   if (gDebug>4) 
375     printf("%f\n", recmatrix[2][2]);
376   
377 }
378
379 Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
380 {
381   //
382   // estimate max
383   Float_t sumteor= 0;
384   Float_t sumamp = 0;
385
386   for (Int_t di = -1;di<=1;di++)
387     for (Int_t dj = -1;dj<=1;dj++){
388       if (vmatrix[2+di][2+dj]>2){
389         Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
390         sumteor += teor*vmatrix[2+di][2+dj];
391         sumamp  += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
392       }
393     }    
394   Float_t max = sumamp/sumteor;
395   return max;
396 }
397
398 void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){
399   //
400   // transform cluster to the global coordinata
401   // add the cluster to the array
402   //
403   Float_t meani = c.GetY()/fPadWidth;
404   Float_t meanj = c.GetZ()/fZWidth;
405
406   Int_t ki = TMath::Nint(meani-3);
407   if (ki<0) ki=0;
408   if (ki>=fMaxPad) ki = fMaxPad-1;
409   Int_t kj = TMath::Nint(meanj-3);
410   if (kj<0) kj=0;
411   if (kj>=fMaxTime-3) kj=fMaxTime-4;
412   // ki and kj shifted to "real" coordinata
413   if (fRowDig) {
414     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
415     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
416     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
417   }
418   
419   
420   Float_t s2 = c.GetSigmaY2();
421   Float_t w=fParam->GetPadPitchWidth(fSector);
422   
423   c.SetSigmaY2(s2*w*w);
424   s2 = c.GetSigmaZ2(); 
425   w=fZWidth;
426   c.SetSigmaZ2(s2*w*w);
427   c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
428   c.SetZ(fZWidth*(meanj-3)); 
429   c.SetZ(c.GetZ()   - 3.*fParam->GetZSigma()); // PASA delay 
430   c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
431   
432   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
433     //c.SetSigmaY2(c.GetSigmaY2()*25.);
434     //c.SetSigmaZ2(c.GetSigmaZ2()*4.);
435     c.SetType(-(c.GetType()+3));  //edge clusters
436   }
437   if (fLoop==2) c.SetType(100);
438
439   TClonesArray * arr = fRowCl->GetArray();
440   // AliTPCclusterMI * cl = 
441   new ((*arr)[fNcluster]) AliTPCclusterMI(c);
442
443   fNcluster++;
444 }
445
446
447 //_____________________________________________________________________________
448 void AliTPCclustererMI::Digits2Clusters()
449 {
450   //-----------------------------------------------------------------
451   // This is a simple cluster finder.
452   //-----------------------------------------------------------------
453
454   if (!fInput) { 
455     Error("Digits2Clusters", "input tree not initialised");
456     return;
457   }
458  
459   if (!fOutput) {
460     Error("Digits2Clusters", "output tree not initialised");
461     return;
462   }
463
464   AliSimDigits digarr, *dummy=&digarr;
465   fRowDig = dummy;
466   fInput->GetBranch("Segment")->SetAddress(&dummy);
467   Stat_t nentries = fInput->GetEntries();
468   
469   fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3   after
470     
471   Int_t nclusters  = 0;
472   
473   for (Int_t n=0; n<nentries; n++) {
474     fInput->GetEvent(n);
475     Int_t row;
476     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,row)) {
477       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
478       continue;
479     }
480         
481     AliTPCClustersRow *clrow= new AliTPCClustersRow();
482     fRowCl = clrow;
483     clrow->SetClass("AliTPCclusterMI");
484     clrow->SetArray(1);
485
486     clrow->SetID(digarr.GetID());
487     fOutput->GetBranch("Segment")->SetAddress(&clrow);
488     fRx=fParam->GetPadRowRadii(fSector,row);
489     
490     
491     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
492     fZWidth = fParam->GetZWidth();
493     if (fSector < kNIS) {
494       fMaxPad = fParam->GetNPadsLow(row);
495       fSign = (fSector < kNIS/2) ? 1 : -1;
496       fPadLength = fParam->GetPadPitchLength(fSector,row);
497       fPadWidth = fParam->GetPadPitchWidth();
498     } else {
499       fMaxPad = fParam->GetNPadsUp(row);
500       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
501       fPadLength = fParam->GetPadPitchLength(fSector,row);
502       fPadWidth  = fParam->GetPadPitchWidth();
503     }
504     
505     
506     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
507     fBins    =new Int_t[fMaxBin];
508     fResBins =new Int_t[fMaxBin];  //fBins with residuals after 1 finder loop 
509     memset(fBins,0,sizeof(Int_t)*fMaxBin);
510     
511     if (digarr.First()) //MI change
512       do {
513         Short_t dig=digarr.CurrentDigit();
514         if (dig<=fParam->GetZeroSup()) continue;
515         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
516         fBins[i*fMaxTime+j]=dig;
517       } while (digarr.Next());
518     digarr.ExpandTrackBuffer();
519
520     FindClusters();
521
522     fOutput->Fill();
523     delete clrow;    
524     nclusters+=fNcluster;    
525     delete[] fBins;      
526     delete[] fResBins;  
527   }  
528
529   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
530 }
531
532 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
533 {
534 //-----------------------------------------------------------------
535 // This is a cluster finder for raw data.
536 //-----------------------------------------------------------------
537
538   if (!fOutput) {
539     Error("Digits2Clusters", "output tree not initialised");
540     return;
541   }
542
543   rawReader->Reset();
544   AliTPCRawStream input(rawReader);
545
546   fRowDig = NULL;
547
548   Int_t nclusters  = 0;
549   
550   fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after
551   const Int_t kNIS = fParam->GetNInnerSector();
552   const Int_t kNOS = fParam->GetNOuterSector();
553   const Int_t kNS = kNIS + kNOS;
554   fZWidth = fParam->GetZWidth();
555   Int_t zeroSup = fParam->GetZeroSup();
556
557   fBins = NULL;
558   Int_t** splitRows = new Int_t* [kNS*2];
559   Int_t** splitRowsRes = new Int_t* [kNS*2];
560   for (Int_t iSector = 0; iSector < kNS*2; iSector++)
561     splitRows[iSector] = NULL;
562   Int_t iSplitRow = -1;
563
564   Bool_t next = kTRUE;
565   while (next) {
566     next = input.Next();
567
568     // when the sector or row number has changed ...
569     if (input.IsNewRow() || !next) {
570
571       // ... find clusters in the previous pad row, and ...
572       if (fBins) {
573         if ((iSplitRow < 0) || splitRows[fSector + kNS*iSplitRow]) {
574           fRowCl = new AliTPCClustersRow;
575           fRowCl->SetClass("AliTPCclusterMI");
576           fRowCl->SetArray(1);
577           fRowCl->SetID(fParam->GetIndex(fSector, input.GetPrevRow()));
578           fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
579
580           FindClusters();
581
582           fOutput->Fill();
583           delete fRowCl;    
584           nclusters += fNcluster;    
585           delete[] fBins;
586           delete[] fResBins;
587           if (iSplitRow >= 0) splitRows[fSector + kNS*iSplitRow] = NULL;
588
589         } else if (iSplitRow >= 0) {
590           splitRows[fSector + kNS*iSplitRow] = fBins;
591           splitRowsRes[fSector + kNS*iSplitRow] = fResBins;
592         }
593       }
594
595       if (!next) break;
596
597       // ... prepare for the next pad row
598       fSector = input.GetSector();
599       Int_t iRow = input.GetRow();
600       fRx = fParam->GetPadRowRadii(fSector, iRow);
601     
602       iSplitRow = -1;
603       if (fSector < kNIS) {
604         fMaxPad = fParam->GetNPadsLow(iRow);
605         fSign = (fSector < kNIS/2) ? 1 : -1;
606         if (iRow == 30) iSplitRow = 0;
607       } else {
608         fMaxPad = fParam->GetNPadsUp(iRow);
609         fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
610         if (iRow == 27) iSplitRow = 0;
611         else if (iRow == 76) iSplitRow = 1;
612       }
613       fPadLength = fParam->GetPadPitchLength(fSector, iRow);
614       fPadWidth  = fParam->GetPadPitchWidth();
615     
616       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
617       if ((iSplitRow < 0) || !splitRows[fSector + kNS*iSplitRow]) {
618         fBins    = new Int_t[fMaxBin];
619         fResBins = new Int_t[fMaxBin];  //fBins with residuals after 1 finder loop 
620         memset(fBins, 0, sizeof(Int_t)*fMaxBin);
621       } else {
622         fBins    = splitRows[fSector + kNS*iSplitRow];
623         fResBins = splitRowsRes[fSector + kNS*iSplitRow];
624       }
625     }
626
627     // fill fBins with digits data
628     if (input.GetSignal() <= zeroSup) continue;
629     Int_t i = input.GetPad() + 3;
630     Int_t j = input.GetTime() + 3;
631     fBins[i*fMaxTime+j] = input.GetSignal();
632   }
633
634   // find clusters in split rows that were skipped until now.
635   // this can happen if the rows were not splitted 
636   for (fSector = 0; fSector < kNS; fSector++)
637     for (Int_t iSplit = 0; iSplit < 2; iSplit++)
638       if (splitRows[fSector + kNS*iSplit]) {
639
640         Int_t iRow = -1;
641         if (fSector < kNIS) {
642           iRow = 30;
643           fMaxPad = fParam->GetNPadsLow(iRow);
644           fSign = (fSector < kNIS/2) ? 1 : -1;
645         } else {
646           if (iSplit == 0) iRow = 27; else iRow = 76;
647           fMaxPad = fParam->GetNPadsUp(iRow);
648           fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
649         }
650         fRx = fParam->GetPadRowRadii(fSector, iRow);
651         fPadLength = fParam->GetPadPitchLength(fSector, iRow);
652         fPadWidth  = fParam->GetPadPitchWidth();
653
654         fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
655         fBins    = splitRows[fSector + kNS*iSplit];
656         fResBins = splitRowsRes[fSector + kNS*iSplit];
657
658         fRowCl = new AliTPCClustersRow;
659         fRowCl->SetClass("AliTPCclusterMI");
660         fRowCl->SetArray(1);
661         fRowCl->SetID(fParam->GetIndex(fSector, iRow));
662         fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
663
664         FindClusters();
665
666         fOutput->Fill();
667         delete fRowCl;    
668         nclusters += fNcluster;    
669         delete[] fBins;
670         delete[] fResBins;
671       }
672
673   delete[] splitRows;
674   delete[] splitRowsRes;
675   Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters);
676 }
677
678 void AliTPCclustererMI::FindClusters()
679 {
680   //add virtual charge at the edge   
681   for (Int_t i=0; i<fMaxTime; i++){
682     Float_t amp1 = fBins[i+3*fMaxTime]; 
683     Float_t amp0 =0;
684     if (amp1>0){
685       Float_t amp2 = fBins[i+4*fMaxTime];
686       if (amp2==0) amp2=0.5;
687       Float_t sigma2 = GetSigmaY2(i);           
688       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);   
689       if (gDebug>4) printf("\n%f\n",amp0);
690     }  
691     fBins[i+2*fMaxTime] = Int_t(amp0);    
692     amp0 = 0;
693     amp1 = fBins[(fMaxPad+2)*fMaxTime+i];
694     if (amp1>0){
695       Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime];
696       if (amp2==0) amp2=0.5;
697       Float_t sigma2 = GetSigmaY2(i);           
698       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);           
699       if (gDebug>4) printf("\n%f\n",amp0);
700     }        
701     fBins[(fMaxPad+3)*fMaxTime+i] = Int_t(amp0);           
702   }
703
704 //  memcpy(fResBins,fBins, fMaxBin*2);
705   memcpy(fResBins,fBins, fMaxBin);
706   //
707   fNcluster=0;
708   //first loop - for "gold cluster" 
709   fLoop=1;
710   Int_t *b=&fBins[-1]+2*fMaxTime;
711   Int_t crtime = Int_t((fParam->GetZLength()-AliTPCReconstructor::GetCtgRange()*fRx)/fZWidth-5);
712
713   for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
714     b++;
715     if (*b<8) continue;   //threshold form maxima
716     if (i%fMaxTime<crtime) {
717       Int_t delta = -(i%fMaxTime)+crtime;
718       b+=delta;
719       i+=delta;
720       continue; 
721     }
722      
723     if (!IsMaximum(*b,fMaxTime,b)) continue;
724     AliTPCclusterMI c;
725     Int_t dummy=0;
726     MakeCluster(i, fMaxTime, fBins, dummy,c);
727     //}
728   }
729   //memcpy(fBins,fResBins, fMaxBin*2);
730   //second  loop - for rest cluster 
731   /*        
732   fLoop=2;
733   b=&fResBins[-1]+2*fMaxTime;
734   for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
735     b++;
736     if (*b<25) continue;   // bigger threshold for maxima
737     if (!IsMaximum(*b,fMaxTime,b)) continue;
738     AliTPCclusterMI c;
739     Int_t dummy;
740     MakeCluster(i, fMaxTime, fResBins, dummy,c);
741     //}
742   }
743   */
744 }