ITS/AliITSclustererV2.cxx

   1 //-------------------------------------------------------------------------
   2 //            Implementation of the ITS clusterer V2 class
   3 //
   4 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
   5 //-------------------------------------------------------------------------
   6 //uncomment the line below for running with V1 cluster finder classes
   7 //#define V1
   8
   9 #include <stdlib.h>
  10
  11 #include "AliRun.h"
  12
  13 #include "AliITSclustererV2.h"
  14 #include "AliITSclusterV2.h"
  15
  16 #include <Riostream.h>
  17 #include <TFile.h>
  18 #include <TTree.h>
  19 #include <TClonesArray.h>
  20 #include "AliITSgeom.h"
  21 #include "AliITSdigit.h"
  22
  23 ClassImp(AliITSclustererV2)
  24
  25 extern AliRun *gAlice;
  26
  27 AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) {
  28   //------------------------------------------------------------
  29   // Standard constructor
  30   //------------------------------------------------------------
  31   AliITSgeom *g=(AliITSgeom*)geom;
  32
  33   fEvent=0;
  34   fI=0;
  35
  36   Int_t mmax=geom->GetIndexMax();
  37   if (mmax>2200) {cerr<<"Too many ITS subdetectors !\n"; exit(1);}
  38   Int_t m;
  39   for (m=0; m<mmax; m++) {
  40      Int_t lay,lad,det; g->GetModuleId(m,lay,lad,det);
  41      Float_t x,y,z;     g->GetTrans(lay,lad,det,x,y,z);
  42      Double_t rot[9];   g->GetRotMatrix(lay,lad,det,rot);
  43      fYshift[m] = x*rot[0] + y*rot[1];
  44      fZshift[m] = (Double_t)z;
  45      fNdet[m] = (lad-1)*g->GetNdetectors(lay) + (det-1);
  46   }
  47
  48   //SPD geometry
  49   fLastSPD1=g->GetModuleIndex(2,1,1)-1;
  50   fNySPD=256; fNzSPD=160;
  51   fYpitchSPD=0.0050;
  52   fZ1pitchSPD=0.0425; fZ2pitchSPD=0.0625;
  53   fHwSPD=0.64; fHlSPD=3.48;
  54   fYSPD[0]=0.5*fYpitchSPD;
  55   for (m=1; m<fNySPD; m++) fYSPD[m]=fYSPD[m-1]+fYpitchSPD;
  56   fZSPD[0]=fZ1pitchSPD;
  57   for (m=1; m<fNzSPD; m++) {
  58     Double_t dz=fZ1pitchSPD;
  59     if (m==31 || m==32 || m==63  || m==64  || m==95 || m==96 ||
  60         m==127 || m==128) dz=fZ2pitchSPD;
  61     fZSPD[m]=fZSPD[m-1]+dz;
  62   }
  63   for (m=0; m<fNzSPD; m++) {
  64     Double_t dz=0.5*fZ1pitchSPD;
  65     if (m==31 || m==32 || m==63  || m==64  || m==95 || m==96 ||
  66         m==127 || m==128) dz=0.5*fZ2pitchSPD;
  67     fZSPD[m]-=dz;
  68   }
  69
  70   //SDD geometry
  71   fNySDD=256; fNzSDD=256;
  72   fYpitchSDD=0.01825;
  73   fZpitchSDD=0.02940;
  74   fHwSDD=3.5085; fHlSDD=3.7632;
  75   fYoffSDD=0.0425;
  76
  77   //SSD geometry
  78   fLastSSD1=g->GetModuleIndex(6,1,1)-1;
  79   fYpitchSSD=0.0095;
  80   fHwSSD=3.65;
  81   fHlSSD=2.00;
  82   fTanP=0.0275;
  83   fTanN=0.0075;
  84 }
  85
  86 void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) {
  87   //------------------------------------------------------------
  88   // This function creates ITS clusters
  89   //------------------------------------------------------------
  90   Int_t ncl=0;
  91   TDirectory *savedir=gDirectory;
  92
  93   if (!out->IsOpen()) {
  94     cerr<<"AliITSclustererV2::Digits2Clusters(): output file not open !\n";
  95     return;
  96   }
  97
  98   Char_t name[100];
  99   sprintf(name,"TreeD%d",fEvent);
 100
 101   //TTree *dTree=(TTree*)((TFile*)in)->Get(name);
 102   TTree *dTree=gAlice->TreeD();
 103   if (!dTree) {
 104     cerr<<"Input tree "<<name<<" not found !\n";
 105     return;
 106   }
 107
 108   TClonesArray *digitsSPD=new TClonesArray("AliITSdigitSPD",3000);
 109   dTree->SetBranchAddress("ITSDigitsSPD",&digitsSPD);
 110   TClonesArray *digitsSDD=new TClonesArray("AliITSdigitSDD",3000);
 111   dTree->SetBranchAddress("ITSDigitsSDD",&digitsSDD);
 112   TClonesArray *digitsSSD=new TClonesArray("AliITSdigitSSD",3000);
 113   dTree->SetBranchAddress("ITSDigitsSSD",&digitsSSD);
 114
 115   Int_t mmax=(Int_t)dTree->GetEntries();
 116
 117   out->cd();
 118
 119   sprintf(name,"TreeC_ITS_%d",fEvent);
 120   TTree cTree(name,"ITS clusters V2");
 121   TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
 122   cTree.Branch("Clusters",&clusters);
 123
 124   for (fI=0; fI<mmax; fI++) {
 125     dTree->GetEvent(fI);
 126
 127     if     (digitsSPD->GetEntriesFast()!=0)
 128                           FindClustersSPD(digitsSPD,clusters);
 129     else if(digitsSDD->GetEntriesFast()!=0)
 130                           FindClustersSDD(digitsSDD,clusters);
 131     else if(digitsSSD->GetEntriesFast()!=0)
 132                           FindClustersSSD(digitsSSD,clusters);
 133
 134     ncl+=clusters->GetEntriesFast();
 135
 136     cTree.Fill();
 137
 138     digitsSPD->Clear();
 139     digitsSDD->Clear();
 140     digitsSSD->Clear();
 141     clusters->Clear();
 142   }
 143   cTree.Write();
 144
 145   delete clusters;
 146
 147   delete digitsSPD;
 148   delete digitsSDD;
 149   delete digitsSSD;
 150
 151   //delete dTree;
 152
 153   cerr<<"Number of found clusters : "<<ncl<<endl;
 154
 155   savedir->cd();
 156 }
 157
 158 //**** Fast clusters *******************************
 159 #include "TParticle.h"
 160
 161 #include "AliITS.h"
 162 #include "AliITSmodule.h"
 163 #include "AliITSRecPoint.h"
 164 #include "AliITSsimulationFastPoints.h"
 165 #include "AliITSRecPoint.h"
 166
 167 static void CheckLabels(Int_t lab[3]) {
 168   //------------------------------------------------------------
 169   // Tries to find mother's labels
 170   //------------------------------------------------------------
 171     Int_t label=lab[0];
 172     if (label>=0) {
 173        TParticle *part=(TParticle*)gAlice->Particle(label);
 174        label=-3;
 175        while (part->P() < 0.005) {
 176           Int_t m=part->GetFirstMother();
 177           if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
 178           if (part->GetStatusCode()>0) {
 179              cerr<<"Primary momentum: "<<part->P()<<endl; break;
 180           }
 181           label=m;
 182           part=(TParticle*)gAlice->Particle(label);
 183         }
 184         if      (lab[1]<0) lab[1]=label;
 185         else if (lab[2]<0) lab[2]=label;
 186         else ;//cerr<<"CheckLabels : No empty labels !\n";
 187     }
 188 }
 189
 190 void AliITSclustererV2::RecPoints2Clusters
 191 (const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
 192   //------------------------------------------------------------
 193   // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS
 194   // subdetector indexed by idx
 195   //------------------------------------------------------------
 196   TClonesArray &cl=*clusters;
 197   Int_t ncl=points->GetEntriesFast();
 198   for (Int_t i=0; i<ncl; i++) {
 199     AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
 200     Float_t lp[5];
 201     lp[0]=-p->GetX()-fYshift[idx]; if (idx<=fLastSPD1) lp[0]*=-1; //SPD1
 202     lp[1]=p->GetZ()+fZshift[idx];
 203     lp[2]=p->GetSigmaX2();
 204     lp[3]=p->GetSigmaZ2();
 205     lp[4]=p->GetQ();
 206     Int_t lab[4];
 207     lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
 208     lab[3]=fNdet[idx];
 209     CheckLabels(lab);
 210     new (cl[i]) AliITSclusterV2(lab,lp);
 211   }
 212 }
 213
 214 void AliITSclustererV2::Hits2Clusters(const TFile *in, TFile *out) {
 215   //------------------------------------------------------------
 216   // This function creates ITS clusters
 217   //------------------------------------------------------------
 218   TDirectory *savedir=gDirectory;
 219
 220   if (!out->IsOpen()) {
 221     cerr<<"AliITSclustererV2::Hits2Clusters: output file not open !\n";
 222     return;
 223   }
 224
 225   if (!gAlice) {
 226      cerr<<"AliITSclustererV2::Hits2Clusters : gAlice==0 !\n";
 227      return;
 228   }
 229
 230   AliITS *its  = (AliITS*)gAlice->GetModule("ITS");
 231   if (!its) {
 232      cerr<<"AliITSclustererV2::Hits2Clusters : Can't find the ITS !\n";
 233      return;
 234   }
 235   AliITSgeom *geom=its->GetITSgeom();
 236   Int_t mmax=geom->GetIndexMax();
 237
 238   its->InitModules(-1,mmax);
 239   its->FillModules(gAlice->TreeH(),0);
 240
 241   out->cd();
 242
 243   Char_t name[100];
 244   sprintf(name,"TreeC_ITS_%d",fEvent);
 245   TTree cTree(name,"ITS clusters V2");
 246   TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
 247   cTree.Branch("Clusters",&clusters);
 248
 249   static TClonesArray *points=its->RecPoints();
 250   AliITSsimulationFastPoints sim;
 251   Int_t ncl=0;
 252   for (Int_t m=0; m<mmax; m++) {
 253     AliITSmodule *mod=its->GetModule(m);
 254     sim.CreateFastRecPoints(mod,m,gRandom);
 255
 256     RecPoints2Clusters(points, m, clusters);
 257     its->ResetRecPoints();
 258
 259     ncl+=clusters->GetEntriesFast();
 260     cTree.Fill();
 261     clusters->Clear();
 262   }
 263   cTree.Write();
 264
 265   cerr<<"Number of found fast clusters : "<<ncl<<endl;
 266
 267   delete clusters;
 268
 269   savedir->cd();
 270 }
 271
 272 //***********************************
 273
 274 #ifndef V1
 275
 276 void AliITSclustererV2::
 277 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
 278   //------------------------------------------------------------
 279   // returns an array of indices of digits belonging to the cluster
 280   // (needed when the segmentation is not regular)
 281   //------------------------------------------------------------
 282   if (n<200) idx[n++]=bins[k].GetIndex();
 283   bins[k].Use();
 284
 285   if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
 286   if (bins[k-1   ].IsNotUsed()) FindCluster(k-1   ,maxz,bins,n,idx);
 287   if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
 288   if (bins[k+1   ].IsNotUsed()) FindCluster(k+1   ,maxz,bins,n,idx);
 289   /*
 290   if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
 291   if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
 292   if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
 293   if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
 294   */
 295 }
 296
 297 void AliITSclustererV2::
 298 FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) {
 299   //------------------------------------------------------------
 300   // Actual SPD cluster finder
 301   //------------------------------------------------------------
 302   const Int_t kMAXBIN=(fNzSPD+2)*(fNySPD+2);
 303
 304   Int_t ndigits=digits->GetEntriesFast();
 305   AliBin *bins=new AliBin[kMAXBIN];
 306
 307   Int_t k;
 308   AliITSdigitSPD *d=0;
 309   for (k=0; k<ndigits; k++) {
 310      d=(AliITSdigitSPD*)digits->UncheckedAt(k);
 311      Int_t i=d->GetCoord2()+1;   //y
 312      Int_t j=d->GetCoord1()+1;
 313      bins[i*fNzSPD+j].SetIndex(k);
 314      bins[i*fNzSPD+j].SetMask(1);
 315   }
 316
 317   Int_t n=0; TClonesArray &cl=*clusters;
 318   for (k=0; k<kMAXBIN; k++) {
 319      if (!bins[k].IsNotUsed()) continue;
 320      Int_t ni=0, idx[200];
 321      FindCluster(k,fNzSPD,bins,ni,idx);
 322      if (ni==200) {cerr<<"SPD: Too big cluster !\n"; continue;}
 323
 324      Int_t lab[4];
 325      lab[0]=-2;
 326      lab[1]=-2;
 327      lab[2]=-2;
 328      lab[3]=fNdet[fI];
 329
 330      d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]);
 331      Int_t ymin=d->GetCoord2(),ymax=ymin;
 332      Int_t zmin=d->GetCoord1(),zmax=zmin;
 333      Float_t y=0.,z=0.,q=0.;
 334      for (Int_t l=0; l<ni; l++) {
 335         d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]);
 336
 337         if (ymin > d->GetCoord2()) ymin=d->GetCoord2();
 338         if (ymax < d->GetCoord2()) ymax=d->GetCoord2();
 339         if (zmin > d->GetCoord1()) zmin=d->GetCoord1();
 340         if (zmax < d->GetCoord1()) zmax=d->GetCoord1();
 341
 342         Int_t lab0=(d->GetTracks())[0];
 343         if (lab0>=0) {
 344           if (lab[0]<0) {
 345              lab[0]=lab0;
 346           } else if (lab[1]<0) {
 347             if (lab0!=lab[0]) lab[1]=lab0;
 348           } else if (lab[2]<0) {
 349             if (lab0!=lab[0])
 350             if (lab0!=lab[1]) lab[2]=lab0;
 351           }
 352         }
 353         Float_t qq=d->GetSignal();
 354         y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq;
 355      }
 356      y/=q; z/=q;
 357      y-=fHwSPD; z-=fHlSPD;
 358
 359      Float_t lp[5];
 360      lp[0]=-y-fYshift[fI]; if (fI<=fLastSPD1) lp[0]=-lp[0];
 361      lp[1]= z+fZshift[fI];
 362      lp[2]= fYpitchSPD*fYpitchSPD/12.;
 363      lp[3]= fZ1pitchSPD*fZ1pitchSPD/12.;
 364      //lp[4]= q;
 365      lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1);
 366
 367      //CheckLabels(lab);
 368      new (cl[n]) AliITSclusterV2(lab,lp); n++;
 369   }
 370
 371   delete bins;
 372 }
 373
 374 Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
 375   //------------------------------------------------------------
 376   //is this a local maximum ?
 377   //------------------------------------------------------------
 378   UShort_t q=bins[k].GetQ();
 379   if (q==1023) return kFALSE;
 380   if (bins[k-max].GetQ() > q) return kFALSE;
 381   if (bins[k-1  ].GetQ() > q) return kFALSE;
 382   if (bins[k+max].GetQ() > q) return kFALSE;
 383   if (bins[k+1  ].GetQ() > q) return kFALSE;
 384   if (bins[k-max-1].GetQ() > q) return kFALSE;
 385   if (bins[k+max-1].GetQ() > q) return kFALSE;
 386   if (bins[k+max+1].GetQ() > q) return kFALSE;
 387   if (bins[k-max+1].GetQ() > q) return kFALSE;
 388   return kTRUE;
 389 }
 390
 391 void AliITSclustererV2::
 392 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
 393   //------------------------------------------------------------
 394   //find local maxima
 395   //------------------------------------------------------------
 396   if (n<31)
 397   if (IsMaximum(k,max,b)) {
 398     idx[n]=k; msk[n]=(2<<n);
 399     n++;
 400   }
 401   b[k].SetMask(0);
 402   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
 403   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
 404   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
 405   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
 406 }
 407
 408 void AliITSclustererV2::
 409 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
 410   //------------------------------------------------------------
 411   //mark this peak
 412   //------------------------------------------------------------
 413   UShort_t q=bins[k].GetQ();
 414
 415   bins[k].SetMask(bins[k].GetMask()|m);
 416
 417   if (bins[k-max].GetQ() <= q)
 418      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
 419   if (bins[k-1  ].GetQ() <= q)
 420      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
 421   if (bins[k+max].GetQ() <= q)
 422      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
 423   if (bins[k+1  ].GetQ() <= q)
 424      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
 425 }
 426
 427 void AliITSclustererV2::
 428 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) {
 429   //------------------------------------------------------------
 430   //make cluster using digits of this peak
 431   //------------------------------------------------------------
 432   Float_t q=(Float_t)bins[k].GetQ();
 433   Int_t i=k/max, j=k-i*max;
 434
 435   c.SetQ(c.GetQ()+q);
 436   c.SetY(c.GetY()+i*q);
 437   c.SetZ(c.GetZ()+j*q);
 438   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
 439   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
 440
 441   bins[k].SetMask(0xFFFFFFFE);
 442
 443   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
 444   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
 445   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
 446   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
 447 }
 448
 449 void AliITSclustererV2::
 450 FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) {
 451   //------------------------------------------------------------
 452   // Actual SDD cluster finder
 453   //------------------------------------------------------------
 454   const Int_t kMAXBIN=(fNzSDD+2)*(fNySDD+2);
 455
 456   AliBin *bins[2];
 457   bins[0]=new AliBin[kMAXBIN];
 458   bins[1]=new AliBin[kMAXBIN];
 459
 460   AliITSdigitSDD *d=0;
 461   Int_t i, ndigits=digits->GetEntriesFast();
 462   for (i=0; i<ndigits; i++) {
 463      d=(AliITSdigitSDD*)digits->UncheckedAt(i);
 464      Int_t y=d->GetCoord2()+1;   //y
 465      Int_t z=d->GetCoord1()+1;   //z
 466      Int_t q=d->GetSignal();
 467      if (z <= fNzSDD) {
 468        bins[0][y*fNzSDD+z].SetQ(q);
 469        bins[0][y*fNzSDD+z].SetMask(1);
 470        bins[0][y*fNzSDD+z].SetIndex(i);
 471      } else {
 472        z-=fNzSDD;
 473        bins[1][y*fNzSDD+z].SetQ(q);
 474        bins[1][y*fNzSDD+z].SetMask(1);
 475        bins[1][y*fNzSDD+z].SetIndex(i);
 476      }
 477   }
 478
 479   Int_t ncl=0; TClonesArray &cl=*clusters;
 480   for (Int_t s=0; s<2; s++)
 481     for (i=0; i<kMAXBIN; i++) {
 482       if (bins[s][i].IsUsed()) continue;
 483       Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
 484       FindPeaks(i, fNzSDD, bins[s], idx, msk, npeaks);
 485
 486       if (npeaks>30) continue;
 487
 488       Int_t k,l;
 489       for (k=0; k<npeaks-1; k++){//mark adjacent peaks
 490         if (idx[k] < 0) continue; //this peak is already removed
 491         for (l=k+1; l<npeaks; l++) {
 492            if (idx[l] < 0) continue; //this peak is already removed
 493            Int_t ki=idx[k]/fNzSDD, kj=idx[k] - ki*fNzSDD;
 494            Int_t li=idx[l]/fNzSDD, lj=idx[l] - li*fNzSDD;
 495            Int_t di=TMath::Abs(ki - li);
 496            Int_t dj=TMath::Abs(kj - lj);
 497            if (di>1 || dj>1) continue;
 498            if (bins[s][idx[k]].GetQ() > bins[s][idx[l]].GetQ()) {
 499               msk[l]=msk[k];
 500               idx[l]*=-1;
 501            } else {
 502               msk[k]=msk[l];
 503               idx[k]*=-1;
 504               break;
 505            }
 506         }
 507       }
 508
 509       for (k=0; k<npeaks; k++) {
 510         MarkPeak(TMath::Abs(idx[k]), fNzSDD, bins[s], msk[k]);
 511       }
 512
 513       for (k=0; k<npeaks; k++) {
 514          if (idx[k] < 0) continue; //removed peak
 515          AliITSclusterV2 c;
 516          MakeCluster(idx[k], fNzSDD, bins[s], msk[k], c);
 517
 518          //if (c.GetQ() < 200) continue; //noise cluster
 519
 520          /*
 521          Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
 522          Float_t w=par->GetPadPitchWidth(sec);
 523          c.SetSigmaY2((s2 + 1./12.)*w*w);
 524          if (s2 != 0.) {
 525            c.SetSigmaY2(c.GetSigmaY2()*0.108);
 526            if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
 527          }
 528
 529          s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
 530          w=par->GetZWidth();
 531          c.SetSigmaZ2((s2 + 1./12.)*w*w);
 532          if (s2 != 0.) {
 533            c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
 534            if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
 535          }
 536          */
 537
 538          c.SetSigmaY2(0.0030*0.0030);
 539          c.SetSigmaZ2(0.0020*0.0020);
 540          c.SetDetectorIndex(fNdet[fI]);
 541
 542          Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ();
 543          y/=q; z/=q;
 544
 545          y=(y-0.5)*fYpitchSDD;
 546          y-=fHwSDD;
 547          y-=fYoffSDD;  //delay ?
 548          if (s) y=-y;
 549
 550          z=(z-0.5)*fZpitchSDD;
 551          z-=fHlSDD;
 552
 553          y=-y-fYshift[fI];
 554          z= z+fZshift[fI];
 555          c.SetY(y);
 556          c.SetZ(z);
 557
 558          c.SetQ(q/20.);  //to be consistent with the SSD charges
 559
 560          AliBin *b=&bins[s][idx[k]];
 561          d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 562          Int_t l0=(d->GetTracks())[0];
 563          if (l0<0) {
 564            b=&bins[s][idx[k]-1];
 565            if (b->GetQ()>0) {
 566              d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 567              l0=(d->GetTracks())[0];
 568            }
 569          }
 570          if (l0<0) {
 571            b=&bins[s][idx[k]+1];
 572            if (b->GetQ()>0) {
 573              d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 574              l0=(d->GetTracks())[0];
 575            }
 576          }
 577          if (l0<0) {
 578            b=&bins[s][idx[k]-fNzSDD];
 579           if (b->GetQ()>0) {
 580              d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 581              l0=(d->GetTracks())[0];
 582            }
 583          }
 584          if (l0<0) {
 585            b=&bins[s][idx[k]+fNzSDD];
 586            if (b->GetQ()>0) {
 587              d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 588              l0=(d->GetTracks())[0];
 589            }
 590          }
 591
 592          if (l0<0) {
 593            b=&bins[s][idx[k]+fNzSDD+1];
 594            if (b->GetQ()>0) {
 595              d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 596              l0=(d->GetTracks())[0];
 597            }
 598          }
 599          if (l0<0) {
 600            b=&bins[s][idx[k]+fNzSDD-1];
 601            if (b->GetQ()>0) {
 602              d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 603              l0=(d->GetTracks())[0];
 604            }
 605          }
 606          if (l0<0) {
 607            b=&bins[s][idx[k]-fNzSDD+1];
 608            if (b->GetQ()>0) {
 609              d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 610              l0=(d->GetTracks())[0];
 611            }
 612          }
 613          if (l0<0) {
 614            b=&bins[s][idx[k]-fNzSDD-1];
 615            if (b->GetQ()>0) {
 616              d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex());
 617              l0=(d->GetTracks())[0];
 618            }
 619          }
 620
 621          {
 622            Int_t lab[3];
 623            lab[0]=(d->GetTracks())[0];
 624            lab[1]=(d->GetTracks())[1];
 625            lab[2]=(d->GetTracks())[2];
 626            //CheckLabels(lab);
 627            c.SetLabel(lab[0],0);
 628            c.SetLabel(lab[1],1);
 629            c.SetLabel(lab[2],2);
 630          }
 631
 632          new (cl[ncl]) AliITSclusterV2(c); ncl++;
 633       }
 634     }
 635
 636   delete[] bins[0];
 637   delete[] bins[1];
 638 }
 639
 640 void AliITSclustererV2::
 641 FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) {
 642   //------------------------------------------------------------
 643   // Actual SSD cluster finder
 644   //------------------------------------------------------------
 645   Int_t smax=digits->GetEntriesFast();
 646   if (smax==0) return;
 647
 648   const Int_t MAX=1000;
 649   Int_t np=0, nn=0;
 650   Ali1Dcluster pos[MAX], neg[MAX];
 651   Float_t y=0., q=0., qmax=0.;
 652   Int_t lab[4]={-2,-2,-2,-2};
 653
 654   TClonesArray &cl=*clusters;
 655
 656   AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0);
 657   q += d->GetSignal();
 658   y += d->GetCoord2()*d->GetSignal();
 659   qmax=d->GetSignal();
 660   lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
 661   Int_t curr=d->GetCoord2();
 662   Int_t flag=d->GetCoord1();
 663   Int_t *n=&nn;
 664   Ali1Dcluster *c=neg;
 665   Int_t nd=0;
 666   for (Int_t s=1; s<smax; s++) {
 667       d=(AliITSdigitSSD*)digits->UncheckedAt(s);
 668       Int_t strip=d->GetCoord2();
 669       if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
 670          c[*n].SetY(y/q);
 671          c[*n].SetQ(q);
 672          c[*n].SetNd(nd);
 673          c[*n].SetLabels(lab);
 674          //Split suspiciously big cluster
 675          if (nd>3) {
 676             c[*n].SetY(y/q-0.5*nd);
 677             c[*n].SetQ(0.5*q);
 678             (*n)++;
 679             if (*n==MAX) {
 680               cerr<<
 681               "AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n";
 682               return;
 683             }
 684             c[*n].SetY(y/q+0.5*nd);
 685             c[*n].SetQ(0.5*q);
 686             c[*n].SetNd(nd);
 687             c[*n].SetLabels(lab);
 688          }
 689          (*n)++;
 690          if (*n==MAX) {
 691           cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n";
 692           return;
 693          }
 694          y=q=qmax=0.;
 695          nd=0;
 696          lab[0]=lab[1]=lab[2]=-2;
 697          if (flag!=d->GetCoord1()) { n=&np; c=pos; }
 698       }
 699       flag=d->GetCoord1();
 700       q += d->GetSignal();
 701       y += d->GetCoord2()*d->GetSignal();
 702       nd++;
 703       if (d->GetSignal()>qmax) {
 704          qmax=d->GetSignal();
 705          lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
 706       }
 707       curr=strip;
 708   }
 709   c[*n].SetY(y/q);
 710   c[*n].SetQ(q);
 711   c[*n].SetNd(nd);
 712   c[*n].SetLabels(lab);
 713   //Split suspiciously big cluster
 714   if (nd>3) {
 715      c[*n].SetY(y/q-0.5*nd);
 716      c[*n].SetQ(0.5*q);
 717      (*n)++;
 718      if (*n==MAX) {
 719         cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n";
 720         return;
 721      }
 722      c[*n].SetY(y/q+0.5*nd);
 723      c[*n].SetQ(0.5*q);
 724      c[*n].SetNd(nd);
 725      c[*n].SetLabels(lab);
 726   }
 727   (*n)++;
 728   if (*n==MAX) {
 729      cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n";
 730      return;
 731   }
 732
 733   Float_t tanp=fTanP, tann=fTanN;
 734   if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;}
 735
 736   Int_t idet=fNdet[fI];
 737   Int_t ncl=0;
 738   for (Int_t i=0; i<np; i++) {
 739     //Float_t dq_min=1.e+33;
 740     Float_t ybest=1000,zbest=1000,qbest=0;
 741     Float_t yp=pos[i].GetY()*fYpitchSSD;
 742     for (Int_t j=0; j<nn; j++) {
 743       //if (pos[i].fTracks[0] != neg[j].fTracks[0]) continue;
 744       Float_t yn=neg[j].GetY()*fYpitchSSD;
 745       Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
 746       Float_t yt=yn + tann*zt;
 747       zt-=fHlSSD; yt-=fHwSSD;
 748       if (TMath::Abs(yt)<fHwSSD+0.01)
 749       if (TMath::Abs(zt)<fHlSSD+0.01) {
 750       //if (TMath::Abs(pos[i].GetQ()-neg[j].GetQ())<dq_min) {
 751         //dq_min=TMath::Abs(pos[i].GetQ()-neg[j].GetQ());
 752         ybest=yt; zbest=zt;
 753         qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
 754
 755         lab[0]=pos[i].GetLabel(0);
 756         lab[1]=pos[i].GetLabel(1);
 757         lab[2]=neg[i].GetLabel(0);
 758         lab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
 759         Float_t lp[5];
 760         lp[0]=-ybest-fYshift[fI];
 761         lp[1]= zbest+fZshift[fI];
 762         lp[2]=0.0025*0.0025;  //SigmaY2
 763         lp[3]=0.110*0.110;  //SigmaZ2
 764         if (pos[i].GetNd()+neg[j].GetNd() > 4) {
 765            lp[2]*=9;
 766            lp[3]*=9;
 767         }
 768         lp[4]=qbest;        //Q
 769
 770         //CheckLabels(lab);
 771         new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++;
 772       }
 773     }
 774     /*
 775     if (ybest<100) {
 776        lab[3]=idet;
 777        Float_t lp[5];
 778        lp[0]=-ybest-fYshift[fI];
 779        lp[1]= zbest+fZshift[fI];
 780        lp[2]=0.002*0.002;  //SigmaY2
 781        lp[3]=0.080*0.080;  //SigmaZ2
 782        lp[4]=qbest;        //Q
 783        //
 784        new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++;
 785     }
 786     */
 787   }
 788 }
 789
 790 #else   //V1
 791
 792 #include "AliITSDetType.h"
 793
 794 #include "AliITSsegmentationSPD.h"
 795 #include "AliITSClusterFinderSPD.h"
 796
 797 #include "AliITSresponseSDD.h"
 798 #include "AliITSsegmentationSDD.h"
 799 #include "AliITSClusterFinderSDD.h"
 800
 801 #include "AliITSsegmentationSSD.h"
 802 #include "AliITSClusterFinderSSD.h"
 803
 804
 805 void AliITSclustererV2::
 806 FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) {
 807   //------------------------------------------------------------
 808   // Actual SPD cluster finding based on AliITSClusterFinderSPD
 809   //------------------------------------------------------------
 810   static AliITS *its=(AliITS*)gAlice->GetModule("ITS");
 811   static TClonesArray *points=its->RecPoints();
 812   static AliITSsegmentationSPD *seg=
 813          (AliITSsegmentationSPD *)its->DetType(0)->GetSegmentationModel();
 814   static AliITSClusterFinderSPD cf(seg, (TClonesArray*)digits, points);
 815
 816   cf.FindRawClusters(fI);
 817   RecPoints2Clusters(points, fI, clusters);
 818   its->ResetRecPoints();
 819
 820 }
 821
 822 void AliITSclustererV2::
 823 FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) {
 824   //------------------------------------------------------------
 825   // Actual SDD cluster finding based on AliITSClusterFinderSDD
 826   //------------------------------------------------------------
 827   static AliITS *its=(AliITS*)gAlice->GetModule("ITS");
 828   static TClonesArray *points=its->RecPoints();
 829   static AliITSresponseSDD *resp=
 830         (AliITSresponseSDD *)its->DetType(1)->GetResponseModel();
 831   static AliITSsegmentationSDD *seg=
 832          (AliITSsegmentationSDD *)its->DetType(1)->GetSegmentationModel();
 833   static AliITSClusterFinderSDD
 834          cf(seg,resp,(TClonesArray*)digits,its->ClustersAddress(1));
 835
 836   cf.FindRawClusters(fI);
 837   Int_t nc=points->GetEntriesFast();
 838   while (nc--) { //To be consistent with the SSD cluster charges
 839      AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(nc);
 840      p->SetQ(p->GetQ/12.);
 841   }
 842   RecPoints2Clusters(points, fI, clusters);
 843   its->ResetClusters(1);
 844   its->ResetRecPoints();
 845
 846 }
 847
 848 void AliITSclustererV2::
 849 FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) {
 850   //------------------------------------------------------------
 851   // Actual SSD cluster finding based on AliITSClusterFinderSSD
 852   //------------------------------------------------------------
 853   static AliITS *its=(AliITS*)gAlice->GetModule("ITS");
 854   static TClonesArray *points=its->RecPoints();
 855   static AliITSsegmentationSSD *seg=
 856          (AliITSsegmentationSSD *)its->DetType(2)->GetSegmentationModel();
 857   static AliITSClusterFinderSSD cf(seg,(TClonesArray*)digits);
 858
 859   cf.FindRawClusters(fI);
 860   RecPoints2Clusters(points, fI, clusters);
 861   its->ResetRecPoints();
 862
 863 }
 864
 865 #endif