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