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