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