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