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