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