TRD/AliTRDclusterizer.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 // TRD cluster finder                                                        //
  21 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TF1.h>
  25 #include <TTree.h>
  26 #include <TH1.h>
  27 #include <TFile.h>
  28 #include <TClonesArray.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 "AliTRDReconstructor.h"
  40 #include "AliTRDgeometry.h"
  41 #include "AliTRDarraySignal.h"
  42 #include "AliTRDarrayDictionary.h"
  43 #include "AliTRDarrayADC.h"
  44 #include "AliTRDdigitsManager.h"
  45 #include "AliTRDrawData.h"
  46 #include "AliTRDcalibDB.h"
  47 #include "AliTRDrecoParam.h"
  48 #include "AliTRDCommonParam.h"
  49 #include "AliTRDtransform.h"
  50 #include "AliTRDSignalIndex.h"
  51 #include "AliTRDrawStreamBase.h"
  52 #include "AliTRDfeeParam.h"
  53
  54 #include "TTreeStream.h"
  55
  56 #include "Cal/AliTRDCalROC.h"
  57 #include "Cal/AliTRDCalDet.h"
  58 #include "Cal/AliTRDCalSingleChamberStatus.h"
  59
  60 ClassImp(AliTRDclusterizer)
  61
  62 //_____________________________________________________________________________
  63 AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
  64   :TNamed()
  65   ,fReconstructor(rec)
  66   ,fRunLoader(NULL)
  67   ,fClusterTree(NULL)
  68   ,fRecPoints(NULL)
  69   ,fTrackletTree(NULL)
  70   ,fDigitsManager(new AliTRDdigitsManager(rec))
  71   ,fTrackletContainer(NULL)
  72   ,fAddLabels(kTRUE)
  73   ,fRawVersion(2)
  74   ,fTransform(new AliTRDtransform(0))
  75   ,fLUTbin(0)
  76   ,fLUT(NULL)
  77   ,fDigitsIn(NULL)
  78   ,fIndexes(NULL)
  79   ,fADCthresh(0)
  80   ,fMaxThresh(0)
  81   ,fSigThresh(0)
  82   ,fMinMaxCutSigma(0)
  83   ,fMinLeftRightCutSigma(0)
  84   ,fLayer(0)
  85   ,fDet(0)
  86   ,fVolid(0)
  87   ,fColMax(0)
  88   ,fTimeTotal(0)
  89   ,fCalGainFactorROC(NULL)
  90   ,fCalGainFactorDetValue(0)
  91   ,fCalNoiseROC(NULL)
  92   ,fCalNoiseDetValue(0)
  93   ,fDigitsOut(NULL)
  94   ,fClusterROC(0)
  95   ,firstClusterROC(0)
  96 {
  97   //
  98   // AliTRDclusterizer default constructor
  99   //
 100
 101   AliTRDcalibDB *trd = 0x0;
 102   if (!(trd = AliTRDcalibDB::Instance())) {
 103     AliFatal("Could not get calibration object");
 104   }
 105
 106   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
 107
 108   // Initialize debug stream
 109   if(fReconstructor){
 110     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 1){
 111       TDirectory *savedir = gDirectory;
 112       //fgGetDebugStream    = new TTreeSRedirector("TRD.ClusterizerDebug.root");
 113       savedir->cd();
 114     }
 115   }
 116
 117 }
 118
 119 //_____________________________________________________________________________
 120 AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, const AliTRDReconstructor *const rec)
 121   :TNamed(name,title)
 122   ,fReconstructor(rec)
 123   ,fRunLoader(NULL)
 124   ,fClusterTree(NULL)
 125   ,fRecPoints(NULL)
 126   ,fTrackletTree(NULL)
 127   ,fDigitsManager(new AliTRDdigitsManager(rec))
 128   ,fTrackletContainer(NULL)
 129   ,fAddLabels(kTRUE)
 130   ,fRawVersion(2)
 131   ,fTransform(new AliTRDtransform(0))
 132   ,fLUTbin(0)
 133   ,fLUT(NULL)
 134   ,fDigitsIn(NULL)
 135   ,fIndexes(NULL)
 136   ,fADCthresh(0)
 137   ,fMaxThresh(0)
 138   ,fSigThresh(0)
 139   ,fMinMaxCutSigma(0)
 140   ,fMinLeftRightCutSigma(0)
 141   ,fLayer(0)
 142   ,fDet(0)
 143   ,fVolid(0)
 144   ,fColMax(0)
 145   ,fTimeTotal(0)
 146   ,fCalGainFactorROC(NULL)
 147   ,fCalGainFactorDetValue(0)
 148   ,fCalNoiseROC(NULL)
 149   ,fCalNoiseDetValue(0)
 150   ,fDigitsOut(NULL)
 151   ,fClusterROC(0)
 152   ,firstClusterROC(0)
 153 {
 154   //
 155   // AliTRDclusterizer constructor
 156   //
 157
 158   AliTRDcalibDB *trd = 0x0;
 159   if (!(trd = AliTRDcalibDB::Instance())) {
 160     AliFatal("Could not get calibration object");
 161   }
 162
 163   fDigitsManager->CreateArrays();
 164
 165   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
 166
 167   FillLUT();
 168
 169 }
 170
 171 //_____________________________________________________________________________
 172 AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
 173   :TNamed(c)
 174   ,fReconstructor(c.fReconstructor)
 175   ,fRunLoader(NULL)
 176   ,fClusterTree(NULL)
 177   ,fRecPoints(NULL)
 178   ,fTrackletTree(NULL)
 179   ,fDigitsManager(NULL)
 180   ,fTrackletContainer(NULL)
 181   ,fAddLabels(kTRUE)
 182   ,fRawVersion(2)
 183   ,fTransform(NULL)
 184   ,fLUTbin(0)
 185   ,fLUT(0)
 186   ,fDigitsIn(NULL)
 187   ,fIndexes(NULL)
 188   ,fADCthresh(0)
 189   ,fMaxThresh(0)
 190   ,fSigThresh(0)
 191   ,fMinMaxCutSigma(0)
 192   ,fMinLeftRightCutSigma(0)
 193   ,fLayer(0)
 194   ,fDet(0)
 195   ,fVolid(0)
 196   ,fColMax(0)
 197   ,fTimeTotal(0)
 198   ,fCalGainFactorROC(NULL)
 199   ,fCalGainFactorDetValue(0)
 200   ,fCalNoiseROC(NULL)
 201   ,fCalNoiseDetValue(0)
 202   ,fDigitsOut(NULL)
 203   ,fClusterROC(0)
 204   ,firstClusterROC(0)
 205 {
 206   //
 207   // AliTRDclusterizer copy constructor
 208   //
 209
 210   FillLUT();
 211
 212 }
 213
 214 //_____________________________________________________________________________
 215 AliTRDclusterizer::~AliTRDclusterizer()
 216 {
 217   //
 218   // AliTRDclusterizer destructor
 219   //
 220
 221   if (fRecPoints/* && IsClustersOwner()*/){
 222     fRecPoints->Delete();
 223     delete fRecPoints;
 224   }
 225
 226   if (fDigitsManager) {
 227     delete fDigitsManager;
 228     fDigitsManager = NULL;
 229   }
 230
 231   if (fTrackletContainer){
 232     delete fTrackletContainer;
 233     fTrackletContainer = NULL;
 234   }
 235
 236   if (fTransform){
 237     delete fTransform;
 238     fTransform     = NULL;
 239   }
 240
 241   if (fLUT) {
 242     delete [] fLUT;
 243     fLUT           = NULL;
 244   }
 245
 246   if (fDigitsOut) {
 247     delete fDigitsOut;
 248     fDigitsOut     = NULL;
 249   }
 250
 251 }
 252
 253 //_____________________________________________________________________________
 254 AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
 255 {
 256   //
 257   // Assignment operator
 258   //
 259
 260   if (this != &c)
 261     {
 262       ((AliTRDclusterizer &) c).Copy(*this);
 263     }
 264
 265   return *this;
 266
 267 }
 268
 269 //_____________________________________________________________________________
 270 void AliTRDclusterizer::Copy(TObject &c) const
 271 {
 272   //
 273   // Copy function
 274   //
 275
 276   ((AliTRDclusterizer &) c).fClusterTree   = NULL;
 277   ((AliTRDclusterizer &) c).fRecPoints     = NULL;
 278   ((AliTRDclusterizer &) c).fTrackletTree  = NULL;
 279   ((AliTRDclusterizer &) c).fDigitsManager = NULL;
 280   ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
 281   ((AliTRDclusterizer &) c).fAddLabels     = fAddLabels;
 282   ((AliTRDclusterizer &) c).fRawVersion    = fRawVersion;
 283   ((AliTRDclusterizer &) c).fTransform     = NULL;
 284   ((AliTRDclusterizer &) c).fLUTbin        = 0;
 285   ((AliTRDclusterizer &) c).fLUT           = NULL;
 286   ((AliTRDclusterizer &) c).fDigitsIn      = NULL;
 287   ((AliTRDclusterizer &) c).fIndexes       = NULL;
 288   ((AliTRDclusterizer &) c).fADCthresh     = 0;
 289   ((AliTRDclusterizer &) c).fMaxThresh     = 0;
 290   ((AliTRDclusterizer &) c).fSigThresh     = 0;
 291   ((AliTRDclusterizer &) c).fMinMaxCutSigma= 0;
 292   ((AliTRDclusterizer &) c).fMinLeftRightCutSigma = 0;
 293   ((AliTRDclusterizer &) c).fLayer         = 0;
 294   ((AliTRDclusterizer &) c).fDet           = 0;
 295   ((AliTRDclusterizer &) c).fVolid         = 0;
 296   ((AliTRDclusterizer &) c).fColMax        = 0;
 297   ((AliTRDclusterizer &) c).fTimeTotal     = 0;
 298   ((AliTRDclusterizer &) c).fCalGainFactorROC = NULL;
 299   ((AliTRDclusterizer &) c).fCalGainFactorDetValue = 0;
 300   ((AliTRDclusterizer &) c).fCalNoiseROC   = NULL;
 301   ((AliTRDclusterizer &) c).fCalNoiseDetValue = 0;
 302   ((AliTRDclusterizer &) c).fDigitsOut     = NULL;
 303   ((AliTRDclusterizer &) c).fClusterROC    = 0;
 304   ((AliTRDclusterizer &) c).firstClusterROC= 0;
 305
 306 }
 307
 308 //_____________________________________________________________________________
 309 Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
 310 {
 311   //
 312   // Opens the AliROOT file. Output and input are in the same file
 313   //
 314
 315   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 316   fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
 317
 318   if (!fRunLoader) {
 319     fRunLoader = AliRunLoader::Open(name);
 320   }
 321
 322   if (!fRunLoader) {
 323     AliError(Form("Can not open session for file %s.",name));
 324     return kFALSE;
 325   }
 326
 327   OpenInput(nEvent);
 328   OpenOutput();
 329
 330   return kTRUE;
 331
 332 }
 333
 334 //_____________________________________________________________________________
 335 Bool_t AliTRDclusterizer::OpenOutput()
 336 {
 337   //
 338   // Open the output file
 339   //
 340
 341   if (!fReconstructor->IsWritingClusters()) return kTRUE;
 342
 343   TObjArray *ioArray = 0x0;
 344
 345   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 346   loader->MakeTree("R");
 347
 348   fClusterTree = loader->TreeR();
 349   fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
 350
 351   return kTRUE;
 352
 353 }
 354
 355 //_____________________________________________________________________________
 356 Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
 357 {
 358   //
 359   // Connect the output tree
 360   //
 361
 362   // clusters writing
 363   if (fReconstructor->IsWritingClusters()){
 364     TObjArray *ioArray = 0x0;
 365     fClusterTree = clusterTree;
 366     fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
 367   }
 368
 369   // tracklet writing
 370   if (fReconstructor->IsWritingTracklets()){
 371     TString evfoldname = AliConfig::GetDefaultEventFolderName();
 372     fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
 373
 374     if (!fRunLoader) {
 375       fRunLoader = AliRunLoader::Open("galice.root");
 376     }
 377     if (!fRunLoader) {
 378       AliError(Form("Can not open session for file galice.root."));
 379       return kFALSE;
 380     }
 381
 382     UInt_t **leaves = new UInt_t *[2];
 383     AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 384     if (!dl) {
 385       AliError("Could not get the tracklets data loader!");
 386       dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
 387       fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl);
 388     }
 389     else {
 390       fTrackletTree = dl->Tree();
 391       if (!fTrackletTree)
 392         {
 393         dl->MakeTree();
 394         fTrackletTree = dl->Tree();
 395         }
 396       TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
 397       if (!trkbranch)
 398         fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
 399     }
 400   }
 401
 402   return kTRUE;
 403
 404 }
 405
 406 //_____________________________________________________________________________
 407 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
 408 {
 409   //
 410   // Opens a ROOT-file with TRD-hits and reads in the digits-tree
 411   //
 412
 413   // Import the Trees for the event nEvent in the file
 414   fRunLoader->GetEvent(nEvent);
 415
 416   return kTRUE;
 417
 418 }
 419
 420 //_____________________________________________________________________________
 421 Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
 422 {
 423   //
 424   // Fills TRDcluster branch in the tree with the clusters
 425   // found in detector = det. For det=-1 writes the tree.
 426   //
 427
 428   if ((det <                      -1) ||
 429       (det >= AliTRDgeometry::Ndet())) {
 430     AliError(Form("Unexpected detector index %d.\n",det));
 431     return kFALSE;
 432   }
 433
 434   TObjArray *ioArray = new TObjArray(400);
 435   TBranch *branch = fClusterTree->GetBranch("TRDcluster");
 436   if (!branch) {
 437     branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
 438   } else branch->SetAddress(&ioArray);
 439
 440   Int_t nRecPoints = RecPoints()->GetEntriesFast();
 441   if(det >= 0){
 442     for (Int_t i = 0; i < nRecPoints; i++) {
 443       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
 444       if(det != c->GetDetector()) continue;
 445       ioArray->AddLast(c);
 446     }
 447     fClusterTree->Fill();
 448   } else {
 449
 450     Int_t detOld = -1;
 451     for (Int_t i = 0; i < nRecPoints; i++) {
 452       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
 453       if(c->GetDetector() != detOld){
 454         fClusterTree->Fill();
 455         ioArray->Clear();
 456         detOld = c->GetDetector();
 457       }
 458       ioArray->AddLast(c);
 459     }
 460   }
 461   delete ioArray;
 462
 463   return kTRUE;
 464
 465 }
 466
 467 //_____________________________________________________________________________
 468 Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
 469 {
 470   //
 471   // Write the raw data tracklets into seperate file
 472   //
 473
 474   UInt_t **leaves = new UInt_t *[2];
 475   for (Int_t i=0; i<2 ;i++){
 476     leaves[i] = new UInt_t[258];
 477     leaves[i][0] = det; // det
 478     leaves[i][1] = i;   // side
 479     memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
 480   }
 481
 482   if (!fTrackletTree){
 483     AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 484     dl->MakeTree();
 485     fTrackletTree = dl->Tree();
 486   }
 487
 488   TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
 489   if (!trkbranch) {
 490     trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
 491   }
 492
 493   for (Int_t i=0; i<2; i++){
 494     if (leaves[i][2]>0) {
 495       trkbranch->SetAddress(leaves[i]);
 496       fTrackletTree->Fill();
 497     }
 498   }
 499
 500   AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 501   dl->WriteData("OVERWRITE");
 502   //dl->Unload();
 503   delete [] leaves;
 504
 505   return kTRUE;
 506
 507 }
 508
 509 //_____________________________________________________________________________
 510 Bool_t AliTRDclusterizer::ReadDigits()
 511 {
 512   //
 513   // Reads the digits arrays from the input aliroot file
 514   //
 515
 516   if (!fRunLoader) {
 517     AliError("No run loader available");
 518     return kFALSE;
 519   }
 520
 521   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 522   if (!loader->TreeD()) {
 523     loader->LoadDigits();
 524   }
 525
 526   // Read in the digit arrays
 527   return (fDigitsManager->ReadDigits(loader->TreeD()));
 528
 529 }
 530
 531 //_____________________________________________________________________________
 532 Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
 533 {
 534   //
 535   // Reads the digits arrays from the input tree
 536   //
 537
 538   // Read in the digit arrays
 539   return (fDigitsManager->ReadDigits(digitsTree));
 540
 541 }
 542
 543 //_____________________________________________________________________________
 544 Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
 545 {
 546   //
 547   // Reads the digits arrays from the ddl file
 548   //
 549
 550   AliTRDrawData raw;
 551   fDigitsManager = raw.Raw2Digits(rawReader);
 552
 553   return kTRUE;
 554
 555 }
 556
 557 //_____________________________________________________________________________
 558 Bool_t AliTRDclusterizer::MakeClusters()
 559 {
 560   //
 561   // Creates clusters from digits
 562   //
 563
 564   // Propagate info from the digits manager
 565   if (fAddLabels == kTRUE){
 566     fAddLabels = fDigitsManager->UsesDictionaries();
 567   }
 568
 569   Bool_t fReturn = kTRUE;
 570   for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
 571
 572     AliTRDarrayADC *digitsIn = (AliTRDarrayADC*) fDigitsManager->GetDigits(i); //mod
 573     // This is to take care of switched off super modules
 574     if (!digitsIn->HasData()) continue;
 575     digitsIn->Expand();
 576     digitsIn->DeleteNegatives();  // Restore digits array to >=0 values
 577     AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
 578     if (indexes->IsAllocated() == kFALSE){
 579       fDigitsManager->BuildIndexes(i);
 580     }
 581
 582     Bool_t fR = kFALSE;
 583     if (indexes->HasEntry()){
 584       if (fAddLabels){
 585         for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
 586           AliTRDarrayDictionary *tracksIn = 0; //mod
 587           tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict);  //mod
 588           tracksIn->Expand();
 589         }
 590       }
 591       fR = MakeClusters(i);
 592       fReturn = fR && fReturn;
 593     }
 594
 595     //if (fR == kFALSE){
 596     //  if(IsWritingClusters()) WriteClusters(i);
 597     //  ResetRecPoints();
 598     //}
 599
 600     // No compress just remove
 601     fDigitsManager->RemoveDigits(i);
 602     fDigitsManager->RemoveDictionaries(i);
 603     fDigitsManager->ClearIndexes(i);
 604   }
 605
 606   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 607
 608   AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
 609
 610   return fReturn;
 611
 612 }
 613
 614 //_____________________________________________________________________________
 615 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
 616 {
 617   //
 618   // Creates clusters from raw data
 619   //
 620
 621   return Raw2ClustersChamber(rawReader);
 622
 623 }
 624
 625 //_____________________________________________________________________________
 626 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
 627 {
 628   //
 629   // Creates clusters from raw data
 630   //
 631
 632   // Create the digits manager
 633   if (!fDigitsManager){
 634     fDigitsManager = new AliTRDdigitsManager(fReconstructor);
 635     fDigitsManager->CreateArrays();
 636   }
 637
 638   fDigitsManager->SetUseDictionaries(fAddLabels);
 639
 640   // tracklet container for raw tracklet writing
 641   if (!fTrackletContainer && fReconstructor->IsWritingTracklets()) {
 642     // maximum tracklets for one HC
 643     const Int_t kTrackletChmb=256;
 644     fTrackletContainer = new UInt_t *[2];
 645     fTrackletContainer[0] = new UInt_t[kTrackletChmb];
 646     fTrackletContainer[1] = new UInt_t[kTrackletChmb];
 647   }
 648
 649   AliTRDrawStreamBase *input = AliTRDrawStreamBase::GetRawStream(rawReader);
 650
 651   AliInfo(Form("Stream version: %s", input->IsA()->GetName()));
 652
 653   Int_t det    = 0;
 654   while ((det = input->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
 655     Bool_t iclusterBranch = kFALSE;
 656     if (fDigitsManager->GetIndexes(det)->HasEntry()){
 657       iclusterBranch = MakeClusters(det);
 658     }
 659
 660     fDigitsManager->ResetArrays(det);
 661
 662     if (!fReconstructor->IsWritingTracklets()) continue;
 663     if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
 664   }
 665
 666   if (fReconstructor->IsWritingTracklets()){
 667     delete [] fTrackletContainer[0];
 668     delete [] fTrackletContainer[1];
 669     delete [] fTrackletContainer;
 670     fTrackletContainer = NULL;
 671   }
 672
 673   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 674
 675   delete fDigitsManager;
 676   fDigitsManager = NULL;
 677
 678   delete input;
 679   input = NULL;
 680
 681   AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
 682   return kTRUE;
 683
 684 }
 685
 686 //_____________________________________________________________________________
 687 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
 688 {
 689   //
 690   // Check if a pad is masked
 691   //
 692
 693   UChar_t status = 0;
 694
 695   if(signal>0 && TESTBIT(signal, 10)){
 696     CLRBIT(signal, 10);
 697     for(int ibit=0; ibit<4; ibit++){
 698       if(TESTBIT(signal, 11+ibit)){
 699         SETBIT(status, ibit);
 700         CLRBIT(signal, 11+ibit);
 701       }
 702     }
 703   }
 704   return status;
 705 }
 706
 707 //_____________________________________________________________________________
 708 void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out){
 709   //
 710   // Set the pad status into out
 711   // First three bits are needed for the position encoding
 712   //
 713   out |= status << 3;
 714 }
 715
 716 //_____________________________________________________________________________
 717 UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
 718   //
 719   // return the staus encoding of the corrupted pad
 720   //
 721   return static_cast<UChar_t>(encoding >> 3);
 722 }
 723
 724 //_____________________________________________________________________________
 725 Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
 726   //
 727   // Return the position of the corruption
 728   //
 729   return encoding & 7;
 730 }
 731
 732 //_____________________________________________________________________________
 733 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 734 {
 735   //
 736   // Generates the cluster.
 737   //
 738
 739   // Get the digits
 740   //   digits should be expanded beforehand!
 741   //   digitsIn->Expand();
 742   fDigitsIn = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
 743
 744   // This is to take care of switched off super modules
 745   if (!fDigitsIn->HasData())
 746     {
 747       return kFALSE;
 748     }
 749
 750   fIndexes = fDigitsManager->GetIndexes(det);
 751   if (fIndexes->IsAllocated() == kFALSE)
 752     {
 753       AliError("Indexes do not exist!");
 754       return kFALSE;
 755     }
 756
 757   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
 758   if (!calibration)
 759     {
 760       AliFatal("No AliTRDcalibDB instance available\n");
 761       return kFALSE;
 762     }
 763
 764   fADCthresh = 0;
 765
 766   if (!fReconstructor){
 767     AliError("Reconstructor not set\n");
 768     return kFALSE;
 769   }
 770
 771   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDReconstructor::kClusterizer);
 772
 773   fMaxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
 774   fSigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
 775   fMinMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
 776   fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
 777
 778   Int_t istack  = fIndexes->GetStack();
 779   fLayer  = fIndexes->GetLayer();
 780   Int_t isector = fIndexes->GetSM();
 781
 782   // Start clustering in the chamber
 783
 784   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
 785   if (fDet != det) {
 786     AliError("Strange Detector number Missmatch!");
 787     return kFALSE;
 788   }
 789
 790   // TRD space point transformation
 791   fTransform->SetDetector(det);
 792
 793   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
 794   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
 795   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
 796
 797   fColMax    = fDigitsIn->GetNcol();
 798   Int_t nRowMax    = fDigitsIn->GetNrow();
 799   fTimeTotal = fDigitsIn->GetNtime();
 800
 801   // Detector wise calibration object for the gain factors
 802   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 803   // Calibration object with pad wise values for the gain factors
 804   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
 805   // Calibration value for chamber wise gain factor
 806   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
 807
 808   // Detector wise calibration object for the noise
 809   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
 810   // Calibration object with pad wise values for the noise
 811   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
 812   // Calibration value for chamber wise noise
 813   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
 814
 815   if(fDigitsOut) delete fDigitsOut;
 816   fDigitsOut = new AliTRDarraySignal(nRowMax, fColMax, fTimeTotal);
 817
 818   firstClusterROC = -1;
 819   fClusterROC     =  0;
 820
 821   // Apply the gain and the tail cancelation via digital filter
 822   TailCancelation();
 823
 824   MaxStruct curr, last;
 825   last.Row = -1;
 826   Int_t nMaximas = 0, nCorrupted = 0;
 827   Float_t Ratio = 1;
 828
 829   // Here the clusterfining is happening
 830
 831   for(curr.Time = 0; curr.Time < fTimeTotal; curr.Time++)
 832     while(fIndexes->NextRCIndex(curr.Row, curr.Col))
 833       if(IsMaximum(curr, curr.padStatus, &curr.Signals[0]))
 834         {
 835           if(last.Row>-1)
 836             {
 837               last.Signals[0] *= Ratio;
 838               if(curr.Row==last.Row && curr.Col==last.Col+2)
 839                 {
 840                   if(IsFivePadCluster(last, curr, Ratio))
 841                     {
 842                       last.Signals[2] *= Ratio;
 843                       Ratio = 1 - Ratio;
 844                     }else Ratio = 1;
 845                 }else Ratio = 1;
 846               CreateCluster(last);
 847             }
 848           last=curr;
 849         }
 850   if(last.Row>-1)
 851     {
 852       last.Signals[0] *= Ratio;
 853       CreateCluster(last);
 854     }
 855
 856   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 2){
 857     (*fDebugStream) << "MakeClusters"
 858                     << "Detector="   << det
 859                     << "NMaxima="    << nMaximas
 860                     << "NClusters="  << fClusterROC
 861                     << "NCorrupted=" << nCorrupted
 862                     << "\n";
 863   }
 864
 865   if (fAddLabels) {
 866     AddLabels(fDet,firstClusterROC,fClusterROC);
 867   }
 868
 869   return kTRUE;
 870
 871 }
 872
 873 //_____________________________________________________________________________
 874 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Float_t *const Signals)
 875 {
 876   //
 877   // Returns true if this row,col,time combination is a maximum.
 878   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
 879   //
 880
 881   Signals[1] = fDigitsOut->GetData(Max.Row, Max.Col, Max.Time);
 882   if(Signals[1] < fMaxThresh) return kFALSE;
 883
 884   Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.Col, Max.Row);
 885   if (Signals[1] < noiseMiddleThresh) return kFALSE;
 886
 887   if (Max.Col + 1 >= fColMax || Max.Col < 1) return kFALSE;
 888
 889   UChar_t status[3]={fDigitsIn->GetPadStatus(Max.Row, Max.Col-1, Max.Time),
 890                      fDigitsIn->GetPadStatus(Max.Row, Max.Col,   Max.Time),
 891                      fDigitsIn->GetPadStatus(Max.Row, Max.Col+1, Max.Time)};
 892
 893   Signals[0] = fDigitsOut->GetData(Max.Row, Max.Col-1, Max.Time);
 894   Signals[2] = fDigitsOut->GetData(Max.Row, Max.Col+1, Max.Time);
 895
 896   if(!(status[0] | status[1] | status[2])) {//all pads are good
 897     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
 898       if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
 899         Float_t  noiseSumThresh = fMinLeftRightCutSigma
 900           * fCalNoiseDetValue
 901           * fCalNoiseROC->GetValue(Max.Col, Max.Row);
 902         if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
 903         padStatus = 0;
 904         return kTRUE;
 905       }
 906     }
 907   }
 908   else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
 909     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) {
 910       fDigitsOut->SetData(Max.Row, Max.Col+1, Max.Time, 0.);
 911       SetPadStatus(status[2], padStatus);
 912       return kTRUE;
 913     }
 914     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
 915       fDigitsOut->SetData(Max.Row, Max.Col-1, Max.Time, 0.);
 916       SetPadStatus(status[0], padStatus);
 917       return kTRUE;
 918     }
 919     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
 920       fDigitsOut->SetData(Max.Row, Max.Col, Max.Time, fMaxThresh);
 921       SetPadStatus(status[1], padStatus);
 922       return kTRUE;
 923     }
 924   }
 925   return kFALSE;
 926 }
 927
 928 //_____________________________________________________________________________
 929 Bool_t AliTRDclusterizer::IsFivePadCluster(const MaxStruct &ThisMax, const MaxStruct &NeighbourMax, Float_t &ratio)
 930 {
 931   //
 932   // Look for 5 pad cluster with minimum in the middle
 933   // Gives back the ratio
 934   //
 935
 936   if (ThisMax.Col >= fColMax - 3) return kFALSE;
 937   if (ThisMax.Col < fColMax - 5){
 938     if (fDigitsOut->GetData(ThisMax.Row, ThisMax.Col+4, ThisMax.Time) >= fSigThresh)
 939       return kFALSE;
 940   }
 941   if (ThisMax.Col > 1) {
 942     if (fDigitsOut->GetData(ThisMax.Row, ThisMax.Col-2, ThisMax.Time) >= fSigThresh)
 943       return kFALSE;
 944   }
 945
 946   //if (fSignalsThisMax[1] >= 0){ //TR: mod
 947
 948   const Float_t kEpsilon = 0.01;
 949   Double_t padSignal[5] = {ThisMax.Signals[0], ThisMax.Signals[1], ThisMax.Signals[2],
 950                            NeighbourMax.Signals[1], NeighbourMax.Signals[2]};
 951
 952   // Unfold the two maxima and set the signal on
 953   // the overlapping pad to the ratio
 954   ratio = Unfold(kEpsilon,fLayer,padSignal);
 955   return kTRUE;
 956 }
 957
 958 //_____________________________________________________________________________
 959 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
 960 {
 961   //
 962   // Creates a cluster at the given position and saves it in fRecPoint
 963   //
 964
 965   const Int_t kNsig = 5;
 966   Double_t padSignal[kNsig];
 967
 968   // The position of the cluster in COL direction relative to the center pad (pad units)
 969   Double_t clusterPosCol = 0.0;
 970   if (fReconstructor->GetRecoParam()->IsLUT()) {
 971     // Calculate the position of the cluster by using the
 972     // lookup table method
 973     clusterPosCol = LUTposition(fLayer,Max.Signals[0]
 974                                 ,Max.Signals[1]
 975                                 ,Max.Signals[2]);
 976   }
 977   else {
 978     // Calculate the position of the cluster by using the
 979     // center of gravity method
 980     padSignal[1] = Max.Signals[0];
 981     padSignal[2] = Max.Signals[1];
 982     padSignal[3] = Max.Signals[2];
 983     if(Max.Col > 2){
 984       padSignal[0] = fDigitsOut->GetData(Max.Row, Max.Col-2, Max.Time);
 985       if(padSignal[0]>= padSignal[1])
 986         padSignal[0] = 0;
 987     }
 988     if(Max.Col < fColMax - 3){
 989       padSignal[4] = fDigitsOut->GetData(Max.Row, Max.Col+2, Max.Time);
 990       if(padSignal[4]>= padSignal[3])
 991         padSignal[4] = 0;
 992     }
 993     clusterPosCol = GetCOG(padSignal);
 994   }
 995
 996   // Count the number of pads in the cluster
 997   Int_t nPadCount = 1;
 998   // Look to the right
 999   Int_t ii = 1;
1000   while (fDigitsOut->GetData(Max.Row, Max.Col-ii, Max.Time) >= fSigThresh) {
1001     nPadCount++;
1002     ii++;
1003     if (Max.Col < ii) break;
1004   }
1005   // Look to the left
1006   ii = 1;
1007   while (fDigitsOut->GetData(Max.Row, Max.Col+ii, Max.Time) >= fSigThresh) {
1008     nPadCount++;
1009     ii++;
1010     if (Max.Col+ii >= fColMax) break;
1011   }
1012
1013   // Store the amplitudes of the pads in the cluster for later analysis
1014   // and check whether one of these pads is masked in the database
1015   Short_t signals[7] = { 0, 0, TMath::Nint(Max.Signals[0]),
1016                                TMath::Nint(Max.Signals[1]),
1017                                TMath::Nint(Max.Signals[2]), 0, 0 };
1018   for(Int_t i = 0; i<2; i++)
1019     {
1020       if(Max.Col+i >= 3)
1021         signals[i] = TMath::Nint(fDigitsOut->GetData(Max.Row, Max.Col-3+i, Max.Time));
1022       if(Max.Col+3-i < fColMax)
1023         signals[6-i] = TMath::Nint(fDigitsOut->GetData(Max.Row, Max.Col+3-i, Max.Time));
1024     }
1025   /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1026     if ((jPad >= 0) && (jPad < fColMax))
1027       signals[jPad-Max.Col+3] = TMath::Nint(fDigitsOut->GetData(Max.Row,jPad,Max.Time));
1028       }*/
1029
1030   // Transform the local cluster coordinates into calibrated
1031   // space point positions defined in the local tracking system.
1032   // Here the calibration for T0, Vdrift and ExB is applied as well.
1033   Double_t clusterXYZ[6];
1034   clusterXYZ[0] = clusterPosCol;
1035   clusterXYZ[1] = Max.Signals[2];
1036   clusterXYZ[2] = Max.Signals[1];
1037   clusterXYZ[3] = Max.Signals[0];
1038   clusterXYZ[4] = 0.0;
1039   clusterXYZ[5] = 0.0;
1040   Int_t    clusterRCT[3];
1041   clusterRCT[0] = Max.Row;
1042   clusterRCT[1] = Max.Col;
1043   clusterRCT[2] = 0;
1044
1045   Bool_t out = kTRUE;
1046   if (fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) Max.Time),out,0)) {
1047
1048     // Add the cluster to the output array
1049     // The track indices will be stored later
1050     Float_t clusterPos[3];
1051     clusterPos[0] = clusterXYZ[0];
1052     clusterPos[1] = clusterXYZ[1];
1053     clusterPos[2] = clusterXYZ[2];
1054     Float_t clusterSig[2];
1055     clusterSig[0] = clusterXYZ[4];
1056     clusterSig[1] = clusterXYZ[5];
1057     Double_t clusterCharge  = clusterXYZ[3];
1058     Char_t   clusterTimeBin = ((Char_t) clusterRCT[2]);
1059
1060     Int_t n = RecPoints()->GetEntriesFast();
1061     AliTRDcluster *cluster = new((*RecPoints())[n]) AliTRDcluster(
1062                                                                   fDet,
1063                                                                   clusterCharge, clusterPos, clusterSig,
1064                                                                   0x0,
1065                                                                   ((Char_t) nPadCount),
1066                                                                   signals,
1067                                                                   ((UChar_t) Max.Col), ((UChar_t) Max.Row), ((UChar_t) Max.Time),
1068                                                                   clusterTimeBin, clusterPosCol,
1069                                                                   fVolid);
1070     cluster->SetInChamber(!out);
1071
1072     UChar_t maskPosition = GetCorruption(Max.padStatus);
1073     if (maskPosition) {
1074       cluster->SetPadMaskedPosition(maskPosition);
1075       cluster->SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1076     }
1077
1078     // Temporarily store the Max.Row, column and time bin of the center pad
1079     // Used to later on assign the track indices
1080     cluster->SetLabel(Max.Row, 0);
1081     cluster->SetLabel(Max.Col, 1);
1082     cluster->SetLabel(Max.Time,2);
1083
1084     // Store the index of the first cluster in the current ROC
1085     if (firstClusterROC < 0) {
1086       firstClusterROC = RecPoints()->GetEntriesFast() - 1;
1087     }
1088
1089     // Count the number of cluster in the current ROC
1090     fClusterROC++;
1091   }
1092
1093 }
1094
1095 //_____________________________________________________________________________
1096 Bool_t AliTRDclusterizer::AddLabels(const Int_t idet, const Int_t firstClusterROC, const Int_t nClusterROC)
1097 {
1098   //
1099   // Add the track indices to the found clusters
1100   //
1101
1102   const Int_t   kNclus  = 3;
1103   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1104   const Int_t   kNtrack = kNdict * kNclus;
1105
1106   Int_t iClusterROC = 0;
1107
1108   Int_t row  = 0;
1109   Int_t col  = 0;
1110   Int_t time = 0;
1111   Int_t iPad = 0;
1112
1113   // Temporary array to collect the track indices
1114   Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
1115
1116   // Loop through the dictionary arrays one-by-one
1117   // to keep memory consumption low
1118   AliTRDarrayDictionary *tracksIn = 0;  //mod
1119   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1120
1121     // tracksIn should be expanded beforehand!
1122     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(idet,iDict);
1123
1124     // Loop though the clusters found in this ROC
1125     for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1126
1127       AliTRDcluster *cluster = (AliTRDcluster *)
1128                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1129       row  = cluster->GetLabel(0);
1130       col  = cluster->GetLabel(1);
1131       time = cluster->GetLabel(2);
1132
1133       for (iPad = 0; iPad < kNclus; iPad++) {
1134         Int_t iPadCol = col - 1 + iPad;
1135         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1136         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1137       }
1138
1139     }
1140
1141   }
1142
1143   // Copy the track indices into the cluster
1144   // Loop though the clusters found in this ROC
1145   for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1146
1147     AliTRDcluster *cluster = (AliTRDcluster *)
1148       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1149     cluster->SetLabel(-9999,0);
1150     cluster->SetLabel(-9999,1);
1151     cluster->SetLabel(-9999,2);
1152
1153     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1154
1155   }
1156
1157   delete [] idxTracks;
1158
1159   return kTRUE;
1160
1161 }
1162
1163 //_____________________________________________________________________________
1164 Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
1165 {
1166   //
1167   // Get COG position
1168   // Used for clusters with more than 3 pads - where LUT not applicable
1169   //
1170
1171   Double_t sum = signal[0]
1172               + signal[1]
1173               + signal[2]
1174               + signal[3]
1175               + signal[4];
1176
1177   // ???????????? CBL
1178   // Go to 3 pad COG ????
1179   // ???????????? CBL
1180   Double_t res = (0.0 * (-signal[0] + signal[4])
1181                       + (-signal[1] + signal[3])) / sum;
1182
1183   return res;
1184
1185 }
1186
1187 //_____________________________________________________________________________
1188 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal) const
1189 {
1190   //
1191   // Method to unfold neighbouring maxima.
1192   // The charge ratio on the overlapping pad is calculated
1193   // until there is no more change within the range given by eps.
1194   // The resulting ratio is then returned to the calling method.
1195   //
1196
1197   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1198   if (!calibration) {
1199     AliError("No AliTRDcalibDB instance available\n");
1200     return kFALSE;
1201   }
1202
1203   Int_t   irc                = 0;
1204   Int_t   itStep             = 0;                 // Count iteration steps
1205
1206   Double_t ratio             = 0.5;               // Start value for ratio
1207   Double_t prevRatio         = 0.0;               // Store previous ratio
1208
1209   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1210   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1211   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1212
1213   // Start the iteration
1214   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1215
1216     itStep++;
1217     prevRatio = ratio;
1218
1219     // Cluster position according to charge ratio
1220     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1221                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1222     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1223                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1224
1225     // Set cluster charge ratio
1226     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1227     Double_t ampLeft  = padSignal[1] / newSignal[1];
1228     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1229     Double_t ampRight = padSignal[3] / newSignal[1];
1230
1231     // Apply pad response to parameters
1232     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1233     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1234
1235     // Calculate new overlapping ratio
1236     ratio = TMath::Min((Double_t) 1.0
1237                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1238
1239   }
1240
1241   return ratio;
1242
1243 }
1244
1245 //_____________________________________________________________________________
1246 void AliTRDclusterizer::TailCancelation()
1247 {
1248   //
1249   // Applies the tail cancelation and gain factors:
1250   // Transform fDigitsIn to fDigitsOut
1251   //
1252
1253   Int_t iRow  = 0;
1254   Int_t iCol  = 0;
1255   Int_t iTime = 0;
1256
1257   Double_t *inADC  = new Double_t[fTimeTotal];  // ADC data before tail cancellation
1258   Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
1259   fIndexes->ResetCounters();
1260   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDReconstructor::kClusterizer);
1261   while(fIndexes->NextRCIndex(iRow, iCol))
1262     {
1263       Float_t  fCalGainFactorROCValue = fCalGainFactorROC->GetValue(iCol,iRow);
1264       Double_t gain                  = fCalGainFactorDetValue
1265                                     * fCalGainFactorROCValue;
1266
1267       Bool_t corrupted = kFALSE;
1268       for (iTime = 0; iTime < fTimeTotal; iTime++)
1269         {
1270           // Apply gain gain factor
1271           inADC[iTime]   = fDigitsIn->GetDataB(iRow,iCol,iTime);
1272           if (fDigitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
1273           inADC[iTime]  /= gain;
1274           outADC[iTime]  = inADC[iTime];
1275           if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 7){
1276             (*fDebugStream) << "TailCancellation"
1277                               << "col="  << iCol
1278                               << "row="  << iRow
1279                               << "time=" << iTime
1280                               << "inADC=" << inADC[iTime]
1281                               << "gain=" << gain
1282                               << "outADC=" << outADC[iTime]
1283                               << "corrupted=" << corrupted
1284                               << "\n";
1285           }
1286         }
1287       if (!corrupted)
1288         {
1289           // Apply the tail cancelation via the digital filter
1290           // (only for non-coorupted pads)
1291           if (fReconstructor->GetRecoParam()->IsTailCancelation())
1292             DeConvExp(inADC,outADC,fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1293         }
1294
1295       for(iTime = 0; iTime < fTimeTotal; iTime++)//while (fIndexes->NextTbinIndex(iTime))
1296         {
1297           // Store the amplitude of the digit if above threshold
1298           if (outADC[iTime] > fADCthresh)
1299             fDigitsOut->SetData(iRow,iCol,iTime,outADC[iTime]);
1300         } // while itime
1301
1302     } // while irow icol
1303
1304   delete [] inADC;
1305   delete [] outADC;
1306
1307   return;
1308
1309 }
1310
1311 //_____________________________________________________________________________
1312 void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1313                                   ,const Int_t n, const Int_t nexp)
1314 {
1315   //
1316   // Tail cancellation by deconvolution for PASA v4 TRF
1317   //
1318
1319   Double_t rates[2];
1320   Double_t coefficients[2];
1321
1322   // Initialization (coefficient = alpha, rates = lambda)
1323   Double_t r1 = 1.0;
1324   Double_t r2 = 1.0;
1325   Double_t c1 = 0.5;
1326   Double_t c2 = 0.5;
1327
1328   if (nexp == 1) {   // 1 Exponentials
1329     r1 = 1.156;
1330     r2 = 0.130;
1331     c1 = 0.066;
1332     c2 = 0.000;
1333   }
1334   if (nexp == 2) {   // 2 Exponentials
1335     Double_t par[4];
1336     fReconstructor->GetTCParams(par);
1337     r1 = par[0];//1.156;
1338     r2 = par[1];//0.130;
1339     c1 = par[2];//0.114;
1340     c2 = par[3];//0.624;
1341   }
1342
1343   coefficients[0] = c1;
1344   coefficients[1] = c2;
1345
1346   Double_t dt = 0.1;
1347
1348   rates[0] = TMath::Exp(-dt/(r1));
1349   rates[1] = TMath::Exp(-dt/(r2));
1350
1351   Int_t i = 0;
1352   Int_t k = 0;
1353
1354   Double_t reminder[2];
1355   Double_t correction = 0.0;
1356   Double_t result     = 0.0;
1357
1358   // Attention: computation order is important
1359   for (k = 0; k < nexp; k++) {
1360     reminder[k] = 0.0;
1361   }
1362
1363   for (i = 0; i < n; i++) {
1364
1365     result    = (source[i] - correction);    // No rescaling
1366     target[i] = result;
1367
1368     for (k = 0; k < nexp; k++) {
1369       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1370     }
1371
1372     correction = 0.0;
1373     for (k = 0; k < nexp; k++) {
1374       correction += reminder[k];
1375     }
1376
1377   }
1378
1379 }
1380
1381 //_____________________________________________________________________________
1382 void AliTRDclusterizer::ResetRecPoints()
1383 {
1384   //
1385   // Resets the list of rec points
1386   //
1387
1388   if (fRecPoints) {
1389     fRecPoints->Delete();
1390     delete fRecPoints;
1391   }
1392 }
1393
1394 //_____________________________________________________________________________
1395 TClonesArray *AliTRDclusterizer::RecPoints()
1396 {
1397   //
1398   // Returns the list of rec points
1399   //
1400
1401   if (!fRecPoints) {
1402     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1403       // determine number of clusters which has to be allocated
1404       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1405       if(fReconstructor->IsHLT()) nclusters /= AliTRDgeometry::kNsector;
1406
1407       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1408     }
1409     //SetClustersOwner(kTRUE);
1410     AliTRDReconstructor::SetClusters(0x0);
1411   }
1412   return fRecPoints;
1413
1414 }
1415
1416 //_____________________________________________________________________________
1417 void AliTRDclusterizer::FillLUT()
1418 {
1419   //
1420   // Create the LUT
1421   //
1422
1423   const Int_t kNlut = 128;
1424
1425   fLUTbin = AliTRDgeometry::kNlayer * kNlut;
1426
1427   // The lookup table from Bogdan
1428   Float_t lut[AliTRDgeometry::kNlayer][kNlut] = {
1429     {
1430       0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611,
1431       0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187,
1432       0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704,
1433       0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160,
1434       0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557,
1435       0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901,
1436       0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207,
1437       0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483,
1438       0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727,
1439       0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952,
1440       0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157,
1441       0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345,
1442       0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520,
1443       0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683,
1444       0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824,
1445       0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978
1446     },
1447     {
1448       0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642,
1449       0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241,
1450       0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770,
1451       0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229,
1452       0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627,
1453       0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968,
1454       0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266,
1455       0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536,
1456       0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772,
1457       0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991,
1458       0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187,
1459       0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369,
1460       0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538,
1461       0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694,
1462       0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832,
1463       0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1464     },
1465     {
1466       0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674,
1467       0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296,
1468       0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839,
1469       0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299,
1470       0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693,
1471       0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033,
1472       0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323,
1473       0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586,
1474       0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815,
1475       0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027,
1476       0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217,
1477       0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393,
1478       0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555,
1479       0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705,
1480       0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839,
1481       0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1482     },
1483     {
1484       0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708,
1485       0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356,
1486       0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910,
1487       0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371,
1488       0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759,
1489       0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095,
1490       0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378,
1491       0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633,
1492       0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857,
1493       0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061,
1494       0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244,
1495       0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414,
1496       0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571,
1497       0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715,
1498       0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846,
1499       0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1500     },
1501     {
1502       0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744,
1503       0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419,
1504       0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984,
1505       0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444,
1506       0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824,
1507       0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152,
1508       0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432,
1509       0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677,
1510       0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896,
1511       0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094,
1512       0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270,
1513       0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435,
1514       0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586,
1515       0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725,
1516       0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852,
1517       0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1518     },
1519     {
1520       0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792,
1521       0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505,
1522       0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081,
1523       0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541,
1524       0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917,
1525       0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235,
1526       0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511,
1527       0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748,
1528       0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962,
1529       0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152,
1530       0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324,
1531       0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483,
1532       0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629,
1533       0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759,
1534       0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859,
1535       0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1536     }
1537   };
1538
1539   if (fLUT) {
1540     delete [] fLUT;
1541   }
1542   fLUT = new Double_t[fLUTbin];
1543
1544   for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
1545     for (Int_t ilut  = 0; ilut  <  kNlut; ilut++  ) {
1546       fLUT[ilayer*kNlut+ilut] = lut[ilayer][ilut];
1547     }
1548   }
1549
1550 }
1551
1552 //_____________________________________________________________________________
1553 Double_t AliTRDclusterizer::LUTposition(Int_t ilayer
1554                                       , Double_t ampL
1555                                       , Double_t ampC
1556                                       , Double_t ampR) const
1557 {
1558   //
1559   // Calculates the cluster position using the lookup table.
1560   // Method provided by Bogdan Vulpescu.
1561   //
1562
1563   const Int_t kNlut = 128;
1564
1565   Double_t pos;
1566   Double_t x    = 0.0;
1567   Double_t xmin;
1568   Double_t xmax;
1569   Double_t xwid;
1570
1571   Int_t    side = 0;
1572   Int_t    ix;
1573
1574   Double_t xMin[AliTRDgeometry::kNlayer] = { 0.006492, 0.006377, 0.006258
1575                                           , 0.006144, 0.006030, 0.005980 };
1576   Double_t xMax[AliTRDgeometry::kNlayer] = { 0.960351, 0.965870, 0.970445
1577                                           , 0.974352, 0.977667, 0.996101 };
1578
1579   if      (ampL > ampR) {
1580     x    = (ampL - ampR) / ampC;
1581     side = -1;
1582   }
1583   else if (ampL < ampR) {
1584     x    = (ampR - ampL) / ampC;
1585     side = +1;
1586   }
1587
1588   if (ampL != ampR) {
1589
1590     xmin = xMin[ilayer] + 0.000005;
1591     xmax = xMax[ilayer] - 0.000005;
1592     xwid = (xmax - xmin) / 127.0;
1593
1594     if      (x < xmin) {
1595       pos = 0.0000;
1596     }
1597     else if (x > xmax) {
1598       pos = side * 0.5000;
1599     }
1600     else {
1601       ix  = (Int_t) ((x - xmin) / xwid);
1602       pos = side * fLUT[ilayer*kNlut+ix];
1603     }
1604
1605   }
1606   else {
1607
1608     pos = 0.0;
1609
1610   }
1611
1612   return pos;
1613
1614 }