TRD/AliTRDdigitizer.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 //  Creates and handles digits from TRD hits                              //
  21 //                                                                        //
  22 //  Authors: C. Blume (blume@ikf.uni-frankfurt.de)                        //
  23 //           C. Lippmann                                                  //
  24 //           B. Vulpescu                                                  //
  25 //                                                                        //
  26 //  The following effects are included:                                   //
  27 //      - Diffusion                                                       //
  28 //      - ExB effects                                                     //
  29 //      - Gas gain including fluctuations                                 //
  30 //      - Pad-response (simple Gaussian approximation)                    //
  31 //      - Time-response                                                   //
  32 //      - Electronics noise                                               //
  33 //      - Electronics gain                                                //
  34 //      - Digitization                                                    //
  35 //      - Zero suppression                                                //
  36 //                                                                        //
  37 ////////////////////////////////////////////////////////////////////////////
  38
  39 #include <TGeoManager.h>
  40 #include <TList.h>
  41 #include <TMath.h>
  42 #include <TRandom.h>
  43 #include <TTree.h>
  44
  45 #include "AliRun.h"
  46 #include "AliMC.h"
  47 #include "AliRunLoader.h"
  48 #include "AliLoader.h"
  49 #include "AliConfig.h"
  50 #include "AliDigitizationInput.h"
  51 #include "AliRunLoader.h"
  52 #include "AliLoader.h"
  53 #include "AliLog.h"
  54
  55 #include "AliTRD.h"
  56 #include "AliTRDhit.h"
  57 #include "AliTRDdigitizer.h"
  58 #include "AliTRDarrayDictionary.h"
  59 #include "AliTRDarrayADC.h"
  60 #include "AliTRDarraySignal.h"
  61 #include "AliTRDdigitsManager.h"
  62 #include "AliTRDgeometry.h"
  63 #include "AliTRDpadPlane.h"
  64 #include "AliTRDcalibDB.h"
  65 #include "AliTRDSimParam.h"
  66 #include "AliTRDCommonParam.h"
  67 #include "AliTRDfeeParam.h"
  68 #include "AliTRDmcmSim.h"
  69 #include "AliTRDdigitsParam.h"
  70
  71 #include "Cal/AliTRDCalROC.h"
  72 #include "Cal/AliTRDCalDet.h"
  73 #include "Cal/AliTRDCalOnlineGainTableROC.h"
  74
  75 ClassImp(AliTRDdigitizer)
  76
  77 //_____________________________________________________________________________
  78 AliTRDdigitizer::AliTRDdigitizer()
  79   :AliDigitizer()
  80   ,fRunLoader(0)
  81   ,fDigitsManager(0)
  82   ,fSDigitsManager(0)
  83   ,fSDigitsManagerList(0)
  84   ,fTRD(0)
  85   ,fGeo(0)
  86   ,fMcmSim(new AliTRDmcmSim)
  87   ,fEvent(0)
  88   ,fMasks(0)
  89   ,fCompress(kTRUE)
  90   ,fSDigits(kFALSE)
  91   ,fMergeSignalOnly(kFALSE)
  92 {
  93   //
  94   // AliTRDdigitizer default constructor
  95   //
  96
  97 }
  98
  99 //_____________________________________________________________________________
 100 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
 101   :AliDigitizer(name,title)
 102   ,fRunLoader(0)
 103   ,fDigitsManager(0)
 104   ,fSDigitsManager(0)
 105   ,fSDigitsManagerList(0)
 106   ,fTRD(0)
 107   ,fGeo(0)
 108   ,fMcmSim(new AliTRDmcmSim)
 109   ,fEvent(0)
 110   ,fMasks(0)
 111   ,fCompress(kTRUE)
 112   ,fSDigits(kFALSE)
 113   ,fMergeSignalOnly(kFALSE)
 114 {
 115   //
 116   // AliTRDdigitizer constructor
 117   //
 118
 119 }
 120
 121 //_____________________________________________________________________________
 122 AliTRDdigitizer::AliTRDdigitizer(AliDigitizationInput* digInput
 123                                , const Text_t *name, const Text_t *title)
 124   :AliDigitizer(digInput,name,title)
 125   ,fRunLoader(0)
 126   ,fDigitsManager(0)
 127   ,fSDigitsManager(0)
 128   ,fSDigitsManagerList(0)
 129   ,fTRD(0)
 130   ,fGeo(0)
 131   ,fMcmSim(new AliTRDmcmSim)
 132   ,fEvent(0)
 133   ,fMasks(0)
 134   ,fCompress(kTRUE)
 135   ,fSDigits(kFALSE)
 136   ,fMergeSignalOnly(kFALSE)
 137 {
 138   //
 139   // AliTRDdigitizer constructor
 140   //
 141
 142 }
 143
 144 //_____________________________________________________________________________
 145 AliTRDdigitizer::AliTRDdigitizer(AliDigitizationInput* digInput)
 146   :AliDigitizer(digInput,"AliTRDdigitizer","TRD digitizer")
 147   ,fRunLoader(0)
 148   ,fDigitsManager(0)
 149   ,fSDigitsManager(0)
 150   ,fSDigitsManagerList(0)
 151   ,fTRD(0)
 152   ,fGeo(0)
 153   ,fMcmSim(new AliTRDmcmSim)
 154   ,fEvent(0)
 155   ,fMasks(0)
 156   ,fCompress(kTRUE)
 157   ,fSDigits(kFALSE)
 158   ,fMergeSignalOnly(kFALSE)
 159 {
 160   //
 161   // AliTRDdigitizer constructor
 162   //
 163
 164 }
 165
 166 //_____________________________________________________________________________
 167 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
 168   :AliDigitizer(d)
 169   ,fRunLoader(0)
 170   ,fDigitsManager(0)
 171   ,fSDigitsManager(0)
 172   ,fSDigitsManagerList(0)
 173   ,fTRD(0)
 174   ,fGeo(0)
 175   ,fMcmSim(new AliTRDmcmSim)
 176   ,fEvent(0)
 177   ,fMasks(0)
 178   ,fCompress(d.fCompress)
 179   ,fSDigits(d.fSDigits)
 180   ,fMergeSignalOnly(d.fMergeSignalOnly)
 181 {
 182   //
 183   // AliTRDdigitizer copy constructor
 184   //
 185
 186 }
 187
 188 //_____________________________________________________________________________
 189 AliTRDdigitizer::~AliTRDdigitizer()
 190 {
 191   //
 192   // AliTRDdigitizer destructor
 193   //
 194
 195   if (fDigitsManager) {
 196     delete fDigitsManager;
 197     fDigitsManager      = 0;
 198   }
 199
 200   if (fSDigitsManager) {
 201     // s-digitsmanager will be deleted via list
 202     fSDigitsManager     = 0;
 203   }
 204
 205   if (fSDigitsManagerList) {
 206     fSDigitsManagerList->Delete();
 207     delete fSDigitsManagerList;
 208     fSDigitsManagerList = 0;
 209   }
 210
 211   if (fMasks) {
 212     delete [] fMasks;
 213     fMasks       = 0;
 214   }
 215
 216   if (fMcmSim) {
 217     delete fMcmSim;
 218     fMcmSim = 0;
 219   }
 220
 221   if (fGeo) {
 222     delete fGeo;
 223     fGeo = 0;
 224   }
 225
 226 }
 227
 228 //_____________________________________________________________________________
 229 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
 230 {
 231   //
 232   // Assignment operator
 233   //
 234
 235   if (this != &d) {
 236     ((AliTRDdigitizer &) d).Copy(*this);
 237   }
 238
 239   return *this;
 240
 241 }
 242
 243 //_____________________________________________________________________________
 244 void AliTRDdigitizer::Copy(TObject &d) const
 245 {
 246   //
 247   // Copy function
 248   //
 249
 250   ((AliTRDdigitizer &) d).fRunLoader          = 0;
 251   ((AliTRDdigitizer &) d).fDigitsManager      = 0;
 252   ((AliTRDdigitizer &) d).fSDigitsManager     = 0;
 253   ((AliTRDdigitizer &) d).fSDigitsManagerList = 0;
 254   ((AliTRDdigitizer &) d).fTRD                = 0;
 255   ((AliTRDdigitizer &) d).fGeo                = 0;
 256   ((AliTRDdigitizer &) d).fEvent              = 0;
 257   ((AliTRDdigitizer &) d).fMasks              = 0;
 258   ((AliTRDdigitizer &) d).fCompress           = fCompress;
 259   ((AliTRDdigitizer &) d).fSDigits            = fSDigits;
 260   ((AliTRDdigitizer &) d).fMergeSignalOnly    = fMergeSignalOnly;
 261
 262 }
 263
 264 //_____________________________________________________________________________
 265 void AliTRDdigitizer::Digitize(const Option_t* option)
 266 {
 267   //
 268   // Executes the merging
 269   //
 270
 271   Int_t iInput;
 272
 273   AliTRDdigitsManager *sdigitsManager;
 274
 275   TString optionString = option;
 276   if (optionString.Contains("deb")) {
 277     AliLog::SetClassDebugLevel("AliTRDdigitizer",1);
 278     AliInfo("Called with debug option");
 279   }
 280
 281   // The AliRoot file is already connected by the manager
 282   AliRunLoader *inrl = 0x0;
 283
 284   if (gAlice) {
 285     AliDebug(1,"AliRun object found on file.");
 286   }
 287   else {
 288     inrl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 289     inrl->LoadgAlice();
 290     gAlice = inrl->GetAliRun();
 291     if (!gAlice) {
 292       AliError("Could not find AliRun object.");
 293       return;
 294     }
 295   }
 296
 297   Int_t nInput = fDigInput->GetNinputs();
 298   fMasks       = new Int_t[nInput];
 299   for (iInput = 0; iInput < nInput; iInput++) {
 300     fMasks[iInput] = fDigInput->GetMask(iInput);
 301   }
 302
 303   //
 304   // Initialization
 305   //
 306
 307   AliRunLoader *orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 308
 309   if (InitDetector()) {
 310
 311     AliLoader *ogime = orl->GetLoader("TRDLoader");
 312
 313     TTree *tree = 0;
 314     if (fSDigits) {
 315       // If we produce SDigits
 316       tree = ogime->TreeS();
 317       if (!tree) {
 318         ogime->MakeTree("S");
 319         tree = ogime->TreeS();
 320       }
 321     }
 322     else {
 323       // If we produce Digits
 324       tree = ogime->TreeD();
 325       if (!tree) {
 326         ogime->MakeTree("D");
 327         tree = ogime->TreeD();
 328       }
 329     }
 330
 331     MakeBranch(tree);
 332
 333   }
 334
 335   for (iInput = 0; iInput < nInput; iInput++) {
 336
 337     AliDebug(1,Form("Add input stream %d",iInput));
 338
 339     // Check if the input tree exists
 340     inrl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(iInput));
 341     AliLoader *gime = inrl->GetLoader("TRDLoader");
 342
 343     TTree *treees = gime->TreeS();
 344     if (treees == 0x0) {
 345       if (gime->LoadSDigits()) {
 346         AliError(Form("Error Occured while loading S. Digits for input %d.",iInput));
 347         return;
 348       }
 349       treees = gime->TreeS();
 350     }
 351
 352     if (treees == 0x0) {
 353       AliError(Form("Input stream %d does not exist",iInput));
 354       return;
 355     }
 356
 357     // Read the s-digits via digits manager
 358     sdigitsManager = new AliTRDdigitsManager();
 359     sdigitsManager->SetSDigits(kTRUE);
 360
 361     AliRunLoader *rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(iInput));
 362     AliLoader *gimme = rl->GetLoader("TRDLoader");
 363     if (!gimme->TreeS())
 364       {
 365         gimme->LoadSDigits();
 366       }
 367
 368     sdigitsManager->ReadDigits(gimme->TreeS());
 369
 370     // Add the s-digits to the input list
 371     AddSDigitsManager(sdigitsManager);
 372
 373   }
 374
 375   // Convert the s-digits to normal digits
 376   AliDebug(1,"Do the conversion");
 377   SDigits2Digits();
 378
 379   // Store the digits
 380   AliDebug(1,"Write the digits");
 381   fDigitsManager->WriteDigits();
 382
 383   // Write parameters
 384   orl->CdGAFile();
 385
 386   // Clean up
 387   DeleteSDigitsManager();
 388
 389   AliDebug(1,"Done");
 390
 391 }
 392
 393 //_____________________________________________________________________________
 394 Bool_t AliTRDdigitizer::Open(const Char_t *file, Int_t nEvent)
 395 {
 396   //
 397   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
 398   //
 399   // Connect the AliRoot file containing Geometry, Kine, and Hits
 400   //
 401
 402   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 403
 404   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
 405   if (!fRunLoader) {
 406     fRunLoader = AliRunLoader::Open(file,evfoldname,"UPDATE");
 407   }
 408   if (!fRunLoader) {
 409     AliError(Form("Can not open session for file %s.",file));
 410     return kFALSE;
 411   }
 412
 413   if (!fRunLoader->GetAliRun()) {
 414     fRunLoader->LoadgAlice();
 415   }
 416   gAlice = fRunLoader->GetAliRun();
 417
 418   if (gAlice) {
 419     AliDebug(1,"AliRun object found on file.");
 420   }
 421   else {
 422     AliError("Could not find AliRun object.");
 423     return kFALSE;
 424   }
 425
 426   fEvent = nEvent;
 427
 428   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
 429   if (!loader) {
 430     AliError("Can not get TRD loader from Run Loader");
 431     return kFALSE;
 432   }
 433
 434   if (InitDetector()) {
 435     TTree *tree = 0;
 436     if (fSDigits) {
 437       // If we produce SDigits
 438       tree = loader->TreeS();
 439       if (!tree) {
 440         loader->MakeTree("S");
 441         tree = loader->TreeS();
 442       }
 443     }
 444     else {
 445       // If we produce Digits
 446       tree = loader->TreeD();
 447       if (!tree) {
 448         loader->MakeTree("D");
 449         tree = loader->TreeD();
 450       }
 451     }
 452     return MakeBranch(tree);
 453   }
 454   else {
 455     return kFALSE;
 456   }
 457
 458 }
 459
 460 //_____________________________________________________________________________
 461 Bool_t AliTRDdigitizer::Open(AliRunLoader * const runLoader, Int_t nEvent)
 462 {
 463   //
 464   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
 465   //
 466   // Connect the AliRoot file containing Geometry, Kine, and Hits
 467   //
 468
 469   fRunLoader = runLoader;
 470   if (!fRunLoader) {
 471     AliError("RunLoader does not exist");
 472     return kFALSE;
 473   }
 474
 475   if (!fRunLoader->GetAliRun()) {
 476     fRunLoader->LoadgAlice();
 477   }
 478   gAlice = fRunLoader->GetAliRun();
 479
 480   if (gAlice) {
 481     AliDebug(1,"AliRun object found on file.");
 482   }
 483   else {
 484     AliError("Could not find AliRun object.");
 485     return kFALSE;
 486   }
 487
 488   fEvent = nEvent;
 489
 490   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
 491   if (!loader) {
 492     AliError("Can not get TRD loader from Run Loader");
 493     return kFALSE;
 494   }
 495
 496   if (InitDetector()) {
 497     TTree *tree = 0;
 498     if (fSDigits) {
 499       // If we produce SDigits
 500       tree = loader->TreeS();
 501       if (!tree) {
 502         loader->MakeTree("S");
 503         tree = loader->TreeS();
 504       }
 505     }
 506     else {
 507       // If we produce Digits
 508       tree = loader->TreeD();
 509       if (!tree) {
 510         loader->MakeTree("D");
 511         tree = loader->TreeD();
 512       }
 513     }
 514     return MakeBranch(tree);
 515   }
 516   else {
 517     return kFALSE;
 518   }
 519
 520 }
 521
 522 //_____________________________________________________________________________
 523 Bool_t AliTRDdigitizer::InitDetector()
 524 {
 525   //
 526   // Sets the pointer to the TRD detector and the geometry
 527   //
 528
 529   // Get the pointer to the detector class and check for version 1
 530   fTRD = (AliTRD *) gAlice->GetDetector("TRD");
 531   if (!fTRD) {
 532     AliFatal("No TRD module found");
 533     exit(1);
 534   }
 535   if (fTRD->IsVersion() != 1) {
 536     AliFatal("TRD must be version 1 (slow simulator)");
 537     exit(1);
 538   }
 539
 540   // Get the geometry
 541   fGeo = new AliTRDgeometry();
 542
 543   // Create a digits manager
 544   if (fDigitsManager) {
 545     delete fDigitsManager;
 546   }
 547   fDigitsManager = new AliTRDdigitsManager();
 548   fDigitsManager->SetSDigits(fSDigits);
 549   fDigitsManager->CreateArrays();
 550   fDigitsManager->SetEvent(fEvent);
 551
 552   // The list for the input s-digits manager to be merged
 553   if (fSDigitsManagerList) {
 554     fSDigitsManagerList->Delete();
 555   }
 556   else {
 557     fSDigitsManagerList = new TList();
 558   }
 559
 560   return kTRUE;
 561
 562 }
 563
 564 //_____________________________________________________________________________
 565 Bool_t AliTRDdigitizer::MakeBranch(TTree *tree) const
 566 {
 567   //
 568   // Create the branches for the digits array
 569   //
 570
 571   return fDigitsManager->MakeBranch(tree);
 572
 573 }
 574
 575 //_____________________________________________________________________________
 576 void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
 577 {
 578   //
 579   // Add a digits manager for s-digits to the input list.
 580   //
 581
 582   fSDigitsManagerList->Add(man);
 583
 584 }
 585
 586 //_____________________________________________________________________________
 587 void AliTRDdigitizer::DeleteSDigitsManager()
 588 {
 589   //
 590   // Removes digits manager from the input list.
 591   //
 592
 593   fSDigitsManagerList->Delete();
 594
 595 }
 596
 597 //_____________________________________________________________________________
 598 Bool_t AliTRDdigitizer::MakeDigits()
 599 {
 600   //
 601   // Creates digits.
 602   //
 603
 604   AliDebug(1,"Start creating digits");
 605
 606   if (!fGeo) {
 607     AliError("No geometry defined");
 608     return kFALSE;
 609   }
 610
 611   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 612   if (!calibration) {
 613     AliFatal("Could not get calibration object");
 614     return kFALSE;
 615   }
 616
 617   const Int_t kNdet  = AliTRDgeometry::Ndet();
 618
 619   Float_t **hits = new Float_t*[kNdet];
 620   Int_t    *nhit = new Int_t[kNdet];
 621   memset(nhit,0,kNdet*sizeof(Int_t));
 622
 623   AliTRDarraySignal *signals = 0x0;
 624
 625   // Check the number of time bins from simParam against OCDB,
 626   // if OCDB value is not supposed to be used.
 627   // As default, the value from OCDB is taken
 628   if (AliTRDSimParam::Instance()->GetNTBoverwriteOCDB()) {
 629     if (calibration->GetNumberOfTimeBinsDCS() != AliTRDSimParam::Instance()->GetNTimeBins()) {
 630       AliWarning(Form("Number of time bins is different to OCDB value [SIM=%d, OCDB=%d]"
 631                      ,AliTRDSimParam::Instance()->GetNTimeBins()
 632                      ,calibration->GetNumberOfTimeBinsDCS()));
 633     }
 634     // Save the values for the raw data headers
 635     fDigitsManager->GetDigitsParam()->SetNTimeBinsAll(AliTRDSimParam::Instance()->GetNTimeBins());
 636   }
 637   else {
 638     // Save the values for the raw data headers
 639     fDigitsManager->GetDigitsParam()->SetNTimeBinsAll(calibration->GetNumberOfTimeBinsDCS());
 640   }
 641
 642   // Save the values for the raw data headers
 643   fDigitsManager->GetDigitsParam()->SetADCbaselineAll(AliTRDSimParam::Instance()->GetADCbaseline());
 644
 645   // Sort all hits according to detector number
 646   if (!SortHits(hits,nhit)) {
 647     AliError("Sorting hits failed");
 648     delete [] hits;
 649     delete [] nhit;
 650     return kFALSE;
 651   }
 652
 653   // Loop through all detectors
 654   for (Int_t det = 0; det < kNdet; det++) {
 655
 656     // Detectors that are switched off, not installed, etc.
 657     if ((!calibration->IsChamberNoData(det))    &&
 658         ( fGeo->ChamberInGeometry(det))         &&
 659         (nhit[det] > 0)) {
 660
 661       signals = new AliTRDarraySignal();
 662
 663       // Convert the hits of the current detector to detector signals
 664       if (!ConvertHits(det,hits[det],nhit[det],signals)) {
 665         AliError(Form("Conversion of hits failed for detector=%d",det));
 666         delete [] hits;
 667         delete [] nhit;
 668         delete signals;
 669         signals = 0x0;
 670         return kFALSE;
 671       }
 672
 673       // Convert the detector signals to digits or s-digits
 674       if (!ConvertSignals(det,signals)) {
 675         AliError(Form("Conversion of signals failed for detector=%d",det));
 676         delete [] hits;
 677         delete [] nhit;
 678         delete signals;
 679         signals = 0x0;
 680         return kFALSE;
 681       }
 682
 683       // Delete the signals array
 684       delete signals;
 685       signals = 0x0;
 686
 687     } // if: detector status
 688
 689     delete [] hits[det];
 690
 691   } // for: detector
 692
 693   if (!fSDigits) {
 694     if (AliDataLoader *trklLoader
 695           = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
 696       if (trklLoader->Tree())
 697         trklLoader->WriteData("OVERWRITE");
 698     }
 699   }
 700
 701   delete [] hits;
 702   delete [] nhit;
 703
 704   return kTRUE;
 705
 706 }
 707
 708 //_____________________________________________________________________________
 709 Bool_t AliTRDdigitizer::SortHits(Float_t **hits, Int_t *nhit)
 710 {
 711   //
 712   // Read all the hits and sorts them according to detector number
 713   // in the output array <hits>.
 714   //
 715
 716   AliDebug(1,"Start sorting hits");
 717
 718   const Int_t kNdet = AliTRDgeometry::Ndet();
 719   // Size of the hit vector
 720   const Int_t kNhit = 6;
 721
 722   Float_t *xyz      = 0;
 723   Int_t    nhitTrk  = 0;
 724
 725   Int_t   *lhit     = new Int_t[kNdet];
 726   memset(lhit,0,kNdet*sizeof(Int_t));
 727
 728   for (Int_t det = 0; det < kNdet; det++) {
 729     hits[det] = 0x0;
 730   }
 731
 732   AliLoader *gimme = fRunLoader->GetLoader("TRDLoader");
 733   if (!gimme->TreeH()) {
 734     gimme->LoadHits();
 735   }
 736   TTree *hitTree = gimme->TreeH();
 737   if (hitTree == 0x0) {
 738     AliError("Can not get TreeH");
 739     delete [] lhit;
 740     return kFALSE;
 741   }
 742   fTRD->SetTreeAddress();
 743
 744   // Get the number of entries in the hit tree
 745   // (Number of primary particles creating a hit somewhere)
 746   Int_t nTrk = (Int_t) hitTree->GetEntries();
 747   AliDebug(1,Form("Found %d tracks",nTrk));
 748
 749   // Loop through all the tracks in the tree
 750   for (Int_t iTrk = 0; iTrk < nTrk; iTrk++) {
 751
 752     gAlice->GetMCApp()->ResetHits();
 753     hitTree->GetEvent(iTrk);
 754
 755     if (!fTRD->Hits()) {
 756       AliError(Form("No hits array for track = %d",iTrk));
 757       continue;
 758     }
 759
 760     // Number of hits for this track
 761     nhitTrk = fTRD->Hits()->GetEntriesFast();
 762
 763     Int_t hitCnt = 0;
 764     // Loop through the TRD hits
 765     AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
 766     while (hit) {
 767
 768       hitCnt++;
 769
 770       // Don't analyze test hits
 771       if (((Int_t) hit->GetCharge()) != 0) {
 772
 773         Int_t   trk  = hit->Track();
 774         Int_t   det  = hit->GetDetector();
 775         Int_t   q    = hit->GetCharge();
 776         Float_t x    = hit->X();
 777         Float_t y    = hit->Y();
 778         Float_t z    = hit->Z();
 779         Float_t time = hit->GetTime();
 780
 781         if (nhit[det] == lhit[det]) {
 782           // Inititialization of new detector
 783           xyz        = new Float_t[kNhit*(nhitTrk+lhit[det])];
 784           if (hits[det]) {
 785             memcpy(xyz,hits[det],sizeof(Float_t)*kNhit*lhit[det]);
 786             delete [] hits[det];
 787           }
 788           lhit[det] += nhitTrk;
 789           hits[det]  = xyz;
 790         }
 791         else {
 792           xyz        = hits[det];
 793         }
 794         xyz[nhit[det]*kNhit+0] = x;
 795         xyz[nhit[det]*kNhit+1] = y;
 796         xyz[nhit[det]*kNhit+2] = z;
 797         xyz[nhit[det]*kNhit+3] = q;
 798         xyz[nhit[det]*kNhit+4] = trk;
 799         xyz[nhit[det]*kNhit+5] = time;
 800         nhit[det]++;
 801
 802       } // if: charge != 0
 803
 804       hit = (AliTRDhit *) fTRD->NextHit();
 805
 806     } // for: hits of one track
 807
 808   } // for: tracks
 809
 810   delete [] lhit;
 811
 812   return kTRUE;
 813
 814 }
 815
 816 //_____________________________________________________________________________
 817 Bool_t AliTRDdigitizer::ConvertHits(Int_t det
 818                                   , const Float_t * const hits
 819                                   , Int_t nhit
 820                                   , AliTRDarraySignal *signals)
 821 {
 822   //
 823   // Converts the detectorwise sorted hits to detector signals
 824   //
 825
 826   AliDebug(1,Form("Start converting hits for detector=%d (nhits=%d)",det,nhit));
 827
 828   // Number of pads included in the pad response
 829   const Int_t kNpad      = 3;
 830   // Number of track dictionary arrays
 831   const Int_t kNdict     = AliTRDdigitsManager::kNDict;
 832   // Size of the hit vector
 833   const Int_t kNhit      = 6;
 834
 835   // Width of the amplification region
 836   const Float_t kAmWidth = AliTRDgeometry::AmThick();
 837   // Width of the drift region
 838   const Float_t kDrWidth = AliTRDgeometry::DrThick();
 839   // Drift + amplification region
 840   const Float_t kDrMin   =          - 0.5 * kAmWidth;
 841   const Float_t kDrMax   = kDrWidth + 0.5 * kAmWidth;
 842
 843   Int_t    iPad          = 0;
 844   Int_t    dict          = 0;
 845   Int_t    timeBinTRFend = 1;
 846
 847   Double_t pos[3];
 848   Double_t loc[3];
 849   Double_t padSignal[kNpad];
 850   Double_t signalOld[kNpad];
 851
 852   AliTRDarrayDictionary *dictionary[kNdict];
 853
 854   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 855   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
 856   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
 857
 858   if (!commonParam) {
 859     AliFatal("Could not get common parameterss");
 860     return kFALSE;
 861   }
 862   if (!simParam) {
 863     AliFatal("Could not get simulation parameters");
 864     return kFALSE;
 865   }
 866   if (!calibration) {
 867     AliFatal("Could not get calibration object");
 868     return kFALSE;
 869   }
 870
 871   // Get the detector wise calibration objects
 872   AliTRDCalROC       *calVdriftROC      = 0;
 873   Float_t             calVdriftDetValue = 0.0;
 874   const AliTRDCalDet *calVdriftDet      = calibration->GetVdriftDet();
 875   AliTRDCalROC       *calT0ROC          = 0;
 876   Float_t             calT0DetValue     = 0.0;
 877   const AliTRDCalDet *calT0Det          = calibration->GetT0Det();
 878   Double_t            calExBDetValue    = 0.0;
 879   const AliTRDCalDet *calExBDet         = calibration->GetExBDet();
 880
 881   if (simParam->TRFOn()) {
 882     timeBinTRFend = ((Int_t) (simParam->GetTRFhi()
 883                   * commonParam->GetSamplingFrequency())) - 1;
 884   }
 885
 886   Int_t   nTimeTotal   = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
 887   Float_t samplingRate = commonParam->GetSamplingFrequency();
 888   Float_t elAttachProp = simParam->GetElAttachProp() / 100.0;
 889
 890   AliTRDpadPlane *padPlane = fGeo->GetPadPlane(det);
 891   Int_t   layer   = fGeo->GetLayer(det);         //update
 892   Float_t row0    = padPlane->GetRow0ROC();
 893   Int_t   nRowMax = padPlane->GetNrows();
 894   Int_t   nColMax = padPlane->GetNcols();
 895
 896   // Create a new array for the signals
 897   signals->Allocate(nRowMax,nColMax,nTimeTotal);
 898
 899   // Create a new array for the dictionary
 900   for (dict = 0; dict < kNdict; dict++) {
 901     dictionary[dict] = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(det,dict);
 902     dictionary[dict]->Allocate(nRowMax,nColMax,nTimeTotal);
 903   }
 904
 905   // Loop through the hits in this detector
 906   for (Int_t hit = 0; hit < nhit; hit++) {
 907
 908     pos[0]          = hits[hit*kNhit+0];
 909     pos[1]          = hits[hit*kNhit+1];
 910     pos[2]          = hits[hit*kNhit+2];
 911     Float_t q       = hits[hit*kNhit+3];
 912     Float_t hittime = hits[hit*kNhit+5];
 913     Int_t   track   = ((Int_t) hits[hit*kNhit+4]);
 914
 915     Int_t   inDrift = 1;
 916
 917     // Find the current volume with the geo manager
 918     gGeoManager->SetCurrentPoint(pos);
 919     gGeoManager->FindNode();
 920     if (strstr(gGeoManager->GetPath(),"/UK")) {
 921       inDrift = 0;
 922     }
 923
 924     // Get the calibration objects
 925     calVdriftROC      = calibration->GetVdriftROC(det);
 926     calVdriftDetValue = calVdriftDet->GetValue(det);
 927     calT0ROC          = calibration->GetT0ROC(det);
 928     calT0DetValue     = calT0Det->GetValue(det);
 929     calExBDetValue    = calExBDet->GetValue(det);
 930
 931     // Go to the local coordinate system:
 932     // loc[0] - col  direction in amplification or driftvolume
 933     // loc[1] - row  direction in amplification or driftvolume
 934     // loc[2] - time direction in amplification or driftvolume
 935     gGeoManager->MasterToLocal(pos,loc);
 936     if (inDrift) {
 937       // Relative to middle of amplification region
 938       loc[2] = loc[2] - kDrWidth/2.0 - kAmWidth/2.0;
 939     }
 940
 941     // The driftlength [cm] (w/o diffusion yet !).
 942     // It is negative if the hit is between pad plane and anode wires.
 943     Double_t driftlength = -1.0 * loc[2];
 944
 945     // Stupid patch to take care of TR photons that are absorbed
 946     // outside the chamber volume. A real fix would actually need
 947     // a more clever implementation of the TR hit generation
 948     if (q < 0.0) {
 949       if ((loc[1] < padPlane->GetRowEndROC()) ||
 950           (loc[1] > padPlane->GetRow0ROC())) {
 951         continue;
 952       }
 953       if ((driftlength < kDrMin) ||
 954           (driftlength > kDrMax)) {
 955         continue;
 956       }
 957     }
 958
 959     // Get row and col of unsmeared electron to retrieve drift velocity
 960     // The pad row (z-direction)
 961     Int_t    rowE         = padPlane->GetPadRowNumberROC(loc[1]);
 962     if (rowE < 0) {
 963       continue;
 964     }
 965     Double_t rowOffset    = padPlane->GetPadRowOffsetROC(rowE,loc[1]);
 966     // The pad column (rphi-direction)
 967     Double_t offsetTilt   = padPlane->GetTiltOffset(rowOffset);
 968     Int_t    colE         = padPlane->GetPadColNumber(loc[0]+offsetTilt);
 969     if (colE < 0) {
 970       continue;
 971     }
 972     Double_t colOffset    = 0.0;
 973
 974     // Normalized drift length
 975     Float_t  driftvelocity  = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
 976     Double_t absdriftlength = TMath::Abs(driftlength);
 977     if (commonParam->ExBOn()) {
 978       absdriftlength /= TMath::Sqrt(1.0 / (1.0 + calExBDetValue*calExBDetValue));
 979     }
 980
 981     // Loop over all electrons of this hit
 982     // TR photons produce hits with negative charge
 983     Int_t nEl = ((Int_t) TMath::Abs(q));
 984     for (Int_t iEl = 0; iEl < nEl; iEl++) {
 985
 986       // Now the real local coordinate system of the ROC
 987       // column direction: locC
 988       // row direction:    locR
 989       // time direction:   locT
 990       // locR and locC are identical to the coordinates of the corresponding
 991       // volumina of the drift or amplification region.
 992       // locT is defined relative to the wire plane (i.e. middle of amplification
 993       // region), meaning locT = 0, and is negative for hits coming from the
 994       // drift region.
 995       Double_t locC = loc[0];
 996       Double_t locR = loc[1];
 997       Double_t locT = loc[2];
 998
 999       // Electron attachment
1000       if (simParam->ElAttachOn()) {
1001         if (gRandom->Rndm() < (absdriftlength * elAttachProp)) {
1002           continue;
1003         }
1004       }
1005
1006       // Apply the diffusion smearing
1007       if (simParam->DiffusionOn()) {
1008         if (!(Diffusion(driftvelocity,absdriftlength,calExBDetValue,locR,locC,locT))) {
1009           continue;
1010         }
1011       }
1012
1013       // Apply E x B effects (depends on drift direction)
1014       if (commonParam->ExBOn()) {
1015         locC = locC + calExBDetValue * driftlength;
1016       }
1017
1018       // The electron position after diffusion and ExB in pad coordinates.
1019       // The pad row (z-direction)
1020       rowE       = padPlane->GetPadRowNumberROC(locR);
1021       if (rowE < 0) continue;
1022       rowOffset  = padPlane->GetPadRowOffsetROC(rowE,locR);
1023
1024       // The pad column (rphi-direction)
1025       offsetTilt = padPlane->GetTiltOffset(rowOffset);
1026       colE       = padPlane->GetPadColNumber(locC+offsetTilt);
1027       if (colE < 0) continue;
1028       colOffset  = padPlane->GetPadColOffset(colE,locC+offsetTilt);
1029
1030       // Also re-retrieve drift velocity because col and row may have changed
1031       driftvelocity = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
1032       Float_t t0    = calT0DetValue     + calT0ROC->GetValue(colE,rowE);
1033
1034       // Convert the position to drift time [mus], using either constant drift velocity or
1035       // time structure of drift cells (non-isochronity, GARFIELD calculation).
1036       // Also add absolute time of hits to take pile-up events into account properly
1037       Double_t drifttime;
1038       if (simParam->TimeStructOn()) {
1039         // Get z-position with respect to anode wire
1040         Double_t zz  =  row0 - locR + padPlane->GetAnodeWireOffset();
1041         zz -= ((Int_t)(2 * zz)) / 2.0;
1042         if (zz > 0.25) {
1043           zz  = 0.5 - zz;
1044         }
1045         // Use drift time map (GARFIELD)
1046         drifttime = commonParam->TimeStruct(driftvelocity,0.5*kAmWidth-1.0*locT,zz)
1047                   + hittime;
1048       }
1049       else {
1050         // Use constant drift velocity
1051         drifttime = TMath::Abs(locT) / driftvelocity
1052                   + hittime;
1053       }
1054
1055       // Apply the gas gain including fluctuations
1056       Double_t ggRndm = 0.0;
1057       do {
1058         ggRndm = gRandom->Rndm();
1059       } while (ggRndm <= 0);
1060       Double_t signal = -(simParam->GetGasGain()) * TMath::Log(ggRndm);
1061
1062       // Apply the pad response
1063       if (simParam->PRFOn()) {
1064         // The distance of the electron to the center of the pad
1065         // in units of pad width
1066         Double_t dist = (colOffset - 0.5*padPlane->GetColSize(colE))
1067                       / padPlane->GetColSize(colE);
1068         // This is a fixed parametrization, i.e. not dependent on
1069         // calibration values !
1070         if (!(calibration->PadResponse(signal,dist,layer,padSignal))) continue;
1071       }
1072       else {
1073         padSignal[0] = 0.0;
1074         padSignal[1] = signal;
1075         padSignal[2] = 0.0;
1076       }
1077
1078       // The time bin (always positive), with t0 distortion
1079       Double_t timeBinIdeal = drifttime * samplingRate + t0;
1080       // Protection
1081       if (TMath::Abs(timeBinIdeal) > 2*nTimeTotal) {
1082         timeBinIdeal = 2 * nTimeTotal;
1083       }
1084       Int_t    timeBinTruncated = ((Int_t) timeBinIdeal);
1085       // The distance of the position to the middle of the timebin
1086       Double_t timeOffset       = ((Float_t) timeBinTruncated
1087                                 + 0.5 - timeBinIdeal) / samplingRate;
1088
1089       // Sample the time response inside the drift region
1090       // + additional time bins before and after.
1091       // The sampling is done always in the middle of the time bin
1092       for (Int_t iTimeBin = TMath::Max(timeBinTruncated,0)
1093           ;iTimeBin < TMath::Min(timeBinTruncated+timeBinTRFend,nTimeTotal)
1094           ;iTimeBin++) {
1095
1096         // Apply the time response
1097         Double_t timeResponse = 1.0;
1098         Double_t crossTalk    = 0.0;
1099         Double_t time         = (iTimeBin - timeBinTruncated) / samplingRate + timeOffset;
1100
1101         if (simParam->TRFOn()) {
1102           timeResponse = simParam->TimeResponse(time);
1103         }
1104         if (simParam->CTOn()) {
1105           crossTalk    = simParam->CrossTalk(time);
1106         }
1107
1108         signalOld[0] = 0.0;
1109         signalOld[1] = 0.0;
1110         signalOld[2] = 0.0;
1111
1112         for (iPad = 0; iPad < kNpad; iPad++) {
1113
1114           Int_t colPos = colE + iPad - 1;
1115           if (colPos <        0) continue;
1116           if (colPos >= nColMax) break;
1117
1118           // Add the signals
1119           signalOld[iPad] = signals->GetData(rowE,colPos,iTimeBin);
1120
1121           if (colPos != colE) {
1122             // Cross talk added to non-central pads
1123             signalOld[iPad] += padSignal[iPad]
1124                              * (timeResponse + crossTalk);
1125           }
1126           else {
1127             // W/o cross talk at central pad
1128             signalOld[iPad] += padSignal[iPad]
1129                              * timeResponse;
1130           }
1131
1132           signals->SetData(rowE,colPos,iTimeBin,signalOld[iPad]);
1133
1134           // Store the track index in the dictionary
1135           // Note: We store index+1 in order to allow the array to be compressed
1136           // Note2: Taking out the +1 in track
1137           if (signalOld[iPad] > 0.0) {
1138             for (dict = 0; dict < kNdict; dict++) {
1139               Int_t oldTrack = dictionary[dict]->GetData(rowE,colPos,iTimeBin);
1140               if (oldTrack == track) break;
1141               if (oldTrack ==  -1 ) {
1142                 dictionary[dict]->SetData(rowE,colPos,iTimeBin,track);
1143                 break;
1144               }
1145             }
1146           }
1147
1148         } // Loop: pads
1149
1150       } // Loop: time bins
1151
1152     } // Loop: electrons of a single hit
1153
1154   } // Loop: hits
1155
1156   AliDebug(2,Form("Finished analyzing %d hits",nhit));
1157
1158   return kTRUE;
1159
1160 }
1161
1162 //_____________________________________________________________________________
1163 Bool_t AliTRDdigitizer::ConvertSignals(Int_t det, AliTRDarraySignal *signals)
1164 {
1165   //
1166   // Convert signals to digits
1167   //
1168
1169   AliDebug(1,Form("Start converting the signals for detector %d",det));
1170
1171   if (fSDigits) {
1172     // Convert the signal array to s-digits
1173     if (!Signal2SDigits(det,signals)) {
1174       return kFALSE;
1175     }
1176   }
1177   else {
1178     // Convert the signal array to digits
1179     if (!Signal2ADC(det,signals)) {
1180       return kFALSE;
1181     }
1182     // Run digital processing for digits
1183     RunDigitalProcessing(det);
1184   }
1185
1186   // Compress the arrays
1187   CompressOutputArrays(det);
1188
1189   return kTRUE;
1190
1191 }
1192
1193 //_____________________________________________________________________________
1194 Bool_t AliTRDdigitizer::Signal2ADC(Int_t det, AliTRDarraySignal *signals)
1195 {
1196   //
1197   // Converts the sampled electron signals to ADC values for a given chamber
1198   //
1199
1200   AliDebug(1,Form("Start converting signals to ADC values for detector=%d",det));
1201
1202   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1203   if (!calibration) {
1204     AliFatal("Could not get calibration object");
1205     return kFALSE;
1206   }
1207
1208   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
1209   if (!simParam) {
1210     AliFatal("Could not get simulation parameters");
1211     return kFALSE;
1212   }
1213
1214   // Converts number of electrons to fC
1215   const Double_t kEl2fC = 1.602e-19 * 1.0e15;
1216
1217   // Coupling factor
1218   Double_t coupling     = simParam->GetPadCoupling()
1219                         * simParam->GetTimeCoupling();
1220   // Electronics conversion factor
1221   Double_t convert      = kEl2fC
1222                         * simParam->GetChipGain();
1223   // ADC conversion factor
1224   Double_t adcConvert   = simParam->GetADCoutRange()
1225                         / simParam->GetADCinRange();
1226   // The electronics baseline in mV
1227   Double_t baseline     = simParam->GetADCbaseline()
1228                         / adcConvert;
1229   // The electronics baseline in electrons
1230   Double_t baselineEl   = baseline
1231                         / convert;
1232
1233   Int_t row  = 0;
1234   Int_t col  = 0;
1235   Int_t time = 0;
1236
1237   Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
1238   Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
1239   Int_t nTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
1240   if (fSDigitsManager->GetDigitsParam()->GetNTimeBins(det)) {
1241     nTimeTotal = fSDigitsManager->GetDigitsParam()->GetNTimeBins(det);
1242   }
1243   else {
1244     AliFatal("Could not get number of time bins");
1245     return kFALSE;
1246   }
1247
1248   // The gain factor calibration objects
1249   const AliTRDCalDet *calGainFactorDet      = calibration->GetGainFactorDet();
1250   AliTRDCalROC       *calGainFactorROC      = 0x0;
1251   Float_t             calGainFactorDetValue = 0.0;
1252
1253   AliTRDarrayADC     *digits = 0x0;
1254
1255   if (!signals) {
1256     AliError(Form("Signals array for detector %d does not exist\n",det));
1257     return kFALSE;
1258   }
1259   if (signals->HasData()) {
1260     // Expand the container if neccessary
1261     signals->Expand();
1262   }
1263   else {
1264     // Create missing containers
1265     signals->Allocate(nRowMax,nColMax,nTimeTotal);
1266   }
1267
1268   // Get the container for the digits of this detector
1269   if (fDigitsManager->HasSDigits()) {
1270     AliError("Digits manager has s-digits");
1271     return kFALSE;
1272   }
1273
1274   digits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det);
1275   // Allocate memory space for the digits buffer
1276   if (!digits->HasData()) {
1277     digits->Allocate(nRowMax,nColMax,nTimeTotal);
1278   }
1279
1280   // Get the calibration objects
1281   calGainFactorROC      = calibration->GetGainFactorROC(det);
1282   calGainFactorDetValue = calGainFactorDet->GetValue(det);
1283
1284   // Get the online gain factors
1285   //AliTRDCalOnlineGainTableROC *onlineGainFactorROC
1286   //  = calibration->GetOnlineGainTableROC(det);
1287
1288   // Create the digits for this chamber
1289   for (row  = 0; row  <  nRowMax; row++ ) {
1290     for (col  = 0; col  <  nColMax; col++ ) {
1291
1292       // Check whether pad is masked
1293       // Bridged pads are not considered yet!!!
1294       if (calibration->IsPadMasked(det,col,row) ||
1295           calibration->IsPadNotConnected(det,col,row)) {
1296         continue;
1297       }
1298
1299       // The gain factors
1300       Float_t padgain = calGainFactorDetValue
1301                       * calGainFactorROC->GetValue(col,row);
1302       if (padgain <= 0) {
1303         AliError(Form("Not a valid gain %f, %d %d %d",padgain,det,col,row));
1304       }
1305
1306       for (time = 0; time < nTimeTotal; time++) {
1307
1308         // Get the signal amplitude
1309         Float_t signalAmp = signals->GetData(row,col,time);
1310         // Pad and time coupling
1311         signalAmp *= coupling;
1312         // Gain factors
1313         signalAmp *= padgain;
1314
1315         // Add the noise, starting from minus ADC baseline in electrons
1316         signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,simParam->GetNoise())
1317                                ,-baselineEl);
1318
1319         // Convert to mV
1320         signalAmp *= convert;
1321         // Add ADC baseline in mV
1322         signalAmp += baseline;
1323
1324         // Convert to ADC counts. Set the overflow-bit fADCoutRange if the
1325         // signal is larger than fADCinRange
1326         Short_t adc  = 0;
1327         if (signalAmp >= simParam->GetADCinRange()) {
1328           adc = ((Short_t) simParam->GetADCoutRange());
1329         }
1330         else {
1331           adc = TMath::Nint(signalAmp * adcConvert);
1332         }
1333
1334         // Saving all digits
1335         digits->SetData(row,col,time,adc);
1336
1337       } // for: time
1338
1339     } // for: col
1340   } // for: row
1341
1342   return kTRUE;
1343
1344 }
1345
1346 //_____________________________________________________________________________
1347 Bool_t AliTRDdigitizer::Signal2SDigits(Int_t det, AliTRDarraySignal *signals)
1348 {
1349   //
1350   // Converts the sampled electron signals to s-digits
1351   //
1352
1353   AliDebug(1,Form("Start converting signals to s-digits for detector=%d",det));
1354
1355   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1356   if (!calibration) {
1357     AliFatal("Could not get calibration object");
1358     return kFALSE;
1359   }
1360
1361   Int_t row  = 0;
1362   Int_t col  = 0;
1363   Int_t time = 0;
1364
1365   Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
1366   Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
1367   Int_t nTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
1368
1369   // Get the container for the digits of this detector
1370   if (!fDigitsManager->HasSDigits()) {
1371     AliError("Digits manager has no s-digits");
1372     return kFALSE;
1373   }
1374
1375   AliTRDarraySignal *digits = (AliTRDarraySignal *) fDigitsManager->GetSDigits(det);
1376   // Allocate memory space for the digits buffer
1377   if (!digits->HasData()) {
1378     digits->Allocate(nRowMax,nColMax,nTimeTotal);
1379   }
1380
1381   // Create the sdigits for this chamber
1382   for (row  = 0; row  <  nRowMax; row++ ) {
1383     for (col  = 0; col  <  nColMax; col++ ) {
1384       for (time = 0; time < nTimeTotal; time++) {
1385         digits->SetData(row,col,time,signals->GetData(row,col,time));
1386       } // for: time
1387     } // for: col
1388   } // for: row
1389
1390   return kTRUE;
1391
1392 }
1393
1394 //_____________________________________________________________________________
1395 Bool_t AliTRDdigitizer::Digits2SDigits(AliTRDdigitsManager * const manDig
1396                                      , AliTRDdigitsManager * const manSDig)
1397 {
1398   //
1399   // Converts digits into s-digits. Needed for embedding into real data.
1400   //
1401
1402   AliDebug(1,"Start converting digits to s-digits");
1403
1404   if (!fGeo) {
1405     fGeo = new AliTRDgeometry();
1406   }
1407
1408   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1409   if (!calibration) {
1410     AliFatal("Could not get calibration object");
1411     return kFALSE;
1412   }
1413
1414   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
1415   if (!simParam) {
1416     AliFatal("Could not get simulation parameters");
1417     return kFALSE;
1418   }
1419
1420   // Converts number of electrons to fC
1421   const Double_t kEl2fC = 1.602e-19 * 1.0e15;
1422
1423   // Coupling factor
1424   Double_t coupling     = simParam->GetPadCoupling()
1425                         * simParam->GetTimeCoupling();
1426   // Electronics conversion factor
1427   Double_t convert      = kEl2fC
1428                         * simParam->GetChipGain();
1429   // ADC conversion factor
1430   Double_t adcConvert   = simParam->GetADCoutRange()
1431                         / simParam->GetADCinRange();
1432   // The electronics baseline in mV
1433   Double_t baseline     = simParam->GetADCbaseline()
1434                         / adcConvert;
1435   // The electronics baseline in electrons
1436   //Double_t baselineEl   = baseline
1437   //                      / convert;
1438
1439   // The gainfactor calibration objects
1440   //const AliTRDCalDet *calGainFactorDet      = calibration->GetGainFactorDet();
1441   //AliTRDCalROC       *calGainFactorROC      = 0;
1442   //Float_t             calGainFactorDetValue = 0.0;
1443
1444   Int_t row  = 0;
1445   Int_t col  = 0;
1446   Int_t time = 0;
1447
1448   for (Int_t det = 0; det < AliTRDgeometry::Ndet(); det++) {
1449
1450     Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
1451     Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
1452     Int_t nTimeTotal = manDig->GetDigitsParam()->GetNTimeBins(det);
1453
1454     // Get the calibration objects
1455     //calGainFactorROC      = calibration->GetGainFactorROC(det);
1456     //calGainFactorDetValue = calGainFactorDet->GetValue(det);
1457
1458     // Get the digits
1459     AliTRDarrayADC *digits = (AliTRDarrayADC *) manDig->GetDigits(det);
1460
1461     if (!manSDig->HasSDigits()) {
1462       AliError("SDigits manager has no s-digits");
1463       return kFALSE;
1464     }
1465     // Get the s-digits
1466     AliTRDarraySignal     *sdigits = (AliTRDarraySignal *)     manSDig->GetSDigits(det);
1467     AliTRDarrayDictionary *tracks0 = (AliTRDarrayDictionary *) manSDig->GetDictionary(det,0);
1468     AliTRDarrayDictionary *tracks1 = (AliTRDarrayDictionary *) manSDig->GetDictionary(det,1);
1469     AliTRDarrayDictionary *tracks2 = (AliTRDarrayDictionary *) manSDig->GetDictionary(det,2);
1470     // Allocate memory space for the digits buffer
1471     sdigits->Allocate(nRowMax,nColMax,nTimeTotal);
1472     tracks0->Allocate(nRowMax,nColMax,nTimeTotal);
1473     tracks1->Allocate(nRowMax,nColMax,nTimeTotal);
1474     tracks2->Allocate(nRowMax,nColMax,nTimeTotal);
1475
1476     // Keep the digits param
1477     manSDig->GetDigitsParam()->SetNTimeBinsAll(manDig->GetDigitsParam()->GetNTimeBins(0));
1478     manSDig->GetDigitsParam()->SetADCbaselineAll(manDig->GetDigitsParam()->GetADCbaseline(0));
1479
1480     if (digits->HasData()) {
1481
1482       digits->Expand();
1483
1484       // Create the sdigits for this chamber
1485       for (row  = 0; row  <  nRowMax; row++ ) {
1486         for (col  = 0; col  <  nColMax; col++ ) {
1487
1488           // The gain factors
1489           //Float_t padgain = calGainFactorDetValue
1490           //                * calGainFactorROC->GetValue(col,row);
1491
1492           for (time = 0; time < nTimeTotal; time++) {
1493
1494             Short_t  adcVal = digits->GetData(row,col,time);
1495             Double_t signal = (Double_t) adcVal;
1496             // ADC -> signal in mV
1497             signal /= adcConvert;
1498             // Subtract baseline in mV
1499             signal -= baseline;
1500             // Signal in mV -> signal in #electrons
1501             signal /= convert;
1502             // Gain factor
1503             //signal /= padgain; // Not needed for real data
1504             // Pad and time coupling
1505             signal /= coupling;
1506
1507             sdigits->SetData(row,col,time,signal);
1508             tracks0->SetData(row,col,time,0);
1509             tracks1->SetData(row,col,time,0);
1510             tracks2->SetData(row,col,time,0);
1511
1512           } // for: time
1513
1514         } // for: col
1515       } // for: row
1516
1517     } // if: has data
1518
1519     sdigits->Compress(0);
1520     tracks0->Compress();
1521     tracks1->Compress();
1522     tracks2->Compress();
1523
1524     // No compress just remove
1525     manDig->RemoveDigits(det);
1526     manDig->RemoveDictionaries(det);
1527
1528   } // for: det
1529
1530   return kTRUE;
1531
1532 }
1533
1534 //_____________________________________________________________________________
1535 Bool_t AliTRDdigitizer::SDigits2Digits()
1536 {
1537   //
1538   // Merges the input s-digits and converts them to normal digits
1539   //
1540
1541   if (!MergeSDigits()) {
1542     return kFALSE;
1543   }
1544
1545   return ConvertSDigits();
1546
1547 }
1548
1549 //_____________________________________________________________________________
1550 Bool_t AliTRDdigitizer::MergeSDigits()
1551 {
1552   //
1553   // Merges the input s-digits:
1554   //   - The amplitude of the different inputs are summed up.
1555   //   - Of the track IDs from the input dictionaries only one is
1556   //     kept for each input. This works for maximal 3 different merged inputs.
1557   //
1558
1559   // Number of track dictionary arrays
1560   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1561
1562   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
1563   if (!simParam) {
1564     AliFatal("Could not get simulation parameters");
1565     return kFALSE;
1566   }
1567
1568   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1569   if (!calibration) {
1570     AliFatal("Could not get calibration object");
1571     return kFALSE;
1572   }
1573
1574   Int_t iDict = 0;
1575   Int_t jDict = 0;
1576
1577   AliTRDarraySignal     *digitsA;
1578   AliTRDarraySignal     *digitsB;
1579   AliTRDarrayDictionary *dictionaryA[kNDict];
1580   AliTRDarrayDictionary *dictionaryB[kNDict];
1581
1582   AliTRDdigitsManager   *mergeSDigitsManager = 0x0;
1583   // Get the first s-digits
1584   fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
1585   if (!fSDigitsManager) {
1586     AliError("No SDigits manager");
1587     return kFALSE;
1588   }
1589
1590   // Loop through the other sets of s-digits
1591   mergeSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->After(fSDigitsManager);
1592
1593   if (mergeSDigitsManager) {
1594     AliDebug(1,Form("Merge %d input files.",fSDigitsManagerList->GetSize()));
1595   }
1596   else {
1597     AliDebug(1,"Only one input file.");
1598   }
1599
1600   Int_t iMerge = 0;
1601
1602   while (mergeSDigitsManager) {
1603
1604     iMerge++;
1605
1606     // Loop through the detectors
1607     for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1608
1609       Int_t nTimeTotal = fSDigitsManager->GetDigitsParam()->GetNTimeBins(iDet);
1610       if (mergeSDigitsManager->GetDigitsParam()->GetNTimeBins(iDet) != nTimeTotal) {
1611         AliError(Form("Mismatch in the number of time bins [%d,%d] in detector %d"
1612                      ,nTimeTotal
1613                      ,mergeSDigitsManager->GetDigitsParam()->GetNTimeBins(iDet)
1614                      ,iDet));
1615         return kFALSE;
1616       }
1617
1618       Int_t nRowMax = fGeo->GetPadPlane(iDet)->GetNrows();
1619       Int_t nColMax = fGeo->GetPadPlane(iDet)->GetNcols();
1620
1621       // Loop through the pixels of one detector and add the signals
1622       digitsA = (AliTRDarraySignal *) fSDigitsManager->GetSDigits(iDet);
1623       digitsB = (AliTRDarraySignal *) mergeSDigitsManager->GetSDigits(iDet);
1624       digitsA->Expand();
1625       if (!digitsA->HasData()) continue;
1626       digitsB->Expand();
1627       if (!digitsB->HasData()) continue;
1628
1629       for (iDict = 0; iDict < kNDict; iDict++) {
1630         dictionaryA[iDict] = (AliTRDarrayDictionary *) fSDigitsManager->GetDictionary(iDet,iDict);
1631         dictionaryB[iDict] = (AliTRDarrayDictionary *) mergeSDigitsManager->GetDictionary(iDet,iDict);
1632         dictionaryA[iDict]->Expand();
1633         dictionaryB[iDict]->Expand();
1634       }
1635
1636       // Merge only detectors that contain a signal
1637       Bool_t doMerge = kTRUE;
1638       if (fMergeSignalOnly) {
1639         if (digitsA->GetOverThreshold(0) == 0) {
1640           doMerge = kFALSE;
1641         }
1642       }
1643
1644       if (doMerge) {
1645
1646         AliDebug(1,Form("Merge detector %d of input no.%d",iDet,iMerge+1));
1647
1648         for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1649           for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1650             for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
1651
1652               // Add the amplitudes of the summable digits
1653               Float_t ampA = digitsA->GetData(iRow,iCol,iTime);
1654               Float_t ampB = digitsB->GetData(iRow,iCol,iTime);
1655               ampA += ampB;
1656               digitsA->SetData(iRow,iCol,iTime,ampA);
1657
1658               // Add the mask to the track id if defined.
1659               for (iDict = 0; iDict < kNDict; iDict++) {
1660                 Int_t trackB = dictionaryB[iDict]->GetData(iRow,iCol,iTime);
1661                 if ((fMasks) && (trackB > 0))  {
1662                   for (jDict = 0; jDict < kNDict; jDict++) {
1663                     Int_t trackA = dictionaryA[iDict]->GetData(iRow,iCol,iTime);
1664                     if (trackA == 0) {
1665                       trackA = trackB + fMasks[iMerge];
1666                       dictionaryA[iDict]->SetData(iRow,iCol,iTime,trackA);
1667                     } // if:  track A == 0
1668                   } // for: jDict
1669                 } // if:  fMasks and trackB > 0
1670               } // for: iDict
1671
1672             } // for: iTime
1673           } // for: iCol
1674         } // for: iRow
1675
1676       } // if:  doMerge
1677
1678       mergeSDigitsManager->RemoveDigits(iDet);
1679       mergeSDigitsManager->RemoveDictionaries(iDet);
1680
1681       if (fCompress) {
1682         digitsA->Compress(0);
1683         for (iDict = 0; iDict < kNDict; iDict++) {
1684           dictionaryA[iDict]->Compress();
1685         }
1686       }
1687
1688     } // for: detectors
1689
1690     // The next set of s-digits
1691     mergeSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->After(mergeSDigitsManager);
1692
1693   } // while: mergeDigitsManagers
1694
1695   return kTRUE;
1696
1697 }
1698
1699 //_____________________________________________________________________________
1700 Bool_t AliTRDdigitizer::ConvertSDigits()
1701 {
1702   //
1703   // Converts s-digits to normal digits
1704   //
1705
1706   AliTRDarraySignal *digitsIn = 0x0;
1707
1708   if (!fSDigitsManager->HasSDigits()) {
1709     AliError("No s-digits in digits manager");
1710     return kFALSE;
1711   }
1712
1713   // Loop through the detectors
1714   for (Int_t det = 0; det < AliTRDgeometry::Ndet(); det++) {
1715
1716     // Get the merged s-digits (signals)
1717     digitsIn = (AliTRDarraySignal *) fSDigitsManager->GetSDigits(det);
1718     if (!digitsIn->HasData()) {
1719       AliDebug(2,Form("No digits for det=%d",det));
1720       continue;
1721     }
1722
1723     // Convert the merged sdigits to digits
1724     if (!Signal2ADC(det,digitsIn)) {
1725       continue;
1726     }
1727
1728     // Copy the dictionary information to the output array
1729     if (!CopyDictionary(det)) {
1730       continue;
1731     }
1732
1733     // Delete
1734     fSDigitsManager->RemoveDigits(det);
1735     fSDigitsManager->RemoveDictionaries(det);
1736
1737     // Run digital processing
1738     RunDigitalProcessing(det);
1739
1740     // Compress the arrays
1741     CompressOutputArrays(det);
1742
1743   } // for: detector numbers
1744
1745   if (AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
1746     if (trklLoader->Tree())
1747       trklLoader->WriteData("OVERWRITE");
1748   }
1749
1750   // Save the values for the raw data headers
1751   if (AliTRDSimParam::Instance()->GetNTBoverwriteOCDB()) {
1752     fDigitsManager->GetDigitsParam()->SetNTimeBinsAll(AliTRDSimParam::Instance()->GetNTimeBins());
1753   }
1754   else {
1755     fDigitsManager->GetDigitsParam()->SetNTimeBinsAll(AliTRDcalibDB::Instance()->GetNumberOfTimeBinsDCS());
1756   }
1757   fDigitsManager->GetDigitsParam()->SetADCbaselineAll(AliTRDSimParam::Instance()->GetADCbaseline());
1758
1759   return kTRUE;
1760
1761 }
1762
1763 //_____________________________________________________________________________
1764 Bool_t AliTRDdigitizer::CopyDictionary(Int_t det)
1765 {
1766   //
1767   // Copies the dictionary information from the s-digits arrays
1768   // to the output arrays
1769   //
1770
1771   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1772   if (!calibration) {
1773     AliFatal("Could not get calibration object");
1774     return kFALSE;
1775   }
1776
1777   AliDebug(1,Form("Start copying dictionaries for detector=%d",det));
1778
1779   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1780   AliTRDarrayDictionary *dictionaryIn[kNDict];
1781   AliTRDarrayDictionary *dictionaryOut[kNDict];
1782
1783   Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
1784   Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
1785   Int_t nTimeTotal = fSDigitsManager->GetDigitsParam()->GetNTimeBins(det);
1786
1787   Int_t row  = 0;
1788   Int_t col  = 0;
1789   Int_t time = 0;
1790   Int_t dict = 0;
1791
1792   for (dict = 0; dict < kNDict; dict++) {
1793
1794     dictionaryIn[dict]  = (AliTRDarrayDictionary *) fSDigitsManager->GetDictionary(det,dict);
1795     dictionaryIn[dict]->Expand();
1796     dictionaryOut[dict] = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(det,dict);
1797     dictionaryOut[dict]->Allocate(nRowMax,nColMax,nTimeTotal);
1798
1799     for (row = 0; row < nRowMax; row++) {
1800       for (col = 0; col < nColMax; col++) {
1801         for (time = 0; time < nTimeTotal; time++) {
1802           Int_t track = dictionaryIn[dict]->GetData(row,col,time);
1803           dictionaryOut[dict]->SetData(row,col,time,track);
1804         } // for: time
1805       } // for: col
1806     } // for: row
1807
1808   } // for: dictionaries
1809
1810   return kTRUE;
1811
1812 }
1813
1814 //_____________________________________________________________________________
1815 void AliTRDdigitizer::CompressOutputArrays(Int_t det)
1816 {
1817   //
1818   // Compress the output arrays
1819   //
1820
1821   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1822   AliTRDarrayDictionary *dictionary = 0x0;
1823
1824   if (fCompress) {
1825
1826     if (!fSDigits) {
1827       AliTRDarrayADC *digits = 0x0;
1828       digits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det);
1829       digits->Compress();
1830     }
1831
1832     if (fSDigits) {
1833       AliTRDarraySignal *digits = 0x0;
1834       digits = (AliTRDarraySignal *) fDigitsManager->GetSDigits(det);
1835       digits->Compress(0);
1836     }
1837
1838     for (Int_t dict = 0; dict < kNDict; dict++) {
1839       dictionary = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(det,dict);
1840       dictionary->Compress();
1841     }
1842
1843   }
1844
1845 }
1846
1847 //_____________________________________________________________________________
1848 Bool_t AliTRDdigitizer::WriteDigits() const
1849 {
1850   //
1851   // Writes out the TRD-digits and the dictionaries
1852   //
1853
1854   // Write parameters
1855   fRunLoader->CdGAFile();
1856
1857   // Store the digits and the dictionary in the tree
1858   return fDigitsManager->WriteDigits();
1859
1860 }
1861
1862 //_____________________________________________________________________________
1863 void AliTRDdigitizer::InitOutput(Int_t iEvent)
1864 {
1865   //
1866   // Initializes the output branches
1867   //
1868
1869   fEvent = iEvent;
1870
1871   if (!fRunLoader) {
1872     AliError("Run Loader is NULL");
1873     return;
1874   }
1875
1876   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
1877   if (!loader) {
1878     AliError("Can not get TRD loader from Run Loader");
1879     return;
1880   }
1881
1882   TTree *tree = 0;
1883
1884   if (fSDigits) {
1885     // If we produce SDigits
1886     tree = loader->TreeS();
1887     if (!tree) {
1888       loader->MakeTree("S");
1889       tree = loader->TreeS();
1890     }
1891   }
1892   else {
1893     // If we produce Digits
1894     tree = loader->TreeD();
1895     if (!tree) {
1896       loader->MakeTree("D");
1897       tree = loader->TreeD();
1898     }
1899   }
1900   fDigitsManager->SetEvent(iEvent);
1901   fDigitsManager->MakeBranch(tree);
1902
1903 }
1904
1905 //_____________________________________________________________________________
1906 Int_t AliTRDdigitizer::Diffusion(Float_t vdrift, Double_t absdriftlength
1907                                , Double_t exbvalue
1908                                , Double_t &lRow, Double_t &lCol, Double_t &lTime)
1909 {
1910   //
1911   // Applies the diffusion smearing to the position of a single electron.
1912   // Depends on absolute drift length.
1913   //
1914
1915   Float_t diffL = 0.0;
1916   Float_t diffT = 0.0;
1917
1918   if (AliTRDCommonParam::Instance()->GetDiffCoeff(diffL,diffT,vdrift)) {
1919
1920     Float_t driftSqrt = TMath::Sqrt(absdriftlength);
1921     Float_t sigmaT    = driftSqrt * diffT;
1922     Float_t sigmaL    = driftSqrt * diffL;
1923     lRow  = gRandom->Gaus(lRow ,sigmaT);
1924     if (AliTRDCommonParam::Instance()->ExBOn()) {
1925       lCol  = gRandom->Gaus(lCol ,sigmaT * 1.0 / (1.0 + exbvalue*exbvalue));
1926       lTime = gRandom->Gaus(lTime,sigmaL * 1.0 / (1.0 + exbvalue*exbvalue));
1927     }
1928     else {
1929       lCol  = gRandom->Gaus(lCol ,sigmaT);
1930       lTime = gRandom->Gaus(lTime,sigmaL);
1931     }
1932
1933     return 1;
1934
1935   }
1936   else {
1937
1938     return 0;
1939
1940   }
1941
1942 }
1943
1944 //_____________________________________________________________________________
1945 void AliTRDdigitizer::RunDigitalProcessing(Int_t det)
1946 {
1947   //
1948   // Run the digital processing in the TRAP
1949   //
1950
1951   AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
1952
1953   AliTRDarrayADC *digits = fDigitsManager->GetDigits(det);
1954   if (!digits)
1955     return;
1956
1957   //Call the methods in the mcm class using the temporary array as input
1958   for(Int_t rob = 0; rob < digits->GetNrow() / 2; rob++)
1959   {
1960     for(Int_t mcm = 0; mcm < 16; mcm++)
1961     {
1962       fMcmSim->Init(det, rob, mcm);
1963       fMcmSim->SetDataByPad(digits, fDigitsManager);
1964       fMcmSim->Filter();
1965       if (feeParam->GetTracklet()) {
1966         fMcmSim->Tracklet();
1967         fMcmSim->StoreTracklets();
1968       }
1969       fMcmSim->ZSMapping();
1970       fMcmSim->WriteData(digits);
1971     }
1972   }
1973
1974 }
1975