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