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