ITS/AliITSClusterFinderSDD.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 #include "AliITSClusterFinderSDD.h"
  18 #include "AliITSMapA2.h"
  19 #include "AliITSMapA1.h"
  20 #include "AliITS.h"
  21 #include "AliRun.h"
  22
  23
  24
  25 ClassImp(AliITSClusterFinderSDD)
  26
  27 //----------------------------------------------------------
  28 AliITSClusterFinderSDD::AliITSClusterFinderSDD
  29 (AliITSsegmentation *seg, AliITSresponse *response, TClonesArray *digits, TClonesArray *recp)
  30 {
  31   // constructor
  32     fSegmentation=seg;
  33     fResponse=response;
  34     fDigits=digits;
  35     fClusters=recp;
  36     fNclusters= fClusters->GetEntriesFast();
  37     SetCutAmplitude();
  38     SetDAnode();
  39     SetDTime();
  40     SetMinPeak();
  41     SetNCells();
  42     SetMap();
  43 }
  44
  45 //_____________________________________________________________________________
  46 AliITSClusterFinderSDD::AliITSClusterFinderSDD()
  47 {
  48   // constructor
  49     fSegmentation=0;
  50     fResponse=0;
  51     fDigits=0;
  52     fClusters=0;
  53     fNclusters=0;
  54     fMap=0;
  55     fMapA2=0;
  56     SetCutAmplitude();
  57     SetDAnode();
  58     SetDTime();
  59     SetMinPeak();
  60     SetNCells();
  61
  62 }
  63
  64 //_____________________________________________________________________________
  65 AliITSClusterFinderSDD::~AliITSClusterFinderSDD()
  66 {
  67     // destructor
  68
  69     if(fMap) delete fMap;
  70     if(fMapA2) delete fMapA2;
  71
  72 }
  73 //__________________________________________________________________________
  74 AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &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->fMapA2 = source.fMapA2 ;
  81   this->fCutAmplitude = source.fCutAmplitude ;
  82   this->fDAnode = source.fDAnode ;
  83   this->fDTime = source.fDTime ;
  84   this->fMinPeak = source.fMinPeak ;
  85   this->fMinNCells = source.fMinNCells ;
  86   return;
  87 }
  88
  89 //_________________________________________________________________________
  90 AliITSClusterFinderSDD&
  91   AliITSClusterFinderSDD::operator=(const AliITSClusterFinderSDD &source) {
  92   //    Assignment operator
  93   if(&source == this) return *this;
  94   this->fClusters = source.fClusters ;
  95   this->fNclusters = source.fNclusters ;
  96   this->fMap = source.fMap ;
  97   this->fMapA2 = source.fMapA2 ;
  98   this->fCutAmplitude = source.fCutAmplitude ;
  99   this->fDAnode = source.fDAnode ;
 100   this->fDTime = source.fDTime ;
 101   this->fMinPeak = source.fMinPeak ;
 102   this->fMinNCells = source.fMinNCells ;
 103   return *this;
 104 }
 105
 106 //_____________________________________________________________________________
 107 void AliITSClusterFinderSDD::SetMap()
 108 {
 109   // set map
 110     if(!fMapA2) fMapA2=new AliITSMapA2(fSegmentation,fDigits,(double)fCutAmplitude);
 111     if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);
 112
 113 }
 114 //_____________________________________________________________________________
 115
 116 void AliITSClusterFinderSDD::Find1DClusters()
 117 {
 118   // find 1D clusters
 119
 120     AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 121
 122    // retrieve the parameters
 123    Int_t fNofMaps = fSegmentation->Npz();
 124    Int_t fMaxNofSamples = fSegmentation->Npx();
 125    Int_t fNofAnodes = fNofMaps/2;
 126    Int_t dummy=0;
 127    Float_t fTimeStep = fSegmentation->Dpx(dummy);
 128    Float_t fSddLength = fSegmentation->Dx();
 129    Float_t fDriftSpeed = fResponse->DriftSpeed();
 130
 131    Float_t anodePitch = fSegmentation->Dpz(dummy);
 132    // map the signal
 133    fMapA2->FillMap();
 134
 135    // Piergiorgio's code - do not subtract baseline since we start
 136    // from digits and do not duplicate arrays, i.e. use fMap instead
 137    // of Float_t fadc[2*fNofAnodes][fMaxNofSamples];
 138
 139
 140   Int_t nofFoundClusters = 0;
 141   Int_t i;
 142   Float_t **dfadc = new Float_t*[fNofMaps];
 143   for(i=0;i<fNofMaps;i++) dfadc[i] = new Float_t[fMaxNofSamples];
 144   Float_t fadc, fadc1, fadc2;
 145   Int_t j,k,idx,l,m;
 146   for(j=0;j<2;j++) {
 147     for(k=0;k<fNofAnodes;k++) {
 148       idx = j*fNofAnodes+k;
 149       // signal (fadc) & derivative (dfadc)
 150       for(l=0; l<fMaxNofSamples; l++) {
 151         fadc2=(Float_t)fMapA2->GetSignal(idx,l);
 152         fadc1=(Float_t)fMapA2->GetSignal(idx,l-1);
 153         if(l>0) dfadc[k][l-1] = fadc2-fadc1;
 154       } // samples
 155     } // anodes
 156
 157     for(k=0;k<fNofAnodes;k++) {
 158       idx = j*fNofAnodes+k;
 159
 160       Int_t imax = 0;
 161       Int_t imaxd = 0;
 162       Int_t it=0;
 163       while(it <= fMaxNofSamples-3) {
 164         //      cout << "sample: " << it << endl;
 165
 166         imax = it;
 167         imaxd = it;
 168         // maximum of signal
 169
 170         Float_t fadcmax = 0.;
 171         Float_t dfadcmax = 0.;
 172         Int_t lthrmina = 1;
 173         //      if(it >= 60) lthrmina = 2;
 174         //      if(it >= 100) lthrmina = 3;
 175         Int_t lthrmint = 2;
 176         //if(it >= 60) lthrmint = 3;
 177         //if(it >= 100) lthrmint = 4;
 178
 179         Int_t lthra = 0;
 180         Int_t lthrt = 0;
 181         for(m=0;m<10;m++) {
 182           Int_t id = it+m;
 183           if(id>=fMaxNofSamples) break;
 184           fadc=fMapA2->GetSignal(idx,id);
 185           if(fadc > fadcmax) {
 186             fadcmax = fadc;
 187             if(fadc > 0) { lthra++; lthrt++; }
 188             imax = id;
 189           }
 190           if(dfadc[k][id] > dfadcmax) {
 191             dfadcmax = dfadc[k][id];
 192             imaxd = id;
 193           }
 194         }
 195         it = imaxd;
 196         // skip if no signal over threshold
 197         if(fMapA2->TestHit(idx,imax) == kEmpty) {it++; continue;}
 198
 199         if(k>0) {
 200           if(fMapA2->TestHit(idx-1,imax) != kEmpty) lthra++;
 201         }
 202         if(k<fNofAnodes-1)
 203           if(fMapA2->TestHit(idx+1,imax) != kEmpty) lthra++;
 204
 205         if(imax>0) {
 206           if(fMapA2->TestHit(idx,imax-1) != kEmpty) lthrt++;
 207         }
 208         if(imax<fMaxNofSamples)
 209           if(fMapA2->TestHit(idx,imax+1) != kEmpty) lthrt++;
 210
 211         // cluster charge
 212         Int_t tstart = it-1;
 213
 214         Bool_t ilcl = 0;
 215         if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
 216
 217         if(ilcl) {
 218           nofFoundClusters++;
 219           Int_t tstop = tstart;
 220           Float_t dfadcmin = 10000.;
 221           Int_t ij;
 222           for(ij=0; ij<10; ij++) {
 223             if(dfadc[k][it+ij] < dfadcmin) {
 224               tstop = it+ij+1;
 225               dfadcmin = dfadc[k][it+ij];
 226             }
 227           }
 228
 229           Float_t clusterCharge = 0.;
 230           Float_t clusterAnode = k+0.5;
 231           Float_t clusterTime = 0.;
 232           Float_t clusterMult = 0.;
 233           Float_t clusterPeakAmplitude = 0.;
 234           Int_t its;
 235           for(its=tstart; its<=tstop; its++) {
 236             fadc=fMapA2->GetSignal(idx,its);
 237             clusterCharge += fadc;
 238             if(fadc > clusterPeakAmplitude) clusterPeakAmplitude = fadc;
 239             clusterTime += fadc*its;
 240             clusterMult++;
 241             if(its == tstop) {
 242                clusterTime /= (clusterCharge/fTimeStep);   // ns
 243                clusterCharge *= (fTimeStep/160.);          // keV
 244                if(clusterTime > 58.2) clusterTime -= 58.2;   // ns
 245               /*
 246               else {
 247                 cout << "Warning: cluster with negative time " << clusterTime << ", peak ampl.: " << clusterPeakAmplitude << ", mult: " << clusterMult << ", charge: " << clusterCharge << endl;
 248                 cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
 249               }
 250               */
 251               }
 252           }
 253           // cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
 254
 255           Float_t clusteranodePath = (0.06 + clusterAnode - fNofAnodes/2)*anodePitch;
 256           Float_t clusterDriftPath = clusterTime*fDriftSpeed;
 257           clusterDriftPath = fSddLength-clusterDriftPath;
 258
 259           if(clusterCharge < 0.) break;
 260
 261           //printf("wing clusterMult clusterAnode clusterTime %d %f %f %f \n",j+1,clusterMult, clusterAnode, clusterTime);
 262
 263           AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,0.,0.,clusterDriftPath,clusteranodePath,clusterMult);
 264           //fClusters->Add(point);
 265           iTS->AddCluster(1,clust);
 266           it = tstop;
 267         } // ilcl
 268
 269         it++;
 270
 271       } // while (samples)
 272     } // anodes
 273   } // detectors (2)
 274
 275
 276   fMapA2->ClearMap();
 277
 278   for(i=0;i<fNofMaps;i++) delete[] dfadc[i];
 279   delete [] dfadc;
 280
 281   return;
 282
 283 }
 284
 285 //_____________________________________________________________________________
 286 void  AliITSClusterFinderSDD::GroupClusters()
 287 {
 288   // group clusters
 289   Int_t dummy=0;
 290   Float_t fTimeStep = fSegmentation->Dpx(dummy);
 291
 292
 293   // get number of clusters for this module
 294   Int_t nofClusters = fClusters->GetEntriesFast();
 295   nofClusters -= fNclusters;
 296
 297   AliITSRawClusterSDD *clusterI;
 298   AliITSRawClusterSDD *clusterJ;
 299
 300   Int_t *label = new Int_t [nofClusters];
 301   Int_t i,j;
 302   for(i=0; i<nofClusters; i++) label[i] = 0;
 303   for(i=0; i<nofClusters; i++) {
 304     if(label[i] != 0) continue;
 305     for(j=i+1; j<nofClusters; j++) {
 306       if(label[j] != 0) continue;
 307       clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
 308       clusterJ = (AliITSRawClusterSDD*) fClusters->At(j);
 309       // 1.3 good
 310       if(clusterI->T() < fTimeStep*60) fDAnode = 3.2;
 311       if(clusterI->T() < fTimeStep*10) fDAnode = 1.2;
 312       Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
 313       if(!pair) continue;
 314       //      clusterI->PrintInfo();
 315       //      clusterJ->PrintInfo();
 316       clusterI->Add(clusterJ);
 317       label[j] = 1;
 318       fClusters->RemoveAt(j);
 319     } // J clusters
 320     label[i] = 1;
 321   } // I clusters
 322   fClusters->Compress();
 323
 324   delete [] label;
 325   return;
 326
 327 }
 328
 329 //_____________________________________________________________________________
 330
 331 void AliITSClusterFinderSDD::SelectClusters()
 332 {
 333   // get number of clusters for this module
 334   Int_t nofClusters = fClusters->GetEntriesFast();
 335   nofClusters -= fNclusters;
 336
 337   Int_t i;
 338   for(i=0; i<nofClusters; i++) {
 339     AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
 340     Int_t rmflg = 0;
 341     Float_t wy = 0.;
 342     if(clusterI->Anodes() != 0.) {
 343       wy = ((Float_t) clusterI->Samples())/clusterI->Anodes();
 344     }
 345     Float_t amp = clusterI->PeakAmpl();
 346     if(amp < fMinPeak) rmflg = 1;
 347     if(wy < fMinNCells) rmflg = 1;
 348     if(rmflg) fClusters->RemoveAt(i);
 349   } // I clusters
 350   fClusters->Compress();
 351   return;
 352
 353 }
 354
 355 //_____________________________________________________________________________
 356
 357 void AliITSClusterFinderSDD::GetRecPoints()
 358 {
 359   // get rec points
 360
 361   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 362
 363   // get number of clusters for this module
 364   Int_t nofClusters = fClusters->GetEntriesFast();
 365   nofClusters -= fNclusters;
 366
 367   const Float_t kconvGeV = 1.e-6; // GeV -> KeV
 368   const Float_t kconv = 1.0e-4;
 369   const Float_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
 370   const Float_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3
 371
 372   fMap->FillMap();
 373
 374   Int_t i, ix, iz;
 375   for(i=0; i<nofClusters; i++) {
 376     AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
 377     fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
 378     AliITSdigitSDD *dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix-1);
 379     AliITSRecPoint rnew;
 380     rnew.SetX(clusterI->X());
 381     rnew.SetZ(clusterI->Z());
 382     rnew.SetQ(clusterI->Q());   // in KeV - should be ADC
 383     rnew.SetdEdX(kconvGeV*clusterI->Q());
 384     rnew.SetSigmaX2(kRMSx*kRMSx);
 385     rnew.SetSigmaZ2(kRMSz*kRMSz);
 386     rnew.fTracks[0]=dig->fTracks[0];
 387     rnew.fTracks[1]=dig->fTracks[1];
 388     rnew.fTracks[2]=dig->fTracks[2];
 389     iTS->AddRecPoint(rnew);
 390   } // I clusters
 391
 392   fMap->ClearMap();
 393 }
 394
 395 //_____________________________________________________________________________
 396
 397 void AliITSClusterFinderSDD::FindRawClusters()
 398 {
 399   // find raw clusters
 400     Find1DClusters();
 401     GroupClusters();
 402     SelectClusters();
 403     GetRecPoints();
 404 }
 405
 406
 407
 408