Correction of the inheritance scheme
[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
44b3710f 131 Int_t i;
b0f5e3fc 132 Int_t fNofMaps = fSegmentation->Npz();
133 Int_t fMaxNofSamples = fSegmentation->Npx();
134 Int_t fNofAnodes = fNofMaps/2;
135 Int_t dummy=0;
136 Float_t fTimeStep = fSegmentation->Dpx(dummy);
137 Float_t fSddLength = fSegmentation->Dx();
138 Float_t fDriftSpeed = fResponse->DriftSpeed();
139
140 Float_t anodePitch = fSegmentation->Dpz(dummy);
141 // map the signal
142 FillMap();
143
144 // Piergiorgio's code - do not subtract baseline since we start
145 // from digits and do not duplicate arrays, i.e. use fMap instead
146 // of Float_t fadc[2*fNofAnodes][fMaxNofSamples];
147
148
149 Int_t nofFoundClusters = 0;
44b3710f 150
151 Float_t **dfadc = new Float_t*[fNofMaps];
152 for(i=0;i<fNofMaps;i++) dfadc[i] = new Float_t[fMaxNofSamples];
153
154
b0f5e3fc 155 Float_t fadc, fadc1, fadc2;
156 Int_t j,k,idx,l,m;
157 for(j=0;j<2;j++) {
158 for(k=0;k<fNofAnodes;k++) {
159 idx = j*fNofAnodes+k;
160 // signal (fadc) & derivative (dfadc)
161 for(l=0; l<fMaxNofSamples; l++) {
162 fadc2=(Float_t)fMap->GetSignal(idx,l);
163 fadc1=(Float_t)fMap->GetSignal(idx,l-1);
164 if(l>0) dfadc[k][l-1] = fadc2-fadc1;
165 } // samples
166 } // anodes
167
168 for(k=0;k<fNofAnodes;k++) {
169 idx = j*fNofAnodes+k;
170
171 Int_t imax = 0;
172 Int_t imaxd = 0;
173 Int_t it=0;
174 while(it <= fMaxNofSamples-3) {
175 // cout << "sample: " << it << endl;
176
177 imax = it;
178 imaxd = it;
179 // maximum of signal
180
181 Float_t fadcmax = 0.;
182 Float_t dfadcmax = 0.;
183 Int_t lthrmina = 1;
184 // if(it >= 60) lthrmina = 2;
185 // if(it >= 100) lthrmina = 3;
186 Int_t lthrmint = 2;
187 //if(it >= 60) lthrmint = 3;
188 //if(it >= 100) lthrmint = 4;
189
190 Int_t lthra = 0;
191 Int_t lthrt = 0;
192 for(m=0;m<10;m++) {
193 Int_t id = it+m;
194 if(id>=fMaxNofSamples) break;
195 fadc=fMap->GetSignal(idx,id);
196 if(fadc > fadcmax) {
197 fadcmax = fadc;
198 if(fadc > fCutAmplitude) { lthra++; lthrt++; }
199 imax = id;
200 }
201 if(dfadc[k][id] > dfadcmax) {
202 dfadcmax = dfadc[k][id];
203 imaxd = id;
204 }
205 }
206 it = imaxd;
207 // skip if no signal over threshold
208 if(fMap->GetSignal(idx,imax) < fCutAmplitude) {it++; continue;}
209
210 if(k>0) {
211 if(fMap->GetSignal(idx-1,imax) > fCutAmplitude) lthra++;
212 }
213 if(k<fNofAnodes-1)
214 if(fMap->GetSignal(idx+1,imax) > fCutAmplitude) lthra++;
215
216 if(imax>0) {
217 if(fMap->GetSignal(idx,imax-1) > fCutAmplitude) lthrt++;
218 }
219 if(imax<fMaxNofSamples)
220 if(fMap->GetSignal(idx,imax+1) > fCutAmplitude) lthrt++;
221
222 // cluster charge
223 Int_t tstart = it-1;
224
225 Bool_t ilcl = 0;
226 if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
227
228 if(ilcl) {
229 nofFoundClusters++;
230 Int_t tstop = tstart;
231 Float_t dfadcmin = 10000.;
232 Int_t ij;
233 for(ij=0; ij<10; ij++) {
234 if(dfadc[k][it+ij] < dfadcmin) {
235 tstop = it+ij+1;
236 dfadcmin = dfadc[k][it+ij];
237 }
238 }
239
240 Float_t clusterCharge = 0.;
241 Float_t clusterAnode = k+0.5;
242 Float_t clusterTime = 0.;
243 Float_t clusterMult = 0.;
244 Float_t clusterPeakAmplitude = 0.;
245 Int_t its;
246 for(its=tstart; its<=tstop; its++) {
247 fadc=fMap->GetSignal(idx,its);
248 clusterCharge += fadc;
249 if(fadc > clusterPeakAmplitude) clusterPeakAmplitude = fadc;
250 clusterTime += fadc*its;
251 clusterMult++;
252 if(its == tstop) {
253 clusterTime /= (clusterCharge/fTimeStep); // ns
254 clusterCharge *= (fTimeStep/160.); // keV
255 if(clusterTime > 58.2) clusterTime -= 58.2; // ns
256 /*
257 else {
258 cout << "Warning: cluster with negative time " << clusterTime << ", peak ampl.: " << clusterPeakAmplitude << ", mult: " << clusterMult << ", charge: " << clusterCharge << endl;
259 cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
260 }
261 */
262 }
263 }
264 // cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
265
266 Float_t clusteranodePath = (0.06 + clusterAnode - fNofAnodes/2)*anodePitch;
267 Float_t clusterDriftPath = clusterTime*fDriftSpeed;
268 clusterDriftPath = fSddLength-clusterDriftPath;
269
270 if(clusterCharge < 0.) break;
271
272 //printf("wing clusterMult clusterAnode clusterTime %d %f %f %f \n",j+1,clusterMult, clusterAnode, clusterTime);
273
274 AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,0.,0.,clusterDriftPath,clusteranodePath,clusterMult);
275 //fClusters->Add(point);
276 iTS->AddCluster(1,clust);
277 it = tstop;
278 } // ilcl
279
280 it++;
281
282 } // while (samples)
283 } // anodes
284 } // detectors (2)
285
286 Int_t nofClusters = fClusters->GetEntriesFast();
287 nofClusters -= fNclusters;
288
289 //printf("SDD- Find1Dclust: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);
290
291 fMap->ClearMap();
44b3710f 292
293 for(i=0;i<fNofMaps;i++) delete[] dfadc[i];
294 delete [] dfadc;
b0f5e3fc 295
296 return;
297
298}
299
300//_____________________________________________________________________________
301void AliITSClusterFinderSDD::GroupClusters()
302{
303 // group clusters
304 Int_t dummy=0;
305 Float_t fTimeStep = fSegmentation->Dpx(dummy);
306
307
308 // get number of clusters for this module
309 Int_t nofClusters = fClusters->GetEntriesFast();
310 nofClusters -= fNclusters;
311
312 //printf("SDD- GroupClusters: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);
313
314 AliITSRawClusterSDD *clusterI;
315 AliITSRawClusterSDD *clusterJ;
316
44b3710f 317 Int_t *label = new Int_t [nofClusters];
b0f5e3fc 318 Int_t i,j;
319 for(i=0; i<nofClusters; i++) label[i] = 0;
320 for(i=0; i<nofClusters; i++) {
321 if(label[i] != 0) continue;
322 for(j=i+1; j<nofClusters; j++) {
323 if(label[j] != 0) continue;
324 clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
325 clusterJ = (AliITSRawClusterSDD*) fClusters->At(j);
326 // 1.3 good
327 if(clusterI->T() < fTimeStep*60) fDAnode = 3.2;
328 if(clusterI->T() < fTimeStep*10) fDAnode = 1.2;
329 Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
330 if(!pair) continue;
331 // clusterI->Print();
332 // clusterJ->Print();
333 clusterI->Add(clusterJ);
334 label[j] = 1;
335 fClusters->RemoveAt(j);
336 } // J clusters
337 label[i] = 1;
338 } // I clusters
339 fClusters->Compress();
44b3710f 340
341 delete [] label;
b0f5e3fc 342 return;
343
344}
345
346//_____________________________________________________________________________
347
348void AliITSClusterFinderSDD::SelectClusters()
349{
350 // get number of clusters for this module
351 Int_t nofClusters = fClusters->GetEntriesFast();
352 nofClusters -= fNclusters;
353
354 //printf("SDD- SelectClusters: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);
355
356 Int_t i;
357 for(i=0; i<nofClusters; i++) {
358 AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
359 Int_t rmflg = 0;
360 Float_t wy = 0.;
361 if(clusterI->Anodes() != 0.) {
362 wy = ((Float_t) clusterI->Samples())/clusterI->Anodes();
363 }
364 Float_t amp = clusterI->PeakAmpl();
365 if(amp < fMinPeak) rmflg = 1;
366 if(wy < fMinNCells) rmflg = 1;
367 if(rmflg) fClusters->RemoveAt(i);
368 } // I clusters
369 fClusters->Compress();
370 return;
371
372}
373
374//_____________________________________________________________________________
375
376void AliITSClusterFinderSDD::GetRecPoints()
377{
378 // get rec points
379 //static Int_t counter=0;
380
381 AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
382
383 // get number of clusters for this module
384 Int_t nofClusters = fClusters->GetEntriesFast();
385 nofClusters -= fNclusters;
386
387 //printf("SDD- GetRecPoints: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);
388
389// const Float_t kdEdXtoQ = 2.778e+2; // KeV -> number of e-hole pairs in Si
390 const Float_t kconvGeV = 1.e-6; // GeV -> KeV
391 const Float_t kconv = 1.0e-4;
392 const Float_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
393 const Float_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3
394
395 Int_t i;
396 for(i=0; i<nofClusters; i++) {
397 AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
398 AliITSRecPoint rnew;
399 rnew.SetX(clusterI->X());
400 rnew.SetZ(clusterI->Z());
401 rnew.SetQ(clusterI->Q()); // in KeV - should be ADC
402 //rnew.SetdEdX((clusterI->Q())/kdEdXtoQ);
403 rnew.SetdEdX(kconvGeV*clusterI->Q());
404 rnew.SetSigmaX2(kRMSx*kRMSx);
405 rnew.SetSigmaZ2(kRMSz*kRMSz);
406 rnew.SetProbability(1.);
407 iTS->AddRecPoint(rnew);
408 //counter++;
409 } // I clusters
410 //printf("counter %d\n",counter);
411
412}
413
414//_____________________________________________________________________________
415
416void AliITSClusterFinderSDD::FindRawClusters()
417{
418 // find raw clusters
419 Find1DClusters();
420 GroupClusters();
421 SelectClusters();
422 GetRecPoints();
423}
424
425
426
427