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