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