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 "AliITSclusterV2.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
  36   //Default constructor
  37
  38   fLastSSD1=fDetTypeRec->GetITSgeom()->GetModuleIndex(6,1,1)-1;
  39   fYpitchSSD=0.0095;
  40   fHwSSD=3.65;
  41   fHlSSD=2.00;
  42   fTanP=0.0275;
  43   fTanN=0.0075;
  44
  45
  46
  47 }
  48
  49
  50 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
  51
  52   //Find clusters V2
  53   SetModule(mod);
  54   FindClustersSSD(fDigits);
  55
  56 }
  57
  58 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
  59   //------------------------------------------------------------
  60   // Actual SSD cluster finder
  61   //------------------------------------------------------------
  62   Int_t smaxall=alldigits->GetEntriesFast();
  63   if (smaxall==0) return;
  64   TObjArray *digits = new TObjArray;
  65   for (Int_t i=0;i<smaxall; i++){
  66     AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
  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 && 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("AliITSclusterV2");
 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();
 288     y += strip * input->GetSignal();
 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   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       AliITSclusterV2 * cl2;
 412       if(clusters) cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
 413       else{
 414         cl2 = new AliITSclusterV2(milab,lp,info);
 415         fDetTypeRec->AddClusterV2(*cl2);
 416       }
 417       ncl++;
 418       cl2->SetChargeRatio(ratio);
 419       cl2->SetType(1);
 420       pairs[ip][j]=1;
 421       if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
 422         cl2->SetType(2);
 423         pairs[ip][j]=2;
 424       }
 425       cused1[ip]++;
 426       cused2[j]++;
 427     }
 428   }
 429
 430   for (Int_t ip=0;ip<np;ip++){
 431     Float_t ybest=1000,zbest=1000,qbest=0;
 432     //
 433     //
 434     // select "silber" cluster
 435     if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
 436       Int_t in  = negativepair[10*ip];
 437       Int_t ip2 = positivepair[10*in];
 438       if (ip2==ip) ip2 =  positivepair[10*in+1];
 439       Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
 440       if (TMath::Abs(pcharge-neg[in].GetQ())<10){
 441         //
 442         // add first pair
 443         if (pairs[ip][in]==100){  //
 444           Float_t yp=pos[ip].GetY()*fYpitchSSD;
 445           Float_t yn=neg[in].GetY()*fYpitchSSD;
 446           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 447           Float_t yt=yn + tann*zt;
 448           zt-=fHlSSD; yt-=fHwSSD;
 449           ybest =yt;  zbest=zt;
 450           qbest =pos[ip].GetQ();
 451           lp[0]=-(-ybest+fYshift[fModule]);
 452           lp[1]=  -zbest+fZshift[fModule];
 453           lp[2]=0.0025*0.0025;  //SigmaY2
 454           lp[3]=0.110*0.110;  //SigmaZ2
 455
 456           lp[4]=qbest;        //Q
 457           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 458           for (Int_t ilab=0;ilab<3;ilab++){
 459             milab[ilab] = pos[ip].GetLabel(ilab);
 460             milab[ilab+3] = neg[in].GetLabel(ilab);
 461           }
 462           //
 463           CheckLabels2(milab);
 464           ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
 465           milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
 466           Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
 467
 468           AliITSclusterV2 * cl2;
 469           if(clusters) cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
 470           else{
 471             cl2 = new AliITSclusterV2(milab,lp,info);
 472             fDetTypeRec->AddClusterV2(*cl2);
 473           }
 474           ncl++;
 475           cl2->SetChargeRatio(ratio);
 476           cl2->SetType(5);
 477           pairs[ip][in] = 5;
 478           if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
 479             cl2->SetType(6);
 480             pairs[ip][in] = 6;
 481           }
 482         }
 483         //
 484         // add second pair
 485
 486         //      if (!(cused1[ip2] || cused2[in])){  //
 487         if (pairs[ip2][in]==100){
 488           Float_t yp=pos[ip2].GetY()*fYpitchSSD;
 489           Float_t yn=neg[in].GetY()*fYpitchSSD;
 490           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 491           Float_t yt=yn + tann*zt;
 492           zt-=fHlSSD; yt-=fHwSSD;
 493           ybest =yt;  zbest=zt;
 494           qbest =pos[ip2].GetQ();
 495           lp[0]=-(-ybest+fYshift[fModule]);
 496           lp[1]=  -zbest+fZshift[fModule];
 497           lp[2]=0.0025*0.0025;  //SigmaY2
 498           lp[3]=0.110*0.110;  //SigmaZ2
 499
 500           lp[4]=qbest;        //Q
 501           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 502           for (Int_t ilab=0;ilab<3;ilab++){
 503             milab[ilab] = pos[ip2].GetLabel(ilab);
 504             milab[ilab+3] = neg[in].GetLabel(ilab);
 505           }
 506           //
 507           CheckLabels2(milab);
 508           ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
 509           milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
 510           Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
 511
 512           AliITSclusterV2 * cl2;
 513           if(clusters) cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
 514           else{
 515             cl2 = new AliITSclusterV2(milab,lp,info);
 516             fDetTypeRec->AddClusterV2(*cl2);
 517           }
 518           ncl++;
 519           cl2->SetChargeRatio(ratio);
 520           cl2->SetType(5);
 521           pairs[ip2][in] =5;
 522           if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
 523             cl2->SetType(6);
 524             pairs[ip2][in] =6;
 525           }
 526         }
 527         cused1[ip]++;
 528         cused1[ip2]++;
 529         cused2[in]++;
 530       }
 531     }
 532   }
 533
 534   //
 535   for (Int_t jn=0;jn<nn;jn++){
 536     if (cused2[jn]) continue;
 537     Float_t ybest=1000,zbest=1000,qbest=0;
 538     // select "silber" cluster
 539     if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
 540       Int_t ip  = positivepair[10*jn];
 541       Int_t jn2 = negativepair[10*ip];
 542       if (jn2==jn) jn2 =  negativepair[10*ip+1];
 543       Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
 544       //
 545       if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
 546         //
 547         // add first pair
 548         //      if (!(cused1[ip]||cused2[jn])){
 549         if (pairs[ip][jn]==100){
 550           Float_t yn=neg[jn].GetY()*fYpitchSSD;
 551           Float_t yp=pos[ip].GetY()*fYpitchSSD;
 552           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 553           Float_t yt=yn + tann*zt;
 554           zt-=fHlSSD; yt-=fHwSSD;
 555           ybest =yt;  zbest=zt;
 556           qbest =neg[jn].GetQ();
 557           lp[0]=-(-ybest+fYshift[fModule]);
 558           lp[1]=  -zbest+fZshift[fModule];
 559           lp[2]=0.0025*0.0025;  //SigmaY2
 560           lp[3]=0.110*0.110;  //SigmaZ2
 561
 562           lp[4]=qbest;        //Q
 563           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 564           for (Int_t ilab=0;ilab<3;ilab++){
 565             milab[ilab] = pos[ip].GetLabel(ilab);
 566             milab[ilab+3] = neg[jn].GetLabel(ilab);
 567           }
 568           //
 569           CheckLabels2(milab);
 570           ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
 571           milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
 572           Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
 573
 574           AliITSclusterV2 * cl2;
 575           if(clusters) cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
 576           else{
 577             cl2 = new AliITSclusterV2(milab,lp,info);
 578             fDetTypeRec->AddClusterV2(*cl2);
 579           }
 580           ncl++;
 581           cl2->SetChargeRatio(ratio);
 582           cl2->SetType(7);
 583           pairs[ip][jn] =7;
 584           if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
 585             cl2->SetType(8);
 586             pairs[ip][jn]=8;
 587           }
 588         }
 589         //
 590         // add second pair
 591         //      if (!(cused1[ip]||cused2[jn2])){
 592         if (pairs[ip][jn2]==100){
 593           Float_t yn=neg[jn2].GetY()*fYpitchSSD;
 594           Double_t yp=pos[ip].GetY()*fYpitchSSD;
 595           Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 596           Double_t yt=yn + tann*zt;
 597           zt-=fHlSSD; yt-=fHwSSD;
 598           ybest =yt;  zbest=zt;
 599           qbest =neg[jn2].GetQ();
 600           lp[0]=-(-ybest+fYshift[fModule]);
 601           lp[1]=  -zbest+fZshift[fModule];
 602           lp[2]=0.0025*0.0025;  //SigmaY2
 603           lp[3]=0.110*0.110;  //SigmaZ2
 604
 605           lp[4]=qbest;        //Q
 606           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 607           for (Int_t ilab=0;ilab<3;ilab++){
 608             milab[ilab] = pos[ip].GetLabel(ilab);
 609             milab[ilab+3] = neg[jn2].GetLabel(ilab);
 610           }
 611           //
 612           CheckLabels2(milab);
 613           ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
 614           milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
 615           Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
 616           AliITSclusterV2 * cl2;
 617           if(clusters) cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
 618           else{
 619             cl2 = new AliITSclusterV2(milab,lp,info);
 620             fDetTypeRec->AddClusterV2(*cl2);
 621           }
 622
 623           ncl++;
 624           cl2->SetChargeRatio(ratio);
 625           pairs[ip][jn2]=7;
 626           cl2->SetType(7);
 627           if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
 628             cl2->SetType(8);
 629             pairs[ip][jn2]=8;
 630           }
 631         }
 632         cused1[ip]++;
 633         cused2[jn]++;
 634         cused2[jn2]++;
 635       }
 636     }
 637   }
 638
 639   for (Int_t ip=0;ip<np;ip++){
 640     Float_t ybest=1000,zbest=1000,qbest=0;
 641     //
 642     // 2x2 clusters
 643     //
 644     if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
 645       Float_t minchargediff =4.;
 646       Int_t j=-1;
 647       for (Int_t di=0;di<cnegative[ip];di++){
 648         Int_t   jc = negativepair[ip*10+di];
 649         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
 650         if (TMath::Abs(chargedif)<minchargediff){
 651           j =jc;
 652           minchargediff = TMath::Abs(chargedif);
 653         }
 654       }
 655       if (j<0) continue;  // not proper cluster
 656       Int_t count =0;
 657       for (Int_t di=0;di<cnegative[ip];di++){
 658         Int_t   jc = negativepair[ip*10+di];
 659         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
 660         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
 661       }
 662       if (count>1) continue;  // more than one "proper" cluster for positive
 663       //
 664       count =0;
 665       for (Int_t dj=0;dj<cpositive[j];dj++){
 666         Int_t   ic  = positivepair[j*10+dj];
 667         Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
 668         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
 669       }
 670       if (count>1) continue;  // more than one "proper" cluster for negative
 671
 672       Int_t jp = 0;
 673
 674       count =0;
 675       for (Int_t dj=0;dj<cnegative[jp];dj++){
 676         Int_t   ic = positivepair[jp*10+dj];
 677         Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
 678         if (TMath::Abs(chargedif)<minchargediff+4.) count++;
 679       }
 680       if (count>1) continue;
 681       if (pairs[ip][j]<100) continue;
 682       //
 683       //almost gold clusters
 684       Float_t yp=pos[ip].GetY()*fYpitchSSD;
 685       Float_t yn=neg[j].GetY()*fYpitchSSD;
 686       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 687       Float_t yt=yn + tann*zt;
 688       zt-=fHlSSD; yt-=fHwSSD;
 689       ybest=yt; zbest=zt;
 690       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
 691       lp[0]=-(-ybest+fYshift[fModule]);
 692       lp[1]=  -zbest+fZshift[fModule];
 693       lp[2]=0.0025*0.0025;  //SigmaY2
 694       lp[3]=0.110*0.110;  //SigmaZ2
 695       lp[4]=qbest;        //Q
 696       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 697       for (Int_t ilab=0;ilab<3;ilab++){
 698         milab[ilab] = pos[ip].GetLabel(ilab);
 699         milab[ilab+3] = neg[j].GetLabel(ilab);
 700       }
 701       //
 702       CheckLabels2(milab);
 703       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
 704       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
 705       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
 706       AliITSclusterV2 * cl2;
 707       if(clusters) cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
 708       else{
 709         cl2 = new AliITSclusterV2(milab,lp,info);
 710             fDetTypeRec->AddClusterV2(*cl2);
 711       }
 712
 713
 714       ncl++;
 715       cl2->SetChargeRatio(ratio);
 716       cl2->SetType(10);
 717       pairs[ip][j]=10;
 718       if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
 719         cl2->SetType(11);
 720         pairs[ip][j]=11;
 721       }
 722       cused1[ip]++;
 723       cused2[j]++;
 724     }
 725
 726   }
 727
 728   //
 729   for (Int_t i=0; i<np; i++) {
 730     Float_t ybest=1000,zbest=1000,qbest=0;
 731     Float_t yp=pos[i].GetY()*fYpitchSSD;
 732     if (pos[i].GetQ()<3) continue;
 733     for (Int_t j=0; j<nn; j++) {
 734     //    for (Int_t di = 0;di<cpositive[i];di++){
 735     //  Int_t j = negativepair[10*i+di];
 736       if (neg[j].GetQ()<3) continue;
 737       if (cused2[j]||cused1[i]) continue;
 738       if (pairs[i][j]>0 &&pairs[i][j]<100) continue;
 739       ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
 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       if (TMath::Abs(yt)<fHwSSD+0.01)
 745       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
 746         ybest=yt; zbest=zt;
 747         qbest=0.5*(pos[i].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
 753         lp[4]=qbest;        //Q
 754         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 755         for (Int_t ilab=0;ilab<3;ilab++){
 756           milab[ilab] = pos[i].GetLabel(ilab);
 757           milab[ilab+3] = neg[j].GetLabel(ilab);
 758         }
 759         //
 760         CheckLabels2(milab);
 761         milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
 762         Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
 763         AliITSclusterV2 * cl2;
 764         if(clusters) cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
 765         else{
 766           cl2 = new AliITSclusterV2(milab,lp,info);
 767           fDetTypeRec->AddClusterV2(*cl2);
 768         }
 769         ncl++;
 770         cl2->SetChargeRatio(ratio);
 771         cl2->SetType(100+cpositive[j]+cnegative[i]);
 772         //cl2->SetType(0);
 773         /*
 774           if (pairs[i][j]<100){
 775           printf("problem:- %d\n", pairs[i][j]);
 776           }
 777           if (cnegative[i]<2&&cpositive[j]<2){
 778           printf("problem:- %d\n", pairs[i][j]);
 779           }
 780         */
 781       }
 782     }
 783   }
 784
 785 //   for (Int_t i=0; i<1000; i++) delete [] pairs[i];
 786 //   delete [] pairs;
 787
 788 }
 789
 790
 791