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