ITS/AliITSClusterFinderV2SSD.cxx

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