TRD/AliTRDclusterizer.cxx

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