1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 #include "AliITSClusterFinderSDD.h"
19 #include "AliITSMapA1.h"
21 #include "AliITSdigit.h"
22 #include "AliITSRawCluster.h"
23 #include "AliITSRecPoint.h"
24 #include "AliITSsegmentation.h"
25 #include "AliITSresponse.h"
30 ClassImp(AliITSClusterFinderSDD)
32 //----------------------------------------------------------
33 AliITSClusterFinderSDD::AliITSClusterFinderSDD
34 (AliITSsegmentation *seg, AliITSresponse *response, TClonesArray *digits, TClonesArray *recp)
42 fNclusters= fClusters->GetEntriesFast();
50 fMap=new AliITSMapA1(fSegmentation,fDigits,fCutAmplitude);
54 //_____________________________________________________________________________
55 AliITSClusterFinderSDD::AliITSClusterFinderSDD()
74 //_____________________________________________________________________________
75 AliITSClusterFinderSDD::~AliITSClusterFinderSDD()
82 //__________________________________________________________________________
83 AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &source){
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 ;
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 ;
118 //_____________________________________________________________________________
120 void AliITSClusterFinderSDD::Find1DClusters()
124 AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
126 // retrieve the parameters
127 Int_t fNofMaps = fSegmentation->Npz();
128 Int_t fMaxNofSamples = fSegmentation->Npx();
129 Int_t fNofAnodes = fNofMaps/2;
131 Float_t fTimeStep = fSegmentation->Dpx(dummy);
132 Float_t fSddLength = fSegmentation->Dx();
133 Float_t fDriftSpeed = fResponse->DriftSpeed();
135 Float_t anodePitch = fSegmentation->Dpz(dummy);
137 fMap->SetThreshold(fCutAmplitude);
142 fResponse->GetNoiseParam(noise,baseline);
144 Float_t maxadc = fResponse->MaxAdc();
145 Float_t topValue = fResponse->MagicValue();
146 Float_t norm = maxadc/topValue;
148 Int_t nofFoundClusters = 0;
150 Float_t **dfadc = new Float_t*[fNofAnodes];
151 for(i=0;i<fNofAnodes;i++) dfadc[i] = new Float_t[fMaxNofSamples];
152 Float_t fadc, fadc1, fadc2;
155 for(k=0;k<fNofAnodes;k++) {
156 idx = j*fNofAnodes+k;
157 // signal (fadc) & derivative (dfadc)
158 for(l=0; l<fMaxNofSamples; l++) {
159 fadc2=(Float_t)fMap->GetSignal(idx,l);
160 if(l>0) fadc1=(Float_t)fMap->GetSignal(idx,l-1);
161 if(l>0) dfadc[k][l-1] = fadc2-fadc1;
165 for(k=0;k<fNofAnodes;k++) {
166 //cout << "Anode: " << k+1 << ", Wing: " << j+1 << endl;
167 idx = j*fNofAnodes+k;
172 while(it <= fMaxNofSamples-3) {
178 Float_t fadcmax = 0.;
179 Float_t dfadcmax = 0.;
181 // if(it >= 60) lthrmina = 2;
182 // if(it >= 100) lthrmina = 3;
184 //if(it >= 60) lthrmint = 3;
185 //if(it >= 100) lthrmint = 4;
192 if(id>=fMaxNofSamples) break;
193 fadc=(float)fMap->GetSignal(idx,id);
194 if(fadc > fadcmax) { fadcmax = fadc; imax = id;}
195 if(fadc > (float)fCutAmplitude) {
199 if(dfadc[k][id] > dfadcmax) {
200 dfadcmax = dfadc[k][id];
206 if(fMap->TestHit(idx,imax) == kEmpty) {it++; continue;}
212 if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
213 //printf("ilcl %d\n",ilcl);
216 fResponse->GetNoiseParam(n,b);
220 Int_t tstop = tstart;
221 Float_t dfadcmin = 10000.;
223 for(ij=0; ij<20; ij++) {
224 if(dfadc[k][it+ij] < dfadcmin) {
226 dfadcmin = dfadc[k][it+ij];
230 Float_t clusterCharge = 0.;
231 Float_t clusterAnode = k+0.5;
232 Float_t clusterTime = 0.;
233 Float_t clusterMult = 0.;
234 Float_t clusterPeakAmplitude = 0.;
235 Int_t its,peakpos=-1;
237 fResponse->GetNoiseParam(n,baseline);
238 for(its=tstart; its<=tstop; its++) {
239 fadc=(float)fMap->GetSignal(idx,its);
244 clusterCharge += fadc;
245 // as a matter of fact we should take the peak pos before FFT
246 // to get the list of tracks !!!
247 if(fadc > clusterPeakAmplitude) {
248 clusterPeakAmplitude = fadc;
249 //peakpos=fMap->GetHitIndex(idx,its);
250 Int_t shift=(int)(fTimeCorr/fTimeStep);
251 if(its>shift && its<(fMaxNofSamples-shift)) peakpos=fMap->GetHitIndex(idx,its+shift);
252 else peakpos=fMap->GetHitIndex(idx,its);
253 if(peakpos<0) peakpos=fMap->GetHitIndex(idx,its);
255 clusterTime += fadc*its;
258 // charge from ADC back to nA
259 //clusterCharge /= norm;
260 if(clusterCharge <= 0.) printf("clusterCharge %f norm %f\n",clusterCharge,norm);
261 clusterTime /= (clusterCharge/fTimeStep); // ns
262 clusterCharge *= (fTimeStep/160.); // keV
263 if(clusterTime > fTimeCorr) clusterTime -= fTimeCorr; // ns
266 // cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
268 Float_t clusteranodePath = (clusterAnode - fNofAnodes/2)*anodePitch;
269 Float_t clusterDriftPath = clusterTime*fDriftSpeed;
270 clusterDriftPath = fSddLength-clusterDriftPath;
272 if(clusterCharge <= 0.) break;
274 AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,peakpos,0.,0.,clusterDriftPath,clusteranodePath,clusterMult);
275 iTS->AddCluster(1,clust);
288 for(i=0;i<fNofAnodes;i++) delete[] dfadc[i];
295 //_____________________________________________________________________________
296 void AliITSClusterFinderSDD::GroupClusters()
300 Float_t fTimeStep = fSegmentation->Dpx(dummy);
303 // get number of clusters for this module
304 Int_t nofClusters = fClusters->GetEntriesFast();
305 nofClusters -= fNclusters;
307 AliITSRawClusterSDD *clusterI;
308 AliITSRawClusterSDD *clusterJ;
310 Int_t *label = new Int_t [nofClusters];
312 for(i=0; i<nofClusters; i++) label[i] = 0;
313 for(i=0; i<nofClusters; i++) {
314 if(label[i] != 0) continue;
315 for(j=i+1; j<nofClusters; j++) {
316 if(label[j] != 0) continue;
317 clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
318 clusterJ = (AliITSRawClusterSDD*) fClusters->At(j);
320 if(clusterI->T() < fTimeStep*60) fDAnode = 3.2;
321 if(clusterI->T() < fTimeStep*10) fDAnode = 1.2;
322 Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
324 // clusterI->PrintInfo();
325 // clusterJ->PrintInfo();
326 clusterI->Add(clusterJ);
328 fClusters->RemoveAt(j);
332 fClusters->Compress();
339 //_____________________________________________________________________________
341 void AliITSClusterFinderSDD::SelectClusters()
343 // get number of clusters for this module
344 Int_t nofClusters = fClusters->GetEntriesFast();
345 nofClusters -= fNclusters;
348 for(i=0; i<nofClusters; i++) {
349 AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
352 if(clusterI->Anodes() != 0.) {
353 wy = ((Float_t) clusterI->Samples())/clusterI->Anodes();
355 Float_t amp = clusterI->PeakAmpl();
356 if(amp < fMinPeak) rmflg = 1;
357 if(wy < fMinNCells) rmflg = 1;
358 //if(wy > fMaxNCells) rmflg = 1;
359 if(rmflg) fClusters->RemoveAt(i);
361 fClusters->Compress();
366 //_____________________________________________________________________________
368 void AliITSClusterFinderSDD::GetRecPoints()
372 AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
374 // get number of clusters for this module
375 Int_t nofClusters = fClusters->GetEntriesFast();
376 nofClusters -= fNclusters;
378 const Float_t kconvGeV = 1.e-6; // GeV -> KeV
379 const Float_t kconv = 1.0e-4;
380 const Float_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
381 const Float_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3
385 Int_t ix, iz, idx=-1;
386 AliITSdigitSDD *dig=0;
387 // Int_t maxt=fSegmentation->Npx();
388 Int_t ndigits=fDigits->GetEntriesFast();
389 for(i=0; i<nofClusters; i++) {
390 AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*)fClusters->At(i);
391 if(!clusterI) Error("SDD: GetRecPoints","i clusterI ",i,clusterI);
392 if(clusterI) idx=clusterI->PeakPos();
393 if(idx>ndigits) Error("SDD: GetRecPoints","idx ndigits",idx,ndigits);
394 // try peak neighbours - to be done
395 if(idx && idx <= ndigits) dig = (AliITSdigitSDD*)fDigits->UncheckedAt(idx);
398 fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
399 dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix-1);
400 // if null try neighbours
401 if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix);
402 if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix+1);
403 if (!dig) printf("SDD: cannot assign the track number!\n");
407 rnew.SetX(clusterI->X());
408 rnew.SetZ(clusterI->Z());
409 rnew.SetQ(clusterI->Q()); // in KeV - should be ADC
410 rnew.SetdEdX(kconvGeV*clusterI->Q());
411 rnew.SetSigmaX2(kRMSx*kRMSx);
412 rnew.SetSigmaZ2(kRMSz*kRMSz);
413 if(dig) rnew.fTracks[0]=dig->fTracks[0];
414 if(dig) rnew.fTracks[1]=dig->fTracks[1];
415 if(dig) rnew.fTracks[2]=dig->fTracks[2];
416 //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());
417 iTS->AddRecPoint(rnew);
423 //_____________________________________________________________________________
425 void AliITSClusterFinderSDD::FindRawClusters()