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 "AliITSMapA1.h"
  19 #include "AliITS.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 ClassImp(AliITSClusterFinderSPD)
  30
  31 //----------------------------------------------------------
  32 AliITSClusterFinderSPD::AliITSClusterFinderSPD
  33 (AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)
  34 {
  35   // constructor
  36     fSegmentation=seg;
  37     fDigits=digits;
  38     fClusters=recp;
  39     fNclusters= fClusters->GetEntriesFast();
  40     SetDx();
  41     SetDz();
  42     fMap=new AliITSMapA1(fSegmentation,fDigits);
  43     SetNCells();
  44 }
  45
  46 //_____________________________________________________________________________
  47 AliITSClusterFinderSPD::AliITSClusterFinderSPD()
  48 {
  49   // constructor
  50   fSegmentation=0;
  51   fDigits=0;
  52   fClusters=0;
  53   fNclusters=0;
  54   fMap=0;
  55   SetDx();
  56   SetDz();
  57   SetNCells();
  58
  59 }
  60
  61 //_____________________________________________________________________________
  62 AliITSClusterFinderSPD::~AliITSClusterFinderSPD()
  63 {
  64   // destructor
  65   if (fMap) delete fMap;
  66
  67
  68 }
  69 //__________________________________________________________________________
  70 AliITSClusterFinderSPD::AliITSClusterFinderSPD(const AliITSClusterFinderSPD &source){
  71   //     Copy Constructor
  72   if(&source == this) return;
  73   this->fClusters = source.fClusters ;
  74   this->fNclusters = source.fNclusters ;
  75   this->fMap = source.fMap ;
  76   this->fDz = source.fDz ;
  77   this->fDx = source.fDx ;
  78   this->fMinNCells = source.fMinNCells ;
  79   return;
  80 }
  81
  82 //_________________________________________________________________________
  83 AliITSClusterFinderSPD&
  84   AliITSClusterFinderSPD::operator=(const AliITSClusterFinderSPD &source) {
  85   //    Assignment operator
  86   if(&source == this) return *this;
  87   this->fClusters = source.fClusters ;
  88   this->fNclusters = source.fNclusters ;
  89   this->fMap = source.fMap ;
  90   this->fDz = source.fDz ;
  91   this->fDx = source.fDx ;
  92   this->fMinNCells = source.fMinNCells ;
  93   return *this;
  94 }
  95
  96 //_____________________________________________________________________________
  97 void AliITSClusterFinderSPD::SetMap()
  98 {
  99   // set map
 100
 101   if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);
 102
 103 }
 104
 105 //_____________________________________________________________________________
 106
 107 void AliITSClusterFinderSPD::Find1DClusters()
 108 {
 109   // Find one dimensional clusters, i.e.
 110   // in r*phi(x) direction for each colunm in z direction
 111
 112   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 113
 114   // retrieve the parameters
 115   Int_t fNofPixels = fSegmentation->Npz();
 116   Int_t fMaxNofSamples = fSegmentation->Npx();
 117
 118   // read in digits -> do not apply threshold
 119   // signal in fired pixels is always 1
 120
 121   fMap->FillMap();
 122
 123   Int_t nofFoundClusters = 0;
 124
 125   Int_t k,it,m;
 126   for(k=0;k<fNofPixels;k++) {
 127
 128     Int_t mmax = 10;  // a size of the window for the cluster finding
 129
 130     for(it=0;it<fMaxNofSamples;it++) {
 131
 132       Int_t lclx = 0;
 133       Int_t xstart = 0;
 134       Int_t xstop = 0;
 135       Int_t id = 0;
 136       Int_t ilcl =0;
 137
 138       for(m=0;m<mmax;m++) {  // find the cluster inside the window
 139         id = it+m;
 140         if(id >= fMaxNofSamples) break;    // ! no possible for the fadc
 141
 142         if(fMap->TestHit(k,id) == kUnused) {   // start of the cluster
 143           lclx += 1;
 144           if(lclx == 1) xstart = id;
 145
 146         }
 147
 148         if(lclx > 0 && fMap->TestHit(k,id) == kEmpty) {
 149           // end of cluster if a gap exists
 150           xstop = id-1;
 151           ilcl = 1;
 152           break;
 153         }
 154
 155       }   //  end of m-loop
 156
 157       if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
 158       // continue the "it" loop
 159
 160       if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
 161
 162       if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {
 163                                    // cluster end is outside of the window,
 164         mmax += 5;                 // increase mmax and repeat the cluster
 165                                    // finding
 166         it -= 1;
 167       }
 168
 169       if(id >= fMaxNofSamples && lclx > 0) {  // the x row finished but
 170         xstop = fMaxNofSamples - 1;           // the end cluster exists
 171         ilcl = 1;
 172       }
 173
 174       // ---  Calculate z and x coordinates for one dimensional clusters
 175
 176       if(ilcl == 1) {         // new cluster exists
 177         it = id;
 178         mmax = 10;
 179             nofFoundClusters++;
 180             Float_t clusterCharge = 0.;
 181             Float_t zpitch = fSegmentation->Dpz(k+1);
 182             Float_t clusterZ, dummyX;
 183             Int_t dummy=0;
 184             fSegmentation->GetPadCxz(dummy,k,dummyX,clusterZ);
 185             Float_t zstart = clusterZ - 0.5*zpitch;
 186             Float_t zstop = clusterZ + 0.5*zpitch;
 187             Float_t clusterX = 0.;
 188             Int_t xstartfull = xstart;
 189             Int_t xstopfull = xstop;
 190             Int_t clusterSizeX = lclx;
 191             Int_t clusterSizeZ = 1;
 192
 193             Int_t its;
 194             for(its=xstart; its<=xstop; its++) {
 195               Int_t firedpixel=0;
 196               if (fMap->GetHitIndex(k,its)>=0) firedpixel=1;
 197               clusterCharge += firedpixel;
 198               clusterX +=its + 0.5;
 199             }
 200             Float_t fRphiPitch = fSegmentation->Dpx(dummy);
 201             clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x
 202
 203
 204             //printf("ClusterZ ClusterX %f %f \n",clusterZ, clusterX);
 205
 206             // Int_t nclusters = fClusters->GetEntriesFast();
 207             //              cout << nclusters << " clusters" << endl;
 208             //            cout<< "Create point"<<endl;
 209
 210             // Write the points (coordinates and some cluster information) to the
 211             // AliITSRawClusterSPD object
 212
 213             AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(clusterZ,clusterX,clusterCharge,clusterSizeZ,clusterSizeX,xstart,xstop,xstartfull,xstopfull,zstart,zstop,k);
 214             // fClusters->Add(point);
 215             iTS->AddCluster(0,clust);
 216             delete clust;
 217             //              cout << "Cluster at Ladder: " << fLadder << ", Detector: " <<fDetector<<endl;
 218
 219             // cout<<" end of cluster finding for Z pixel "<<endl;
 220
 221       }    // new cluster (ilcl=1)
 222     } // X direction loop (it)
 223   } // Z direction loop (k)
 224
 225   //fMap->ClearMap();
 226   return;
 227
 228 }
 229
 230 //_____________________________________________________________________________
 231 void  AliITSClusterFinderSPD::GroupClusters()
 232 {
 233   // Find two dimensional clusters, i.e. group one dimensional clusters
 234   // into two dimensional ones (go both in x and z directions).
 235
 236   // get number of clusters for this module
 237   Int_t nofClusters = fClusters->GetEntriesFast();
 238   nofClusters -= fNclusters;
 239
 240   AliITSRawClusterSPD *clusterI;
 241   AliITSRawClusterSPD *clusterJ;
 242
 243   Int_t *label=new Int_t[nofClusters];
 244   Int_t i,j;
 245   for(i=0; i<nofClusters; i++) label[i] = 0;
 246   for(i=0; i<nofClusters; i++) {
 247     if(label[i] != 0) continue;
 248     for(j=i+1; j<nofClusters; j++) {
 249       if(label[j] != 0) continue;
 250       clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
 251       clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
 252       Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
 253       if(pair) {
 254
 255         //    if((clusterI->XStop() == clusterJ->XStart()-1)||(clusterI->XStart()==clusterJ->XStop()+1)) cout<<"!! Diagonal cluster"<<endl;
 256         /*
 257       cout << "clusters " << i << "," << j << " before grouping" << endl;
 258               clusterI->PrintInfo();
 259               clusterJ->PrintInfo();
 260         */
 261         clusterI->Add(clusterJ);
 262         //        cout << "remove cluster " << j << endl;
 263         label[j] = 1;
 264         fClusters->RemoveAt(j);
 265         /*
 266           cout << "cluster  " << i << " after grouping" << endl;
 267           clusterI->PrintInfo();
 268         */
 269       }  // pair
 270     } // J clusters
 271     label[i] = 1;
 272   } // I clusters
 273   fClusters->Compress();
 274   //  Int_t totalNofClusters = fClusters->GetEntriesFast();
 275   //  cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
 276
 277   delete [] label;
 278
 279   return;
 280
 281
 282 }
 283 //_____________________________________________________________________________
 284
 285 void AliITSClusterFinderSPD::TracksInCluster()
 286 {
 287
 288   // Find tracks creating one cluster
 289
 290   // get number of clusters for this module
 291   Int_t nofClusters = fClusters->GetEntriesFast();
 292   nofClusters -= fNclusters;
 293
 294   //printf("nofClusters fNclusters %d %d\n",nofClusters,fNclusters);   //commented 18-09-00
 295
 296   Int_t i, ix, iz, jx, jz, xstart, xstop, zstart, zstop, nclx, nclz;
 297   //  Int_t signal, track0, track1, track2;
 298   const Int_t trmax = 100;
 299   Int_t cltracks[trmax], itr, tracki, ii, is, js, ie, ntr, tr0, tr1, tr2;
 300
 301   for(i=0; i<nofClusters; i++) {
 302     ii = 0;
 303     memset(cltracks,-1,sizeof(int)*trmax);
 304     tr0=tr1=tr2=-3;
 305
 306     AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
 307     //printf("clusterI %p\n",clusterI);               //commented 18-09-00
 308
 309     nclx = clusterI->NclX();
 310     nclz = clusterI->NclZ();
 311     xstart = clusterI->XStartf();
 312     xstop = clusterI->XStopf();
 313     zstart = clusterI->Zend()-nclz+1;
 314     zstop = clusterI->Zend();
 315
 316     Int_t ind;
 317
 318      for(iz=0; iz<nclz; iz++) {
 319          jz = zstart + iz;
 320        for(ix=0; ix<nclx; ix++) {
 321          jx = xstart + ix;
 322          ind = fMap->GetHitIndex(jz,jx);
 323          //printf("ind %d\n",ind);                              // commented 18-09-00
 324          // this part must be wrong - the index of the digit
 325          // can be zero but not -1 !!! comment out this part
 326          // and replace it - see below
 327          /*
 328          if(ind == 0 && iz >= 0 && ix > 0) {
 329           continue;
 330          }
 331          if(ind == 0 && iz > 0 && ix >= 0) {
 332           continue;
 333          }
 334          if(ind == 0 && iz == 0 && ix == 0 && i > 0) {
 335           continue;
 336          }
 337          */
 338          // change the conditions above to :
 339          if (ind < 0) continue;
 340
 341         AliITSdigitSPD *dig = (AliITSdigitSPD*)fMap->GetHit(jz,jx);
 342         if (!dig) printf("SPD: something wrong ! dig %p\n",dig);
 343         /*
 344          signal=dig->fSignal;
 345          track0=dig->fTracks[0];
 346          track1=dig->fTracks[1];
 347          track2=dig->fTracks[2];
 348         */
 349           for(itr=0; itr<3; itr++) {
 350             tracki = dig->fTracks[itr];
 351             //printf("tracki %d\n",tracki);                        //commented 18-09-00
 352             if(tracki >= 0) {
 353               ii += 1;
 354              cltracks[ii-1] = tracki;
 355             }
 356           }
 357        } // ix pixel
 358      }  // iz pixel
 359
 360      for(is=0; is<trmax; is++) {
 361          if(cltracks[is]<0) continue;
 362        for(js=is+1; js<trmax; js++) {
 363          if(cltracks[js]<0) continue;
 364          if(cltracks[js]==cltracks[is]) cltracks[js]=-5;
 365        }
 366      }
 367
 368      ntr = 0;
 369      for(ie=0; ie<trmax; ie++) {
 370        if(cltracks[ie] >= 0) {
 371         ntr=ntr+1;
 372         if(ntr==1) tr0=cltracks[ie];
 373         if(ntr==2) tr1=cltracks[ie];
 374         if(ntr==3) tr2=cltracks[ie];
 375        }
 376      }
 377      // if delta ray only
 378      if(ntr == 0) ntr = 1;
 379
 380      clusterI->SetNTracks(ntr);
 381      clusterI->SetTracks(tr0,tr1,tr2);
 382
 383   } // I cluster
 384
 385 }
 386 //_____________________________________________________________________________
 387
 388 void AliITSClusterFinderSPD::GetRecPoints()
 389 {
 390   // get rec points
 391   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 392
 393   // get number of clusters for this module
 394   Int_t nofClusters = fClusters->GetEntriesFast();
 395   nofClusters -= fNclusters;
 396   const Float_t kconv = 1.0e-4;
 397   const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
 398   const Float_t kRMSz = 120.0*kconv; // resolution for 425 micron pixel
 399
 400   Float_t spdLength = fSegmentation->Dz();
 401   Float_t spdWidth = fSegmentation->Dx();
 402
 403   Int_t i;
 404   Int_t track0, track1, track2;
 405
 406   for(i=0; i<nofClusters; i++) {
 407
 408     AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
 409     clusterI->GetTracks(track0, track1, track2);
 410     AliITSRecPoint rnew;
 411
 412     rnew.SetX((clusterI->X() - spdWidth/2)*kconv);
 413     rnew.SetZ((clusterI->Z() - spdLength/2)*kconv);
 414     rnew.SetQ(1.);
 415     rnew.SetdEdX(0.);
 416     rnew.SetSigmaX2(kRMSx*kRMSx);
 417     rnew.SetSigmaZ2(kRMSz*kRMSz);
 418     rnew.fTracks[0]=track0;
 419     rnew.fTracks[1]=track1;
 420     rnew.fTracks[2]=track2;
 421     iTS->AddRecPoint(rnew);
 422   } // I clusters
 423
 424   fMap->ClearMap();
 425
 426 }
 427 //_____________________________________________________________________________
 428
 429 void AliITSClusterFinderSPD::FindRawClusters()
 430 {
 431   // find raw clusters
 432   Find1DClusters();
 433   GroupClusters();
 434   TracksInCluster();
 435   GetRecPoints();
 436
 437 }
 438