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