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