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 /*
  17 $Log$
  18 Revision 1.3  2006/08/11 17:58:05  cblume
  19 Next round of effc++ changes
  20
  21 Revision 1.2  2006/04/05 12:45:40  hristov
  22 Updated TRD trigger code. Now the trigger code can run:
  23
  24 - in the simulation step, through the central trigger interface, to
  25 produce and store "tracklets" and to produce and analyze tracks, with
  26 inputs to the central trigger
  27
  28 - in the reconstruction step: if the tracklets have already been produced
  29 in the simulation step (by the trigger option), then only the tracks are
  30 produced and stored in the ESD event; if not, the trigger start by
  31 producing the tracklets, etc.
  32
  33 Bogdan
  34
  35 Revision 1.1.1.1  2004/08/19 14:58:11  vulpescu
  36 CVS head
  37
  38 Revision 1.1.1.1  2004/08/18 07:47:17  vulpescu
  39 test
  40
  41 */
  42
  43 ///////////////////////////////////////////////////////////////////////////////
  44 //                                                                           //
  45 //                                                                           //
  46 //  Multi Chip Module exponential filter and tracklet finder                 //
  47 //                                                                           //
  48 //                                                                           //
  49 ///////////////////////////////////////////////////////////////////////////////
  50
  51 #include <TMath.h>
  52
  53 #include "AliLog.h"
  54
  55 #include "AliTRDmcm.h"
  56 #include "AliTRDtrigParam.h"
  57
  58 ClassImp(AliTRDmcm)
  59
  60 //_____________________________________________________________________________
  61 AliTRDmcm::AliTRDmcm()
  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)
  81 {
  82   //
  83   // AliTRDmcm default constructor
  84   //
  85
  86   Int_t i = 0;
  87   Int_t j = 0;
  88
  89   for (i = 0; i < kMaxTrackletsPerMCM; i++) {
  90     fTrkIndex[i] = 0;
  91     fSeedCol[i]  = -1;
  92   }
  93   for (i = 0; i < kMcmCol; i++) {
  94     fPadHits[i]  = 0;
  95     for (j = 0; j < kMcmTBmax; j++) {
  96       fADC[i][j]    = 0.0;
  97       fIsClus[i][j] = kFALSE;
  98     }
  99   }
 100
 101 }
 102
 103 //_____________________________________________________________________________
 104 AliTRDmcm::AliTRDmcm(AliTRDtrigParam *trigp, const Int_t id)
 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)
 124 {
 125   //
 126   // AliTRDmcm constructor
 127   //
 128
 129   Int_t i = 0;
 130   Int_t j = 0;
 131
 132   for (i = 0; i < kMaxTrackletsPerMCM; i++) {
 133     fTrkIndex[i] = 0;
 134     fSeedCol[i]  = -1;
 135   }
 136   for (i = 0; i < kMcmCol; i++) {
 137     fPadHits[i]  = 0;
 138     for (j = 0; j < kMcmTBmax; j++) {
 139       fADC[i][j]    = 0.0;
 140       fIsClus[i][j] = kFALSE;
 141     }
 142   }
 143
 144   fTrigParam->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
 145
 146 }
 147
 148 //_____________________________________________________________________________
 149 AliTRDmcm::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   }
 188
 189 }
 190
 191 //_____________________________________________________________________________
 192 AliTRDmcm::~AliTRDmcm()
 193 {
 194   //
 195   // AliTRDmcm destructor
 196   //
 197
 198 }
 199
 200 //_____________________________________________________________________________
 201 AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
 202 {
 203   //
 204   // Assignment operator
 205   //
 206
 207   if (this != &m) ((AliTRDmcm &) m).Copy(*this);
 208
 209   return *this;
 210
 211 }
 212
 213 //_____________________________________________________________________________
 214 void AliTRDmcm::Copy(TObject &m) const
 215 {
 216   //
 217   // Copy function
 218   //
 219
 220   Int_t i = 0;
 221   Int_t j = 0;
 222
 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
 242   for (i = 0; i < kMaxTrackletsPerMCM; i++) {
 243     ((AliTRDmcm &) m).fTrkIndex[i] = 0;
 244     ((AliTRDmcm &) m).fSeedCol[i]  = -1;
 245   }
 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;
 250       ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
 251     }
 252   }
 253
 254 }
 255
 256 //_____________________________________________________________________________
 257 void 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 //_____________________________________________________________________________
 271 void AliTRDmcm::Reset()
 272 {
 273   //
 274   // Reset MCM data
 275   //
 276
 277   Int_t i = 0;
 278   Int_t j = 0;
 279
 280   for (i = 0; i < kMcmCol; i++) {
 281     fPadHits[i] = 0;
 282     for (j = 0; j < kMcmTBmax; j++) {
 283       fADC[i][j]    = 0.0;
 284       fIsClus[i][j] = kFALSE;
 285     }
 286   }
 287   for (i = 0; i < kMaxTrackletsPerMCM; i++) {
 288     fSeedCol[i] = -1;
 289   }
 290
 291 }
 292
 293 //_____________________________________________________________________________
 294 Bool_t AliTRDmcm::Run()
 295 {
 296   //
 297   // Run MCM
 298   //
 299
 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;
 308
 309   Float_t amp[3] = { 0.0, 0.0, 0.0 };
 310   Int_t   nClus;
 311   Int_t   clusCol[kMcmCol/2];
 312   Float_t clusAmp[kMcmCol/2];
 313   Float_t veryLarge;
 314   Int_t   clusMin = -1;
 315
 316   // Main TB loop
 317   for (iTime = fTime1; iTime <= fTime2; iTime++) {
 318
 319     // Find clusters...
 320     nClus = 0;
 321     for (iCol = 1; iCol < (kMcmCol-1); iCol++) {
 322       amp[0] = fADC[iCol-1][iTime];
 323       amp[1] = fADC[iCol  ][iTime];
 324       amp[2] = fADC[iCol+1][iTime];
 325       if (IsCluster(amp)) {
 326         fIsClus[iCol][iTime] = kTRUE;
 327         clusCol[nClus]       = iCol;
 328         clusAmp[nClus]       = amp[0]+amp[1]+amp[2];
 329         nClus++;
 330         if (nClus == kMcmCol/2) {
 331           AliWarning(Form("Too many clusters in time bin %2d MCM %d...\n",iTime,Id()));
 332           break;
 333         }
 334       }
 335     }
 336
 337     // ...but no more than six...
 338     if (nClus > (Int_t) kSelClus) {
 339       for (j = kSelClus/2; j < nClus-kSelClus/2; j++) {
 340         fIsClus[clusCol[j]][iTime] = kFALSE;
 341       }
 342     }
 343
 344     // ...and take the largest four.
 345     Int_t nClusPlus = nClus - kMaxClus;
 346     for (iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
 347       veryLarge = 1.E+10;
 348       for (i = 0; i < nClus; i++) {
 349         if (fIsClus[clusCol[i]][iTime]) {
 350           if (clusAmp[i] <= veryLarge) {
 351             veryLarge = clusAmp[i];
 352             clusMin = i;
 353           }
 354         }
 355       }
 356       fIsClus[clusCol[clusMin]][iTime] = kFALSE;
 357     }
 358
 359     AddTimeBin(iTime);
 360
 361   }  // end main TB loop
 362
 363   if ((fNtrkSeeds = CreateSeeds())) {
 364     return kTRUE;
 365   }
 366
 367   return kFALSE;
 368
 369 }
 370
 371 //_____________________________________________________________________________
 372 Int_t AliTRDmcm::CreateSeeds()
 373 {
 374   //
 375   // Make column seeds (from Falk Lesser, ex KIP)
 376   //
 377
 378   Int_t i      = 0;
 379   Int_t j      = 0;
 380   Int_t nSeeds = 0;
 381
 382   AliDebug(2,Form("AliTRDmcm::CreateSeeds MCM %d \n",Id()));
 383
 384   // Working array for hit sums
 385   Int_t fHit2padSum[2][kMcmCol];
 386
 387   // Initialize the array
 388   for (i = 0; i < 2; i++) {
 389     for (j = 0; j < kMcmCol; j++ ) {
 390       if (i == 0) {
 391         fHit2padSum[i][j] = j;
 392       } else {
 393         fHit2padSum[i][j] = -1;
 394       }
 395     }
 396   }
 397
 398   Int_t sum10 = fTrigParam->GetSum10();
 399   Int_t sum12 = fTrigParam->GetSum12();
 400
 401   // Build the 2padSum
 402   Int_t nsum2seed = 0;
 403   for (i = 0; i < kMcmCol; i++) {
 404     if (i < (kMcmCol-1)) {
 405       if ((fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12)) {
 406         fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1];
 407       }
 408       else {
 409         fHit2padSum[1][i] = -1;
 410       }
 411     }
 412     else {
 413       if (fPadHits[i] >= sum12) {
 414         fHit2padSum[1][i] = fPadHits[i];
 415       }
 416       else {
 417         fHit2padSum[1][i] = -1;
 418       }
 419     }
 420     if (fHit2padSum[1][i] > 0) {
 421       nsum2seed++;
 422     }
 423   }
 424
 425   // sort the sums in decreasing order of the amplitude
 426   Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
 427
 428   // arrange (maximum number of) candidates in increasing order of the column number
 429   nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
 430   Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
 431
 432   for (i = 0; i < nSeeds; i++) {
 433     fSeedCol[i] = fHit2padSum[0][i];
 434   }
 435
 436   // reject multiple found tracklets
 437   Int_t imax = nSeeds - 1;
 438   for (i = 0; i < imax; i++) {
 439
 440     if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
 441       nSeeds--;
 442       if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
 443         if (fTrigParam->GetDebugLevel() > 1) {
 444           AliWarning(Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]));
 445         }
 446         fSeedCol[i+1] = -1;
 447       }
 448       else {
 449         if (fTrigParam->GetDebugLevel() > 1) {
 450           AliWarning(Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]));
 451         }
 452         fSeedCol[i]   = -1;
 453       }
 454     }
 455
 456   }
 457
 458   return nSeeds;
 459
 460 }
 461
 462 //_____________________________________________________________________________
 463 void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir)
 464 {
 465   //
 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
 469   //
 470
 471   Int_t  i = 0;
 472   Bool_t sort;
 473   Int_t  tmp1;
 474   Int_t  tmp2;
 475
 476   if (dir == 0) {
 477
 478     do {
 479       sort = kTRUE;
 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];
 484           x2[i+1] = tmp2;
 485           tmp1    = x1[i];
 486           x1[i]   = x1[i+1];
 487           x1[i+1] = tmp1;
 488           sort    = kFALSE;
 489         }
 490       }
 491     } while (sort == kFALSE);
 492
 493   }
 494
 495   if (dir == 1) {
 496
 497     do {
 498       sort = kTRUE;
 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];
 503           x2[i+1] = tmp2;
 504           tmp1    = x1[i];
 505           x1[i]   = x1[i+1];
 506           x1[i+1] = tmp1;
 507           sort    = kFALSE;
 508         }
 509       }
 510     } while (sort == kFALSE);
 511
 512   }
 513
 514 }
 515
 516 //_____________________________________________________________________________
 517 void 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 //_____________________________________________________________________________
 532 Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
 533 {
 534   //
 535   // Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
 536   //
 537
 538   // -> shape
 539   if (amp[0] > amp[1] || amp[2] > amp[1]) {
 540     return kFALSE;
 541   }
 542
 543   // -> cluster amplitude
 544   if ((amp[0]+amp[1]+amp[2]) < fClusThr) {
 545     return kFALSE;
 546   }
 547
 548   // -> pad amplitude
 549   if (amp[0] < fPadThr && amp[2] < fPadThr) {
 550     return kFALSE;
 551   }
 552
 553   return kTRUE;
 554
 555 }
 556
 557 //_____________________________________________________________________________
 558 void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
 559 {
 560   //
 561   // Exponential filter
 562   //
 563
 564   Int_t iCol  = 0;
 565   Int_t iTime = 0;
 566
 567   Double_t sour[kMcmTBmax];
 568   Double_t dtarg[kMcmTBmax];
 569   Int_t    itarg[kMcmTBmax];
 570
 571   switch(ftype) {
 572
 573    case 0:
 574
 575     for (iCol = 0; iCol < kMcmCol; iCol++) {
 576       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
 577         sour[iTime] = fADC[iCol][iTime];
 578       }
 579       DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
 580       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
 581         fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
 582       }
 583     }
 584     break;
 585
 586   case 1:
 587
 588     for (iCol = 0; iCol < kMcmCol; iCol++) {
 589       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
 590         sour[iTime] = fADC[iCol][iTime];
 591       }
 592       DeConvExpD(sour,itarg,kMcmTBmax,nexp);
 593       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
 594         fADC[iCol][iTime] = itarg[iTime];
 595       }
 596     }
 597     break;
 598
 599   case 2:
 600
 601     for (iCol = 0; iCol < kMcmCol; iCol++) {
 602       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
 603         sour[iTime] = fADC[iCol][iTime];
 604       }
 605       DeConvExpMI(sour,dtarg,kMcmTBmax);
 606       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
 607         fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
 608       }
 609     }
 610     break;
 611
 612   default:
 613
 614     AliError(Form("Invalid filter type %d ! \n",ftype));
 615     return;
 616
 617   }
 618
 619 }
 620
 621 //_____________________________________________________________________________
 622 void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp)
 623 {
 624   //
 625   // Exponential filter "analog"
 626   //
 627
 628   Int_t    i = 0;
 629   Int_t    k = 0;
 630   Double_t reminder[2];
 631   Double_t correction;
 632   Double_t result;
 633   Double_t rates[2];
 634   Double_t coefficients[2];
 635
 636   // Initialize (coefficient = alpha, rates = lambda)
 637
 638   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 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;
 644
 645   coefficients[0] = c1;
 646   coefficients[1] = c2;
 647
 648   Double_t dt = 0.1;
 649   rates[0] = TMath::Exp(-dt/(r1));
 650   rates[1] = TMath::Exp(-dt/(r2));
 651
 652   // Attention: computation order is important
 653   correction = 0.0;
 654   for (k = 0; k < nexp; k++) {
 655     reminder[k] = 0.0;
 656   }
 657
 658   for (i = 0; i < n; i++) {
 659
 660     result    = (source[i] - correction);    // no rescaling
 661     target[i] = result;
 662
 663     for (k = 0; k < nexp; k++) {
 664       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
 665     }
 666
 667     correction = 0.0;
 668     for (k = 0; k < nexp; k++) {
 669       correction += reminder[k];
 670     }
 671
 672   }
 673
 674 }
 675
 676 //_____________________________________________________________________________
 677 void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp)
 678 {
 679   //
 680   // Exponential filter "digital"
 681   //
 682
 683   Int_t i = 0;
 684
 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;
 695
 696   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 697   // initialize (coefficient = alpha, rates = lambda)
 698
 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);
 721   iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600;     //  9 bit paramter + fixed bits
 722   iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200;     //  9 bit paramter + fixed bits
 723
 724   if (nexp == 1) {
 725     fAlphaL = (Int_t) (c1 * 2048.0);
 726     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
 727   }
 728   if (nexp == 2) {
 729     fAlphaL = (Int_t) (c1 * 2048.0);
 730     fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
 731     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
 732     iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;              // 10 bit paramter + fixed bits
 733   }
 734
 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
 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
 739
 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));
 751
 752   // further initialization
 753   if ((rem1 + rem2) > 0x0FFF) {
 754     correction = 0x0FFF;
 755   }
 756   else {
 757     correction = (rem1 + rem2) & 0x0FFF;
 758   }
 759
 760   fTailPed = iFactor - correction;
 761
 762   for (i = 0; i < n; i++) {
 763
 764     result = ((Int_t)source[i] - correction);
 765     if (result < 0) {
 766       result = 0;
 767     }
 768
 769     target[i] = result;
 770
 771     h1 = (rem1 + ((iAlphaL * result) >> 11));
 772     if (h1 > 0x0FFF) {
 773       h1 = 0x0FFF;
 774     }
 775     else {
 776       h1 &= 0x0FFF;
 777     }
 778
 779     h2 = (rem2 + ((iAlphaS * result) >> 11));
 780     if (h2 > 0x0FFF) {
 781       h2 = 0x0FFF;
 782     }
 783     else {
 784       h2 &= 0x0FFF;
 785     }
 786
 787     rem1 = (iLambdaL * h1 ) >> 11;
 788     rem2 = (iLambdaS * h2 ) >> 11;
 789
 790     if ((rem1 + rem2) > 0x0FFF) {
 791       correction = 0x0FFF;
 792     }
 793     else {
 794       correction = (rem1 + rem2) & 0x0FFF;
 795     }
 796
 797   }
 798
 799 }
 800
 801 //_____________________________________________________________________________
 802 void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n)
 803 {
 804   //
 805   // Exponential filter (M. Ivanov)
 806   //
 807
 808   Int_t i = 0;
 809
 810   Double_t sig1[100];
 811   Double_t sig2[100];
 812   Double_t sig3[100];
 813
 814   for (i = 0; i < n; i++) {
 815     sig1[i] = source[i];
 816   }
 817
 818   Float_t dt = 0.1;
 819
 820   Float_t lambda0 = (1.0 / fR2) * dt;
 821   Float_t lambda1 = (1.0 / fR1) * dt;
 822
 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   }
 829
 830 }
 831
 832 //______________________________________________________________________________
 833 void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
 834 {
 835   //
 836   // Special filter (M. Ivanov)
 837   //
 838
 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];
 853   in[1] = (ampin[0] + ampin[1]) * 0.5;
 854
 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
 868   Double_t temp;
 869   out[2*n]    = in[2*n];
 870   temp        = 0;
 871   for (i = 2*n; i >= 0; i--) {
 872     out[i]    = in[i] + temp;
 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   }
 880
 881 }
 882
 883 //______________________________________________________________________________
 884 void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout
 885                                 , Double_t norm, Double_t lambda, Int_t n)
 886 {
 887   //
 888   // Special filter (M. Ivanov)
 889   //
 890
 891   Int_t    i = 0;
 892
 893   Double_t l = TMath::Exp(-lambda*0.5);
 894   Double_t k = l*(1.0 - norm*lambda*0.5);
 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];
 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;
 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
 931     //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
 932     ampout[i] = TMath::Max(out[2*i],0.0);
 933   }
 934
 935 }
 936