Variable name fixed.
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSDD.cxx
CommitLineData
b0f5e3fc 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
22ClassImp(AliITSClusterFinderSDD)
23
24//----------------------------------------------------------
25AliITSClusterFinderSDD::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//_____________________________________________________________________________
44AliITSClusterFinderSDD::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//__________________________________________________________________________
62AliITSClusterFinderSDD::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//_________________________________________________________________________
77AliITSClusterFinderSDD&
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//_____________________________________________________________________________
93void AliITSClusterFinderSDD::SetMap()
94{
95 // set map
96 if(!fMap) fMap=new AliITSMapA2(fSegmentation);
97
98}
99//_____________________________________________________________________________
100void 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
124void 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//_____________________________________________________________________________
293void 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
338void 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
366void 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
406void AliITSClusterFinderSDD::FindRawClusters()
407{
408 // find raw clusters
409 Find1DClusters();
410 GroupClusters();
411 SelectClusters();
412 GetRecPoints();
413}
414
415
416
417