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