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