ITS/AliITSClusterFinderSPD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 //
  16 //  Cluster finder
  17 //  for Silicon pixels
  18 //
  19 //
  20 #include "AliITSClusterFinderSPD.h"
  21 #include "AliITS.h"
  22 #include "AliITSdigitSPD.h"
  23 #include "AliITSRawClusterSPD.h"
  24 #include "AliITSRecPoint.h"
  25 #include "AliITSresponseSPD.h"
  26 #include "AliITSsegmentationSPD.h"
  27 #include "AliRun.h"
  28
  29 //#define DEBUG
  30
  31 ClassImp(AliITSClusterFinderSPD)
  32
  33 //______________________________________________________________________
  34 AliITSClusterFinderSPD::AliITSClusterFinderSPD():AliITSClusterFinder(),
  35 fDz(0.0),
  36 fDx(0.0),
  37 fMinNCells(0){
  38     // constructor
  39 }
  40 //----------------------------------------------------------
  41 AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSsegmentation *seg,
  42                                                AliITSresponse *res):
  43 AliITSClusterFinder(seg,res),
  44 fDz(0.0),
  45 fDx(0.0),
  46 fMinNCells(0){
  47     // constructor
  48
  49     SetDx();
  50     SetDz();
  51 }
  52 //----------------------------------------------------------
  53 AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSsegmentation *seg,
  54                                                TClonesArray *digits,
  55                                                TClonesArray *recp):
  56 AliITSClusterFinder(seg,0),
  57 fDz(0.0),
  58 fDx(0.0),
  59 fMinNCells(0){
  60     // constructor
  61
  62     SetDigits(digits);
  63     SetClusters(recp);
  64     SetDx();
  65     SetDz();
  66 }
  67 //_____________________________________________________________________
  68 AliITSClusterFinderSPD::AliITSClusterFinderSPD(
  69          const AliITSClusterFinderSPD &source): AliITSClusterFinder(source){
  70     //     Copy Constructor
  71
  72     if(&source == this) return;
  73     this->fDz        = source.fDz;
  74     this->fDx        = source.fDx;
  75     this->fMinNCells = source.fMinNCells;
  76     return;
  77 }
  78 //______________________________________________________________________
  79 AliITSClusterFinderSPD& AliITSClusterFinderSPD::operator=(
  80                                        const AliITSClusterFinderSPD &source) {
  81     //    Assignment operator
  82
  83     if(&source == this) return *this;
  84     this->fDz        = source.fDz;
  85     this->fDx        = source.fDx;
  86     this->fMinNCells = source.fMinNCells;
  87     return *this;
  88 }
  89 //______________________________________________________________________
  90 void AliITSClusterFinderSPD::FindRawClusters(Int_t module){
  91     // input of Cluster Finder
  92     Int_t   digitcount  = 0;
  93     Int_t   numberd     = 100000;
  94     Int_t   *digx       = new Int_t[numberd];
  95     Int_t   *digz       = new Int_t[numberd];
  96     Int_t   *digtr1     = new Int_t[numberd];
  97     Int_t   *digtr2     = new Int_t[numberd];
  98     Int_t   *digtr3     = new Int_t[numberd];
  99     Int_t   *digtr4     = new Int_t[numberd];
 100     //  output of Cluster Finder
 101     Int_t   numberc     = 10000;
 102     Double_t *xcenterl   = new Double_t[numberc];
 103     Double_t *zcenterl   = new Double_t[numberc];
 104     Double_t *errxcenter = new Double_t[numberc];
 105     Double_t *errzcenter = new Double_t[numberc];
 106     Int_t   *tr1clus    = new Int_t[numberc];
 107     Int_t   *tr2clus    = new Int_t[numberc];
 108     Int_t   *tr3clus    = new Int_t[numberc];
 109     Int_t   nclus;
 110
 111     SetModule(module);
 112     digitcount=0;
 113     Int_t ndigits = Digits()->GetEntriesFast();
 114     if (!ndigits) return;
 115
 116     AliITSdigitSPD *dig;
 117     Int_t ndig=0,i;
 118     /*
 119     if(GetDebug(4)){
 120         cout << "SPD - FindRawclusters"<<endl;
 121         scanf("%d",&ndig);
 122     } // end if GetDebug
 123     */
 124     for(ndig=0; ndig<ndigits; ndig++) {
 125         dig= (AliITSdigitSPD*)GetDigit(ndig);
 126         digx[digitcount] = dig->GetCoord2()+1;  //starts at 1
 127         digz[digitcount] = dig->GetCoord1()+1;  //starts at 1
 128         digtr1[digitcount] = dig->GetTrack(0);
 129         digtr2[digitcount] = -3;
 130         digtr3[digitcount] = -3;
 131         if(GetDebug(5)){
 132             cout << "digtr1["<<digitcount <<"]="<<digtr1[digitcount];
 133             cout << " fTracks["<<0<<"]="<<dig->GetTrack(0)<<": ";
 134         } // end if GetDebug
 135         i=1;
 136         while(digtr1[digitcount]==dig->GetTrack(i) && i<dig->GetNTracks()) i++;
 137         if(GetDebug(5)) cout << " fTracks["<<i<<"]="<<dig->GetTrack(i);
 138         if(i<dig->GetNTracks()){
 139             digtr2[digitcount] = dig->GetTrack(i);
 140             if(GetDebug(5)) cout<<"digtr2["<<digitcount <<"]="
 141                                 <<digtr2[digitcount]<<": ";
 142             while((digtr1[digitcount]==dig->GetTrack(i) ||
 143                    digtr2[digitcount]==dig->GetTrack(i))&&
 144                   i<=dig->GetNTracks()) i++;
 145             if(i<dig->GetNTracks()) digtr3[digitcount] = dig->GetTrack(i);
 146             if(GetDebug(5)){
 147                 cout << " fTracks["<<i<<"]=";
 148                 if(i<dig->GetNTracks()) cout <<dig->GetTrack(i);
 149                 cout << "digtr3["<<digitcount <<"]="<<digtr3[digitcount];
 150             } // end if GetDebug
 151         } // end if
 152         if(GetDebug(4)) cout<<endl;
 153         digtr4[digitcount] = dig->GetSignal();
 154         digitcount++;
 155     } // end for ndig
 156     ClusterFinder(digitcount,digx,digz,digtr1,digtr2,digtr3,digtr4,
 157                   nclus,xcenterl,zcenterl,errxcenter,errzcenter,
 158                   tr1clus, tr2clus, tr3clus);
 159     DigitToPoint(nclus,xcenterl,zcenterl,errxcenter,errzcenter,
 160                  tr1clus, tr2clus, tr3clus);
 161     delete[] digx;
 162     delete[] digz;
 163     delete[] digtr1;
 164     delete[] digtr2;
 165     delete[] digtr3;
 166     delete[] digtr4;
 167     delete[] xcenterl;
 168     delete[] zcenterl;
 169     delete[] errxcenter;
 170     delete[] errzcenter;
 171     delete[] tr1clus;
 172     delete[] tr2clus;
 173     delete[] tr3clus;
 174 }
 175 //----------------------------------------------------------------------
 176 void AliITSClusterFinderSPD::ClusterFinder(Int_t ndigits,Int_t digx[],
 177                                            Int_t digz[],Int_t digtr1[],
 178                                            Int_t digtr2[],Int_t digtr3[],
 179                                            Int_t digtr4[],Int_t &nclus,
 180                                            Double_t xcenter[],Double_t zcenter[],
 181                                            Double_t errxcenter[],
 182                                            Double_t errzcenter[],
 183                                            Int_t tr1clus[],Int_t tr2clus[],
 184                                            Int_t tr3clus[]){
 185     // Search for clusters of fired pixels (digits). Two digits are linked
 186     // inside a cluster if they are countiguous both in row or column
 187     // direction.  Diagonal digits are not linked.
 188     // xcenter, ycenter, zcenter are the coordinates of the center
 189     // of each found cluster, calculated from the averaging the corresponding
 190     // coordinate of the center of the linked digits. The coordinates are
 191     // given in the local reference sistem.
 192     // errxcenter, errycenter, errzcenter are the errors associated to
 193     // the corresponding average.
 194     Int_t   if1, min, max, nd;
 195     Int_t   x1, z1, t1, t2, t3, t4;
 196     Int_t   ndx, ndz, ndxmin=0, ndxmax=0, ndzmin=0, ndzmax=0;
 197     Double_t dx, dz;
 198     Int_t   i,k,ipos=0;
 199     Float_t xdum, zdum;
 200     Int_t   kmax, sigmax;
 201     Double_t deltax, deltaz;
 202     Double_t ndig;
 203     Int_t   numberd = 10000;
 204     Int_t   *ifpad  = new Int_t[numberd];
 205     Int_t   *xpad   = new Int_t[numberd];
 206     Int_t   *zpad   = new Int_t[numberd];
 207     Int_t   *tr1pad = new Int_t[numberd];
 208     Int_t   *tr2pad = new Int_t[numberd];
 209     Int_t   *tr3pad = new Int_t[numberd];
 210     Int_t   *tr4pad = new Int_t[numberd];
 211     Int_t   *iclus  = new Int_t[numberd];
 212
 213     nclus=1;
 214     for (i=0; i < ndigits ; i++){
 215         ifpad[i] = -1;
 216         iclus[i] = 0;
 217     } // end for i
 218
 219     ifpad[0]=0;
 220     for (i=0; i < ndigits-1 ; i++) {
 221         if ( ifpad[i] == -1 ) {
 222             nclus++;
 223             ipos++;
 224             ifpad[i]=nclus-1;
 225         } // end if ipad[i]
 226         for (Int_t j=i+1 ; j < ndigits ; j++)  {
 227             if (ifpad[j]== -1 ) {
 228                 dx = TMath::Abs(digx[i]-digx[j]);
 229                 dz = TMath::Abs(digz[i]-digz[j]);
 230                 // if ( ( dx+dz )==1 )  //clusters are not diagonal
 231                 if(( dx+dz )==1 || (dx==1 && dz==1)){
 232                     //diagonal clusters allowed
 233                     ipos++;
 234                     ifpad[j] = ifpad[i];
 235
 236                     x1         = digx[j];
 237                     z1         = digz[j];
 238                     digx[j]    = digx[ipos];
 239                     digz[j]    = digz[ipos];
 240                     digx[ipos] = x1;
 241                     digz[ipos] = z1;
 242
 243                     t1 = digtr1[j];
 244                     t2 = digtr2[j];
 245                     t3 = digtr3[j];
 246                     t4 = digtr4[j];
 247                     digtr1[j] = digtr1[ipos];
 248                     digtr2[j] = digtr2[ipos];
 249                     digtr3[j] = digtr3[ipos];
 250                     digtr4[j] = digtr4[ipos];
 251                     digtr1[ipos] = t1;
 252                     digtr2[ipos] = t2;
 253                     digtr3[ipos] = t3;
 254                     digtr4[ipos] = t4;
 255
 256                     if1 = ifpad[j];
 257                     ifpad[j] = ifpad[ipos];
 258                     ifpad[ipos] = if1;
 259                 } // end dx+dx...
 260             }// end if ifpad[j]== -1
 261         } // end for j
 262     }//end loop on digits
 263     if ( ifpad[ndigits-1] == -1 ) {
 264         nclus++;
 265         ifpad[ndigits-1]=nclus-1;
 266     } // end if ifpad[ndigits-1] == -1
 267
 268     for (i=0 ; i < ndigits ; i++) iclus[ifpad[i]]++;
 269
 270     min=0;
 271     max=0;
 272     // loop on found clusters
 273     for (i=0 ; i < nclus ; i++){
 274         min = max;
 275         max += iclus[i];
 276         deltax = GetSeg()->Dpx(0);
 277         if (iclus[i]!=1){
 278             //cluster with more than one digit
 279             nd=iclus[i];
 280             ndig=(Double_t) nd;
 281             Int_t count=0;
 282             for (k=min;k<min+nd;k++){
 283                 xpad[count] = digx[k];
 284                 zpad[count] = digz[k];
 285
 286                 tr1pad[count] = digtr1[k];
 287                 tr2pad[count] = digtr2[k];
 288                 tr3pad[count] = digtr3[k];
 289                 tr4pad[count] = digtr4[k];
 290
 291                 count++;
 292             } // end for k
 293             ndxmin = xpad[TMath::LocMin(nd,xpad)];
 294             ndxmax = xpad[TMath::LocMax(nd,xpad)];
 295             ndzmin = zpad[TMath::LocMin(nd,zpad)];
 296             ndzmax = zpad[TMath::LocMax(nd,zpad)];
 297             ndx = ndxmax - ndxmin+1;
 298             ndz = ndzmax - ndzmin+1;
 299
 300             // calculate x and z coordinates of the center of the cluster
 301             GetSeg()->GetPadCxz(digx[min],digz[min]-1,xdum, zdum);
 302
 303             if (ndx == 1) {
 304                 xcenter[i] = xdum;
 305             }else{
 306                 xcenter[i] = 0.;
 307                 for (k=0;k<nd;k++) {
 308                     GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
 309                     xcenter[i] += (xdum / nd);
 310                 } // end for k
 311             } // end if ndx
 312
 313             if (ndz == 1) {
 314                 zcenter[i] = zdum;
 315             } else {
 316                 zcenter[i] = 0.;
 317                 for (k=0;k<nd;k++) {
 318                     GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
 319                     zcenter[i] += (zdum / nd);
 320                 } // end for k
 321             } // end if ndz
 322
 323             // error on points in x and z directions
 324
 325             if (ndx == 1) {
 326                 errxcenter[i] = deltax / TMath::Sqrt(12.);
 327             } else {
 328                 errxcenter[i] = 0.;
 329                 for (k=0;k<nd;k++){
 330                     GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
 331                     errxcenter[i] += ((xdum-xcenter[i])*(xdum-xcenter[i]))/
 332                         (nd*(nd-1));
 333                 } // end for k
 334                 errxcenter[i] = TMath::Sqrt(errxcenter[i]);
 335             } // end if ndx
 336             if (ndz == 1) {
 337                 deltaz = GetSeg()->Dpz(digz[min]);
 338                 errzcenter[i] = deltaz / TMath::Sqrt(12.);
 339             } else {
 340                 errzcenter[i] = 0.;
 341                 for (k=0;k<nd;k++){
 342                     GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
 343                     errzcenter[i] += ((zdum-zcenter[i])*(zdum-zcenter[i]))/
 344                         (nd*(nd-1));
 345                 } // end for k
 346                 errzcenter[i] = TMath::Sqrt(errzcenter[i]);
 347             } // end if ndz
 348             // take three track numbers for the cluster
 349             // choose the track numbers of the digit with higher signal
 350             kmax = 0;
 351             sigmax = 0;
 352             for (k=0;k<nd;k++){
 353                 if(tr4pad[k] > sigmax){
 354                     sigmax = tr4pad[k];
 355                     kmax   = k;
 356                 } // end if tr4pad[k]
 357             } // end for k
 358             if(sigmax != 0) {
 359                 tr1clus[i]= tr1pad[kmax];
 360                 tr2clus[i]= tr2pad[kmax];
 361                 tr3clus[i]= tr3pad[kmax];
 362             } else {
 363                 tr1clus[i]= -2;
 364                 tr2clus[i]= -2;
 365                 tr3clus[i]= -2;
 366             } // end if sigmax
 367         } else {
 368             // cluster with single digit
 369             ndig= 1.;
 370             ndx = 1;
 371             ndz = 1;
 372             GetSeg()->GetPadCxz(digx[min],digz[min]-1,xdum,zdum);
 373             xcenter[i] = xdum;
 374             zcenter[i] = zdum;
 375             tr1clus[i]=digtr1[min];
 376             tr2clus[i]=digtr2[min];
 377             tr3clus[i]=digtr3[min];
 378             deltaz = GetSeg()->Dpz(digz[min]);
 379             errxcenter[i] = deltax / TMath::Sqrt(12.);
 380             errzcenter[i] = deltaz / TMath::Sqrt(12.);
 381         } // end if iclus[i]
 382
 383         // store the cluster information to the AliITSRawCLusterSPD object
 384         static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 385
 386         //put the cluster center in local reference frame of the detector
 387         // and in microns
 388         xcenter[i] = xcenter[i] - GetSeg()->Dx()/2.;
 389         zcenter[i] = zcenter[i] - GetSeg()->Dz()/2.;
 390
 391         AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i], //f
 392                                                              xcenter[i], //f
 393                                                              ndig, //f
 394                                                              ndz,ndx, //ii
 395                                                              ndxmin,ndxmax,//ii
 396                                                              (Double_t) ndzmin,
 397                                                              (Double_t) ndzmax,
 398                                                              0,GetModule());
 399         iTS->AddCluster(0,clust);
 400         delete clust;
 401     }//end loop on clusters
 402     delete[] ifpad;
 403     delete[] xpad ;
 404     delete[] zpad ;
 405     delete[] iclus;
 406     delete[] tr1pad;
 407     delete[] tr2pad;
 408     delete[] tr3pad;
 409     delete[] tr4pad;
 410 }
 411 //______________________________________________________----------------
 412 void AliITSClusterFinderSPD::DigitToPoint(Int_t nclus,
 413                                           Double_t *xcenter,Double_t *zcenter,
 414                                           Double_t *errxcenter,
 415                                           Double_t *errzcenter,
 416                                           Int_t *tr1clus, Int_t *tr2clus,
 417                                           Int_t *tr3clus){
 418     // A point is associated to each cluster of SPD digits. The points
 419     // and their associated errors are stored in the file galiceSP.root.
 420     Double_t l[3],xg,zg;
 421     const Double_t kconv = 1.0e-4; // micron -> cm
 422
 423     // get rec points
 424     static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 425     for (Int_t i=0; i<nclus; i++){
 426         l[0] = kconv*xcenter[i];
 427         l[1] = kconv*GetSeg()->Dy()/2.;
 428         l[2] = kconv*zcenter[i];
 429
 430         xg = l[0];
 431         zg = l[2];
 432
 433         Double_t sigma2x = (kconv*errxcenter[i]) * (kconv*errxcenter[i]);
 434         Double_t sigma2z = (kconv*errzcenter[i]) * (kconv*errzcenter[i]);
 435         AliITSRecPoint rnew;
 436         rnew.SetX(xg);
 437         rnew.SetZ(zg);
 438         rnew.SetQ(1.);
 439         rnew.SetdEdX(0.);
 440         rnew.SetSigmaX2(sigma2x);
 441         rnew.SetSigmaZ2(sigma2z);
 442         rnew.fTracks[0]=tr1clus[i];
 443         rnew.fTracks[1]=tr2clus[i];
 444         rnew.fTracks[2]=tr3clus[i];
 445         iTS->AddRecPoint(rnew);
 446     } // end for i
 447 }