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   const Int_t iEveryNTB = recoParam->GetRecEveryNTB();
 787
 788   Int_t istack  = fIndexes->GetStack();
 789   fLayer  = fIndexes->GetLayer();
 790   Int_t isector = fIndexes->GetSM();
 791
 792   // Start clustering in the chamber
 793
 794   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
 795   if (fDet != det) {
 796     AliError("Strange Detector number Missmatch!");
 797     return kFALSE;
 798   }
 799
 800   AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
 801
 802   // TRD space point transformation
 803   fTransform->SetDetector(det);
 804
 805   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
 806   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
 807   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
 808
 809   if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
 810     AddTrackletsToArray();
 811
 812   fColMax    = fDigits->GetNcol();
 813   fTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
 814
 815   // Check consistency between OCDB and raw data
 816   Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS();
 817   if(TestBit(kHLT)){
 818     if((nTimeOCDB > -1) && (fTimeTotal != nTimeOCDB)){
 819       AliWarning(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d]), using raw value"
 820                       ,fTimeTotal,nTimeOCDB));
 821     }
 822   }else{
 823     if(nTimeOCDB == -1){
 824       AliWarning("Undefined number of timebins in OCDB, using value from raw data.");
 825       if(!fTimeTotal>0){
 826         AliError("Number of timebins in raw data is negative, skipping chamber!");
 827         return kFALSE;
 828       }
 829     }else if(nTimeOCDB == -2){
 830       AliError("Mixed number of timebins in OCDB, no reconstruction of TRD data!");
 831       return kFALSE;
 832     }else if(fTimeTotal != nTimeOCDB){
 833       AliError(Form("Number of timebins in raw data does not match OCDB value (RAW[%d] OCDB[%d]), skipping chamber!"
 834                     ,fTimeTotal,nTimeOCDB));
 835       return kFALSE;
 836     }
 837   }
 838
 839   // Detector wise calibration object for the gain factors
 840   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 841   // Calibration object with pad wise values for the gain factors
 842   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
 843   // Calibration value for chamber wise gain factor
 844   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
 845
 846   // Detector wise calibration object for the noise
 847   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
 848   // Calibration object with pad wise values for the noise
 849   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
 850   // Calibration value for chamber wise noise
 851   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
 852
 853   // Calibration object with the pad status
 854   fCalPadStatusROC       = calibration->GetPadStatusROC(fDet);
 855
 856   firstClusterROC = -1;
 857   fClusterROC     =  0;
 858
 859   SetBit(kLUT, recoParam->UseLUT());
 860   SetBit(kGAUS, recoParam->UseGAUS());
 861
 862   // Apply the gain and the tail cancelation via digital filter
 863   if(recoParam->UseTailCancelation()) TailCancelation(recoParam);
 864
 865   MaxStruct curr, last;
 866   Int_t nMaximas = 0, nCorrupted = 0;
 867
 868   // Here the clusterfining is happening
 869
 870   for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
 871     while(fIndexes->NextRCIndex(curr.row, curr.col)){
 872       if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
 873         if(last.row>-1){
 874           if(curr.time==last.time && curr.row==last.row && curr.col==last.col+2) FivePadCluster(last, curr);
 875           CreateCluster(last);
 876         }
 877         last=curr; curr.fivePad=kFALSE;
 878       }
 879     }
 880   }
 881   if(last.row>-1) CreateCluster(last);
 882
 883   if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
 884     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
 885     (*fDebugStream) << "MakeClusters"
 886       << "Detector="   << det
 887       << "NMaxima="    << nMaximas
 888       << "NClusters="  << fClusterROC
 889       << "NCorrupted=" << nCorrupted
 890       << "\n";
 891   }
 892   if (TestBit(kLabels)) AddLabels();
 893
 894   return kTRUE;
 895
 896 }
 897
 898 //_____________________________________________________________________________
 899 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
 900 {
 901   //
 902   // Returns true if this row,col,time combination is a maximum.
 903   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
 904   //
 905
 906   Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
 907   Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
 908   if(Signals[1] <= fMaxThresh) return kFALSE;
 909
 910   Short_t noiseMiddleThresh = (Short_t)(fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row));
 911   if (Signals[1] <= noiseMiddleThresh) return kFALSE;
 912
 913   if (Max.col + 1 >= fColMax || Max.col < 1) return kFALSE;
 914
 915   UChar_t status[3]={
 916     fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
 917    ,fCalPadStatusROC->GetStatus(Max.col,   Max.row)
 918    ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
 919   };
 920
 921   gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
 922   Signals[0] = (Short_t)((fDigits->GetData(Max.row, Max.col-1, Max.time) - fBaseline) / gain + 0.5f);
 923   gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
 924   Signals[2] = (Short_t)((fDigits->GetData(Max.row, Max.col+1, Max.time) - fBaseline) / gain + 0.5f);
 925
 926   if(!(status[0] | status[1] | status[2])) {//all pads are good
 927     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
 928       if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) {
 929         if(Signals[0]<0)Signals[0]=0;
 930         if(Signals[2]<0)Signals[2]=0;
 931         Short_t noiseSumThresh = (Short_t)(fMinLeftRightCutSigma * fCalNoiseDetValue
 932                                            * fCalNoiseROC->GetValue(Max.col, Max.row));
 933         if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
 934         padStatus = 0;
 935         return kTRUE;
 936       }
 937     }
 938   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
 939     if(Signals[0]<0)Signals[0]=0;
 940     if(Signals[2]<0)Signals[2]=0;
 941     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) {
 942       Signals[2]=0;
 943       SetPadStatus(status[2], padStatus);
 944       return kTRUE;
 945     }
 946     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
 947       Signals[0]=0;
 948       SetPadStatus(status[0], padStatus);
 949       return kTRUE;
 950     }
 951     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
 952       Signals[1] = fMaxThresh;
 953       SetPadStatus(status[1], padStatus);
 954       return kTRUE;
 955     }
 956   }
 957   return kFALSE;
 958 }
 959
 960 //_____________________________________________________________________________
 961 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
 962 {
 963   //
 964   // Look for 5 pad cluster with minimum in the middle
 965   // Gives back the ratio
 966   //
 967
 968   if (ThisMax.col >= fColMax - 3) return kFALSE;
 969   Float_t gain;
 970   if (ThisMax.col < fColMax - 5){
 971     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
 972     if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
 973       return kFALSE;
 974   }
 975   if (ThisMax.col > 1) {
 976     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
 977     if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
 978       return kFALSE;
 979   }
 980
 981   const Float_t kEpsilon = 0.01;
 982   Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
 983       NeighbourMax.signals[1], NeighbourMax.signals[2]};
 984
 985   // Unfold the two maxima and set the signal on
 986   // the overlapping pad to the ratio
 987   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
 988   ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
 989   NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
 990   ThisMax.fivePad=kTRUE;
 991   NeighbourMax.fivePad=kTRUE;
 992   return kTRUE;
 993
 994 }
 995
 996 //_____________________________________________________________________________
 997 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
 998 {
 999   //
1000   // Creates a cluster at the given position and saves it in fRecPoints
1001   //
1002
1003   Int_t nPadCount = 1;
1004   Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
1005   if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
1006
1007   AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
1008   cluster.SetNPads(nPadCount);
1009   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1010   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1011   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
1012
1013   cluster.SetFivePad(Max.fivePad);
1014   // set pads status for the cluster
1015   UChar_t maskPosition = GetCorruption(Max.padStatus);
1016   if (maskPosition) {
1017     cluster.SetPadMaskedPosition(maskPosition);
1018     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1019   }
1020
1021   // Transform the local cluster coordinates into calibrated
1022   // space point positions defined in the local tracking system.
1023   // Here the calibration for T0, Vdrift and ExB is applied as well.
1024   if(!fTransform->Transform(&cluster)) return;
1025   // Temporarily store the Max.Row, column and time bin of the center pad
1026   // Used to later on assign the track indices
1027   cluster.SetLabel(Max.row, 0);
1028   cluster.SetLabel(Max.col, 1);
1029   cluster.SetLabel(Max.time,2);
1030
1031   //needed for HLT reconstruction
1032   AddClusterToArray(&cluster);
1033
1034   // Store the index of the first cluster in the current ROC
1035   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1036
1037   fNoOfClusters++;
1038   fClusterROC++;
1039 }
1040
1041 //_____________________________________________________________________________
1042 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1043 {
1044   // Look to the right
1045   Int_t ii = 1;
1046   while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
1047     nPadCount++;
1048     ii++;
1049     if (Max.col < ii) break;
1050   }
1051   // Look to the left
1052   ii = 1;
1053   while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
1054     nPadCount++;
1055     ii++;
1056     if (Max.col+ii >= fColMax) break;
1057   }
1058
1059   // Store the amplitudes of the pads in the cluster for later analysis
1060   // and check whether one of these pads is masked in the database
1061   signals[2]=Max.signals[0];
1062   signals[3]=Max.signals[1];
1063   signals[4]=Max.signals[2];
1064   Float_t gain;
1065   for(Int_t i = 0; i<2; i++)
1066     {
1067       if(Max.col+i >= 3){
1068         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
1069         signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
1070       }
1071       if(Max.col+3-i < fColMax){
1072         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
1073         signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
1074       }
1075     }
1076   /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1077     if ((jPad >= 0) && (jPad < fColMax))
1078       signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1079       }*/
1080 }
1081
1082 //_____________________________________________________________________________
1083 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1084 {
1085   //
1086   // Add a cluster to the array
1087   //
1088
1089   Int_t n = RecPoints()->GetEntriesFast();
1090   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1091   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1092 }
1093
1094 //_____________________________________________________________________________
1095 void AliTRDclusterizer::AddTrackletsToArray()
1096 {
1097   //
1098   // Add the online tracklets of this chamber to the array
1099   //
1100
1101   UInt_t* trackletword;
1102   for(Int_t side=0; side<2; side++)
1103     {
1104       Int_t trkl=0;
1105       trackletword=fTrackletContainer[side];
1106       while(trackletword[trkl]>0){
1107         Int_t n = TrackletsArray()->GetEntriesFast();
1108         AliTRDtrackletWord tmp(trackletword[trkl]);
1109         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1110         trkl++;
1111       }
1112     }
1113 }
1114
1115 //_____________________________________________________________________________
1116 Bool_t AliTRDclusterizer::AddLabels()
1117 {
1118   //
1119   // Add the track indices to the found clusters
1120   //
1121
1122   const Int_t   kNclus  = 3;
1123   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1124   const Int_t   kNtrack = kNdict * kNclus;
1125
1126   Int_t iClusterROC = 0;
1127
1128   Int_t row  = 0;
1129   Int_t col  = 0;
1130   Int_t time = 0;
1131   Int_t iPad = 0;
1132
1133   // Temporary array to collect the track indices
1134   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1135
1136   // Loop through the dictionary arrays one-by-one
1137   // to keep memory consumption low
1138   AliTRDarrayDictionary *tracksIn = 0;  //mod
1139   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1140
1141     // tracksIn should be expanded beforehand!
1142     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1143
1144     // Loop though the clusters found in this ROC
1145     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1146
1147       AliTRDcluster *cluster = (AliTRDcluster *)
1148                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1149       row  = cluster->GetLabel(0);
1150       col  = cluster->GetLabel(1);
1151       time = cluster->GetLabel(2);
1152
1153       for (iPad = 0; iPad < kNclus; iPad++) {
1154         Int_t iPadCol = col - 1 + iPad;
1155         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1156         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1157       }
1158
1159     }
1160
1161   }
1162
1163   // Copy the track indices into the cluster
1164   // Loop though the clusters found in this ROC
1165   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1166
1167     AliTRDcluster *cluster = (AliTRDcluster *)
1168       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1169     cluster->SetLabel(-9999,0);
1170     cluster->SetLabel(-9999,1);
1171     cluster->SetLabel(-9999,2);
1172
1173     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1174
1175   }
1176
1177   delete [] idxTracks;
1178
1179   return kTRUE;
1180
1181 }
1182
1183 //_____________________________________________________________________________
1184 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1185 {
1186   //
1187   // Method to unfold neighbouring maxima.
1188   // The charge ratio on the overlapping pad is calculated
1189   // until there is no more change within the range given by eps.
1190   // The resulting ratio is then returned to the calling method.
1191   //
1192
1193   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1194   if (!calibration) {
1195     AliError("No AliTRDcalibDB instance available\n");
1196     return kFALSE;
1197   }
1198
1199   Int_t   irc                = 0;
1200   Int_t   itStep             = 0;                 // Count iteration steps
1201
1202   Double_t ratio             = 0.5;               // Start value for ratio
1203   Double_t prevRatio         = 0.0;               // Store previous ratio
1204
1205   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1206   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1207   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1208
1209   // Start the iteration
1210   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1211
1212     itStep++;
1213     prevRatio = ratio;
1214
1215     // Cluster position according to charge ratio
1216     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1217                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1218     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1219                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1220
1221     // Set cluster charge ratio
1222     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1223     Double_t ampLeft  = padSignal[1] / newSignal[1];
1224     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1225     Double_t ampRight = padSignal[3] / newSignal[1];
1226
1227     // Apply pad response to parameters
1228     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1229     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1230
1231     // Calculate new overlapping ratio
1232     ratio = TMath::Min((Double_t) 1.0
1233                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1234
1235   }
1236
1237   return ratio;
1238
1239 }
1240
1241 //_____________________________________________________________________________
1242 void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
1243 {
1244   //
1245   // Applies the tail cancelation
1246   //
1247
1248   Int_t iRow  = 0;
1249   Int_t iCol  = 0;
1250   Int_t iTime = 0;
1251
1252   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1253   Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1254   Int_t nexp = recoParam->GetTCnexp();
1255   while(fIndexes->NextRCIndex(iRow, iCol))
1256     {
1257       // if corrupted then don't make the tail cancallation
1258       if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
1259
1260       if(debugStreaming){
1261         for (iTime = 0; iTime < fTimeTotal; iTime++)
1262           (*fDebugStream) << "TailCancellation"
1263                           << "col="  << iCol
1264                           << "row="  << iRow
1265                           << "\n";
1266       }
1267
1268       // Apply the tail cancelation via the digital filter
1269       DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1270
1271     } // while irow icol
1272
1273   return;
1274
1275 }
1276
1277 //_____________________________________________________________________________
1278 void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1279 {
1280   //
1281   // Tail cancellation by deconvolution for PASA v4 TRF
1282   //
1283
1284   Float_t rates[2];
1285   Float_t coefficients[2];
1286
1287   // Initialization (coefficient = alpha, rates = lambda)
1288   Float_t r1 = 1.0;
1289   Float_t r2 = 1.0;
1290   Float_t c1 = 0.5;
1291   Float_t c2 = 0.5;
1292
1293   if (nexp == 1) {   // 1 Exponentials
1294     r1 = 1.156;
1295     r2 = 0.130;
1296     c1 = 0.066;
1297     c2 = 0.000;
1298   }
1299   if (nexp == 2) {   // 2 Exponentials
1300     Double_t par[4];
1301     fReconstructor->GetRecoParam()->GetTCParams(par);
1302     r1 = par[0];//1.156;
1303     r2 = par[1];//0.130;
1304     c1 = par[2];//0.114;
1305     c2 = par[3];//0.624;
1306   }
1307
1308   coefficients[0] = c1;
1309   coefficients[1] = c2;
1310
1311   Double_t dt = 0.1;
1312
1313   rates[0] = TMath::Exp(-dt/(r1));
1314   rates[1] = TMath::Exp(-dt/(r2));
1315
1316   Int_t i = 0;
1317   Int_t k = 0;
1318
1319   Float_t reminder[2];
1320   Float_t correction = 0.0;
1321   Float_t result     = 0.0;
1322
1323   // Attention: computation order is important
1324   for (k = 0; k < nexp; k++) {
1325     reminder[k] = 0.0;
1326   }
1327
1328   for (i = 0; i < nTime; i++) {
1329
1330     result = arr[i] - correction - fBaseline;    // No rescaling
1331     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1332
1333     for (k = 0; k < nexp; k++) {
1334       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1335     }
1336
1337     correction = 0.0;
1338     for (k = 0; k < nexp; k++) {
1339       correction += reminder[k];
1340     }
1341
1342   }
1343
1344 }
1345
1346 //_____________________________________________________________________________
1347 void AliTRDclusterizer::ResetRecPoints()
1348 {
1349   //
1350   // Resets the list of rec points
1351   //
1352
1353   if (fRecPoints) {
1354     fRecPoints->Clear();
1355     fNoOfClusters = 0;
1356     //    delete fRecPoints;
1357   }
1358 }
1359
1360 //_____________________________________________________________________________
1361 TClonesArray *AliTRDclusterizer::RecPoints()
1362 {
1363   //
1364   // Returns the list of rec points
1365   //
1366
1367   if (!fRecPoints) {
1368     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1369       // determine number of clusters which has to be allocated
1370       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1371
1372       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1373     }
1374     //SetClustersOwner(kTRUE);
1375     AliTRDReconstructor::SetClusters(0x0);
1376   }
1377   return fRecPoints;
1378
1379 }
1380
1381 //_____________________________________________________________________________
1382 TClonesArray *AliTRDclusterizer::TrackletsArray()
1383 {
1384   //
1385   // Returns the list of rec points
1386   //
1387
1388   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1389     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1390     //SetClustersOwner(kTRUE);
1391     //AliTRDReconstructor::SetTracklets(0x0);
1392   }
1393   return fTracklets;
1394
1395 }
1396