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