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