ITS/AliITSClusterFinderV2SSD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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   /*
 223   const UInt_t *evid; evid = rawReader->GetEventId();
 224   cout<<"Event="<<evid[0]<<endl;
 225   */
 226   AliITSRawStreamSSD inputSSD(rawReader);
 227   //  rawReader->SelectEquipment(-1,0,15);
 228   FindClustersSSD(&inputSSD,clusters);
 229
 230 }
 231
 232 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input,
 233                                         TClonesArray** clusters)
 234 {
 235   //------------------------------------------------------------
 236   // Actual SSD cluster finder for raw data
 237   //------------------------------------------------------------
 238   Int_t nClustersSSD = 0;
 239   const Int_t kMax = 1000;
 240   Ali1Dcluster clusters1D[2][kMax];
 241   Int_t nClusters[2] = {0, 0};
 242   Int_t lab[3]={-2,-2,-2};
 243   Float_t q = 0.;
 244   Float_t y = 0.;
 245   Int_t nDigits = 0;
 246   Float_t gain=0;
 247   Float_t noise=0.;
 248   Float_t pedestal=0.;
 249   Float_t oldnoise=0.;
 250   AliITSCalibrationSSD* cal=NULL;
 251
 252   Int_t matrix[12][1536];
 253   Int_t iddl=-1;
 254   Int_t iad=-1;
 255   Int_t oddl = -1;
 256   Int_t oad = -1;
 257   Int_t oadc = -1;
 258   Int_t ostrip = -1;
 259   Int_t osignal = 65535;
 260   Int_t n=0;
 261   Bool_t next=0;
 262
 263   // read raw data input stream
 264   while (kTRUE) {
 265
 266     // reset signal matrix
 267     for(Int_t i=0; i<12; i++) { for(Int_t j=0; j<1536; j++) { matrix[i][j] = 65535;} }
 268
 269     if(osignal!=65535) {
 270       n++;
 271       matrix[oadc][ostrip] = osignal; // recover data from previous occurence of input->Next()
 272     }
 273
 274     // buffer data for ddl=iddl and ad=iad
 275     while(kTRUE) {
 276
 277       next = input->Next();
 278       if((!next)&&(input->flag)) continue;
 279       Int_t ddl=input->GetDDL();
 280       Int_t ad=input->GetAD();
 281       Int_t adc = input->GetADC(); adc = (adc<6)? adc : adc - 2;
 282       Int_t strip = input->GetStrip();
 283       if(input->GetSideFlag()) strip=1535-strip;
 284       Int_t signal = input->GetSignal();
 285       //cout<<ddl<<" "<<ad<<" "<<adc<<" "<<strip<<" "<<signal<<endl;
 286
 287       if((ddl==iddl)&&(ad==iad)) {n++; matrix[adc][strip] = signal;}
 288       else {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
 289
 290       if(!next)  {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
 291       //break;
 292     }
 293
 294     // No SSD data
 295     if(!next && oddl<0) break;
 296
 297     if(n==0) continue; // first occurence
 298     n=0; osignal=0;
 299
 300     // fill 1Dclusters
 301     for(Int_t iadc=0; iadc<12; iadc++) {  // loop over ADC index for ddl=oddl and ad=oad
 302
 303       Int_t iimod = (oad - 1)  * 12 + iadc;
 304       Int_t iModule = AliITSRawStreamSSD::GetModuleNumber(oddl,iimod);
 305       if(iModule==-1) continue;
 306       //      cout<<"ddl="<<oddl<<" ad"<<oad<<" module="<<iModule<<endl;
 307       cal = (AliITSCalibrationSSD*)GetResp(iModule);
 308
 309       Bool_t first = 0;
 310
 311       for(Int_t istrip=0; istrip<768; istrip++) { // P-side
 312         Int_t signal = matrix[iadc][istrip];
 313         pedestal = cal->GetPedestalP(istrip);
 314         matrix[iadc][istrip]=signal-pedestal;
 315       }
 316
 317       Float_t cmode=0;
 318       for(Int_t l=0; l<6; l++) {
 319         cmode=0;
 320         for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
 321         cmode/=88.;
 322         for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=cmode;
 323
 324       }
 325
 326       for(Int_t istrip=0; istrip<768; istrip++) { // P-side
 327
 328         Int_t signal = TMath::Abs(matrix[iadc][istrip]);
 329
 330         oldnoise = noise;
 331         noise = cal->GetNoiseP(istrip); if(noise<1.) signal = 65535;
 332         if(signal<5*noise) signal = 65535; // in case ZS was not done in hw do it now
 333         if( (signal<30.) || (istrip<10) || (istrip>758) ) signal=65535;
 334
 335         if (signal!=65535) {
 336           gain = cal->GetGainP(istrip);
 337           signal = (Int_t) ( signal * gain ); // signal is corrected for gain
 338           signal = (Int_t) cal->ADCToKeV( signal ); // signal is  converted in KeV
 339
 340           q += signal;    // add digit to current cluster
 341           y += istrip * signal;
 342           nDigits++;
 343           first=1;
 344         }
 345
 346         else if(first) {
 347
 348           if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
 349
 350             Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
 351             cluster.SetY(y/q);
 352             cluster.SetQ(q);
 353             cluster.SetNd(nDigits);
 354             cluster.SetLabels(lab);
 355
 356             //Split suspiciously big cluster
 357             if (nDigits > 4&&nDigits < 25) {
 358               cluster.SetY(y/q - 0.25*nDigits);
 359               cluster.SetQ(0.5*q);
 360               if (nClusters[0] == kMax) {
 361                 Error("FindClustersSSD", "Too many 1D clusters !");
 362                 return;
 363               }
 364               Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
 365               cluster2.SetY(y/q + 0.25*nDigits);
 366               cluster2.SetQ(0.5*q);
 367               cluster2.SetNd(nDigits);
 368               cluster2.SetLabels(lab);
 369             }
 370
 371           }
 372
 373           y = q = 0.;
 374           nDigits = 0;
 375           first=0;
 376         }
 377
 378       } // loop over strip on P-side
 379
 380       // if last strip does have signal
 381       if(first) {
 382
 383           if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
 384
 385             Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
 386             cluster.SetY(y/q);
 387             cluster.SetQ(q);
 388             cluster.SetNd(nDigits);
 389             cluster.SetLabels(lab);
 390
 391             //Split suspiciously big cluster
 392             if (nDigits > 4&&nDigits < 25) {
 393               cluster.SetY(y/q - 0.25*nDigits);
 394               cluster.SetQ(0.5*q);
 395               if (nClusters[0] == kMax) {
 396                 Error("FindClustersSSD", "Too many 1D clusters !");
 397                 return;
 398               }
 399               Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
 400               cluster2.SetY(y/q + 0.25*nDigits);
 401               cluster2.SetQ(0.5*q);
 402               cluster2.SetNd(nDigits);
 403               cluster2.SetLabels(lab);
 404             }
 405
 406           }
 407           y = q = 0.;
 408           nDigits = 0;
 409           first=0;
 410       }
 411
 412       for(Int_t istrip=768; istrip<1536; istrip++) { // P-side
 413         Int_t signal = matrix[iadc][istrip];
 414         pedestal = cal->GetPedestalN(1535-istrip);
 415         matrix[iadc][istrip]=signal-pedestal;
 416       }
 417
 418       for(Int_t l=6; l<12; l++) {
 419         Float_t cmode=0;
 420         for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
 421         cmode/=88.;
 422         for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=cmode;
 423       }
 424
 425       oldnoise = 0.;
 426       noise = 0.;
 427       for(Int_t istrip=768; istrip<1536; istrip++) { // N-side
 428
 429         Int_t signal = TMath::Abs(matrix[iadc][istrip]);
 430         //cout<<"####"<<" "<<oddl<<" "<<oad<<" "<<iadc<<" "<<istrip<<" "<<signal<<endl;
 431
 432         Int_t strip = 1535-istrip;
 433
 434         oldnoise = noise;
 435         noise = cal->GetNoiseN(strip); if(noise<1.) signal=65535;
 436         if(signal<5*noise) signal = 65535; // in case ZS was not done in hw do it now
 437         if( (signal<30.) || (istrip<778) || (istrip>1526) ) signal=65535;
 438
 439         if (signal!=65535) {
 440           //      cout<<"ddl="<<oddl<<" ad"<<oad<<" module="<<iModule<<" strip= "<<istrip<<
 441           //  " sig="<<signal<<" "<<cal->GetPedestalN(strip)<<endl;
 442           gain = cal->GetGainN(strip);
 443           signal = (Int_t) ( signal * gain); // signal is corrected for gain
 444           signal = (Int_t) cal->ADCToKeV( signal ); // signal is  converted in KeV
 445
 446           // add digit to current cluster
 447           q += signal;
 448           y += strip * signal;
 449           nDigits++;
 450           first=1;
 451         }
 452
 453         else if(first) {
 454
 455           if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
 456
 457             Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
 458             cluster.SetY(y/q);
 459             cluster.SetQ(q);
 460             cluster.SetNd(nDigits);
 461             cluster.SetLabels(lab);
 462
 463             //Split suspiciously big cluster
 464             if (nDigits > 4&&nDigits < 25) {
 465               cluster.SetY(y/q - 0.25*nDigits);
 466               cluster.SetQ(0.5*q);
 467               if (nClusters[1] == kMax) {
 468                 Error("FindClustersSSD", "Too many 1D clusters !");
 469                 return;
 470               }
 471               Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
 472               cluster2.SetY(y/q + 0.25*nDigits);
 473               cluster2.SetQ(0.5*q);
 474               cluster2.SetNd(nDigits);
 475               cluster2.SetLabels(lab);
 476             }
 477
 478           }
 479           y = q = 0.;
 480           nDigits = 0;
 481           first=0;
 482         }
 483
 484       } // loop over strips on N-side
 485
 486       if(first) {
 487
 488           if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
 489
 490             Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
 491             cluster.SetY(y/q);
 492             cluster.SetQ(q);
 493             cluster.SetNd(nDigits);
 494             cluster.SetLabels(lab);
 495
 496             //Split suspiciously big cluster
 497             if (nDigits > 4&&nDigits < 25) {
 498               cluster.SetY(y/q - 0.25*nDigits);
 499               cluster.SetQ(0.5*q);
 500               if (nClusters[1] == kMax) {
 501                 Error("FindClustersSSD", "Too many 1D clusters !");
 502                 return;
 503               }
 504               Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
 505               cluster2.SetY(y/q + 0.25*nDigits);
 506               cluster2.SetQ(0.5*q);
 507               cluster2.SetNd(nDigits);
 508               cluster2.SetLabels(lab);
 509             }
 510
 511           }
 512           y = q = 0.;
 513           nDigits = 0;
 514           first=0;
 515       }
 516
 517       // create recpoints
 518       if((nClusters[0])&&(nClusters[1])) {
 519
 520         //cout<<"creating recpoint for module="<<iModule<<" "<<nClusters[0]<<" "<<nClusters[1]<<endl;
 521         clusters[iModule] = new TClonesArray("AliITSRecPoint");
 522         fModule = iModule;
 523         //      fModule = 500;
 524         FindClustersSSD(&clusters1D[0][0], nClusters[0],
 525                         &clusters1D[1][0], nClusters[1], clusters[iModule]);
 526         Int_t nClusters = clusters[iModule]->GetEntriesFast();
 527         nClustersSSD += nClusters;
 528       }
 529
 530       nClusters[0] = nClusters[1] = 0;
 531       y = q = 0.;
 532       nDigits = 0;
 533
 534     } // loop over adc
 535
 536     if(!next) break;
 537   }
 538
 539   Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
 540 }
 541
 542 void AliITSClusterFinderV2SSD::
 543 FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
 544                 Ali1Dcluster* pos, Int_t np,
 545                 TClonesArray *clusters) {
 546   //------------------------------------------------------------
 547   // Actual SSD cluster finder
 548   //------------------------------------------------------------
 549
 550   //  Float_t xyz[3];
 551
 552   const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
 553
 554   TClonesArray &cl=*clusters;
 555   //
 556   Float_t tanp=fTanP, tann=fTanN;
 557   if (fModule>fLastSSD1) {tann=fTanP; tanp=fTanN;}
 558   Int_t idet=fNdet[fModule];
 559   Int_t ncl=0;
 560   //
 561   Int_t negativepair[30000];
 562   Int_t cnegative[3000];
 563   Int_t cused1[3000];
 564   Int_t positivepair[30000];
 565   Int_t cpositive[3000];
 566   Int_t cused2[3000];
 567   for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
 568   for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
 569   for (Int_t i=0;i<30000;i++) {negativepair[i]=0; positivepair[i]=0;}
 570
 571   if ((np*nn) > fgPairsSize) {
 572     if (fgPairs) delete [] fgPairs;
 573     fgPairsSize = 4*np*nn;
 574     fgPairs = new Short_t[fgPairsSize];
 575   }
 576   memset(fgPairs,0,sizeof(Short_t)*np*nn);
 577
 578   //
 579   // find available pairs
 580   //
 581   for (Int_t i=0; i<np; i++) {
 582     Float_t yp=pos[i].GetY()*fYpitchSSD;
 583     if (pos[i].GetQ()<3) continue;
 584     for (Int_t j=0; j<nn; j++) {
 585       if (neg[j].GetQ()<3) continue;
 586       Float_t yn=neg[j].GetY()*fYpitchSSD;
 587       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 588       Float_t yt=yn + tann*zt;
 589       zt-=fHlSSD; yt-=fHwSSD;
 590       if (TMath::Abs(yt)<fHwSSD+0.01)
 591       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
 592         negativepair[i*10+cnegative[i]] =j;  //index
 593         positivepair[j*10+cpositive[j]] =i;
 594         cnegative[i]++;  //counters
 595         cpositive[j]++;
 596         fgPairs[i*nn+j]=100;
 597       }
 598     }
 599   }
 600
 601   //
 602   // try to recover points out of but close to the module boundaries
 603   //
 604   for (Int_t i=0; i<np; i++) {
 605     Float_t yp=pos[i].GetY()*fYpitchSSD;
 606     if (pos[i].GetQ()<3) continue;
 607     for (Int_t j=0; j<nn; j++) {
 608       if (neg[j].GetQ()<3) continue;
 609       // if both 1Dclusters have an other cross continue
 610       if (cpositive[j]&&cnegative[i]) continue;
 611       Float_t yn=neg[j].GetY()*fYpitchSSD;
 612       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 613       Float_t yt=yn + tann*zt;
 614       zt-=fHlSSD; yt-=fHwSSD;
 615       if (TMath::Abs(yt)<fHwSSD+0.1)
 616       if (TMath::Abs(zt)<fHlSSD+0.15) {
 617         // tag 1Dcluster (eventually will produce low quality recpoint)
 618         if (cnegative[i]==0) pos[i].SetNd(100);  // not available pair
 619         if (cpositive[j]==0) neg[j].SetNd(100);  // not available pair
 620         negativepair[i*10+cnegative[i]] =j;  //index
 621         positivepair[j*10+cpositive[j]] =i;
 622         cnegative[i]++;  //counters
 623         cpositive[j]++;
 624         fgPairs[i*nn+j]=100;
 625       }
 626     }
 627   }
 628
 629   //
 630   Float_t lp[5];
 631   Int_t milab[10];
 632   Double_t ratio;
 633
 634   //
 635   // sign gold tracks
 636   //
 637   for (Int_t ip=0;ip<np;ip++){
 638     Float_t ybest=1000,zbest=1000,qbest=0;
 639     //
 640     // select gold clusters
 641     if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
 642       Float_t yp=pos[ip].GetY()*fYpitchSSD;
 643       Int_t j = negativepair[10*ip];
 644       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
 645       //
 646       Float_t yn=neg[j].GetY()*fYpitchSSD;
 647       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 648       Float_t yt=yn + tann*zt;
 649       zt-=fHlSSD; yt-=fHwSSD;
 650       ybest=yt; zbest=zt;
 651       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
 652
 653       //cout<<yt<<" "<<zt<<" "<<qbest<<endl;
 654
 655       {
 656       Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
 657       mT2L->MasterToLocal(loc,trk);
 658       lp[0]=trk[1];
 659       lp[1]=trk[2];
 660       }
 661       lp[2]=0.0025*0.0025;  //SigmaY2
 662       lp[3]=0.110*0.110;  //SigmaZ2
 663
 664       lp[4]=qbest;        //Q
 665       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 666       for (Int_t ilab=0;ilab<3;ilab++){
 667         milab[ilab] = pos[ip].GetLabel(ilab);
 668         milab[ilab+3] = neg[j].GetLabel(ilab);
 669       }
 670       //
 671       CheckLabels2(milab);
 672       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
 673       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
 674       AliITSRecPoint * cl2;
 675
 676       if(clusters){  // Note clusters != 0 when method is called for rawdata
 677
 678
 679         cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
 680
 681         //      cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
 682
 683         cl2->SetChargeRatio(ratio);
 684         cl2->SetType(1);
 685         fgPairs[ip*nn+j]=1;
 686         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
 687           cl2->SetType(2);
 688           fgPairs[ip*nn+j]=2;
 689         }
 690         cused1[ip]++;
 691         cused2[j]++;
 692
 693       }
 694       else{ // Note clusters == 0 when method is called for digits
 695
 696         cl2 = new AliITSRecPoint(milab,lp,info);
 697
 698         //      cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
 699
 700         cl2->SetChargeRatio(ratio);
 701         cl2->SetType(1);
 702         fgPairs[ip*nn+j]=1;
 703         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
 704           cl2->SetType(2);
 705           fgPairs[ip*nn+j]=2;
 706         }
 707         cused1[ip]++;
 708         cused2[j]++;
 709         //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" gold"<<endl;
 710         fDetTypeRec->AddRecPoint(*cl2);
 711       }
 712       ncl++;
 713     }
 714   }
 715
 716   for (Int_t ip=0;ip<np;ip++){
 717     Float_t ybest=1000,zbest=1000,qbest=0;
 718     //
 719     //
 720     // select "silber" cluster
 721     if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
 722       Int_t in  = negativepair[10*ip];
 723       Int_t ip2 = positivepair[10*in];
 724       if (ip2==ip) ip2 =  positivepair[10*in+1];
 725       Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
 726       if (TMath::Abs(pcharge-neg[in].GetQ())<10){
 727         //
 728         // add first pair
 729         if (fgPairs[ip*nn+in]==100){  //
 730           Float_t yp=pos[ip].GetY()*fYpitchSSD;
 731           Float_t yn=neg[in].GetY()*fYpitchSSD;
 732           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 733           Float_t yt=yn + tann*zt;
 734           zt-=fHlSSD; yt-=fHwSSD;
 735           ybest =yt;  zbest=zt;
 736           qbest =pos[ip].GetQ();
 737           {
 738           Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
 739           mT2L->MasterToLocal(loc,trk);
 740           lp[0]=trk[1];
 741           lp[1]=trk[2];
 742           }
 743           lp[2]=0.0025*0.0025;  //SigmaY2
 744           lp[3]=0.110*0.110;  //SigmaZ2
 745
 746           lp[4]=qbest;        //Q
 747           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 748           for (Int_t ilab=0;ilab<3;ilab++){
 749             milab[ilab] = pos[ip].GetLabel(ilab);
 750             milab[ilab+3] = neg[in].GetLabel(ilab);
 751           }
 752           //
 753           CheckLabels2(milab);
 754           ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
 755           milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
 756           Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
 757
 758           AliITSRecPoint * cl2;
 759           if(clusters){
 760
 761             cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
 762
 763             //  cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
 764
 765             cl2->SetChargeRatio(ratio);
 766             cl2->SetType(5);
 767             fgPairs[ip*nn+in] = 5;
 768             if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
 769               cl2->SetType(6);
 770               fgPairs[ip*nn+in] = 6;
 771             }
 772           }
 773           else{
 774             cl2 = new AliITSRecPoint(milab,lp,info);
 775             cl2->SetChargeRatio(ratio);
 776             cl2->SetType(5);
 777             fgPairs[ip*nn+in] = 5;
 778             if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
 779               cl2->SetType(6);
 780               fgPairs[ip*nn+in] = 6;
 781             }
 782             //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silver1"<<endl;
 783
 784             fDetTypeRec->AddRecPoint(*cl2);
 785           }
 786           ncl++;
 787         }
 788
 789         //
 790         // add second pair
 791
 792       //        if (!(cused1[ip2] || cused2[in])){  //
 793         if (fgPairs[ip2*nn+in]==100){
 794           Float_t yp=pos[ip2].GetY()*fYpitchSSD;
 795           Float_t yn=neg[in].GetY()*fYpitchSSD;
 796           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 797           Float_t yt=yn + tann*zt;
 798           zt-=fHlSSD; yt-=fHwSSD;
 799           ybest =yt;  zbest=zt;
 800           qbest =pos[ip2].GetQ();
 801           {
 802           Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
 803           mT2L->MasterToLocal(loc,trk);
 804           lp[0]=trk[1];
 805           lp[1]=trk[2];
 806           }
 807           lp[2]=0.0025*0.0025;  //SigmaY2
 808           lp[3]=0.110*0.110;  //SigmaZ2
 809
 810           lp[4]=qbest;        //Q
 811           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 812           for (Int_t ilab=0;ilab<3;ilab++){
 813             milab[ilab] = pos[ip2].GetLabel(ilab);
 814             milab[ilab+3] = neg[in].GetLabel(ilab);
 815           }
 816           //
 817           CheckLabels2(milab);
 818           ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
 819           milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
 820           Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
 821
 822           AliITSRecPoint * cl2;
 823           if(clusters){
 824             cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
 825
 826             //  cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
 827
 828             cl2->SetChargeRatio(ratio);
 829             cl2->SetType(5);
 830             fgPairs[ip2*nn+in] =5;
 831             if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
 832               cl2->SetType(6);
 833               fgPairs[ip2*nn+in] =6;
 834             }
 835           }
 836           else{
 837             cl2 = new AliITSRecPoint(milab,lp,info);
 838             cl2->SetChargeRatio(ratio);
 839             cl2->SetType(5);
 840             fgPairs[ip2*nn+in] =5;
 841             if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
 842               cl2->SetType(6);
 843               fgPairs[ip2*nn+in] =6;
 844             }
 845
 846             //   cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silver2"<<endl;
 847             fDetTypeRec->AddRecPoint(*cl2);
 848           }
 849           ncl++;
 850         }
 851         cused1[ip]++;
 852         cused1[ip2]++;
 853         cused2[in]++;
 854       }
 855     }
 856   }
 857
 858   //
 859   for (Int_t jn=0;jn<nn;jn++){
 860     if (cused2[jn]) continue;
 861     Float_t ybest=1000,zbest=1000,qbest=0;
 862     // select "silber" cluster
 863     if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
 864       Int_t ip  = positivepair[10*jn];
 865       Int_t jn2 = negativepair[10*ip];
 866       if (jn2==jn) jn2 =  negativepair[10*ip+1];
 867       Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
 868       //
 869       if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
 870         //
 871         // add first pair
 872         //      if (!(cused1[ip]||cused2[jn])){
 873         if (fgPairs[ip*nn+jn]==100){
 874           Float_t yn=neg[jn].GetY()*fYpitchSSD;
 875           Float_t yp=pos[ip].GetY()*fYpitchSSD;
 876           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 877           Float_t yt=yn + tann*zt;
 878           zt-=fHlSSD; yt-=fHwSSD;
 879           ybest =yt;  zbest=zt;
 880           qbest =neg[jn].GetQ();
 881           {
 882           Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
 883           mT2L->MasterToLocal(loc,trk);
 884           lp[0]=trk[1];
 885           lp[1]=trk[2];
 886           }
 887           lp[2]=0.0025*0.0025;  //SigmaY2
 888           lp[3]=0.110*0.110;  //SigmaZ2
 889
 890           lp[4]=qbest;        //Q
 891           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 892           for (Int_t ilab=0;ilab<3;ilab++){
 893             milab[ilab] = pos[ip].GetLabel(ilab);
 894             milab[ilab+3] = neg[jn].GetLabel(ilab);
 895           }
 896           //
 897           CheckLabels2(milab);
 898           ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
 899           milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
 900           Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
 901
 902           AliITSRecPoint * cl2;
 903           if(clusters){
 904             cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
 905
 906             //  cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
 907
 908             cl2->SetChargeRatio(ratio);
 909             cl2->SetType(7);
 910             fgPairs[ip*nn+jn] =7;
 911             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
 912               cl2->SetType(8);
 913               fgPairs[ip*nn+jn]=8;
 914             }
 915
 916           }
 917           else{
 918             cl2 = new AliITSRecPoint(milab,lp,info);
 919             cl2->SetChargeRatio(ratio);
 920             cl2->SetType(7);
 921             fgPairs[ip*nn+jn] =7;
 922             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
 923               cl2->SetType(8);
 924               fgPairs[ip*nn+jn]=8;
 925             }
 926             //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silverN1"<<endl;
 927
 928             fDetTypeRec->AddRecPoint(*cl2);
 929           }
 930           ncl++;
 931         }
 932         //
 933         // add second pair
 934         //      if (!(cused1[ip]||cused2[jn2])){
 935         if (fgPairs[ip*nn+jn2]==100){
 936           Float_t yn=neg[jn2].GetY()*fYpitchSSD;
 937           Double_t yp=pos[ip].GetY()*fYpitchSSD;
 938           Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 939           Double_t yt=yn + tann*zt;
 940           zt-=fHlSSD; yt-=fHwSSD;
 941           ybest =yt;  zbest=zt;
 942           qbest =neg[jn2].GetQ();
 943           {
 944           Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
 945           mT2L->MasterToLocal(loc,trk);
 946           lp[0]=trk[1];
 947           lp[1]=trk[2];
 948           }
 949           lp[2]=0.0025*0.0025;  //SigmaY2
 950           lp[3]=0.110*0.110;  //SigmaZ2
 951
 952           lp[4]=qbest;        //Q
 953           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 954           for (Int_t ilab=0;ilab<3;ilab++){
 955             milab[ilab] = pos[ip].GetLabel(ilab);
 956             milab[ilab+3] = neg[jn2].GetLabel(ilab);
 957           }
 958           //
 959           CheckLabels2(milab);
 960           ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
 961           milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
 962           Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
 963           AliITSRecPoint * cl2;
 964           if(clusters){
 965             cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
 966
 967             //  cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
 968
 969             cl2->SetChargeRatio(ratio);
 970             fgPairs[ip*nn+jn2]=7;
 971             cl2->SetType(7);
 972             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
 973               cl2->SetType(8);
 974               fgPairs[ip*nn+jn2]=8;
 975             }
 976
 977           }
 978           else{
 979             cl2 = new AliITSRecPoint(milab,lp,info);
 980             cl2->SetChargeRatio(ratio);
 981             fgPairs[ip*nn+jn2]=7;
 982             cl2->SetType(7);
 983             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
 984               cl2->SetType(8);
 985               fgPairs[ip*nn+jn2]=8;
 986             }
 987             //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silverN2"<<endl;
 988
 989             fDetTypeRec->AddRecPoint(*cl2);
 990           }
 991
 992           ncl++;
 993         }
 994         cused1[ip]++;
 995         cused2[jn]++;
 996         cused2[jn2]++;
 997       }
 998     }
 999   }
