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