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