TRD/AliTRDclusterizerV2.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                                                        //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 #include <TF1.h>
  26 #include <TTree.h>
  27 #include <TH1.h>
  28 #include <TFile.h>
  29
  30 #include "AliRunLoader.h"
  31 #include "AliLoader.h"
  32 #include "AliRawReader.h"
  33 #include "AliLog.h"
  34 #include "AliAlignObj.h"
  35
  36 #include "AliTRDclusterizerV2.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 "AliTRDcluster.h"
  47
  48 #include "Cal/AliTRDCalROC.h"
  49 #include "Cal/AliTRDCalDet.h"
  50
  51 #include "AliTRDSignalIndex.h"
  52 #include "AliTRDRawStream.h"
  53 #include "AliTRDRawStreamV2.h"
  54
  55 #include "AliTRDfeeParam.h"
  56
  57 ClassImp(AliTRDclusterizerV2)
  58
  59 //_____________________________________________________________________________
  60 AliTRDclusterizerV2::AliTRDclusterizerV2()
  61   :AliTRDclusterizer()
  62   ,fDigitsManager(NULL)
  63   ,fGeometry(NULL)
  64   ,fAddLabels(kTRUE)
  65   ,fRawVersion(2)
  66   ,fIndexesOut(NULL)
  67   ,fIndexesMaxima(NULL)
  68 {
  69   //
  70   // AliTRDclusterizerV2 default constructor
  71   //
  72
  73   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
  74 }
  75
  76 //_____________________________________________________________________________
  77 AliTRDclusterizerV2::AliTRDclusterizerV2(const Text_t *name, const Text_t *title)
  78   :AliTRDclusterizer(name,title)
  79   ,fDigitsManager(new AliTRDdigitsManager())
  80   ,fGeometry(NULL)
  81   ,fAddLabels(kTRUE)
  82   ,fRawVersion(2)
  83   ,fIndexesOut(NULL)
  84   ,fIndexesMaxima(NULL)
  85 {
  86   //
  87   // AliTRDclusterizerV2 constructor
  88   //
  89
  90   fDigitsManager->CreateArrays();
  91   fGeometry = new AliTRDgeometry;
  92
  93   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
  94 }
  95
  96 //_____________________________________________________________________________
  97 AliTRDclusterizerV2::AliTRDclusterizerV2(const AliTRDclusterizerV2 &c)
  98   :AliTRDclusterizer(c)
  99   ,fDigitsManager(NULL)
 100   ,fGeometry(NULL)
 101   ,fAddLabels(kTRUE)
 102   ,fRawVersion(2)
 103   ,fIndexesOut(NULL)
 104   ,fIndexesMaxima(NULL)
 105 {
 106   //
 107   // AliTRDclusterizerV2 copy constructor
 108   //
 109
 110 }
 111
 112 //_____________________________________________________________________________
 113 AliTRDclusterizerV2::~AliTRDclusterizerV2()
 114 {
 115   //
 116   // AliTRDclusterizerV2 destructor
 117   //
 118
 119   if (fDigitsManager) {
 120     delete fDigitsManager;
 121     fDigitsManager = NULL;
 122   }
 123
 124   if (fGeometry)
 125     {
 126       delete fGeometry;
 127       fGeometry = NULL;
 128     }
 129
 130   if (fIndexesOut)
 131     {
 132       delete fIndexesOut;
 133       fIndexesOut = NULL;
 134     }
 135
 136   if (fIndexesMaxima)
 137     {
 138       delete fIndexesMaxima;
 139       fIndexesMaxima = NULL;
 140     }
 141 }
 142
 143 //_____________________________________________________________________________
 144 AliTRDclusterizerV2 &AliTRDclusterizerV2::operator=(const AliTRDclusterizerV2 &c)
 145 {
 146   //
 147   // Assignment operator
 148   //
 149
 150   if (this != &c) ((AliTRDclusterizerV2 &) c).Copy(*this);
 151   return *this;
 152
 153 }
 154
 155 //_____________________________________________________________________________
 156 void AliTRDclusterizerV2::Copy(TObject &c) const
 157 {
 158   //
 159   // Copy function
 160   //
 161
 162   ((AliTRDclusterizerV2 &) c).fDigitsManager = NULL;
 163   ((AliTRDclusterizerV2 &) c).fGeometry = NULL;
 164   ((AliTRDclusterizerV2 &) c).fAddLabels = fAddLabels;
 165   ((AliTRDclusterizerV2 &) c).fRawVersion = fRawVersion;
 166   ((AliTRDclusterizerV2 &) c).fIndexesOut = NULL;
 167   ((AliTRDclusterizerV2 &) c).fIndexesMaxima = NULL;
 168
 169   AliTRDclusterizer::Copy(c);
 170
 171 }
 172
 173 //_____________________________________________________________________________
 174 void AliTRDclusterizerV2::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
 175 {
 176   //
 177   // Reset the helper indexes
 178   //
 179   if (fIndexesOut)
 180     {
 181       // carefull here - we assume that only row number may change - most probable
 182       if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
 183         fIndexesOut->ResetContent();
 184       else
 185         fIndexesOut->ResetContentConditional(indexesIn->GetNrow(), indexesIn->GetNcol(), indexesIn->GetNtime());
 186     }
 187   else
 188     {
 189       fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow(), indexesIn->GetNcol(), indexesIn->GetNtime());
 190     }
 191
 192   if (fIndexesMaxima)
 193     {
 194       // carefull here - we assume that only row number may change - most probable
 195       if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
 196         fIndexesMaxima->ResetContent();
 197       else
 198         fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow(), indexesIn->GetNcol(), indexesIn->GetNtime());
 199     }
 200   else
 201     {
 202       fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow(), indexesIn->GetNcol(), indexesIn->GetNtime());
 203     }
 204 }
 205
 206 //_____________________________________________________________________________
 207 Bool_t AliTRDclusterizerV2::ReadDigits()
 208 {
 209   //
 210   // Reads the digits arrays from the input aliroot file
 211   //
 212
 213   if (!fRunLoader) {
 214     AliError("No run loader available");
 215     return kFALSE;
 216   }
 217
 218   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 219   if (!loader->TreeD()) {
 220     loader->LoadDigits();
 221   }
 222
 223   // Read in the digit arrays
 224   return (fDigitsManager->ReadDigits(loader->TreeD()));
 225
 226 }
 227
 228 //_____________________________________________________________________________
 229 Bool_t AliTRDclusterizerV2::ReadDigits(TTree *digitsTree)
 230 {
 231   //
 232   // Reads the digits arrays from the input tree
 233   //
 234
 235   // Read in the digit arrays
 236   return (fDigitsManager->ReadDigits(digitsTree));
 237
 238 }
 239
 240 //_____________________________________________________________________________
 241 Bool_t AliTRDclusterizerV2::ReadDigits(AliRawReader *rawReader)
 242 {
 243   //
 244   // Reads the digits arrays from the ddl file
 245   //
 246
 247   AliTRDrawData raw;
 248   fDigitsManager = raw.Raw2Digits(rawReader);
 249
 250   return kTRUE;
 251
 252 }
 253
 254 //_____________________________________________________________________________
 255 Bool_t AliTRDclusterizerV2::MakeClusters()
 256 {
 257   //
 258   // Creates clusters from digits
 259   //
 260
 261   //propagate info from the digits manager
 262   if (fAddLabels == kTRUE)
 263     fAddLabels = fDigitsManager->UsesDictionaries();
 264
 265   Bool_t fReturn = kTRUE;
 266   for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++)
 267     {
 268       AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(i);
 269       // This is to take care of switched off super modules
 270       if (digitsIn->GetNtime() == 0) {
 271         continue;
 272       }
 273       //AliInfo(Form("digitsIn->Expand() 0x%x", digitsIn));
 274       digitsIn->Expand();
 275       AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
 276       if (indexes->IsAllocated() == kFALSE)
 277         {
 278           fDigitsManager->BuildIndexes(i);
 279         }
 280
 281       Bool_t fR = kFALSE;
 282       if (indexes->HasEntry())
 283         {
 284           if (fAddLabels)
 285             {
 286               for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++)
 287                 {
 288                   AliTRDdataArrayI *tracksIn = 0;
 289                   tracksIn = fDigitsManager->GetDictionary(i,iDict);
 290                   tracksIn->Expand();
 291                 }
 292             }
 293           fR = MakeClusters(i);
 294           fReturn = fR && fReturn;
 295         }
 296
 297       if (fR == kFALSE)
 298         {
 299           WriteClusters(i);
 300           ResetRecPoints();
 301         }
 302       //digitsIn->Compress(1,0);
 303       // no compress just remove
 304       fDigitsManager->RemoveDigits(i);
 305       fDigitsManager->RemoveDictionaries(i);
 306       fDigitsManager->ClearIndexes(i);
 307     }
 308
 309   return fReturn;
 310 }
 311
 312
 313 //_____________________________________________________________________________
 314 Bool_t AliTRDclusterizerV2::Raw2Clusters(AliRawReader *rawReader)
 315 {
 316   //
 317   // Creates clusters from raw data
 318   //
 319
 320   //AliDebug(2, "Read all");
 321
 322   AliTRDdataArrayI *digits = 0;
 323   AliTRDdataArrayI *track0 = 0;
 324   AliTRDdataArrayI *track1 = 0;
 325   AliTRDdataArrayI *track2 = 0;
 326
 327   AliTRDSignalIndex *indexes = 0;
 328   // Create the digits manager
 329   if (!fDigitsManager)
 330     {
 331       fDigitsManager = new AliTRDdigitsManager();
 332       fDigitsManager->CreateArrays();
 333     }
 334
 335   //AliTRDRawStream input(rawReader);
 336   AliTRDRawStreamV2 input(rawReader);
 337   input.SetRawVersion( fRawVersion );
 338   input.Init();
 339
 340   AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
 341
 342   // Loop through the digits
 343   Int_t lastdet = -1;
 344   Int_t det    = 0;
 345   Int_t it = 0;
 346   while (input.Next())
 347     {
 348
 349       det    = input.GetDet();
 350
 351       if (det != lastdet)
 352         {
 353           if (lastdet != -1)
 354             {
 355               digits = fDigitsManager->GetDigits(lastdet);
 356               Bool_t iclusterBranch = kFALSE;
 357               if (indexes->HasEntry())
 358                 iclusterBranch = MakeClusters(lastdet);
 359               if (iclusterBranch == kFALSE)
 360                 {
 361                   WriteClusters(lastdet);
 362                   ResetRecPoints();
 363                 }
 364             }
 365
 366           if (digits)
 367             {
 368               fDigitsManager->RemoveDigits(lastdet);
 369               fDigitsManager->RemoveDictionaries(lastdet);
 370               fDigitsManager->ClearIndexes(lastdet);
 371             }
 372
 373           lastdet = det;
 374
 375           // Add a container for the digits of this detector
 376           digits = fDigitsManager->GetDigits(det);
 377           track0 = fDigitsManager->GetDictionary(det,0);
 378           track1 = fDigitsManager->GetDictionary(det,1);
 379           track2 = fDigitsManager->GetDictionary(det,2);
 380
 381           // Allocate memory space for the digits buffer
 382           if (digits->GetNtime() == 0)
 383             {
 384               //AliDebug(5, Form("Alloc digits for det %d", det));
 385               digits->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
 386               track0->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
 387               track1->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
 388               track2->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
 389             }
 390
 391           indexes = fDigitsManager->GetIndexes(det);
 392           indexes->SetSM(input.GetSM());
 393           indexes->SetStack(input.GetStack());
 394           indexes->SetLayer(input.GetLayer());
 395           indexes->SetDetNumber(det);
 396           if (indexes->IsAllocated() == kFALSE)
 397             indexes->Allocate(input.GetMaxRow(), input.GetMaxCol(), input.GetNumberOfTimeBins());
 398         }
 399
 400       for (it = 0; it < 3; it++)
 401         {
 402           if ( input.GetTimeBin() + it < input.GetNumberOfTimeBins() )
 403             {
 404               if (input.GetSignals()[it] > 0)
 405                 {
 406                   digits->SetDataUnchecked(input.GetRow(), input.GetCol(),
 407                                            input.GetTimeBin() + it, input.GetSignals()[it]);
 408
 409                   indexes->AddIndexTBin(input.GetRow(), input.GetCol(),
 410                                         input.GetTimeBin() + it);
 411                   track0->SetDataUnchecked(input.GetRow(), input.GetCol(),
 412                                            input.GetTimeBin() + it, 0);
 413                   track1->SetDataUnchecked(input.GetRow(), input.GetCol(),
 414                                            input.GetTimeBin() + it, 0);
 415                   track2->SetDataUnchecked(input.GetRow(), input.GetCol(),
 416                                            input.GetTimeBin() + it, 0);
 417                 }
 418             }
 419         }
 420   }
 421
 422   if (lastdet != -1)
 423     {
 424       Bool_t iclusterBranch = kFALSE;
 425       if (indexes->HasEntry())
 426         iclusterBranch = MakeClusters(lastdet);
 427       if (iclusterBranch == kFALSE)
 428         {
 429           WriteClusters(lastdet);
 430           ResetRecPoints();
 431         }
 432       //MakeClusters(lastdet);
 433       if (digits)
 434         {
 435           fDigitsManager->RemoveDigits(lastdet);
 436           fDigitsManager->RemoveDictionaries(lastdet);
 437           fDigitsManager->ClearIndexes(lastdet);
 438         }
 439     }
 440
 441   delete fDigitsManager;
 442   fDigitsManager = NULL;
 443   return kTRUE;
 444 }
 445
 446 //_____________________________________________________________________________
 447 Bool_t AliTRDclusterizerV2::Raw2ClustersChamber(AliRawReader *rawReader)
 448 {
 449   //
 450   // Creates clusters from raw data
 451   //
 452
 453   //AliDebug(1, "Raw2ClustersChamber");
 454
 455   // Create the digits manager
 456   // Create the digits manager
 457   if (!fDigitsManager)
 458     {
 459       fDigitsManager = new AliTRDdigitsManager();
 460       fDigitsManager->CreateArrays();
 461     }
 462
 463   fDigitsManager->SetUseDictionaries(fAddLabels);
 464
 465   //AliTRDRawStream input(rawReader);
 466   AliTRDRawStreamV2 input(rawReader);
 467   input.SetRawVersion( fRawVersion );
 468   input.Init();
 469
 470   AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
 471
 472   Int_t det    = 0;
 473   while ((det = input.NextChamber(fDigitsManager)) >= 0)
 474     {
 475       Bool_t iclusterBranch = kFALSE;
 476       if (fDigitsManager->GetIndexes(det)->HasEntry())
 477         iclusterBranch = MakeClusters(det);
 478       if (iclusterBranch == kFALSE)
 479         {
 480           WriteClusters(det);
 481           ResetRecPoints();
 482         }
 483       fDigitsManager->RemoveDigits(det);
 484       fDigitsManager->RemoveDictionaries(det);
 485       fDigitsManager->ClearIndexes(det);
 486     }
 487
 488   delete fDigitsManager;
 489   fDigitsManager = NULL;
 490   return kTRUE;
 491 }
 492
 493 //_____________________________________________________________________________
 494 Bool_t AliTRDclusterizerV2::MakeClusters(Int_t det)
 495 {
 496   //
 497   // Generates the cluster.
 498   //
 499
 500   // Get the digits
 501   //   digits should be expanded beforehand!
 502   //   digitsIn->Expand();
 503   AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(det);
 504   // This is to take care of switched off super modules
 505   if (digitsIn->GetNtime() == 0)
 506     {
 507       //AliDebug(5, Form("digitsIn->GetNtime() == 0 [%d]", det));
 508       return kFALSE;
 509     }
 510
 511   AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
 512   if (indexesIn->IsAllocated() == kFALSE)
 513     {
 514       AliError("Indexes do not exist!");
 515       return kFALSE;
 516     }
 517
 518   AliTRDcalibDB  *calibration    = AliTRDcalibDB::Instance();
 519   if (!calibration) {
 520     AliFatal("No AliTRDcalibDB instance available\n");
 521     return kFALSE;
 522   }
 523
 524   AliTRDSimParam *simParam       = AliTRDSimParam::Instance();
 525   if (!simParam) {
 526     AliError("No AliTRDSimParam instance available\n");
 527     return kFALSE;
 528   }
 529
 530   AliTRDRecParam *recParam       = AliTRDRecParam::Instance();
 531   if (!recParam) {
 532     AliError("No AliTRDRecParam instance available\n");
 533     return kFALSE;
 534   }
 535
 536   // ADC thresholds
 537   //  Float_t ADCthreshold   = simParam->GetADCthreshold();
 538   Float_t ADCthreshold   = 0; // There is no ADC threshold anymore, and simParam should not be used ni clusterizer. KO
 539
 540   // Threshold value for the maximum
 541   Float_t maxThresh      = recParam->GetClusMaxThresh();
 542   // Threshold value for the digit signal
 543   Float_t sigThresh      = recParam->GetClusSigThresh();
 544
 545   // Detector wise calibration object for t0
 546   const AliTRDCalDet *calT0Det         = calibration->GetT0Det();
 547   // Detector wise calibration object for the gain factors
 548   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 549
 550   // Iteration limit for unfolding procedure
 551   const Float_t kEpsilon = 0.01;
 552   const Int_t   kNclus   = 3;
 553   const Int_t   kNsig    = 5;
 554
 555   Int_t    iUnfold       = 0;
 556   Double_t ratioLeft     = 1.0;
 557   Double_t ratioRight    = 1.0;
 558
 559   Double_t padSignal[kNsig];
 560   Double_t clusterSignal[kNclus];
 561   Double_t clusterPads[kNclus];
 562
 563   Int_t icham = indexesIn->GetChamber();
 564   Int_t iplan = indexesIn->GetPlane();
 565   Int_t isect = indexesIn->GetSM();
 566
 567   // Start clustering in the chamber
 568
 569   Int_t    idet    = fGeometry->GetDetector(iplan,icham,isect);
 570   if (idet != det)
 571     {
 572       AliError("Strange Detector number Missmatch!");
 573       return kFALSE;
 574     }
 575
 576   Int_t    ilayer  = AliGeomManager::kTRD1 + iplan;
 577   Int_t    imodule = icham + AliTRDgeometry::Ncham() * isect;
 578   UShort_t volid   = AliGeomManager::LayerToVolUID(ilayer,imodule);
 579
 580   //Int_t nRowMax = digitsIn->GetNrow();
 581   Int_t nColMax = digitsIn->GetNcol();
 582   Int_t nTimeTotal = digitsIn->GetNtime();
 583
 584   AliTRDpadPlane *padPlane = fGeometry->GetPadPlane(iplan,icham);
 585
 586   // Calibration object with pad wise values for t0
 587   AliTRDCalROC *calT0ROC              = calibration->GetT0ROC(idet);
 588   // Calibration object with pad wise values for the gain factors
 589   AliTRDCalROC *calGainFactorROC      = calibration->GetGainFactorROC(idet);
 590   // Calibration value for chamber wise t0
 591   Float_t       calT0DetValue         = calT0Det->GetValue(idet);
 592   // Calibration value for chamber wise gain factor
 593   Float_t       calGainFactorDetValue = calGainFactorDet->GetValue(idet);
 594
 595   Int_t nClusters      = 0;
 596
 597   // Apply the gain and the tail cancelation via digital filter
 598   AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow()
 599                                                      ,digitsIn->GetNcol()
 600                                                      ,digitsIn->GetNtime());
 601
 602 //   AliInfo(Form("nrows %d cols %d time %d",
 603 //             digitsIn->GetNrow()
 604 //             ,digitsIn->GetNcol()
 605 //             ,digitsIn->GetNtime()));
 606
 607   ResetHelperIndexes(indexesIn);
 608
 609   Transform(digitsIn
 610             ,digitsOut
 611             ,indexesIn
 612             ,fIndexesOut
 613             ,nTimeTotal
 614             ,ADCthreshold
 615             ,calGainFactorROC
 616             ,calGainFactorDetValue);
 617
 618   Int_t row   = 0;
 619   Int_t col   = 0;
 620   Int_t time  = 0;
 621   Int_t iPad  = 0;
 622
 623   fIndexesOut->ResetCounters();
 624   while (fIndexesOut->NextRCTbinIndex(row, col, time))
 625     {
 626       Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
 627
 628       // Look for the maximum
 629       if (signalM >= maxThresh)
 630         {
 631
 632           if (col + 1 >= nColMax || col-1 < 0)
 633             continue;
 634
 635           Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1  ,time));
 636
 637           Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
 638
 639           if ((TMath::Abs(signalL) <= signalM) &&
 640               (TMath::Abs(signalR) <  signalM))
 641             {
 642               if ((TMath::Abs(signalL) >= sigThresh) ||
 643                   (TMath::Abs(signalR) >= sigThresh))
 644                 {
 645                   // Maximum found, mark the position by a negative signal
 646                   digitsOut->SetDataUnchecked(row,col,time,-signalM);
 647                   fIndexesMaxima->AddIndexTBin(row,col,time);
 648                 }
 649             }
 650         }
 651     }
 652
 653   // The index to the first cluster of a given ROC
 654   Int_t firstClusterROC = -1;
 655   // The number of cluster in a given ROC
 656   Int_t nClusterROC     =  0;
 657
 658   // Now check the maxima and calculate the cluster position
 659   fIndexesMaxima->ResetCounters();
 660   while (fIndexesMaxima->NextRCTbinIndex(row, col, time))
 661     {
 662       // Maximum found ?
 663       if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) {
 664
 665         for (iPad = 0; iPad < kNclus; iPad++) {
 666           Int_t iPadCol = col - 1 + iPad;
 667           clusterSignal[iPad] =
 668             TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
 669         }
 670
 671         // Count the number of pads in the cluster
 672         Int_t nPadCount = 0;
 673         Int_t ii;
 674         // Look to the left
 675         ii = 0;
 676         while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii  ,time)) >= sigThresh) {
 677           nPadCount++;
 678           ii++;
 679           if (col-ii   <        0) break;
 680         }
 681         // Look to the right
 682         ii = 0;
 683         while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh) {
 684           nPadCount++;
 685           ii++;
 686           if (col+ii+1 >= nColMax) break;
 687         }
 688         nClusters++;
 689
 690         // Look for 5 pad cluster with minimum in the middle
 691         Bool_t fivePadCluster = kFALSE;
 692         if (col < (nColMax - 3)) {
 693           if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
 694             fivePadCluster = kTRUE;
 695           }
 696           if ((fivePadCluster) && (col < (nColMax - 5))) {
 697             if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) {
 698               fivePadCluster = kFALSE;
 699             }
 700           }
 701           if ((fivePadCluster) && (col >             1)) {
 702             if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) {
 703               fivePadCluster = kFALSE;
 704             }
 705           }
 706         }
 707
 708         // 5 pad cluster
 709         // Modify the signal of the overlapping pad for the left part
 710         // of the cluster which remains from a previous unfolding
 711         if (iUnfold) {
 712           clusterSignal[0] *= ratioLeft;
 713           iUnfold = 0;
 714         }
 715
 716         // Unfold the 5 pad cluster
 717         if (fivePadCluster) {
 718           for (iPad = 0; iPad < kNsig; iPad++) {
 719             padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
 720                                                                      ,col-1+iPad
 721                                                                      ,time));
 722           }
 723           // Unfold the two maxima and set the signal on
 724           // the overlapping pad to the ratio
 725           ratioRight        = Unfold(kEpsilon,iplan,padSignal);
 726           ratioLeft         = 1.0 - ratioRight;
 727           clusterSignal[2] *= ratioRight;
 728           iUnfold = 1;
 729         }
 730
 731         Double_t clusterCharge = clusterSignal[0]
 732           + clusterSignal[1]
 733           + clusterSignal[2];
 734
 735         // The position of the cluster
 736         clusterPads[0] =  row + 0.5;
 737         // Take the shift of the additional time bins into account
 738         clusterPads[2] = time + 0.5;
 739
 740         if (recParam->LUTOn()) {
 741           // Calculate the position of the cluster by using the
 742           // lookup table method
 743           clusterPads[1] = recParam->LUTposition(iplan,clusterSignal[0]
 744                                                  ,clusterSignal[1]
 745                                                  ,clusterSignal[2]);
 746         }
 747         else {
 748           // Calculate the position of the cluster by using the
 749           // center of gravity method
 750           for (Int_t i = 0; i < kNsig; i++) {
 751             padSignal[i] = 0.0;
 752           }
 753           padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col  ,time)); // Central pad
 754           padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left    pad
 755           padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right   pad
 756           if ((col >           2) &&
 757               (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) {
 758             padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
 759           }
 760           if ((col < nColMax - 3) &&
 761               (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])) {
 762             padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
 763           }
 764           clusterPads[1] = GetCOG(padSignal);
 765         }
 766
 767         Double_t q0 = clusterSignal[0];
 768         Double_t q1 = clusterSignal[1];
 769         Double_t q2 = clusterSignal[2];
 770         Double_t clusterSigmaY2 = (q1 * (q0 + q2) + 4.0 * q0 * q2)
 771           / (clusterCharge*clusterCharge);
 772
 773         //
 774         // Calculate the position and the error
 775         //
 776
 777         // Correct for t0 (sum of chamber and pad wise values !!!)
 778         Float_t  calT0ROCValue  = calT0ROC->GetValue(col,row);
 779         Char_t   clusterTimeBin = ((Char_t) TMath::Nint(time - (calT0DetValue + calT0ROCValue)));
 780         Double_t colSize        = padPlane->GetColSize(col);
 781         Double_t rowSize        = padPlane->GetRowSize(row);
 782
 783         Float_t clusterPos[3];
 784         clusterPos[0] = padPlane->GetColPos(col) - (clusterPads[1] + 0.5) * colSize;
 785         clusterPos[1] = padPlane->GetRowPos(row) - 0.5                    * rowSize;
 786         clusterPos[2] = CalcXposFromTimebin(clusterPads[2],idet,col,row);
 787         Float_t clusterSig[2];
 788         clusterSig[0] = (clusterSigmaY2 + 1.0/12.0) * colSize*colSize;
 789         clusterSig[1] = rowSize * rowSize / 12.0;
 790
 791         // Store the amplitudes of the pads in the cluster for later analysis
 792         Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
 793         for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
 794           if ((jPad <          0) ||
 795               (jPad >= nColMax-1)) {
 796             continue;
 797           }
 798           signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
 799         }
 800
 801         // Add the cluster to the output array
 802         // The track indices will be stored later
 803         AliTRDcluster *cluster = new AliTRDcluster(idet
 804                                                    ,clusterCharge
 805                                                    ,clusterPos
 806                                                    ,clusterSig
 807                                                    ,0x0
 808                                                    ,((Char_t) nPadCount)
 809                                                    ,signals
 810                                                    ,((UChar_t) col)
 811                                                    ,clusterTimeBin
 812                                                    ,clusterPads[1]
 813                                                    ,volid);
 814         // Temporarily store the row, column and time bin of the center pad
 815         // Used to later on assign the track indices
 816         cluster->SetLabel( row,0);
 817         cluster->SetLabel( col,1);
 818         cluster->SetLabel(time,2);
 819         RecPoints()->Add(cluster);
 820
 821         // Store the index of the first cluster in the current ROC
 822         if (firstClusterROC < 0) {
 823           firstClusterROC = RecPoints()->GetEntriesFast() - 1;
 824         }
 825         // Count the number of cluster in the current ROC
 826         nClusterROC++;
 827
 828       } // if: Maximum found ?
 829
 830     }
 831   delete digitsOut;
 832   //delete fIndexesOut;
 833   //delete fIndexesMaxima;
 834
 835   if (fAddLabels)
 836     AddLabels(idet, firstClusterROC, nClusterROC);
 837
 838   // Write the cluster and reset the array
 839   WriteClusters(idet);
 840   ResetRecPoints();
 841
 842   return kTRUE;
 843
 844 }
 845
 846 //_____________________________________________________________________________
 847 Bool_t AliTRDclusterizerV2::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
 848 {
 849   //
 850   // Add the track indices to the found clusters
 851   //
 852
 853   const Int_t   kNclus   = 3;
 854   const Int_t   kNdict   = AliTRDdigitsManager::kNDict;
 855   const Int_t   kNtrack  = kNdict * kNclus;
 856
 857   Int_t    iClusterROC   = 0;
 858
 859   Int_t row   = 0;
 860   Int_t col   = 0;
 861   Int_t time  = 0;
 862   Int_t iPad  = 0;
 863
 864   // Temporary array to collect the track indices
 865   Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
 866
 867   // Loop through the dictionary arrays one-by-one
 868   // to keep memory consumption low
 869   AliTRDdataArrayI *tracksIn = 0;
 870   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
 871
 872     tracksIn = fDigitsManager->GetDictionary(idet,iDict);
 873     // tracksIn should be expanded beforehand!
 874
 875     // Loop though the clusters found in this ROC
 876     for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
 877
 878       AliTRDcluster *cluster = (AliTRDcluster *)
 879         RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
 880       row  = cluster->GetLabel(0);
 881       col  = cluster->GetLabel(1);
 882       time = cluster->GetLabel(2);
 883
 884       for (iPad = 0; iPad < kNclus; iPad++) {
 885         Int_t iPadCol = col - 1 + iPad;
 886         Int_t index   = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
 887         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
 888       }
 889
 890     }
 891
 892     // Compress the arrays
 893     // no do not compress - we will delete them when we are done with the detector
 894     //tracksIn->Compress(1,0);
 895
 896   }
 897
 898   // Copy the track indices into the cluster
 899   // Loop though the clusters found in this ROC
 900   for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
 901
 902     AliTRDcluster *cluster = (AliTRDcluster *)
 903       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
 904     cluster->SetLabel(-9999,0);
 905     cluster->SetLabel(-9999,1);
 906     cluster->SetLabel(-9999,2);
 907
 908     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
 909
 910   }
 911
 912   delete [] idxTracks;
 913
 914   return kTRUE;
 915 }
 916
 917 //_____________________________________________________________________________
 918 Double_t AliTRDclusterizerV2::GetCOG(Double_t signal[5])
 919 {
 920   //
 921   // Get COG position
 922   // Used for clusters with more than 3 pads - where LUT not applicable
 923   //
 924
 925   Double_t sum = signal[0]
 926                + signal[1]
 927                + signal[2]
 928                + signal[3]
 929                + signal[4];
 930
 931   Double_t res = (0.0 * (-signal[0] + signal[4])
 932                       + (-signal[1] + signal[3])) / sum;
 933
 934   return res;
 935
 936 }
 937
 938 //_____________________________________________________________________________
 939 Double_t AliTRDclusterizerV2::Unfold(Double_t eps, Int_t plane, Double_t *padSignal)
 940 {
 941   //
 942   // Method to unfold neighbouring maxima.
 943   // The charge ratio on the overlapping pad is calculated
 944   // until there is no more change within the range given by eps.
 945   // The resulting ratio is then returned to the calling method.
 946   //
 947
 948   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 949   if (!calibration) {
 950     AliError("No AliTRDcalibDB instance available\n");
 951     return kFALSE;
 952   }
 953
 954   Int_t   irc                = 0;
 955   Int_t   itStep             = 0;                 // Count iteration steps
 956
 957   Double_t ratio             = 0.5;               // Start value for ratio
 958   Double_t prevRatio         = 0.0;               // Store previous ratio
 959
 960   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
 961   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
 962   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
 963
 964   // Start the iteration
 965   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
 966
 967     itStep++;
 968     prevRatio = ratio;
 969
 970     // Cluster position according to charge ratio
 971     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
 972                       / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
 973     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
 974                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
 975
 976     // Set cluster charge ratio
 977     irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
 978     Double_t ampLeft  = padSignal[1] / newSignal[1];
 979     irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
 980     Double_t ampRight = padSignal[3] / newSignal[1];
 981
 982     // Apply pad response to parameters
 983     irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
 984     irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
 985
 986     // Calculate new overlapping ratio
 987     ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] /
 988                                     (newLeftSignal[2] + newRightSignal[0]));
 989
 990   }
 991
 992   return ratio;
 993
 994 }
 995
 996 //_____________________________________________________________________________
 997 void AliTRDclusterizerV2::Transform(AliTRDdataArrayI *digitsIn
 998                                     , AliTRDdataArrayF *digitsOut
 999                                     , AliTRDSignalIndex *indexesIn
1000                                     , AliTRDSignalIndex *indexesOut
1001                                     , Int_t nTimeTotal
1002                                     , Float_t ADCthreshold
1003                                     , AliTRDCalROC *calGainFactorROC
1004                                     , Float_t calGainFactorDetValue)
1005 {
1006   //
1007   // Apply gain factor
1008   // Apply tail cancelation: Transform digitsIn to digitsOut
1009   //
1010
1011   Int_t iRow  = 0;
1012   Int_t iCol  = 0;
1013   Int_t iTime = 0;
1014
1015   AliTRDRecParam *recParam = AliTRDRecParam::Instance();
1016   if (!recParam) {
1017     AliError("No AliTRDRecParam instance available\n");
1018     return;
1019   }
1020
1021   Double_t *inADC  = new Double_t[nTimeTotal];  // ADC data before tail cancellation
1022   Double_t *outADC = new Double_t[nTimeTotal];  // ADC data after tail cancellation
1023   indexesIn->ResetCounters();
1024   while (indexesIn->NextRCIndex(iRow, iCol))
1025     {
1026       Float_t  calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1027       Double_t gain                  = calGainFactorDetValue
1028                                      * calGainFactorROCValue;
1029
1030       for (iTime = 0; iTime < nTimeTotal; iTime++)
1031         {
1032           //
1033           // Add gain
1034           //
1035           inADC[iTime]   = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1036           inADC[iTime]  /= gain;
1037           outADC[iTime]  = inADC[iTime];
1038         }
1039
1040       // Apply the tail cancelation via the digital filter
1041       if (recParam->TCOn()) {
1042         DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp());
1043       }
1044
1045       indexesIn->ResetTbinCounter();
1046       while (indexesIn->NextTbinIndex(iTime))
1047         {
1048           // Store the amplitude of the digit if above threshold
1049           if (outADC[iTime] > ADCthreshold)
1050             {
1051               digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1052               //AliDebug(5, Form("add index %d", indexesIn->GetDetNumber()));
1053               indexesOut->AddIndexTBin(iRow,iCol,iTime);
1054             }
1055         } //while itime
1056     }//while irow icol
1057
1058   delete [] inADC;
1059   delete [] outADC;
1060
1061   //AliDebug(5, Form("Stop %d", indexesIn->GetDetNumber()));
1062
1063   return;
1064
1065 }
1066
1067 //_____________________________________________________________________________
1068 void AliTRDclusterizerV2::DeConvExp(Double_t *source, Double_t *target
1069                                   , Int_t n, Int_t nexp)
1070 {
1071   //
1072   // Tail cancellation by deconvolution for PASA v4 TRF
1073   //
1074
1075   Double_t rates[2];
1076   Double_t coefficients[2];
1077
1078   // Initialization (coefficient = alpha, rates = lambda)
1079   Double_t R1 = 1.0;
1080   Double_t R2 = 1.0;
1081   Double_t C1 = 0.5;
1082   Double_t C2 = 0.5;
1083
1084   if (nexp == 1) {   // 1 Exponentials
1085     R1 = 1.156;
1086     R2 = 0.130;
1087     C1 = 0.066;
1088     C2 = 0.000;
1089   }
1090   if (nexp == 2) {   // 2 Exponentials
1091     R1 = 1.156;
1092     R2 = 0.130;
1093     C1 = 0.114;
1094     C2 = 0.624;
1095   }
1096
1097   coefficients[0] = C1;
1098   coefficients[1] = C2;
1099
1100   Double_t Dt = 0.1;
1101
1102   rates[0] = TMath::Exp(-Dt/(R1));
1103   rates[1] = TMath::Exp(-Dt/(R2));
1104
1105   Int_t i = 0;
1106   Int_t k = 0;
1107
1108   Double_t reminder[2];
1109   Double_t correction;
1110   Double_t result;
1111
1112   // Attention: computation order is important
1113   correction = 0.0;
1114   for (k = 0; k < nexp; k++) {
1115     reminder[k] = 0.0;
1116   }
1117   for (i = 0; i < n; i++) {
1118     result    = (source[i] - correction);    // No rescaling
1119     target[i] = result;
1120
1121     for (k = 0; k < nexp; k++) {
1122       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1123     }
1124     correction = 0.0;
1125     for (k = 0; k < nexp; k++) {
1126       correction += reminder[k];
1127     }
1128   }
1129
1130 }