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