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