ITS/AliITSClusterFinderV2SSD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2003, 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 //            Implementation of the ITS clusterer V2 class                //
  17 //                                                                        //
  18 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch            //
  19 //                                                                        //
  20 ///////////////////////////////////////////////////////////////////////////
  21
  22
  23 #include "AliITSClusterFinderV2SSD.h"
  24 #include "AliITSRecPoint.h"
  25 #include "AliITSDetTypeRec.h"
  26 #include "AliRawReader.h"
  27 #include "AliITSRawStreamSSD.h"
  28 #include <TClonesArray.h>
  29 #include "AliITSdigitSSD.h"
  30
  31 ClassImp(AliITSClusterFinderV2SSD)
  32
  33
  34 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinderV2(dettyp),
  35 fLastSSD1(0),
  36 fYpitchSSD(0.0095),
  37 fHwSSD(3.65),
  38 fHlSSD(2.00),
  39 fTanP(0.0275),
  40 fTanN(0.0075){
  41
  42   //Default constructor
  43
  44   fLastSSD1=fDetTypeRec->GetITSgeom()->GetModuleIndex(6,1,1)-1;
  45
  46 }
  47
  48
  49 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
  50
  51   //Find clusters V2
  52   SetModule(mod);
  53   FindClustersSSD(fDigits);
  54
  55 }
  56
  57 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
  58   //------------------------------------------------------------
  59   // Actual SSD cluster finder
  60   //------------------------------------------------------------
  61   Int_t smaxall=alldigits->GetEntriesFast();
  62   if (smaxall==0) return;
  63   TObjArray *digits = new TObjArray;
  64   for (Int_t i=0;i<smaxall; i++){
  65     AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
  66     if (d->GetSignal()<3) continue;
  67     digits->AddLast(d);
  68   }
  69   Int_t smax = digits->GetEntriesFast();
  70   if (smax==0) return;
  71
  72   const Int_t kMax=1000;
  73   Int_t np=0, nn=0;
  74   Ali1Dcluster pos[kMax], neg[kMax];
  75   Float_t y=0., q=0., qmax=0.;
  76   Int_t lab[4]={-2,-2,-2,-2};
  77
  78   AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0);
  79   q += d->GetSignal();
  80   y += d->GetCoord2()*d->GetSignal();
  81   qmax=d->GetSignal();
  82   lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
  83   Int_t curr=d->GetCoord2();
  84   Int_t flag=d->GetCoord1();
  85   Int_t *n=&nn;
  86   Ali1Dcluster *c=neg;
  87   Int_t nd=1;
  88   Int_t milab[10];
  89   for (Int_t ilab=0;ilab<10;ilab++){
  90     milab[ilab]=-2;
  91   }
  92   milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
  93
  94   for (Int_t s=1; s<smax; s++) {
  95       d=(AliITSdigitSSD*)digits->UncheckedAt(s);
  96       Int_t strip=d->GetCoord2();
  97       if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
  98          c[*n].SetY(y/q);
  99          c[*n].SetQ(q);
 100          c[*n].SetNd(nd);
 101          CheckLabels2(milab);
 102          c[*n].SetLabels(milab);
 103          //Split suspiciously big cluster
 104          /*
 105          if (nd>10&&nd<16){
 106            c[*n].SetY(y/q-0.3*nd);
 107            c[*n].SetQ(0.5*q);
 108            (*n)++;
 109            if (*n==MAX) {
 110              Error("FindClustersSSD","Too many 1D clusters !");
 111               return;
 112            }
 113            c[*n].SetY(y/q-0.0*nd);
 114            c[*n].SetQ(0.5*q);
 115            c[*n].SetNd(nd);
 116            (*n)++;
 117            if (*n==MAX) {
 118              Error("FindClustersSSD","Too many 1D clusters !");
 119               return;
 120            }
 121            //
 122            c[*n].SetY(y/q+0.3*nd);
 123            c[*n].SetQ(0.5*q);
 124            c[*n].SetNd(nd);
 125            c[*n].SetLabels(milab);
 126          }
 127          else{
 128          */
 129          if (nd>4&&nd<25) {
 130            c[*n].SetY(y/q-0.25*nd);
 131            c[*n].SetQ(0.5*q);
 132            (*n)++;
 133            if (*n==kMax) {
 134              Error("FindClustersSSD","Too many 1D clusters !");
 135              return;
 136            }
 137            c[*n].SetY(y/q+0.25*nd);
 138            c[*n].SetQ(0.5*q);
 139            c[*n].SetNd(nd);
 140            c[*n].SetLabels(milab);
 141          }
 142          (*n)++;
 143          if (*n==kMax) {
 144           Error("FindClustersSSD","Too many 1D clusters !");
 145           return;
 146          }
 147          y=q=qmax=0.;
 148          nd=0;
 149          lab[0]=lab[1]=lab[2]=-2;
 150          //
 151          for (Int_t ilab=0;ilab<10;ilab++){
 152            milab[ilab]=-2;
 153          }
 154          //
 155          if (flag!=d->GetCoord1()) { n=&np; c=pos; }
 156       }
 157       flag=d->GetCoord1();
 158       q += d->GetSignal();
 159       y += d->GetCoord2()*d->GetSignal();
 160       nd++;
 161       if (d->GetSignal()>qmax) {
 162          qmax=d->GetSignal();
 163          lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
 164       }
 165       for (Int_t ilab=0;ilab<10;ilab++) {
 166         if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
 167       }
 168       curr=strip;
 169   }
 170   c[*n].SetY(y/q);
 171   c[*n].SetQ(q);
 172   c[*n].SetNd(nd);
 173   c[*n].SetLabels(lab);
 174   //Split suspiciously big cluster
 175   if (nd>4 && nd<25) {
 176      c[*n].SetY(y/q-0.25*nd);
 177      c[*n].SetQ(0.5*q);
 178      (*n)++;
 179      if (*n==kMax) {
 180         Error("FindClustersSSD","Too many 1D clusters !");
 181         return;
 182      }
 183      c[*n].SetY(y/q+0.25*nd);
 184      c[*n].SetQ(0.5*q);
 185      c[*n].SetNd(nd);
 186      c[*n].SetLabels(lab);
 187   }
 188   (*n)++;
 189   if (*n==kMax) {
 190      Error("FindClustersSSD","Too many 1D clusters !");
 191      return;
 192   }
 193
 194   FindClustersSSD(neg, nn, pos, np);
 195 }
 196
 197
 198 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
 199
 200     //------------------------------------------------------------
 201   // This function creates ITS clusters from raw data
 202   //------------------------------------------------------------
 203   rawReader->Reset();
 204   AliITSRawStreamSSD inputSSD(rawReader);
 205   FindClustersSSD(&inputSSD,clusters);
 206
 207 }
 208
 209 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStream* input,
 210                                         TClonesArray** clusters)
 211 {
 212   //------------------------------------------------------------
 213   // Actual SSD cluster finder for raw data
 214   //------------------------------------------------------------
 215   Int_t nClustersSSD = 0;
 216   const Int_t kMax = 1000;
 217   Ali1Dcluster clusters1D[2][kMax];
 218   Int_t nClusters[2] = {0, 0};
 219   Int_t lab[3]={-2,-2,-2};
 220   Float_t q = 0.;
 221   Float_t y = 0.;
 222   Int_t nDigits = 0;
 223   Int_t prevStrip = -1;
 224   Int_t prevFlag = -1;
 225   Int_t prevModule = -1;
 226
 227   // read raw data input stream
 228   while (kTRUE) {
 229     Bool_t next = input->Next();
 230
 231     if(input->GetSignal()<3 && next) continue;
 232     // check if a new cluster starts
 233     Int_t strip = input->GetCoord2();
 234     Int_t flag = input->GetCoord1();
 235     if ((!next || (input->GetModuleID() != prevModule)||
 236          (strip-prevStrip > 1) || (flag != prevFlag)) &&
 237         (nDigits > 0)) {
 238       if (nClusters[prevFlag] == kMax) {
 239         Error("FindClustersSSD", "Too many 1D clusters !");
 240         return;
 241       }
 242       Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++];
 243       cluster.SetY(y/q);
 244       cluster.SetQ(q);
 245       cluster.SetNd(nDigits);
 246       cluster.SetLabels(lab);
 247
 248       //Split suspiciously big cluster
 249       if (nDigits > 4&&nDigits < 25) {
 250         cluster.SetY(y/q - 0.25*nDigits);
 251         cluster.SetQ(0.5*q);
 252         if (nClusters[prevFlag] == kMax) {
 253           Error("FindClustersSSD", "Too many 1D clusters !");
 254           return;
 255         }
 256         Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++];
 257         cluster2.SetY(y/q + 0.25*nDigits);
 258         cluster2.SetQ(0.5*q);
 259         cluster2.SetNd(nDigits);
 260         cluster2.SetLabels(lab);
 261       }
 262       y = q = 0.;
 263       nDigits = 0;
 264     }
 265
 266     if (!next || (input->GetModuleID() != prevModule)) {
 267       Int_t iModule = prevModule;
 268
 269       // when all data from a module was read, search for clusters
 270       if (prevFlag >= 0) {
 271         clusters[iModule] = new TClonesArray("AliITSRecPoint");
 272         fModule = iModule;
 273         FindClustersSSD(&clusters1D[0][0], nClusters[0],
 274                         &clusters1D[1][0], nClusters[1], clusters[iModule]);
 275         Int_t nClusters = clusters[iModule]->GetEntriesFast();
 276         nClustersSSD += nClusters;
 277       }
 278
 279       if (!next) break;
 280       nClusters[0] = nClusters[1] = 0;
 281       y = q = 0.;
 282       nDigits = 0;
 283     }
 284
 285     // add digit to current cluster
 286     q += input->GetSignal();
 287     y += strip * input->GetSignal();
 288     nDigits++;
 289     prevStrip = strip;
 290     prevFlag = flag;
 291     prevModule = input->GetModuleID();
 292
 293   }
 294
 295   Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
 296 }
 297
 298 void AliITSClusterFinderV2SSD::
 299 FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
 300                 Ali1Dcluster* pos, Int_t np,
 301                 TClonesArray *clusters) {
 302   //------------------------------------------------------------
 303   // Actual SSD cluster finder
 304   //------------------------------------------------------------
 305   TClonesArray &cl=*clusters;
 306   //
 307   Float_t tanp=fTanP, tann=fTanN;
 308   if (fModule>fLastSSD1) {tann=fTanP; tanp=fTanN;}
 309   Int_t idet=fNdet[fModule];
 310   Int_t ncl=0;
 311   //
 312   Int_t negativepair[30000];
 313   Int_t cnegative[3000];
 314   Int_t cused1[3000];
 315   Int_t positivepair[30000];
 316   Int_t cpositive[3000];
 317   Int_t cused2[3000];
 318   for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
 319   for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
 320   for (Int_t i=0;i<30000;i++) {negativepair[i]=0; positivepair[i]=0;}
 321   static Short_t pairs[1000][1000];
 322   memset(pairs,0,sizeof(Short_t)*1000000);
 323 //   Short_t ** pairs = new Short_t*[1000];
 324 //   for (Int_t i=0; i<1000; i++) {
 325 //     pairs[i] = new Short_t[1000];
 326 //     memset(pairs[i],0,sizeof(Short_t)*1000);
 327 //   }
 328   //
 329   // find available pairs
 330   //
 331   for (Int_t i=0; i<np; i++) {
 332     Float_t yp=pos[i].GetY()*fYpitchSSD;
 333     if (pos[i].GetQ()<3) continue;
 334     for (Int_t j=0; j<nn; j++) {
 335       if (neg[j].GetQ()<3) continue;
 336       Float_t yn=neg[j].GetY()*fYpitchSSD;
 337       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 338       Float_t yt=yn + tann*zt;
 339       zt-=fHlSSD; yt-=fHwSSD;
 340       if (TMath::Abs(yt)<fHwSSD+0.01)
 341       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
 342         negativepair[i*10+cnegative[i]] =j;  //index
 343         positivepair[j*10+cpositive[j]] =i;
 344         cnegative[i]++;  //counters
 345         cpositive[j]++;
 346         pairs[i][j]=100;
 347       }
 348     }
 349   }
 350   //
 351   for (Int_t i=0; i<np; i++) {
 352     Float_t yp=pos[i].GetY()*fYpitchSSD;
 353     if (pos[i].GetQ()<3) continue;
 354     for (Int_t j=0; j<nn; j++) {
 355       if (neg[j].GetQ()<3) continue;
 356       if (cpositive[j]&&cnegative[i]) continue;
 357       Float_t yn=neg[j].GetY()*fYpitchSSD;
 358       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 359       Float_t yt=yn + tann*zt;
 360       zt-=fHlSSD; yt-=fHwSSD;
 361       if (TMath::Abs(yt)<fHwSSD+0.1)
 362       if (TMath::Abs(zt)<fHlSSD+0.15) {
 363         if (cnegative[i]==0) pos[i].SetNd(100);  // not available pair
 364         if (cpositive[j]==0) neg[j].SetNd(100);  // not available pair
 365         negativepair[i*10+cnegative[i]] =j;  //index
 366         positivepair[j*10+cpositive[j]] =i;
 367         cnegative[i]++;  //counters
 368         cpositive[j]++;
 369         pairs[i][j]=100;
 370       }
 371     }
 372   }
 373   //
 374   Float_t lp[5];
 375   Int_t milab[10];
 376   Double_t ratio;
 377
 378   //
 379   // sign gold tracks
 380   //
 381   for (Int_t ip=0;ip<np;ip++){
 382     Float_t ybest=1000,zbest=1000,qbest=0;
 383     //
 384     // select gold clusters
 385     if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
 386       Float_t yp=pos[ip].GetY()*fYpitchSSD;
 387       Int_t j = negativepair[10*ip];
 388       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
 389       //
 390       Float_t yn=neg[j].GetY()*fYpitchSSD;
 391       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 392       Float_t yt=yn + tann*zt;
 393       zt-=fHlSSD; yt-=fHwSSD;
 394       ybest=yt; zbest=zt;
 395       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
 396       lp[0]=-(-ybest+fYshift[fModule]);
 397       lp[1]=  -zbest+fZshift[fModule];
 398       lp[2]=0.0025*0.0025;  //SigmaY2
 399       lp[3]=0.110*0.110;  //SigmaZ2
 400
 401       lp[4]=qbest;        //Q
 402       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 403       for (Int_t ilab=0;ilab<3;ilab++){
 404         milab[ilab] = pos[ip].GetLabel(ilab);
 405         milab[ilab+3] = neg[j].GetLabel(ilab);
 406       }
 407       //
 408       CheckLabels2(milab);
 409       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
 410       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
 411       AliITSRecPoint * cl2;
 412       if(clusters){
 413         cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 414         cl2->SetChargeRatio(ratio);
 415         cl2->SetType(1);
 416         pairs[ip][j]=1;
 417         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
 418           cl2->SetType(2);
 419           pairs[ip][j]=2;
 420         }
 421         cused1[ip]++;
 422         cused2[j]++;
 423
 424       }
 425       else{
 426         cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 427         cl2->SetChargeRatio(ratio);
 428         cl2->SetType(1);
 429         pairs[ip][j]=1;
 430         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
 431           cl2->SetType(2);
 432           pairs[ip][j]=2;
 433         }
 434         cused1[ip]++;
 435         cused2[j]++;
 436         fDetTypeRec->AddRecPoint(*cl2);
 437       }
 438       ncl++;
 439     }
 440   }
 441
 442   for (Int_t ip=0;ip<np;ip++){
 443     Float_t ybest=1000,zbest=1000,qbest=0;
 444     //
 445     //
 446     // select "silber" cluster
 447     if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
 448       Int_t in  = negativepair[10*ip];
 449       Int_t ip2 = positivepair[10*in];
 450       if (ip2==ip) ip2 =  positivepair[10*in+1];
 451       Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
 452       if (TMath::Abs(pcharge-neg[in].GetQ())<10){
 453         //
 454         // add first pair
 455         if (pairs[ip][in]==100){  //
 456           Float_t yp=pos[ip].GetY()*fYpitchSSD;
 457           Float_t yn=neg[in].GetY()*fYpitchSSD;
 458           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 459           Float_t yt=yn + tann*zt;
 460           zt-=fHlSSD; yt-=fHwSSD;
 461           ybest =yt;  zbest=zt;
 462           qbest =pos[ip].GetQ();
 463           lp[0]=-(-ybest+fYshift[fModule]);
 464           lp[1]=  -zbest+fZshift[fModule];
 465           lp[2]=0.0025*0.0025;  //SigmaY2
 466           lp[3]=0.110*0.110;  //SigmaZ2
 467
 468           lp[4]=qbest;        //Q
 469           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 470           for (Int_t ilab=0;ilab<3;ilab++){
 471             milab[ilab] = pos[ip].GetLabel(ilab);
 472             milab[ilab+3] = neg[in].GetLabel(ilab);
 473           }
 474           //
 475           CheckLabels2(milab);
 476           ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
 477           milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
 478           Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
 479
 480           AliITSRecPoint * cl2;
 481           if(clusters){
 482             cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 483             cl2->SetChargeRatio(ratio);
 484             cl2->SetType(5);
 485             pairs[ip][in] = 5;
 486             if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
 487               cl2->SetType(6);
 488               pairs[ip][in] = 6;
 489             }
 490           }
 491           else{
 492             cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 493             cl2->SetChargeRatio(ratio);
 494             cl2->SetType(5);
 495             pairs[ip][in] = 5;
 496             if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
 497               cl2->SetType(6);
 498               pairs[ip][in] = 6;
 499             }
 500
 501             fDetTypeRec->AddRecPoint(*cl2);
 502           }
 503           ncl++;
 504         }
 505
 506         //
 507         // add second pair
 508
 509       //        if (!(cused1[ip2] || cused2[in])){  //
 510         if (pairs[ip2][in]==100){
 511           Float_t yp=pos[ip2].GetY()*fYpitchSSD;
 512           Float_t yn=neg[in].GetY()*fYpitchSSD;
 513           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 514           Float_t yt=yn + tann*zt;
 515           zt-=fHlSSD; yt-=fHwSSD;
 516           ybest =yt;  zbest=zt;
 517           qbest =pos[ip2].GetQ();
 518           lp[0]=-(-ybest+fYshift[fModule]);
 519           lp[1]=  -zbest+fZshift[fModule];
 520           lp[2]=0.0025*0.0025;  //SigmaY2
 521           lp[3]=0.110*0.110;  //SigmaZ2
 522
 523           lp[4]=qbest;        //Q
 524           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 525           for (Int_t ilab=0;ilab<3;ilab++){
 526             milab[ilab] = pos[ip2].GetLabel(ilab);
 527             milab[ilab+3] = neg[in].GetLabel(ilab);
 528           }
 529           //
 530           CheckLabels2(milab);
 531           ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
 532           milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
 533           Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
 534
 535           AliITSRecPoint * cl2;
 536           if(clusters){
 537             cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 538             cl2->SetChargeRatio(ratio);
 539             cl2->SetType(5);
 540             pairs[ip2][in] =5;
 541             if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
 542               cl2->SetType(6);
 543               pairs[ip2][in] =6;
 544             }
 545           }
 546           else{
 547             cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 548             cl2->SetChargeRatio(ratio);
 549             cl2->SetType(5);
 550             pairs[ip2][in] =5;
 551             if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
 552               cl2->SetType(6);
 553               pairs[ip2][in] =6;
 554             }
 555
 556             fDetTypeRec->AddRecPoint(*cl2);
 557           }
 558           ncl++;
 559         }
 560         cused1[ip]++;
 561         cused1[ip2]++;
 562         cused2[in]++;
 563       }
 564     }
 565   }
 566
 567   //
 568   for (Int_t jn=0;jn<nn;jn++){
 569     if (cused2[jn]) continue;
 570     Float_t ybest=1000,zbest=1000,qbest=0;
 571     // select "silber" cluster
 572     if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
 573       Int_t ip  = positivepair[10*jn];
 574       Int_t jn2 = negativepair[10*ip];
 575       if (jn2==jn) jn2 =  negativepair[10*ip+1];
 576       Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
 577       //
 578       if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
 579         //
 580         // add first pair
 581         //      if (!(cused1[ip]||cused2[jn])){
 582         if (pairs[ip][jn]==100){
 583           Float_t yn=neg[jn].GetY()*fYpitchSSD;
 584           Float_t yp=pos[ip].GetY()*fYpitchSSD;
 585           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 586           Float_t yt=yn + tann*zt;
 587           zt-=fHlSSD; yt-=fHwSSD;
 588           ybest =yt;  zbest=zt;
 589           qbest =neg[jn].GetQ();
 590           lp[0]=-(-ybest+fYshift[fModule]);
 591           lp[1]=  -zbest+fZshift[fModule];
 592           lp[2]=0.0025*0.0025;  //SigmaY2
 593           lp[3]=0.110*0.110;  //SigmaZ2
 594
 595           lp[4]=qbest;        //Q
 596           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 597           for (Int_t ilab=0;ilab<3;ilab++){
 598             milab[ilab] = pos[ip].GetLabel(ilab);
 599             milab[ilab+3] = neg[jn].GetLabel(ilab);
 600           }
 601           //
 602           CheckLabels2(milab);
 603           ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
 604           milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
 605           Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
 606
 607           AliITSRecPoint * cl2;
 608           if(clusters){
 609             cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 610             cl2->SetChargeRatio(ratio);
 611             cl2->SetType(7);
 612             pairs[ip][jn] =7;
 613             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
 614               cl2->SetType(8);
 615               pairs[ip][jn]=8;
 616             }
 617
 618           }
 619           else{
 620             cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 621             cl2->SetChargeRatio(ratio);
 622             cl2->SetType(7);
 623             pairs[ip][jn] =7;
 624             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
 625               cl2->SetType(8);
 626               pairs[ip][jn]=8;
 627             }
 628
 629             fDetTypeRec->AddRecPoint(*cl2);
 630           }
 631           ncl++;
 632         }
 633         //
 634         // add second pair
 635         //      if (!(cused1[ip]||cused2[jn2])){
 636         if (pairs[ip][jn2]==100){
 637           Float_t yn=neg[jn2].GetY()*fYpitchSSD;
 638           Double_t yp=pos[ip].GetY()*fYpitchSSD;
 639           Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 640           Double_t yt=yn + tann*zt;
 641           zt-=fHlSSD; yt-=fHwSSD;
 642           ybest =yt;  zbest=zt;
 643           qbest =neg[jn2].GetQ();
 644           lp[0]=-(-ybest+fYshift[fModule]);
 645           lp[1]=  -zbest+fZshift[fModule];
 646           lp[2]=0.0025*0.0025;  //SigmaY2
 647           lp[3]=0.110*0.110;  //SigmaZ2
 648
 649           lp[4]=qbest;        //Q
 650           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 651           for (Int_t ilab=0;ilab<3;ilab++){
 652             milab[ilab] = pos[ip].GetLabel(ilab);
 653             milab[ilab+3] = neg[jn2].GetLabel(ilab);
 654           }
 655           //
 656           CheckLabels2(milab);
 657           ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
 658           milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
 659           Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
 660           AliITSRecPoint * cl2;
 661           if(clusters){
 662             cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 663             cl2->SetChargeRatio(ratio);
 664             pairs[ip][jn2]=7;
 665             cl2->SetType(7);
 666             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
 667               cl2->SetType(8);
 668               pairs[ip][jn2]=8;
 669             }
 670
 671           }
 672           else{
 673             cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 674             cl2->SetChargeRatio(ratio);
 675             pairs[ip][jn2]=7;
 676             cl2->SetType(7);
 677             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
 678               cl2->SetType(8);
 679               pairs[ip][jn2]=8;
 680             }
 681
 682             fDetTypeRec->AddRecPoint(*cl2);
 683           }
 684
 685           ncl++;
 686         }
 687         cused1[ip]++;
 688         cused2[jn]++;
 689         cused2[jn2]++;
 690       }
 691     }
 692   }
 693
 694   for (Int_t ip=0;ip<np;ip++){
 695     Float_t ybest=1000,zbest=1000,qbest=0;
 696     //
 697     // 2x2 clusters
 698     //
 699     if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
 700       Float_t minchargediff =4.;
 701       Int_t j=-1;
 702       for (Int_t di=0;di<cnegative[ip];di++){
 703         Int_t   jc = negativepair[ip*10+di];
 704         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
 705         if (TMath::Abs(chargedif)<minchargediff){
 706           j =jc;
 707           minchargediff = TMath::Abs(chargedif);
 708         }
 709       }
 710       if (j<0) continue;  // not proper cluster
 711       Int_t count =0;
 712       for (Int_t di=0;di<cnegative[ip];di++){
 713         Int_t   jc = negativepair[ip*10+di];
 714         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
 715         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
 716       }
 717       if (count>1) continue;  // more than one "proper" cluster for positive
 718       //
 719       count =0;
 720       for (Int_t dj=0;dj<cpositive[j];dj++){
 721         Int_t   ic  = positivepair[j*10+dj];
 722         Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
 723         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
 724       }
 725       if (count>1) continue;  // more than one "proper" cluster for negative
 726
 727       Int_t jp = 0;
 728
 729       count =0;
 730       for (Int_t dj=0;dj<cnegative[jp];dj++){
 731         Int_t   ic = positivepair[jp*10+dj];
 732         Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
 733         if (TMath::Abs(chargedif)<minchargediff+4.) count++;
 734       }
 735       if (count>1) continue;
 736       if (pairs[ip][j]<100) continue;
 737       //
 738       //almost gold clusters
 739       Float_t yp=pos[ip].GetY()*fYpitchSSD;
 740       Float_t yn=neg[j].GetY()*fYpitchSSD;
 741       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 742       Float_t yt=yn + tann*zt;
 743       zt-=fHlSSD; yt-=fHwSSD;
 744       ybest=yt; zbest=zt;
 745       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
 746       lp[0]=-(-ybest+fYshift[fModule]);
 747       lp[1]=  -zbest+fZshift[fModule];
 748       lp[2]=0.0025*0.0025;  //SigmaY2
 749       lp[3]=0.110*0.110;  //SigmaZ2
 750       lp[4]=qbest;        //Q
 751       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 752       for (Int_t ilab=0;ilab<3;ilab++){
 753         milab[ilab] = pos[ip].GetLabel(ilab);
 754         milab[ilab+3] = neg[j].GetLabel(ilab);
 755       }
 756       //
 757       CheckLabels2(milab);
 758       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
 759       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
 760       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
 761       AliITSRecPoint * cl2;
 762       if(clusters){
 763         cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 764         cl2->SetChargeRatio(ratio);
 765         cl2->SetType(10);
 766         pairs[ip][j]=10;
 767         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
 768           cl2->SetType(11);
 769           pairs[ip][j]=11;
 770         }
 771         cused1[ip]++;
 772         cused2[j]++;
 773       }
 774       else{
 775         cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 776         cl2->SetChargeRatio(ratio);
 777         cl2->SetType(10);
 778         pairs[ip][j]=10;
 779         if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
 780           cl2->SetType(11);
 781           pairs[ip][j]=11;
 782         }
 783         cused1[ip]++;
 784         cused2[j]++;
 785
 786         fDetTypeRec->AddRecPoint(*cl2);
 787       }
 788       ncl++;
 789     }
 790
 791   }
 792
 793   //
 794   for (Int_t i=0; i<np; i++) {
 795     Float_t ybest=1000,zbest=1000,qbest=0;
 796     Float_t yp=pos[i].GetY()*fYpitchSSD;
 797     if (pos[i].GetQ()<3) continue;
 798     for (Int_t j=0; j<nn; j++) {
 799     //    for (Int_t di = 0;di<cpositive[i];di++){
 800     //  Int_t j = negativepair[10*i+di];
 801       if (neg[j].GetQ()<3) continue;
 802       if (cused2[j]||cused1[i]) continue;
 803       if (pairs[i][j]>0 &&pairs[i][j]<100) continue;
 804       ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
 805       Float_t yn=neg[j].GetY()*fYpitchSSD;
 806       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 807       Float_t yt=yn + tann*zt;
 808       zt-=fHlSSD; yt-=fHwSSD;
 809       if (TMath::Abs(yt)<fHwSSD+0.01)
 810       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
 811         ybest=yt; zbest=zt;
 812         qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
 813         lp[0]=-(-ybest+fYshift[fModule]);
 814         lp[1]=  -zbest+fZshift[fModule];
 815         lp[2]=0.0025*0.0025;  //SigmaY2
 816         lp[3]=0.110*0.110;  //SigmaZ2
 817
 818         lp[4]=qbest;        //Q
 819         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 820         for (Int_t ilab=0;ilab<3;ilab++){
 821           milab[ilab] = pos[i].GetLabel(ilab);
 822           milab[ilab+3] = neg[j].GetLabel(ilab);
 823         }
 824         //
 825         CheckLabels2(milab);
 826         milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
 827         Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
 828         AliITSRecPoint * cl2;
 829         if(clusters){
 830           cl2 = new (cl[ncl]) AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 831           cl2->SetChargeRatio(ratio);
 832           cl2->SetType(100+cpositive[j]+cnegative[i]);
 833         }
 834         else{
 835           cl2 = new AliITSRecPoint(fModule,fDetTypeRec->GetITSgeom(),milab,lp,info);
 836           cl2->SetChargeRatio(ratio);
 837           cl2->SetType(100+cpositive[j]+cnegative[i]);
 838           fDetTypeRec->AddRecPoint(*cl2);
 839         }
 840         ncl++;
 841         //cl2->SetType(0);
 842         /*
 843           if (pairs[i][j]<100){
 844           printf("problem:- %d\n", pairs[i][j]);
 845           }
 846           if (cnegative[i]<2&&cpositive[j]<2){
 847           printf("problem:- %d\n", pairs[i][j]);
 848           }
 849         */
 850       }
 851     }
 852   }
 853
 854 //   for (Int_t i=0; i<1000; i++) delete [] pairs[i];
 855 //   delete [] pairs;
 856
 857 }
 858
 859
 860