TRD/AliTRDclusterizerV1.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /* $Id$ */
  18
  19 ///////////////////////////////////////////////////////////////////////////////
  20 //                                                                           //
  21 // TRD cluster finder for the slow simulator.
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 #include <TF1.h>
  26 #include <TTree.h>
  27 #include <TH1.h>
  28 #include <TFile.h>
  29
  30 #include "AliRun.h"
  31 #include "AliRunLoader.h"
  32 #include "AliLoader.h"
  33 #include "AliRawReader.h"
  34
  35 #include "AliTRDclusterizerV1.h"
  36 #include "AliTRDmatrix.h"
  37 #include "AliTRDgeometry.h"
  38 #include "AliTRDdataArrayF.h"
  39 #include "AliTRDdataArrayI.h"
  40 #include "AliTRDdigitsManager.h"
  41 #include "AliTRDpadPlane.h"
  42 #include "AliTRDrawData.h"
  43 #include "AliTRDcalibDB.h"
  44 #include "AliTRDSimParam.h"
  45 #include "AliTRDRecParam.h"
  46 #include "AliTRDCommonParam.h"
  47 #include "AliTRDcluster.h"
  48
  49 ClassImp(AliTRDclusterizerV1)
  50
  51 //_____________________________________________________________________________
  52 AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
  53 {
  54   //
  55   // AliTRDclusterizerV1 default constructor
  56   //
  57
  58   fDigitsManager = 0;
  59
  60 }
  61
  62 //_____________________________________________________________________________
  63 AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
  64                     :AliTRDclusterizer(name,title)
  65 {
  66   //
  67   // AliTRDclusterizerV1 default constructor
  68   //
  69
  70   fDigitsManager = new AliTRDdigitsManager();
  71   fDigitsManager->CreateArrays();
  72
  73 }
  74
  75 //_____________________________________________________________________________
  76 AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
  77 :AliTRDclusterizer(c)
  78 {
  79   //
  80   // AliTRDclusterizerV1 copy constructor
  81   //
  82
  83   ((AliTRDclusterizerV1 &) c).Copy(*this);
  84
  85 }
  86
  87 //_____________________________________________________________________________
  88 AliTRDclusterizerV1::~AliTRDclusterizerV1()
  89 {
  90   //
  91   // AliTRDclusterizerV1 destructor
  92   //
  93
  94   if (fDigitsManager) {
  95     delete fDigitsManager;
  96     fDigitsManager = NULL;
  97   }
  98
  99 }
 100
 101 //_____________________________________________________________________________
 102 AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
 103 {
 104   //
 105   // Assignment operator
 106   //
 107
 108   if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
 109   return *this;
 110
 111 }
 112
 113 //_____________________________________________________________________________
 114 void AliTRDclusterizerV1::Copy(TObject &c) const
 115 {
 116   //
 117   // Copy function
 118   //
 119
 120   ((AliTRDclusterizerV1 &) c).fDigitsManager = 0;
 121
 122   AliTRDclusterizer::Copy(c);
 123
 124 }
 125
 126 //_____________________________________________________________________________
 127 Bool_t AliTRDclusterizerV1::ReadDigits()
 128 {
 129   //
 130   // Reads the digits arrays from the input aliroot file
 131   //
 132
 133   if (!fRunLoader) {
 134     printf("<AliTRDclusterizerV1::ReadDigits> ");
 135     printf("No input file open\n");
 136     return kFALSE;
 137   }
 138   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 139   if (!loader->TreeD()) loader->LoadDigits();
 140
 141   // Read in the digit arrays
 142   return (fDigitsManager->ReadDigits(loader->TreeD()));
 143
 144 }
 145
 146 //_____________________________________________________________________________
 147 Bool_t AliTRDclusterizerV1::ReadDigits(AliRawReader* rawReader)
 148 {
 149   //
 150   // Reads the digits arrays from the ddl file
 151   //
 152
 153   AliTRDrawData *raw = new AliTRDrawData();
 154   raw->SetDebug(1);
 155
 156   fDigitsManager = raw->Raw2Digits(rawReader);
 157
 158   return kTRUE;
 159
 160 }
 161
 162 //_____________________________________________________________________________
 163 Bool_t AliTRDclusterizerV1::MakeClusters()
 164 {
 165   //
 166   // Generates the cluster.
 167   //
 168
 169   Int_t row, col, time;
 170
 171   /*
 172   if (fTRD->IsVersion() != 1) {
 173     printf("<AliTRDclusterizerV1::MakeCluster> ");
 174     printf("TRD must be version 1 (slow simulator).\n");
 175     return kFALSE;
 176   }
 177   */
 178
 179   // Get the geometry
 180   AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader);
 181   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
 182   if (!calibration)
 183   {
 184     printf("<AliTRDclusterizerV1::MakeCluster> ");
 185     printf("ERROR getting instance of AliTRDcalibDB");
 186     return kFALSE;
 187   }
 188
 189   AliTRDSimParam* simParam = AliTRDSimParam::Instance();
 190   if (!simParam)
 191   {
 192     printf("<AliTRDclusterizerV1::MakeCluster> ");
 193     printf("ERROR getting instance of AliTRDSimParam");
 194     return kFALSE;
 195   }
 196
 197   AliTRDRecParam* recParam = AliTRDRecParam::Instance();
 198   if (!recParam)
 199   {
 200     printf("<AliTRDclusterizerV1::MakeCluster> ");
 201     printf("ERROR getting instance of AliTRDRecParam");
 202     return kFALSE;
 203   }
 204
 205   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
 206   if (!commonParam)
 207   {
 208     printf("<AliTRDclusterizerV1::MakeDigits> ");
 209     printf("Could not get common params\n");
 210     return kFALSE;
 211   }
 212
 213   Float_t ADCthreshold = simParam->GetADCthreshold();
 214
 215   if (fVerbose > 0) {
 216     //printf("<AliTRDclusterizerV1::MakeCluster> ");
 217     //printf("OmegaTau = %f \n",omegaTau);
 218     printf("<AliTRDclusterizerV1::MakeCluster> ");
 219     printf("Start creating clusters.\n");
 220   }
 221
 222   AliTRDdataArrayI *digitsIn;
 223   AliTRDdataArrayI *track0;
 224   AliTRDdataArrayI *track1;
 225   AliTRDdataArrayI *track2;
 226
 227   // Threshold value for the maximum
 228   Float_t maxThresh = recParam->GetClusMaxThresh();
 229   // Threshold value for the digit signal
 230   Float_t sigThresh = recParam->GetClusSigThresh();
 231   // Iteration limit for unfolding procedure
 232   const Float_t kEpsilon = 0.01;
 233
 234   const Int_t   kNclus   = 3;
 235   const Int_t   kNsig    = 5;
 236   const Int_t   kNtrack  = 3 * kNclus;
 237
 238   Int_t    iType         = 0;
 239   Int_t    iUnfold       = 0;
 240   Double_t ratioLeft     = 1.0;
 241   Double_t ratioRight    = 1.0;
 242
 243   //
 244   Double_t padSignal[kNsig];
 245   Double_t clusterSignal[kNclus];
 246   Double_t clusterPads[kNclus];
 247   Int_t    clusterTracks[kNtrack];
 248
 249   Int_t    chamBeg = 0;
 250   Int_t    chamEnd = AliTRDgeometry::Ncham();
 251   Int_t    planBeg = 0;
 252   Int_t    planEnd = AliTRDgeometry::Nplan();
 253   Int_t    sectBeg = 0;
 254   Int_t    sectEnd = AliTRDgeometry::Nsect();
 255
 256   Int_t    nTimeTotal  = calibration->GetNumberOfTimeBins();
 257
 258   if (fVerbose > 0) {
 259     printf("<AliTRDclusterizerV1::MakeCluster> ");
 260     printf("Number of Time Bins = %d.\n",nTimeTotal);
 261   }
 262
 263   // Start clustering in every chamber
 264   for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
 265     for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
 266       for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
 267
 268         Int_t idet    = geo->GetDetector(iplan,icham,isect);
 269
 270         Int_t nRowMax = commonParam->GetRowMax(iplan,icham,isect);
 271         Int_t nColMax = commonParam->GetColMax(iplan);
 272
 273         Int_t nClusters      = 0;
 274         Int_t nClusters2pad  = 0;
 275         Int_t nClusters3pad  = 0;
 276         Int_t nClusters4pad  = 0;
 277         Int_t nClusters5pad  = 0;
 278         Int_t nClustersLarge = 0;
 279
 280         if (fVerbose > 0) {
 281           printf("<AliTRDclusterizerV1::MakeCluster> ");
 282           printf("Analyzing chamber %d, plane %d, sector %d.\n"
 283                 ,icham,iplan,isect);
 284         }
 285
 286         AliTRDpadPlane *padPlane = commonParam->GetPadPlane(iplan,icham);
 287
 288         // Get the digits
 289         digitsIn = fDigitsManager->GetDigits(idet);
 290         digitsIn->Expand();
 291         AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow(), digitsIn->GetNcol(), digitsIn->GetNtime());
 292
 293         Transform(digitsIn, digitsOut, idet, nRowMax, nColMax, nTimeTotal, ADCthreshold);
 294
 295         track0 = fDigitsManager->GetDictionary(idet,0);
 296         track0->Expand();
 297         track1 = fDigitsManager->GetDictionary(idet,1);
 298         track1->Expand();
 299         track2 = fDigitsManager->GetDictionary(idet,2);
 300         track2->Expand();
 301
 302         // Loop through the chamber and find the maxima
 303         for ( row = 0;  row <  nRowMax;    row++) {
 304           for ( col = 2;  col <  nColMax;    col++) {
 305             //for ( col = 4;  col <  nColMax-2;    col++) {
 306             for (time = 0; time < nTimeTotal; time++) {
 307
 308               Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col  ,time));
 309               Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
 310               Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
 311
 312 //            // Lonok for the maximum
 313 //               if (signalM >= maxThresh) {
 314 //                 if (((signalL >= sigThresh) &&
 315 //                      (signalL <  signalM))  ||
 316 //                     ((signalR >= sigThresh) &&
 317 //                      (signalR <  signalM))) {
 318 //                   // Maximum found, mark the position by a negative signal
 319 //                   digitsOut->SetDataUnchecked(row,col-1,time,-signalM);
 320 //              }
 321 //            }
 322               // Look for the maximum
 323               if (signalM >= maxThresh) {
 324                 if ( (TMath::Abs(signalL)<=signalM) && (TMath::Abs(signalR)<=signalM) &&
 325                      (TMath::Abs(signalL)+TMath::Abs(signalR))>sigThresh ) {
 326                   // Maximum found, mark the position by a negative signal
 327                   digitsOut->SetDataUnchecked(row,col-1,time,-signalM);
 328                 }
 329               }
 330             }
 331           }
 332         }
 333
 334         // Now check the maxima and calculate the cluster position
 335         for ( row = 0;  row <  nRowMax  ;  row++) {
 336           for (time = 0; time < nTimeTotal; time++) {
 337             for ( col = 1;  col <  nColMax-1;  col++) {
 338
 339               // Maximum found ?
 340               if (digitsOut->GetDataUnchecked(row,col,time) < 0) {
 341
 342                 Int_t iPad;
 343                 for (iPad = 0; iPad < kNclus; iPad++) {
 344                   Int_t iPadCol = col - 1 + iPad;
 345                   clusterSignal[iPad]     = TMath::Abs(digitsOut->GetDataUnchecked(row
 346                                                                                ,iPadCol
 347                                                                                ,time));
 348                   clusterTracks[3*iPad  ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
 349                   clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
 350                   clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
 351                 }
 352
 353                 // Count the number of pads in the cluster
 354                 Int_t nPadCount = 0;
 355                 Int_t ii        = 0;
 356                 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii  ,time))
 357                                                                   >= sigThresh) {
 358                   nPadCount++;
 359                   ii++;
 360                   if (col-ii   <        0) break;
 361                 }
 362                 ii = 0;
 363                 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time))
 364                                                                   >= sigThresh) {
 365                   nPadCount++;
 366                   ii++;
 367                   if (col+ii+1 >= nColMax) break;
 368                 }
 369
 370                 nClusters++;
 371                 switch (nPadCount) {
 372                 case 2:
 373                   iType = 0;
 374                   nClusters2pad++;
 375                   break;
 376                 case 3:
 377                   iType = 1;
 378                   nClusters3pad++;
 379                   break;
 380                 case 4:
 381                   iType = 2;
 382                   nClusters4pad++;
 383                   break;
 384                 case 5:
 385                   iType = 3;
 386                   nClusters5pad++;
 387                   break;
 388                 default:
 389                   iType = 4;
 390                   nClustersLarge++;
 391                   break;
 392                 };
 393
 394                  // Look for 5 pad cluster with minimum in the middle
 395                 Bool_t fivePadCluster = kFALSE;
 396                 if (col < nColMax-3) {
 397                   if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
 398                     fivePadCluster = kTRUE;
 399                   }
 400                   if ((fivePadCluster) && (col < nColMax-5)) {
 401                     if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) {
 402                       fivePadCluster = kFALSE;
 403                     }
 404                   }
 405                   if ((fivePadCluster) && (col >         1)) {
 406                     if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) {
 407                       fivePadCluster = kFALSE;
 408                     }
 409                   }
 410                 }
 411
 412                 // 5 pad cluster
 413                 // Modify the signal of the overlapping pad for the left part
 414                 // of the cluster which remains from a previous unfolding
 415                 if (iUnfold) {
 416                   clusterSignal[0] *= ratioLeft;
 417                   iType   = 5;
 418                   iUnfold = 0;
 419                 }
 420
 421                 // Unfold the 5 pad cluster
 422                 if (fivePadCluster) {
 423                   for (iPad = 0; iPad < kNsig; iPad++) {
 424                     padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
 425                                                                              ,col-1+iPad
 426                                                                              ,time));
 427                   }
 428                   // Unfold the two maxima and set the signal on
 429                   // the overlapping pad to the ratio
 430                   ratioRight        = Unfold(kEpsilon,iplan,padSignal);
 431                   ratioLeft         = 1.0 - ratioRight;
 432                   clusterSignal[2] *= ratioRight;
 433                   iType   = 5;
 434                   iUnfold = 1;
 435                 }
 436
 437                 Double_t clusterCharge = clusterSignal[0]
 438                                        + clusterSignal[1]
 439                                        + clusterSignal[2];
 440
 441                 // The position of the cluster
 442                 clusterPads[0] = row + 0.5;
 443                 // Take the shift of the additional time bins into account
 444                 clusterPads[2] = time + 0.5;
 445
 446
 447                 if (recParam->LUTOn()) {
 448                   // Calculate the position of the cluster by using the
 449                   // lookup table method
 450                   clusterPads[1] = recParam->LUTposition(iplan,clusterSignal[0]
 451                                                          ,clusterSignal[1]
 452                                                          ,clusterSignal[2]);
 453                 }
 454                 else {
 455                   // Calculate the position of the cluster by using the
 456                   // center of gravity method
 457                   for (Int_t i=0;i<5;i++) padSignal[i]=0;
 458                   padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));   // central  pad
 459                   padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // left     pad
 460                   padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // right    pad
 461                   if (col>2 &&TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)<padSignal[1])){
 462                     padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
 463                   }
 464                   if (col<nColMax-3 &&TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)<padSignal[3])){
 465                     padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
 466                   }
 467                   clusterPads[1] =  GetCOG(padSignal);
 468
 469                 }
 470
 471                 Double_t q0 = clusterSignal[0];
 472                 Double_t q1 = clusterSignal[1];
 473                 Double_t q2 = clusterSignal[2];
 474                 Double_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
 475                                           (clusterCharge*clusterCharge);
 476
 477
 478
 479                 // Calculate the position and the error
 480
 481                 // correct for t0
 482                 Int_t clusterTimeBin = TMath::Nint(time - calibration->GetT0(idet, col, row));
 483
 484                 Double_t colSize = padPlane->GetColSize(col);
 485                 Double_t rowSize = padPlane->GetRowSize(row);
 486                 Double_t clusterPos[3];
 487                 clusterPos[0] = padPlane->GetColPos(col) - (clusterPads[1]+0.5)*colSize;  // MI change
 488                 clusterPos[1] = padPlane->GetRowPos(row) - 0.5*rowSize; //MI change
 489                 clusterPos[2] = CalcXposFromTimebin(clusterPads[2], idet, col, row);
 490                 Double_t clusterSig[2];
 491                 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
 492                 clusterSig[1] = rowSize * rowSize / 12.;
 493
 494
 495                 // Add the cluster to the output array
 496                 AliTRDcluster * cluster = AddCluster(clusterPos
 497                           ,clusterTimeBin
 498                           ,idet
 499                           ,clusterCharge
 500                           ,clusterTracks
 501                           ,clusterSig
 502                           ,iType,clusterPads[1]);
 503                 //
 504                 //
 505                 Short_t signals[7]={0,0,0,0,0,0,0};
 506                 for (Int_t jPad = col-3;jPad<=col+3;jPad++){
 507                   if (jPad<0 ||jPad>=nColMax-1) continue;
 508                   signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
 509                 }
 510                 cluster->SetSignals(signals);
 511               }
 512             }
 513           }
 514         }
 515
 516         delete digitsOut;
 517
 518         // Compress the arrays
 519         track0->Compress(1,0);
 520         track1->Compress(1,0);
 521         track2->Compress(1,0);
 522
 523         // Write the cluster and reset the array
 524         WriteClusters(idet);
 525         ResetRecPoints();
 526       }
 527     }
 528   }
 529
 530   if (fVerbose > 0) {
 531     printf("<AliTRDclusterizerV1::MakeCluster> ");
 532     printf("Done.\n");
 533   }
 534
 535   //delete digitsIn;
 536
 537   return kTRUE;
 538
 539 }
 540
 541 //_____________________________________________________________________________
 542 Double_t AliTRDclusterizerV1::GetCOG(Double_t signal[5])
 543 {
 544   //
 545   // get COG position
 546   // used for clusters with more than 3 pads - where LUT not applicable
 547   Double_t sum = signal[0]+signal[1]+signal[2]+signal[3]+signal[4];
 548   Double_t res = (0.0*(-signal[0]+signal[4])+(-signal[1]+signal[3]))/sum;
 549   return res;
 550 }
 551
 552 //_____________________________________________________________________________
 553 Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSignal)
 554 {
 555   //
 556   // Method to unfold neighbouring maxima.
 557   // The charge ratio on the overlapping pad is calculated
 558   // until there is no more change within the range given by eps.
 559   // The resulting ratio is then returned to the calling method.
 560   //
 561
 562   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
 563   if (!calibration)
 564   {
 565     printf("<AliTRDclusterizerMI::Unfold> ");
 566     printf("ERROR getting instance of AliTRDcalibDB");
 567     return kFALSE;
 568   }
 569
 570   Int_t   irc                = 0;
 571   Int_t   itStep             = 0;      // Count iteration steps
 572
 573   Double_t ratio             = 0.5;    // Start value for ratio
 574   Double_t prevRatio         = 0;      // Store previous ratio
 575
 576   Double_t newLeftSignal[3]  = {0};    // Array to store left cluster signal
 577   Double_t newRightSignal[3] = {0};    // Array to store right cluster signal
 578   Double_t newSignal[3]      = {0};
 579
 580   // Start the iteration
 581   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
 582
 583     itStep++;
 584     prevRatio = ratio;
 585
 586     // Cluster position according to charge ratio
 587     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
 588                       / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
 589     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
 590                       / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
 591
 592     // Set cluster charge ratio
 593     irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
 594     Double_t ampLeft  = padSignal[1] / newSignal[1];
 595     irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
 596     Double_t ampRight = padSignal[3] / newSignal[1];
 597
 598     // Apply pad response to parameters
 599     irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
 600     irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
 601
 602     // Calculate new overlapping ratio
 603     ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] /
 604                           (newLeftSignal[2] + newRightSignal[0]));
 605
 606   }
 607
 608   return ratio;
 609
 610 }
 611
 612 //_____________________________________________________________________________
 613 void AliTRDclusterizerV1::Transform(AliTRDdataArrayI* digitsIn,
 614                                     AliTRDdataArrayF* digitsOut,
 615                                     Int_t idet, Int_t nRowMax,
 616                                     Int_t nColMax, Int_t nTimeTotal,
 617                                     Float_t ADCthreshold)
 618 {
 619
 620   //
 621   // Apply gain factor
 622   // Apply tail cancellation: Transform digitsIn to digitsOut
 623   //
 624
 625
 626   AliTRDRecParam* recParam = AliTRDRecParam::Instance();
 627   if (!recParam)
 628   {
 629     printf("<AliTRDclusterizerV1::Transform> ");
 630     printf("ERROR getting instance of AliTRDRecParam");
 631     return;
 632   }
 633   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
 634
 635   Double_t *inADC  = new Double_t[nTimeTotal];  // adc data before tail cancellation
 636   Double_t *outADC = new Double_t[nTimeTotal];  // adc data after tail cancellation
 637
 638   if (fVerbose > 0) {
 639     printf("<AliTRDclusterizerV1::Transform> ");
 640     printf("Tail cancellation (nExp = %d) for detector %d.\n",
 641            recParam->GetTCnexp(),idet);
 642   }
 643
 644   for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
 645     for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
 646       for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
 647         //
 648         // add gain
 649         //
 650         Double_t gain = calibration->GetGainFactor(idet, iCol, iRow);
 651         if (gain==0) {
 652           AliError("Not a valid gain\n");
 653         }
 654         inADC[iTime]  = digitsIn->GetDataUnchecked(iRow, iCol, iTime);
 655
 656         inADC[iTime] /= gain;
 657         outADC[iTime] = inADC[iTime];
 658       }
 659
 660       // Apply the tail cancelation via the digital filter
 661       if (recParam->TCOn())
 662       {
 663         DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp());
 664       }
 665
 666       for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
 667         // Store the amplitude of the digit if above threshold
 668         if (outADC[iTime] > ADCthreshold) {
 669           if (fVerbose > 1)
 670           {
 671             printf("  iRow = %d, iCol = %d, iTime = %d, adc = %f\n"
 672                    ,iRow,iCol,iTime,outADC[iTime]);
 673           }
 674           digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
 675         }
 676
 677       }
 678
 679     }
 680
 681   }
 682
 683   delete [] inADC;
 684   delete [] outADC;
 685
 686   return;
 687
 688 }
 689
 690
 691 //_____________________________________________________________________________
 692 void AliTRDclusterizerV1::DeConvExp(Double_t *source, Double_t *target,
 693                                     Int_t n, Int_t nexp)
 694 {
 695   //
 696   // Tail Cancellation by Deconvolution for PASA v4 TRF
 697   //
 698
 699   Double_t rates[2];
 700   Double_t coefficients[2];
 701
 702   // initialize (coefficient = alpha, rates = lambda)
 703
 704   Double_t R1 = 1.0;
 705   Double_t R2 = 1.0;
 706   Double_t C1 = 0.5;
 707   Double_t C2 = 0.5;
 708
 709   if (nexp == 1) {   // 1 Exponentials
 710     R1 = 1.156;
 711     R2 = 0.130;
 712     C1 = 0.066;
 713     C2 = 0.000;
 714   }
 715   if (nexp == 2) {   // 2 Exponentials
 716     R1 = 1.156;
 717     R2 = 0.130;
 718     C1 = 0.114;
 719     C2 = 0.624;
 720   }
 721
 722   coefficients[0] = C1;
 723   coefficients[1] = C2;
 724
 725   Double_t Dt = 0.100;
 726
 727   rates[0] = TMath::Exp(-Dt/(R1));
 728   rates[1] = TMath::Exp(-Dt/(R2));
 729
 730   Int_t i, k;
 731   Double_t reminder[2];
 732   Double_t correction, result;
 733
 734   /* attention: computation order is important */
 735   correction=0.0;
 736
 737   for ( k=0; k<nexp; k++ ) reminder[k]=0.0;
 738
 739   for ( i=0; i<n; i++ ) {
 740     result = ( source[i] - correction );    // no rescaling
 741     target[i] = result;
 742
 743     for ( k=0; k<nexp; k++ ) reminder[k] = rates[k] *
 744                                ( reminder[k] + coefficients[k] * result);
 745     correction=0.0;
 746     for ( k=0; k<nexp; k++ ) correction += reminder[k];
 747   }
 748
 749 }