1000
1001   for (Int_t ip=0;ip<np;ip++){
1002     Float_t ybest=1000,zbest=1000,qbest=0;
1003     //
1004     // 2x2 clusters
1005     //
1006     if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1007       Float_t minchargediff =4.;
1008       Int_t j=-1;
1009       for (Int_t di=0;di<cnegative[ip];di++){
1010         Int_t   jc = negativepair[ip*10+di];
1011         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1012         if (TMath::Abs(chargedif)<minchargediff){
1013           j =jc;
1014           minchargediff = TMath::Abs(chargedif);
1015         }
1016       }
1017       if (j<0) continue;  // not proper cluster
1018
1019       Int_t count =0;
1020       for (Int_t di=0;di<cnegative[ip];di++){
1021         Int_t   jc = negativepair[ip*10+di];
1022         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1023         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1024       }
1025       if (count>1) continue;  // more than one "proper" cluster for positive
1026       //
1027       count =0;
1028       for (Int_t dj=0;dj<cpositive[j];dj++){
1029         Int_t   ic  = positivepair[j*10+dj];
1030         Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1031         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1032       }
1033       if (count>1) continue;  // more than one "proper" cluster for negative
1034
1035       Int_t jp = 0;
1036
1037       count =0;
1038       for (Int_t dj=0;dj<cnegative[jp];dj++){
1039         Int_t   ic = positivepair[jp*10+dj];
1040         Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1041         if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1042       }
1043       if (count>1) continue;
1044       if (fgPairs[ip*nn+j]<100) continue;
1045       //
1046       //almost gold clusters
1047       Float_t yp=pos[ip].GetY()*fYpitchSSD;
1048       Float_t yn=neg[j].GetY()*fYpitchSSD;
1049       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1050       Float_t yt=yn + tann*zt;
1051       zt-=fHlSSD; yt-=fHwSSD;
1052       ybest=yt; zbest=zt;
1053       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1054       {
1055       Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1056       mT2L->MasterToLocal(loc,trk);
1057       lp[0]=trk[1];
1058       lp[1]=trk[2];
1059       }
1060       lp[2]=0.0025*0.0025;  //SigmaY2
1061       lp[3]=0.110*0.110;  //SigmaZ2
1062       lp[4]=qbest;        //Q
1063       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1064       for (Int_t ilab=0;ilab<3;ilab++){
1065         milab[ilab] = pos[ip].GetLabel(ilab);
1066         milab[ilab+3] = neg[j].GetLabel(ilab);
1067       }
1068       //
1069       CheckLabels2(milab);
1070       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1071       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1072       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1073       AliITSRecPoint * cl2;
1074       if(clusters){
1075         cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1076
1077         //      cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
1078
1079         cl2->SetChargeRatio(ratio);
1080         cl2->SetType(10);
1081         fgPairs[ip*nn+j]=10;
1082         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1083           cl2->SetType(11);
1084           fgPairs[ip*nn+j]=11;
1085         }
1086         cused1[ip]++;
1087         cused2[j]++;
1088       }
1089       else{
1090         cl2 = new AliITSRecPoint(milab,lp,info);
1091         cl2->SetChargeRatio(ratio);
1092         cl2->SetType(10);
1093         fgPairs[ip*nn+j]=10;
1094         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1095           cl2->SetType(11);
1096           fgPairs[ip*nn+j]=11;
1097         }
1098         cused1[ip]++;
1099         cused2[j]++;
1100
1101         //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" 2x2"<<endl;
1102
1103         fDetTypeRec->AddRecPoint(*cl2);
1104       }
1105       ncl++;
1106     }
1107
1108   }
1109
1110   //
1111   for (Int_t i=0; i<np; i++) {
1112     Float_t ybest=1000,zbest=1000,qbest=0;
1113     Float_t yp=pos[i].GetY()*fYpitchSSD;
1114     if (pos[i].GetQ()<3) continue;
1115     for (Int_t j=0; j<nn; j++) {
1116     //    for (Int_t di = 0;di<cpositive[i];di++){
1117     //  Int_t j = negativepair[10*i+di];
1118       if (neg[j].GetQ()<3) continue;
1119       if (cused2[j]||cused1[i]) continue;
1120       if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1121       ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1122       Float_t yn=neg[j].GetY()*fYpitchSSD;
1123       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1124       Float_t yt=yn + tann*zt;
1125       zt-=fHlSSD; yt-=fHwSSD;
1126       if (TMath::Abs(yt)<fHwSSD+0.01)
1127       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1128         ybest=yt; zbest=zt;
1129         qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1130         {
1131         Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1132         mT2L->MasterToLocal(loc,trk);
1133         lp[0]=trk[1];
1134         lp[1]=trk[2];
1135         }
1136         lp[2]=0.0025*0.0025;  //SigmaY2
1137         lp[3]=0.110*0.110;  //SigmaZ2
1138
1139         lp[4]=qbest;        //Q
1140         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1141         for (Int_t ilab=0;ilab<3;ilab++){
1142           milab[ilab] = pos[i].GetLabel(ilab);
1143           milab[ilab+3] = neg[j].GetLabel(ilab);
1144         }
1145         //
1146         CheckLabels2(milab);
1147         milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1148         Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1149         AliITSRecPoint * cl2;
1150         if(clusters){
1151           cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1152
1153           //    cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
1154
1155           cl2->SetChargeRatio(ratio);
1156           cl2->SetType(100+cpositive[j]+cnegative[i]);
1157         }
1158         else{
1159           cl2 = new AliITSRecPoint(milab,lp,info);
1160           cl2->SetChargeRatio(ratio);
1161           cl2->SetType(100+cpositive[j]+cnegative[i]);
1162
1163           //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" other"<<endl;
1164
1165           fDetTypeRec->AddRecPoint(*cl2);
1166         }
1167         ncl++;
1168         //cl2->SetType(0);
1169         /*
1170           if (fgPairs[i*nn+j]<100){
1171           printf("problem:- %d\n", fgPairs[i*nn+j]);
1172           }
1173           if (cnegative[i]<2&&cpositive[j]<2){
1174           printf("problem:- %d\n", fgPairs[i*nn+j]);
1175           }
1176         */
1177       }
1178     }
1179   }
1180
1181 }
1182
1183