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