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