5e040ef41d2f9aefbe5656982a1eb86a69d4fb42
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderV2SSD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, 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 ITS clusterer V2 class                //
20 //                                                                        //
21 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch            //
22 //          Revised: Enrico Fragiacomo, enrico.fragiacomo@ts.infn.it      //
23 //                                                                        //
24 ///////////////////////////////////////////////////////////////////////////
25
26 #include <Riostream.h>
27
28
29 #include "AliITSClusterFinderV2SSD.h"
30 #include "AliITSRecPoint.h"
31 #include "AliITSgeomTGeo.h"
32 #include "AliITSDetTypeRec.h"
33 #include "AliRawReader.h"
34 #include "AliITSRawStreamSSD.h"
35 #include <TClonesArray.h>
36 #include "AliITSdigitSSD.h"
37 #include "AliITSCalibrationSSD.h"
38
39 Short_t *AliITSClusterFinderV2SSD::fgPairs = 0x0;
40 Int_t    AliITSClusterFinderV2SSD::fgPairsSize = 0;
41
42 ClassImp(AliITSClusterFinderV2SSD)
43
44
45 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinderV2(dettyp),
46 fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1),
47 fYpitchSSD(0.0095),
48 fHwSSD(3.65),
49 fHlSSD(2.00),
50 fTanP(0.0275),
51 fTanN(0.0075){
52
53   //Default constructor
54
55 }
56  
57 //______________________________________________________________________
58 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinderV2(cf),                                             fLastSSD1(cf.fLastSSD1),
59 fYpitchSSD(cf.fYpitchSSD),
60 fHwSSD(cf.fHwSSD),
61 fHlSSD(cf.fHlSSD),
62 fTanP(cf.fTanP),
63 fTanN(cf.fTanN)
64 {
65   // Copy constructor
66
67 }
68
69 //______________________________________________________________________
70 AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD&  cf ){
71   // Assignment operator
72
73   this->~AliITSClusterFinderV2SSD();
74   new(this) AliITSClusterFinderV2SSD(cf);
75   return *this;
76 }
77
78
79 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
80
81   //Find clusters V2
82   SetModule(mod);
83   FindClustersSSD(fDigits);
84
85 }
86
87 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
88   //------------------------------------------------------------
89   // Actual SSD cluster finder
90   //------------------------------------------------------------
91   AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
92   Float_t gain=0;
93
94   Int_t smaxall=alldigits->GetEntriesFast();
95   if (smaxall==0) return;
96   //  TObjArray *digits = new TObjArray;
97   TObjArray digits;
98   for (Int_t i=0;i<smaxall; i++){
99     AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
100
101     if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());  
102     else gain = cal->GetGainN(d->GetStripNumber());
103
104     Float_t q=gain*d->GetSignal(); // calibration brings mip peaks around 120 (in ADC units)
105     q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
106     //Float_t q=d->GetSignal()/4.29;// temp. fix (for PID purposed - normalis. to be checked)
107     d->SetSignal(Int_t(q));
108
109     if (d->GetSignal()<3) continue;
110     digits.AddLast(d);
111   }
112   Int_t smax = digits.GetEntriesFast();
113   if (smax==0) return;
114   
115   const Int_t kMax=1000;
116   Int_t np=0, nn=0; 
117   Ali1Dcluster pos[kMax], neg[kMax];
118   Float_t y=0., q=0., qmax=0.; 
119   Int_t lab[4]={-2,-2,-2,-2};
120   
121   AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
122   q += d->GetSignal();
123   y += d->GetCoord2()*d->GetSignal();
124   qmax=d->GetSignal();
125   lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
126   Int_t curr=d->GetCoord2();
127   Int_t flag=d->GetCoord1();
128   Int_t *n=&nn;
129   Ali1Dcluster *c=neg;
130   Int_t nd=1;
131   Int_t milab[10];
132   for (Int_t ilab=0;ilab<10;ilab++){
133     milab[ilab]=-2;
134   }
135   milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
136
137   for (Int_t s=1; s<smax; s++) {
138       d=(AliITSdigitSSD*)digits.UncheckedAt(s);      
139       Int_t strip=d->GetCoord2();
140       if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
141          c[*n].SetY(y/q);
142          c[*n].SetQ(q);
143          c[*n].SetNd(nd);
144          CheckLabels2(milab);
145          c[*n].SetLabels(milab);
146          //Split suspiciously big cluster
147          if (nd>4&&nd<25) {
148            c[*n].SetY(y/q-0.25*nd);
149            c[*n].SetQ(0.5*q);
150            (*n)++;
151            if (*n==kMax) {
152              Error("FindClustersSSD","Too many 1D clusters !");
153              return;
154            }
155            c[*n].SetY(y/q+0.25*nd);
156            c[*n].SetQ(0.5*q);
157            c[*n].SetNd(nd);
158            c[*n].SetLabels(milab);
159          }       
160          (*n)++;
161          if (*n==kMax) {
162           Error("FindClustersSSD","Too many 1D clusters !");
163           return;
164          }
165          y=q=qmax=0.;
166          nd=0;
167          lab[0]=lab[1]=lab[2]=-2;
168          //
169          for (Int_t ilab=0;ilab<10;ilab++){
170            milab[ilab]=-2;
171          }
172          //
173          if (flag!=d->GetCoord1()) { n=&np; c=pos; }
174       }
175       flag=d->GetCoord1();
176       q += d->GetSignal();
177       y += d->GetCoord2()*d->GetSignal();
178       nd++;
179       if (d->GetSignal()>qmax) {
180          qmax=d->GetSignal();
181          lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
182       }
183       for (Int_t ilab=0;ilab<10;ilab++) {
184         if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab))); 
185       }
186       curr=strip;
187   }
188   c[*n].SetY(y/q);
189   c[*n].SetQ(q);
190   c[*n].SetNd(nd);
191   c[*n].SetLabels(lab);
192   //Split suspiciously big cluster
193   if (nd>4 && nd<25) {
194      c[*n].SetY(y/q-0.25*nd);
195      c[*n].SetQ(0.5*q);
196      (*n)++;
197      if (*n==kMax) {
198         Error("FindClustersSSD","Too many 1D clusters !");
199         return;
200      }
201      c[*n].SetY(y/q+0.25*nd);
202      c[*n].SetQ(0.5*q);
203      c[*n].SetNd(nd);
204      c[*n].SetLabels(lab);
205   }
206   (*n)++;
207   if (*n==kMax) {
208      Error("FindClustersSSD","Too many 1D clusters !");
209      return;
210   }
211
212   FindClustersSSD(neg, nn, pos, np);
213 }
214
215
216 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
217
218     //------------------------------------------------------------
219   // This function creates ITS clusters from raw data
220   //------------------------------------------------------------
221   rawReader->Reset();
222   AliITSRawStreamSSD inputSSD(rawReader);
223   //  rawReader->SelectEquipment(-1,0,15);
224   FindClustersSSD(&inputSSD,clusters);
225   
226 }
227
228 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input, 
229                                         TClonesArray** clusters) 
230 {
231   //------------------------------------------------------------
232   // Actual SSD cluster finder for raw data
233   //------------------------------------------------------------
234   Int_t nClustersSSD = 0;
235   const Int_t kMax = 1000;
236   Ali1Dcluster clusters1D[2][kMax];
237   Int_t nClusters[2] = {0, 0};
238   Int_t lab[3]={-2,-2,-2};
239   Float_t q = 0.;
240   Float_t y = 0.;
241   Int_t nDigits = 0;
242   Float_t gain=0;
243   Float_t noise=0.;
244   Float_t oldnoise=0.;
245   AliITSCalibrationSSD* cal=NULL;
246
247   Int_t matrix[12][1536];
248   Int_t iddl=-1;
249   Int_t iad=-1;
250   Int_t oddl = -1;
251   Int_t oad = -1;
252   Int_t oadc = -1;
253   Int_t ostrip = -1;
254   Int_t osignal = 65535;
255   Int_t n=0;
256   Bool_t next=0;
257
258   // read raw data input stream
259   while (kTRUE) {
260     
261     // reset signal matrix
262     for(Int_t i=0; i<12; i++) { for(Int_t j=0; j<1536; j++) { matrix[i][j] = 65535;} }
263
264     if(osignal!=65535) { 
265       n++;
266       matrix[oadc][ostrip] = osignal; // recover data from previous occurence of input->Next() 
267     }
268
269     // buffer data for ddl=iddl and ad=iad
270     while(kTRUE) {
271
272       next = input->Next();
273       if((!next)&&(input->flag)) continue;
274       Int_t ddl=input->GetDDL(); 
275       Int_t ad=input->GetAD();
276       Int_t adc = input->GetADC(); adc = (adc<6)? adc : adc - 2;
277       Int_t strip = input->GetStrip();
278       if(input->GetSideFlag()) strip=1535-strip;
279       Int_t signal = input->GetSignal();
280       //cout<<ddl<<" "<<ad<<" "<<adc<<" "<<strip<<" "<<signal<<endl;
281
282       if((ddl==iddl)&&(ad==iad)) {n++; matrix[adc][strip] = signal;}
283       else {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
284       
285       if(!next) break;
286     }
287
288     // No SDD data
289     if(!next && oddl<0) break;
290
291     if(n==0) continue; // first occurence
292     n=0; osignal=0;
293
294     // fill 1Dclusters
295     for(Int_t iadc=0; iadc<12; iadc++) {  // loop over ADC index for ddl=oddl and ad=oad
296
297       Int_t iimod = (oad - 1)  * 12 + iadc;
298       Int_t iModule = AliITSRawStreamSSD::GetModuleNumber(oddl,iimod);
299       if(iModule==-1) continue;
300       //cout<<"ddl="<<oddl<<" ad"<<oad<<" module="<<iModule<<endl;
301       cal = (AliITSCalibrationSSD*)GetResp(iModule);
302
303       Bool_t first = 0;
304
305       for(Int_t istrip=0; istrip<768; istrip++) { // P-side
306
307         Int_t signal = matrix[iadc][istrip];
308         oldnoise = noise;
309         noise = cal->GetNoiseP(istrip);
310         if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
311
312         if (signal!=65535) {
313           gain = cal->GetGainP(istrip);
314           signal = (Int_t) ( signal * gain ); // signal is corrected for gain
315           //cout<<"ddl="<<oddl<<" ad"<<oad<<" module="<<iModule<<" strip= "<<istrip<<
316           //        " sig="<<signal<<endl;
317           signal = (Int_t) cal->ADCToKeV( signal ); // signal is  converted in KeV 
318          
319           q += signal;    // add digit to current cluster
320           y += istrip * signal;   
321           nDigits++;
322           first=1;
323         }
324
325         else if(first) {
326
327           if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
328
329             Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
330             cluster.SetY(y/q);
331             cluster.SetQ(q);
332             cluster.SetNd(nDigits);
333             cluster.SetLabels(lab);
334             
335             //Split suspiciously big cluster
336             if (nDigits > 4&&nDigits < 25) {
337               cluster.SetY(y/q - 0.25*nDigits);
338               cluster.SetQ(0.5*q);
339               if (nClusters[0] == kMax) {
340                 Error("FindClustersSSD", "Too many 1D clusters !");
341                 return;
342               }
343               Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
344               cluster2.SetY(y/q + 0.25*nDigits);
345               cluster2.SetQ(0.5*q);
346               cluster2.SetNd(nDigits);
347               cluster2.SetLabels(lab);
348             }
349             
350           }
351           
352           y = q = 0.;
353           nDigits = 0;
354           first=0;
355         }
356
357       } // loop over strip on P-side
358
359       // if last strip does have signal
360       if(first) {
361         
362           if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
363
364             Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
365             cluster.SetY(y/q);
366             cluster.SetQ(q);
367             cluster.SetNd(nDigits);
368             cluster.SetLabels(lab);
369             
370             //Split suspiciously big cluster
371             if (nDigits > 4&&nDigits < 25) {
372               cluster.SetY(y/q - 0.25*nDigits);
373               cluster.SetQ(0.5*q);
374               if (nClusters[0] == kMax) {
375                 Error("FindClustersSSD", "Too many 1D clusters !");
376                 return;
377               }
378               Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
379               cluster2.SetY(y/q + 0.25*nDigits);
380               cluster2.SetQ(0.5*q);
381               cluster2.SetNd(nDigits);
382               cluster2.SetLabels(lab);
383             }
384             
385           }
386           y = q = 0.;
387           nDigits = 0;
388           first=0;
389       }
390       
391       oldnoise = 0.;
392       noise = 0.;
393       for(Int_t istrip=768; istrip<1536; istrip++) { // N-side
394         
395         Int_t signal = matrix[iadc][istrip];
396         Int_t strip = 1535-istrip;
397         oldnoise = noise;
398         noise = cal->GetNoiseN(strip);
399         if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
400
401         if (signal!=65535) {
402           gain = cal->GetGainN(strip);
403           signal = (Int_t) ( signal * gain); // signal is corrected for gain
404           signal = (Int_t) cal->ADCToKeV( signal ); // signal is  converted in KeV 
405           
406           // add digit to current cluster
407           q += signal;
408           y += strip * signal;
409           nDigits++;
410           first=1;
411         }
412
413         else if(first) {
414
415           if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
416             
417             Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
418             cluster.SetY(y/q);
419             cluster.SetQ(q);
420             cluster.SetNd(nDigits);
421             cluster.SetLabels(lab);
422             
423             //Split suspiciously big cluster
424             if (nDigits > 4&&nDigits < 25) {
425               cluster.SetY(y/q - 0.25*nDigits);
426               cluster.SetQ(0.5*q);
427               if (nClusters[1] == kMax) {
428                 Error("FindClustersSSD", "Too many 1D clusters !");
429                 return;
430               }
431               Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
432               cluster2.SetY(y/q + 0.25*nDigits);
433               cluster2.SetQ(0.5*q);
434               cluster2.SetNd(nDigits);
435               cluster2.SetLabels(lab);
436             }
437             
438           }
439           y = q = 0.;
440           nDigits = 0;
441           first=0;        
442         }
443         
444       } // loop over strips on N-side
445
446       if(first) {
447         
448           if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
449
450             Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
451             cluster.SetY(y/q);
452             cluster.SetQ(q);
453             cluster.SetNd(nDigits);
454             cluster.SetLabels(lab);
455             
456             //Split suspiciously big cluster
457             if (nDigits > 4&&nDigits < 25) {
458               cluster.SetY(y/q - 0.25*nDigits);
459               cluster.SetQ(0.5*q);
460               if (nClusters[1] == kMax) {
461                 Error("FindClustersSSD", "Too many 1D clusters !");
462                 return;
463               }
464               Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
465               cluster2.SetY(y/q + 0.25*nDigits);
466               cluster2.SetQ(0.5*q);
467               cluster2.SetNd(nDigits);
468               cluster2.SetLabels(lab);
469             }
470             
471           }
472           y = q = 0.;
473           nDigits = 0;
474           first=0;        
475       }
476       
477       // create recpoints
478       if((nClusters[0])&&(nClusters[1])) {
479  
480         //cout<<"creating recpoint for module="<<iModule<<" "<<nClusters[0]<<" "<<nClusters[1]<<endl;
481         clusters[iModule] = new TClonesArray("AliITSRecPoint");
482         fModule = iModule;
483         //      fModule = 500;
484         FindClustersSSD(&clusters1D[0][0], nClusters[0], 
485                         &clusters1D[1][0], nClusters[1], clusters[iModule]);
486         Int_t nClusters = clusters[iModule]->GetEntriesFast();
487         nClustersSSD += nClusters;
488       }
489
490       nClusters[0] = nClusters[1] = 0;
491       y = q = 0.;
492       nDigits = 0;
493
494     } // loop over adc
495
496     if(!next) break;
497   }
498   
499   Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
500 }
501
502 void AliITSClusterFinderV2SSD::
503 FindClustersSSD(Ali1Dcluster* neg, Int_t nn, 
504                 Ali1Dcluster* pos, Int_t np,
505                 TClonesArray *clusters) {
506   //------------------------------------------------------------
507   // Actual SSD cluster finder
508   //------------------------------------------------------------
509
510   //  Float_t xyz[3];
511
512   const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
513
514   TClonesArray &cl=*clusters;
515   //
516   Float_t tanp=fTanP, tann=fTanN;
517   if (fModule>fLastSSD1) {tann=fTanP; tanp=fTanN;}
518   Int_t idet=fNdet[fModule];
519   Int_t ncl=0;
520   //
521   Int_t negativepair[30000];
522   Int_t cnegative[3000];  
523   Int_t cused1[3000];
524   Int_t positivepair[30000];
525   Int_t cpositive[3000];
526   Int_t cused2[3000];
527   for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
528   for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
529   for (Int_t i=0;i<30000;i++) {negativepair[i]=0; positivepair[i]=0;}
530
531   if ((np*nn) > fgPairsSize) {
532     if (fgPairs) delete [] fgPairs;
533     fgPairsSize = 4*np*nn;
534     fgPairs = new Short_t[fgPairsSize];
535   }
536   memset(fgPairs,0,sizeof(Short_t)*np*nn);
537
538   //
539   // find available pairs
540   //
541   for (Int_t i=0; i<np; i++) {
542     Float_t yp=pos[i].GetY()*fYpitchSSD; 
543     if (pos[i].GetQ()<3) continue;
544     for (Int_t j=0; j<nn; j++) {
545       if (neg[j].GetQ()<3) continue;
546       Float_t yn=neg[j].GetY()*fYpitchSSD;
547       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
548       Float_t yt=yn + tann*zt;
549       zt-=fHlSSD; yt-=fHwSSD;
550       if (TMath::Abs(yt)<fHwSSD+0.01)
551       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
552         negativepair[i*10+cnegative[i]] =j;  //index
553         positivepair[j*10+cpositive[j]] =i;
554         cnegative[i]++;  //counters
555         cpositive[j]++; 
556         fgPairs[i*nn+j]=100;
557       }
558     }
559   }
560
561   //
562   // try to recover points out of but close to the module boundaries 
563   //
564   for (Int_t i=0; i<np; i++) {
565     Float_t yp=pos[i].GetY()*fYpitchSSD; 
566     if (pos[i].GetQ()<3) continue;
567     for (Int_t j=0; j<nn; j++) {
568       if (neg[j].GetQ()<3) continue;
569       // if both 1Dclusters have an other cross continue
570       if (cpositive[j]&&cnegative[i]) continue;
571       Float_t yn=neg[j].GetY()*fYpitchSSD;
572       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
573       Float_t yt=yn + tann*zt;
574       zt-=fHlSSD; yt-=fHwSSD;
575       if (TMath::Abs(yt)<fHwSSD+0.1)
576       if (TMath::Abs(zt)<fHlSSD+0.15) {
577         // tag 1Dcluster (eventually will produce low quality recpoint)
578         if (cnegative[i]==0) pos[i].SetNd(100);  // not available pair
579         if (cpositive[j]==0) neg[j].SetNd(100);  // not available pair
580         negativepair[i*10+cnegative[i]] =j;  //index
581         positivepair[j*10+cpositive[j]] =i;
582         cnegative[i]++;  //counters
583         cpositive[j]++; 
584         fgPairs[i*nn+j]=100;
585       }
586     }
587   }
588
589   //
590   Float_t lp[5];
591   Int_t milab[10];
592   Double_t ratio;
593   
594   //
595   // sign gold tracks
596   //
597   for (Int_t ip=0;ip<np;ip++){
598     Float_t ybest=1000,zbest=1000,qbest=0;
599     //
600     // select gold clusters
601     if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){ 
602       Float_t yp=pos[ip].GetY()*fYpitchSSD; 
603       Int_t j = negativepair[10*ip];      
604       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
605       //
606       Float_t yn=neg[j].GetY()*fYpitchSSD;
607       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
608       Float_t yt=yn + tann*zt;
609       zt-=fHlSSD; yt-=fHwSSD;
610       ybest=yt; zbest=zt; 
611       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
612       {
613       Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
614       mT2L->MasterToLocal(loc,trk);
615       lp[0]=trk[1];
616       lp[1]=trk[2];
617       }
618       lp[2]=0.0025*0.0025;  //SigmaY2
619       lp[3]=0.110*0.110;  //SigmaZ2
620       
621       lp[4]=qbest;        //Q
622       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
623       for (Int_t ilab=0;ilab<3;ilab++){
624         milab[ilab] = pos[ip].GetLabel(ilab);
625         milab[ilab+3] = neg[j].GetLabel(ilab);
626       }
627       //
628       CheckLabels2(milab);
629       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
630       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
631       AliITSRecPoint * cl2;
632     
633       if(clusters){  // Note clusters != 0 when method is called for rawdata
634
635
636         cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
637
638         //      cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
639
640         cl2->SetChargeRatio(ratio);     
641         cl2->SetType(1);
642         fgPairs[ip*nn+j]=1;
643         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
644           cl2->SetType(2);
645           fgPairs[ip*nn+j]=2;
646         }
647         cused1[ip]++;
648         cused2[j]++;
649         
650       }
651       else{ // Note clusters == 0 when method is called for digits
652
653         cl2 = new AliITSRecPoint(milab,lp,info);        
654
655         //      cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
656
657         cl2->SetChargeRatio(ratio);     
658         cl2->SetType(1);
659         fgPairs[ip*nn+j]=1;
660         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
661           cl2->SetType(2);
662           fgPairs[ip*nn+j]=2;
663         }
664         cused1[ip]++;
665         cused2[j]++;
666         //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" gold"<<endl;
667         fDetTypeRec->AddRecPoint(*cl2);
668       }
669       ncl++;
670     }
671   }
672     
673   for (Int_t ip=0;ip<np;ip++){
674     Float_t ybest=1000,zbest=1000,qbest=0;
675     //
676     //
677     // select "silber" cluster
678     if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
679       Int_t in  = negativepair[10*ip];
680       Int_t ip2 = positivepair[10*in];
681       if (ip2==ip) ip2 =  positivepair[10*in+1];
682       Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
683       if (TMath::Abs(pcharge-neg[in].GetQ())<10){
684         //
685         // add first pair
686         if (fgPairs[ip*nn+in]==100){  //
687           Float_t yp=pos[ip].GetY()*fYpitchSSD; 
688           Float_t yn=neg[in].GetY()*fYpitchSSD;
689           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
690           Float_t yt=yn + tann*zt;
691           zt-=fHlSSD; yt-=fHwSSD;
692           ybest =yt;  zbest=zt; 
693           qbest =pos[ip].GetQ();
694           {
695           Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
696           mT2L->MasterToLocal(loc,trk);
697           lp[0]=trk[1];
698           lp[1]=trk[2];
699           }
700           lp[2]=0.0025*0.0025;  //SigmaY2
701           lp[3]=0.110*0.110;  //SigmaZ2
702           
703           lp[4]=qbest;        //Q
704           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
705           for (Int_t ilab=0;ilab<3;ilab++){
706             milab[ilab] = pos[ip].GetLabel(ilab);
707             milab[ilab+3] = neg[in].GetLabel(ilab);
708           }
709           //
710           CheckLabels2(milab);
711           ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
712           milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
713           Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
714
715           AliITSRecPoint * cl2;
716           if(clusters){
717
718             cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
719
720             //  cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
721
722             cl2->SetChargeRatio(ratio);         
723             cl2->SetType(5);
724             fgPairs[ip*nn+in] = 5;
725             if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
726               cl2->SetType(6);
727               fgPairs[ip*nn+in] = 6;
728             }       
729           }
730           else{
731             cl2 = new AliITSRecPoint(milab,lp,info);
732             cl2->SetChargeRatio(ratio);         
733             cl2->SetType(5);
734             fgPairs[ip*nn+in] = 5;
735             if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
736               cl2->SetType(6);
737               fgPairs[ip*nn+in] = 6;
738             }
739             //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silver1"<<endl;
740
741             fDetTypeRec->AddRecPoint(*cl2);
742           }
743           ncl++;
744         }
745         
746         //
747         // add second pair
748         
749       //        if (!(cused1[ip2] || cused2[in])){  //
750         if (fgPairs[ip2*nn+in]==100){
751           Float_t yp=pos[ip2].GetY()*fYpitchSSD;
752           Float_t yn=neg[in].GetY()*fYpitchSSD;
753           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
754           Float_t yt=yn + tann*zt;
755           zt-=fHlSSD; yt-=fHwSSD;
756           ybest =yt;  zbest=zt; 
757           qbest =pos[ip2].GetQ();
758           {
759           Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
760           mT2L->MasterToLocal(loc,trk);
761           lp[0]=trk[1];
762           lp[1]=trk[2];
763           }
764           lp[2]=0.0025*0.0025;  //SigmaY2
765           lp[3]=0.110*0.110;  //SigmaZ2
766           
767           lp[4]=qbest;        //Q
768           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
769           for (Int_t ilab=0;ilab<3;ilab++){
770             milab[ilab] = pos[ip2].GetLabel(ilab);
771             milab[ilab+3] = neg[in].GetLabel(ilab);
772           }
773           //
774           CheckLabels2(milab);
775           ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
776           milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
777           Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
778           
779           AliITSRecPoint * cl2;
780           if(clusters){
781             cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
782
783             //  cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
784
785             cl2->SetChargeRatio(ratio);         
786             cl2->SetType(5);
787             fgPairs[ip2*nn+in] =5;
788             if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
789               cl2->SetType(6);
790               fgPairs[ip2*nn+in] =6;
791             }
792           }
793           else{
794             cl2 = new AliITSRecPoint(milab,lp,info);
795             cl2->SetChargeRatio(ratio);         
796             cl2->SetType(5);
797             fgPairs[ip2*nn+in] =5;
798             if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
799               cl2->SetType(6);
800               fgPairs[ip2*nn+in] =6;
801             }
802             
803             //   cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silver2"<<endl;
804             fDetTypeRec->AddRecPoint(*cl2);
805           }
806           ncl++;
807         }       
808         cused1[ip]++;
809         cused1[ip2]++;
810         cused2[in]++;
811       }
812     }    
813   }
814
815   //  
816   for (Int_t jn=0;jn<nn;jn++){
817     if (cused2[jn]) continue;
818     Float_t ybest=1000,zbest=1000,qbest=0;
819     // select "silber" cluster
820     if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
821       Int_t ip  = positivepair[10*jn];
822       Int_t jn2 = negativepair[10*ip];
823       if (jn2==jn) jn2 =  negativepair[10*ip+1];
824       Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
825       //
826       if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
827         //
828         // add first pair
829         //      if (!(cused1[ip]||cused2[jn])){
830         if (fgPairs[ip*nn+jn]==100){
831           Float_t yn=neg[jn].GetY()*fYpitchSSD; 
832           Float_t yp=pos[ip].GetY()*fYpitchSSD;
833           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
834           Float_t yt=yn + tann*zt;
835           zt-=fHlSSD; yt-=fHwSSD;
836           ybest =yt;  zbest=zt; 
837           qbest =neg[jn].GetQ();
838           {
839           Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
840           mT2L->MasterToLocal(loc,trk);
841           lp[0]=trk[1];
842           lp[1]=trk[2];
843           }
844           lp[2]=0.0025*0.0025;  //SigmaY2
845           lp[3]=0.110*0.110;  //SigmaZ2
846           
847           lp[4]=qbest;        //Q
848           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
849           for (Int_t ilab=0;ilab<3;ilab++){
850             milab[ilab] = pos[ip].GetLabel(ilab);
851             milab[ilab+3] = neg[jn].GetLabel(ilab);
852           }
853           //
854           CheckLabels2(milab);
855           ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
856           milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
857           Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
858
859           AliITSRecPoint * cl2;
860           if(clusters){
861             cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
862
863             //  cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
864
865             cl2->SetChargeRatio(ratio);         
866             cl2->SetType(7);
867             fgPairs[ip*nn+jn] =7;
868             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
869               cl2->SetType(8);
870               fgPairs[ip*nn+jn]=8;
871             }
872
873           }
874           else{
875             cl2 = new AliITSRecPoint(milab,lp,info);
876             cl2->SetChargeRatio(ratio);         
877             cl2->SetType(7);
878             fgPairs[ip*nn+jn] =7;
879             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
880               cl2->SetType(8);
881               fgPairs[ip*nn+jn]=8;
882             }
883             //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silverN1"<<endl;
884
885             fDetTypeRec->AddRecPoint(*cl2);
886           }
887           ncl++;
888         }
889         //
890         // add second pair
891         //      if (!(cused1[ip]||cused2[jn2])){
892         if (fgPairs[ip*nn+jn2]==100){
893           Float_t yn=neg[jn2].GetY()*fYpitchSSD; 
894           Double_t yp=pos[ip].GetY()*fYpitchSSD; 
895           Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
896           Double_t yt=yn + tann*zt;
897           zt-=fHlSSD; yt-=fHwSSD;
898           ybest =yt;  zbest=zt; 
899           qbest =neg[jn2].GetQ();
900           {
901           Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
902           mT2L->MasterToLocal(loc,trk);
903           lp[0]=trk[1];
904           lp[1]=trk[2];
905           }
906           lp[2]=0.0025*0.0025;  //SigmaY2
907           lp[3]=0.110*0.110;  //SigmaZ2
908           
909           lp[4]=qbest;        //Q
910           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
911           for (Int_t ilab=0;ilab<3;ilab++){
912             milab[ilab] = pos[ip].GetLabel(ilab);
913             milab[ilab+3] = neg[jn2].GetLabel(ilab);
914           }
915           //
916           CheckLabels2(milab);
917           ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
918           milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
919           Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
920           AliITSRecPoint * cl2;
921           if(clusters){
922             cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
923
924             //  cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
925
926             cl2->SetChargeRatio(ratio);         
927             fgPairs[ip*nn+jn2]=7;
928             cl2->SetType(7);
929             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
930               cl2->SetType(8);
931               fgPairs[ip*nn+jn2]=8;
932             }
933             
934           }
935           else{
936             cl2 = new AliITSRecPoint(milab,lp,info);
937             cl2->SetChargeRatio(ratio);         
938             fgPairs[ip*nn+jn2]=7;
939             cl2->SetType(7);
940             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
941               cl2->SetType(8);
942               fgPairs[ip*nn+jn2]=8;
943             }
944             //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silverN2"<<endl;
945
946             fDetTypeRec->AddRecPoint(*cl2);
947           }
948
949           ncl++;
950         }
951         cused1[ip]++;
952         cused2[jn]++;
953         cused2[jn2]++;
954       }
955     }    
956   }
957   
958   for (Int_t ip=0;ip<np;ip++){
959     Float_t ybest=1000,zbest=1000,qbest=0;
960     //
961     // 2x2 clusters
962     //
963     if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){ 
964       Float_t minchargediff =4.;
965       Int_t j=-1;
966       for (Int_t di=0;di<cnegative[ip];di++){
967         Int_t   jc = negativepair[ip*10+di];
968         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
969         if (TMath::Abs(chargedif)<minchargediff){
970           j =jc;
971           minchargediff = TMath::Abs(chargedif);
972         }
973       }
974       if (j<0) continue;  // not proper cluster      
975
976       Int_t count =0;
977       for (Int_t di=0;di<cnegative[ip];di++){
978         Int_t   jc = negativepair[ip*10+di];
979         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
980         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
981       }
982       if (count>1) continue;  // more than one "proper" cluster for positive
983       //
984       count =0;
985       for (Int_t dj=0;dj<cpositive[j];dj++){
986         Int_t   ic  = positivepair[j*10+dj];
987         Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
988         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
989       }
990       if (count>1) continue;  // more than one "proper" cluster for negative
991       
992       Int_t jp = 0;
993       
994       count =0;
995       for (Int_t dj=0;dj<cnegative[jp];dj++){
996         Int_t   ic = positivepair[jp*10+dj];
997         Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
998         if (TMath::Abs(chargedif)<minchargediff+4.) count++;
999       }
1000       if (count>1) continue;   
1001       if (fgPairs[ip*nn+j]<100) continue;
1002       //
1003       //almost gold clusters
1004       Float_t yp=pos[ip].GetY()*fYpitchSSD; 
1005       Float_t yn=neg[j].GetY()*fYpitchSSD;
1006       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1007       Float_t yt=yn + tann*zt;
1008       zt-=fHlSSD; yt-=fHwSSD;
1009       ybest=yt; zbest=zt; 
1010       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1011       {
1012       Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1013       mT2L->MasterToLocal(loc,trk);
1014       lp[0]=trk[1];
1015       lp[1]=trk[2];
1016       }
1017       lp[2]=0.0025*0.0025;  //SigmaY2
1018       lp[3]=0.110*0.110;  //SigmaZ2     
1019       lp[4]=qbest;        //Q
1020       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1021       for (Int_t ilab=0;ilab<3;ilab++){
1022         milab[ilab] = pos[ip].GetLabel(ilab);
1023         milab[ilab+3] = neg[j].GetLabel(ilab);
1024       }
1025       //
1026       CheckLabels2(milab);
1027       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1028       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1029       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1030       AliITSRecPoint * cl2;
1031       if(clusters){
1032         cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1033
1034         //      cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
1035
1036         cl2->SetChargeRatio(ratio);     
1037         cl2->SetType(10);
1038         fgPairs[ip*nn+j]=10;
1039         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1040           cl2->SetType(11);
1041           fgPairs[ip*nn+j]=11;
1042         }
1043         cused1[ip]++;
1044         cused2[j]++;      
1045       }
1046       else{
1047         cl2 = new AliITSRecPoint(milab,lp,info);
1048         cl2->SetChargeRatio(ratio);     
1049         cl2->SetType(10);
1050         fgPairs[ip*nn+j]=10;
1051         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1052           cl2->SetType(11);
1053           fgPairs[ip*nn+j]=11;
1054         }
1055         cused1[ip]++;
1056         cused2[j]++;      
1057         
1058         //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" 2x2"<<endl;
1059
1060         fDetTypeRec->AddRecPoint(*cl2);
1061       }      
1062       ncl++;
1063     }
1064
1065   }
1066   
1067   //  
1068   for (Int_t i=0; i<np; i++) {
1069     Float_t ybest=1000,zbest=1000,qbest=0;
1070     Float_t yp=pos[i].GetY()*fYpitchSSD; 
1071     if (pos[i].GetQ()<3) continue;
1072     for (Int_t j=0; j<nn; j++) {
1073     //    for (Int_t di = 0;di<cpositive[i];di++){
1074     //  Int_t j = negativepair[10*i+di];
1075       if (neg[j].GetQ()<3) continue;
1076       if (cused2[j]||cused1[i]) continue;      
1077       if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1078       ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());      
1079       Float_t yn=neg[j].GetY()*fYpitchSSD;
1080       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1081       Float_t yt=yn + tann*zt;
1082       zt-=fHlSSD; yt-=fHwSSD;
1083       if (TMath::Abs(yt)<fHwSSD+0.01)
1084       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1085         ybest=yt; zbest=zt; 
1086         qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1087         {
1088         Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1089         mT2L->MasterToLocal(loc,trk);
1090         lp[0]=trk[1];
1091         lp[1]=trk[2];
1092         }
1093         lp[2]=0.0025*0.0025;  //SigmaY2
1094         lp[3]=0.110*0.110;  //SigmaZ2
1095
1096         lp[4]=qbest;        //Q
1097         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1098         for (Int_t ilab=0;ilab<3;ilab++){
1099           milab[ilab] = pos[i].GetLabel(ilab);
1100           milab[ilab+3] = neg[j].GetLabel(ilab);
1101         }
1102         //
1103         CheckLabels2(milab);
1104         milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1105         Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1106         AliITSRecPoint * cl2;
1107         if(clusters){
1108           cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1109
1110           //    cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
1111
1112           cl2->SetChargeRatio(ratio);
1113           cl2->SetType(100+cpositive[j]+cnegative[i]);    
1114         }
1115         else{
1116           cl2 = new AliITSRecPoint(milab,lp,info);
1117           cl2->SetChargeRatio(ratio);
1118           cl2->SetType(100+cpositive[j]+cnegative[i]);
1119           
1120           //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" other"<<endl;
1121
1122           fDetTypeRec->AddRecPoint(*cl2);
1123         }
1124         ncl++;
1125         //cl2->SetType(0);
1126         /*
1127           if (fgPairs[i*nn+j]<100){
1128           printf("problem:- %d\n", fgPairs[i*nn+j]);
1129           }
1130           if (cnegative[i]<2&&cpositive[j]<2){
1131           printf("problem:- %d\n", fgPairs[i*nn+j]);
1132           }
1133         */
1134       }
1135     }
1136   }
1137
1138 }
1139
1140