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