ITS/AliITSClusterFinderV2SSD.cxx

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