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