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 <TClonesArray.h>
  25 #include <TObjArray.h>
  26
  27 #include "AliRunLoader.h"
  28 #include "AliLoader.h"
  29 #include "AliAlignObj.h"
  30
  31 #include "AliTRDclusterizer.h"
  32 #include "AliTRDcluster.h"
  33 #include "AliTRDReconstructor.h"
  34 #include "AliTRDgeometry.h"
  35 #include "AliTRDarrayDictionary.h"
  36 #include "AliTRDarrayADC.h"
  37 #include "AliTRDdigitsManager.h"
  38 #include "AliTRDdigitsParam.h"
  39 #include "AliTRDrawData.h"
  40 #include "AliTRDcalibDB.h"
  41 #include "AliTRDtransform.h"
  42 #include "AliTRDSignalIndex.h"
  43 #include "AliTRDrawStreamBase.h"
  44 #include "AliTRDfeeParam.h"
  45 #include "AliTRDtrackletWord.h"
  46
  47 #include "TTreeStream.h"
  48
  49 #include "Cal/AliTRDCalROC.h"
  50 #include "Cal/AliTRDCalDet.h"
  51 #include "Cal/AliTRDCalSingleChamberStatus.h"
  52
  53 ClassImp(AliTRDclusterizer)
  54
  55 //_____________________________________________________________________________
  56 AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
  57   :TNamed()
  58   ,fReconstructor(rec)
  59   ,fRunLoader(NULL)
  60   ,fClusterTree(NULL)
  61   ,fRecPoints(NULL)
  62   ,fTracklets(NULL)
  63   ,fTrackletTree(NULL)
  64   ,fDigitsManager(new AliTRDdigitsManager())
  65   ,fTrackletContainer(NULL)
  66   ,fRawVersion(2)
  67   ,fTransform(new AliTRDtransform(0))
  68   ,fDigits(NULL)
  69   ,fIndexes(NULL)
  70   ,fMaxThresh(0)
  71   ,fSigThresh(0)
  72   ,fMinMaxCutSigma(0)
  73   ,fMinLeftRightCutSigma(0)
  74   ,fLayer(0)
  75   ,fDet(0)
  76   ,fVolid(0)
  77   ,fColMax(0)
  78   ,fTimeTotal(0)
  79   ,fCalGainFactorROC(NULL)
  80   ,fCalGainFactorDetValue(0)
  81   ,fCalNoiseROC(NULL)
  82   ,fCalNoiseDetValue(0)
  83   ,fCalPadStatusROC(NULL)
  84   ,fClusterROC(0)
  85   ,firstClusterROC(0)
  86   ,fNoOfClusters(0)
  87   ,fBaseline(0)
  88   ,fRawStream(NULL)
  89 {
  90   //
  91   // AliTRDclusterizer default constructor
  92   //
  93
  94   SetBit(kLabels, kTRUE);
  95   SetBit(knewDM, kFALSE);
  96
  97   AliTRDcalibDB *trd = 0x0;
  98   if (!(trd = AliTRDcalibDB::Instance())) {
  99     AliFatal("Could not get calibration object");
 100   }
 101
 102   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
 103
 104   // Initialize debug stream
 105   if(fReconstructor){
 106     if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 1){
 107       TDirectory *savedir = gDirectory;
 108       //fgGetDebugStream    = new TTreeSRedirector("TRD.ClusterizerDebug.root");
 109       savedir->cd();
 110     }
 111   }
 112
 113 }
 114
 115 //_____________________________________________________________________________
 116 AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, const AliTRDReconstructor *const rec)
 117   :TNamed(name,title)
 118   ,fReconstructor(rec)
 119   ,fRunLoader(NULL)
 120   ,fClusterTree(NULL)
 121   ,fRecPoints(NULL)
 122   ,fTracklets(NULL)
 123   ,fTrackletTree(NULL)
 124   ,fDigitsManager(new AliTRDdigitsManager())
 125   ,fTrackletContainer(NULL)
 126   ,fRawVersion(2)
 127   ,fTransform(new AliTRDtransform(0))
 128   ,fDigits(NULL)
 129   ,fIndexes(NULL)
 130   ,fMaxThresh(0)
 131   ,fSigThresh(0)
 132   ,fMinMaxCutSigma(0)
 133   ,fMinLeftRightCutSigma(0)
 134   ,fLayer(0)
 135   ,fDet(0)
 136   ,fVolid(0)
 137   ,fColMax(0)
 138   ,fTimeTotal(0)
 139   ,fCalGainFactorROC(NULL)
 140   ,fCalGainFactorDetValue(0)
 141   ,fCalNoiseROC(NULL)
 142   ,fCalNoiseDetValue(0)
 143   ,fCalPadStatusROC(NULL)
 144   ,fClusterROC(0)
 145   ,firstClusterROC(0)
 146   ,fNoOfClusters(0)
 147   ,fBaseline(0)
 148   ,fRawStream(NULL)
 149 {
 150   //
 151   // AliTRDclusterizer constructor
 152   //
 153
 154   SetBit(kLabels, kTRUE);
 155   SetBit(knewDM, kFALSE);
 156
 157   AliTRDcalibDB *trd = 0x0;
 158   if (!(trd = AliTRDcalibDB::Instance())) {
 159     AliFatal("Could not get calibration object");
 160   }
 161
 162   fDigitsManager->CreateArrays();
 163
 164   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
 165
 166   //FillLUT();
 167
 168 }
 169
 170 //_____________________________________________________________________________
 171 AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
 172   :TNamed(c)
 173   ,fReconstructor(c.fReconstructor)
 174   ,fRunLoader(NULL)
 175   ,fClusterTree(NULL)
 176   ,fRecPoints(NULL)
 177   ,fTracklets(NULL)
 178   ,fTrackletTree(NULL)
 179   ,fDigitsManager(NULL)
 180   ,fTrackletContainer(NULL)
 181   ,fRawVersion(2)
 182   ,fTransform(NULL)
 183   ,fDigits(NULL)
 184   ,fIndexes(NULL)
 185   ,fMaxThresh(0)
 186   ,fSigThresh(0)
 187   ,fMinMaxCutSigma(0)
 188   ,fMinLeftRightCutSigma(0)
 189   ,fLayer(0)
 190   ,fDet(0)
 191   ,fVolid(0)
 192   ,fColMax(0)
 193   ,fTimeTotal(0)
 194   ,fCalGainFactorROC(NULL)
 195   ,fCalGainFactorDetValue(0)
 196   ,fCalNoiseROC(NULL)
 197   ,fCalNoiseDetValue(0)
 198   ,fCalPadStatusROC(NULL)
 199   ,fClusterROC(0)
 200   ,firstClusterROC(0)
 201   ,fNoOfClusters(0)
 202   ,fBaseline(0)
 203   ,fRawStream(NULL)
 204 {
 205   //
 206   // AliTRDclusterizer copy constructor
 207   //
 208
 209   SetBit(kLabels, kTRUE);
 210   SetBit(knewDM, kFALSE);
 211
 212   //FillLUT();
 213
 214 }
 215
 216 //_____________________________________________________________________________
 217 AliTRDclusterizer::~AliTRDclusterizer()
 218 {
 219   //
 220   // AliTRDclusterizer destructor
 221   //
 222
 223   if (fRecPoints/* && IsClustersOwner()*/){
 224     fRecPoints->Delete();
 225     delete fRecPoints;
 226   }
 227
 228   if (fTracklets){
 229     fTracklets->Delete();
 230     delete fTracklets;
 231   }
 232
 233   if (fDigitsManager) {
 234     delete fDigitsManager;
 235     fDigitsManager = NULL;
 236   }
 237
 238   if (fTrackletContainer){
 239     delete [] fTrackletContainer[0];
 240     delete [] fTrackletContainer[1];
 241     delete [] fTrackletContainer;
 242     fTrackletContainer = NULL;
 243   }
 244
 245   if (fTransform){
 246     delete fTransform;
 247     fTransform = NULL;
 248   }
 249
 250   if (fRawStream){
 251     delete fRawStream;
 252     fRawStream = NULL;
 253   }
 254
 255 }
 256
 257 //_____________________________________________________________________________
 258 AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
 259 {
 260   //
 261   // Assignment operator
 262   //
 263
 264   if (this != &c)
 265     {
 266       ((AliTRDclusterizer &) c).Copy(*this);
 267     }
 268
 269   return *this;
 270
 271 }
 272
 273 //_____________________________________________________________________________
 274 void AliTRDclusterizer::Copy(TObject &c) const
 275 {
 276   //
 277   // Copy function
 278   //
 279
 280   ((AliTRDclusterizer &) c).fClusterTree   = NULL;
 281   ((AliTRDclusterizer &) c).fRecPoints     = NULL;
 282   ((AliTRDclusterizer &) c).fTrackletTree  = NULL;
 283   ((AliTRDclusterizer &) c).fDigitsManager = NULL;
 284   ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
 285   ((AliTRDclusterizer &) c).fRawVersion    = fRawVersion;
 286   ((AliTRDclusterizer &) c).fTransform     = NULL;
 287   ((AliTRDclusterizer &) c).fDigits      = NULL;
 288   ((AliTRDclusterizer &) c).fIndexes       = NULL;
 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).fCalPadStatusROC = NULL;
 303   ((AliTRDclusterizer &) c).fClusterROC    = 0;
 304   ((AliTRDclusterizer &) c).firstClusterROC= 0;
 305   ((AliTRDclusterizer &) c).fNoOfClusters  = 0;
 306   ((AliTRDclusterizer &) c).fBaseline      = 0;
 307   ((AliTRDclusterizer &) c).fRawStream     = NULL;
 308
 309 }
 310
 311 //_____________________________________________________________________________
 312 Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
 313 {
 314   //
 315   // Opens the AliROOT file. Output and input are in the same file
 316   //
 317
 318   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 319   fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
 320
 321   if (!fRunLoader) {
 322     fRunLoader = AliRunLoader::Open(name);
 323   }
 324
 325   if (!fRunLoader) {
 326     AliError(Form("Can not open session for file %s.",name));
 327     return kFALSE;
 328   }
 329
 330   OpenInput(nEvent);
 331   OpenOutput();
 332
 333   return kTRUE;
 334
 335 }
 336
 337 //_____________________________________________________________________________
 338 Bool_t AliTRDclusterizer::OpenOutput()
 339 {
 340   //
 341   // Open the output file
 342   //
 343
 344   if (!fReconstructor->IsWritingClusters()) return kTRUE;
 345
 346   TObjArray *ioArray = 0x0;
 347
 348   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 349   loader->MakeTree("R");
 350
 351   fClusterTree = loader->TreeR();
 352   fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
 353
 354   return kTRUE;
 355
 356 }
 357
 358 //_____________________________________________________________________________
 359 Bool_t AliTRDclusterizer::OpenOutput(TTree *const clusterTree)
 360 {
 361   //
 362   // Connect the output tree
 363   //
 364
 365   // clusters writing
 366   if (fReconstructor->IsWritingClusters()){
 367     TObjArray *ioArray = 0x0;
 368     fClusterTree = clusterTree;
 369     fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
 370   }
 371   return kTRUE;
 372 }
 373
 374 //_____________________________________________________________________________
 375 Bool_t AliTRDclusterizer::OpenTrackletOutput()
 376 {
 377   //
 378   // Tracklet writing
 379   //
 380
 381   if (fReconstructor->IsWritingTracklets()){
 382     TString evfoldname = AliConfig::GetDefaultEventFolderName();
 383     fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
 384
 385     if (!fRunLoader) {
 386       fRunLoader = AliRunLoader::Open("galice.root");
 387     }
 388     if (!fRunLoader) {
 389       AliError(Form("Can not open session for file galice.root."));
 390       return kFALSE;
 391     }
 392
 393     UInt_t **leaves = new UInt_t *[2];
 394     AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 395     if (!dl) {
 396       AliError("Could not get the tracklets data loader!");
 397       dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
 398       fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl);
 399     }
 400     fTrackletTree = dl->Tree();
 401     if (!fTrackletTree)
 402       {
 403        dl->MakeTree();
 404        fTrackletTree = dl->Tree();
 405       }
 406     TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
 407     if (!trkbranch)
 408       fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
 409   }
 410
 411   return kTRUE;
 412 }
 413
 414 //_____________________________________________________________________________
 415 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
 416 {
 417   //
 418   // Opens a ROOT-file with TRD-hits and reads in the digits-tree
 419   //
 420
 421   // Import the Trees for the event nEvent in the file
 422   fRunLoader->GetEvent(nEvent);
 423
 424   return kTRUE;
 425
 426 }
 427
 428 //_____________________________________________________________________________
 429 Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
 430 {
 431   //
 432   // Fills TRDcluster branch in the tree with the clusters
 433   // found in detector = det. For det=-1 writes the tree.
 434   //
 435
 436   if ((det <                      -1) ||
 437       (det >= AliTRDgeometry::Ndet())) {
 438     AliError(Form("Unexpected detector index %d.\n",det));
 439     return kFALSE;
 440   }
 441
 442   TObjArray *ioArray = new TObjArray(400);
 443   TBranch *branch = fClusterTree->GetBranch("TRDcluster");
 444   if (!branch) {
 445     branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
 446   } else branch->SetAddress(&ioArray);
 447
 448   Int_t nRecPoints = RecPoints()->GetEntriesFast();
 449   if(det >= 0){
 450     for (Int_t i = 0; i < nRecPoints; i++) {
 451       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
 452       if(det != c->GetDetector()) continue;
 453       ioArray->AddLast(c);
 454     }
 455     fClusterTree->Fill();
 456     ioArray->Clear();
 457   } else {
 458     Int_t detOld = -1, nw(0);
 459     for (Int_t i = 0; i < nRecPoints; i++) {
 460       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
 461       if(c->GetDetector() != detOld){
 462         nw += ioArray->GetEntriesFast();
 463         fClusterTree->Fill();
 464         ioArray->Clear();
 465         detOld = c->GetDetector();
 466       }
 467       ioArray->AddLast(c);
 468     }
 469     if(ioArray->GetEntriesFast()){
 470       nw += ioArray->GetEntriesFast();
 471       fClusterTree->Fill();
 472       ioArray->Clear();
 473     }
 474     AliDebug(2, Form("Clusters FOUND[%d] WRITTEN[%d] STATUS[%s]", nRecPoints, nw, nw==nRecPoints?"OK":"FAILED"));
 475   }
 476   delete ioArray;
 477
 478   return kTRUE;
 479 }
 480
 481 //_____________________________________________________________________________
 482 Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
 483 {
 484   //
 485   // Write the raw data tracklets into seperate file
 486   //
 487
 488   UInt_t **leaves = new UInt_t *[2];
 489   for (Int_t i=0; i<2 ;i++){
 490     leaves[i] = new UInt_t[258];
 491     leaves[i][0] = det; // det
 492     leaves[i][1] = i;   // side
 493     memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
 494   }
 495
 496   if (!fTrackletTree){
 497     AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 498     dl->MakeTree();
 499     fTrackletTree = dl->Tree();
 500   }
 501
 502   TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
 503   if (!trkbranch) {
 504     trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
 505   }
 506
 507   for (Int_t i=0; i<2; i++){
 508     if (leaves[i][2]>0) {
 509       trkbranch->SetAddress(leaves[i]);
 510       fTrackletTree->Fill();
 511     }
 512   }
 513
 514   AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 515   dl->WriteData("OVERWRITE");
 516   //dl->Unload();
 517   delete [] leaves;
 518
 519   return kTRUE;
 520
 521 }
 522
 523 //_____________________________________________________________________________
 524 Bool_t AliTRDclusterizer::ReadDigits()
 525 {
 526   //
 527   // Reads the digits arrays from the input aliroot file
 528   //
 529
 530   if (!fRunLoader) {
 531     AliError("No run loader available");
 532     return kFALSE;
 533   }
 534
 535   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 536   if (!loader->TreeD()) {
 537     loader->LoadDigits();
 538   }
 539
 540   // Read in the digit arrays
 541   return (fDigitsManager->ReadDigits(loader->TreeD()));
 542
 543 }
 544
 545 //_____________________________________________________________________________
 546 Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
 547 {
 548   //
 549   // Reads the digits arrays from the input tree
 550   //
 551
 552   // Read in the digit arrays
 553   return (fDigitsManager->ReadDigits(digitsTree));
 554
 555 }
 556
 557 //_____________________________________________________________________________
 558 Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
 559 {
 560   //
 561   // Reads the digits arrays from the ddl file
 562   //
 563
 564   AliTRDrawData raw;
 565   fDigitsManager = raw.Raw2Digits(rawReader);
 566
 567   return kTRUE;
 568
 569 }
 570
 571 //_____________________________________________________________________________
 572 Bool_t AliTRDclusterizer::MakeClusters()
 573 {
 574   //
 575   // Creates clusters from digits
 576   //
 577
 578   // Propagate info from the digits manager
 579   if (TestBit(kLabels)){
 580     SetBit(kLabels, fDigitsManager->UsesDictionaries());
 581   }
 582
 583   Bool_t fReturn = kTRUE;
 584   for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
 585
 586     AliTRDarrayADC *digitsIn = (AliTRDarrayADC*) fDigitsManager->GetDigits(i); //mod
 587     // This is to take care of switched off super modules
 588     if (!digitsIn->HasData()) continue;
 589     digitsIn->Expand();
 590     digitsIn->DeleteNegatives();  // Restore digits array to >=0 values
 591     AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
 592     if (indexes->IsAllocated() == kFALSE){
 593       fDigitsManager->BuildIndexes(i);
 594     }
 595
 596     Bool_t fR = kFALSE;
 597     if (indexes->HasEntry()){
 598       if (TestBit(kLabels)){
 599         for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
 600           AliTRDarrayDictionary *tracksIn = 0; //mod
 601           tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict);  //mod
 602           tracksIn->Expand();
 603         }
 604       }
 605       fR = MakeClusters(i);
 606       fReturn = fR && fReturn;
 607     }
 608
 609     //if (fR == kFALSE){
 610     //  if(IsWritingClusters()) WriteClusters(i);
 611     //  ResetRecPoints();
 612     //}
 613
 614     // No compress just remove
 615     fDigitsManager->RemoveDigits(i);
 616     fDigitsManager->RemoveDictionaries(i);
 617     fDigitsManager->ClearIndexes(i);
 618   }
 619   fReconstructor->SetDigitsParam(fDigitsManager->GetDigitsParam());
 620
 621   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 622
 623   AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
 624
 625   return fReturn;
 626
 627 }
 628
 629 //_____________________________________________________________________________
 630 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
 631 {
 632   //
 633   // Creates clusters from raw data
 634   //
 635
 636   return Raw2ClustersChamber(rawReader);
 637
 638 }
 639
 640 //_____________________________________________________________________________
 641 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
 642 {
 643   //
 644   // Creates clusters from raw data
 645   //
 646
 647   // Create the digits manager
 648   if (!fDigitsManager){
 649     SetBit(knewDM, kTRUE);
 650     fDigitsManager = new AliTRDdigitsManager(kTRUE);
 651     fDigitsManager->CreateArrays();
 652   }
 653
 654   fDigitsManager->SetUseDictionaries(TestBit(kLabels));
 655
 656   // tracklet container for raw tracklet writing
 657   if (!fTrackletContainer && ( fReconstructor->IsWritingTracklets() || fReconstructor->IsProcessingTracklets() )) {
 658     // maximum tracklets for one HC
 659     const Int_t kTrackletChmb=256;
 660     fTrackletContainer = new UInt_t *[2];
 661     fTrackletContainer[0] = new UInt_t[kTrackletChmb];
 662     fTrackletContainer[1] = new UInt_t[kTrackletChmb];
 663   }
 664
 665   if(!fRawStream)
 666     fRawStream = AliTRDrawStreamBase::GetRawStream(rawReader);
 667   else
 668     fRawStream->SetReader(rawReader);
 669
 670   if(fReconstructor->IsHLT()){
 671     fRawStream->SetSharedPadReadout(kFALSE);
 672     fRawStream->SetNoErrorWarning();
 673   }
 674
 675   AliDebug(1,Form("Stream version: %s", fRawStream->IsA()->GetName()));
 676
 677   Int_t det    = 0;
 678   while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
 679     Bool_t iclusterBranch = kFALSE;
 680     if (fDigitsManager->GetIndexes(det)->HasEntry()){
 681       iclusterBranch = MakeClusters(det);
 682     }
 683
 684     fDigitsManager->ClearArrays(det);
 685
 686     if (!fReconstructor->IsWritingTracklets()) continue;
 687     if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
 688   }
 689   fReconstructor->SetDigitsParam(fDigitsManager->GetDigitsParam());
 690
 691   if (fTrackletContainer){
 692     delete [] fTrackletContainer[0];
 693     delete [] fTrackletContainer[1];
 694     delete [] fTrackletContainer;
 695     fTrackletContainer = NULL;
 696   }
 697
 698   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 699
 700   if(!TestBit(knewDM)){
 701     delete fDigitsManager;
 702     fDigitsManager = NULL;
 703   }
 704
 705   AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
 706   return kTRUE;
 707
 708 }
 709
 710 //_____________________________________________________________________________
 711 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
 712 {
 713   //
 714   // Check if a pad is masked
 715   //
 716
 717   UChar_t status = 0;
 718
 719   if(signal>0 && TESTBIT(signal, 10)){
 720     CLRBIT(signal, 10);
 721     for(int ibit=0; ibit<4; ibit++){
 722       if(TESTBIT(signal, 11+ibit)){
 723         SETBIT(status, ibit);
 724         CLRBIT(signal, 11+ibit);
 725       }
 726     }
 727   }
 728   return status;
 729 }
 730
 731 //_____________________________________________________________________________
 732 void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
 733   //
 734   // Set the pad status into out
 735   // First three bits are needed for the position encoding
 736   //
 737   out |= status << 3;
 738 }
 739
 740 //_____________________________________________________________________________
 741 UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
 742   //
 743   // return the staus encoding of the corrupted pad
 744   //
 745   return static_cast<UChar_t>(encoding >> 3);
 746 }
 747
 748 //_____________________________________________________________________________
 749 Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
 750   //
 751   // Return the position of the corruption
 752   //
 753   return encoding & 7;
 754 }
 755
 756 //_____________________________________________________________________________
 757 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 758 {
 759   //
 760   // Generates the cluster.
 761   //
 762
 763   // Get the digits
 764   fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
 765   fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline();
 766
 767   // This is to take care of switched off super modules
 768   if (!fDigits->HasData()) return kFALSE;
 769
 770   fIndexes = fDigitsManager->GetIndexes(det);
 771   if (fIndexes->IsAllocated() == kFALSE) {
 772     AliError("Indexes do not exist!");
 773     return kFALSE;
 774   }
 775
 776   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
 777   if (!calibration) {
 778     AliFatal("No AliTRDcalibDB instance available\n");
 779     return kFALSE;
 780   }
 781
 782   if (!fReconstructor){
 783     AliError("Reconstructor not set\n");
 784     return kFALSE;
 785   }
 786
 787
 788   fMaxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
 789   fSigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
 790   fMinMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
 791   fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
 792
 793   Int_t istack  = fIndexes->GetStack();
 794   fLayer  = fIndexes->GetLayer();
 795   Int_t isector = fIndexes->GetSM();
 796
 797   // Start clustering in the chamber
 798
 799   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
 800   if (fDet != det) {
 801     AliError("Strange Detector number Missmatch!");
 802     return kFALSE;
 803   }
 804
 805   AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
 806
 807   // TRD space point transformation
 808   fTransform->SetDetector(det);
 809
 810   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
 811   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
 812   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
 813
 814   if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
 815     AddTrackletsToArray();
 816
 817   fColMax    = fDigits->GetNcol();
 818   //Int_t nRowMax    = fDigits->GetNrow();
 819   fTimeTotal = fDigits->GetNtime();
 820
 821   // Check consistency between OCDB and raw data
 822   Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS();
 823   if ((nTimeOCDB  >         -1) &&
 824       (fTimeTotal != nTimeOCDB)) {
 825     AliError(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d])"
 826                 ,fTimeTotal,calibration->GetNumberOfTimeBinsDCS()));
 827   }
 828
 829   // Detector wise calibration object for the gain factors
 830   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 831   // Calibration object with pad wise values for the gain factors
 832   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
 833   // Calibration value for chamber wise gain factor
 834   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
 835
 836   // Detector wise calibration object for the noise
 837   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
 838   // Calibration object with pad wise values for the noise
 839   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
 840   // Calibration value for chamber wise noise
 841   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
 842
 843   // Calibration object with the pad status
 844   fCalPadStatusROC       = calibration->GetPadStatusROC(fDet);
 845
 846   SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT());
 847   SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS());
 848   SetBit(kHLT, fReconstructor->IsHLT());
 849
 850   firstClusterROC = -1;
 851   fClusterROC     =  0;
 852
 853   // Apply the gain and the tail cancelation via digital filter
 854   if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation();
 855
 856   MaxStruct curr, last;
 857   Int_t nMaximas = 0, nCorrupted = 0;
 858
 859   // Here the clusterfining is happening
 860
 861   for(curr.time = 0; curr.time < fTimeTotal; curr.time++){
 862     while(fIndexes->NextRCIndex(curr.row, curr.col)){
 863       if(IsMaximum(curr, curr.padStatus, &curr.signals[0])){
 864         if(last.row>-1){
 865           if(curr.time==last.time && curr.row==last.row && curr.col==last.col+2) FivePadCluster(last, curr);
 866           CreateCluster(last);
 867         }
 868         last=curr; curr.fivePad=kFALSE;
 869       }
 870     }
 871   }
 872   if(last.row>-1) CreateCluster(last);
 873
 874   if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
 875     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
 876     (*fDebugStream) << "MakeClusters"
 877       << "Detector="   << det
 878       << "NMaxima="    << nMaximas
 879       << "NClusters="  << fClusterROC
 880       << "NCorrupted=" << nCorrupted
 881       << "\n";
 882   }
 883   if (TestBit(kLabels)) AddLabels();
 884
 885   return kTRUE;
 886
 887 }
 888
 889 //_____________________________________________________________________________
 890 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
 891 {
 892   //
 893   // Returns true if this row,col,time combination is a maximum.
 894   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
 895   //
 896
 897   Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
 898   Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
 899   if(Signals[1] < fMaxThresh) return kFALSE;
 900
 901   Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
 902   if (Signals[1] < noiseMiddleThresh) return kFALSE;
 903
 904   if (Max.col + 1 >= fColMax || Max.col < 1) return kFALSE;
 905
 906   UChar_t status[3]={
 907     fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
 908    ,fCalPadStatusROC->GetStatus(Max.col,   Max.row)
 909    ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
 910   };
 911
 912   gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
 913   Signals[0] = (Short_t)((fDigits->GetData(Max.row, Max.col-1, Max.time) - fBaseline) / gain + 0.5f);
 914   gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
 915   Signals[2] = (Short_t)((fDigits->GetData(Max.row, Max.col+1, Max.time) - fBaseline) / gain + 0.5f);
 916
 917   if(!(status[0] | status[1] | status[2])) {//all pads are good
 918     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
 919       if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
 920         if(Signals[0]<0)Signals[0]=0;
 921         if(Signals[2]<0)Signals[2]=0;
 922         Float_t  noiseSumThresh = fMinLeftRightCutSigma
 923           * fCalNoiseDetValue
 924           * fCalNoiseROC->GetValue(Max.col, Max.row);
 925         if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
 926         padStatus = 0;
 927         return kTRUE;
 928       }
 929     }
 930   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
 931     if(Signals[0]<0)Signals[0]=0;
 932     if(Signals[2]<0)Signals[2]=0;
 933     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) {
 934       Signals[2]=0;
 935       SetPadStatus(status[2], padStatus);
 936       return kTRUE;
 937     }
 938     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
 939       Signals[0]=0;
 940       SetPadStatus(status[0], padStatus);
 941       return kTRUE;
 942     }
 943     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
 944       Signals[1] = (Short_t)(fMaxThresh + 0.5f);
 945       SetPadStatus(status[1], padStatus);
 946       return kTRUE;
 947     }
 948   }
 949   return kFALSE;
 950 }
 951
 952 //_____________________________________________________________________________
 953 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
 954 {
 955   //
 956   // Look for 5 pad cluster with minimum in the middle
 957   // Gives back the ratio
 958   //
 959
 960   if (ThisMax.col >= fColMax - 3) return kFALSE;
 961   Float_t gain;
 962   if (ThisMax.col < fColMax - 5){
 963     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
 964     if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
 965       return kFALSE;
 966   }
 967   if (ThisMax.col > 1) {
 968     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
 969     if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
 970       return kFALSE;
 971   }
 972
 973   const Float_t kEpsilon = 0.01;
 974   Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
 975       NeighbourMax.signals[1], NeighbourMax.signals[2]};
 976
 977   // Unfold the two maxima and set the signal on
 978   // the overlapping pad to the ratio
 979   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
 980   ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
 981   NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
 982   ThisMax.fivePad=kTRUE;
 983   NeighbourMax.fivePad=kTRUE;
 984   return kTRUE;
 985
 986 }
 987
 988 //_____________________________________________________________________________
 989 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
 990 {
 991   //
 992   // Creates a cluster at the given position and saves it in fRecPoints
 993   //
 994
 995   Int_t nPadCount = 1;
 996   Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
 997   if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
 998
 999   AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
1000   cluster.SetNPads(nPadCount);
1001   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1002   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1003   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
1004
1005   cluster.SetFivePad(Max.fivePad);
1006   // set pads status for the cluster
1007   UChar_t maskPosition = GetCorruption(Max.padStatus);
1008   if (maskPosition) {
1009     cluster.SetPadMaskedPosition(maskPosition);
1010     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1011   }
1012
1013   // Transform the local cluster coordinates into calibrated
1014   // space point positions defined in the local tracking system.
1015   // Here the calibration for T0, Vdrift and ExB is applied as well.
1016   if(!fTransform->Transform(&cluster)) return;
1017   // Temporarily store the Max.Row, column and time bin of the center pad
1018   // Used to later on assign the track indices
1019   cluster.SetLabel(Max.row, 0);
1020   cluster.SetLabel(Max.col, 1);
1021   cluster.SetLabel(Max.time,2);
1022
1023   //needed for HLT reconstruction
1024   AddClusterToArray(&cluster);
1025
1026   // Store the index of the first cluster in the current ROC
1027   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1028
1029   fNoOfClusters++;
1030   fClusterROC++;
1031 }
1032
1033 //_____________________________________________________________________________
1034 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1035 {
1036   // Look to the right
1037   Int_t ii = 1;
1038   while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
1039     nPadCount++;
1040     ii++;
1041     if (Max.col < ii) break;
1042   }
1043   // Look to the left
1044   ii = 1;
1045   while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
1046     nPadCount++;
1047     ii++;
1048     if (Max.col+ii >= fColMax) break;
1049   }
1050
1051   // Store the amplitudes of the pads in the cluster for later analysis
1052   // and check whether one of these pads is masked in the database
1053   signals[2]=Max.signals[0];
1054   signals[3]=Max.signals[1];
1055   signals[4]=Max.signals[2];
1056   Float_t gain;
1057   for(Int_t i = 0; i<2; i++)
1058     {
1059       if(Max.col+i >= 3){
1060         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
1061         signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
1062       }
1063       if(Max.col+3-i < fColMax){
1064         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
1065         signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
1066       }
1067     }
1068   /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1069     if ((jPad >= 0) && (jPad < fColMax))
1070       signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1071       }*/
1072 }
1073
1074 //_____________________________________________________________________________
1075 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1076 {
1077   //
1078   // Add a cluster to the array
1079   //
1080
1081   Int_t n = RecPoints()->GetEntriesFast();
1082   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1083   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1084 }
1085
1086 //_____________________________________________________________________________
1087 void AliTRDclusterizer::AddTrackletsToArray()
1088 {
1089   //
1090   // Add the online tracklets of this chamber to the array
1091   //
1092
1093   UInt_t* trackletword;
1094   for(Int_t side=0; side<2; side++)
1095     {
1096       Int_t trkl=0;
1097       trackletword=fTrackletContainer[side];
1098       while(trackletword[trkl]>0){
1099         Int_t n = TrackletsArray()->GetEntriesFast();
1100         AliTRDtrackletWord tmp(trackletword[trkl]);
1101         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1102         trkl++;
1103       }
1104     }
1105 }
1106
1107 //_____________________________________________________________________________
1108 Bool_t AliTRDclusterizer::AddLabels()
1109 {
1110   //
1111   // Add the track indices to the found clusters
1112   //
1113
1114   const Int_t   kNclus  = 3;
1115   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1116   const Int_t   kNtrack = kNdict * kNclus;
1117
1118   Int_t iClusterROC = 0;
1119
1120   Int_t row  = 0;
1121   Int_t col  = 0;
1122   Int_t time = 0;
1123   Int_t iPad = 0;
1124
1125   // Temporary array to collect the track indices
1126   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1127
1128   // Loop through the dictionary arrays one-by-one
1129   // to keep memory consumption low
1130   AliTRDarrayDictionary *tracksIn = 0;  //mod
1131   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1132
1133     // tracksIn should be expanded beforehand!
1134     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1135
1136     // Loop though the clusters found in this ROC
1137     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1138
1139       AliTRDcluster *cluster = (AliTRDcluster *)
1140                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1141       row  = cluster->GetLabel(0);
1142       col  = cluster->GetLabel(1);
1143       time = cluster->GetLabel(2);
1144
1145       for (iPad = 0; iPad < kNclus; iPad++) {
1146         Int_t iPadCol = col - 1 + iPad;
1147         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1148         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1149       }
1150
1151     }
1152
1153   }
1154
1155   // Copy the track indices into the cluster
1156   // Loop though the clusters found in this ROC
1157   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1158
1159     AliTRDcluster *cluster = (AliTRDcluster *)
1160       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1161     cluster->SetLabel(-9999,0);
1162     cluster->SetLabel(-9999,1);
1163     cluster->SetLabel(-9999,2);
1164
1165     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1166
1167   }
1168
1169   delete [] idxTracks;
1170
1171   return kTRUE;
1172
1173 }
1174
1175 //_____________________________________________________________________________
1176 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1177 {
1178   //
1179   // Method to unfold neighbouring maxima.
1180   // The charge ratio on the overlapping pad is calculated
1181   // until there is no more change within the range given by eps.
1182   // The resulting ratio is then returned to the calling method.
1183   //
1184
1185   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1186   if (!calibration) {
1187     AliError("No AliTRDcalibDB instance available\n");
1188     return kFALSE;
1189   }
1190
1191   Int_t   irc                = 0;
1192   Int_t   itStep             = 0;                 // Count iteration steps
1193
1194   Double_t ratio             = 0.5;               // Start value for ratio
1195   Double_t prevRatio         = 0.0;               // Store previous ratio
1196
1197   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1198   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1199   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1200
1201   // Start the iteration
1202   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1203
1204     itStep++;
1205     prevRatio = ratio;
1206
1207     // Cluster position according to charge ratio
1208     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1209                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1210     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1211                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1212
1213     // Set cluster charge ratio
1214     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1215     Double_t ampLeft  = padSignal[1] / newSignal[1];
1216     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1217     Double_t ampRight = padSignal[3] / newSignal[1];
1218
1219     // Apply pad response to parameters
1220     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1221     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1222
1223     // Calculate new overlapping ratio
1224     ratio = TMath::Min((Double_t) 1.0
1225                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1226
1227   }
1228
1229   return ratio;
1230
1231 }
1232
1233 //_____________________________________________________________________________
1234 void AliTRDclusterizer::TailCancelation()
1235 {
1236   //
1237   // Applies the tail cancelation and gain factors:
1238   // Transform fDigits to fDigits
1239   //
1240
1241   Int_t iRow  = 0;
1242   Int_t iCol  = 0;
1243   Int_t iTime = 0;
1244
1245   Double_t *inADC = new Double_t[fTimeTotal];  // ADC data before tail cancellation
1246   Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
1247
1248   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1249   Bool_t debugStreaming = fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1250   while(fIndexes->NextRCIndex(iRow, iCol))
1251     {
1252       Bool_t corrupted = kFALSE;
1253       if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
1254
1255       // Save data into the temporary processing array and substract the baseline,
1256       // since DeConvExp does not expect a baseline
1257       for (iTime = 0; iTime < fTimeTotal; iTime++)
1258         inADC[iTime]   = fDigits->GetData(iRow,iCol,iTime)-fBaseline;
1259
1260       if(debugStreaming){
1261         for (iTime = 0; iTime < fTimeTotal; iTime++)
1262           (*fDebugStream) << "TailCancellation"
1263                           << "col="  << iCol
1264                           << "row="  << iRow
1265                           << "time=" << iTime
1266                           << "inADC=" << inADC[iTime]
1267                           << "outADC=" << outADC[iTime]
1268                           << "corrupted=" << corrupted
1269                           << "\n";
1270       }
1271
1272       if (!corrupted)
1273         {
1274           // Apply the tail cancelation via the digital filter
1275           // (only for non-coorupted pads)
1276           DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1277         }
1278       else memcpy(&outADC[0],&inADC[0],fTimeTotal*sizeof(inADC[0]));
1279
1280       // Save tailcancalled data and add the baseline
1281       for(iTime = 0; iTime < fTimeTotal; iTime++)
1282         fDigits->SetData(iRow,iCol,iTime,(Short_t)(outADC[iTime] + fBaseline + 0.5));
1283
1284     } // while irow icol
1285
1286   delete [] inADC;
1287   delete [] outADC;
1288
1289   return;
1290
1291 }
1292
1293 //_____________________________________________________________________________
1294 void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1295                                   ,const Int_t n, const Int_t nexp)
1296 {
1297   //
1298   // Tail cancellation by deconvolution for PASA v4 TRF
1299   //
1300
1301   Double_t rates[2];
1302   Double_t coefficients[2];
1303
1304   // Initialization (coefficient = alpha, rates = lambda)
1305   Double_t r1 = 1.0;
1306   Double_t r2 = 1.0;
1307   Double_t c1 = 0.5;
1308   Double_t c2 = 0.5;
1309
1310   if (nexp == 1) {   // 1 Exponentials
1311     r1 = 1.156;
1312     r2 = 0.130;
1313     c1 = 0.066;
1314     c2 = 0.000;
1315   }
1316   if (nexp == 2) {   // 2 Exponentials
1317     Double_t par[4];
1318     fReconstructor->GetRecoParam()->GetTCParams(par);
1319     r1 = par[0];//1.156;
1320     r2 = par[1];//0.130;
1321     c1 = par[2];//0.114;
1322     c2 = par[3];//0.624;
1323   }
1324
1325   coefficients[0] = c1;
1326   coefficients[1] = c2;
1327
1328   Double_t dt = 0.1;
1329
1330   rates[0] = TMath::Exp(-dt/(r1));
1331   rates[1] = TMath::Exp(-dt/(r2));
1332
1333   Int_t i = 0;
1334   Int_t k = 0;
1335
1336   Double_t reminder[2];
1337   Double_t correction = 0.0;
1338   Double_t result     = 0.0;
1339
1340   // Attention: computation order is important
1341   for (k = 0; k < nexp; k++) {
1342     reminder[k] = 0.0;
1343   }
1344
1345   for (i = 0; i < n; i++) {
1346
1347     result    = (source[i] - correction);    // No rescaling
1348     target[i] = result;
1349
1350     for (k = 0; k < nexp; k++) {
1351       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1352     }
1353
1354     correction = 0.0;
1355     for (k = 0; k < nexp; k++) {
1356       correction += reminder[k];
1357     }
1358
1359   }
1360
1361 }
1362
1363 //_____________________________________________________________________________
1364 void AliTRDclusterizer::ResetRecPoints()
1365 {
1366   //
1367   // Resets the list of rec points
1368   //
1369
1370   if (fRecPoints) {
1371     fRecPoints->Delete();
1372     delete fRecPoints;
1373   }
1374 }
1375
1376 //_____________________________________________________________________________
1377 TClonesArray *AliTRDclusterizer::RecPoints()
1378 {
1379   //
1380   // Returns the list of rec points
1381   //
1382
1383   if (!fRecPoints) {
1384     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1385       // determine number of clusters which has to be allocated
1386       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1387
1388       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1389     }
1390     //SetClustersOwner(kTRUE);
1391     AliTRDReconstructor::SetClusters(0x0);
1392   }
1393   return fRecPoints;
1394
1395 }
1396
1397 //_____________________________________________________________________________
1398 TClonesArray *AliTRDclusterizer::TrackletsArray()
1399 {
1400   //
1401   // Returns the list of rec points
1402   //
1403
1404   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1405     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1406     //SetClustersOwner(kTRUE);
1407     //AliTRDReconstructor::SetTracklets(0x0);
1408   }
1409   return fTracklets;
1410
1411 }
1412