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           Float_t n, baseline;
 236           fResponse->GetNoiseParam(n,baseline);
 237           for(its=tstart; its<=tstop; its++) {
 238             fadc=fMapA2->GetSignal(idx,its)-baseline;
 239             clusterCharge += fadc;
 240             if(fadc > clusterPeakAmplitude) clusterPeakAmplitude = fadc;
 241             clusterTime += fadc*its;
 242             clusterMult++;
 243             if(its == tstop) {
 244                clusterTime /= (clusterCharge/fTimeStep);   // ns
 245                clusterCharge *= (fTimeStep/160.);          // keV
 246                if(clusterTime > 58.2) clusterTime -= 58.2;   // ns
 247               /*
 248               else {
 249                 cout << "Warning: cluster with negative time " << clusterTime << ", peak ampl.: " << clusterPeakAmplitude << ", mult: " << clusterMult << ", charge: " << clusterCharge << endl;
 250                 cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
 251               }
 252               */
 253               }
 254           }
 255           // cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
 256
 257           Float_t clusteranodePath = (0.06 + clusterAnode - fNofAnodes/2)*anodePitch;
 258           Float_t clusterDriftPath = clusterTime*fDriftSpeed;
 259           clusterDriftPath = fSddLength-clusterDriftPath;
 260
 261           if(clusterCharge < 0.) break;
 262
 263           //printf("wing clusterMult clusterAnode clusterTime %d %f %f %f \n",j+1,clusterMult, clusterAnode, clusterTime);
 264
 265           AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,0.,0.,clusterDriftPath,clusteranodePath,clusterMult);
 266           //fClusters->Add(point);
 267           iTS->AddCluster(1,clust);
 268           it = tstop;
 269         } // ilcl
 270
 271         it++;
 272
 273       } // while (samples)
 274     } // anodes
 275   } // detectors (2)
 276
 277
 278   fMapA2->ClearMap();
 279
 280   for(i=0;i<fNofMaps;i++) delete[] dfadc[i];
 281   delete [] dfadc;
 282
 283   return;
 284
 285 }
 286
 287 //_____________________________________________________________________________
 288 void  AliITSClusterFinderSDD::GroupClusters()
 289 {
 290   // group clusters
 291   Int_t dummy=0;
 292   Float_t fTimeStep = fSegmentation->Dpx(dummy);
 293
 294
 295   // get number of clusters for this module
 296   Int_t nofClusters = fClusters->GetEntriesFast();
 297   nofClusters -= fNclusters;
 298
 299   AliITSRawClusterSDD *clusterI;
 300   AliITSRawClusterSDD *clusterJ;
 301
 302   Int_t *label = new Int_t [nofClusters];
 303   Int_t i,j;
 304   for(i=0; i<nofClusters; i++) label[i] = 0;
 305   for(i=0; i<nofClusters; i++) {
 306     if(label[i] != 0) continue;
 307     for(j=i+1; j<nofClusters; j++) {
 308       if(label[j] != 0) continue;
 309       clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
 310       clusterJ = (AliITSRawClusterSDD*) fClusters->At(j);
 311       // 1.3 good
 312       if(clusterI->T() < fTimeStep*60) fDAnode = 3.2;
 313       if(clusterI->T() < fTimeStep*10) fDAnode = 1.2;
 314       Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
 315       if(!pair) continue;
 316       //      clusterI->PrintInfo();
 317       //      clusterJ->PrintInfo();
 318       clusterI->Add(clusterJ);
 319       label[j] = 1;
 320       fClusters->RemoveAt(j);
 321     } // J clusters
 322     label[i] = 1;
 323   } // I clusters
 324   fClusters->Compress();
 325
 326   delete [] label;
 327   return;
 328
 329 }
 330
 331 //_____________________________________________________________________________
 332
 333 void AliITSClusterFinderSDD::SelectClusters()
 334 {
 335   // get number of clusters for this module
 336   Int_t nofClusters = fClusters->GetEntriesFast();
 337   nofClusters -= fNclusters;
 338
 339   Int_t i;
 340   for(i=0; i<nofClusters; i++) {
 341     AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
 342     Int_t rmflg = 0;
 343     Float_t wy = 0.;
 344     if(clusterI->Anodes() != 0.) {
 345       wy = ((Float_t) clusterI->Samples())/clusterI->Anodes();
 346     }
 347     Float_t amp = clusterI->PeakAmpl();
 348     if(amp < fMinPeak) rmflg = 1;
 349     if(wy < fMinNCells) rmflg = 1;
 350     if(rmflg) fClusters->RemoveAt(i);
 351   } // I clusters
 352   fClusters->Compress();
 353   return;
 354
 355 }
 356
 357 //_____________________________________________________________________________
 358
 359 void AliITSClusterFinderSDD::GetRecPoints()
 360 {
 361   // get rec points
 362
 363   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
 364
 365   // get number of clusters for this module
 366   Int_t nofClusters = fClusters->GetEntriesFast();
 367   nofClusters -= fNclusters;
 368
 369   const Float_t kconvGeV = 1.e-6; // GeV -> KeV
 370   const Float_t kconv = 1.0e-4;
 371   const Float_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
 372   const Float_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3
 373
 374   fMap->FillMap();
 375
 376   Int_t i, ix, iz;
 377   for(i=0; i<nofClusters; i++) {
 378     AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
 379     fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
 380     AliITSdigitSDD *dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix-1);
 381     AliITSRecPoint rnew;
 382     rnew.SetX(clusterI->X());
 383     rnew.SetZ(clusterI->Z());
 384     rnew.SetQ(clusterI->Q());   // in KeV - should be ADC
 385     rnew.SetdEdX(kconvGeV*clusterI->Q());
 386     rnew.SetSigmaX2(kRMSx*kRMSx);
 387     rnew.SetSigmaZ2(kRMSz*kRMSz);
 388     rnew.fTracks[0]=dig->fTracks[0];
 389     rnew.fTracks[1]=dig->fTracks[1];
 390     rnew.fTracks[2]=dig->fTracks[2];
 391     iTS->AddRecPoint(rnew);
 392   } // I clusters
 393
 394   fMap->ClearMap();
 395 }
 396
 397 //_____________________________________________________________________________
 398
 399 void AliITSClusterFinderSDD::FindRawClusters()
 400 {
 401   // find raw clusters
 402     Find1DClusters();
 403     GroupClusters();
 404     SelectClusters();
 405     GetRecPoints();
 406 }
 407
 408
 409
 410