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     for(ndig=0; ndig<ndigits; ndig++) {
 107         dig= (AliITSdigitSPD*)fDigits->UncheckedAt(ndig);
 108         digx[digitcount] = dig->fCoord2+1;  //starts at 1
 109         digz[digitcount] = dig->fCoord1+1;  //starts at 1
 110         digtr1[digitcount] = dig->fTracks[0];
 111         digtr2[digitcount] = -3;
 112         digtr3[digitcount] = -3;
 113         i=1;
 114         while(digtr1[digitcount]==dig->fTracks[i] && i<dig->GetNTracks()) i++;
 115         if(i<dig->GetNTracks()){
 116             digtr2[digitcount] = dig->fTracks[i];
 117             while(digtr1[digitcount]==dig->fTracks[i] &&
 118                   i<dig->GetNTracks()) i++;
 119             if(i<dig->GetNTracks()) digtr3[digitcount] = dig->fTracks[i];
 120         } // end if
 121         digtr4[digitcount] = dig->fSignal;
 122         digitcount++;
 123     } // end for ndig
 124     ClusterFinder(digitcount,digx,digz,digtr1,digtr2,digtr3,digtr4,
 125                   nclus,xcenterl,zcenterl,errxcenter,errzcenter,
 126                   tr1clus, tr2clus, tr3clus,module);
 127     DigitToPoint(nclus,xcenterl,zcenterl,errxcenter,errzcenter,
 128                  tr1clus, tr2clus, tr3clus);
 129     delete[] digx;
 130     delete[] digz;
 131     delete[] digtr1;
 132     delete[] digtr2;
 133     delete[] digtr3;
 134     delete[] digtr4;
 135     delete[] xcenterl;
 136     delete[] zcenterl;
 137     delete[] errxcenter;
 138     delete[] errzcenter;
 139     delete[] tr1clus;
 140     delete[] tr2clus;
 141     delete[] tr3clus;
 142 }
 143 //----------------------------------------------------------------------
 144 void AliITSClusterFinderSPD::ClusterFinder(Int_t ndigits,Int_t digx[],
 145                                            Int_t digz[],Int_t digtr1[],
 146                                            Int_t digtr2[],Int_t digtr3[],
 147                                            Int_t digtr4[],Int_t &nclus,
 148                                            Float_t xcenter[],Float_t zcenter[],
 149                                            Float_t errxcenter[],
 150                                            Float_t errzcenter[],
 151                                            Int_t tr1clus[],Int_t tr2clus[],
 152                                            Int_t tr3clus[],Int_t module){
 153     // Search for clusters of fired pixels (digits). Two digits are linked
 154     // inside a cluster if they are countiguous both in row or column
 155     // direction.  Diagonal digits are not linked.
 156     // xcenter, ycenter, zcenter are the coordinates of the center
 157     // of each found cluster, calculated from the averaging the corresponding
 158     // coordinate of the center of the linked digits. The coordinates are
 159     // given in the local reference sistem.
 160     // errxcenter, errycenter, errzcenter are the errors associated to
 161     // the corresponding average.
 162     Int_t   if1, min, max, nd;
 163     Int_t   x1, z1, t1, t2, t3, t4;
 164     Int_t   ndx, ndz, ndxmin=0, ndxmax=0, ndzmin=0, ndzmax=0;
 165     Float_t dx, dz;
 166     Int_t   i,k,ipos=0;
 167     Float_t xdum, zdum;
 168     Int_t   kmax, sigmax;
 169     Float_t deltax, deltaz;
 170     Float_t ndig;
 171     Int_t   numberd = 10000;
 172     Int_t   *ifpad  = new Int_t[numberd];
 173     Int_t   *xpad   = new Int_t[numberd];
 174     Int_t   *zpad   = new Int_t[numberd];
 175     Int_t   *tr1pad = new Int_t[numberd];
 176     Int_t   *tr2pad = new Int_t[numberd];
 177     Int_t   *tr3pad = new Int_t[numberd];
 178     Int_t   *tr4pad = new Int_t[numberd];
 179     Int_t   *iclus  = new Int_t[numberd];
 180
 181     nclus=1;
 182     for (i=0; i < ndigits ; i++){
 183         ifpad[i] = -1;
 184         iclus[i] = 0;
 185     } // end for i
 186
 187     ifpad[0]=0;
 188     for (i=0; i < ndigits-1 ; i++) {
 189         if ( ifpad[i] == -1 ) {
 190             nclus++;
 191             ipos++;
 192             ifpad[i]=nclus-1;
 193         } // end if ipad[i]
 194         for (Int_t j=i+1 ; j < ndigits ; j++)  {
 195             if (ifpad[j]== -1 ) {
 196                 dx = TMath::Abs(digx[i]-digx[j]);
 197                 dz = TMath::Abs(digz[i]-digz[j]);
 198                 // if ( ( dx+dz )==1 )  //clusters are not diagonal
 199                 if(( dx+dz )==1 || (dx==1 && dz==1)){
 200                     //diagonal clusters allowed
 201                     ipos++;
 202                     ifpad[j] = ifpad[i];
 203
 204                     x1         = digx[j];
 205                     z1         = digz[j];
 206                     digx[j]    = digx[ipos];
 207                     digz[j]    = digz[ipos];
 208                     digx[ipos] = x1;
 209                     digz[ipos] = z1;
 210
 211                     t1 = digtr1[j];
 212                     t2 = digtr2[j];
 213                     t3 = digtr3[j];
 214                     t4 = digtr4[j];
 215                     digtr1[j] = digtr1[ipos];
 216                     digtr2[j] = digtr2[ipos];
 217                     digtr3[j] = digtr3[ipos];
 218                     digtr4[j] = digtr4[ipos];
 219                     digtr1[ipos] = t1;
 220                     digtr2[ipos] = t2;
 221                     digtr3[ipos] = t3;
 222                     digtr4[ipos] = t4;
 223
 224                     if1 = ifpad[j];
 225                     ifpad[j] = ifpad[ipos];
 226                     ifpad[ipos] = if1;
 227                 } // end dx+dx...
 228             }// end if ifpad[j]== -1
 229         } // end for j
 230     }//end loop on digits
 231     if ( ifpad[ndigits-1] == -1 ) {
 232         nclus++;
 233         ifpad[ndigits-1]=nclus-1;
 234     } // end if ifpad[ndigits-1] == -1
 235
 236     for (i=0 ; i < ndigits ; i++) iclus[ifpad[i]]++;
 237
 238     min=0;
 239     max=0;
 240     // loop on found clusters
 241     for (i=0 ; i < nclus ; i++){
 242         min = max;
 243         max += iclus[i];
 244         deltax = fSegmentation->Dpx(0);
 245         if (iclus[i]!=1){
 246             //cluster with more than one digit
 247             nd=iclus[i];
 248             ndig=(Float_t) nd;
 249             Int_t count=0;
 250             for (k=min;k<min+nd;k++){
 251                 xpad[count] = digx[k];
 252                 zpad[count] = digz[k];
 253
 254                 tr1pad[count] = digtr1[k];
 255                 tr2pad[count] = digtr2[k];
 256                 tr3pad[count] = digtr3[k];
 257                 tr4pad[count] = digtr4[k];
 258
 259                 count++;
 260             } // end for k
 261             ndxmin = xpad[TMath::LocMin(nd,xpad)];
 262             ndxmax = xpad[TMath::LocMax(nd,xpad)];
 263             ndzmin = zpad[TMath::LocMin(nd,zpad)];
 264             ndzmax = zpad[TMath::LocMax(nd,zpad)];
 265             ndx = ndxmax - ndxmin+1;
 266             ndz = ndzmax - ndzmin+1;
 267
 268             // calculate x and z coordinates of the center of the cluster
 269             fSegmentation->GetPadCxz(digx[min],digz[min]-1,xdum, zdum);
 270
 271             if (ndx == 1) {
 272                 xcenter[i] = xdum;
 273             }else{
 274                 xcenter[i] = 0.;
 275                 for (k=0;k<nd;k++) {
 276                     fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
 277                     xcenter[i] += (xdum / nd);
 278                 } // end for k
 279             } // end if ndx
 280
 281             if (ndz == 1) {
 282                 zcenter[i] = zdum;
 283             } else {
 284                 zcenter[i] = 0.;
 285                 for (k=0;k<nd;k++) {
 286                     fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
 287                     zcenter[i] += (zdum / nd);
 288                 } // end for k
 289             } // end if ndz
 290
 291             // error on points in x and z directions
 292
 293             if (ndx == 1) {
 294                 errxcenter[i] = deltax / TMath::Sqrt(12.);
 295             } else {
 296                 errxcenter[i] = 0.;
 297                 for (k=0;k<nd;k++){
 298                     fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
 299                     errxcenter[i] += ((xdum-xcenter[i])*(xdum-xcenter[i]))/
 300                                      (nd*(nd-1));
 301                 } // end for k
 302                 errxcenter[i] = TMath::Sqrt(errxcenter[i]);
 303             } // end if ndx
 304             if (ndz == 1) {
 305                 deltaz = fSegmentation->Dpz(digz[min]);
 306                 errzcenter[i] = deltaz / TMath::Sqrt(12.);
 307             } else {
 308                 errzcenter[i] = 0.;
 309                 for (k=0;k<nd;k++){
 310                     fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
 311                     errzcenter[i] += ((zdum-zcenter[i])*(zdum-zcenter[i]))/
 312                                      (nd*(nd-1));
 313                 } // end for k
 314                 errzcenter[i] = TMath::Sqrt(errzcenter[i]);
 315             } // end if ndz
 316             // take three track numbers for the cluster
 317             // choose the track numbers of the digit with higher signal
 318             kmax = 0;
 319             sigmax = 0;
 320             for (k=0;k<nd;k++){
 321                 if(tr4pad[k] > sigmax){
 322                     sigmax = tr4pad[k];
 323                     kmax   = k;
 324                 } // end if tr4pad[k]
 325             } // end for k
 326             if(sigmax != 0) {
 327                 tr1clus[i]= tr1pad[kmax];
 328                 tr2clus[i]= tr2pad[kmax];
 329                 tr3clus[i]= tr3pad[kmax];
 330             } else {
 331                 tr1clus[i]= -2;
 332                 tr2clus[i]= -2;
 333                 tr3clus[i]= -2;
 334             } // end if sigmax
 335         } else {
 336             // cluster with single digit
 337             ndig= 1.;
 338             ndx = 1;
 339             ndz = 1;
 340             fSegmentation->GetPadCxz(digx[min],digz[min]-1,xdum,zdum);
 341             xcenter[i] = xdum;
 342             zcenter[i] = zdum;
 343             tr1clus[i]=digtr1[min];
 344             tr2clus[i]=digtr2[min];
 345             tr3clus[i]=digtr3[min];
 346             deltaz = fSegmentation->Dpz(digz[min]);
 347             errxcenter[i] = deltax / TMath::Sqrt(12.);
 348             errzcenter[i] = deltaz / TMath::Sqrt(12.);
 349         } // end if iclus[i]
 350
 351         // store the cluster information to the AliITSRawCLusterSPD object
 352         static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 353
 354         //put the cluster center in local reference frame of the detector
 355         // and in microns
 356         xcenter[i] = xcenter[i] - fSegmentation->Dx()/2.;
 357         zcenter[i] = zcenter[i] - fSegmentation->Dz()/2.;
 358
 359         AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i], //f
 360                                                              xcenter[i], //f
 361                                                              ndig, //f
 362                                                              ndz,ndx, //ii
 363                                                              ndxmin,ndxmax,//ii
 364                                                              (Float_t) ndzmin,
 365                                                              (Float_t) ndzmax,
 366                                                              0,module); //ii
 367         iTS->AddCluster(0,clust);
 368         delete clust;
 369     }//end loop on clusters
 370     delete[] ifpad;
 371     delete[] xpad ;
 372     delete[] zpad ;
 373     delete[] iclus;
 374     delete[] tr1pad;
 375     delete[] tr2pad;
 376     delete[] tr3pad;
 377     delete[] tr4pad;
 378 }
 379 //______________________________________________________----------------
 380 void AliITSClusterFinderSPD::DigitToPoint(Int_t nclus,
 381                                           Float_t *xcenter,Float_t *zcenter,
 382                                           Float_t *errxcenter,
 383                                           Float_t *errzcenter,
 384                                           Int_t *tr1clus, Int_t *tr2clus,
 385                                           Int_t *tr3clus){
 386     // A point is associated to each cluster of SPD digits. The points
 387     // and their associated errors are stored in the file galiceSP.root.
 388     Float_t l[3],xg,zg;
 389     const Float_t kconv = 1.0e-4; // micron -> cm
 390
 391     // get rec points
 392     static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 393     for (Int_t i=0; i<nclus; i++){
 394         l[0] = kconv*xcenter[i];
 395         l[1] = kconv*fSegmentation->Dy()/2.;
 396         l[2] = kconv*zcenter[i];
 397
 398         xg = l[0];
 399         zg = l[2];
 400
 401         Float_t sigma2x = (kconv*errxcenter[i]) * (kconv*errxcenter[i]);
 402         Float_t sigma2z = (kconv*errzcenter[i]) * (kconv*errzcenter[i]);
 403         AliITSRecPoint rnew;
 404         rnew.SetX(xg);
 405         rnew.SetZ(zg);
 406         rnew.SetQ(1.);
 407         rnew.SetdEdX(0.);
 408         rnew.SetSigmaX2(sigma2x);
 409         rnew.SetSigmaZ2(sigma2z);
 410         rnew.fTracks[0]=tr1clus[i];
 411         rnew.fTracks[1]=tr2clus[i];
 412         rnew.fTracks[2]=tr3clus[i];
 413         iTS->AddRecPoint(rnew);
 414     } // end for i
 415 }