Fixes for Coverity warnings
[u/mrichter/AliRoot.git] / ITS / AliITSclustererV2.cxx
1 //-------------------------------------------------------------------------
2 //            Implementation of the ITS clusterer V2 class
3 //
4 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
5 //-------------------------------------------------------------------------
6
7
8 #include "AliRun.h"
9
10 #include "AliITSclustererV2.h"
11 #include "AliITSclusterV2.h"
12 #include "AliITSDetTypeRec.h"
13 #include "AliRawReader.h"
14 #include "AliITSRawStreamSPD.h"
15 #include "AliITSRawStreamSDD.h"
16 #include "AliITSRawStreamSSD.h"
17
18 #include <TFile.h>
19 #include <TTree.h>
20 #include <TClonesArray.h>
21 #include "AliITSgeomTGeo.h"
22 #include "AliITSdigitSPD.h"
23 #include "AliITSdigitSDD.h"
24 #include "AliITSdigitSSD.h"
25 #include "AliMC.h"
26
27 ClassImp(AliITSclustererV2)
28
29 extern AliRun *gAlice;
30
31 AliITSclustererV2::AliITSclustererV2():
32 fNModules(0),
33 fEvent(0),
34 fI(0),
35 fLastSPD1(0),
36 fNySPD(0),
37 fNzSPD(0),
38 fYpitchSPD(0),
39 fZ1pitchSPD(0),
40 fZ2pitchSPD(0),
41 fHwSPD(0),
42 fHlSPD(0),
43 fNySDD(0),
44 fNzSDD(0),
45 fYpitchSDD(0),
46 fZpitchSDD(0),
47 fHwSDD(0),
48 fHlSDD(0),
49 fYoffSDD(0),
50 fLastSSD1(0),
51 fYpitchSSD(0),
52 fHwSSD(0),
53 fHlSSD(0),
54 fTanP(0),
55 fTanN(0){
56    //default constructor
57   for(Int_t i=0;i<260;i++)fYSPD[i]=0.;
58   for(Int_t i=0;i<170;i++)fZSPD[i]=0.;
59   for(Int_t i=0;i<2200;i++){
60     fYshift[i]=0.;
61     fZshift[i]=0.;
62     fNdet[i]=0;
63     fNlayer[i]=0;
64   }
65
66 }
67 AliITSclustererV2::AliITSclustererV2(const Char_t *geom):
68 fNModules(AliITSgeomTGeo::GetNModules()),
69 fEvent(0),
70 fI(0),
71 fLastSPD1(AliITSgeomTGeo::GetModuleIndex(2,1,1)-1),
72 fNySPD(256),
73 fNzSPD(160),
74 fYpitchSPD(0.0050),
75 fZ1pitchSPD(0.0425),
76 fZ2pitchSPD(0.0625),
77 fHwSPD(0.64),
78 fHlSPD(3.48),
79 fNySDD(256),
80 fNzSDD(256),
81 fYpitchSDD(0.01825),
82 fZpitchSDD(0.02940),
83 fHwSDD(3.5085),
84 fHlSDD(3.7632),
85 fYoffSDD(0.0425),
86 fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1),
87 fYpitchSSD(0.0095),
88 fHwSSD(3.65),
89 fHlSSD(2.00),
90 fTanP(0.0275),
91 fTanN(0.0075) {
92   //------------------------------------------------------------
93   // Standard constructor
94   //------------------------------------------------------------
95   if (geom) {
96     AliWarning("\"geom\" is actually a dummy argument !");
97   }
98   Int_t m;
99   for (m=0; m<fNModules; m++) {
100      Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
101      const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(m);
102      fYshift[m] = (tm->Inverse()).GetTranslation()[1];
103      fZshift[m] = (tm->Inverse()).GetTranslation()[2];
104      fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
105      fNlayer[m] = lay-1;
106   }
107
108   //SPD geometry  
109   fYSPD[0]=0.5*fYpitchSPD;
110   for (m=1; m<fNySPD; m++) fYSPD[m]=fYSPD[m-1]+fYpitchSPD; 
111   fZSPD[0]=fZ1pitchSPD;
112   for (m=1; m<fNzSPD; m++) {
113     Double_t dz=fZ1pitchSPD;
114     if (m==31 || m==32 || m==63  || m==64  || m==95 || m==96 || 
115         m==127 || m==128) dz=fZ2pitchSPD; 
116     fZSPD[m]=fZSPD[m-1]+dz;
117   }
118   for (m=0; m<fNzSPD; m++) {
119     Double_t dz=0.5*fZ1pitchSPD;
120     if (m==31 || m==32 || m==63  || m==64  || m==95 || m==96 || 
121         m==127 || m==128) dz=0.5*fZ2pitchSPD; 
122     fZSPD[m]-=dz;
123   }
124
125 }
126
127
128 Int_t AliITSclustererV2::Digits2Clusters(TTree *dTree, TTree *cTree) {
129   //------------------------------------------------------------
130   // This function creates ITS clusters
131   //------------------------------------------------------------
132   Int_t ncl=0;
133
134   if (!dTree) {
135     Error("Digits2Clusters","Can't get the tree with digits !");
136     return 1;
137   }
138
139   TClonesArray *digitsSPD=new TClonesArray("AliITSdigitSPD",3000);
140   dTree->SetBranchAddress("ITSDigitsSPD",&digitsSPD);
141   TClonesArray *digitsSDD=new TClonesArray("AliITSdigitSDD",3000);
142   dTree->SetBranchAddress("ITSDigitsSDD",&digitsSDD);
143   TClonesArray *digitsSSD=new TClonesArray("AliITSdigitSSD",3000);
144   dTree->SetBranchAddress("ITSDigitsSSD",&digitsSSD);
145
146   TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
147   TBranch *branch=cTree->GetBranch("Clusters");
148   if (!branch) cTree->Branch("Clusters",&clusters);
149   else branch->SetAddress(&clusters);
150
151   Int_t mmax=(Int_t)dTree->GetEntries();
152   if (mmax!=fNModules) {
153     Error("Digits2Clusters","Number of entries != number of modules !");
154     return 1;
155   }
156
157   for (fI=0; fI<mmax; fI++) {
158     dTree->GetEvent(fI);
159
160     if     (digitsSPD->GetEntriesFast()!=0) 
161       FindClustersSPD(digitsSPD,clusters);
162     else 
163       if(digitsSDD->GetEntriesFast()!=0) 
164         FindClustersSDD(digitsSDD,clusters);
165       else if(digitsSSD->GetEntriesFast()!=0) 
166         FindClustersSSD(digitsSSD,clusters);
167     
168     ncl+=clusters->GetEntriesFast();
169
170     cTree->Fill();
171
172     digitsSPD->Clear();
173     digitsSDD->Clear();
174     digitsSSD->Clear();
175     clusters->Clear();
176   }
177
178   //cTree->Write();
179
180   delete clusters;
181
182   delete digitsSPD;
183   delete digitsSDD;
184   delete digitsSSD;
185
186   //delete dTree;
187
188   Info("Digits2Clusters","Number of found clusters : %d",ncl);
189
190   return 0;
191 }
192
193 void AliITSclustererV2::Digits2Clusters(AliRawReader* rawReader) {
194   //------------------------------------------------------------
195   // This function creates ITS clusters from raw data
196   //------------------------------------------------------------
197   AliRunLoader* runLoader = AliRunLoader::Instance();
198   if (!runLoader) {
199     Error("Digits2Clusters", "no run loader found");
200     return;
201   }
202   runLoader->LoadKinematics();
203   AliLoader* itsLoader = runLoader->GetLoader("ITSLoader");
204   if (!itsLoader) {
205     Error("Digits2Clusters", "no loader for ITS found");
206     return;
207   }
208   if (!itsLoader->TreeR()) itsLoader->MakeTree("R");
209   TTree* cTree = itsLoader->TreeR();
210
211   TClonesArray *array=new TClonesArray("AliITSclusterV2",1000);
212   cTree->Branch("Clusters",&array);
213   delete array;
214
215   TClonesArray** clusters = new TClonesArray*[fNModules]; 
216   for (Int_t iModule = 0; iModule < fNModules; iModule++) {
217     clusters[iModule] = NULL;
218   }
219   // one TClonesArray per module
220
221   rawReader->Reset();
222   AliITSRawStreamSPD inputSPD(rawReader);
223   FindClustersSPD(&inputSPD, clusters);
224
225   rawReader->Reset();
226   AliITSRawStreamSDD inputSDD(rawReader);
227   FindClustersSDD(&inputSDD, clusters);
228
229   rawReader->Reset();
230   AliITSRawStreamSSD inputSSD(rawReader);
231   FindClustersSSD(&inputSSD, clusters);
232
233   // write all clusters to the tree
234   Int_t nClusters = 0;
235   TClonesArray *emptyArray=new TClonesArray("AliITSclusterV2");
236   for (Int_t iModule = 0; iModule < fNModules; iModule++) {
237     array = clusters[iModule];
238     if (!array) {
239       Error("Digits2Clusters", "data for module %d missing!", iModule);
240       array = emptyArray;
241     }
242     cTree->SetBranchAddress("Clusters", &array);
243     cTree->Fill();
244     nClusters += array->GetEntriesFast();
245   }
246   delete emptyArray;
247
248   itsLoader->WriteRecPoints("OVERWRITE");
249
250   delete[] clusters;
251
252   Info("Digits2Clusters", "total number of found clusters in ITS: %d\n", 
253        nClusters);
254 }
255
256 //**** Fast clusters *******************************
257 #include "TParticle.h"
258
259 //#include "AliITS.h"
260 #include "AliITSmodule.h"
261 #include "AliITSRecPoint.h"
262 #include "AliITSsimulationFastPoints.h"
263 #include "AliITSRecPoint.h"
264
265
266 static void CheckLabels(Int_t lab[3]) {
267   //------------------------------------------------------------
268   // Tries to find mother's labels
269   //------------------------------------------------------------
270   AliRunLoader *rl = AliRunLoader::Instance();
271   TTree *trK=(TTree*)rl->TreeK();
272
273   if(trK){
274     Int_t label=lab[0];
275     if (label>=0) {
276       TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
277       label=-3;
278       while (part->P() < 0.005) {
279         Int_t m=part->GetFirstMother();
280         if (m<0) {
281           Info("CheckLabels","Primary momentum: %f",part->P()); 
282           break;
283         }
284         if (part->GetStatusCode()>0) {
285           Info("CheckLabels","Primary momentum: %f",part->P()); 
286           break;
287         }
288         label=m;
289         part=(TParticle*)gAlice->GetMCApp()->Particle(label);
290       }
291       if(lab[1]<0){
292         lab[1]=label;
293       }
294       else if (lab[2]<0) {
295         lab[2]=label;
296       }
297       else {
298         //      cerr<<"CheckLabels : No empty labels !\n";
299       }
300     }
301   }
302 }
303
304 /*
305 static void CheckLabels(Int_t lab[3]) {
306   //------------------------------------------------------------
307   // Tries to find mother's labels
308   //------------------------------------------------------------
309
310   if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit
311
312   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
313   for (Int_t i=0;i<3;i++){
314     Int_t label = lab[i];
315     if (label>=0 && label<ntracks) {
316       TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
317       if (part->P() < 0.005) {
318         Int_t m=part->GetFirstMother();
319         if (m<0) {      
320           continue;
321         }
322         if (part->GetStatusCode()>0) {
323           continue;
324         }
325         lab[i]=m;       
326       }
327     }    
328   }
329   
330 }
331 */
332 static void CheckLabels2(Int_t lab[10]) {
333   //------------------------------------------------------------
334   // Tries to find mother's labels
335   //------------------------------------------------------------
336   AliRunLoader *rl = AliRunLoader::Instance();
337   TTree *trK=(TTree*)rl->TreeK();
338
339   if(trK){
340     Int_t nlabels =0; 
341     for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
342     if(nlabels == 0) return; // In case of no labels just exit
343
344     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
345
346     for (Int_t i=0;i<10;i++){
347       Int_t label = lab[i];
348       if (label>=0 && label<ntracks) {
349         TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
350         if (part->P() < 0.02) {
351           Int_t m=part->GetFirstMother();
352           if (m<0) {    
353             continue;
354           }
355           if (part->GetStatusCode()>0) {
356             continue;
357           }
358           lab[i]=m;       
359         }
360         else
361           if (part->P() < 0.12 && nlabels>3) {
362             lab[i]=-2;
363             nlabels--;
364           } 
365       }
366       else{
367         if ( (label>ntracks||label <0) && nlabels>3) {
368           lab[i]=-2;
369           nlabels--;
370         } 
371       }
372     }  
373     if (nlabels>3){
374       for (Int_t i=0;i<10;i++){
375         if (nlabels>3){
376           Int_t label = lab[i];
377           if (label>=0 && label<ntracks) {
378             TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
379             if (part->P() < 0.1) {
380               lab[i]=-2;
381               nlabels--;
382             }
383           }
384         }
385       }
386     }
387
388     //compress labels -- if multi-times the same
389     Int_t lab2[10];
390     for (Int_t i=0;i<10;i++) lab2[i]=-2;
391     for (Int_t i=0;i<10  ;i++){
392       if (lab[i]<0) continue;
393       for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
394         if (lab2[j]<0) {
395           lab2[j]= lab[i];
396           break;
397         }
398       }
399     }
400     for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
401   }
402 }
403
404 static void AddLabel(Int_t lab[10], Int_t label) {
405 // add label of the particle in the kine tree which originates this cluster
406 // (used for reconstruction of MC data only, for comparison purposes)
407   AliRunLoader *rl = AliRunLoader::Instance();
408   TTree *trK=(TTree*)rl->TreeK();
409
410   if(trK){
411     if(label<0) return; // In case of no label just exit
412
413     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
414     if (label>ntracks) return;
415     for (Int_t i=0;i<10;i++){
416       //    if (label<0) break;
417       if (lab[i]==label) break;
418       if (lab[i]<0) {
419         lab[i]= label;
420         break;
421       }
422     }
423   }
424 }
425
426 void AliITSclustererV2::RecPoints2Clusters
427 (const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
428   //------------------------------------------------------------
429   // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS 
430   // subdetector indexed by idx 
431   //------------------------------------------------------------
432   TClonesArray &cl=*clusters;
433   Int_t ncl=points->GetEntriesFast();
434   for (Int_t i=0; i<ncl; i++) {
435     AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
436     Float_t lp[5];
437     lp[0]=-(-p->GetDetLocalX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1
438     lp[1]=  -p->GetZ()+fZshift[idx];
439     lp[2]=p->GetSigmaDetLocX2();
440     lp[3]=p->GetSigmaZ2();
441     lp[4]=p->GetQ()*36./23333.;  //electrons -> ADC
442     Int_t lab[4]; 
443     lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
444     lab[3]=fNdet[idx];
445     CheckLabels(lab);
446     Int_t dummy[3]={0,0,0};
447     new (cl[i]) AliITSclusterV2(lab,lp, dummy);
448   }  
449
450
451 //***********************************
452
453
454 void AliITSclustererV2:: 
455 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
456   //------------------------------------------------------------
457   // returns an array of indices of digits belonging to the cluster
458   // (needed when the segmentation is not regular) 
459   //------------------------------------------------------------
460   if (n<200) idx[n++]=bins[k].GetIndex();
461   bins[k].Use();
462
463   if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
464   if (bins[k-1   ].IsNotUsed()) FindCluster(k-1   ,maxz,bins,n,idx);
465   if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
466   if (bins[k+1   ].IsNotUsed()) FindCluster(k+1   ,maxz,bins,n,idx);
467   /*
468   if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
469   if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
470   if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
471   if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
472   */
473 }
474
475 void AliITSclustererV2::
476 FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) {
477   //------------------------------------------------------------
478   // Actual SPD cluster finder
479   //------------------------------------------------------------
480   Int_t kNzBins = fNzSPD + 2;
481   const Int_t kMAXBIN=kNzBins*(fNySPD+2);
482
483   Int_t ndigits=digits->GetEntriesFast();
484   AliBin *bins=new AliBin[kMAXBIN];
485
486   Int_t k;
487   AliITSdigitSPD *d=0;
488   for (k=0; k<ndigits; k++) {
489      d=(AliITSdigitSPD*)digits->UncheckedAt(k);
490      Int_t i=d->GetCoord2()+1;   //y
491      Int_t j=d->GetCoord1()+1;
492      bins[i*kNzBins+j].SetIndex(k);
493      bins[i*kNzBins+j].SetMask(1);
494   }
495    
496   Int_t n=0; TClonesArray &cl=*clusters;
497   for (k=0; k<kMAXBIN; k++) {
498      if (!bins[k].IsNotUsed()) continue;
499      Int_t ni=0, idx[200];
500      FindCluster(k,kNzBins,bins,ni,idx);
501      if (ni==200) {
502         Info("FindClustersSPD","Too big cluster !"); 
503         continue;
504      }
505      Int_t milab[10];
506      for (Int_t ilab=0;ilab<10;ilab++){
507        milab[ilab]=-2;
508      }
509
510      d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]);
511      Int_t ymin=d->GetCoord2(),ymax=ymin;
512      Int_t zmin=d->GetCoord1(),zmax=zmin;
513
514      for (Int_t l=0; l<ni; l++) {
515         d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]);
516
517         if (ymin > d->GetCoord2()) ymin=d->GetCoord2();
518         if (ymax < d->GetCoord2()) ymax=d->GetCoord2();
519         if (zmin > d->GetCoord1()) zmin=d->GetCoord1();
520         if (zmax < d->GetCoord1()) zmax=d->GetCoord1();
521         // MI addition - find all labels in cluster
522         for (Int_t dlab=0;dlab<10;dlab++){
523           Int_t digitlab = (d->GetTracks())[dlab];
524           if (digitlab<0) continue;
525           AddLabel(milab,digitlab);       
526         }
527         if (milab[9]>0) CheckLabels2(milab);
528      }
529      CheckLabels2(milab);
530      //
531      //Int_t idy = (fNlayer[fI]==0)? 2:3; 
532      //for (Int_t iz=zmin; iz<=zmax;iz+=2)
533      //Int_t idy = (ymax-ymin)/4.; // max 2 clusters
534      Int_t idy = 0; // max 2 clusters
535      if (fNlayer[fI]==0 &&idy<3) idy=3;
536      if (fNlayer[fI]==1 &&idy<4) idy=4; 
537      Int_t idz =3;
538      for (Int_t iz=zmin; iz<=zmax;iz+=idz)
539        for (Int_t iy=ymin; iy<=ymax;iy+=idy){
540          //
541          Int_t nodigits =0;
542          Float_t y=0.,z=0.,q=0.;         
543          for (Int_t l=0; l<ni; l++) {
544            d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]);
545            if (zmax-zmin>=idz || ymax-ymin>=idy){
546              if (TMath::Abs( d->GetCoord2()-iy)>0.75*idy) continue;
547              if (TMath::Abs( d->GetCoord1()-iz)>0.75*idz) continue;
548            }
549            nodigits++;
550            Float_t qq=d->GetSignal();
551            y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq;   
552           
553          }     
554          if (nodigits==0) continue;
555          y/=q; z/=q;
556          y-=fHwSPD; z-=fHlSPD;
557          
558          Float_t lp[5];
559          lp[0]=-(-y+fYshift[fI]); if (fI<=fLastSPD1) lp[0]=-lp[0];
560          lp[1]=  -z+fZshift[fI];
561          // Float_t factor=TMath::Max(double(ni-3.),1.5);
562          Float_t factory=TMath::Max(ymax-ymin,1);
563          Float_t factorz=TMath::Max(zmax-zmin,1);
564          factory*= factory;
565          factorz*= factorz;     
566          //lp[2]= (fYpitchSPD*fYpitchSPD/12.)*factory;
567          //lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.)*factorz;
568          lp[2]= (fYpitchSPD*fYpitchSPD/12.);
569          lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.);
570          //lp[4]= q;
571          lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1);
572          
573          milab[3]=fNdet[fI];
574          Int_t info[3] = {ymax-ymin+1,zmax-zmin+1,fNlayer[fI]};
575          new (cl[n]) AliITSclusterV2(milab,lp,info); n++;        
576        }
577   }
578   
579   delete [] bins;
580 }
581
582 void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, 
583                                         TClonesArray** clusters) 
584 {
585   //------------------------------------------------------------
586   // Actual SPD cluster finder for raw data
587   //------------------------------------------------------------
588
589   Int_t nClustersSPD = 0;
590   Int_t kNzBins = fNzSPD + 2;
591   Int_t kNyBins = fNySPD + 2;
592   Int_t kMaxBin = kNzBins * kNyBins;
593   AliBin *binsSPD = new AliBin[kMaxBin];
594   AliBin *binsSPDInit = new AliBin[kMaxBin];
595   AliBin *bins = NULL;
596
597   // read raw data input stream
598   while (kTRUE) {
599     Bool_t next = input->Next();
600     if (!next || input->IsNewModule()) {
601       Int_t iModule = input->GetPrevModuleID();
602
603       // when all data from a module was read, search for clusters
604       if (bins) { 
605         clusters[iModule] = new TClonesArray("AliITSclusterV2");
606         Int_t nClusters = 0;
607
608         for (Int_t iBin = 0; iBin < kMaxBin; iBin++) {
609           if (bins[iBin].IsUsed()) continue;
610           Int_t nBins = 0;
611           Int_t idxBins[200];
612           FindCluster(iBin, kNzBins, bins, nBins, idxBins);
613           if (nBins == 200) {
614             Error("FindClustersSPD", "SPD: Too big cluster !\n"); 
615             continue;
616           }
617
618           Int_t label[4]; 
619           label[0] = -2;
620           label[1] = -2;
621           label[2] = -2;
622 //        label[3] = iModule;
623           label[3] = fNdet[iModule];
624
625           Int_t ymin = (idxBins[0] / kNzBins) - 1;
626           Int_t ymax = ymin;
627           Int_t zmin = (idxBins[0] % kNzBins) - 1;
628           Int_t zmax = zmin;
629           for (Int_t idx = 0; idx < nBins; idx++) {
630             Int_t iy = (idxBins[idx] / kNzBins) - 1;
631             Int_t iz = (idxBins[idx] % kNzBins) - 1;
632             if (ymin > iy) ymin = iy;
633             if (ymax < iy) ymax = iy;
634             if (zmin > iz) zmin = iz;
635             if (zmax < iz) zmax = iz;
636           }
637
638           Int_t idy = 0; // max 2 clusters
639           if ((iModule <= fLastSPD1) &&idy<3) idy=3;
640           if ((iModule > fLastSPD1) &&idy<4) idy=4; 
641           Int_t idz =3;
642           for (Int_t iiz=zmin; iiz<=zmax;iiz+=idz)
643             for (Int_t iiy=ymin; iiy<=ymax;iiy+=idy){
644               //
645               Int_t ndigits =0;
646               Float_t y=0.,z=0.,q=0.;    
647               for (Int_t idx = 0; idx < nBins; idx++) {
648                 Int_t iy = (idxBins[idx] / kNzBins) - 1;
649                 Int_t iz = (idxBins[idx] % kNzBins) - 1;
650                 if (zmax-zmin>=idz || ymax-ymin>=idy){
651                   if (TMath::Abs(iy-iiy)>0.75*idy) continue;
652                   if (TMath::Abs(iz-iiz)>0.75*idz) continue;
653                 }
654                 ndigits++;
655                 Float_t qBin = bins[idxBins[idx]].GetQ();
656                 y += qBin * fYSPD[iy]; 
657                 z += qBin * fZSPD[iz]; 
658                 q += qBin;   
659               }
660               if (ndigits==0) continue;
661               y /= q; 
662               z /= q;
663               y -= fHwSPD; 
664               z -= fHlSPD;
665
666               Float_t hit[5];  // y, z, sigma(y)^2, sigma(z)^2, charge
667               hit[0] = -(-y+fYshift[iModule]); 
668               if (iModule <= fLastSPD1) hit[0] = -hit[0];
669               hit[1] = -z+fZshift[iModule];
670               hit[2] = fYpitchSPD*fYpitchSPD/12.;
671               hit[3] = fZ1pitchSPD*fZ1pitchSPD/12.;
672               //          hit[4] = q;
673               hit[4] = (zmax-zmin+1)*100 + (ymax-ymin+1);
674               //          CheckLabels(label);
675               Int_t info[3]={ymax-ymin+1,zmax-zmin+1,fNlayer[iModule]};
676               new (clusters[iModule]->AddrAt(nClusters)) 
677                 AliITSclusterV2(label, hit,info); 
678               nClusters++;
679             }
680         }
681
682         nClustersSPD += nClusters;
683         bins = NULL;
684       }
685
686       if (!next) break;
687       bins = binsSPD;
688       memcpy(binsSPD,binsSPDInit,sizeof(AliBin)*kMaxBin);
689     }
690
691     // fill the current digit into the bins array
692     if(bins){
693       Int_t index = (input->GetCoord2()+1) * kNzBins + (input->GetCoord1()+1);
694       bins[index].SetIndex(index);
695       bins[index].SetMask(1);
696       bins[index].SetQ(1);
697     }
698   }
699
700   delete [] binsSPDInit;
701   delete [] binsSPD;
702
703   Info("FindClustersSPD", "found clusters in ITS SPD: %d", nClustersSPD);
704 }
705
706
707 Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
708   //------------------------------------------------------------
709   //is this a local maximum ?
710   //------------------------------------------------------------
711   UShort_t q=bins[k].GetQ();
712   if (q==1023) return kFALSE;
713   if (bins[k-max].GetQ() > q) return kFALSE;
714   if (bins[k-1  ].GetQ() > q) return kFALSE; 
715   if (bins[k+max].GetQ() > q) return kFALSE; 
716   if (bins[k+1  ].GetQ() > q) return kFALSE; 
717   if (bins[k-max-1].GetQ() > q) return kFALSE;
718   if (bins[k+max-1].GetQ() > q) return kFALSE; 
719   if (bins[k+max+1].GetQ() > q) return kFALSE; 
720   if (bins[k-max+1].GetQ() > q) return kFALSE;
721   return kTRUE; 
722 }
723
724 void AliITSclustererV2::
725 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
726   //------------------------------------------------------------
727   //find local maxima
728   //------------------------------------------------------------
729   if (n<31)
730   if (IsMaximum(k,max,b)) {
731     idx[n]=k; msk[n]=(2<<n);
732     n++;
733   }
734   b[k].SetMask(0);
735   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
736   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
737   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
738   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
739 }
740
741 void AliITSclustererV2::
742 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
743   //------------------------------------------------------------
744   //mark this peak
745   //------------------------------------------------------------
746   UShort_t q=bins[k].GetQ();
747
748   bins[k].SetMask(bins[k].GetMask()|m); 
749
750   if (bins[k-max].GetQ() <= q)
751      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
752   if (bins[k-1  ].GetQ() <= q)
753      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
754   if (bins[k+max].GetQ() <= q)
755      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
756   if (bins[k+1  ].GetQ() <= q)
757      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
758 }
759
760 void AliITSclustererV2::
761 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) {
762   //------------------------------------------------------------
763   //make cluster using digits of this peak
764   //------------------------------------------------------------
765   Float_t q=(Float_t)bins[k].GetQ();
766   Int_t i=k/max, j=k-i*max;
767
768   c.SetQ(c.GetQ()+q);
769   c.SetY(c.GetY()+i*q); 
770   c.SetZ(c.GetZ()+j*q); 
771   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
772   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
773
774   bins[k].SetMask(0xFFFFFFFE);
775   
776   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
777   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
778   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
779   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
780 }
781
782 void AliITSclustererV2::
783 FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, 
784                 const TClonesArray *digits, TClonesArray *clusters) {
785   //------------------------------------------------------------
786   // Actual SDD cluster finder
787   //------------------------------------------------------------
788   Int_t ncl=0; TClonesArray &cl=*clusters;
789   for (Int_t s=0; s<2; s++)
790     for (Int_t i=0; i<nMaxBin; i++) {
791       if (bins[s][i].IsUsed()) continue;
792       Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
793       FindPeaks(i, nzBins, bins[s], idx, msk, npeaks);
794
795       if (npeaks>30) continue;
796       if (npeaks==0) continue;
797
798       Int_t k,l;
799       for (k=0; k<npeaks-1; k++){//mark adjacent peaks
800         if (idx[k] < 0) continue; //this peak is already removed
801         for (l=k+1; l<npeaks; l++) {
802            if (idx[l] < 0) continue; //this peak is already removed
803            Int_t ki=idx[k]/nzBins, kj=idx[k] - ki*nzBins;
804            Int_t li=idx[l]/nzBins, lj=idx[l] - li*nzBins;
805            Int_t di=TMath::Abs(ki - li);
806            Int_t dj=TMath::Abs(kj - lj);
807            if (di>1 || dj>1) continue;
808            if (bins[s][idx[k]].GetQ() > bins[s][idx[l]].GetQ()) {
809               msk[l]=msk[k];
810               idx[l]*=-1;
811            } else {
812               msk[k]=msk[l];
813               idx[k]*=-1;
814               break;
815            } 
816         }
817       }
818
819       for (k=0; k<npeaks; k++) {
820         MarkPeak(TMath::Abs(idx[k]), nzBins, bins[s], msk[k]);
821       }
822         
823       for (k=0; k<npeaks; k++) {
824          if (idx[k] < 0) continue; //removed peak
825          AliITSclusterV2 c;
826          MakeCluster(idx[k], nzBins, bins[s], msk[k], c);
827          //mi change
828          Int_t milab[10];
829          for (Int_t ilab=0;ilab<10;ilab++){
830            milab[ilab]=-2;
831          }
832          Int_t maxi=0,mini=0,maxj=0,minj=0;
833          //AliBin *bmax=&bins[s][idx[k]];
834          //Float_t max = TMath::Max(TMath::Abs(bmax->GetQ())/5.,3.);
835          Float_t max=3;
836          for (Int_t di=-2; di<=2;di++)
837            for (Int_t dj=-3;dj<=3;dj++){
838              Int_t index = idx[k]+di+dj*nzBins;
839              if (index<0) continue;
840              if (index>=nMaxBin) continue;
841              AliBin *b=&bins[s][index];
842              if (TMath::Abs(b->GetQ())>max){
843                if (di>maxi) maxi=di;
844                if (di<mini) mini=di;
845                if (dj>maxj) maxj=dj;
846                if (dj<minj) minj=dj;
847                //
848                if(digits) {
849                  if (TMath::Abs(di)<2&&TMath::Abs(dj)<2){
850                    AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
851                    for (Int_t itrack=0;itrack<10;itrack++){
852                      Int_t track = (d->GetTracks())[itrack];
853                      if (track>=0) {
854                        AddLabel(milab, track); 
855                      }
856                    }
857                  }
858                }
859              }
860            }
861          
862          /* 
863             Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
864             Float_t w=par->GetPadPitchWidth(sec);
865             c.SetSigmaY2(s2);
866             if (s2 != 0.) {
867             c.SetSigmaY2(c.GetSigmaY2()*0.108);
868             if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
869             }    
870             s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
871             w=par->GetZWidth();
872             c.SetSigmaZ2(s2);
873             
874             if (s2 != 0.) {
875             c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
876             if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
877             }
878          */
879
880          Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ();
881          y/=q; z/=q;
882          //
883          //Float_t s2 = c.GetSigmaY2()/c.GetQ() - y*y;
884          // c.SetSigmaY2(s2);
885          //s2 = c.GetSigmaZ2()/c.GetQ() - z*z;
886          //c.SetSigmaZ2(s2);
887          //
888          y=(y-0.5)*fYpitchSDD;
889          y-=fHwSDD;
890          y-=fYoffSDD;  //delay ?
891          if (s) y=-y;
892
893          z=(z-0.5)*fZpitchSDD;
894          z-=fHlSDD;
895
896          y=-(-y+fYshift[fI]);
897          z=  -z+fZshift[fI];
898
899          q/=12.7;  //this WAS consistent with SSD. To be reassessed 
900                    // 23-MAR-2007
901          Float_t hit[5] = {y, z, 0.0030*0.0030, 0.0020*0.0020, q};
902          Int_t  info[3] = {maxj-minj+1, maxi-mini+1, fNlayer[fI]};
903
904          if (c.GetQ() < 20.) continue; //noise cluster
905          
906          if (digits) {    
907            //      AliBin *b=&bins[s][idx[k]];
908            //      AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
909            {
910              //Int_t lab[3];
911              //lab[0]=(d->GetTracks())[0];
912              //lab[1]=(d->GetTracks())[1];
913              //lab[2]=(d->GetTracks())[2];
914              //CheckLabels(lab);
915              CheckLabels2(milab); 
916            }
917          }
918          milab[3]=fNdet[fI];
919
920          AliITSclusterV2 *cc = new (cl[ncl]) AliITSclusterV2(milab,hit,info);
921          cc->SetType(npeaks); 
922          ncl++;
923       }
924     }
925 }
926
927 void AliITSclustererV2::
928 FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) {
929   //------------------------------------------------------------
930   // Actual SDD cluster finder
931   //------------------------------------------------------------
932   Int_t kNzBins = fNzSDD + 2;
933   const Int_t kMAXBIN=kNzBins*(fNySDD+2);
934
935   AliBin *bins[2];
936   bins[0]=new AliBin[kMAXBIN];
937   bins[1]=new AliBin[kMAXBIN];
938
939   AliITSdigitSDD *d=0;
940   Int_t i, ndigits=digits->GetEntriesFast();
941   for (i=0; i<ndigits; i++) {
942      d=(AliITSdigitSDD*)digits->UncheckedAt(i);
943      Int_t y=d->GetCoord2()+1;   //y
944      Int_t z=d->GetCoord1()+1;   //z
945      Int_t q=d->GetSignal();
946      if (q<3) continue;
947
948      if (z <= fNzSDD) {
949        bins[0][y*kNzBins+z].SetQ(q);
950        bins[0][y*kNzBins+z].SetMask(1);
951        bins[0][y*kNzBins+z].SetIndex(i);
952      } else {
953        z-=fNzSDD; 
954        bins[1][y*kNzBins+z].SetQ(q);
955        bins[1][y*kNzBins+z].SetMask(1);
956        bins[1][y*kNzBins+z].SetIndex(i);
957      }
958   }
959   
960   FindClustersSDD(bins, kMAXBIN, kNzBins, digits, clusters);
961
962   delete[] bins[0];
963   delete[] bins[1];
964
965 }
966
967 void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, 
968                                         TClonesArray** clusters) 
969 {
970   //------------------------------------------------------------
971   // Actual SDD cluster finder for raw data
972   //------------------------------------------------------------
973   Int_t nClustersSDD = 0;
974   Int_t kNzBins = fNzSDD + 2;
975   Int_t kMaxBin = kNzBins * (fNySDD+2);
976   AliBin *binsSDDInit = new AliBin[kMaxBin];
977   AliBin *binsSDD1 = new AliBin[kMaxBin];
978   AliBin *binsSDD2 = new AliBin[kMaxBin];
979   AliBin *bins[2] = {NULL, NULL};
980
981   // read raw data input stream
982   while (kTRUE) {
983     Bool_t next = input->Next();
984     if (!next || input->IsNewModule()) {
985       Int_t iModule = input->GetPrevModuleID();
986
987       // when all data from a module was read, search for clusters
988       if (bins[0]) { 
989         clusters[iModule] = new TClonesArray("AliITSclusterV2");
990         fI = iModule;
991         FindClustersSDD(bins, kMaxBin, kNzBins, NULL, clusters[iModule]);
992         Int_t nClusters = clusters[iModule]->GetEntriesFast();
993         nClustersSDD += nClusters;
994         bins[0] = bins[1] = NULL;
995       }
996
997       if (!next) break;
998       bins[0]=binsSDD1;
999       bins[1]=binsSDD2;
1000       memcpy(binsSDD1,binsSDDInit,sizeof(AliBin)*kMaxBin);
1001       memcpy(binsSDD2,binsSDDInit,sizeof(AliBin)*kMaxBin);
1002     }
1003
1004     // fill the current digit into the bins array
1005     if(input->GetSignal()>=3) {
1006       Int_t iz = input->GetCoord1()+1;
1007       Int_t side = ((iz <= fNzSDD) ? 0 : 1);
1008       iz -= side*fNzSDD;
1009       Int_t index = (input->GetCoord2()+1) * kNzBins + iz;
1010       bins[side][index].SetQ(input->GetSignal());
1011       bins[side][index].SetMask(1);
1012       bins[side][index].SetIndex(index);
1013     }
1014   }
1015
1016   delete[] binsSDD1;
1017   delete[] binsSDD2;
1018   delete[] binsSDDInit;
1019
1020   Info("FindClustersSDD", "found clusters in ITS SDD: %d", nClustersSDD);
1021 }
1022
1023
1024
1025 void AliITSclustererV2::
1026 FindClustersSSD(Ali1Dcluster* neg, Int_t nn, 
1027                 Ali1Dcluster* pos, Int_t np,
1028                 TClonesArray *clusters) {
1029   //------------------------------------------------------------
1030   // Actual SSD cluster finder
1031   //------------------------------------------------------------
1032   TClonesArray &cl=*clusters;
1033   //
1034   Float_t tanp=fTanP, tann=fTanN;
1035   if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;}
1036   Int_t idet=fNdet[fI];
1037   Int_t ncl=0;
1038   //
1039   Int_t negativepair[30000];
1040   Int_t cnegative[3000];  
1041   Int_t cused1[3000];
1042   Int_t positivepair[30000];
1043   Int_t cpositive[3000];
1044   Int_t cused2[3000];
1045   for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
1046   for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
1047   for (Int_t i=0;i<30000;i++){negativepair[i]=positivepair[i]=0;}
1048   static Short_t pairs[1000][1000];
1049   memset(pairs,0,sizeof(Short_t)*1000000);
1050 //   Short_t ** pairs = new Short_t*[1000];
1051 //   for (Int_t i=0; i<1000; i++) {
1052 //     pairs[i] = new Short_t[1000];
1053 //     memset(pairs[i],0,sizeof(Short_t)*1000);
1054 //   }  
1055   //
1056   // find available pairs
1057   //
1058   for (Int_t i=0; i<np; i++) {
1059     Float_t yp=pos[i].GetY()*fYpitchSSD; 
1060     if (pos[i].GetQ()<3) continue;
1061     for (Int_t j=0; j<nn; j++) {
1062       if (neg[j].GetQ()<3) continue;
1063       Float_t yn=neg[j].GetY()*fYpitchSSD;
1064       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1065       Float_t yt=yn + tann*zt;
1066       zt-=fHlSSD; yt-=fHwSSD;
1067       if (TMath::Abs(yt)<fHwSSD+0.01)
1068       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1069         negativepair[i*10+cnegative[i]] =j;  //index
1070         positivepair[j*10+cpositive[j]] =i;
1071         cnegative[i]++;  //counters
1072         cpositive[j]++; 
1073         pairs[i][j]=100;
1074       }
1075     }
1076   }
1077   //
1078   for (Int_t i=0; i<np; i++) {
1079     Float_t yp=pos[i].GetY()*fYpitchSSD; 
1080     if (pos[i].GetQ()<3) continue;
1081     for (Int_t j=0; j<nn; j++) {
1082       if (neg[j].GetQ()<3) continue;
1083       if (cpositive[j]&&cnegative[i]) continue;
1084       Float_t yn=neg[j].GetY()*fYpitchSSD;
1085       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1086       Float_t yt=yn + tann*zt;
1087       zt-=fHlSSD; yt-=fHwSSD;
1088       if (TMath::Abs(yt)<fHwSSD+0.1)
1089       if (TMath::Abs(zt)<fHlSSD+0.15) {
1090         if (cnegative[i]==0) pos[i].SetNd(100);  // not available pair
1091         if (cpositive[j]==0) neg[j].SetNd(100);  // not available pair
1092         negativepair[i*10+cnegative[i]] =j;  //index
1093         positivepair[j*10+cpositive[j]] =i;
1094         cnegative[i]++;  //counters
1095         cpositive[j]++; 
1096         pairs[i][j]=100;
1097       }
1098     }
1099   }
1100   //
1101   Float_t lp[5];
1102   Int_t milab[10];
1103   Double_t ratio;
1104   
1105   //
1106   // sign gold tracks
1107   //
1108   for (Int_t ip=0;ip<np;ip++){
1109     Float_t ybest=1000,zbest=1000,qbest=0;
1110     //
1111     // select gold clusters
1112     if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){ 
1113       Float_t yp=pos[ip].GetY()*fYpitchSSD; 
1114       Int_t j = negativepair[10*ip];      
1115       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1116       //
1117       Float_t yn=neg[j].GetY()*fYpitchSSD;
1118       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1119       Float_t yt=yn + tann*zt;
1120       zt-=fHlSSD; yt-=fHwSSD;
1121       ybest=yt; zbest=zt; 
1122       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1123       lp[0]=-(-ybest+fYshift[fI]);
1124       lp[1]=  -zbest+fZshift[fI];
1125       lp[2]=0.0025*0.0025;  //SigmaY2
1126       lp[3]=0.110*0.110;  //SigmaZ2
1127       
1128       lp[4]=qbest;        //Q
1129       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1130       for (Int_t ilab=0;ilab<3;ilab++){
1131         milab[ilab] = pos[ip].GetLabel(ilab);
1132         milab[ilab+3] = neg[j].GetLabel(ilab);
1133       }
1134       //
1135       CheckLabels2(milab);
1136       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1137       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]};
1138       AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1139       ncl++;
1140       cl2->SetChargeRatio(ratio);       
1141       cl2->SetType(1);
1142       pairs[ip][j]=1;
1143       if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1144         cl2->SetType(2);
1145         pairs[ip][j]=2;
1146       }
1147       cused1[ip]++;
1148       cused2[j]++;
1149     }
1150   }
1151     
1152   for (Int_t ip=0;ip<np;ip++){
1153     Float_t ybest=1000,zbest=1000,qbest=0;
1154     //
1155     //
1156     // select "silber" cluster
1157     if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
1158       Int_t in  = negativepair[10*ip];
1159       Int_t ip2 = positivepair[10*in];
1160       if (ip2==ip) ip2 =  positivepair[10*in+1];
1161       Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
1162       if (TMath::Abs(pcharge-neg[in].GetQ())<10){
1163         //
1164         // add first pair
1165         if (pairs[ip][in]==100){  //
1166           Float_t yp=pos[ip].GetY()*fYpitchSSD; 
1167           Float_t yn=neg[in].GetY()*fYpitchSSD;
1168           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1169           Float_t yt=yn + tann*zt;
1170           zt-=fHlSSD; yt-=fHwSSD;
1171           ybest =yt;  zbest=zt; 
1172           qbest =pos[ip].GetQ();
1173           lp[0]=-(-ybest+fYshift[fI]);
1174           lp[1]=  -zbest+fZshift[fI];
1175           lp[2]=0.0025*0.0025;  //SigmaY2
1176           lp[3]=0.110*0.110;  //SigmaZ2
1177           
1178           lp[4]=qbest;        //Q
1179           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1180           for (Int_t ilab=0;ilab<3;ilab++){
1181             milab[ilab] = pos[ip].GetLabel(ilab);
1182             milab[ilab+3] = neg[in].GetLabel(ilab);
1183           }
1184           //
1185           CheckLabels2(milab);
1186           ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
1187           milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
1188           Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fI]};
1189           AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1190           ncl++;
1191           cl2->SetChargeRatio(ratio);           
1192           cl2->SetType(5);
1193           pairs[ip][in] = 5;
1194           if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1195             cl2->SetType(6);
1196             pairs[ip][in] = 6;
1197           }
1198         }
1199         //
1200         // add second pair
1201         
1202         //      if (!(cused1[ip2] || cused2[in])){  //
1203         if (pairs[ip2][in]==100){
1204           Float_t yp=pos[ip2].GetY()*fYpitchSSD;
1205           Float_t yn=neg[in].GetY()*fYpitchSSD;
1206           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1207           Float_t yt=yn + tann*zt;
1208           zt-=fHlSSD; yt-=fHwSSD;
1209           ybest =yt;  zbest=zt; 
1210           qbest =pos[ip2].GetQ();
1211           lp[0]=-(-ybest+fYshift[fI]);
1212           lp[1]=  -zbest+fZshift[fI];
1213           lp[2]=0.0025*0.0025;  //SigmaY2
1214           lp[3]=0.110*0.110;  //SigmaZ2
1215           
1216           lp[4]=qbest;        //Q
1217           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1218           for (Int_t ilab=0;ilab<3;ilab++){
1219             milab[ilab] = pos[ip2].GetLabel(ilab);
1220             milab[ilab+3] = neg[in].GetLabel(ilab);
1221           }
1222           //
1223           CheckLabels2(milab);
1224           ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
1225           milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
1226           Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fI]};
1227           AliITSclusterV2 *cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1228           ncl++;
1229           cl2->SetChargeRatio(ratio);           
1230           cl2->SetType(5);
1231           pairs[ip2][in] =5;
1232           if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1233             cl2->SetType(6);
1234             pairs[ip2][in] =6;
1235           }
1236         }       
1237         cused1[ip]++;
1238         cused1[ip2]++;
1239         cused2[in]++;
1240       }
1241     }    
1242   }
1243   
1244   //  
1245   for (Int_t jn=0;jn<nn;jn++){
1246     if (cused2[jn]) continue;
1247     Float_t ybest=1000,zbest=1000,qbest=0;
1248     // select "silber" cluster
1249     if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
1250       Int_t ip  = positivepair[10*jn];
1251       Int_t jn2 = negativepair[10*ip];
1252       if (jn2==jn) jn2 =  negativepair[10*ip+1];
1253       Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
1254       //
1255       if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
1256         //
1257         // add first pair
1258         //      if (!(cused1[ip]||cused2[jn])){
1259         if (pairs[ip][jn]==100){
1260           Float_t yn=neg[jn].GetY()*fYpitchSSD; 
1261           Float_t yp=pos[ip].GetY()*fYpitchSSD;
1262           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1263           Float_t yt=yn + tann*zt;
1264           zt-=fHlSSD; yt-=fHwSSD;
1265           ybest =yt;  zbest=zt; 
1266           qbest =neg[jn].GetQ();
1267           lp[0]=-(-ybest+fYshift[fI]);
1268           lp[1]=  -zbest+fZshift[fI];
1269           lp[2]=0.0025*0.0025;  //SigmaY2
1270           lp[3]=0.110*0.110;  //SigmaZ2
1271           
1272           lp[4]=qbest;        //Q
1273           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1274           for (Int_t ilab=0;ilab<3;ilab++){
1275             milab[ilab] = pos[ip].GetLabel(ilab);
1276             milab[ilab+3] = neg[jn].GetLabel(ilab);
1277           }
1278           //
1279           CheckLabels2(milab);
1280           ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1281           milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1282           Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fI]};
1283           AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1284           ncl++;
1285           cl2->SetChargeRatio(ratio);           
1286           cl2->SetType(7);
1287           pairs[ip][jn] =7;
1288           if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1289             cl2->SetType(8);
1290             pairs[ip][jn]=8;
1291           }
1292         }
1293         //
1294         // add second pair
1295         //      if (!(cused1[ip]||cused2[jn2])){
1296         if (pairs[ip][jn2]==100){
1297           Float_t yn=neg[jn2].GetY()*fYpitchSSD; 
1298           Double_t yp=pos[ip].GetY()*fYpitchSSD; 
1299           Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1300           Double_t yt=yn + tann*zt;
1301           zt-=fHlSSD; yt-=fHwSSD;
1302           ybest =yt;  zbest=zt; 
1303           qbest =neg[jn2].GetQ();
1304           lp[0]=-(-ybest+fYshift[fI]);
1305           lp[1]=  -zbest+fZshift[fI];
1306           lp[2]=0.0025*0.0025;  //SigmaY2
1307           lp[3]=0.110*0.110;  //SigmaZ2
1308           
1309           lp[4]=qbest;        //Q
1310           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1311           for (Int_t ilab=0;ilab<3;ilab++){
1312             milab[ilab] = pos[ip].GetLabel(ilab);
1313             milab[ilab+3] = neg[jn2].GetLabel(ilab);
1314           }
1315           //
1316           CheckLabels2(milab);
1317           ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1318           milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1319           Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fI]};
1320           AliITSclusterV2* cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1321           ncl++;
1322           cl2->SetChargeRatio(ratio);           
1323           pairs[ip][jn2]=7;
1324           cl2->SetType(7);
1325           if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1326             cl2->SetType(8);
1327             pairs[ip][jn2]=8;
1328           }
1329         }
1330         cused1[ip]++;
1331         cused2[jn]++;
1332         cused2[jn2]++;
1333       }
1334     }    
1335   }
1336   
1337   for (Int_t ip=0;ip<np;ip++){
1338     Float_t ybest=1000,zbest=1000,qbest=0;
1339     //
1340     // 2x2 clusters
1341     //
1342     if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){ 
1343       Float_t minchargediff =4.;
1344       Int_t j=-1;
1345       for (Int_t di=0;di<cnegative[ip];di++){
1346         Int_t   jc = negativepair[ip*10+di];
1347         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1348         if (TMath::Abs(chargedif)<minchargediff){
1349           j =jc;
1350           minchargediff = TMath::Abs(chargedif);
1351         }
1352       }
1353       if (j<0) continue;  // not proper cluster      
1354       Int_t count =0;
1355       for (Int_t di=0;di<cnegative[ip];di++){
1356         Int_t   jc = negativepair[ip*10+di];
1357         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1358         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1359       }
1360       if (count>1) continue;  // more than one "proper" cluster for positive
1361       //
1362       count =0;
1363       for (Int_t dj=0;dj<cpositive[j];dj++){
1364         Int_t   ic  = positivepair[j*10+dj];
1365         Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1366         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1367       }
1368       if (count>1) continue;  // more than one "proper" cluster for negative
1369       
1370       Int_t jp = 0;
1371       
1372       count =0;
1373       for (Int_t dj=0;dj<cnegative[jp];dj++){
1374         Int_t   ic = positivepair[jp*10+dj];
1375         Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1376         if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1377       }
1378       if (count>1) continue;   
1379       if (pairs[ip][j]<100) continue;
1380       //
1381       //almost gold clusters
1382       Float_t yp=pos[ip].GetY()*fYpitchSSD; 
1383       Float_t yn=neg[j].GetY()*fYpitchSSD;
1384       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1385       Float_t yt=yn + tann*zt;
1386       zt-=fHlSSD; yt-=fHwSSD;
1387       ybest=yt; zbest=zt; 
1388       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1389       lp[0]=-(-ybest+fYshift[fI]);
1390       lp[1]=  -zbest+fZshift[fI];
1391       lp[2]=0.0025*0.0025;  //SigmaY2
1392       lp[3]=0.110*0.110;  //SigmaZ2     
1393       lp[4]=qbest;        //Q
1394       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1395       for (Int_t ilab=0;ilab<3;ilab++){
1396         milab[ilab] = pos[ip].GetLabel(ilab);
1397         milab[ilab+3] = neg[j].GetLabel(ilab);
1398       }
1399       //
1400       CheckLabels2(milab);
1401       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1402       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1403       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]};
1404       AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1405       ncl++;
1406       cl2->SetChargeRatio(ratio);       
1407       cl2->SetType(10);
1408       pairs[ip][j]=10;
1409       if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1410         cl2->SetType(11);
1411         pairs[ip][j]=11;
1412       }
1413       cused1[ip]++;
1414       cused2[j]++;      
1415     }
1416
1417   }
1418   
1419   //  
1420   for (Int_t i=0; i<np; i++) {
1421     Float_t ybest=1000,zbest=1000,qbest=0;
1422     Float_t yp=pos[i].GetY()*fYpitchSSD; 
1423     if (pos[i].GetQ()<3) continue;
1424     for (Int_t j=0; j<nn; j++) {
1425     //    for (Int_t di = 0;di<cpositive[i];di++){
1426     //  Int_t j = negativepair[10*i+di];
1427       if (neg[j].GetQ()<3) continue;
1428       if (cused2[j]||cused1[i]) continue;      
1429       if (pairs[i][j]>0 &&pairs[i][j]<100) continue;
1430       ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());      
1431       Float_t yn=neg[j].GetY()*fYpitchSSD;
1432       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1433       Float_t yt=yn + tann*zt;
1434       zt-=fHlSSD; yt-=fHwSSD;
1435       if (TMath::Abs(yt)<fHwSSD+0.01)
1436       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1437         ybest=yt; zbest=zt; 
1438         qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1439         lp[0]=-(-ybest+fYshift[fI]);
1440         lp[1]=  -zbest+fZshift[fI];
1441         lp[2]=0.0025*0.0025;  //SigmaY2
1442         lp[3]=0.110*0.110;  //SigmaZ2
1443
1444         lp[4]=qbest;        //Q
1445         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1446         for (Int_t ilab=0;ilab<3;ilab++){
1447           milab[ilab] = pos[i].GetLabel(ilab);
1448           milab[ilab+3] = neg[j].GetLabel(ilab);
1449         }
1450         //
1451         CheckLabels2(milab);
1452         milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1453         Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fI]};
1454         AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); 
1455         ncl++;
1456         cl2->SetChargeRatio(ratio);
1457         cl2->SetType(100+cpositive[j]+cnegative[i]);
1458         //cl2->SetType(0);
1459         /*
1460           if (pairs[i][j]<100){
1461           printf("problem:- %d\n", pairs[i][j]);
1462           }
1463           if (cnegative[i]<2&&cpositive[j]<2){
1464           printf("problem:- %d\n", pairs[i][j]);
1465           }
1466         */
1467       }
1468     }
1469   }
1470
1471 //   for (Int_t i=0; i<1000; i++) delete [] pairs[i];
1472 //   delete [] pairs;
1473
1474 }
1475
1476
1477 void AliITSclustererV2::
1478 FindClustersSSD(const TClonesArray *alldigits, TClonesArray *clusters) {
1479   //------------------------------------------------------------
1480   // Actual SSD cluster finder
1481   //------------------------------------------------------------
1482   Int_t smaxall=alldigits->GetEntriesFast();
1483   if (smaxall==0) return;
1484   TObjArray *digits = new TObjArray;
1485   for (Int_t i=0;i<smaxall; i++){
1486     AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
1487     if (d->GetSignal()<3) continue;
1488     digits->AddLast(d);
1489   }
1490   Int_t smax = digits->GetEntriesFast();
1491   if (smax==0) return;
1492   
1493   const Int_t MAX=1000;
1494   Int_t np=0, nn=0; 
1495   Ali1Dcluster pos[MAX], neg[MAX];
1496   Float_t y=0., q=0., qmax=0.; 
1497   Int_t lab[4]={-2,-2,-2,-2};
1498   
1499   AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0);
1500   q += d->GetSignal();
1501   y += d->GetCoord2()*d->GetSignal();
1502   qmax=d->GetSignal();
1503   lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
1504   Int_t curr=d->GetCoord2();
1505   Int_t flag=d->GetCoord1();
1506   Int_t *n=&nn;
1507   Ali1Dcluster *c=neg;
1508   Int_t nd=1;
1509   Int_t milab[10];
1510   for (Int_t ilab=0;ilab<10;ilab++){
1511     milab[ilab]=-2;
1512   }
1513   milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
1514
1515   for (Int_t s=1; s<smax; s++) {
1516       d=(AliITSdigitSSD*)digits->UncheckedAt(s);      
1517       Int_t strip=d->GetCoord2();
1518       if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
1519          c[*n].SetY(y/q);
1520          c[*n].SetQ(q);
1521          c[*n].SetNd(nd);
1522          CheckLabels2(milab);
1523          c[*n].SetLabels(milab);
1524          //Split suspiciously big cluster
1525          /*
1526          if (nd>10&&nd<16){
1527            c[*n].SetY(y/q-0.3*nd);
1528            c[*n].SetQ(0.5*q);
1529            (*n)++;
1530            if (*n==MAX) {
1531              Error("FindClustersSSD","Too many 1D clusters !");
1532               return;
1533            }
1534            c[*n].SetY(y/q-0.0*nd);
1535            c[*n].SetQ(0.5*q);
1536            c[*n].SetNd(nd);
1537            (*n)++;
1538            if (*n==MAX) {
1539              Error("FindClustersSSD","Too many 1D clusters !");
1540               return;
1541            }
1542            //
1543            c[*n].SetY(y/q+0.3*nd);
1544            c[*n].SetQ(0.5*q);
1545            c[*n].SetNd(nd);
1546            c[*n].SetLabels(milab);
1547          }
1548          else{
1549          */
1550          if (nd>4&&nd<25) {
1551            c[*n].SetY(y/q-0.25*nd);
1552            c[*n].SetQ(0.5*q);
1553            (*n)++;
1554            if (*n==MAX) {
1555              Error("FindClustersSSD","Too many 1D clusters !");
1556              return;
1557            }
1558            c[*n].SetY(y/q+0.25*nd);
1559            c[*n].SetQ(0.5*q);
1560            c[*n].SetNd(nd);
1561            c[*n].SetLabels(milab);
1562          }       
1563          (*n)++;
1564          if (*n==MAX) {
1565           Error("FindClustersSSD","Too many 1D clusters !");
1566           return;
1567          }
1568          y=q=qmax=0.;
1569          nd=0;
1570          lab[0]=lab[1]=lab[2]=-2;
1571          //
1572          for (Int_t ilab=0;ilab<10;ilab++){
1573            milab[ilab]=-2;
1574          }
1575          //
1576          if (flag!=d->GetCoord1()) { n=&np; c=pos; }
1577       }
1578       flag=d->GetCoord1();
1579       q += d->GetSignal();
1580       y += d->GetCoord2()*d->GetSignal();
1581       nd++;
1582       if (d->GetSignal()>qmax) {
1583          qmax=d->GetSignal();
1584          lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
1585       }
1586       for (Int_t ilab=0;ilab<10;ilab++) {
1587         if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab))); 
1588       }
1589       curr=strip;
1590   }
1591   c[*n].SetY(y/q);
1592   c[*n].SetQ(q);
1593   c[*n].SetNd(nd);
1594   c[*n].SetLabels(lab);
1595   //Split suspiciously big cluster
1596   if (nd>4 && nd<25) {
1597      c[*n].SetY(y/q-0.25*nd);
1598      c[*n].SetQ(0.5*q);
1599      (*n)++;
1600      if (*n==MAX) {
1601         Error("FindClustersSSD","Too many 1D clusters !");
1602         return;
1603      }
1604      c[*n].SetY(y/q+0.25*nd);
1605      c[*n].SetQ(0.5*q);
1606      c[*n].SetNd(nd);
1607      c[*n].SetLabels(lab);
1608   }
1609   (*n)++;
1610   if (*n==MAX) {
1611      Error("FindClustersSSD","Too many 1D clusters !");
1612      return;
1613   }
1614
1615   FindClustersSSD(neg, nn, pos, np, clusters);
1616 }
1617
1618 void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, 
1619                                         TClonesArray** clusters) 
1620 {
1621   //------------------------------------------------------------
1622   // Actual SSD cluster finder for raw data
1623   //------------------------------------------------------------
1624   Int_t nClustersSSD = 0;
1625   const Int_t MAX = 1000;
1626   Ali1Dcluster clusters1D[2][MAX];
1627   Int_t nClusters[2] = {0, 0};
1628   Int_t lab[3]={-2,-2,-2};
1629   Float_t q = 0.;
1630   Float_t y = 0.;
1631   Int_t nDigits = 0;
1632   Int_t prevStrip = -1;
1633   Int_t prevFlag = -1;
1634   Int_t prevModule = -1;
1635
1636   // read raw data input stream
1637   while (kTRUE) {
1638     Bool_t next = input->Next();
1639
1640     if(input->GetSignal()<3 && next) continue;
1641     // check if a new cluster starts
1642     Int_t strip = input->GetCoord2();
1643     Int_t flag = input->GetCoord1();
1644     if ((!next || (input->GetModuleID() != prevModule)||
1645          (strip-prevStrip > 1) || (flag != prevFlag)) &&
1646         (nDigits > 0)) {
1647       if (nClusters[prevFlag] == MAX) {
1648         Error("FindClustersSSD", "Too many 1D clusters !");
1649         return;
1650       }
1651       Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++];
1652       cluster.SetY(y/q);
1653       cluster.SetQ(q);
1654       cluster.SetNd(nDigits);
1655       cluster.SetLabels(lab);
1656
1657       //Split suspiciously big cluster
1658       if (nDigits > 4&&nDigits < 25) {
1659         cluster.SetY(y/q - 0.25*nDigits);
1660         cluster.SetQ(0.5*q);
1661         if (nClusters[prevFlag] == MAX) {
1662           Error("FindClustersSSD", "Too many 1D clusters !");
1663           return;
1664         }
1665         Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++];
1666         cluster2.SetY(y/q + 0.25*nDigits);
1667         cluster2.SetQ(0.5*q);
1668         cluster2.SetNd(nDigits);
1669         cluster2.SetLabels(lab);
1670       }
1671       y = q = 0.;
1672       nDigits = 0;
1673     }
1674
1675     if (!next || (input->GetModuleID() != prevModule)) {
1676       Int_t iModule = prevModule;
1677
1678       // when all data from a module was read, search for clusters
1679       if (prevFlag >= 0) {
1680         clusters[iModule] = new TClonesArray("AliITSclusterV2");
1681         fI = iModule;
1682         FindClustersSSD(&clusters1D[0][0], nClusters[0], 
1683                         &clusters1D[1][0], nClusters[1], clusters[iModule]);
1684         Int_t noClusters = clusters[iModule]->GetEntriesFast();
1685         nClustersSSD += noClusters;
1686       }
1687
1688       if (!next) break;
1689       nClusters[0] = nClusters[1] = 0;
1690       y = q = 0.;
1691       nDigits = 0;
1692     }
1693
1694     // add digit to current cluster
1695     q += input->GetSignal();
1696     y += strip * input->GetSignal();
1697     nDigits++;
1698     prevStrip = strip;
1699     prevFlag = flag;
1700     prevModule = input->GetModuleID();
1701
1702   }
1703
1704   Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
1705 }
1706