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