]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDmcm.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated error definitions.
[u/mrichter/AliRoot.git] / TRD / AliTRDmcm.cxx
CommitLineData
0ee00e25 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
acc49af9 16/* $Id$ */
f162af62 17
acc49af9 18////////////////////////////////////////////////////////////////////////////
19// //
20// //
21// Multi Chip Module exponential filter and tracklet finder //
22// //
23// //
24////////////////////////////////////////////////////////////////////////////
e3b2b5e5 25
0ee00e25 26#include <TMath.h>
0ee00e25 27
6d50f529 28#include "AliLog.h"
29
0ee00e25 30#include "AliTRDmcm.h"
31#include "AliTRDtrigParam.h"
32
33ClassImp(AliTRDmcm)
34
35//_____________________________________________________________________________
36AliTRDmcm::AliTRDmcm()
6d50f529 37 :TObject()
6d50f529 38 ,fNtrk(0)
39 ,fRobId(0)
40 ,fChaId(0)
41 ,fRow(0)
42 ,fColFirst(0)
43 ,fColLast(0)
44 ,fTime1(0)
45 ,fTime2(0)
46 ,fClusThr(0)
47 ,fPadThr(0)
48 ,fNtrkSeeds(0)
49 ,fR1(0)
50 ,fR2(0)
51 ,fC1(0)
52 ,fC2(0)
53 ,fPedestal(0)
54 ,fId(0)
0ee00e25 55{
0ee00e25 56 //
57 // AliTRDmcm default constructor
58 //
59
6d50f529 60 Int_t i = 0;
61 Int_t j = 0;
0ee00e25 62
6d50f529 63 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
0ee00e25 64 fTrkIndex[i] = 0;
6d50f529 65 fSeedCol[i] = -1;
0ee00e25 66 }
6d50f529 67 for (i = 0; i < kMcmCol; i++) {
68 fPadHits[i] = 0;
69 for (j = 0; j < kMcmTBmax; j++) {
70 fADC[i][j] = 0.0;
0ee00e25 71 fIsClus[i][j] = kFALSE;
72 }
73 }
0ee00e25 74
75}
76
77//_____________________________________________________________________________
f162af62 78AliTRDmcm::AliTRDmcm(const Int_t id)
6d50f529 79 :TObject()
6d50f529 80 ,fNtrk(0)
81 ,fRobId(0)
82 ,fChaId(0)
83 ,fRow(0)
84 ,fColFirst(0)
85 ,fColLast(0)
f162af62 86 ,fTime1(0)
87 ,fTime2(0)
88 ,fClusThr(0)
89 ,fPadThr(0)
6d50f529 90 ,fNtrkSeeds(0)
91 ,fR1(0)
92 ,fR2(0)
93 ,fC1(0)
94 ,fC2(0)
95 ,fPedestal(0)
96 ,fId(id)
0ee00e25 97{
98 //
99 // AliTRDmcm constructor
100 //
101
6d50f529 102 Int_t i = 0;
103 Int_t j = 0;
0ee00e25 104
6d50f529 105 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
0ee00e25 106 fTrkIndex[i] = 0;
6d50f529 107 fSeedCol[i] = -1;
0ee00e25 108 }
6d50f529 109 for (i = 0; i < kMcmCol; i++) {
110 fPadHits[i] = 0;
111 for (j = 0; j < kMcmTBmax; j++) {
112 fADC[i][j] = 0.0;
0ee00e25 113 fIsClus[i][j] = kFALSE;
114 }
115 }
f162af62 116
117 fTime1 = AliTRDtrigParam::Instance()->GetTime1();
118 fTime2 = AliTRDtrigParam::Instance()->GetTime2();
119 fClusThr = AliTRDtrigParam::Instance()->GetClusThr();
120 fPadThr = AliTRDtrigParam::Instance()->GetPadThr();
0ee00e25 121
f162af62 122 AliTRDtrigParam::Instance()->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
0ee00e25 123
6d50f529 124}
125
126//_____________________________________________________________________________
127AliTRDmcm::AliTRDmcm(const AliTRDmcm &m)
128 :TObject(m)
6d50f529 129 ,fNtrk(m.fNtrk)
130 ,fRobId(m.fRobId)
131 ,fChaId(m.fChaId)
132 ,fRow(m.fRow)
133 ,fColFirst(m.fColFirst)
134 ,fColLast(m.fColLast)
135 ,fTime1(m.fTime1)
136 ,fTime2(m.fTime2)
137 ,fClusThr(m.fClusThr)
138 ,fPadThr(m.fPadThr)
139 ,fNtrkSeeds(m.fNtrkSeeds)
140 ,fR1(m.fR1)
141 ,fR2(m.fR2)
142 ,fC1(m.fC1)
143 ,fC2(m.fC2)
144 ,fPedestal(m.fPedestal)
145 ,fId(m.fId)
146{
147 //
148 // AliTRDmcm copy constructor
149 //
150
151 Int_t i = 0;
152 Int_t j = 0;
153
154 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
155 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
156 ((AliTRDmcm &) m).fSeedCol[i] = -1;
157 }
158 for (i = 0; i < kMcmCol; i++) {
159 ((AliTRDmcm &) m).fPadHits[i] = 0;
160 for (j = 0; j < kMcmTBmax; j++) {
161 ((AliTRDmcm &) m).fADC[i][j] = 0.0;
162 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
163 }
164 }
0ee00e25 165
166}
167
168//_____________________________________________________________________________
169AliTRDmcm::~AliTRDmcm()
170{
0ee00e25 171 //
172 // AliTRDmcm destructor
173 //
174
175}
176
e3b2b5e5 177//_____________________________________________________________________________
178AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
179{
180 //
6d50f529 181 // Assignment operator
e3b2b5e5 182 //
183
184 if (this != &m) ((AliTRDmcm &) m).Copy(*this);
04eeac11 185
e3b2b5e5 186 return *this;
187
188}
189
190//_____________________________________________________________________________
191void AliTRDmcm::Copy(TObject &m) const
192{
193 //
6d50f529 194 // Copy function
e3b2b5e5 195 //
196
6d50f529 197 Int_t i = 0;
198 Int_t j = 0;
199
e3b2b5e5 200 ((AliTRDmcm &) m).fNtrk = fNtrk;
201 ((AliTRDmcm &) m).fRobId = fRobId;
202 ((AliTRDmcm &) m).fChaId = fChaId;
203 ((AliTRDmcm &) m).fRow = fRow;
204 ((AliTRDmcm &) m).fColFirst = fColFirst;
205 ((AliTRDmcm &) m).fColLast = fColLast;
206 ((AliTRDmcm &) m).fTime1 = fTime1;
207 ((AliTRDmcm &) m).fTime2 = fTime2;
208 ((AliTRDmcm &) m).fClusThr = fClusThr;
209 ((AliTRDmcm &) m).fPadThr = fPadThr;
210 ((AliTRDmcm &) m).fNtrkSeeds = fNtrkSeeds;
211 ((AliTRDmcm &) m).fR1 = fR1;
212 ((AliTRDmcm &) m).fR2 = fR2;
213 ((AliTRDmcm &) m).fC1 = fC1;
214 ((AliTRDmcm &) m).fC2 = fC2;
215 ((AliTRDmcm &) m).fPedestal = fPedestal;
216 ((AliTRDmcm &) m).fId = fId;
217
6d50f529 218 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
e3b2b5e5 219 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
6d50f529 220 ((AliTRDmcm &) m).fSeedCol[i] = -1;
e3b2b5e5 221 }
6d50f529 222 for (i = 0; i < kMcmCol; i++) {
223 ((AliTRDmcm &) m).fPadHits[i] = 0;
224 for (j = 0; j < kMcmTBmax; j++) {
225 ((AliTRDmcm &) m).fADC[i][j] = 0.0;
e3b2b5e5 226 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
227 }
228 }
e3b2b5e5 229
230}
231
0ee00e25 232//_____________________________________________________________________________
233void AliTRDmcm::AddTrk(const Int_t id)
234{
235 //
236 // Add a tracklet index
237 //
238
239 fTrkIndex[fNtrk] = id;
240 fNtrk++;
241
242 return;
243
244}
245
246//_____________________________________________________________________________
247void AliTRDmcm::Reset()
248{
249 //
250 // Reset MCM data
251 //
252
6d50f529 253 Int_t i = 0;
254 Int_t j = 0;
255
256 for (i = 0; i < kMcmCol; i++) {
0ee00e25 257 fPadHits[i] = 0;
6d50f529 258 for (j = 0; j < kMcmTBmax; j++) {
259 fADC[i][j] = 0.0;
0ee00e25 260 fIsClus[i][j] = kFALSE;
261 }
262 }
6d50f529 263 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
0ee00e25 264 fSeedCol[i] = -1;
265 }
266
267}
268
269//_____________________________________________________________________________
270Bool_t AliTRDmcm::Run()
271{
272 //
273 // Run MCM
274 //
275
6d50f529 276 AliDebug(2,(Form("Run MCM %d\n",Id())));
277
278 Int_t iTime = 0;
279 Int_t iCol = 0;
6d50f529 280 Int_t iPlus = 0;
281 Int_t i = 0;
282 Int_t j = 0;
0ee00e25 283
6d50f529 284 Float_t amp[3] = { 0.0, 0.0, 0.0 };
0ee00e25 285 Int_t nClus;
e3b2b5e5 286 Int_t clusCol[kMcmCol/2];
287 Float_t clusAmp[kMcmCol/2];
288 Float_t veryLarge;
289 Int_t clusMin = -1;
0ee00e25 290
6d50f529 291 // Main TB loop
292 for (iTime = fTime1; iTime <= fTime2; iTime++) {
0ee00e25 293
6d50f529 294 // Find clusters...
0ee00e25 295 nClus = 0;
6d50f529 296 for (iCol = 1; iCol < (kMcmCol-1); iCol++) {
e3b2b5e5 297 amp[0] = fADC[iCol-1][iTime];
298 amp[1] = fADC[iCol ][iTime];
299 amp[2] = fADC[iCol+1][iTime];
300 if (IsCluster(amp)) {
0ee00e25 301 fIsClus[iCol][iTime] = kTRUE;
6d50f529 302 clusCol[nClus] = iCol;
303 clusAmp[nClus] = amp[0]+amp[1]+amp[2];
0ee00e25 304 nClus++;
305 if (nClus == kMcmCol/2) {
6d50f529 306 AliWarning(Form("Too many clusters in time bin %2d MCM %d...\n",iTime,Id()));
0ee00e25 307 break;
308 }
309 }
310 }
311
312 // ...but no more than six...
6d50f529 313 if (nClus > (Int_t) kSelClus) {
314 for (j = kSelClus/2; j < nClus-kSelClus/2; j++) {
e3b2b5e5 315 fIsClus[clusCol[j]][iTime] = kFALSE;
0ee00e25 316 }
317 }
318
319 // ...and take the largest four.
0ee00e25 320 Int_t nClusPlus = nClus - kMaxClus;
6d50f529 321 for (iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
e3b2b5e5 322 veryLarge = 1.E+10;
6d50f529 323 for (i = 0; i < nClus; i++) {
e3b2b5e5 324 if (fIsClus[clusCol[i]][iTime]) {
325 if (clusAmp[i] <= veryLarge) {
326 veryLarge = clusAmp[i];
327 clusMin = i;
0ee00e25 328 }
329 }
330 }
e3b2b5e5 331 fIsClus[clusCol[clusMin]][iTime] = kFALSE;
0ee00e25 332 }
333
334 AddTimeBin(iTime);
335
336 } // end main TB loop
04eeac11 337
0ee00e25 338 if ((fNtrkSeeds = CreateSeeds())) {
0ee00e25 339 return kTRUE;
0ee00e25 340 }
341
342 return kFALSE;
343
344}
6d50f529 345
0ee00e25 346//_____________________________________________________________________________
347Int_t AliTRDmcm::CreateSeeds()
348{
349 //
350 // Make column seeds (from Falk Lesser, ex KIP)
351 //
352
6d50f529 353 Int_t i = 0;
354 Int_t j = 0;
0ee00e25 355 Int_t nSeeds = 0;
356
6d50f529 357 AliDebug(2,Form("AliTRDmcm::CreateSeeds MCM %d \n",Id()));
358
359 // Working array for hit sums
0ee00e25 360 Int_t fHit2padSum[2][kMcmCol];
361
6d50f529 362 // Initialize the array
363 for (i = 0; i < 2; i++) {
364 for (j = 0; j < kMcmCol; j++ ) {
365 if (i == 0) {
0ee00e25 366 fHit2padSum[i][j] = j;
367 } else {
368 fHit2padSum[i][j] = -1;
369 }
370 }
371 }
372
f162af62 373 Int_t sum10 = AliTRDtrigParam::Instance()->GetSum10();
374 Int_t sum12 = AliTRDtrigParam::Instance()->GetSum12();
0ee00e25 375
6d50f529 376 // Build the 2padSum
e3b2b5e5 377 Int_t nsum2seed = 0;
6d50f529 378 for (i = 0; i < kMcmCol; i++) {
379 if (i < (kMcmCol-1)) {
380 if ((fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12)) {
0ee00e25 381 fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1];
6d50f529 382 }
383 else {
0ee00e25 384 fHit2padSum[1][i] = -1;
385 }
6d50f529 386 }
387 else {
388 if (fPadHits[i] >= sum12) {
0ee00e25 389 fHit2padSum[1][i] = fPadHits[i];
6d50f529 390 }
391 else {
0ee00e25 392 fHit2padSum[1][i] = -1;
393 }
394 }
04eeac11 395 if (fHit2padSum[1][i] > 0) {
396 nsum2seed++;
0ee00e25 397 }
0ee00e25 398 }
399
400 // sort the sums in decreasing order of the amplitude
401 Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
402
0ee00e25 403 // arrange (maximum number of) candidates in increasing order of the column number
e3b2b5e5 404 nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
0ee00e25 405 Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
406
6d50f529 407 for (i = 0; i < nSeeds; i++) {
0ee00e25 408 fSeedCol[i] = fHit2padSum[0][i];
409 }
410
0ee00e25 411 // reject multiple found tracklets
04eeac11 412 Int_t imax = nSeeds - 1;
6d50f529 413 for (i = 0; i < imax; i++) {
0ee00e25 414
415 if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
416 nSeeds--;
417 if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
f162af62 418 AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]));
0ee00e25 419 fSeedCol[i+1] = -1;
04eeac11 420 }
421 else {
f162af62 422 AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]));
6d50f529 423 fSeedCol[i] = -1;
0ee00e25 424 }
425 }
426
427 }
428
0ee00e25 429 return nSeeds;
430
431}
432
433//_____________________________________________________________________________
acc49af9 434void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir) const
0ee00e25 435{
e3b2b5e5 436 //
0ee00e25 437 // Sort two parallel vectors (x1[nel], x2[nel]) after the second one (x2)
438 // in the direction: dir = 0 ascending order
439 // dir = 1 descending order
e3b2b5e5 440 //
0ee00e25 441
6d50f529 442 Int_t i = 0;
0ee00e25 443 Bool_t sort;
6d50f529 444 Int_t tmp1;
445 Int_t tmp2;
0ee00e25 446
6d50f529 447 if (dir == 0) {
0ee00e25 448
449 do {
450 sort = kTRUE;
6d50f529 451 for (i = 0; i < (nel-1); i++) {
452 if (x2[i+1] < x2[i]) {
453 tmp2 = x2[i];
454 x2[i] = x2[i+1];
0ee00e25 455 x2[i+1] = tmp2;
6d50f529 456 tmp1 = x1[i];
457 x1[i] = x1[i+1];
0ee00e25 458 x1[i+1] = tmp1;
6d50f529 459 sort = kFALSE;
0ee00e25 460 }
6d50f529 461 }
462 } while (sort == kFALSE);
0ee00e25 463
464 }
465
6d50f529 466 if (dir == 1) {
0ee00e25 467
468 do {
469 sort = kTRUE;
6d50f529 470 for (i = 0; i < (nel-1); i++) {
471 if (x2[i+1] > x2[i]) {
472 tmp2 = x2[i];
473 x2[i] = x2[i+1];
0ee00e25 474 x2[i+1] = tmp2;
6d50f529 475 tmp1 = x1[i];
476 x1[i] = x1[i+1];
0ee00e25 477 x1[i+1] = tmp1;
6d50f529 478 sort = kFALSE;
0ee00e25 479 }
6d50f529 480 }
481 } while (sort == kFALSE);
0ee00e25 482
483 }
484
485}
486
487//_____________________________________________________________________________
488void AliTRDmcm::AddTimeBin(const Int_t iTime)
489{
490 //
491 // Build column seeds
492 //
493
494 for (Int_t iPad = 1; iPad < (kMcmCol-1); iPad++) {
495 if (fIsClus[iPad][iTime]) {
496 fPadHits[iPad]++;
497 }
498 }
499
500}
501
502//_____________________________________________________________________________
e3b2b5e5 503Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
0ee00e25 504{
505 //
506 // Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
507 //
508
509 // -> shape
6d50f529 510 if (amp[0] > amp[1] || amp[2] > amp[1]) {
511 return kFALSE;
512 }
0ee00e25 513
514 // -> cluster amplitude
6d50f529 515 if ((amp[0]+amp[1]+amp[2]) < fClusThr) {
516 return kFALSE;
517 }
0ee00e25 518
519 // -> pad amplitude
6d50f529 520 if (amp[0] < fPadThr && amp[2] < fPadThr) {
521 return kFALSE;
522 }
0ee00e25 523
524 return kTRUE;
525
526}
527
528//_____________________________________________________________________________
529void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
530{
531 //
6d50f529 532 // Exponential filter
0ee00e25 533 //
534
6d50f529 535 Int_t iCol = 0;
536 Int_t iTime = 0;
537
0ee00e25 538 Double_t sour[kMcmTBmax];
e3b2b5e5 539 Double_t dtarg[kMcmTBmax];
540 Int_t itarg[kMcmTBmax];
0ee00e25 541
542 switch(ftype) {
543
6d50f529 544 case 0:
0ee00e25 545
6d50f529 546 for (iCol = 0; iCol < kMcmCol; iCol++) {
547 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
0ee00e25 548 sour[iTime] = fADC[iCol][iTime];
549 }
e3b2b5e5 550 DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
6d50f529 551 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
e3b2b5e5 552 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
0ee00e25 553 }
554 }
555 break;
556
557 case 1:
558
6d50f529 559 for (iCol = 0; iCol < kMcmCol; iCol++) {
560 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
0ee00e25 561 sour[iTime] = fADC[iCol][iTime];
562 }
e3b2b5e5 563 DeConvExpD(sour,itarg,kMcmTBmax,nexp);
6d50f529 564 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
e3b2b5e5 565 fADC[iCol][iTime] = itarg[iTime];
0ee00e25 566 }
567 }
568 break;
569
570 case 2:
571
6d50f529 572 for (iCol = 0; iCol < kMcmCol; iCol++) {
573 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
0ee00e25 574 sour[iTime] = fADC[iCol][iTime];
575 }
e3b2b5e5 576 DeConvExpMI(sour,dtarg,kMcmTBmax);
6d50f529 577 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
e3b2b5e5 578 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
0ee00e25 579 }
580 }
581 break;
582
583 default:
584
6d50f529 585 AliError(Form("Invalid filter type %d ! \n",ftype));
0ee00e25 586 return;
587
588 }
589
590}
591
592//_____________________________________________________________________________
593void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp)
594{
e3b2b5e5 595 //
596 // Exponential filter "analog"
597 //
0ee00e25 598
6d50f529 599 Int_t i = 0;
600 Int_t k = 0;
601 Double_t reminder[2];
602 Double_t correction;
603 Double_t result;
0ee00e25 604 Double_t rates[2];
605 Double_t coefficients[2];
606
6d50f529 607 // Initialize (coefficient = alpha, rates = lambda)
0ee00e25 608
609 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
e3b2b5e5 610 Double_t r1, r2, c1, c2;
611 r1 = (Double_t)fR1;
612 r2 = (Double_t)fR2;
613 c1 = (Double_t)fC1;
614 c2 = (Double_t)fC2;
0ee00e25 615
e3b2b5e5 616 coefficients[0] = c1;
617 coefficients[1] = c2;
0ee00e25 618
04eeac11 619 Double_t dt = 0.1;
e3b2b5e5 620 rates[0] = TMath::Exp(-dt/(r1));
621 rates[1] = TMath::Exp(-dt/(r2));
0ee00e25 622
6d50f529 623 // Attention: computation order is important
624 correction = 0.0;
625 for (k = 0; k < nexp; k++) {
626 reminder[k] = 0.0;
627 }
0ee00e25 628
6d50f529 629 for (i = 0; i < n; i++) {
630
631 result = (source[i] - correction); // no rescaling
0ee00e25 632 target[i] = result;
633
6d50f529 634 for (k = 0; k < nexp; k++) {
635 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
636 }
0ee00e25 637
6d50f529 638 correction = 0.0;
639 for (k = 0; k < nexp; k++) {
640 correction += reminder[k];
641 }
642
0ee00e25 643 }
644
645}
646
647//_____________________________________________________________________________
e3b2b5e5 648void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp)
649{
650 //
651 // Exponential filter "digital"
652 //
0ee00e25 653
6d50f529 654 Int_t i = 0;
655
04eeac11 656 Int_t fAlphaL;
657 Int_t fAlphaS;
658 Int_t fLambdaL;
659 Int_t fLambdaS;
660 Int_t fTailPed;
661
662 Int_t iAlphaL;
663 Int_t iAlphaS;
664 Int_t iLambdaL;
665 Int_t iLambdaS;
0ee00e25 666
667 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
668 // initialize (coefficient = alpha, rates = lambda)
669
04eeac11 670 fLambdaL = 0;
671 fAlphaL = 0;
672 fLambdaS = 0;
673 fAlphaS = 0;
674 iLambdaL = 0;
675 iAlphaL = 0;
676 iLambdaS = 0;
677 iAlphaS = 0;
678
679 Double_t dt = 0.1;
680
681 Double_t r1;
682 Double_t r2;
683 Double_t c1;
684 Double_t c2;
685 r1 = (Double_t) fR1;
686 r2 = (Double_t) fR2;
687 c1 = (Double_t) fC1;
688 c2 = (Double_t) fC2;
689
690 fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
691 fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
e3b2b5e5 692 iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
693 iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
0ee00e25 694
695 if (nexp == 1) {
04eeac11 696 fAlphaL = (Int_t) (c1 * 2048.0);
e3b2b5e5 697 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
0ee00e25 698 }
699 if (nexp == 2) {
04eeac11 700 fAlphaL = (Int_t) (c1 * 2048.0);
701 fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
e3b2b5e5 702 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
703 iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
0ee00e25 704 }
705
04eeac11 706 Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers
707 Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers
e3b2b5e5 708 Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers
709 Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers
0ee00e25 710
04eeac11 711 Int_t h1;
712 Int_t h2;
713 Int_t rem1;
714 Int_t rem2;
715 Int_t correction;
716 Int_t result;
717 Int_t iFactor = ((Int_t) fPedestal) << 2;
718
719 Double_t xi = 1 - (iLl*iAs + iLs*iAl); // Calculation of equilibrium values of the
720 rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
721 rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
0ee00e25 722
723 // further initialization
6d50f529 724 if ((rem1 + rem2) > 0x0FFF) {
0ee00e25 725 correction = 0x0FFF;
6d50f529 726 }
727 else {
0ee00e25 728 correction = (rem1 + rem2) & 0x0FFF;
729 }
730
731 fTailPed = iFactor - correction;
732
6d50f529 733 for (i = 0; i < n; i++) {
0ee00e25 734
6d50f529 735 result = ((Int_t)source[i] - correction);
736 if (result < 0) {
0ee00e25 737 result = 0;
738 }
739
740 target[i] = result;
741
6d50f529 742 h1 = (rem1 + ((iAlphaL * result) >> 11));
743 if (h1 > 0x0FFF) {
0ee00e25 744 h1 = 0x0FFF;
6d50f529 745 }
746 else {
0ee00e25 747 h1 &= 0x0FFF;
748 }
749
6d50f529 750 h2 = (rem2 + ((iAlphaS * result) >> 11));
751 if (h2 > 0x0FFF) {
0ee00e25 752 h2 = 0x0FFF;
6d50f529 753 }
754 else {
0ee00e25 755 h2 &= 0x0FFF;
756 }
757
e3b2b5e5 758 rem1 = (iLambdaL * h1 ) >> 11;
759 rem2 = (iLambdaS * h2 ) >> 11;
0ee00e25 760
6d50f529 761 if ((rem1 + rem2) > 0x0FFF) {
0ee00e25 762 correction = 0x0FFF;
6d50f529 763 }
764 else {
0ee00e25 765 correction = (rem1 + rem2) & 0x0FFF;
766 }
04eeac11 767
0ee00e25 768 }
769
770}
771
772//_____________________________________________________________________________
e3b2b5e5 773void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n)
774{
775 //
776 // Exponential filter (M. Ivanov)
777 //
0ee00e25 778
6d50f529 779 Int_t i = 0;
0ee00e25 780
6d50f529 781 Double_t sig1[100];
782 Double_t sig2[100];
783 Double_t sig3[100];
0ee00e25 784
6d50f529 785 for (i = 0; i < n; i++) {
786 sig1[i] = source[i];
787 }
0ee00e25 788
04eeac11 789 Float_t dt = 0.1;
0ee00e25 790
04eeac11 791 Float_t lambda0 = (1.0 / fR2) * dt;
792 Float_t lambda1 = (1.0 / fR1) * dt;
0ee00e25 793
6d50f529 794 TailMakerSpline(sig1,sig2,lambda0,n);
795 TailCancelationMI(sig2,sig3,0.7,lambda1,n);
796
797 for (i = 0; i < n; i++) {
798 target[i] = sig3[i];
799 }
0ee00e25 800
801}
802
803//______________________________________________________________________________
e3b2b5e5 804void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
805{
806 //
807 // Special filter (M. Ivanov)
808 //
0ee00e25 809
6d50f529 810 Int_t i = 0;
811
812 Double_t l = TMath::Exp(-lambda*0.5);
813 Double_t in[1000];
814 Double_t out[1000];
815
816 // Initialize in[] and out[] goes 0 ... 2*n+19
817 for (i = 0; i < n*2+20; i++) {
818 in[i] = 0;
819 out[i] = 0;
820 }
821
822 // in[] goes 0, 1
823 in[0] = ampin[0];
04eeac11 824 in[1] = (ampin[0] + ampin[1]) * 0.5;
0ee00e25 825
6d50f529 826 // Add charge to the end
827 for (i = 0; i < 22; i++) {
828 // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
829 in[2*(n-1)+i] = ampin[n-1];
830 }
831
832 // Use arithmetic mean
833 for (i = 1; i < n-1; i++) {
834 // in[] goes 2, 3, ... , 2*n-4, 2*n-3
835 in[2*i] = ampin[i];
836 in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
837 }
838
6d50f529 839 Double_t temp;
840 out[2*n] = in[2*n];
841 temp = 0;
842 for (i = 2*n; i >= 0; i--) {
04eeac11 843 out[i] = in[i] + temp;
6d50f529 844 temp = l*(temp+in[i]);
845 }
846
847 for (i = 0; i < n; i++){
848 //ampout[i] = out[2*i+1]; // org
849 ampout[i] = out[2*i];
850 }
0ee00e25 851
852}
853
854//______________________________________________________________________________
6d50f529 855void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout
856 , Double_t norm, Double_t lambda, Int_t n)
e3b2b5e5 857{
858 //
859 // Special filter (M. Ivanov)
860 //
0ee00e25 861
6d50f529 862 Int_t i = 0;
863
864 Double_t l = TMath::Exp(-lambda*0.5);
04eeac11 865 Double_t k = l*(1.0 - norm*lambda*0.5);
6d50f529 866 Double_t in[1000];
867 Double_t out[1000];
868
869 // Initialize in[] and out[] goes 0 ... 2*n+19
870 for (i = 0; i < n*2+20; i++) {
871 in[i] = 0;
872 out[i] = 0;
873 }
874
875 // in[] goes 0, 1
876 in[0] = ampin[0];
877 in[1] = (ampin[0]+ampin[1])*0.5;
878
879 // Add charge to the end
880 for (i =-2; i < 22; i++) {
881 // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
882 in[2*(n-1)+i] = ampin[n-1];
883 }
884
885 for (i = 1; i < n-2; i++) {
886 // in[] goes 2, 3, ... , 2*n-6, 2*n-5
887 in[2*i] = ampin[i];
04eeac11 888 in[2*i+1] = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
889 //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.0;
6d50f529 890 }
891
892 Double_t temp;
893 out[0] = in[0];
894 temp = in[0];
895 for (i = 1; i <= 2*n; i++) {
896 out[i] = in[i] + (k-l)*temp;
897 temp = in[i] + k *temp;
898 }
899
900 for (i = 0; i < n; i++) {
901 //ampout[i] = out[2*i+1]; // org
04eeac11 902 //ampout[i] = TMath::Max(out[2*i+1],0.0); // org
903 ampout[i] = TMath::Max(out[2*i],0.0);
6d50f529 904 }
0ee00e25 905
906}
907