ITS/AliITSClusterFinderSPDdubna.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 /* $Id$ */
  17
  18 //  Cluster Finder for
  19 //  Silicon
  20 //  Pixel
  21 //  developed by the Dubna group
  22 //
  23 #include <Riostream.h>
  24 #include "AliITSClusterFinderSPDdubna.h"
  25 #include "AliITSMapA1.h"
  26 #include "AliITS.h"
  27 #include "AliITSdigit.h"
  28 #include "AliITSRawCluster.h"
  29 #include "AliITSRecPoint.h"
  30 #include "AliITSsegmentation.h"
  31 #include "AliRun.h"
  32
  33 //#define DEBUG
  34
  35 ClassImp(AliITSClusterFinderSPDdubna)
  36
  37 //----------------------------------------------------------------------
  38 AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna(
  39            AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp){
  40     // constructor
  41
  42     fSegmentation = seg;
  43     fDigits       = digits;
  44     fClusters     = recp;
  45     fNclusters    = 0;
  46     fMap          = new AliITSMapA1(fSegmentation,fDigits);
  47     SetDx();
  48     SetDz();
  49     SetNCells();
  50 }
  51 //_______________________________________________________________________
  52 AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna(){
  53     // constructor
  54
  55     fSegmentation = 0;
  56     fDigits       = 0;
  57     fClusters     = 0;
  58     fNclusters    = 0;
  59     fMap          = 0;
  60     SetDx();
  61     SetDz();
  62     SetNCells();
  63 }
  64 //______________________________________________________________________
  65 AliITSClusterFinderSPDdubna &AliITSClusterFinderSPDdubna::operator=(const
  66                                           AliITSClusterFinderSPDdubna &s){
  67     // The = operator for the class AliITSClusterFinderSPDdugna
  68     // Inputs:
  69     //    AliITSClusterFinderSPDdubna @s  The AliITSClusterFinderSPDdubna
  70     //                                    class to be copy constructed.
  71     // Outputs:
  72     //    none.
  73     // Returned:
  74     //    this, a copy of the class s.
  75
  76     this->fNclusters = s.fNclusters;
  77     this->fDz        = s.fDz;
  78     this->fDx        = s.fDx;
  79     this->fMinNCells = s.fMinNCells;
  80     this->fClusters  = new TClonesArray(*(s.fClusters));
  81     return *this;
  82 }
  83 //______________________________________________________________________
  84 AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna(
  85       const AliITSClusterFinderSPDdubna &s): AliITSClusterFinder(s){
  86     // The Copy constructortor the class AliITSClusterFinderSPDdugna
  87     // It calles the = operator.
  88     // Inputs:
  89     //    AliITSClusterFinderSPDdubna @s  The AliITSClusterFinderSPDdubna
  90     //                                    class to be copy constructed.
  91     // Outputs:
  92     //    none.
  93     // Returned:
  94     //    none.
  95
  96     *this = s;
  97     return;
  98 }
  99 //______________________________________________________________________
 100 AliITSClusterFinderSPDdubna::~AliITSClusterFinderSPDdubna(){
 101     // destructor
 102
 103     if (fMap) delete fMap;
 104 }
 105 //______________________________________________________________________
 106 void AliITSClusterFinderSPDdubna::Find1DClusters(Int_t mod){
 107     // Find one dimensional clusters, i.e.
 108     // in r*phi(x) direction for each colunm in z direction
 109
 110     static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 111
 112     // retrieve the parameters
 113     Int_t fNofPixels     = fSegmentation->Npz();
 114     Int_t fMaxNofSamples = fSegmentation->Npx();
 115
 116     // read in digits -> do not apply threshold
 117     // signal in fired pixels is always 1
 118     fMap->FillMap();
 119
 120     Int_t nofFoundClusters = 0;
 121
 122     Int_t k,it,m;
 123     for(k=0;k<fNofPixels;k++) {
 124         Int_t mmax = 10;  // a size of the window for the cluster finding
 125         for(it=0;it<fMaxNofSamples;it++) {
 126             Int_t lclx   = 0;
 127             Int_t xstart = 0;
 128             Int_t xstop  = 0;
 129             Int_t id     = 0;
 130             Int_t ilcl   = 0;
 131             for(m=0;m<mmax;m++) {  // find the cluster inside the window
 132                 id = it+m;
 133                 if(id >= fMaxNofSamples) break;    // ! no possible for the fadc
 134                 if(fMap->TestHit(k,id) == kUnused) {   // start of the cluster
 135                     lclx += 1;
 136                     if(lclx == 1) xstart = id;
 137                 } // end if fMap
 138                 if(lclx > 0 && fMap->TestHit(k,id) == kEmpty) {
 139                     // end of cluster if a gap exists
 140                     xstop = id-1;
 141                     ilcl = 1;
 142                     break;
 143                 } // end if lclx >0 ...
 144             }   //  end of m-loop
 145             if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
 146                                                 // continue the "it" loop
 147             if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
 148             if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {
 149                 // cluster end is outside of the window,
 150                 mmax += 5;     // increase mmax and repeat the cluster finding
 151                 it -= 1;
 152             } // end if id<
 153             if(id >= fMaxNofSamples && lclx > 0) {  // the x row finished but
 154                 xstop = fMaxNofSamples - 1;           // the end cluster exists
 155                 ilcl = 1;
 156             }  // end if id>=
 157             // ---  Calculate z and x coordinates for one dimensional clusters
 158             if(ilcl == 1) {         // new cluster exists
 159                 it = id;
 160                 mmax = 10;
 161                 nofFoundClusters++;
 162                 Float_t clusterCharge = 0.;
 163                 Float_t zpitch        = fSegmentation->Dpz(k+1);
 164                 Int_t   dummy         = 0;
 165                 Float_t clusterZ, dummyX;
 166                 fSegmentation->GetPadCxz(dummy,k,dummyX,clusterZ);
 167                 Float_t zstart        = clusterZ - 0.5*zpitch;
 168                 Float_t zstop         = clusterZ + 0.5*zpitch;
 169                 Float_t clusterX      = 0.;
 170                 Int_t   clusterSizeX  = lclx;
 171                 Int_t   clusterSizeZ  = 1;
 172                 Int_t   its;
 173                 for(its=xstart; its<=xstop; its++) {
 174                     Int_t firedpixel=0;
 175                     if (fMap->GetHitIndex(k,its)>=0) firedpixel=1;
 176                     clusterCharge += firedpixel;
 177                     clusterX      += its + 0.5;
 178                 } // end for its
 179                 Float_t fRphiPitch = fSegmentation->Dpx(dummy);
 180                 clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x
 181                 // Write the points (coordinates and some cluster information)
 182                 // to the AliITSRawClusterSPD object
 183
 184                 AliITSRawClusterSPD clust(clusterZ,clusterX,clusterCharge,
 185                                           clusterSizeZ,clusterSizeX,
 186                                           xstart,xstop,zstart,zstop,k,mod);
 187                 iTS->AddCluster(0,&clust);
 188             }    // new cluster (ilcl=1)
 189         } // X direction loop (it)
 190     } // Z direction loop (k)
 191     //fMap->ClearMap();
 192     return;
 193 }
 194 //______________________________________________________________________
 195 void  AliITSClusterFinderSPDdubna::GroupClusters(){
 196     // Find two dimensional clusters, i.e. group one dimensional clusters
 197     // into two dimensional ones (go both in x and z directions).
 198     // get number of clusters for this module
 199
 200     Int_t nofClusters = fClusters->GetEntriesFast();
 201     nofClusters -= fNclusters;
 202     AliITSRawClusterSPD *clusterI;
 203     AliITSRawClusterSPD *clusterJ;
 204     Int_t *label=new Int_t[nofClusters];
 205     Int_t i,j;
 206     for(i=0; i<nofClusters; i++) label[i] = 0;
 207     for(i=0; i<nofClusters; i++) {
 208         if(label[i] != 0) continue;
 209         for(j=i+1; j<nofClusters; j++) {
 210             if(label[j] != 0) continue;
 211             clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
 212             clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
 213             Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
 214             if(pair) {
 215 #ifndef DEBUG
 216                 cout << "clusters " << i << "," << j << " before grouping" << endl;
 217                 clusterI->PrintInfo();
 218                 clusterJ->PrintInfo();
 219 #endif
 220                 clusterI->Add(clusterJ);
 221 #ifndef DEBUG
 222                 cout << "remove cluster " << j << endl;
 223 #endif
 224                 label[j] = 1;
 225                 fClusters->RemoveAt(j);
 226 #ifndef DEBUG
 227                 cout << "cluster  " << i << " after grouping" << endl;
 228                 clusterI->PrintInfo();
 229 #endif
 230             }  // pair
 231         } // J clusters
 232         label[i] = 1;
 233     } // I clusters
 234     fClusters->Compress();
 235 #ifndef DEBUG
 236     Int_t totalNofClusters = fClusters->GetEntriesFast();
 237     cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
 238 #endif
 239     delete [] label;
 240     return;
 241 }
 242 //______________________________________________________________________
 243 void AliITSClusterFinderSPDdubna::TracksInCluster(){
 244     // Find tracks creating one cluster
 245     // get number of clusters for this module
 246     Int_t nofClusters = fClusters->GetEntriesFast();
 247     Int_t i, ix, iz, jx, jz, xstart, xstop, zstart, zstop, nclx, nclz;
 248     const Int_t kTrmax = 100;
 249     Int_t cltracks[kTrmax], itr, tracki, ii, is, js, ie, ntr, tr0, tr1, tr2;
 250
 251     nofClusters -= fNclusters;
 252     for(i=0; i<nofClusters; i++) {
 253         ii = 0;
 254         memset(cltracks,-1,sizeof(int)*kTrmax);
 255         tr0=tr1=tr2=-1;
 256         AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
 257         nclx = clusterI->NclX();
 258         nclz = clusterI->NclZ();
 259         xstart = clusterI->XStartf();
 260         xstop = clusterI->XStopf();
 261         zstart = clusterI->Zend()-nclz+1;
 262         zstop = clusterI->Zend();
 263         Int_t ind;
 264         for(iz=0; iz<nclz; iz++) {
 265             jz = zstart + iz;
 266             for(ix=0; ix<nclx; ix++) {
 267                 jx = xstart + ix;
 268                 ind = fMap->GetHitIndex(jz,jx);
 269                 if(ind < 0) continue;
 270                 if(ind == 0 && iz >= 0 && ix > 0) continue;
 271                 if(ind == 0 && iz > 0 && ix >= 0) continue;
 272                 if(ind == 0 && iz == 0 && ix == 0 && i > 0) continue;
 273                 AliITSdigitSPD *dig = (AliITSdigitSPD*)fMap->GetHit(jz,jx);
 274                 for(itr=0; itr<3; itr++) {
 275                     tracki = dig->GetTrack(i);
 276                     if(tracki >= 0) {
 277                         ii += 1;
 278                         if(ii < 99) cltracks[ii-1] = tracki;
 279                     } // end if tracki >-0
 280                 } // end for itr
 281             } // ix pixel
 282         }  // iz pixel
 283         for(is=0; is<kTrmax; is++) {
 284             if(cltracks[is]<0) continue;
 285             for(js=is+1; js<kTrmax; js++) {
 286                 if(cltracks[js]<0) continue;
 287                 if(cltracks[js]==cltracks[is]) cltracks[js]=-5;
 288             } // end for js
 289         } // end for is
 290         ntr = 0;
 291         for(ie=0; ie<kTrmax; ie++) {
 292             if(cltracks[ie] >= 0) {
 293                 ntr=ntr+1;
 294                 if(ntr==1) tr0=cltracks[ie];
 295                 if(ntr==2) tr1=cltracks[ie];
 296                 if(ntr==3) tr2=cltracks[ie];
 297             } // end if cltracks[ie] >-=
 298         } // end for ie
 299         // if delta ray only
 300         if(ntr == 0) ntr = 1;
 301         clusterI->SetNTracks(ntr);
 302         clusterI->SetTracks(tr0,tr1,tr2);
 303     } // I cluster
 304 }
 305 //______________________________________________________________________
 306 void AliITSClusterFinderSPDdubna::GetRecPoints(){
 307     // get rec points
 308     static AliITS *iTS = (AliITS*)gAlice->GetModule("ITS");
 309     // get number of clusters for this module
 310     Int_t nofClusters   = fClusters->GetEntriesFast();
 311     const Float_t kconv = 1.0e-4;
 312     const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
 313     const Float_t kRMSz = 70.0*kconv; // microns -> cm ITS TDR Table 1.3
 314     Float_t spdLength   = fSegmentation->Dz();
 315     Float_t spdWidth    = fSegmentation->Dx();
 316     Int_t   dummy       = 0;
 317     Float_t xpitch      = fSegmentation->Dpx(dummy);
 318     Int_t i;
 319     Int_t track0, track1, track2;
 320
 321     nofClusters -= fNclusters;
 322     for(i=0; i<nofClusters; i++) {
 323         AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
 324         Int_t clustersizex = clusterI->NclX();
 325         Int_t clustersizez = clusterI->NclZ();
 326         Int_t xstart = clusterI->XStartf();
 327         Float_t clusterx = 0;
 328         Float_t clusterz = clusterI->Z();
 329         for(Int_t ii=0; ii<clustersizex; ii++) clusterx += (xstart+0.5+ii)*xpitch;
 330         clusterx /= clustersizex;
 331         clusterz /= clustersizez;
 332         clusterI->GetTracks(track0, track1, track2);
 333         AliITSRecPoint rnew;
 334         rnew.SetX((clusterx - spdWidth/2)*kconv);
 335         rnew.SetZ((clusterz - spdLength/2)*kconv);
 336         rnew.SetQ(1.);
 337         rnew.SetdEdX(0.);
 338         rnew.SetSigmaX2(kRMSx*kRMSx);
 339         rnew.SetSigmaZ2(kRMSz*kRMSz);
 340         rnew.fTracks[0]=track0;
 341         rnew.fTracks[1]=track1;
 342         rnew.fTracks[2]=track2;
 343         iTS->AddRecPoint(rnew);
 344     } // I clusters
 345     fMap->ClearMap();
 346 }
 347 //______________________________________________________________________
 348 void AliITSClusterFinderSPDdubna::FindRawClusters(Int_t mod){
 349     // find raw clusters
 350
 351 #ifndef DEBUG
 352     cout<<"SPDdubna: module ="<<mod<<endl;
 353 #endif
 354     Find1DClusters(mod);
 355     GroupClusters();
 356     TracksInCluster();
 357     GetRecPoints();
 358 }