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