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