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.1.1.1  2004/08/19 14:58:11  vulpescu
  19 CVS head
  20
  21 Revision 1.1.1.1  2004/08/18 07:47:17  vulpescu
  22 test
  23
  24 */
  25
  26 ///////////////////////////////////////////////////////////////////////////////
  27 //                                                                           //
  28 //                                                                           //
  29 //  Multi Chip Module exponential filter and tracklet finder                 //
  30 //                                                                           //
  31 //                                                                           //
  32 ///////////////////////////////////////////////////////////////////////////////
  33
  34 #include <TMath.h>
  35
  36 #include "AliTRDmcm.h"
  37 #include "AliTRDtrigParam.h"
  38
  39 ClassImp(AliTRDmcm)
  40
  41 //_____________________________________________________________________________
  42 AliTRDmcm::AliTRDmcm()
  43 {
  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 //_____________________________________________________________________________
  86 AliTRDmcm::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 //_____________________________________________________________________________
 132 AliTRDmcm::~AliTRDmcm()
 133 {
 134   //
 135   // AliTRDmcm destructor
 136   //
 137
 138 }
 139
 140 //_____________________________________________________________________________
 141 AliTRDmcm &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 //_____________________________________________________________________________
 153 void 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
 194 //_____________________________________________________________________________
 195 void 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 //_____________________________________________________________________________
 209 void 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 //_____________________________________________________________________________
 229 Bool_t AliTRDmcm::Run()
 230 {
 231   //
 232   // Run MCM
 233   //
 234
 235   if ( fTrigParam->GetDebugLevel() > 1 ) printf("AliTRDmcm::Run MCM %d\n",Id());
 236
 237   Float_t amp[3] = {0.0, 0.0, 0.0};
 238   Int_t   nClus;
 239   Int_t   clusCol[kMcmCol/2];
 240   Float_t clusAmp[kMcmCol/2];
 241   Float_t veryLarge;
 242   Int_t   clusMin = -1;
 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++) {
 249       amp[0] = fADC[iCol-1][iTime];
 250       amp[1] = fADC[iCol  ][iTime];
 251       amp[2] = fADC[iCol+1][iTime];
 252       if (IsCluster(amp)) {
 253         fIsClus[iCol][iTime] = kTRUE;
 254         clusCol[nClus] = iCol;
 255         clusAmp[nClus] = amp[0]+amp[1]+amp[2];
 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++) {
 268         fIsClus[clusCol[j]][iTime] = kFALSE;
 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++ ) {
 276       veryLarge = 1.E+10;
 277       for (Int_t i = 0; i < nClus; i++) {
 278         if (fIsClus[clusCol[i]][iTime]) {
 279           if (clusAmp[i] <= veryLarge) {
 280             veryLarge = clusAmp[i];
 281             clusMin = i;
 282           }
 283         }
 284       }
 285       fIsClus[clusCol[clusMin]][iTime] = kFALSE;
 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 //_____________________________________________________________________________
 317 Int_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
 341   Int_t sum10 = fTrigParam->GetSum10();
 342   Int_t sum12 = fTrigParam->GetSum12();
 343
 344   // build the 2padSum
 345   Int_t nsum2seed = 0;
 346   for( Int_t i = 0; i < kMcmCol; i++ ) {
 347     if( i < (kMcmCol-1) ) {
 348       if( (fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12) ) {
 349         fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1];
 350       } else {
 351         fHit2padSum[1][i] = -1;
 352       }
 353     } else {
 354       if ( fPadHits[i] >= sum12 ) {
 355         fHit2padSum[1][i] = fPadHits[i];
 356       } else {
 357         fHit2padSum[1][i] = -1;
 358       }
 359     }
 360     if (fHit2padSum[1][i] > 0) nsum2seed++;
 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
 393   nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
 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 //_____________________________________________________________________________
 446 void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir)
 447 {
 448   //
 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
 452   //
 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 //_____________________________________________________________________________
 496 void 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 //_____________________________________________________________________________
 511 Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
 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 //_____________________________________________________________________________
 531 void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
 532 {
 533   //
 534   // exponential filter
 535   //
 536
 537   Double_t sour[kMcmTBmax];
 538   Double_t dtarg[kMcmTBmax];
 539   Int_t    itarg[kMcmTBmax];
 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       }
 549       DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
 550       for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
 551         //fADC[iCol][iTime] = (Int_t)TMath::Max(0.0,dtarg[iTime]);
 552         fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
 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       }
 563       DeConvExpD(sour,itarg,kMcmTBmax,nexp);
 564       for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
 565         fADC[iCol][iTime] = itarg[iTime];
 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       }
 576       DeConvExpMI(sour,dtarg,kMcmTBmax);
 577       for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
 578         //fADC[iCol][iTime] = (Int_t)TMath::Max(0.0,dtarg[iTime]);
 579         fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
 580       }
 581     }
 582     break;
 583
 584   default:
 585
 586     printf("Invalid filter type %d ! \n",ftype);
 587     return;
 588
 589   }
 590
 591 }
 592
 593 //_____________________________________________________________________________
 594 void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp)
 595 {
 596   //
 597   // Exponential filter "analog"
 598   //
 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)
 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;
 611
 612   coefficients[0] = c1;
 613   coefficients[1] = c2;
 614
 615   Double_t dt = 0.100;
 616   rates[0] = TMath::Exp(-dt/(r1));
 617   rates[1] = TMath::Exp(-dt/(r2));
 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 //_____________________________________________________________________________
 640 void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp)
 641 {
 642   //
 643   // Exponential filter "digital"
 644   //
 645
 646   Int_t fAlphaL, fAlphaS, fLambdaL, fLambdaS, fTailPed;
 647   Int_t iAlphaL, iAlphaS, iLambdaL, iLambdaS;
 648
 649   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 650   // initialize (coefficient = alpha, rates = lambda)
 651
 652   fLambdaL = 0; fAlphaL = 0; fLambdaS = 0; fAlphaS = 0;
 653   iLambdaL = 0; iAlphaL = 0; iLambdaS = 0; iAlphaS = 0;
 654
 655   Double_t dt = 0.100;
 656
 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;
 662
 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
 667
 668   if (nexp == 1) {
 669     fAlphaL = (Int_t)(c1*2048.0);
 670     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
 671   }
 672   if (nexp == 2) {
 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
 677   }
 678
 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
 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
 711     h1 = ( rem1 + ((iAlphaL * result) >> 11) );
 712     if ( h1 > 0x0FFF ) {
 713       h1 = 0x0FFF;
 714     } else {
 715       h1 &= 0x0FFF;
 716     }
 717
 718     h2 = ( rem2 + ((iAlphaS * result) >> 11));
 719     if ( h2 > 0x0FFF ) {
 720       h2 = 0x0FFF;
 721     } else {
 722       h2 &= 0x0FFF;
 723     }
 724
 725     rem1 = (iLambdaL * h1 ) >> 11;
 726     rem2 = (iLambdaS * h2 ) >> 11;
 727
 728     if ( (rem1 + rem2) > 0x0FFF ) {
 729       correction = 0x0FFF;
 730     } else {
 731       correction = (rem1 + rem2) & 0x0FFF;
 732     }
 733   }
 734
 735 }
 736
 737 //_____________________________________________________________________________
 738 void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n)
 739 {
 740   //
 741   // Exponential filter (M. Ivanov)
 742   //
 743
 744    Double_t sig1[100], sig2[100], sig3[100];//, sig4[100];
 745    for (Int_t i = 0; i < n; i++) sig1[i] = source[i];
 746
 747    Float_t dt = 0.100;
 748
 749    //Float_t lambda0 = 9.8016*dt;  // short
 750    //Float_t lambda1 = 1.0778*dt;  // long
 751
 752    Float_t lambda0 = (1.0/fR2)*dt;
 753    Float_t lambda1 = (1.0/fR1)*dt;
 754
 755    TailMakerSpline(sig1,sig2,lambda0,n);
 756    TailCancelationMI(sig2,sig3,0.7,lambda1,n);
 757
 758    for (Int_t i = 0; i < n; i++) target[i] = sig3[i];
 759
 760 }
 761
 762 //______________________________________________________________________________
 763 void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
 764 {
 765   //
 766   // Special filter (M. Ivanov)
 767   //
 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 //______________________________________________________________________________
 828 void 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   //
 833
 834    Double_t l = TMath::Exp(-lambda*0.5);
 835    Double_t k = l*(1.-norm*lambda*0.5);
 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++){
 869      out[i]      = in[i]   + (k-l)*temp;
 870      temp        = in[i]   +  k*temp;
 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