]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/attic/AliTRDmcm.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding inteface for High voltage
[u/mrichter/AliRoot.git] / TRD / attic / AliTRDmcm.cxx
CommitLineData
7aa4d4fd 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/* $Id: AliTRDmcm.cxx 29514 2008-10-26 10:24:38Z hristov $ */
17
18////////////////////////////////////////////////////////////////////////////
19// //
20// //
21// Multi Chip Module exponential filter and tracklet finder //
22// //
23// //
24////////////////////////////////////////////////////////////////////////////
25
26#include <TMath.h>
27
28#include "AliLog.h"
29
30#include "AliTRDmcm.h"
31#include "AliTRDtrigParam.h"
32
33ClassImp(AliTRDmcm)
34
35//_____________________________________________________________________________
36AliTRDmcm::AliTRDmcm()
37 :TObject()
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)
55{
56 //
57 // AliTRDmcm default constructor
58 //
59
60 Int_t i = 0;
61 Int_t j = 0;
62
63 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
64 fTrkIndex[i] = 0;
65 fSeedCol[i] = -1;
66 }
67 for (i = 0; i < kMcmCol; i++) {
68 fPadHits[i] = 0;
69 for (j = 0; j < kMcmTBmax; j++) {
70 fADC[i][j] = 0.0;
71 fIsClus[i][j] = kFALSE;
72 }
73 }
74
75}
76
77//_____________________________________________________________________________
78AliTRDmcm::AliTRDmcm(Int_t id)
79 :TObject()
80 ,fNtrk(0)
81 ,fRobId(0)
82 ,fChaId(0)
83 ,fRow(0)
84 ,fColFirst(0)
85 ,fColLast(0)
86 ,fTime1(0)
87 ,fTime2(0)
88 ,fClusThr(0)
89 ,fPadThr(0)
90 ,fNtrkSeeds(0)
91 ,fR1(0)
92 ,fR2(0)
93 ,fC1(0)
94 ,fC2(0)
95 ,fPedestal(0)
96 ,fId(id)
97{
98 //
99 // AliTRDmcm constructor
100 //
101
102 Int_t i = 0;
103 Int_t j = 0;
104
105 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
106 fTrkIndex[i] = 0;
107 fSeedCol[i] = -1;
108 }
109 for (i = 0; i < kMcmCol; i++) {
110 fPadHits[i] = 0;
111 for (j = 0; j < kMcmTBmax; j++) {
112 fADC[i][j] = 0.0;
113 fIsClus[i][j] = kFALSE;
114 }
115 }
116
117 fTime1 = AliTRDtrigParam::Instance()->GetTime1();
118 fTime2 = AliTRDtrigParam::Instance()->GetTime2();
119 fClusThr = AliTRDtrigParam::Instance()->GetClusThr();
120 fPadThr = AliTRDtrigParam::Instance()->GetPadThr();
121
122 AliTRDtrigParam::Instance()->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
123
124}
125
126//_____________________________________________________________________________
127AliTRDmcm::AliTRDmcm(const AliTRDmcm &m)
128 :TObject(m)
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 }
165
166}
167
168//_____________________________________________________________________________
169AliTRDmcm::~AliTRDmcm()
170{
171 //
172 // AliTRDmcm destructor
173 //
174
175}
176
177//_____________________________________________________________________________
178AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
179{
180 //
181 // Assignment operator
182 //
183
184 if (this != &m) ((AliTRDmcm &) m).Copy(*this);
185
186 return *this;
187
188}
189
190//_____________________________________________________________________________
191void AliTRDmcm::Copy(TObject &m) const
192{
193 //
194 // Copy function
195 //
196
197 Int_t i = 0;
198 Int_t j = 0;
199
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
218 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
219 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
220 ((AliTRDmcm &) m).fSeedCol[i] = -1;
221 }
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;
226 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
227 }
228 }
229
230}
231
232//_____________________________________________________________________________
233void AliTRDmcm::AddTrk(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
253 Int_t i = 0;
254 Int_t j = 0;
255
256 for (i = 0; i < kMcmCol; i++) {
257 fPadHits[i] = 0;
258 for (j = 0; j < kMcmTBmax; j++) {
259 fADC[i][j] = 0.0;
260 fIsClus[i][j] = kFALSE;
261 }
262 }
263 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
264 fSeedCol[i] = -1;
265 }
266
267}
268
269//_____________________________________________________________________________
270Bool_t AliTRDmcm::Run()
271{
272 //
273 // Run MCM
274 //
275
276 AliDebug(2,(Form("Run MCM %d\n",Id())));
277
278 Int_t iTime = 0;
279 Int_t iCol = 0;
280 Int_t iPlus = 0;
281 Int_t i = 0;
282 Int_t j = 0;
283
284 Float_t amp[3] = { 0.0, 0.0, 0.0 };
285 Int_t nClus;
286 Int_t clusCol[kMcmCol/2];
287 Float_t clusAmp[kMcmCol/2];
288 Float_t veryLarge;
289 Int_t clusMin = -1;
290
291 // Main TB loop
292 for (iTime = fTime1; iTime <= fTime2; iTime++) {
293
294 // Find clusters...
295 nClus = 0;
296 for (iCol = 1; iCol < (kMcmCol-1); iCol++) {
297 amp[0] = fADC[iCol-1][iTime];
298 amp[1] = fADC[iCol ][iTime];
299 amp[2] = fADC[iCol+1][iTime];
300 if (IsCluster(amp)) {
301 fIsClus[iCol][iTime] = kTRUE;
302 clusCol[nClus] = iCol;
303 clusAmp[nClus] = amp[0]+amp[1]+amp[2];
304 nClus++;
305 if (nClus == kMcmCol/2) {
306 AliWarning(Form("Too many clusters in time bin %2d MCM %d...\n",iTime,Id()));
307 break;
308 }
309 }
310 }
311
312 // ...but no more than six...
313 if (nClus > (Int_t) kSelClus) {
314 for (j = kSelClus/2; j < nClus-kSelClus/2; j++) {
315 fIsClus[clusCol[j]][iTime] = kFALSE;
316 }
317 }
318
319 // ...and take the largest four.
320 Int_t nClusPlus = nClus - kMaxClus;
321 for (iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
322 veryLarge = 1.E+10;
323 for (i = 0; i < nClus; i++) {
324 if (fIsClus[clusCol[i]][iTime]) {
325 if (clusAmp[i] <= veryLarge) {
326 veryLarge = clusAmp[i];
327 clusMin = i;
328 }
329 }
330 }
331 fIsClus[clusCol[clusMin]][iTime] = kFALSE;
332 }
333
334 AddTimeBin(iTime);
335
336 } // end main TB loop
337
338 if ((fNtrkSeeds = CreateSeeds())) {
339 return kTRUE;
340 }
341
342 return kFALSE;
343
344}
345
346//_____________________________________________________________________________
347Int_t AliTRDmcm::CreateSeeds()
348{
349 //
350 // Make column seeds (from Falk Lesser, ex KIP)
351 //
352
353 Int_t i = 0;
354 Int_t j = 0;
355 Int_t nSeeds = 0;
356
357 AliDebug(2,Form("AliTRDmcm::CreateSeeds MCM %d \n",Id()));
358
359 // Working array for hit sums
360 Int_t fHit2padSum[2][kMcmCol];
361
362 // Initialize the array
363 for (i = 0; i < 2; i++) {
364 for (j = 0; j < kMcmCol; j++ ) {
365 if (i == 0) {
366 fHit2padSum[i][j] = j;
367 } else {
368 fHit2padSum[i][j] = -1;
369 }
370 }
371 }
372
373 Int_t sum10 = AliTRDtrigParam::Instance()->GetSum10();
374 Int_t sum12 = AliTRDtrigParam::Instance()->GetSum12();
375
376 // Build the 2padSum
377 Int_t nsum2seed = 0;
378 for (i = 0; i < kMcmCol; i++) {
379 if (i < (kMcmCol-1)) {
380 if ((fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12)) {
381 fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1];
382 }
383 else {
384 fHit2padSum[1][i] = -1;
385 }
386 }
387 else {
388 if (fPadHits[i] >= sum12) {
389 fHit2padSum[1][i] = fPadHits[i];
390 }
391 else {
392 fHit2padSum[1][i] = -1;
393 }
394 }
395 if (fHit2padSum[1][i] > 0) {
396 nsum2seed++;
397 }
398 }
399
400 // sort the sums in decreasing order of the amplitude
401 Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
402
403 // arrange (maximum number of) candidates in increasing order of the column number
404 nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
405 Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
406
407 for (i = 0; i < nSeeds; i++) {
408 fSeedCol[i] = fHit2padSum[0][i];
409 }
410
411 // reject multiple found tracklets
412 Int_t imax = nSeeds - 1;
413 for (i = 0; i < imax; i++) {
414
415 if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
416 nSeeds--;
417 if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
418 AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]));
419 fSeedCol[i+1] = -1;
420 }
421 else {
422 AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]));
423 fSeedCol[i] = -1;
424 }
425 }
426
427 }
428
429 return nSeeds;
430
431}
432
433//_____________________________________________________________________________
434void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir) const
435{
436 //
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
440 //
441
442 Int_t i = 0;
443 Bool_t sort;
444 Int_t tmp1;
445 Int_t tmp2;
446
447 if (dir == 0) {
448
449 do {
450 sort = kTRUE;
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];
455 x2[i+1] = tmp2;
456 tmp1 = x1[i];
457 x1[i] = x1[i+1];
458 x1[i+1] = tmp1;
459 sort = kFALSE;
460 }
461 }
462 } while (sort == kFALSE);
463
464 }
465
466 if (dir == 1) {
467
468 do {
469 sort = kTRUE;
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];
474 x2[i+1] = tmp2;
475 tmp1 = x1[i];
476 x1[i] = x1[i+1];
477 x1[i+1] = tmp1;
478 sort = kFALSE;
479 }
480 }
481 } while (sort == kFALSE);
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//_____________________________________________________________________________
503Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
504{
505 //
506 // Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
507 //
508
509 // -> shape
510 if (amp[0] > amp[1] || amp[2] > amp[1]) {
511 return kFALSE;
512 }
513
514 // -> cluster amplitude
515 if ((amp[0]+amp[1]+amp[2]) < fClusThr) {
516 return kFALSE;
517 }
518
519 // -> pad amplitude
520 if (amp[0] < fPadThr && amp[2] < fPadThr) {
521 return kFALSE;
522 }
523
524 return kTRUE;
525
526}
527
528//_____________________________________________________________________________
529void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
530{
531 //
532 // Exponential filter
533 //
534
535 Int_t iCol = 0;
536 Int_t iTime = 0;
537
538 Double_t sour[kMcmTBmax];
539 Double_t dtarg[kMcmTBmax];
540 Int_t itarg[kMcmTBmax];
541
542 switch(ftype) {
543
544 case 0:
545
546 for (iCol = 0; iCol < kMcmCol; iCol++) {
547 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
548 sour[iTime] = fADC[iCol][iTime];
549 }
550 DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
551 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
552 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
553 }
554 }
555 break;
556
557 case 1:
558
559 for (iCol = 0; iCol < kMcmCol; iCol++) {
560 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
561 sour[iTime] = fADC[iCol][iTime];
562 }
563 DeConvExpD(sour,itarg,kMcmTBmax,nexp);
564 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
565 fADC[iCol][iTime] = itarg[iTime];
566 }
567 }
568 break;
569
570 case 2:
571
572 for (iCol = 0; iCol < kMcmCol; iCol++) {
573 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
574 sour[iTime] = fADC[iCol][iTime];
575 }
576 DeConvExpMI(sour,dtarg,kMcmTBmax);
577 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
578 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
579 }
580 }
581 break;
582
583 default:
584
585 AliError(Form("Invalid filter type %d ! \n",ftype));
586 return;
587
588 }
589
590}
591
592//_____________________________________________________________________________
593void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp)
594{
595 //
596 // Exponential filter "analog"
597 //
598
599 Int_t i = 0;
600 Int_t k = 0;
601 Double_t reminder[2];
602 Double_t correction;
603 Double_t result;
604 Double_t rates[2];
605 Double_t coefficients[2];
606
607 // Initialize (coefficient = alpha, rates = lambda)
608
609 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
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;
615
616 coefficients[0] = c1;
617 coefficients[1] = c2;
618
619 Double_t dt = 0.1;
620 rates[0] = TMath::Exp(-dt/(r1));
621 rates[1] = TMath::Exp(-dt/(r2));
622
623 // Attention: computation order is important
624 correction = 0.0;
625 for (k = 0; k < nexp; k++) {
626 reminder[k] = 0.0;
627 }
628
629 for (i = 0; i < n; i++) {
630
631 result = (source[i] - correction); // no rescaling
632 target[i] = result;
633
634 for (k = 0; k < nexp; k++) {
635 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
636 }
637
638 correction = 0.0;
639 for (k = 0; k < nexp; k++) {
640 correction += reminder[k];
641 }
642
643 }
644
645}
646
647//_____________________________________________________________________________
648void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp)
649{
650 //
651 // Exponential filter "digital"
652 //
653
654 Int_t i = 0;
655
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;
666
667 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
668 // initialize (coefficient = alpha, rates = lambda)
669
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);
692 iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
693 iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
694
695 if (nexp == 1) {
696 fAlphaL = (Int_t) (c1 * 2048.0);
697 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
698 }
699 if (nexp == 2) {
700 fAlphaL = (Int_t) (c1 * 2048.0);
701 fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
702 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
703 iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
704 }
705
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
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
710
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));
722
723 // further initialization
724 if ((rem1 + rem2) > 0x0FFF) {
725 correction = 0x0FFF;
726 }
727 else {
728 correction = (rem1 + rem2) & 0x0FFF;
729 }
730
731 fTailPed = iFactor - correction;
732
733 for (i = 0; i < n; i++) {
734
735 result = ((Int_t)source[i] - correction);
736 if (result < 0) {
737 result = 0;
738 }
739
740 target[i] = result;
741
742 h1 = (rem1 + ((iAlphaL * result) >> 11));
743 if (h1 > 0x0FFF) {
744 h1 = 0x0FFF;
745 }
746 else {
747 h1 &= 0x0FFF;
748 }
749
750 h2 = (rem2 + ((iAlphaS * result) >> 11));
751 if (h2 > 0x0FFF) {
752 h2 = 0x0FFF;
753 }
754 else {
755 h2 &= 0x0FFF;
756 }
757
758 rem1 = (iLambdaL * h1 ) >> 11;
759 rem2 = (iLambdaS * h2 ) >> 11;
760
761 if ((rem1 + rem2) > 0x0FFF) {
762 correction = 0x0FFF;
763 }
764 else {
765 correction = (rem1 + rem2) & 0x0FFF;
766 }
767
768 }
769
770}
771
772//_____________________________________________________________________________
773void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n)
774{
775 //
776 // Exponential filter (M. Ivanov)
777 //
778
779 Int_t i = 0;
780
781 Double_t sig1[100];
782 Double_t sig2[100];
783 Double_t sig3[100];
784
785 for (i = 0; i < n; i++) {
786 sig1[i] = source[i];
787 }
788
789 Float_t dt = 0.1;
790
791 Float_t lambda0 = (1.0 / fR2) * dt;
792 Float_t lambda1 = (1.0 / fR1) * dt;
793
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 }
800
801}
802
803//______________________________________________________________________________
804void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
805{
806 //
807 // Special filter (M. Ivanov)
808 //
809
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];
824 in[1] = (ampin[0] + ampin[1]) * 0.5;
825
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
839 Double_t temp;
840 out[2*n] = in[2*n];
841 temp = 0;
842 for (i = 2*n; i >= 0; i--) {
843 out[i] = in[i] + temp;
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 }
851
852}
853
854//______________________________________________________________________________
855void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout
856 , Double_t norm, Double_t lambda, Int_t n)
857{
858 //
859 // Special filter (M. Ivanov)
860 //
861
862 Int_t i = 0;
863
864 Double_t l = TMath::Exp(-lambda*0.5);
865 Double_t k = l*(1.0 - norm*lambda*0.5);
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];
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;
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
902 //ampout[i] = TMath::Max(out[2*i+1],0.0); // org
903 ampout[i] = TMath::Max(out[2*i],0.0);
904 }
905
906}
907