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