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