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