ITS/AliITSClusterFinderV2SSD.cxx

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