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