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