TRD/AliTRDmcm.cxx

   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$ */
  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
  33 ClassImp(AliTRDmcm)
  34
  35 //_____________________________________________________________________________
  36 AliTRDmcm::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 //_____________________________________________________________________________
  78 AliTRDmcm::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 //_____________________________________________________________________________
 127 AliTRDmcm::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 //_____________________________________________________________________________
 169 AliTRDmcm::~AliTRDmcm()
 170 {
 171   //
 172   // AliTRDmcm destructor
 173   //
 174
 175 }
 176
 177 //_____________________________________________________________________________
 178 AliTRDmcm &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 //_____________________________________________________________________________
 191 void 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 //_____________________________________________________________________________
 233 void 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 //_____________________________________________________________________________
 247 void 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 //_____________________________________________________________________________
 270 Bool_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 //_____________________________________________________________________________
 347 Int_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 //_____________________________________________________________________________
 434 void 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 //_____________________________________________________________________________
 488 void 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 //_____________________________________________________________________________
 503 Bool_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 //_____________________________________________________________________________
 529 void 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 //_____________________________________________________________________________
 593 void 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 //_____________________________________________________________________________
 648 void 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 //_____________________________________________________________________________
 773 void 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 //______________________________________________________________________________
 804 void 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 //______________________________________________________________________________
 855 void 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