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