ITS/AliITSclustererV2.cxx

   1 //-------------------------------------------------------------------------
   2 //            Implementation of the ITS clusterer V2 class
   3 //
   4 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
   5 //-------------------------------------------------------------------------
   6
   7
   8 #include "AliRun.h"
   9
  10 #include "AliITSclustererV2.h"
  11 #include "AliITSclusterV2.h"
  12 #include "AliITSDetTypeRec.h"
  13 #include "AliRawReader.h"
  14 #include "AliITSRawStreamSPD.h"
  15 #include "AliITSRawStreamSDD.h"
  16 #include "AliITSRawStreamSSD.h"
  17
  18 #include <TFile.h>
  19 #include <TTree.h>
  20 #include <TClonesArray.h>
  21 #include "AliITSgeomTGeo.h"
  22 #include "AliITSdigitSPD.h"
  23 #include "AliITSdigitSDD.h"
  24 #include "AliITSdigitSSD.h"
  25 #include "AliMC.h"
  26
  27 ClassImp(AliITSclustererV2)
  28
  29 extern AliRun *gAlice;
  30
  31 AliITSclustererV2::AliITSclustererV2():
  32 fNModules(0),
  33 fEvent(0),
  34 fI(0),
  35 fLastSPD1(0),
  36 fNySPD(0),
  37 fNzSPD(0),
  38 fYpitchSPD(0),
  39 fZ1pitchSPD(0),
  40 fZ2pitchSPD(0),
  41 fHwSPD(0),
  42 fHlSPD(0),
  43 fNySDD(0),
  44 fNzSDD(0),
  45 fYpitchSDD(0),
  46 fZpitchSDD(0),
  47 fHwSDD(0),
  48 fHlSDD(0),
  49 fYoffSDD(0),
  50 fLastSSD1(0),
  51 fYpitchSSD(0),
  52 fHwSSD(0),
  53 fHlSSD(0),
  54 fTanP(0),
  55 fTanN(0){
  56    //default constructor
  57   for(Int_t i=0;i<260;i++)fYSPD[i]=0.;
  58   for(Int_t i=0;i<170;i++)fZSPD[i]=0.;
  59   for(Int_t i=0;i<2200;i++){
  60     fYshift[i]=0.;
  61     fZshift[i]=0.;
  62     fNdet[i]=0;
  63     fNlayer[i]=0;
  64   }
  65
  66 }
  67 AliITSclustererV2::AliITSclustererV2(const Char_t *geom):
  68 fNModules(AliITSgeomTGeo::GetNModules()),
  69 fEvent(0),
  70 fI(0),
  71 fLastSPD1(AliITSgeomTGeo::GetModuleIndex(2,1,1)-1),
  72 fNySPD(256),
  73 fNzSPD(160),
  74 fYpitchSPD(0.0050),
  75 fZ1pitchSPD(0.0425),
  76 fZ2pitchSPD(0.0625),
  77 fHwSPD(0.64),
  78 fHlSPD(3.48),
  79 fNySDD(256),
  80 fNzSDD(256),
  81 fYpitchSDD(0.01825),
  82 fZpitchSDD(0.02940),
  83 fHwSDD(3.5085),
  84 fHlSDD(3.7632),
  85 fYoffSDD(0.0425),
  86 fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1),
  87 fYpitchSSD(0.0095),
  88 fHwSSD(3.65),
  89 fHlSSD(2.00),
  90 fTanP(0.0275),
  91 fTanN(0.0075) {
  92   //------------------------------------------------------------
  93   // Standard constructor
  94   //------------------------------------------------------------
  95   if (geom) {
  96     AliWarning("\"geom\" is actually a dummy argument !");
  97   }
  98   Int_t m;
  99   for (m=0; m<fNModules; m++) {
 100      Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
 101      const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(m);
 102      fYshift[m] = (tm->Inverse()).GetTranslation()[1];
 103      fZshift[m] = (tm->Inverse()).GetTranslation()[2];
 104      fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
 105      fNlayer[m] = lay-1;
 106   }
 107
 108   //SPD geometry
 109   fYSPD[0]=0.5*fYpitchSPD;
 110   for (m=1; m<fNySPD; m++) fYSPD[m]=fYSPD[m-1]+fYpitchSPD;
 111   fZSPD[0]=fZ1pitchSPD;
 112   for (m=1; m<fNzSPD; m++) {
 113     Double_t dz=fZ1pitchSPD;
 114     if (m==31 || m==32 || m==63  || m==64  || m==95 || m==96 ||
 115         m==127 || m==128) dz=fZ2pitchSPD;
 116     fZSPD[m]=fZSPD[m-1]+dz;
 117   }
 118   for (m=0; m<fNzSPD; m++) {
 119     Double_t dz=0.5*fZ1pitchSPD;
 120     if (m==31 || m==32 || m==63  || m==64  || m==95 || m==96 ||
 121         m==127 || m==128) dz=0.5*fZ2pitchSPD;
 122     fZSPD[m]-=dz;
 123   }
 124
 125 }
 126
 127
 128 Int_t AliITSclustererV2::Digits2Clusters(TTree *dTree, TTree *cTree) {
 129   //------------------------------------------------------------
 130   // This function creates ITS clusters
 131   //------------------------------------------------------------
 132   Int_t ncl=0;
 133
 134   if (!dTree) {
 135     Error("Digits2Clusters","Can't get the tree with digits !");
 136     return 1;
 137   }
 138
 139   TClonesArray *digitsSPD=new TClonesArray("AliITSdigitSPD",3000);
 140   dTree->SetBranchAddress("ITSDigitsSPD",&digitsSPD);
 141   TClonesArray *digitsSDD=new TClonesArray("AliITSdigitSDD",3000);
 142   dTree->SetBranchAddress("ITSDigitsSDD",&digitsSDD);
 143   TClonesArray *digitsSSD=new TClonesArray("AliITSdigitSSD",3000);
 144   dTree->SetBranchAddress("ITSDigitsSSD",&digitsSSD);
 145
 146   TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
 147   TBranch *branch=cTree->GetBranch("Clusters");
 148   if (!branch) cTree->Branch("Clusters",&clusters);
 149   else branch->SetAddress(&clusters);
 150
 151   Int_t mmax=(Int_t)dTree->GetEntries();
 152   if (mmax!=fNModules) {
 153     Error("Digits2Clusters","Number of entries != number of modules !");
 154     return 1;
 155   }
 156
 157   for (fI=0; fI<mmax; fI++) {
 158     dTree->GetEvent(fI);
 159
 160     if     (digitsSPD->GetEntriesFast()!=0)
 161       FindClustersSPD(digitsSPD,clusters);
 162     else
 163       if(digitsSDD->GetEntriesFast()!=0)
 164         FindClustersSDD(digitsSDD,clusters);
 165       else if(digitsSSD->GetEntriesFast()!=0)
 166         FindClustersSSD(digitsSSD,clusters);
 167
 168     ncl+=clusters->GetEntriesFast();
 169
 170     cTree->Fill();
 171
 172     digitsSPD->Clear();
 173     digitsSDD->Clear();
 174     digitsSSD->Clear();
 175     clusters->Clear();
 176   }
 177
 178   //cTree->Write();
 179
 180   delete clusters;
 181
 182   delete digitsSPD;
 183   delete digitsSDD;
 184   delete digitsSSD;
 185
 186   //delete dTree;
 187
 188   Info("Digits2Clusters","Number of found clusters : %d",ncl);
 189
 190   return 0;
 191 }
 192
 193 void AliITSclustererV2::Digits2Clusters(AliRawReader* rawReader) {
 194   //------------------------------------------------------------
 195   // This function creates ITS clusters from raw data
 196   //------------------------------------------------------------
 197   AliRunLoader* runLoader = AliRunLoader::Instance();
 198   if (!runLoader) {
 199     Error("Digits2Clusters", "no run loader found");
 200     return;
 201   }
 202   runLoader->LoadKinematics();
 203   AliLoader* itsLoader = runLoader->GetLoader("ITSLoader");
 204   if (!itsLoader) {
 205     Error("Digits2Clusters", "no loader for ITS found");
 206     return;
 207   }
 208   if (!itsLoader->TreeR()) itsLoader->MakeTree("R");
 209   TTree* cTree = itsLoader->TreeR();
 210
 211   TClonesArray *array=new TClonesArray("AliITSclusterV2",1000);
 212   cTree->Branch("Clusters",&array);
 213   delete array;
 214
 215   TClonesArray** clusters = new TClonesArray*[fNModules];
 216   for (Int_t iModule = 0; iModule < fNModules; iModule++) {
 217     clusters[iModule] = NULL;
 218   }
 219   // one TClonesArray per module
 220
 221   rawReader->Reset();
 222   AliITSRawStreamSPD inputSPD(rawReader);
 223   FindClustersSPD(&inputSPD, clusters);
 224
 225   rawReader->Reset();
 226   AliITSRawStreamSDD inputSDD(rawReader);
 227   FindClustersSDD(&inputSDD, clusters);
 228
 229   rawReader->Reset();
 230   AliITSRawStreamSSD inputSSD(rawReader);
 231   FindClustersSSD(&inputSSD, clusters);
 232
 233   // write all clusters to the tree
 234   Int_t nClusters = 0;
 235   TClonesArray *emptyArray=new TClonesArray("AliITSclusterV2");
 236   for (Int_t iModule = 0; iModule < fNModules; iModule++) {
 237     array = clusters[iModule];
 238     if (!array) {
 239       Error("Digits2Clusters", "data for module %d missing!", iModule);
 240       array = emptyArray;
 241     }
 242     cTree->SetBranchAddress("Clusters", &array);
 243     cTree->Fill();
 244     nClusters += array->GetEntriesFast();
 245   }
 246   delete emptyArray;
 247
 248   itsLoader->WriteRecPoints("OVERWRITE");
 249
 250   delete[] clusters;
 251
 252   Info("Digits2Clusters", "total number of found clusters in ITS: %d\n",
 253        nClusters);
 254 }
 255
 256 //**** Fast clusters *******************************
 257 #include "TParticle.h"
 258
 259 //#include "AliITS.h"
 260 #include "AliITSmodule.h"
 261 #include "AliITSRecPoint.h"
 262 #include "AliITSsimulationFastPoints.h"
 263 #include "AliITSRecPoint.h"
 264
 265
 266 static void CheckLabels(Int_t lab[3]) {
 267   //------------------------------------------------------------
 268   // Tries to find mother's labels
 269   //------------------------------------------------------------
 270   AliRunLoader *rl = AliRunLoader::Instance();
 271   TTree *trK=(TTree*)rl->TreeK();
 272
 273   if(trK){
 274     Int_t label=lab[0];
 275     if (label>=0) {
 276       TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 277       label=-3;
 278       while (part->P() < 0.005) {
 279         Int_t m=part->GetFirstMother();
 280         if (m<0) {
 281           Info("CheckLabels","Primary momentum: %f",part->P());
 282           break;
 283         }
 284         if (part->GetStatusCode()>0) {
 285           Info("CheckLabels","Primary momentum: %f",part->P());
 286           break;
 287         }
 288         label=m;
 289         part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 290       }
 291       if(lab[1]<0){
 292         lab[1]=label;
 293       }
 294       else if (lab[2]<0) {
 295         lab[2]=label;
 296       }
 297       else {
 298         //      cerr<<"CheckLabels : No empty labels !\n";
 299       }
 300     }
 301   }
 302 }
 303
 304 /*
 305 static void CheckLabels(Int_t lab[3]) {
 306   //------------------------------------------------------------
 307   // Tries to find mother's labels
 308   //------------------------------------------------------------
 309
 310   if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit
 311
 312   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 313   for (Int_t i=0;i<3;i++){
 314     Int_t label = lab[i];
 315     if (label>=0 && label<ntracks) {
 316       TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 317       if (part->P() < 0.005) {
 318         Int_t m=part->GetFirstMother();
 319         if (m<0) {
 320           continue;
 321         }
 322         if (part->GetStatusCode()>0) {
 323           continue;
 324         }
 325         lab[i]=m;
 326       }
 327     }
 328   }
 329
 330 }
 331 */
 332 static void CheckLabels2(Int_t lab[10]) {
 333   //------------------------------------------------------------
 334   // Tries to find mother's labels
 335   //------------------------------------------------------------
 336   AliRunLoader *rl = AliRunLoader::Instance();
 337   TTree *trK=(TTree*)rl->TreeK();
 338
 339   if(trK){
 340     Int_t nlabels =0;
 341     for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
 342     if(nlabels == 0) return; // In case of no labels just exit
 343
 344     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 345
 346     for (Int_t i=0;i<10;i++){
 347       Int_t label = lab[i];
 348       if (label>=0 && label<ntracks) {
 349         TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 350         if (part->P() < 0.02) {
 351           Int_t m=part->GetFirstMother();
 352           if (m<0) {
 353             continue;
 354           }
 355           if (part->GetStatusCode()>0) {
 356             continue;
 357           }
 358           lab[i]=m;
 359         }
 360         else
 361           if (part->P() < 0.12 && nlabels>3) {
 362             lab[i]=-2;
 363             nlabels--;
 364           }
 365       }
 366       else{
 367         if ( (label>ntracks||label <0) && nlabels>3) {
 368           lab[i]=-2;
 369           nlabels--;
 370         }
 371       }
 372     }
 373     if (nlabels>3){
 374       for (Int_t i=0;i<10;i++){
 375         if (nlabels>3){
 376           Int_t label = lab[i];
 377           if (label>=0 && label<ntracks) {
 378             TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 379             if (part->P() < 0.1) {
 380               lab[i]=-2;
 381               nlabels--;
 382             }
 383           }
 384         }
 385       }
 386     }
 387
 388     //compress labels -- if multi-times the same
 389     Int_t lab2[10];
 390     for (Int_t i=0;i<10;i++) lab2[i]=-2;
 391     for (Int_t i=0;i<10  ;i++){
 392       if (lab[i]<0) continue;
 393       for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
 394         if (lab2[j]<0) {
 395           lab2[j]= lab[i];
 396           break;
 397         }
 398       }
 399     }
 400     for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
 401   }
 402 }
 403
 404 static void AddLabel(Int_t lab[10], Int_t label) {
 405 // add label of the particle in the kine tree which originates this cluster
 406 // (used for reconstruction of MC data only, for comparison purposes)
 407   AliRunLoader *rl = AliRunLoader::Instance();
 408   TTree *trK=(TTree*)rl->TreeK();
 409
 410   if(trK){
 411     if(label<0) return; // In case of no label just exit
 412
 413     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 414     if (label>ntracks) return;
 415     for (Int_t i=0;i<10;i++){
 416       //    if (label<0) break;
 417       if (lab[i]==label) break;
 418       if (lab[i]<0) {
 419         lab[i]= label;
 420         break;
 421       }
 422     }
 423   }
 424 }
 425
 426 void AliITSclustererV2::RecPoints2Clusters
 427 (const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
 428   //------------------------------------------------------------
 429   // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS
 430   // subdetector indexed by idx
 431   //------------------------------------------------------------
 432   TClonesArray &cl=*clusters;
 433   Int_t ncl=points->GetEntriesFast();
 434   for (Int_t i=0; i<ncl; i++) {
 435     AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
 436     Float_t lp[5];
 437     lp[0]=-(-p->GetDetLocalX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1
 438     lp[1]=  -p->GetZ()+fZshift[idx];
 439     lp[2]=p->GetSigmaDetLocX2();
 440     lp[3]=p->GetSigmaZ2();
 441     lp[4]=p->GetQ()*36./23333.;  //electrons -> ADC
 442     Int_t lab[4];
 443     lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
 444     lab[3]=fNdet[idx];
 445     CheckLabels(lab);
 446     Int_t dummy[3]={0,0,0};
 447     new (cl[i]) AliITSclusterV2(lab,lp, dummy);
 448   }
 449 }
 450
 451 //***********************************
 452
 453
 454 void AliITSclustererV2::
 455 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
 456   //------------------------------------------------------------
 457   // returns an array of indices of digits belonging to the cluster
 458   // (needed when the segmentation is not regular)
 459   //------------------------------------------------------------
 460   if (n<200) idx[n++]=bins[k].GetIndex();
 461   bins[k].Use();
 462
 463   if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
 464   if (bins[k-1   ].IsNotUsed()) FindCluster(k-1   ,maxz,bins,n,idx);
 465   if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
 466   if (bins[k+1   ].IsNotUsed()) FindCluster(k+1   ,maxz,bins,n,idx);
 467   /*
 468   if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
 469   if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
 470   if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
 471   if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
 472   */
 473 }
 474
 475 void AliITSclustererV2::
 476 FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) {
 477   //------------------------------------------------------------
 478   // Actual SPD cluster finder
 479   //------------------------------------------------------------
 480   Int_t kNzBins = fNzSPD + 2;
 481   const Int_t kMAXBIN=kNzBins*(fNySPD+2);
 482
 483   Int_t ndigits=digits->GetEntriesFast();
 484   AliBin *bins=new AliBin[kMAXBIN];
 485
 486   Int_t k;
 487   AliITSdigitSPD *d=0;
 488   for (k=0; k<ndigits; k++) {
 489      d=(AliITSdigitSPD*)digits->UncheckedAt(k);
 490      Int_t i=d->GetCoord2()+1;   //y
 491      Int_t j=d->GetCoord1()+1;
 492      bins[i*kNzBins+j].SetIndex(k);
 493      bins[i*kNzBins+j].SetMask(1);
 494   }
 495
 496   Int_t n=0; TClonesArray &cl=*clusters;
 497   for (k=0; k<kMAXBIN; k++) {
 498      if (!bins[k].IsNotUsed()) continue;
 499      Int_t ni=0, idx[200];
 500      FindCluster(k,kNzBins,bins,ni,idx);
 501      if (ni==200) {
 502         Info("FindClustersSPD","Too big cluster !");
 503         continue;
 504      }
 505      Int_t milab[10];
 506      for (Int_t ilab=0;ilab<10;ilab++){
 507        milab[ilab]=-2;
 508      }
 509
 510      d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]);
 511      Int_t ymin=d->GetCoord2(),ymax=ymin;
 512      Int_t zmin=d->GetCoord1(),zmax=zmin;
 513
 514      for (Int_t l=0; l<ni; l++) {
 515         d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]);
 516
 517         if (ymin > d->GetCoord2()) ymin=d->GetCoord2();
 518         if (ymax < d->GetCoord2()) ymax=d->GetCoord2();
 519         if (zmin > d->GetCoord1()) zmin=d->GetCoord1();
 520         if (zmax < d->GetCoord1()) zmax=d->GetCoord1();
 521         // MI addition - find all labels in cluster
 522         for (Int_t dlab=0;dlab<10;dlab++){
 523           Int_t digitlab = (d->GetTracks())[dlab];
 524           if (digitlab<0) continue;
 525           AddLabel(milab,digitlab);
 526         }
 527         if (milab[9]>0) CheckLabels2(milab);
 528      }
 529      CheckLabels2(milab);
 530      //
 531      //Int_t idy = (fNlayer[fI]==0)? 2:3;
 532      //for (Int_t iz=zmin; iz<=zmax;iz+=2)
 533      //Int_t idy = (ymax-ymin)/4.; // max 2 clusters
 534      Int_t idy = 0; // max 2 clusters
 535      if (fNlayer[fI]==0 &&idy<3) idy=3;
 536      if (fNlayer[fI]==1 &&idy<4) idy=4;
 537      Int_t idz =3;
 538      for (Int_t iz=zmin; iz<=zmax;iz+=idz)
 539        for (Int_t iy=ymin; iy<=ymax;iy+=idy){
 540          //
 541          Int_t nodigits =0;
 542          Float_t y=0.,z=0.,q=0.;
 543          for (Int_t l=0; l<ni; l++) {
 544            d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]);
 545            if (zmax-zmin>=idz || ymax-ymin>=idy){
 546              if (TMath::Abs( d->GetCoord2()-iy)>0.75*idy) continue;
 547              if (TMath::Abs( d->GetCoord1()-iz)>0.75*idz) continue;
 548            }
 549            nodigits++;
 550            Float_t qq=d->GetSignal();
 551            y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq;
 552
 553          }
 554          if (nodigits==0) continue;
 555          y/=q; z/=q;
 556          y-=fHwSPD; z-=fHlSPD;
 557
 558          Float_t lp[5];
 559          lp[0]=-(-y+fYshift[fI]); if (fI<=fLastSPD1) lp[0]=-lp[0];
 560          lp[1]=  -z+fZshift[fI];
 561          // Float_t factor=TMath::Max(double(ni-3.),1.5);
 562          Float_t factory=TMath::Max(ymax-ymin,1);
 563          Float_t factorz=TMath::Max(zmax-zmin,1);
 564          factory*= factory;
 565          factorz*= factorz;
 566          //lp[2]= (fYpitchSPD*fYpitchSPD/12.)*factory;
 567          //lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.)*factorz;
 568          lp[2]= (fYpitchSPD*fYpitchSPD/12.);
 569          lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.);
 570          //lp[4]= q;
 571          lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1);
 572
 573          milab[3]=fNdet[fI];
 574          Int_t info[3] = {ymax-ymin+1,zmax-zmin+1,fNlayer[fI]};
 575          new (cl[n]) AliITSclusterV2(milab,lp,info); n++;
 576        }
 577   }
 578
 579   delete [] bins;
 580 }
 581
 582 void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input,
 583                                         TClonesArray** clusters)
 584 {
 585   //------------------------------------------------------------
 586   // Actual SPD cluster finder for raw data
 587   //------------------------------------------------------------
 588
 589   Int_t nClustersSPD = 0;
 590   Int_t kNzBins = fNzSPD + 2;
 591   Int_t kNyBins = fNySPD + 2;
 592   Int_t kMaxBin = kNzBins * kNyBins;
 593   AliBin *binsSPD = new AliBin[kMaxBin];
 594   AliBin *binsSPDInit = new AliBin[kMaxBin];
 595   AliBin *bins = NULL;
 596
 597   // read raw data input stream
 598   while (kTRUE) {
 599     Bool_t next = input->Next();
 600     if (!next || input->IsNewModule()) {
 601       Int_t iModule = input->GetPrevModuleID();
 602
 603       // when all data from a module was read, search for clusters
 604       if (bins) {
 605         clusters[iModule] = new TClonesArray("AliITSclusterV2");
 606         Int_t nClusters = 0;
 607
 608         for (Int_t iBin = 0; iBin < kMaxBin; iBin++) {
 609           if (bins[iBin].IsUsed()) continue;
 610           Int_t nBins = 0;
 611           Int_t idxBins[200];
 612           FindCluster(iBin, kNzBins, bins, nBins, idxBins);
 613           if (nBins == 200) {
 614             Error("FindClustersSPD", "SPD: Too big cluster !\n");
 615             continue;
 616           }
 617
 618           Int_t label[4];
 619           label[0] = -2;
 620           label[1] = -2;
 621           label[2] = -2;
 622 //        label[3] = iModule;
 623           label[3] = fNdet[iModule];
 624
 625           Int_t ymin = (idxBins[0] / kNzBins) - 1;
 626           Int_t ymax = ymin;
 627           Int_t zmin = (idxBins[0] % kNzBins) - 1;
 628           Int_t zmax = zmin;
 629           for (Int_t idx = 0; idx < nBins; idx++) {
 630             Int_t iy = (idxBins[idx] / kNzBins) - 1;
 631             Int_t iz = (idxBins[idx] % kNzBins) - 1;
 632             if (ymin > iy) ymin = iy;
 633             if (ymax < iy) ymax = iy;
 634             if (zmin > iz) zmin = iz;
 635             if (zmax < iz) zmax = iz;
 636           }
 637
 638           Int_t idy = 0; // max 2 clusters
 639           if ((iModule <= fLastSPD1) &&idy<3) idy=3;
 640           if ((iModule > fLastSPD1) &&idy<4) idy=4;
 641           Int_t idz =3;
 642           for (Int_t iiz=zmin; iiz<=zmax;iiz+=idz)
 643             for (Int_t iiy=ymin; iiy<=ymax;iiy+=idy){
 644               //
 645               Int_t ndigits =0;
 646               Float_t y=0.,z=0.,q=0.;
 647               for (Int_t idx = 0; idx < nBins; idx++) {
 648                 Int_t iy = (idxBins[idx] / kNzBins) - 1;
 649                 Int_t iz = (idxBins[idx] % kNzBins) - 1;
 650                 if (zmax-zmin>=idz || ymax-ymin>=idy){
 651                   if (TMath::Abs(iy-iiy)>0.75*idy) continue;
 652                   if (TMath::Abs(iz-iiz)>0.75*idz) continue;
 653                 }
 654                 ndigits++;
 655                 Float_t qBin = bins[idxBins[idx]].GetQ();
 656                 y += qBin * fYSPD[iy];
 657                 z += qBin * fZSPD[iz];
 658                 q += qBin;
 659               }
 660               if (ndigits==0) continue;
 661               y /= q;
 662               z /= q;
 663               y -= fHwSPD;
 664               z -= fHlSPD;
 665
 666               Float_t hit[5];  // y, z, sigma(y)^2, sigma(z)^2, charge
 667               hit[0] = -(-y+fYshift[iModule]);
 668               if (iModule <= fLastSPD1) hit[0] = -hit[0];
 669               hit[1] = -z+fZshift[iModule];
 670               hit[2] = fYpitchSPD*fYpitchSPD/12.;
 671               hit[3] = fZ1pitchSPD*fZ1pitchSPD/12.;
 672               //          hit[4] = q;
 673               hit[4] = (zmax-zmin+1)*100 + (ymax-ymin+1);
 674               //          CheckLabels(label);
 675               Int_t info[3]={ymax-ymin+1,zmax-zmin+1,fNlayer[iModule]};
 676               new (clusters[iModule]->AddrAt(nClusters))
 677                 AliITSclusterV2(label, hit,info);
 678               nClusters++;
 679             }
 680         }
 681
 682         nClustersSPD += nClusters;
 683         bins = NULL;
 684       }
 685
 686       if (!next) break;
 687       bins = binsSPD;
 688       memcpy(binsSPD,binsSPDInit,sizeof(AliBin)*kMaxBin);
 689     }
 690
 691     // fill the current digit into the bins array
 692     if(bins){
 693       Int_t index = (input->GetCoord2()+1) * kNzBins + (input->GetCoord1()+1);
 694       bins[index].SetIndex(index);
 695       bins[index].SetMask(1);
 696       bins[index].SetQ(1);
 697     }
 698   }
 699
 700   delete [] binsSPDInit;
 701   delete [] binsSPD;
 702
 703   Info("FindClustersSPD", "found clusters in ITS SPD: %d", nClustersSPD);
 704 }
 705
 706
 707 Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
 708   //------------------------------------------------------------
 709   //is this a local maximum ?
 710   //------------------------------------------------------------
 711   UShort_t q=bins[k].GetQ();
 712   if (q==1023) return kFALSE;
 713   if (bins[k-max].GetQ() > q) return kFALSE;
 714   if (bins[k-1  ].GetQ() > q) return kFALSE;
 715   if (bins[k+max].GetQ() > q) return kFALSE;
 716   if (bins[k+1  ].GetQ() > q) return kFALSE;
 717   if (bins[k-max-1].GetQ() > q) return kFALSE;
 718   if (bins[k+max-1].GetQ() > q) return kFALSE;
 719   if (bins[k+max+1].GetQ() > q) return kFALSE;
 720   if (bins[k-max+1].GetQ() > q) return kFALSE;
 721   return kTRUE;
 722 }
 723
 724 void AliITSclustererV2::
 725 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
 726   //------------------------------------------------------------
 727   //find local maxima
 728   //------------------------------------------------------------
 729   if (n<31)
 730   if (IsMaximum(k,max,b)) {
 731     idx[n]=k; msk[n]=(2<<n);
 732     n++;
 733   }
 734   b[k].SetMask(0);
 735   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
 736   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
 737   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
 738   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
 739 }
 740
 741 void AliITSclustererV2::
 742 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
 743   //------------------------------------------------------------
 744   //mark this peak
 745   //------------------------------------------------------------
 746   UShort_t q=bins[k].GetQ();
 747
 748   bins[k].SetMask(bins[k].GetMask()|m);
 749
 750   if (bins[k-max].GetQ() <= q)
 751      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
 752   if (bins[k-1  ].GetQ() <= q)
 753      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
 754   if (bins[k+max].GetQ() <= q)
 755      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
 756   if (bins[k+1  ].GetQ() <= q)
 757      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
 758 }
 759
 760 void AliITSclustererV2::
 761 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) {
 762   //------------------------------------------------------------
 763   //make cluster using digits of this peak
 764   //------------------------------------------------------------
 765   Float_t q=(Float_t)bins[k].GetQ();
 766   Int_t i=k/max, j=k-i*max;
 767
 768   c.SetQ(c.GetQ()+q);
 769   c.SetY(c.GetY()+i*q);
 770   c.SetZ(c.GetZ()+j*q);
 771   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
 772   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
 773
 774   bins[k].SetMask(0xFFFFFFFE);
 775
 776   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
 777   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
 778   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
 779   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
 780 }
 781
 782 void AliITSclustererV2::
 783 FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins,
 784                 const TClonesArray *digits, TClonesArray *clusters) {
 785   //------------------------------------------------------------
 786   // Actual SDD cluster finder
 787   //------------------------------------------------------------
 788   Int_t ncl=0; TClonesArray &cl=*clusters;
 789   for (Int_t s=0; s<2; s++)
 790     for (Int_t i=0; i<nMaxBin; i++) {
 791       if (bins[s][i].IsUsed()) continue;
 792       Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
 793       FindPeaks(i, nzBins, bins[s], idx, msk, npeaks);
 794
 795       if (npeaks>30) continue;
 796       if (npeaks==0) continue;
 797
 798       Int_t k,l;
 799       for (k=0; k<npeaks-1; k++){//mark adjacent peaks
 800         if (idx[k] < 0) continue; //this peak is already removed
 801         for (l=k+1; l<npeaks; l++) {
 802            if (idx[l] < 0) continue; //this peak is already removed
 803            Int_t ki=idx[k]/nzBins, kj=idx[k] - ki*nzBins;
 804            Int_t li=idx[l]/nzBins, lj=idx[l] - li*nzBins;
 805            Int_t di=TMath::Abs(ki - li);
 806            Int_t dj=TMath::Abs(kj - lj);
 807            if (di>1 || dj>1) continue;
 808            if (bins[s][idx[k]].GetQ() > bins[s][idx[l]].GetQ()) {
 809               msk[l]=msk[k];
 810               idx[l]*=-1;
 811            } else {
 812               msk[k]=msk[l];
 813               idx[k]*=-1;
 814               break;
 815            }
 816         }
 817       }
 818
 819       for (k=0; k<npeaks; k++) {
 820         MarkPeak(TMath::Abs(idx[k]), nzBins, bins[s], msk[k]);
 821       }
 822
 823       for (k=0; k<npeaks; k++) {
 824          if (idx[k] < 0) continue; //removed peak
 825          AliITSclusterV2 c;
 826          MakeCluster(idx[k], nzBins, bins[s], msk[k], c);
 827          //mi change
 828          Int_t milab[10];
 829          for (Int_t ilab=0;ilab<10;ilab++){
 830            milab[ilab]=-2;
 831          }
 832          Int_t maxi=0,mini=0,maxj=0,minj=0;
 833          //AliBin *bmax=&bins[s][idx[k]];
 834          //Float_t max = TMath::Max(TMath::Abs(bmax->GetQ())/5.,3.);
 835          Float_t max=3;
 836          for (Int_t di=-2; di<=2;di++)
 837            for (Int_t dj=-3;dj<=3;dj++){
 838              Int_t index = idx[k]+di+dj*nzBins;
 839              if (index<0) continue;
 840              if (index>=nMaxBin) continue;
 841              AliBin *b=&bins[s][index];
 842              if (TMath::Abs(b->GetQ())>max){
 843                if (di>maxi) maxi=di;
 844                if (di<mini) mini=di;
 845                if (dj>maxj) maxj=dj;
 846                if (dj<minj) minj=dj;
 847                //
 848                if(digits) {
 849                  if (TMath::Abs(di)<2&&TMath::Abs(dj)<2){
 850                    AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 851                    for (Int_t itrack=0;itrack<10;itrack++){
 852                      Int_t track = (d->GetTracks())[itrack];
 853                      if (track>=0) {
 854                        AddLabel(milab, track);
 855                      }
 856                    }
 857                  }
 858                }
 859              }
 860            }
 861
 862          /*
 863             Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
 864             Float_t w=par->GetPadPitchWidth(sec);
 865             c.SetSigmaY2(s2);
 866             if (s2 != 0.) {
 867             c.SetSigmaY2(c.GetSigmaY2()*0.108);
 868             if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
 869             }
 870             s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
 871             w=par->GetZWidth();
 872             c.SetSigmaZ2(s2);
 873
 874             if (s2 != 0.) {
 875             c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
 876             if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
 877             }
 878          */
 879
 880          Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ();
 881          y/=q; z/=q;
 882          //
 883          //Float_t s2 = c.GetSigmaY2()/c.GetQ() - y*y;
 884          // c.SetSigmaY2(s2);
 885          //s2 = c.GetSigmaZ2()/c.GetQ() - z*z;
 886          //c.SetSigmaZ2(s2);
 887          //
 888          y=(y-0.5)*fYpitchSDD;
 889          y-=fHwSDD;
 890          y-=fYoffSDD;  //delay ?
 891          if (s) y=-y;
 892
 893          z=(z-0.5)*fZpitchSDD;
 894          z-=fHlSDD;
 895
 896          y=-(-y+fYshift[fI]);
 897          z=  -z+fZshift[fI];
 898
 899          q/=12.7;  //this WAS consistent with SSD. To be reassessed
 900                    // 23-MAR-2007
 901          Float_t hit[5] = {y, z, 0.0030*0.0030, 0.0020*0.0020, q};
 902          Int_t  info[3] = {maxj-minj+1, maxi-mini+1, fNlayer[fI]};
 903
 904          if (c.GetQ() < 20.) continue; //noise cluster
 905
 906          if (digits) {
 907            //      AliBin *b=&bins[s][idx[k]];
 908            //      AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 909            {
 910              //Int_t lab[3];
 911              //lab[0]=(d->GetTracks())[0];
 912              //lab[1]=(d->GetTracks())[1];
 913              //lab[2]=(d->GetTracks())[2];
 914              //CheckLabels(lab);
 915              CheckLabels2(milab);
 916            }
 917          }
 918          milab[3]=fNdet[fI];
 919
 920          AliITSclusterV2 *cc = new (cl[ncl]) AliITSclusterV2(milab,hit,info);
 921          cc->SetType(npeaks);
 922          ncl++;
 923       }
 924     }
 925 }
 926
 927 void AliITSclustererV2::
 928 FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) {
 929   //------------------------------------------------------------
 930   // Actual SDD cluster finder
 931   //------------------------------------------------------------
 932   Int_t kNzBins = fNzSDD + 2;
 933   const Int_t kMAXBIN=kNzBins*(fNySDD+2);
 934
 935   AliBin *bins[2];
 936   bins[0]=new AliBin[kMAXBIN];
 937   bins[1]=new AliBin[kMAXBIN];
 938
 939   AliITSdigitSDD *d=0;
 940   Int_t i, ndigits=digits->GetEntriesFast();
 941   for (i=0; i<ndigits; i++) {
 942      d=(AliITSdigitSDD*)digits->UncheckedAt(i);
 943      Int_t y=d->GetCoord2()+1;   //y
 944      Int_t z=d->GetCoord1()+1;   //z
 945      Int_t q=d->GetSignal();
 946      if (q<3) continue;
 947
 948      if (z <= fNzSDD) {
 949        bins[0][y*kNzBins+z].SetQ(q);
 950        bins[0][y*kNzBins+z].SetMask(1);
 951        bins[0][y*kNzBins+z].SetIndex(i);
 952      } else {
 953        z-=fNzSDD;
 954        bins[1][y*kNzBins+z].SetQ(q);
 955        bins[1][y*kNzBins+z].SetMask(1);
 956        bins[1][y*kNzBins+z].SetIndex(i);
 957      }
 958   }
 959
 960   FindClustersSDD(bins, kMAXBIN, kNzBins, digits, clusters);
 961
 962   delete[] bins[0];
 963   delete[] bins[1];
 964
 965 }
 966
 967 void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input,
 968                                         TClonesArray** clusters)
 969 {
 970   //------------------------------------------------------------
 971   // Actual SDD cluster finder for raw data
 972   //------------------------------------------------------------
 973   Int_t nClustersSDD = 0;
 974   Int_t kNzBins = fNzSDD + 2;
 975   Int_t kMaxBin = kNzBins * (fNySDD+2);
 976   AliBin *binsSDDInit = new AliBin[kMaxBin];
 977   AliBin *binsSDD1 = new AliBin[kMaxBin];
 978   AliBin *binsSDD2 = new AliBin[kMaxBin];
 979   AliBin *bins[2] = {NULL, NULL};
 980
 981   // read raw data input stream
 982   while (kTRUE) {
 983     Bool_t next = input->Next();
 984     if (!next || input->IsNewModule()) {
 985       Int_t iModule = input->GetPrevModuleID();
 986
 987       // when all data from a module was read, search for clusters
 988       if (bins[0]) {
 989         clusters[iModule] = new TClonesArray("AliITSclusterV2");
 990         fI = iModule;
 991         FindClustersSDD(bins, kMaxBin, kNzBins, NULL, clusters[iModule]);
 992         Int_t nClusters = clusters[iModule]->GetEntriesFast();
 993         nClustersSDD += nClusters;
 994         bins[0] = bins[1] = NULL;
 995       }
 996
 997       if (!next) break;
 998       bins[0]=binsSDD1;
 999       bins[1]=binsSDD2;
1000       memcpy(binsSDD1,binsSDDInit,sizeof(AliBin)*kMaxBin);
1001       memcpy(binsSDD2,binsSDDInit,sizeof(AliBin)*kMaxBin);
1002     }
1003
1004     // fill the current digit into the bins array
1005     if(input->GetSignal()>=3) {
1006       Int_t iz = input->GetCoord1()+1;
1007       Int_t side = ((iz <= fNzSDD) ? 0 : 1);
1008       iz -= side*fNzSDD;
1009       Int_t index = (input->GetCoord2()+1) * kNzBins + iz;
1010       bins[side][index].SetQ(input->GetSignal());
1011       bins[side][index].SetMask(1);
1012       bins[side][index].SetIndex(index);
1013     }
1014   }
1015
1016   delete[] binsSDD1;
1017   delete[] binsSDD2;
1018   delete[] binsSDDInit;
1019
1020   Info("FindClustersSDD", "found clusters in ITS SDD: %d", nClustersSDD);
1021 }
1022
1023
1024
1025 void AliITSclustererV2::
1026 FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
1027                 Ali1Dcluster* pos, Int_t np,
1028                 TClonesArray *clusters) {
1029   //------------------------------------------------------------
1030   // Actual SSD cluster finder
1031   //------------------------------------------------------------
1032   TClonesArray &cl=*clusters;
1033   //
1034   Float_t tanp=fTanP, tann=fTanN;
1035   if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;}
1036   Int_t idet=fNdet[fI];
1037   Int_t ncl=0;
1038   //
1039   Int_t negativepair[30000];
1040   Int_t cnegative[3000];
1041   Int_t cused1[3000];
1042   Int_t positivepair[30000];
1043   Int_t cpositive[3000];
1044   Int_t cused2[3000];
1045   for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
1046   for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
1047   for (Int_t i=0;i<30000;i++){negativepair[i]=positivepair[i]=0;}
1048   static Short_t pairs[1000][1000];
1049   memset(pairs,0,sizeof(Short_t)*1000000);
1050 //   Short_t ** pairs = new Short_t*[1000];
1051 //   for (Int_t i=0; i<1000; i++) {
1052 //     pairs[i] = new Short_t[1000];
1053 //     memset(pairs[i],0,sizeof(Short_t)*1000);
1054 //   }
1055   //
1056   // find available pairs
1057   //
1058   for (Int_t i=0; i<np; i++) {
1059     Float_t yp=pos[i].GetY()*fYpitchSSD;
1060     if (pos[i].GetQ()<3) continue;
1061     for (Int_t j=0; j<nn; j++) {
1062       if (neg[j].GetQ()<3) continue;
1063       Float_t yn=neg[j].GetY()*fYpitchSSD;
1064       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1065       Float_t yt=yn + tann*zt;
1066       zt-=fHlSSD; yt-=fHwSSD;
1067       if (TMath::Abs(yt)<fHwSSD+0.01)
1068       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1069         negativepair[i*10+cnegative[i]] =j;  //index
1070         positivepair[j*10+cpositive[j]] =i;
1071         cnegative[i]++;  //counters
1072         cpositive[j]++;
1073         pairs[i][j]=100;
1074       }
1075     }
1076   }
1077   //
1078   for (Int_t i=0; i<np; i++) {
1079     Float_t yp=pos[i].GetY()*fYpitchSSD;
1080     if (pos[i].GetQ()<3) continue;
1081     for (Int_t j=0; j<nn; j++) {
1082       if (neg[j].GetQ()<3) continue;
1083       if (cpositive[j]&&cnegative[i]) continue;
1084       Float_t yn=neg[j].GetY()*fYpitchSSD;
1085       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1086       Float_t yt=yn + tann*zt;
1087       zt-=fHlSSD; yt-=fHwSSD;
1088       if (TMath::Abs(yt)<fHwSSD+0.1)
1089       if (TMath::Abs(zt)<fHlSSD+0.15) {
1090         if (cnegative[i]==0) pos[i].SetNd(100);  // not available pair
1091         if (cpositive[j]==0) neg[j].SetNd(100);  // not available pair
1092         negativepair[i*10+cnegative[i]] =j;  //index
1093         positivepair[j*10+cpositive[j]] =i;
1094         cnegative[i]++;  //counters
1095         cpositive[j]++;
1096         pairs[i][j]=100;
1097       }
1098     }
1099   }
1100   //
1101   Float_t lp[5];
1102   Int_t milab[10];
1103   Double_t ratio;
1104
1105   //
1106   // sign gold tracks
1107   //
1108   for (Int_t ip=0;ip<np;ip++){
1109     Float_t ybest=1000,zbest=1000,qbest=0;
1110     //
1111     // select gold clusters
1112     if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
1113       Float_t yp=pos[ip].GetY()*fYpitchSSD;
1114       Int_t j = negativepair[10*ip];
1115       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1116       //
1117       Float_t yn=neg[j].GetY()*fYpitchSSD;
1118       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1119       Float_t yt=yn + tann*zt;
1120       zt-=fHlSSD; yt-=fHwSSD;
1121       ybest=yt; zbest=zt;
1122       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1123       lp[0]=-(-ybest+fYshift[fI]);
1124       lp[1]=  -zbest+fZshift[fI];
1125       lp[2]=0.0025*0.0025;  //SigmaY2
1126       lp[3]=0.110*0.110;  //SigmaZ2
1127
1128       lp[4]=qbest;        //Q
1129       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1130       for (Int_t ilab=0;ilab<3;ilab++){
1131         milab[ilab] = pos[ip].GetLabel(ilab);
1132         milab[ilab+3] = neg[j].GetLabel(ilab);
1133       }
1134       //
1135       CheckLabels2(milab);
1136       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1137       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]};
1138       AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1139       ncl++;
1140       cl2->SetChargeRatio(ratio);
1141       cl2->SetType(1);
1142       pairs[ip][j]=1;
1143       if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1144         cl2->SetType(2);
1145         pairs[ip][j]=2;
1146       }
1147       cused1[ip]++;
1148       cused2[j]++;
1149     }
1150   }
1151
1152   for (Int_t ip=0;ip<np;ip++){
1153     Float_t ybest=1000,zbest=1000,qbest=0;
1154     //
1155     //
1156     // select "silber" cluster
1157     if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
1158       Int_t in  = negativepair[10*ip];
1159       Int_t ip2 = positivepair[10*in];
1160       if (ip2==ip) ip2 =  positivepair[10*in+1];
1161       Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
1162       if (TMath::Abs(pcharge-neg[in].GetQ())<10){
1163         //
1164         // add first pair
1165         if (pairs[ip][in]==100){  //
1166           Float_t yp=pos[ip].GetY()*fYpitchSSD;
1167           Float_t yn=neg[in].GetY()*fYpitchSSD;
1168           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1169           Float_t yt=yn + tann*zt;
1170           zt-=fHlSSD; yt-=fHwSSD;
1171           ybest =yt;  zbest=zt;
1172           qbest =pos[ip].GetQ();
1173           lp[0]=-(-ybest+fYshift[fI]);
1174           lp[1]=  -zbest+fZshift[fI];
1175           lp[2]=0.0025*0.0025;  //SigmaY2
1176           lp[3]=0.110*0.110;  //SigmaZ2
1177
1178           lp[4]=qbest;        //Q
1179           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1180           for (Int_t ilab=0;ilab<3;ilab++){
1181             milab[ilab] = pos[ip].GetLabel(ilab);
1182             milab[ilab+3] = neg[in].GetLabel(ilab);
1183           }
1184           //
1185           CheckLabels2(milab);
1186           ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
1187           milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
1188           Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fI]};
1189           AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1190           ncl++;
1191           cl2->SetChargeRatio(ratio);
1192           cl2->SetType(5);
1193           pairs[ip][in] = 5;
1194           if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1195             cl2->SetType(6);
1196             pairs[ip][in] = 6;
1197           }
1198         }
1199         //
1200         // add second pair
1201
1202         //      if (!(cused1[ip2] || cused2[in])){  //
1203         if (pairs[ip2][in]==100){
1204           Float_t yp=pos[ip2].GetY()*fYpitchSSD;
1205           Float_t yn=neg[in].GetY()*fYpitchSSD;
1206           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1207           Float_t yt=yn + tann*zt;
1208           zt-=fHlSSD; yt-=fHwSSD;
1209           ybest =yt;  zbest=zt;
1210           qbest =pos[ip2].GetQ();
1211           lp[0]=-(-ybest+fYshift[fI]);
1212           lp[1]=  -zbest+fZshift[fI];
1213           lp[2]=0.0025*0.0025;  //SigmaY2
1214           lp[3]=0.110*0.110;  //SigmaZ2
1215
1216           lp[4]=qbest;        //Q
1217           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1218           for (Int_t ilab=0;ilab<3;ilab++){
1219             milab[ilab] = pos[ip2].GetLabel(ilab);
1220             milab[ilab+3] = neg[in].GetLabel(ilab);
1221           }
1222           //
1223           CheckLabels2(milab);
1224           ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
1225           milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
1226           Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fI]};
1227           AliITSclusterV2 *cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1228           ncl++;
1229           cl2->SetChargeRatio(ratio);
1230           cl2->SetType(5);
1231           pairs[ip2][in] =5;
1232           if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1233             cl2->SetType(6);
1234             pairs[ip2][in] =6;
1235           }
1236         }
1237         cused1[ip]++;
1238         cused1[ip2]++;
1239         cused2[in]++;
1240       }
1241     }
1242   }
1243
1244   //
1245   for (Int_t jn=0;jn<nn;jn++){
1246     if (cused2[jn]) continue;
1247     Float_t ybest=1000,zbest=1000,qbest=0;
1248     // select "silber" cluster
1249     if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
1250       Int_t ip  = positivepair[10*jn];
1251       Int_t jn2 = negativepair[10*ip];
1252       if (jn2==jn) jn2 =  negativepair[10*ip+1];
1253       Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
1254       //
1255       if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
1256         //
1257         // add first pair
1258         //      if (!(cused1[ip]||cused2[jn])){
1259         if (pairs[ip][jn]==100){
1260           Float_t yn=neg[jn].GetY()*fYpitchSSD;
1261           Float_t yp=pos[ip].GetY()*fYpitchSSD;
1262           Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1263           Float_t yt=yn + tann*zt;
1264           zt-=fHlSSD; yt-=fHwSSD;
1265           ybest =yt;  zbest=zt;
1266           qbest =neg[jn].GetQ();
1267           lp[0]=-(-ybest+fYshift[fI]);
1268           lp[1]=  -zbest+fZshift[fI];
1269           lp[2]=0.0025*0.0025;  //SigmaY2
1270           lp[3]=0.110*0.110;  //SigmaZ2
1271
1272           lp[4]=qbest;        //Q
1273           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1274           for (Int_t ilab=0;ilab<3;ilab++){
1275             milab[ilab] = pos[ip].GetLabel(ilab);
1276             milab[ilab+3] = neg[jn].GetLabel(ilab);
1277           }
1278           //
1279           CheckLabels2(milab);
1280           ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1281           milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1282           Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fI]};
1283           AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1284           ncl++;
1285           cl2->SetChargeRatio(ratio);
1286           cl2->SetType(7);
1287           pairs[ip][jn] =7;
1288           if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1289             cl2->SetType(8);
1290             pairs[ip][jn]=8;
1291           }
1292         }
1293         //
1294         // add second pair
1295         //      if (!(cused1[ip]||cused2[jn2])){
1296         if (pairs[ip][jn2]==100){
1297           Float_t yn=neg[jn2].GetY()*fYpitchSSD;
1298           Double_t yp=pos[ip].GetY()*fYpitchSSD;
1299           Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1300           Double_t yt=yn + tann*zt;
1301           zt-=fHlSSD; yt-=fHwSSD;
1302           ybest =yt;  zbest=zt;
1303           qbest =neg[jn2].GetQ();
1304           lp[0]=-(-ybest+fYshift[fI]);
1305           lp[1]=  -zbest+fZshift[fI];
1306           lp[2]=0.0025*0.0025;  //SigmaY2
1307           lp[3]=0.110*0.110;  //SigmaZ2
1308
1309           lp[4]=qbest;        //Q
1310           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1311           for (Int_t ilab=0;ilab<3;ilab++){
1312             milab[ilab] = pos[ip].GetLabel(ilab);
1313             milab[ilab+3] = neg[jn2].GetLabel(ilab);
1314           }
1315           //
1316           CheckLabels2(milab);
1317           ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1318           milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1319           Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fI]};
1320           AliITSclusterV2* cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1321           ncl++;
1322           cl2->SetChargeRatio(ratio);
1323           pairs[ip][jn2]=7;
1324           cl2->SetType(7);
1325           if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1326             cl2->SetType(8);
1327             pairs[ip][jn2]=8;
1328           }
1329         }
1330         cused1[ip]++;
1331         cused2[jn]++;
1332         cused2[jn2]++;
1333       }
1334     }
1335   }
1336
1337   for (Int_t ip=0;ip<np;ip++){
1338     Float_t ybest=1000,zbest=1000,qbest=0;
1339     //
1340     // 2x2 clusters
1341     //
1342     if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1343       Float_t minchargediff =4.;
1344       Int_t j=-1;
1345       for (Int_t di=0;di<cnegative[ip];di++){
1346         Int_t   jc = negativepair[ip*10+di];
1347         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1348         if (TMath::Abs(chargedif)<minchargediff){
1349           j =jc;
1350           minchargediff = TMath::Abs(chargedif);
1351         }
1352       }
1353       if (j<0) continue;  // not proper cluster
1354       Int_t count =0;
1355       for (Int_t di=0;di<cnegative[ip];di++){
1356         Int_t   jc = negativepair[ip*10+di];
1357         Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1358         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1359       }
1360       if (count>1) continue;  // more than one "proper" cluster for positive
1361       //
1362       count =0;
1363       for (Int_t dj=0;dj<cpositive[j];dj++){
1364         Int_t   ic  = positivepair[j*10+dj];
1365         Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1366         if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1367       }
1368       if (count>1) continue;  // more than one "proper" cluster for negative
1369
1370       Int_t jp = 0;
1371
1372       count =0;
1373       for (Int_t dj=0;dj<cnegative[jp];dj++){
1374         Int_t   ic = positivepair[jp*10+dj];
1375         Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1376         if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1377       }
1378       if (count>1) continue;
1379       if (pairs[ip][j]<100) continue;
1380       //
1381       //almost gold clusters
1382       Float_t yp=pos[ip].GetY()*fYpitchSSD;
1383       Float_t yn=neg[j].GetY()*fYpitchSSD;
1384       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1385       Float_t yt=yn + tann*zt;
1386       zt-=fHlSSD; yt-=fHwSSD;
1387       ybest=yt; zbest=zt;
1388       qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1389       lp[0]=-(-ybest+fYshift[fI]);
1390       lp[1]=  -zbest+fZshift[fI];
1391       lp[2]=0.0025*0.0025;  //SigmaY2
1392       lp[3]=0.110*0.110;  //SigmaZ2
1393       lp[4]=qbest;        //Q
1394       for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1395       for (Int_t ilab=0;ilab<3;ilab++){
1396         milab[ilab] = pos[ip].GetLabel(ilab);
1397         milab[ilab+3] = neg[j].GetLabel(ilab);
1398       }
1399       //
1400       CheckLabels2(milab);
1401       ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1402       milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1403       Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]};
1404       AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1405       ncl++;
1406       cl2->SetChargeRatio(ratio);
1407       cl2->SetType(10);
1408       pairs[ip][j]=10;
1409       if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1410         cl2->SetType(11);
1411         pairs[ip][j]=11;
1412       }
1413       cused1[ip]++;
1414       cused2[j]++;
1415     }
1416
1417   }
1418
1419   //
1420   for (Int_t i=0; i<np; i++) {
1421     Float_t ybest=1000,zbest=1000,qbest=0;
1422     Float_t yp=pos[i].GetY()*fYpitchSSD;
1423     if (pos[i].GetQ()<3) continue;
1424     for (Int_t j=0; j<nn; j++) {
1425     //    for (Int_t di = 0;di<cpositive[i];di++){
1426     //  Int_t j = negativepair[10*i+di];
1427       if (neg[j].GetQ()<3) continue;
1428       if (cused2[j]||cused1[i]) continue;
1429       if (pairs[i][j]>0 &&pairs[i][j]<100) continue;
1430       ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1431       Float_t yn=neg[j].GetY()*fYpitchSSD;
1432       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1433       Float_t yt=yn + tann*zt;
1434       zt-=fHlSSD; yt-=fHwSSD;
1435       if (TMath::Abs(yt)<fHwSSD+0.01)
1436       if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1437         ybest=yt; zbest=zt;
1438         qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1439         lp[0]=-(-ybest+fYshift[fI]);
1440         lp[1]=  -zbest+fZshift[fI];
1441         lp[2]=0.0025*0.0025;  //SigmaY2
1442         lp[3]=0.110*0.110;  //SigmaZ2
1443
1444         lp[4]=qbest;        //Q
1445         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1446         for (Int_t ilab=0;ilab<3;ilab++){
1447           milab[ilab] = pos[i].GetLabel(ilab);
1448           milab[ilab+3] = neg[j].GetLabel(ilab);
1449         }
1450         //
1451         CheckLabels2(milab);
1452         milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1453         Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fI]};
1454         AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
1455         ncl++;
1456         cl2->SetChargeRatio(ratio);
1457         cl2->SetType(100+cpositive[j]+cnegative[i]);
1458         //cl2->SetType(0);
1459         /*
1460           if (pairs[i][j]<100){
1461           printf("problem:- %d\n", pairs[i][j]);
1462           }
1463           if (cnegative[i]<2&&cpositive[j]<2){
1464           printf("problem:- %d\n", pairs[i][j]);
1465           }
1466         */
1467       }
1468     }
1469   }
1470
1471 //   for (Int_t i=0; i<1000; i++) delete [] pairs[i];
1472 //   delete [] pairs;
1473
1474 }
1475
1476
1477 void AliITSclustererV2::
1478 FindClustersSSD(const TClonesArray *alldigits, TClonesArray *clusters) {
1479   //------------------------------------------------------------
1480   // Actual SSD cluster finder
1481   //------------------------------------------------------------
1482   Int_t smaxall=alldigits->GetEntriesFast();
1483   if (smaxall==0) return;
1484   TObjArray *digits = new TObjArray;
1485   for (Int_t i=0;i<smaxall; i++){
1486     AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
1487     if (d->GetSignal()<3) continue;
1488     digits->AddLast(d);
1489   }
1490   Int_t smax = digits->GetEntriesFast();
1491   if (smax==0) return;
1492
1493   const Int_t MAX=1000;
1494   Int_t np=0, nn=0;
1495   Ali1Dcluster pos[MAX], neg[MAX];
1496   Float_t y=0., q=0., qmax=0.;
1497   Int_t lab[4]={-2,-2,-2,-2};
1498
1499   AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0);
1500   q += d->GetSignal();
1501   y += d->GetCoord2()*d->GetSignal();
1502   qmax=d->GetSignal();
1503   lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
1504   Int_t curr=d->GetCoord2();
1505   Int_t flag=d->GetCoord1();
1506   Int_t *n=&nn;
1507   Ali1Dcluster *c=neg;
1508   Int_t nd=1;
1509   Int_t milab[10];
1510   for (Int_t ilab=0;ilab<10;ilab++){
1511     milab[ilab]=-2;
1512   }
1513   milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
1514
1515   for (Int_t s=1; s<smax; s++) {
1516       d=(AliITSdigitSSD*)digits->UncheckedAt(s);
1517       Int_t strip=d->GetCoord2();
1518       if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
1519          c[*n].SetY(y/q);
1520          c[*n].SetQ(q);
1521          c[*n].SetNd(nd);
1522          CheckLabels2(milab);
1523          c[*n].SetLabels(milab);
1524          //Split suspiciously big cluster
1525          /*
1526          if (nd>10&&nd<16){
1527            c[*n].SetY(y/q-0.3*nd);
1528            c[*n].SetQ(0.5*q);
1529            (*n)++;
1530            if (*n==MAX) {
1531              Error("FindClustersSSD","Too many 1D clusters !");
1532               return;
1533            }
1534            c[*n].SetY(y/q-0.0*nd);
1535            c[*n].SetQ(0.5*q);
1536            c[*n].SetNd(nd);
1537            (*n)++;
1538            if (*n==MAX) {
1539              Error("FindClustersSSD","Too many 1D clusters !");
1540               return;
1541            }
1542            //
1543            c[*n].SetY(y/q+0.3*nd);
1544            c[*n].SetQ(0.5*q);
1545            c[*n].SetNd(nd);
1546            c[*n].SetLabels(milab);
1547          }
1548          else{
1549          */
1550          if (nd>4&&nd<25) {
1551            c[*n].SetY(y/q-0.25*nd);
1552            c[*n].SetQ(0.5*q);
1553            (*n)++;
1554            if (*n==MAX) {
1555              Error("FindClustersSSD","Too many 1D clusters !");
1556              return;
1557            }
1558            c[*n].SetY(y/q+0.25*nd);
1559            c[*n].SetQ(0.5*q);
1560            c[*n].SetNd(nd);
1561            c[*n].SetLabels(milab);
1562          }
1563          (*n)++;
1564          if (*n==MAX) {
1565           Error("FindClustersSSD","Too many 1D clusters !");
1566           return;
1567          }
1568          y=q=qmax=0.;
1569          nd=0;
1570          lab[0]=lab[1]=lab[2]=-2;
1571          //
1572          for (Int_t ilab=0;ilab<10;ilab++){
1573            milab[ilab]=-2;
1574          }
1575          //
1576          if (flag!=d->GetCoord1()) { n=&np; c=pos; }
1577       }
1578       flag=d->GetCoord1();
1579       q += d->GetSignal();
1580       y += d->GetCoord2()*d->GetSignal();
1581       nd++;
1582       if (d->GetSignal()>qmax) {
1583          qmax=d->GetSignal();
1584          lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
1585       }
1586       for (Int_t ilab=0;ilab<10;ilab++) {
1587         if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
1588       }
1589       curr=strip;
1590   }
1591   c[*n].SetY(y/q);
1592   c[*n].SetQ(q);
1593   c[*n].SetNd(nd);
1594   c[*n].SetLabels(lab);
1595   //Split suspiciously big cluster
1596   if (nd>4 && nd<25) {
1597      c[*n].SetY(y/q-0.25*nd);
1598      c[*n].SetQ(0.5*q);
1599      (*n)++;
1600      if (*n==MAX) {
1601         Error("FindClustersSSD","Too many 1D clusters !");
1602         return;
1603      }
1604      c[*n].SetY(y/q+0.25*nd);
1605      c[*n].SetQ(0.5*q);
1606      c[*n].SetNd(nd);
1607      c[*n].SetLabels(lab);
1608   }
1609   (*n)++;
1610   if (*n==MAX) {
1611      Error("FindClustersSSD","Too many 1D clusters !");
1612      return;
1613   }
1614
1615   FindClustersSSD(neg, nn, pos, np, clusters);
1616 }
1617
1618 void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input,
1619                                         TClonesArray** clusters)
1620 {
1621   //------------------------------------------------------------
1622   // Actual SSD cluster finder for raw data
1623   //------------------------------------------------------------
1624   Int_t nClustersSSD = 0;
1625   const Int_t MAX = 1000;
1626   Ali1Dcluster clusters1D[2][MAX];
1627   Int_t nClusters[2] = {0, 0};
1628   Int_t lab[3]={-2,-2,-2};
1629   Float_t q = 0.;
1630   Float_t y = 0.;
1631   Int_t nDigits = 0;
1632   Int_t prevStrip = -1;
1633   Int_t prevFlag = -1;
1634   Int_t prevModule = -1;
1635
1636   // read raw data input stream
1637   while (kTRUE) {
1638     Bool_t next = input->Next();
1639
1640     if(input->GetSignal()<3 && next) continue;
1641     // check if a new cluster starts
1642     Int_t strip = input->GetCoord2();
1643     Int_t flag = input->GetCoord1();
1644     if ((!next || (input->GetModuleID() != prevModule)||
1645          (strip-prevStrip > 1) || (flag != prevFlag)) &&
1646         (nDigits > 0)) {
1647       if (nClusters[prevFlag] == MAX) {
1648         Error("FindClustersSSD", "Too many 1D clusters !");
1649         return;
1650       }
1651       Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++];
1652       cluster.SetY(y/q);
1653       cluster.SetQ(q);
1654       cluster.SetNd(nDigits);
1655       cluster.SetLabels(lab);
1656
1657       //Split suspiciously big cluster
1658       if (nDigits > 4&&nDigits < 25) {
1659         cluster.SetY(y/q - 0.25*nDigits);
1660         cluster.SetQ(0.5*q);
1661         if (nClusters[prevFlag] == MAX) {
1662           Error("FindClustersSSD", "Too many 1D clusters !");
1663           return;
1664         }
1665         Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++];
1666         cluster2.SetY(y/q + 0.25*nDigits);
1667         cluster2.SetQ(0.5*q);
1668         cluster2.SetNd(nDigits);
1669         cluster2.SetLabels(lab);
1670       }
1671       y = q = 0.;
1672       nDigits = 0;
1673     }
1674
1675     if (!next || (input->GetModuleID() != prevModule)) {
1676       Int_t iModule = prevModule;
1677
1678       // when all data from a module was read, search for clusters
1679       if (prevFlag >= 0) {
1680         clusters[iModule] = new TClonesArray("AliITSclusterV2");
1681         fI = iModule;
1682         FindClustersSSD(&clusters1D[0][0], nClusters[0],
1683                         &clusters1D[1][0], nClusters[1], clusters[iModule]);
1684         Int_t noClusters = clusters[iModule]->GetEntriesFast();
1685         nClustersSSD += noClusters;
1686       }
1687
1688       if (!next) break;
1689       nClusters[0] = nClusters[1] = 0;
1690       y = q = 0.;
1691       nDigits = 0;
1692     }
1693
1694     // add digit to current cluster
1695     q += input->GetSignal();
1696     y += strip * input->GetSignal();
1697     nDigits++;
1698     prevStrip = strip;
1699     prevFlag = flag;
1700     prevModule = input->GetModuleID();
1701
1702   }
1703
1704   Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
1705 }
1706