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 //  Authors: C. Blume (blume@ikf.uni-frankfurt.de)                           //
  22 //           C. Lippmann                                                     //
  23 //           B. Vulpescu                                                     //
  24 //                                                                           //
  25 //  The following effects are included:                                      //
  26 //      - Diffusion                                                          //
  27 //      - ExB effects                                                        //
  28 //      - Gas gain including fluctuations                                    //
  29 //      - Pad-response (simple Gaussian approximation)                       //
  30 //      - Time-response                                                      //
  31 //      - Electronics noise                                                  //
  32 //      - Electronics gain                                                   //
  33 //      - Digitization                                                       //
  34 //      - ADC threshold                                                      //
  35 //  The corresponding parameter can be adjusted via the various              //
  36 //  Set-functions. If these parameters are not explicitly set, default       //
  37 //  values are used (see Init-function).                                     //
  38 //                                                                           //
  39 ///////////////////////////////////////////////////////////////////////////////
  40
  41 #include <stdlib.h>
  42
  43 #include <TMath.h>
  44 #include <TVector.h>
  45 #include <TRandom.h>
  46 #include <TROOT.h>
  47 #include <TTree.h>
  48 #include <TFile.h>
  49 #include <TF1.h>
  50 #include <TList.h>
  51 #include <TTask.h>
  52 #include <TGeoManager.h>
  53
  54 #include "AliRun.h"
  55 #include "AliRunLoader.h"
  56 #include "AliLoader.h"
  57 #include "AliConfig.h"
  58 #include "AliMagF.h"
  59 #include "AliRunDigitizer.h"
  60 #include "AliRunLoader.h"
  61 #include "AliLoader.h"
  62 #include "AliLog.h"
  63
  64 #include "AliTRD.h"
  65 #include "AliTRDhit.h"
  66 #include "AliTRDdigitizer.h"
  67 #include "AliTRDdataArrayI.h"
  68 #include "AliTRDdataArrayF.h"
  69 #include "AliTRDsegmentArray.h"
  70 #include "AliTRDdigitsManager.h"
  71 #include "AliTRDgeometry.h"
  72 #include "AliTRDpadPlane.h"
  73 #include "AliTRDcalibDB.h"
  74 #include "AliTRDSimParam.h"
  75 #include "AliTRDCommonParam.h"
  76
  77 #include "Cal/AliTRDCalROC.h"
  78 #include "Cal/AliTRDCalDet.h"
  79
  80 ClassImp(AliTRDdigitizer)
  81
  82 //_____________________________________________________________________________
  83 AliTRDdigitizer::AliTRDdigitizer()
  84   :AliDigitizer()
  85   ,fRunLoader(0)
  86   ,fDigitsManager(0)
  87   ,fSDigitsManager(0)
  88   ,fSDigitsManagerList(0)
  89   ,fTRD(0)
  90   ,fGeo(0)
  91   ,fEvent(0)
  92   ,fMasks(0)
  93   ,fCompress(kTRUE)
  94   ,fSDigits(kFALSE)
  95   ,fSDigitsScale(0.0)
  96   ,fMergeSignalOnly(kFALSE)
  97   ,fDiffLastVdrift(0)
  98   ,fDiffusionT(0.0)
  99   ,fDiffusionL(0.0)
 100   ,fOmegaTau(0.0)
 101   ,fLorentzFactor(0.0)
 102   ,fTimeLastVdrift(0)
 103   ,fTimeStruct1(0)
 104   ,fTimeStruct2(0)
 105   ,fVDlo(0)
 106   ,fVDhi(0)
 107 {
 108   //
 109   // AliTRDdigitizer default constructor
 110   //
 111
 112   Init();
 113
 114 }
 115
 116 //_____________________________________________________________________________
 117 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
 118   :AliDigitizer(name,title)
 119   ,fRunLoader(0)
 120   ,fDigitsManager(0)
 121   ,fSDigitsManager(0)
 122   ,fSDigitsManagerList(0)
 123   ,fTRD(0)
 124   ,fGeo(0)
 125   ,fEvent(0)
 126   ,fMasks(0)
 127   ,fCompress(kTRUE)
 128   ,fSDigits(kFALSE)
 129   ,fSDigitsScale(0.0)
 130   ,fMergeSignalOnly(kFALSE)
 131   ,fDiffLastVdrift(0)
 132   ,fDiffusionT(0.0)
 133   ,fDiffusionL(0.0)
 134   ,fOmegaTau(0.0)
 135   ,fLorentzFactor(0.0)
 136   ,fTimeLastVdrift(0)
 137   ,fTimeStruct1(0)
 138   ,fTimeStruct2(0)
 139   ,fVDlo(0)
 140   ,fVDhi(0)
 141 {
 142   //
 143   // AliTRDdigitizer constructor
 144   //
 145
 146   Init();
 147
 148 }
 149
 150 //_____________________________________________________________________________
 151 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager
 152                                , const Text_t *name, const Text_t *title)
 153   :AliDigitizer(manager,name,title)
 154   ,fRunLoader(0)
 155   ,fDigitsManager(0)
 156   ,fSDigitsManager(0)
 157   ,fSDigitsManagerList(0)
 158   ,fTRD(0)
 159   ,fGeo(0)
 160   ,fEvent(0)
 161   ,fMasks(0)
 162   ,fCompress(kTRUE)
 163   ,fSDigits(kFALSE)
 164   ,fSDigitsScale(0.0)
 165   ,fMergeSignalOnly(kFALSE)
 166   ,fDiffLastVdrift(0)
 167   ,fDiffusionT(0.0)
 168   ,fDiffusionL(0.0)
 169   ,fOmegaTau(0.0)
 170   ,fLorentzFactor(0.0)
 171   ,fTimeLastVdrift(0)
 172   ,fTimeStruct1(0)
 173   ,fTimeStruct2(0)
 174   ,fVDlo(0)
 175   ,fVDhi(0)
 176 {
 177   //
 178   // AliTRDdigitizer constructor
 179   //
 180
 181   Init();
 182
 183 }
 184
 185 //_____________________________________________________________________________
 186 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager)
 187   :AliDigitizer(manager,"AliTRDdigitizer","TRD digitizer")
 188   ,fRunLoader(0)
 189   ,fDigitsManager(0)
 190   ,fSDigitsManager(0)
 191   ,fSDigitsManagerList(0)
 192   ,fTRD(0)
 193   ,fGeo(0)
 194   ,fEvent(0)
 195   ,fMasks(0)
 196   ,fCompress(kTRUE)
 197   ,fSDigits(kFALSE)
 198   ,fSDigitsScale(0.0)
 199   ,fMergeSignalOnly(kFALSE)
 200   ,fDiffLastVdrift(0)
 201   ,fDiffusionT(0.0)
 202   ,fDiffusionL(0.0)
 203   ,fOmegaTau(0.0)
 204   ,fLorentzFactor(0.0)
 205   ,fTimeLastVdrift(0)
 206   ,fTimeStruct1(0)
 207   ,fTimeStruct2(0)
 208   ,fVDlo(0)
 209   ,fVDhi(0)
 210 {
 211   //
 212   // AliTRDdigitizer constructor
 213   //
 214
 215   Init();
 216
 217 }
 218
 219 //_____________________________________________________________________________
 220 Bool_t AliTRDdigitizer::Init()
 221 {
 222   //
 223   // Initialize the digitizer with default values
 224   //
 225
 226   fRunLoader          = 0;
 227   fDigitsManager      = 0;
 228   fSDigitsManager     = 0;
 229   fSDigitsManagerList = 0;
 230   fTRD                = 0;
 231   fGeo                = 0;
 232
 233   fEvent              = 0;
 234   fMasks              = 0;
 235   fCompress           = kTRUE;
 236   fSDigits            = kFALSE;
 237   fSDigitsScale       = 100.0;
 238   fMergeSignalOnly    = kFALSE;
 239
 240   fTimeLastVdrift     = -1;
 241   fTimeStruct1        =  0;
 242   fTimeStruct2        =  0;
 243   fVDlo               =  0;
 244   fVDhi               =  0;
 245
 246   fDiffLastVdrift     = -1;
 247   fDiffusionT         =  0.0;
 248   fDiffusionL         =  0.0;
 249   fLorentzFactor      =  0.0;
 250
 251   return AliDigitizer::Init();
 252
 253 }
 254
 255 //_____________________________________________________________________________
 256 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
 257   :AliDigitizer(d)
 258   ,fRunLoader(0)
 259   ,fDigitsManager(0)
 260   ,fSDigitsManager(0)
 261   ,fSDigitsManagerList(0)
 262   ,fTRD(0)
 263   ,fGeo(0)
 264   ,fEvent(0)
 265   ,fMasks(0)
 266   ,fCompress(d.fCompress)
 267   ,fSDigits(d.fSDigits)
 268   ,fSDigitsScale(d.fSDigitsScale)
 269   ,fMergeSignalOnly(d.fMergeSignalOnly)
 270   ,fDiffLastVdrift(-1)
 271   ,fDiffusionT(0.0)
 272   ,fDiffusionL(0.0)
 273   ,fOmegaTau(0.0)
 274   ,fLorentzFactor(0.0)
 275   ,fTimeLastVdrift(-1)
 276   ,fTimeStruct1(0)
 277   ,fTimeStruct2(0)
 278   ,fVDlo(0)
 279   ,fVDhi(0)
 280 {
 281   //
 282   // AliTRDdigitizer copy constructor
 283   //
 284
 285   // Do not copy timestructs, just invalidate lastvdrift.
 286   // Next time they are requested, they get recalculated
 287   if (((AliTRDdigitizer &) d).fTimeStruct1) {
 288     delete [] ((AliTRDdigitizer &) d).fTimeStruct1;
 289     ((AliTRDdigitizer &) d).fTimeStruct1 = 0;
 290   }
 291   if (((AliTRDdigitizer &) d).fTimeStruct2) {
 292     delete [] ((AliTRDdigitizer &) d).fTimeStruct2;
 293     ((AliTRDdigitizer &) d).fTimeStruct2 = 0;
 294   }
 295
 296 }
 297
 298 //_____________________________________________________________________________
 299 AliTRDdigitizer::~AliTRDdigitizer()
 300 {
 301   //
 302   // AliTRDdigitizer destructor
 303   //
 304
 305   if (fDigitsManager) {
 306     delete fDigitsManager;
 307     fDigitsManager      = 0;
 308   }
 309
 310   if (fDigitsManager) {
 311     delete fSDigitsManager;
 312     fSDigitsManager     = 0;
 313   }
 314
 315   if (fSDigitsManagerList) {
 316     fSDigitsManagerList->Delete();
 317     delete fSDigitsManagerList;
 318     fSDigitsManagerList = 0;
 319   }
 320
 321   if (fMasks) {
 322     delete [] fMasks;
 323     fMasks       = 0;
 324   }
 325
 326   if (fTimeStruct1) {
 327     delete [] fTimeStruct1;
 328     fTimeStruct1 = 0;
 329   }
 330
 331   if (fTimeStruct2) {
 332     delete [] fTimeStruct2;
 333     fTimeStruct2 = 0;
 334   }
 335
 336 }
 337
 338 //_____________________________________________________________________________
 339 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
 340 {
 341   //
 342   // Assignment operator
 343   //
 344
 345   if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
 346
 347   return *this;
 348
 349 }
 350
 351 //_____________________________________________________________________________
 352 void AliTRDdigitizer::Copy(TObject &d) const
 353 {
 354   //
 355   // Copy function
 356   //
 357
 358   ((AliTRDdigitizer &) d).fRunLoader          = 0;
 359   ((AliTRDdigitizer &) d).fDigitsManager      = 0;
 360   ((AliTRDdigitizer &) d).fSDigitsManager     = 0;
 361   ((AliTRDdigitizer &) d).fSDigitsManagerList = 0;
 362   ((AliTRDdigitizer &) d).fTRD                = 0;
 363   ((AliTRDdigitizer &) d).fGeo                = 0;
 364   ((AliTRDdigitizer &) d).fEvent              = 0;
 365   ((AliTRDdigitizer &) d).fMasks              = 0;
 366   ((AliTRDdigitizer &) d).fCompress           = fCompress;
 367   ((AliTRDdigitizer &) d).fSDigits            = fSDigits;
 368   ((AliTRDdigitizer &) d).fSDigitsScale       = fSDigitsScale;
 369   ((AliTRDdigitizer &) d).fMergeSignalOnly    = fMergeSignalOnly;
 370
 371   AliTRDdigitizer& target = (AliTRDdigitizer &) d;
 372
 373   // Do not copy timestructs, just invalidate lastvdrift.
 374   // Next time they are requested, they get recalculated
 375   if (target.fTimeStruct1) {
 376     delete [] target.fTimeStruct1;
 377     target.fTimeStruct1 = 0;
 378   }
 379   if (target.fTimeStruct2) {
 380     delete [] target.fTimeStruct2;
 381     target.fTimeStruct2 = 0;
 382   }
 383   target.fTimeLastVdrift = -1;
 384
 385 }
 386
 387 //_____________________________________________________________________________
 388 void AliTRDdigitizer::Exec(Option_t *option)
 389 {
 390   //
 391   // Executes the merging
 392   //
 393
 394   Int_t iInput;
 395
 396   AliTRDdigitsManager *sdigitsManager;
 397
 398   TString optionString = option;
 399   if (optionString.Contains("deb")) {
 400     AliLog::SetClassDebugLevel("AliTRDdigitizer",1);
 401     AliInfo("Called with debug option");
 402   }
 403
 404   // The AliRoot file is already connected by the manager
 405   AliRunLoader *inrl;
 406
 407   if (gAlice) {
 408     AliDebug(1,"AliRun object found on file.");
 409   }
 410   else {
 411     inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
 412     inrl->LoadgAlice();
 413     gAlice = inrl->GetAliRun();
 414     if (!gAlice) {
 415       AliError("Could not find AliRun object.")
 416       return;
 417     }
 418   }
 419
 420   Int_t nInput = fManager->GetNinputs();
 421   fMasks       = new Int_t[nInput];
 422   for (iInput = 0; iInput < nInput; iInput++) {
 423     fMasks[iInput] = fManager->GetMask(iInput);
 424   }
 425
 426   //
 427   // Initialization
 428   //
 429
 430   AliRunLoader *orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
 431
 432   if (InitDetector()) {
 433
 434     AliLoader *ogime = orl->GetLoader("TRDLoader");
 435
 436     TTree *tree = 0;
 437     if (fSDigits) {
 438       // If we produce SDigits
 439       tree = ogime->TreeS();
 440       if (!tree) {
 441         ogime->MakeTree("S");
 442         tree = ogime->TreeS();
 443       }
 444     }
 445     else {
 446       // If we produce Digits
 447       tree = ogime->TreeD();
 448       if (!tree) {
 449         ogime->MakeTree("D");
 450         tree = ogime->TreeD();
 451       }
 452     }
 453
 454     MakeBranch(tree);
 455
 456   }
 457
 458   for (iInput = 0; iInput < nInput; iInput++) {
 459
 460     AliDebug(1,Form("Add input stream %d",iInput));
 461
 462     // Check if the input tree exists
 463     inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
 464     AliLoader *gime = inrl->GetLoader("TRDLoader");
 465
 466     TTree *treees = gime->TreeS();
 467     if (treees == 0x0) {
 468       if (gime->LoadSDigits()) {
 469         AliError(Form("Error Occured while loading S. Digits for input %d.",iInput));
 470         return;
 471       }
 472       treees = gime->TreeS();
 473     }
 474
 475     if (treees == 0x0) {
 476       AliError(Form("Input stream %d does not exist",iInput));
 477       return;
 478     }
 479
 480     // Read the s-digits via digits manager
 481     sdigitsManager = new AliTRDdigitsManager();
 482     sdigitsManager->SetSDigits(kTRUE);
 483
 484     AliRunLoader *rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
 485     AliLoader *gimme = rl->GetLoader("TRDLoader");
 486     if (!gimme->TreeS()) {
 487       gimme->LoadSDigits();
 488     }
 489     sdigitsManager->ReadDigits(gimme->TreeS());
 490
 491     // Add the s-digits to the input list
 492     AddSDigitsManager(sdigitsManager);
 493
 494   }
 495
 496   // Convert the s-digits to normal digits
 497   AliDebug(1,"Do the conversion");
 498   SDigits2Digits();
 499
 500   // Store the digits
 501   AliDebug(1,"Write the digits");
 502   fDigitsManager->WriteDigits();
 503
 504   // Write parameters
 505   orl->CdGAFile();
 506
 507   AliDebug(1,"Done");
 508
 509   DeleteSDigitsManager();
 510
 511 }
 512
 513 //_____________________________________________________________________________
 514 Bool_t AliTRDdigitizer::Open(const Char_t *file, Int_t nEvent)
 515 {
 516   //
 517   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
 518   //
 519   // Connect the AliRoot file containing Geometry, Kine, and Hits
 520   //
 521
 522   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 523
 524   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
 525   if (!fRunLoader) {
 526     fRunLoader = AliRunLoader::Open(file,evfoldname,"UPDATE");
 527   }
 528   if (!fRunLoader) {
 529     AliError(Form("Can not open session for file %s.",file));
 530     return kFALSE;
 531   }
 532
 533   if (!fRunLoader->GetAliRun()) {
 534     fRunLoader->LoadgAlice();
 535   }
 536   gAlice = fRunLoader->GetAliRun();
 537
 538   if (gAlice) {
 539     AliDebug(1,"AliRun object found on file.");
 540   }
 541   else {
 542     AliError("Could not find AliRun object.");
 543     return kFALSE;
 544   }
 545
 546   fEvent = nEvent;
 547
 548   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
 549   if (!loader) {
 550     AliError("Can not get TRD loader from Run Loader");
 551     return kFALSE;
 552   }
 553
 554   if (InitDetector()) {
 555     TTree *tree = 0;
 556     if (fSDigits) {
 557       // If we produce SDigits
 558       tree = loader->TreeS();
 559       if (!tree) {
 560         loader->MakeTree("S");
 561         tree = loader->TreeS();
 562       }
 563     }
 564     else {
 565       // If we produce Digits
 566       tree = loader->TreeD();
 567       if (!tree) {
 568         loader->MakeTree("D");
 569         tree = loader->TreeD();
 570       }
 571     }
 572     return MakeBranch(tree);
 573   }
 574   else {
 575     return kFALSE;
 576   }
 577
 578 }
 579
 580 //_____________________________________________________________________________
 581 Bool_t AliTRDdigitizer::Open(AliRunLoader *runLoader, Int_t nEvent)
 582 {
 583   //
 584   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
 585   //
 586   // Connect the AliRoot file containing Geometry, Kine, and Hits
 587   //
 588
 589   fRunLoader = runLoader;
 590   if (!fRunLoader) {
 591     AliError("RunLoader does not exist");
 592     return kFALSE;
 593   }
 594
 595   if (!fRunLoader->GetAliRun()) {
 596     fRunLoader->LoadgAlice();
 597   }
 598   gAlice = fRunLoader->GetAliRun();
 599
 600   if (gAlice) {
 601     AliDebug(1,"AliRun object found on file.");
 602   }
 603   else {
 604     AliError("Could not find AliRun object.");
 605     return kFALSE;
 606   }
 607
 608   fEvent = nEvent;
 609
 610   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
 611   if (!loader) {
 612     AliError("Can not get TRD loader from Run Loader");
 613     return kFALSE;
 614   }
 615
 616   if (InitDetector()) {
 617     TTree *tree = 0;
 618     if (fSDigits) {
 619       // If we produce SDigits
 620       tree = loader->TreeS();
 621       if (!tree) {
 622         loader->MakeTree("S");
 623         tree = loader->TreeS();
 624       }
 625     }
 626     else {
 627       // If we produce Digits
 628       tree = loader->TreeD();
 629       if (!tree) {
 630         loader->MakeTree("D");
 631         tree = loader->TreeD();
 632       }
 633     }
 634     return MakeBranch(tree);
 635   }
 636   else {
 637     return kFALSE;
 638   }
 639
 640 }
 641
 642 //_____________________________________________________________________________
 643 Bool_t AliTRDdigitizer::InitDetector()
 644 {
 645   //
 646   // Sets the pointer to the TRD detector and the geometry
 647   //
 648
 649   // Get the pointer to the detector class and check for version 1
 650   fTRD = (AliTRD *) gAlice->GetDetector("TRD");
 651   if (!fTRD) {
 652     AliFatal("No TRD module found");
 653     exit(1);
 654   }
 655   if (fTRD->IsVersion() != 1) {
 656     AliFatal("TRD must be version 1 (slow simulator)");
 657     exit(1);
 658   }
 659
 660   // Get the geometry
 661   fGeo = fTRD->GetGeometry();
 662   AliDebug(1,Form("Geometry version %d",fGeo->IsVersion()));
 663
 664   // Create a digits manager
 665   delete fDigitsManager;
 666   fDigitsManager = new AliTRDdigitsManager();
 667   fDigitsManager->SetSDigits(fSDigits);
 668   fDigitsManager->CreateArrays();
 669   fDigitsManager->SetEvent(fEvent);
 670
 671   // The list for the input s-digits manager to be merged
 672   if (fSDigitsManagerList) {
 673     fSDigitsManagerList->Delete();
 674   }
 675   else {
 676     fSDigitsManagerList = new TList();
 677   }
 678
 679   return kTRUE;
 680
 681 }
 682
 683 //_____________________________________________________________________________
 684 Bool_t AliTRDdigitizer::MakeBranch(TTree *tree) const
 685 {
 686   //
 687   // Create the branches for the digits array
 688   //
 689
 690   return fDigitsManager->MakeBranch(tree);
 691
 692 }
 693
 694 //_____________________________________________________________________________
 695 Bool_t AliTRDdigitizer::MakeDigits()
 696 {
 697   //
 698   // Creates digits.
 699   //
 700
 701   ///////////////////////////////////////////////////////////////
 702   // Parameter
 703   ///////////////////////////////////////////////////////////////
 704
 705   // Converts number of electrons to fC
 706   const Double_t kEl2fC  = 1.602e-19 * 1.0e15;
 707
 708   ///////////////////////////////////////////////////////////////
 709
 710   // Number of pads included in the pad response
 711   const Int_t kNpad      = 3;
 712
 713   // Number of track dictionary arrays
 714   const Int_t kNDict     = AliTRDdigitsManager::kNDict;
 715
 716   // Width of the amplification region
 717   const Float_t kAmWidth = AliTRDgeometry::AmThick();
 718   // Width of the drift region
 719   const Float_t kDrWidth = AliTRDgeometry::DrThick();
 720   // Drift + amplification region
 721   const Float_t kDrMin   =          - 0.5 * kAmWidth;
 722   const Float_t kDrMax   = kDrWidth + 0.5 * kAmWidth;
 723
 724   Int_t    iRow;
 725   Int_t    iCol;
 726   Int_t    iTime;
 727   Int_t    iPad;
 728   Int_t    iDict  = 0;
 729   Int_t    nBytes = 0;
 730
 731   Int_t    totalSizeDigits = 0;
 732   Int_t    totalSizeDict0  = 0;
 733   Int_t    totalSizeDict1  = 0;
 734   Int_t    totalSizeDict2  = 0;
 735
 736   Int_t    timeBinTRFend   = 1;
 737
 738   Double_t pos[3];
 739   Double_t loc[3];
 740   Double_t padSignal[kNpad];
 741   Double_t signalOld[kNpad];
 742
 743   AliTRDdataArrayF *signals  = 0;
 744   AliTRDdataArrayI *digits   = 0;
 745   AliTRDdataArrayI *dictionary[kNDict];
 746
 747   AliTRDpadPlane   *padPlane = 0;
 748
 749   AliTRDCalROC     *calVdriftROC          = 0;
 750   Float_t           calVdriftDetValue     = 0.0;
 751   AliTRDCalROC     *calT0ROC              = 0;
 752   Float_t           calT0DetValue         = 0.0;
 753   AliTRDCalROC     *calGainFactorROC      = 0;
 754   Float_t           calGainFactorDetValue = 0.0;
 755
 756   if (!gGeoManager) {
 757     AliFatal("No geometry manager!");
 758   }
 759
 760   if (!fGeo) {
 761     AliError("No geometry defined");
 762     return kFALSE;
 763   }
 764   fGeo->ReadGeoMatrices();
 765
 766   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 767   if (!simParam) {
 768     AliFatal("Could not get simulation parameters");
 769     return kFALSE;
 770   }
 771
 772   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
 773   if (!commonParam) {
 774     AliFatal("Could not get common parameterss");
 775     return kFALSE;
 776   }
 777
 778   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
 779   if (!calibration) {
 780     AliFatal("Could not get calibration object");
 781     return kFALSE;
 782   }
 783
 784   AliDebug(1,"Start creating digits");
 785
 786   // Create a container for the amplitudes
 787   AliTRDsegmentArray *signalsArray = new AliTRDsegmentArray("AliTRDdataArrayF"
 788                                                            ,AliTRDgeometry::Ndet());
 789
 790   AliLoader *gimme = fRunLoader->GetLoader("TRDLoader");
 791   if (!gimme->TreeH()) {
 792     gimme->LoadHits();
 793   }
 794   TTree *hitTree = gimme->TreeH();
 795   if (hitTree == 0x0) {
 796     AliError("Can not get TreeH");
 797     return kFALSE;
 798   }
 799   fTRD->SetTreeAddress();
 800
 801   // Get the number of entries in the hit tree
 802   // (Number of primary particles creating a hit somewhere)
 803   Int_t nTrack = (Int_t) hitTree->GetEntries();
 804   AliDebug(1,Form("Found %d primary particles",nTrack));
 805   AliDebug(1,Form("Sampling = %.0fMHz"        ,commonParam->GetSamplingFrequency()));
 806   AliDebug(1,Form("Gain     = %d"             ,((Int_t) simParam->GetGasGain())));
 807   AliDebug(1,Form("Noise    = %d"             ,((Int_t) simParam->GetNoise())));
 808   if (simParam->TimeStructOn()) {
 809     AliDebug(1,"Time Structure of drift cells implemented.");
 810   }
 811   else {
 812     AliDebug(1,"Constant drift velocity in drift cells.");
 813   }
 814   if (simParam->TRFOn()) {
 815     timeBinTRFend = ((Int_t) (simParam->GetTRFhi()
 816                   * commonParam->GetSamplingFrequency())) - 1;
 817     AliDebug(1,Form("Sample the TRF up to bin %d",timeBinTRFend));
 818   }
 819
 820   // Get the detector wise calibration objects
 821   const AliTRDCalDet *calVdriftDet     = calibration->GetVdriftDet();
 822   const AliTRDCalDet *calT0Det         = calibration->GetT0Det();
 823   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 824
 825   Int_t   detectorOld  = -1;
 826   Int_t   countHits    =  0;
 827
 828   Int_t   nTimeTotal   = calibration->GetNumberOfTimeBins();
 829   Float_t samplingRate = commonParam->GetSamplingFrequency();
 830   Float_t elAttachProp = simParam->GetElAttachProp() / 100.0;
 831
 832   // Loop through all entries in the tree
 833   for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
 834
 835     gAlice->ResetHits();
 836     nBytes += hitTree->GetEvent(iTrack);
 837
 838     // Loop through the TRD hits
 839     Int_t iHit = 0;
 840     AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
 841     while (hit) {
 842
 843       countHits++;
 844       iHit++;
 845
 846       Float_t q        = hit->GetCharge();
 847       // Don't analyze test hits
 848       if (((Int_t) q) == 0) {
 849         hit = (AliTRDhit *) fTRD->NextHit();
 850         continue;
 851       }
 852
 853       pos[0]           = hit->X();
 854       pos[1]           = hit->Y();
 855       pos[2]           = hit->Z();
 856       Int_t   track    = hit->Track();
 857       Int_t   detector = hit->GetDetector();
 858       Float_t hittime  = hit->GetTime();
 859       Int_t   plane    = fGeo->GetPlane(detector);
 860       Int_t   chamber  = fGeo->GetChamber(detector);
 861       padPlane         = commonParam->GetPadPlane(plane,chamber);
 862       Float_t row0     = padPlane->GetRow0ROC();
 863       Int_t   nRowMax  = padPlane->GetNrows();
 864       Int_t   nColMax  = padPlane->GetNcols();
 865       Int_t   inDrift  = 1;
 866
 867       // Find the current volume with the geo manager
 868       gGeoManager->SetCurrentPoint(pos);
 869       gGeoManager->FindNode();
 870       if (strstr(gGeoManager->GetPath(),"/UK")) {
 871         inDrift = 0;
 872       }
 873
 874       if (detector != detectorOld) {
 875
 876         // Compress the old one if enabled
 877         if ((fCompress) && (detectorOld > -1)) {
 878           signals->Compress(1,0);
 879           for (iDict = 0; iDict < kNDict; iDict++) {
 880             dictionary[iDict]->Compress(1,0);
 881           }
 882         }
 883         // Get the new container
 884         signals = (AliTRDdataArrayF *) signalsArray->At(detector);
 885         if (signals->GetNtime() == 0) {
 886           // Allocate a new one if not yet existing
 887           signals->Allocate(nRowMax,nColMax,nTimeTotal);
 888         }
 889         else {
 890           // Expand an existing one
 891           if (fCompress) {
 892             signals->Expand();
 893           }
 894         }
 895         // The same for the dictionary
 896         for (iDict = 0; iDict < kNDict; iDict++) {
 897           dictionary[iDict] = fDigitsManager->GetDictionary(detector,iDict);
 898           if (dictionary[iDict]->GetNtime() == 0) {
 899             dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
 900           }
 901           else {
 902             if (fCompress) dictionary[iDict]->Expand();
 903           }
 904         }
 905
 906         // Get the calibration objects
 907         calVdriftROC      = calibration->GetVdriftROC(detector);
 908         calVdriftDetValue = calVdriftDet->GetValue(detector);
 909         calT0ROC          = calibration->GetT0ROC(detector);
 910         calT0DetValue     = calT0Det->GetValue(detector);
 911
 912         detectorOld = detector;
 913
 914       }
 915
 916       // Go to the local coordinate system:
 917       // loc[0] - col  direction in amplification or driftvolume
 918       // loc[1] - row  direction in amplification or driftvolume
 919       // loc[2] - time direction in amplification or driftvolume
 920       gGeoManager->MasterToLocal(pos,loc);
 921       if (inDrift) {
 922         // Relative to middle of amplification region
 923         loc[2] = loc[2] - kDrWidth/2.0 - kAmWidth/2.0;
 924       }
 925
 926       // The driftlength [cm] (w/o diffusion yet !).
 927       // It is negative if the hit is between pad plane and anode wires.
 928       Double_t driftlength = -1.0 * loc[2];
 929
 930       // Stupid patch to take care of TR photons that are absorbed
 931       // outside the chamber volume. A real fix would actually need
 932       // a more clever implementation of the TR hit generation
 933       if (q < 0.0) {
 934         if ((loc[1] < padPlane->GetRowEndROC()) ||
 935             (loc[1] > padPlane->GetRow0ROC())) {
 936           AliDebug(2,Form("Hit outside of sensitive volume, row (z=%f, row0=%f, rowE=%f)\n"
 937                          ,loc[1],padPlane->GetRow0ROC(),padPlane->GetRowEndROC()));
 938           hit = (AliTRDhit *) fTRD->NextHit();
 939           continue;
 940         }
 941         if ((driftlength < kDrMin) ||
 942             (driftlength > kDrMax)) {
 943           AliDebug(2,Form("Hit outside of sensitive volume, time (Q = %d)\n"
 944                          ,((Int_t) q)));
 945           hit = (AliTRDhit *) fTRD->NextHit();
 946           continue;
 947         }
 948       }
 949
 950       // Get row and col of unsmeared electron to retrieve drift velocity
 951       // The pad row (z-direction)
 952       Int_t    rowE         = padPlane->GetPadRowNumberROC(loc[1]);
 953       if (rowE < 0) {
 954         hit = (AliTRDhit *) fTRD->NextHit();
 955         continue;
 956       }
 957       Double_t rowOffset    = padPlane->GetPadRowOffsetROC(rowE,loc[1]);
 958
 959       // The pad column (rphi-direction)
 960       Double_t offsetTilt   = padPlane->GetTiltOffset(rowOffset);
 961       Int_t    colE         = padPlane->GetPadColNumber(loc[0]+offsetTilt);
 962       if (colE < 0) {
 963         hit = (AliTRDhit *) fTRD->NextHit();
 964         continue;
 965       }
 966       Double_t colOffset    = padPlane->GetPadColOffset(colE,loc[0]+offsetTilt);
 967
 968       Float_t driftvelocity = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
 969
 970       // Normalized drift length
 971       Double_t absdriftlength = TMath::Abs(driftlength);
 972       if (commonParam->ExBOn()) {
 973         absdriftlength /= TMath::Sqrt(GetLorentzFactor(driftvelocity));
 974       }
 975
 976       // Loop over all electrons of this hit
 977       // TR photons produce hits with negative charge
 978       Int_t nEl = ((Int_t) TMath::Abs(q));
 979       for (Int_t iEl = 0; iEl < nEl; iEl++) {
 980
 981         // Now the real local coordinate system of the ROC
 982         // column direction: locC
 983         // row direction:    locR
 984         // time direction:   locT
 985         // locR and locC are identical to the coordinates of the corresponding
 986         // volumina of the drift or amplification region.
 987         // locT is defined relative to the wire plane (i.e. middle of amplification
 988         // region), meaming locT = 0, and is negative for hits coming from the
 989         // drift region.
 990         Double_t locC = loc[0];
 991         Double_t locR = loc[1];
 992         Double_t locT = loc[2];
 993
 994         // Electron attachment
 995         if (simParam->ElAttachOn()) {
 996           if (gRandom->Rndm() < (absdriftlength * elAttachProp)) {
 997             continue;
 998           }
 999         }
1000
1001         // Apply the diffusion smearing
1002         if (simParam->DiffusionOn()) {
1003           if (!(Diffusion(driftvelocity,absdriftlength,locR,locC,locT))) {
1004             continue;
1005           }
1006         }
1007
1008         // Apply E x B effects (depends on drift direction)
1009         if (commonParam->ExBOn()) {
1010           if (!(ExB(driftvelocity,driftlength,locC))) {
1011             continue;
1012           }
1013         }
1014
1015         // The electron position after diffusion and ExB in pad coordinates.
1016         // The pad row (z-direction)
1017         rowE       = padPlane->GetPadRowNumberROC(locR);
1018         if (rowE < 0) continue;
1019         rowOffset  = padPlane->GetPadRowOffsetROC(rowE,locR);
1020
1021         // The pad column (rphi-direction)
1022         offsetTilt = padPlane->GetTiltOffset(rowOffset);
1023         colE       = padPlane->GetPadColNumber(locC+offsetTilt);
1024         if (colE < 0) continue;
1025         colOffset  = padPlane->GetPadColOffset(colE,locC+offsetTilt);
1026
1027         // Also re-retrieve drift velocity because col and row may have changed
1028         driftvelocity = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
1029         Float_t t0    = calT0DetValue     + calT0ROC->GetValue(colE,rowE);
1030
1031         // Convert the position to drift time [mus], using either constant drift velocity or
1032         // time structure of drift cells (non-isochronity, GARFIELD calculation).
1033         // Also add absolute time of hits to take pile-up events into account properly
1034         Double_t drifttime;
1035         if (simParam->TimeStructOn()) {
1036           // Get z-position with respect to anode wire
1037           Double_t zz  =  row0 - locR + simParam->GetAnodeWireOffset();
1038           zz -= ((Int_t)(2 * zz)) / 2.0;
1039           if (zz > 0.25) {
1040             zz  = 0.5 - zz;
1041           }
1042           // Use drift time map (GARFIELD)
1043           drifttime = TimeStruct(driftvelocity,0.5*kAmWidth-1.0*locT,zz)
1044                     + hittime;
1045         }
1046         else {
1047           // Use constant drift velocity
1048           drifttime = -1.0 * locT / driftvelocity
1049                     + hittime;
1050         }
1051
1052         // Apply the gas gain including fluctuations
1053         Double_t ggRndm = 0.0;
1054         do {
1055           ggRndm = gRandom->Rndm();
1056         } while (ggRndm <= 0);
1057         Int_t signal = (Int_t) (-(simParam->GetGasGain()) * TMath::Log(ggRndm));
1058
1059         // Apply the pad response
1060         if (simParam->PRFOn()) {
1061           // The distance of the electron to the center of the pad
1062           // in units of pad width
1063           Double_t dist = (colOffset - 0.5*padPlane->GetColSize(colE))
1064                         / padPlane->GetColSize(colE);
1065           // This is still the fixed parametrization, i.e. not dependent on
1066           // calibration values !!!!
1067           if (!(calibration->PadResponse(signal,dist,plane,padSignal))) continue;
1068         }
1069         else {
1070           padSignal[0] = 0.0;
1071           padSignal[1] = signal;
1072           padSignal[2] = 0.0;
1073         }
1074
1075         // The time bin (always positive), with t0 distortion
1076         Double_t timeBinIdeal = drifttime * samplingRate + t0;
1077         // Protection
1078         if (TMath::Abs(timeBinIdeal) > 2*nTimeTotal) {
1079           timeBinIdeal = 2 * nTimeTotal;
1080         }
1081         Int_t    timeBinTruncated = (Int_t) timeBinIdeal;
1082         // The distance of the position to the middle of the timebin
1083         Double_t timeOffset       = ((Float_t) timeBinTruncated
1084                                   + 0.5 - timeBinIdeal) / samplingRate;
1085
1086         // Sample the time response inside the drift region
1087         // + additional time bins before and after.
1088         // The sampling is done always in the middle of the time bin
1089         for (Int_t iTimeBin = TMath::Max(timeBinTruncated,0)
1090             ;iTimeBin < TMath::Min(timeBinTruncated+timeBinTRFend,nTimeTotal)
1091             ;iTimeBin++) {
1092
1093           // Apply the time response
1094           Double_t timeResponse = 1.0;
1095           Double_t crossTalk    = 0.0;
1096           Double_t time         = (iTimeBin - timeBinTruncated) / samplingRate + timeOffset;
1097           if (simParam->TRFOn()) {
1098             timeResponse = simParam->TimeResponse(time);
1099           }
1100           if (simParam->CTOn()) {
1101             crossTalk    = simParam->CrossTalk(time);
1102           }
1103
1104           signalOld[0] = 0.0;
1105           signalOld[1] = 0.0;
1106           signalOld[2] = 0.0;
1107
1108           for (iPad = 0; iPad < kNpad; iPad++) {
1109
1110             Int_t colPos = colE + iPad - 1;
1111             if (colPos <        0) continue;
1112             if (colPos >= nColMax) break;
1113
1114             // Add the signals
1115             Int_t iCurrentTimeBin = iTimeBin;
1116             signalOld[iPad]  = signals->GetDataUnchecked(rowE,colPos,iCurrentTimeBin);
1117             if (colPos != colE) {
1118               signalOld[iPad] += padSignal[iPad] * (timeResponse + crossTalk);
1119             }
1120             else {
1121               signalOld[iPad] += padSignal[iPad] * timeResponse;
1122             }
1123             signals->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,signalOld[iPad]);
1124
1125             // Store the track index in the dictionary
1126             // Note: We store index+1 in order to allow the array to be compressed
1127             if (signalOld[iPad] > 0) {
1128               for (iDict = 0; iDict < kNDict; iDict++) {
1129                 Int_t oldTrack = dictionary[iDict]->GetDataUnchecked(rowE
1130                                                                     ,colPos
1131                                                                     ,iCurrentTimeBin);
1132                 if (oldTrack == track+1) break;
1133                 if (oldTrack ==       0) {
1134                   dictionary[iDict]->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,track+1);
1135                   break;
1136                 }
1137               }
1138             }
1139
1140           } // Loop: pads
1141
1142         } // Loop: time bins
1143
1144       } // Loop: electrons of a single hit
1145
1146       hit = (AliTRDhit *) fTRD->NextHit();
1147
1148     } // Loop: hits of one primary track
1149
1150   } // Loop: primary tracks
1151
1152   AliDebug(1,Form("Finished analyzing %d hits",countHits));
1153
1154   //____________________________________________________________________
1155   //
1156   // Create (s)digits
1157   //____________________________________________________________________
1158   //
1159
1160   // The coupling factor
1161   Double_t coupling = simParam->GetPadCoupling()
1162                     * simParam->GetTimeCoupling();
1163
1164   // The conversion factor
1165   Double_t convert  = kEl2fC
1166                     * simParam->GetChipGain();
1167
1168   // Loop through all chambers to finalize the digits
1169   Int_t iDetBeg = 0;
1170   Int_t iDetEnd = AliTRDgeometry::Ndet();
1171   for (Int_t iDet = iDetBeg; iDet < iDetEnd; iDet++) {
1172
1173     Int_t plane      = fGeo->GetPlane(iDet);
1174     Int_t sector     = fGeo->GetSector(iDet);
1175     Int_t chamber    = fGeo->GetChamber(iDet);
1176     Int_t nRowMax    = commonParam->GetRowMax(plane,chamber,sector);
1177     Int_t nColMax    = commonParam->GetColMax(plane);
1178
1179     Double_t *inADC  = new Double_t[nTimeTotal];
1180     Double_t *outADC = new Double_t[nTimeTotal];
1181
1182     // Add a container for the digits of this detector
1183     digits = fDigitsManager->GetDigits(iDet);
1184     // Allocate memory space for the digits buffer
1185     if (digits->GetNtime() == 0) {
1186       digits->Allocate(nRowMax,nColMax,nTimeTotal);
1187     }
1188
1189     // Get the signal container
1190     signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
1191     if (signals->GetNtime() == 0) {
1192       // Create missing containers
1193       signals->Allocate(nRowMax,nColMax,nTimeTotal);
1194     }
1195     else {
1196       // Expand the container if neccessary
1197       if (fCompress) signals->Expand();
1198     }
1199     // Create the missing dictionary containers
1200     for (iDict = 0; iDict < kNDict; iDict++) {
1201       dictionary[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1202       if (dictionary[iDict]->GetNtime() == 0) {
1203         dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1204       }
1205     }
1206
1207     Int_t nDigits = 0;
1208
1209     // Don't create noise in detectors that are switched off / not installed, etc.
1210     if (( calibration->IsChamberInstalled(iDet)) &&
1211         (!calibration->IsChamberMasked(iDet))    &&
1212         ( fGeo->GetSMstatus(sector))) {
1213
1214       // Get the calibration objects
1215       calGainFactorROC      = calibration->GetGainFactorROC(iDet);
1216       calGainFactorDetValue = calGainFactorDet->GetValue(iDet);
1217
1218       // Create the digits for this chamber
1219       for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1220         for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1221
1222           // Check whether pad is masked
1223           // Bridged pads are not considered yet!!!
1224           if (calibration->IsPadMasked(iDet,iCol,iRow)) {
1225             continue;
1226           }
1227
1228           // Create summable digits
1229           if (fSDigits) {
1230
1231             for (iTime = 0; iTime < nTimeTotal; iTime++) {
1232               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1233               signalAmp *= fSDigitsScale;
1234               signalAmp  = TMath::Min(signalAmp,(Float_t) 1.0e9);
1235               Int_t adc  = (Int_t) signalAmp;
1236               if (adc > 0) nDigits++;
1237               digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1238             }
1239
1240           }
1241           // Create normal digits
1242           else {
1243
1244             Float_t padgain = calGainFactorDetValue
1245                             * calGainFactorROC->GetValue(iCol,iRow);
1246             if (padgain <= 0) {
1247               AliError(Form("Not a valid gain %f, %d %d %d",padgain,iDet,iCol,iRow));
1248             }
1249
1250             for (iTime = 0; iTime < nTimeTotal; iTime++) {
1251
1252               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1253
1254               // Pad and time coupling
1255               signalAmp *= coupling;
1256               // Gain factors
1257               signalAmp *= padgain;
1258
1259               // Add the noise, starting from minus ADC baseline in electrons
1260               Double_t baselineEl = simParam->GetADCbaseline() * (simParam->GetADCinRange()
1261                                                                / simParam->GetADCoutRange())
1262                                                                / convert;
1263               signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,simParam->GetNoise())
1264                                      ,-baselineEl);
1265               // Convert to mV
1266               signalAmp *= convert;
1267               // Add ADC baseline in mV
1268               signalAmp += simParam->GetADCbaseline() * (simParam->GetADCinRange()
1269                                                       / simParam->GetADCoutRange());
1270               // Convert to ADC counts. Set the overflow-bit fADCoutRange if the
1271               // signal is larger than fADCinRange
1272               Int_t adc  = 0;
1273               if (signalAmp >= simParam->GetADCinRange()) {
1274                 adc = ((Int_t) simParam->GetADCoutRange());
1275               }
1276               else {
1277                 adc = TMath::Nint(signalAmp * (simParam->GetADCoutRange()
1278                                             / simParam->GetADCinRange()));
1279               }
1280
1281               inADC[iTime]  = adc;
1282               outADC[iTime] = adc;
1283
1284             }
1285
1286             for (iTime = 0; iTime < nTimeTotal; iTime++) {
1287               // Store the amplitude of the digit if above threshold
1288               if ( outADC[iTime] > ( simParam->GetADCbaseline() + simParam->GetADCthreshold() ) ) {
1289                 nDigits++;
1290                 digits->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
1291               }
1292             }
1293
1294           }
1295
1296         }
1297       }
1298
1299     }
1300
1301     // Compress the arrays
1302     digits->Compress(1,0);
1303     for (iDict = 0; iDict < kNDict; iDict++) {
1304       dictionary[iDict]->Compress(1,0);
1305     }
1306
1307     totalSizeDigits += digits->GetSize();
1308     totalSizeDict0  += dictionary[0]->GetSize();
1309     totalSizeDict1  += dictionary[1]->GetSize();
1310     totalSizeDict2  += dictionary[2]->GetSize();
1311
1312     if (nDigits > 0) {
1313       Float_t nPixel = nRowMax * nColMax * nTimeTotal;
1314       AliDebug(1,Form("Found %d digits in detector %d (%3.0f)."
1315                      ,nDigits,iDet
1316                      ,100.0 * ((Float_t) nDigits) / nPixel));
1317     }
1318
1319     if (fCompress) {
1320       signals->Compress(1,0);
1321     }
1322
1323     delete [] inADC;
1324     delete [] outADC;
1325
1326   }
1327
1328   if (signalsArray) {
1329     delete signalsArray;
1330     signalsArray = 0;
1331   }
1332
1333   AliDebug(1,Form("Total number of analyzed hits = %d",countHits));
1334   AliDebug(1,Form("Total digits data size = %d, %d, %d, %d",totalSizeDigits
1335                                                            ,totalSizeDict0
1336                                                            ,totalSizeDict1
1337                                                            ,totalSizeDict2));
1338
1339   return kTRUE;
1340
1341 }
1342
1343 //_____________________________________________________________________________
1344 void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
1345 {
1346   //
1347   // Add a digits manager for s-digits to the input list.
1348   //
1349
1350   fSDigitsManagerList->Add(man);
1351
1352 }
1353
1354 //_____________________________________________________________________________
1355 void AliTRDdigitizer::DeleteSDigitsManager()
1356 {
1357   //
1358   // Removes digits manager from the input list.
1359   //
1360
1361   fSDigitsManagerList->Delete();
1362
1363 }
1364
1365 //_____________________________________________________________________________
1366 Bool_t AliTRDdigitizer::ConvertSDigits()
1367 {
1368   //
1369   // Converts s-digits to normal digits
1370   //
1371
1372   // Number of track dictionary arrays
1373   const Int_t    kNDict = AliTRDdigitsManager::kNDict;
1374
1375   // Converts number of electrons to fC
1376   const Double_t kEl2fC = 1.602e-19 * 1.0e15;
1377
1378   Int_t iDict = 0;
1379   Int_t iRow;
1380   Int_t iCol;
1381   Int_t iTime;
1382
1383   AliTRDCalROC *calGainFactorROC      = 0;
1384   Float_t       calGainFactorDetValue = 0.0;
1385
1386   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
1387   if (!simParam) {
1388     AliFatal("Could not get simulation parameters");
1389     return kFALSE;
1390   }
1391
1392   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
1393   if (!commonParam) {
1394     AliFatal("Could not get common parameters");
1395     return kFALSE;
1396   }
1397
1398   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1399   if (!calibration) {
1400     AliFatal("Could not get calibration object");
1401     return kFALSE;
1402   }
1403
1404   Double_t sDigitsScale = 1.0 / GetSDigitsScale();
1405   Double_t noise        = simParam->GetNoise();
1406   Double_t padCoupling  = simParam->GetPadCoupling();
1407   Double_t timeCoupling = simParam->GetTimeCoupling();
1408   Double_t chipGain     = simParam->GetChipGain();
1409   Double_t coupling     = padCoupling * timeCoupling;
1410   Double_t convert      = kEl2fC * chipGain;
1411   Double_t adcInRange   = simParam->GetADCinRange();
1412   Double_t adcOutRange  = simParam->GetADCoutRange();
1413   Int_t    adcThreshold = simParam->GetADCthreshold();
1414   Int_t    adcBaseline  = simParam->GetADCbaseline();
1415   Int_t    nTimeTotal   = calibration->GetNumberOfTimeBins();
1416
1417   AliTRDdataArrayI *digitsIn;
1418   AliTRDdataArrayI *digitsOut;
1419   AliTRDdataArrayI *dictionaryIn[kNDict];
1420   AliTRDdataArrayI *dictionaryOut[kNDict];
1421
1422   // Get the detector wise calibration objects
1423   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
1424
1425   // Loop through the detectors
1426   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1427
1428     Int_t plane      = fGeo->GetPlane(iDet);
1429     Int_t sector     = fGeo->GetSector(iDet);
1430     Int_t chamber    = fGeo->GetChamber(iDet);
1431     Int_t nRowMax    = commonParam->GetRowMax(plane,chamber,sector);
1432     Int_t nColMax    = commonParam->GetColMax(plane);
1433
1434     Double_t *inADC  = new Double_t[nTimeTotal];
1435     Double_t *outADC = new Double_t[nTimeTotal];
1436
1437     digitsIn  = fSDigitsManager->GetDigits(iDet);
1438     digitsIn->Expand();
1439     digitsOut = fDigitsManager->GetDigits(iDet);
1440     digitsOut->Allocate(nRowMax,nColMax,nTimeTotal);
1441     for (iDict = 0; iDict < kNDict; iDict++) {
1442       dictionaryIn[iDict]  = fSDigitsManager->GetDictionary(iDet,iDict);
1443       dictionaryIn[iDict]->Expand();
1444       dictionaryOut[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1445       dictionaryOut[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1446     }
1447
1448     // Don't create noise in detectors that are switched off / not installed, etc.
1449     if (( calibration->IsChamberInstalled(iDet)) &&
1450         (!calibration->IsChamberMasked(iDet))    &&
1451         ( fGeo->GetSMstatus(sector))) {
1452
1453       // Get the calibration objects
1454       calGainFactorROC      = calibration->GetGainFactorROC(iDet);
1455       calGainFactorDetValue = calGainFactorDet->GetValue(iDet);
1456
1457       for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1458         for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1459
1460           // Check whether pad is masked
1461           // Bridged pads are not considered yet!!!
1462           if (calibration->IsPadMasked(iDet,iCol,iRow)) {
1463             continue;
1464           }
1465
1466           Float_t padgain = calGainFactorDetValue
1467                           * calGainFactorROC->GetValue(iCol,iRow);
1468           if (padgain <= 0) {
1469             AliError(Form("Not a valid gain %f, %d %d %d",padgain,iDet,iCol,iRow));
1470           }
1471
1472           for (iTime = 0; iTime < nTimeTotal; iTime++) {
1473
1474             // Scale s-digits to normal digits
1475             Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1476             signal         *= sDigitsScale;
1477             // Pad and time coupling
1478             signal *= coupling;
1479             // Gain factors
1480             signal *= padgain;
1481             // Add the noise, starting from minus ADC baseline in electrons
1482             Double_t baselineEl = adcBaseline * (adcInRange / adcOutRange) / convert;
1483             signal  = TMath::Max((Double_t) gRandom->Gaus(signal,noise),-baselineEl);
1484             // Convert to mV
1485             signal *= convert;
1486             // add ADC baseline in mV
1487             signal += adcBaseline * (adcInRange / adcOutRange);
1488             // Convert to ADC counts. Set the overflow-bit adcOutRange if the
1489             // signal is larger than adcInRange
1490             Int_t adc  = 0;
1491             if (signal >= adcInRange) {
1492               adc = ((Int_t) adcOutRange);
1493             }
1494             else {
1495               adc = TMath::Nint(signal * (adcOutRange / adcInRange));
1496             }
1497             inADC[iTime]  = adc;
1498             outADC[iTime] = adc;
1499
1500           }
1501
1502           for (iTime = 0; iTime < nTimeTotal; iTime++) {
1503             // Store the amplitude of the digit if above threshold
1504             if (outADC[iTime] > (adcBaseline + adcThreshold)) {
1505               digitsOut->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
1506               // Copy the dictionary
1507               for (iDict = 0; iDict < kNDict; iDict++) {
1508                 Int_t track = dictionaryIn[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1509                 dictionaryOut[iDict]->SetDataUnchecked(iRow,iCol,iTime,track);
1510               }
1511             }
1512           }
1513
1514         }
1515       }
1516
1517     }
1518
1519     if (fCompress) {
1520       digitsIn->Compress(1,0);
1521       digitsOut->Compress(1,0);
1522       for (iDict = 0; iDict < kNDict; iDict++) {
1523         dictionaryIn[iDict]->Compress(1,0);
1524         dictionaryOut[iDict]->Compress(1,0);
1525       }
1526     }
1527
1528     delete [] inADC;
1529     delete [] outADC;
1530
1531   }
1532
1533   return kTRUE;
1534
1535 }
1536
1537 //_____________________________________________________________________________
1538 Bool_t AliTRDdigitizer::MergeSDigits()
1539 {
1540   //
1541   // Merges the input s-digits:
1542   //   - The amplitude of the different inputs are summed up.
1543   //   - Of the track IDs from the input dictionaries only one is
1544   //     kept for each input. This works for maximal 3 different merged inputs.
1545   //
1546
1547   // Number of track dictionary arrays
1548   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1549
1550   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
1551   if (!simParam) {
1552     AliFatal("Could not get simulation parameters");
1553     return kFALSE;
1554   }
1555
1556   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
1557   if (!commonParam) {
1558     AliFatal("Could not get common parameters");
1559     return kFALSE;
1560   }
1561
1562   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1563   if (!calibration) {
1564     AliFatal("Could not get calibration object");
1565     return kFALSE;
1566   }
1567
1568   Int_t iDict = 0;
1569   Int_t jDict = 0;
1570
1571   AliTRDdataArrayI *digitsA;
1572   AliTRDdataArrayI *digitsB;
1573   AliTRDdataArrayI *dictionaryA[kNDict];
1574   AliTRDdataArrayI *dictionaryB[kNDict];
1575
1576   // Get the first s-digits
1577   fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
1578   if (!fSDigitsManager) {
1579     AliError("No SDigits manager");
1580     return kFALSE;
1581   }
1582
1583   // Loop through the other sets of s-digits
1584   AliTRDdigitsManager *mergeSDigitsManager;
1585   mergeSDigitsManager = (AliTRDdigitsManager *)
1586                         fSDigitsManagerList->After(fSDigitsManager);
1587
1588   if (mergeSDigitsManager) {
1589     AliDebug(1,Form("Merge %d input files.",fSDigitsManagerList->GetSize()));
1590   }
1591   else {
1592     AliDebug(1,"Only one input file.");
1593   }
1594
1595   Int_t nTimeTotal = calibration->GetNumberOfTimeBins();
1596
1597   Int_t iMerge = 0;
1598   while (mergeSDigitsManager) {
1599
1600     iMerge++;
1601
1602     // Loop through the detectors
1603     for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1604
1605       Int_t plane   = fGeo->GetPlane(iDet);
1606       Int_t sector  = fGeo->GetSector(iDet);
1607       Int_t chamber = fGeo->GetChamber(iDet);
1608       Int_t nRowMax = commonParam->GetRowMax(plane,chamber,sector);
1609       Int_t nColMax = commonParam->GetColMax(plane);
1610
1611       // Loop through the pixels of one detector and add the signals
1612       digitsA = fSDigitsManager->GetDigits(iDet);
1613       digitsB = mergeSDigitsManager->GetDigits(iDet);
1614       digitsA->Expand();
1615       digitsB->Expand();
1616       for (iDict = 0; iDict < kNDict; iDict++) {
1617         dictionaryA[iDict] = fSDigitsManager->GetDictionary(iDet,iDict);
1618         dictionaryB[iDict] = mergeSDigitsManager->GetDictionary(iDet,iDict);
1619         dictionaryA[iDict]->Expand();
1620         dictionaryB[iDict]->Expand();
1621       }
1622
1623       // Merge only detectors that contain a signal
1624       Bool_t doMerge = kTRUE;
1625       if (fMergeSignalOnly) {
1626         if (digitsA->GetOverThreshold(0) == 0) {
1627           doMerge = kFALSE;
1628         }
1629       }
1630
1631       if (doMerge) {
1632
1633         AliDebug(1,Form("Merge detector %d of input no.%d",iDet,iMerge+1));
1634
1635         for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1636           for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1637             for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
1638
1639               // Add the amplitudes of the summable digits
1640               Int_t ampA = digitsA->GetDataUnchecked(iRow,iCol,iTime);
1641               Int_t ampB = digitsB->GetDataUnchecked(iRow,iCol,iTime);
1642               ampA += ampB;
1643               digitsA->SetDataUnchecked(iRow,iCol,iTime,ampA);
1644
1645               // Add the mask to the track id if defined.
1646               for (iDict = 0; iDict < kNDict; iDict++) {
1647                 Int_t trackB = dictionaryB[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1648                 if ((fMasks) && (trackB > 0)) {
1649                   for (jDict = 0; jDict < kNDict; jDict++) {
1650                     Int_t trackA = dictionaryA[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1651                     if (trackA == 0) {
1652                       trackA = trackB + fMasks[iMerge];
1653                       dictionaryA[iDict]->SetDataUnchecked(iRow,iCol,iTime,trackA);
1654                     }
1655                   }
1656                 }
1657               }
1658
1659             }
1660           }
1661         }
1662
1663       }
1664
1665       if (fCompress) {
1666         digitsA->Compress(1,0);
1667         digitsB->Compress(1,0);
1668         for (iDict = 0; iDict < kNDict; iDict++) {
1669           dictionaryA[iDict]->Compress(1,0);
1670           dictionaryB[iDict]->Compress(1,0);
1671         }
1672       }
1673
1674     }
1675
1676     // The next set of s-digits
1677     mergeSDigitsManager = (AliTRDdigitsManager *)
1678                           fSDigitsManagerList->After(mergeSDigitsManager);
1679
1680   }
1681
1682   return kTRUE;
1683
1684 }
1685
1686 //_____________________________________________________________________________
1687 Bool_t AliTRDdigitizer::SDigits2Digits()
1688 {
1689   //
1690   // Merges the input s-digits and converts them to normal digits
1691   //
1692
1693   if (!MergeSDigits()) {
1694     return kFALSE;
1695   }
1696
1697   return ConvertSDigits();
1698
1699 }
1700
1701 //_____________________________________________________________________________
1702 Bool_t AliTRDdigitizer::WriteDigits() const
1703 {
1704   //
1705   // Writes out the TRD-digits and the dictionaries
1706   //
1707
1708   // Write parameters
1709   fRunLoader->CdGAFile();
1710
1711   // Store the digits and the dictionary in the tree
1712   return fDigitsManager->WriteDigits();
1713
1714 }
1715
1716 //_____________________________________________________________________________
1717 void AliTRDdigitizer::InitOutput(Int_t iEvent)
1718 {
1719   //
1720   // Initializes the output branches
1721   //
1722
1723   fEvent = iEvent;
1724
1725   if (!fRunLoader) {
1726     AliError("Run Loader is NULL");
1727     return;
1728   }
1729
1730   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
1731   if (!loader) {
1732     AliError("Can not get TRD loader from Run Loader");
1733     return;
1734   }
1735
1736   TTree *tree = 0;
1737
1738   if (fSDigits) {
1739     // If we produce SDigits
1740     tree = loader->TreeS();
1741     if (!tree) {
1742       loader->MakeTree("S");
1743       tree = loader->TreeS();
1744     }
1745   }
1746   else {
1747     // If we produce Digits
1748     tree = loader->TreeD();
1749     if (!tree) {
1750       loader->MakeTree("D");
1751       tree = loader->TreeD();
1752     }
1753   }
1754
1755   fDigitsManager->SetEvent(iEvent);
1756   fDigitsManager->MakeBranch(tree);
1757
1758 }
1759
1760 //_____________________________________________________________________________
1761 Double_t AliTRDdigitizer::TimeStruct(Float_t vdrift, Double_t dist, Double_t z)
1762 {
1763   //
1764   // Applies the time structure of the drift cells (by C.Lippmann).
1765   // The drift time of electrons to the anode wires depends on the
1766   // distance to the wire (z) and on the position in the drift region.
1767   //
1768   // input :
1769   // dist = radial distance from (cathode) pad plane [cm]
1770   // z    = distance from anode wire (parallel to cathode planes) [cm]
1771   //
1772   // output :
1773   // tdrift = the drift time of an electron at the given position
1774   //
1775   // We interpolate between the drift time values at the two drift
1776   // velocities fVDlo and fVDhi, being smaller and larger than
1777   // fDriftVelocity. We use the two stored drift time maps fTimeStruct1
1778   // and fTimeStruct2, calculated for the two mentioned drift velocities.
1779   //
1780
1781   SampleTimeStruct(vdrift);
1782
1783   // Indices:
1784   Int_t r1 = (Int_t)(10 * dist);
1785   if (r1 <  0) r1 =  0;
1786   if (r1 > 37) r1 = 37;
1787   Int_t r2 = r1 + 1;
1788   if (r2 <  0) r2 =  0;
1789   if (r2 > 37) r2 = 37;
1790   const Int_t kz1 = ((Int_t)(100 * z / 2.5));
1791   const Int_t kz2 = kz1 + 1;
1792
1793   if ((r1  <  0) ||
1794       (r1  > 37) ||
1795       (kz1 <  0) ||
1796       (kz1 > 10)) {
1797     AliWarning(Form("Indices out of range: dist=%.2f, z=%.2f, r1=%d, kz1=%d"
1798                    ,dist,z,r1,kz1));
1799   }
1800
1801   const Float_t ky111 = fTimeStruct1[r1+38*kz1];
1802   const Float_t ky221 = ((r2 <= 37) && (kz2 <= 10))
1803                       ? fTimeStruct1[r2+38*kz2]
1804                       : fTimeStruct1[37+38*10];
1805   const Float_t ky121 = (kz2 <= 10)
1806                       ? fTimeStruct1[r1+38*kz2]
1807                       : fTimeStruct1[r1+38*10];
1808   const Float_t ky211 = (r2 <= 37)
1809                       ? fTimeStruct1[r2+38*kz1]
1810                       : fTimeStruct1[37+38*kz1];
1811
1812   // 2D Interpolation, lower drift time map
1813   const Float_t ky11  = (ky211-ky111)*10*dist + ky111 - (ky211-ky111)*r1;
1814   const Float_t ky21  = (ky221-ky121)*10*dist + ky121 - (ky221-ky121)*r1;
1815
1816   const Float_t ky112 = fTimeStruct2[r1+38*kz1];
1817   const Float_t ky222 = ((r2 <= 37) && (kz2 <= 10))
1818                       ? fTimeStruct2[r2+38*kz2]
1819                       : fTimeStruct2[37+38*10];
1820   const Float_t ky122 = (kz2 <= 10)
1821                       ? fTimeStruct2[r1+38*kz2]
1822                       : fTimeStruct2[r1+38*10];
1823   const Float_t ky212 = (r2 <= 37)
1824                       ? fTimeStruct2[r2+38*kz1]
1825                       : fTimeStruct2[37+38*kz1];
1826
1827   // 2D Interpolation, larger drift time map
1828   const Float_t ky12  = (ky212-ky112)*10*dist + ky112 - (ky212-ky112)*r1;
1829   const Float_t ky22  = (ky222-ky122)*10*dist + ky122 - (ky222-ky122)*r1;
1830
1831   // Dist now is the drift distance to the anode wires (negative if electrons are
1832   // between anode wire plane and cathode pad plane)
1833   dist -= AliTRDgeometry::AmThick() / 2.0;
1834
1835   // Get the drift times for the drift velocities fVDlo and fVDhi
1836   const Float_t ktdrift1 = ((TMath::Abs(dist) > 0.005) || (z > 0.005))
1837                          ? (ky21 - ky11) * 100 * z / 2.5 + ky11 - (ky21 - ky11) * kz1
1838                          : 0.0;
1839   const Float_t ktdrift2 = ((TMath::Abs(dist) > 0.005) || (z > 0.005))
1840                          ? (ky22 - ky12) * 100 * z / 2.5 + ky12 - (ky22 - ky12) * kz1
1841                          : 0.0;
1842
1843   // 1D Interpolation between the values at fVDlo and fVDhi
1844   Float_t a = (ktdrift2 - ktdrift1)
1845             / (fVDhi - fVDlo);
1846   Float_t b = ktdrift2 - a * fVDhi;
1847
1848   return a * vdrift + b;
1849
1850 }
1851
1852 //_____________________________________________________________________________
1853 void AliTRDdigitizer::SampleTimeStruct(Float_t vdrift)
1854 {
1855   //
1856   // Samples the timing structure of a drift cell
1857   // Drift Time data calculated with Garfield (by C.Lippmann)
1858   //
1859
1860   // Nothing to do
1861   if (vdrift == fTimeLastVdrift) {
1862     return;
1863   }
1864   fTimeLastVdrift = vdrift;
1865
1866   // Drift time maps are saved for some drift velocity values (in drift region):
1867   Float_t  fVDsmp[8];
1868   fVDsmp[0] = 1.032;
1869   fVDsmp[1] = 1.158;
1870   fVDsmp[2] = 1.299;
1871   fVDsmp[3] = 1.450;
1872   fVDsmp[4] = 1.610;
1873   fVDsmp[5] = 1.783;
1874   fVDsmp[6] = 1.959;
1875   fVDsmp[7] = 2.134;
1876
1877   if      (vdrift < fVDsmp[0]) {
1878     AliWarning(Form("Drift Velocity too small (%.3f<%.3f)",vdrift,fVDsmp[0]));
1879     vdrift = fVDsmp[0];
1880   }
1881   else if (vdrift > fVDsmp[7]) {
1882     AliWarning(Form("Drift Velocity too large (%.3f>%.3f)",vdrift,fVDsmp[6]));
1883     vdrift = fVDsmp[7];
1884   }
1885
1886   const Int_t ktimebin  = 38;
1887   const Int_t kZbin     = 11;
1888
1889   // Garfield simulation at UD = -1500V; vd = 1.032cm/microsec, <driftfield> = 525V/cm
1890   Float_t time1500[ktimebin][kZbin] =
1891       {{0.09170, 0.09205, 0.09306, 0.09475, 0.09716, 0.10035,
1892         0.10445, 0.10993, 0.11838, 0.13986, 0.37858},
1893        {0.06588, 0.06626, 0.06739, 0.06926, 0.07186, 0.07524,
1894         0.07951, 0.08515, 0.09381, 0.11601, 0.35673},
1895        {0.03946, 0.04003, 0.04171, 0.04435, 0.04780, 0.05193,
1896         0.05680, 0.06306, 0.07290, 0.09873, 0.34748},
1897        {0.01151, 0.01283, 0.01718, 0.02282, 0.02880, 0.03479,
1898         0.04098, 0.04910, 0.06413, 0.10567, 0.36897},
1899        {0.01116, 0.01290, 0.01721, 0.02299, 0.02863, 0.03447,
1900         0.04074, 0.04984, 0.06839, 0.11625, 0.37277},
1901        {0.03919, 0.03974, 0.04131, 0.04380, 0.04703, 0.05102,
1902         0.05602, 0.06309, 0.07651, 0.10938, 0.36838},
1903        {0.06493, 0.06560, 0.06640, 0.06802, 0.07051, 0.07392,
1904         0.07853, 0.08510, 0.09690, 0.12621, 0.38058},
1905        {0.09174, 0.09186, 0.09225, 0.09303, 0.09477, 0.00000,
1906         0.11205, 0.11952, 0.13461, 0.16984, 0.43017},
1907        {0.14356, 0.14494, 0.14959, 0.16002, 0.18328, 0.27981,
1908         0.22785, 0.21240, 0.21948, 0.25965, 0.52392},
1909        {0.23120, 0.23366, 0.24046, 0.25422, 0.28071, 0.36914,
1910         0.32999, 0.31208, 0.31772, 0.35804, 0.62249},
1911        {0.32686, 0.32916, 0.33646, 0.35053, 0.37710, 0.46292,
1912         0.42773, 0.40948, 0.41497, 0.45527, 0.71955},
1913        {0.42353, 0.42583, 0.43317, 0.44727, 0.47380, 0.55884,
1914         0.52479, 0.50650, 0.51194, 0.55225, 0.81658},
1915        {0.52038, 0.52271, 0.53000, 0.54415, 0.57064, 0.65545,
1916         0.62172, 0.60341, 0.60885, 0.64915, 0.91339},
1917        {0.61724, 0.61953, 0.62694, 0.64098, 0.66756, 0.75226,
1918         0.71862, 0.70030, 0.70575, 0.74604, 1.01035},
1919        {0.71460, 0.71678, 0.72376, 0.73786, 0.76447, 0.84913,
1920         0.81551, 0.79720, 0.80264, 0.84292, 1.10723},
1921        {0.81101, 0.81334, 0.82066, 0.83475, 0.86127, 0.94599,
1922         0.91240, 0.89408, 0.89952, 0.93981, 1.20409},
1923        {0.90788, 0.91023, 0.91752, 0.93163, 0.95815, 1.04293,
1924         1.00929, 0.99096, 0.99640, 1.03669, 1.30106},
1925        {1.00477, 1.00707, 1.01449, 1.02852, 1.05504, 1.13976,
1926         1.10617, 1.08784, 1.09329, 1.13358, 1.39796},
1927        {1.10166, 1.10397, 1.11130, 1.12541, 1.15257, 1.23672,
1928         1.20307, 1.18472, 1.19018, 1.23046, 1.49486},
1929        {1.19854, 1.20084, 1.20818, 1.22235, 1.24885, 1.33355,
1930         1.29992, 1.28156, 1.28707, 1.32735, 1.59177},
1931        {1.29544, 1.29780, 1.30507, 1.31917, 1.34575, 1.43073,
1932         1.39681, 1.37851, 1.38396, 1.42377, 1.68854},
1933        {1.39236, 1.39462, 1.40205, 1.41607, 1.44259, 1.52745,
1934         1.49368, 1.47541, 1.48083, 1.52112, 1.78546},
1935        {1.49314, 1.49149, 1.49885, 1.51297, 1.53949, 1.62420,
1936         1.59016, 1.57231, 1.57772, 1.61800, 1.88048},
1937        {1.58610, 1.58839, 1.59572, 1.60983, 1.63635, 1.72109,
1938         1.68651, 1.66921, 1.67463, 1.71489, 1.97918},
1939        {1.68400, 1.68529, 1.69261, 1.70671, 1.73331, 1.81830,
1940         1.78341, 1.76605, 1.77150, 1.81179, 2.07608},
1941        {1.77991, 1.78215, 1.78952, 1.80385, 1.83014, 1.91486,
1942         1.88128, 1.86215, 1.86837, 1.90865, 2.17304},
1943        {1.87674, 1.87904, 1.88647, 1.90052, 1.92712, 2.01173,
1944         1.97812, 1.95905, 1.96527, 2.00710, 2.26979},
1945        {1.97369, 1.97594, 1.98326, 1.99869, 2.02442, 2.10865,
1946         2.07501, 2.05666, 2.06214, 2.10243, 2.36669},
1947        {2.07052, 2.07281, 2.08016, 2.09425, 2.12132, 2.20555,
1948         2.17182, 2.15341, 2.15904, 2.19933, 2.46363},
1949        {2.16742, 2.16971, 2.17707, 2.19114, 2.21766, 2.30240,
1950         2.26877, 2.25015, 2.25573, 2.29586, 2.56060},
1951        {2.26423, 2.26659, 2.27396, 2.28803, 2.31456, 2.40828,
1952         2.36567, 2.34705, 2.35282, 2.39765, 2.65744},
1953        {2.36153, 2.36349, 2.37330, 2.38501, 2.41159, 2.49940,
1954         2.46257, 2.44420, 2.44843, 2.48987, 2.75431},
1955        {2.46558, 2.46035, 2.46822, 2.48181, 2.50849, 2.59630,
1956         2.55947, 2.54112, 2.54513, 2.58677, 2.85094},
1957        {2.56248, 2.55723, 2.56486, 2.57871, 2.60520, 2.68998,
1958         2.65626, 2.63790, 2.64316, 2.68360, 2.94813},
1959        {2.65178, 2.65441, 2.66153, 2.67556, 2.70210, 2.78687,
1960         2.75319, 2.73463, 2.74032, 2.78060, 3.04503},
1961        {2.74868, 2.75131, 2.75870, 2.77245, 2.79385, 2.88700,
1962         2.85009, 2.83177, 2.83723, 2.87750, 3.14193},
1963        {2.84574, 2.84789, 2.85560, 2.86935, 2.89075, 2.98060,
1964         2.94576, 2.92868, 2.93356, 2.97436, 3.23868},
1965        {2.94239, 2.94469, 2.95221, 2.96625, 2.99345, 3.07747,
1966         3.04266, 3.02545, 3.03051, 3.07118, 3.33555}};
1967
1968   // Garfield simulation at UD = -1600V; vd = 1.158cm/microsec, <driftfield> = 558V/cm
1969   Float_t time1600[ktimebin][kZbin] =
1970       {{0.09169, 0.09203, 0.09305, 0.09473, 0.09714, 0.10032,
1971         0.10441, 0.10990, 0.11835, 0.13986, 0.37845},
1972        {0.06589, 0.06626, 0.06738, 0.06924, 0.07184, 0.07521,
1973         0.07947, 0.08512, 0.09379, 0.11603, 0.35648},
1974        {0.03947, 0.04003, 0.04171, 0.04434, 0.04778, 0.05190,
1975         0.05678, 0.06304, 0.07292, 0.09876, 0.34759},
1976        {0.01111, 0.01281, 0.01718, 0.02281, 0.02879, 0.03477,
1977         0.04097, 0.04910, 0.06415, 0.10573, 0.36896},
1978        {0.01120, 0.01311, 0.01721, 0.02279, 0.02862, 0.03446,
1979         0.04074, 0.04981, 0.06825, 0.11595, 0.37255},
1980        {0.03919, 0.03980, 0.04132, 0.04380, 0.04704, 0.05102,
1981         0.05602, 0.06302, 0.07633, 0.10896, 0.36743},
1982        {0.06531, 0.06528, 0.06631, 0.06805, 0.07053, 0.07392,
1983         0.07853, 0.08505, 0.09669, 0.12578, 0.37967},
1984        {0.09157, 0.09171, 0.09216, 0.09301, 0.09475, 0.00000,
1985         0.11152, 0.11879, 0.13352, 0.16802, 0.42750},
1986        {0.13977, 0.14095, 0.14509, 0.15433, 0.17534, 0.26406,
1987         0.21660, 0.20345, 0.21113, 0.25067, 0.51434},
1988        {0.21816, 0.22041, 0.22631, 0.23850, 0.26208, 0.34340,
1989         0.30755, 0.29237, 0.29878, 0.33863, 0.60258},
1990        {0.30344, 0.30547, 0.31241, 0.32444, 0.34809, 0.42696,
1991         0.39464, 0.37919, 0.38546, 0.42530, 0.68926},
1992        {0.38969, 0.39164, 0.39810, 0.41059, 0.43441, 0.51246,
1993         0.48112, 0.46562, 0.47191, 0.51172, 0.77558},
1994        {0.47592, 0.47799, 0.48442, 0.49689, 0.52061, 0.59855,
1995         0.56752, 0.55201, 0.55826, 0.59808, 0.86202},
1996        {0.56226, 0.56428, 0.57074, 0.58324, 0.60696, 0.68483,
1997         0.65388, 0.63837, 0.64461, 0.68445, 0.94830},
1998        {0.64881, 0.65063, 0.65709, 0.66958, 0.69331, 0.77117,
1999         0.74023, 0.72472, 0.73098, 0.77079, 1.03486},
2000        {0.73506, 0.73698, 0.74344, 0.75596, 0.77964, 0.85751,
2001         0.82658, 0.81107, 0.81731, 0.85712, 1.12106},
2002        {0.82132, 0.82333, 0.82979, 0.84228, 0.86608, 0.94386,
2003         0.91293, 0.89742, 0.90367, 0.94335, 1.20737},
2004        {0.90767, 0.90968, 0.91614, 0.92863, 0.95236, 1.03021,
2005         0.99928, 0.98377, 0.99001, 1.02984, 1.29371},
2006        {0.99410, 0.99602, 1.00257, 1.01498, 1.03869, 1.11720,
2007         1.08563, 1.07011, 1.07637, 1.11621, 1.37873},
2008        {1.08036, 1.08240, 1.08884, 1.10138, 1.12504, 1.20301,
2009         1.17198, 1.15647, 1.16272, 1.20255, 1.46651},
2010        {1.16670, 1.16872, 1.17525, 1.18783, 1.21139, 1.28934,
2011         1.25833, 1.24281, 1.24909, 1.28889, 1.55275},
2012        {1.25307, 1.25510, 1.26153, 1.27404, 1.29773, 1.37584,
2013         1.34469, 1.32916, 1.33536, 1.37524, 1.63915},
2014        {1.33942, 1.34146, 1.34788, 1.36040, 1.38410, 1.46438,
2015         1.43105, 1.41537, 1.42176, 1.46158, 1.72538},
2016        {1.42579, 1.42782, 1.43458, 1.44674, 1.47043, 1.55085,
2017         1.51675, 1.50168, 1.50810, 1.54793, 1.81174},
2018        {1.51207, 1.51454, 1.52060, 1.53307, 1.55684, 1.63478,
2019         1.60336, 1.58820, 1.59446, 1.63414, 1.89814},
2020        {1.59856, 1.60047, 1.60693, 1.61942, 1.64317, 1.72257,
2021         1.69008, 1.67454, 1.68080, 1.72063, 1.98433},
2022        {1.68481, 1.68682, 1.69330, 1.70584, 1.72949, 1.80752,
2023         1.77643, 1.76089, 1.76716, 1.80692, 2.07069},
2024        {1.77117, 1.77319, 1.77969, 1.79260, 1.81583, 1.89376,
2025         1.86226, 1.84720, 1.85355, 1.89256, 2.15343},
2026        {1.85748, 1.85967, 1.86605, 1.87848, 1.90222, 1.98010,
2027         1.94913, 1.93271, 1.93981, 1.97968, 2.24355},
2028        {1.94386, 1.94587, 1.95233, 1.96484, 1.98854, 2.06646,
2029         2.03542, 2.01755, 2.02617, 2.06604, 2.32993},
2030        {2.03022, 2.03230, 2.03868, 2.05134, 2.07488, 2.15367,
2031         2.12178, 2.10391, 2.11252, 2.15432, 2.41623},
2032        {2.11656, 2.11857, 2.12505, 2.13772, 2.16147, 2.23919,
2033         2.20817, 2.19265, 2.20744, 2.23872, 2.49996},
2034        {2.20291, 2.20611, 2.21137, 2.22387, 2.24758, 2.32563,
2035         2.29450, 2.27901, 2.28525, 2.32507, 2.58897},
2036        {2.28922, 2.29172, 2.29774, 2.31345, 2.33400, 2.41287,
2037         2.38086, 2.36535, 2.37160, 2.40869, 2.67113},
2038        {2.37572, 2.37764, 2.38410, 2.39803, 2.42046, 2.49817,
2039         2.46721, 2.45171, 2.45794, 2.49505, 2.76061},
2040        {2.46190, 2.46396, 2.47043, 2.48340, 2.50665, 2.58453,
2041         2.55357, 2.53728, 2.54430, 2.58407, 2.84816},
2042        {2.54833, 2.55032, 2.55679, 2.56976, 2.59312, 2.67103,
2043         2.63993, 2.62364, 2.63062, 2.67040, 2.93444},
2044        {2.63456, 2.63660, 2.64304, 2.65555, 2.67938, 2.75739,
2045         2.72629, 2.71064, 2.71688, 2.75671, 3.01886}};
2046
2047   // Garfield simulation at UD = -1700V; vd = 1.299cm/microsec, <driftfield> = 590V/cm
2048   Float_t time1700[ktimebin][kZbin] =
2049       {{0.09167, 0.09201, 0.09302, 0.09471, 0.09712, 0.10029,
2050         0.10438, 0.10986, 0.11832, 0.13986, 0.37824},
2051        {0.06591, 0.06626, 0.06736, 0.06923, 0.07183, 0.07519,
2052         0.07944, 0.08511, 0.09378, 0.11603, 0.35625},
2053        {0.03946, 0.04003, 0.04170, 0.04433, 0.04777, 0.05189,
2054         0.05676, 0.06301, 0.07291, 0.09880, 0.34724},
2055        {0.01110, 0.01281, 0.01718, 0.02280, 0.02879, 0.03476,
2056         0.04096, 0.04910, 0.06417, 0.10582, 0.36861},
2057        {0.01115, 0.01294, 0.01721, 0.02276, 0.02862, 0.03447,
2058         0.04074, 0.04980, 0.06812, 0.11565, 0.37231},
2059        {0.03920, 0.03974, 0.04133, 0.04381, 0.04706, 0.05105,
2060         0.05603, 0.06299, 0.07618, 0.10860, 0.36646},
2061        {0.06498, 0.06529, 0.06634, 0.06808, 0.07055, 0.07395,
2062         0.07852, 0.08500, 0.09650, 0.12532, 0.37850},
2063        {0.09143, 0.09159, 0.09207, 0.09297, 0.09473, 0.00000,
2064         0.11102, 0.11809, 0.13245, 0.16627, 0.42496},
2065        {0.13646, 0.13750, 0.14112, 0.14926, 0.16806, 0.24960,
2066         0.20627, 0.19536, 0.20366, 0.24256, 0.50557},
2067        {0.20678, 0.20848, 0.21384, 0.22450, 0.24552, 0.32018,
2068         0.28740, 0.27477, 0.28196, 0.32128, 0.58475},
2069        {0.28287, 0.28461, 0.29020, 0.30108, 0.32224, 0.39467,
2070         0.36500, 0.35217, 0.35926, 0.39860, 0.66194},
2071        {0.35972, 0.36145, 0.36713, 0.37797, 0.39912, 0.47091,
2072         0.44212, 0.42925, 0.43632, 0.47563, 0.73892},
2073        {0.43667, 0.43841, 0.44413, 0.45494, 0.47607, 0.54780,
2074         0.51912, 0.50627, 0.51334, 0.55254, 0.81595},
2075        {0.51365, 0.51540, 0.52101, 0.53193, 0.55305, 0.62463,
2076         0.59617, 0.58328, 0.59035, 0.62965, 0.89303},
2077        {0.59064, 0.59240, 0.59801, 0.60893, 0.63009, 0.70169,
2078         0.67317, 0.66028, 0.66735, 0.70666, 0.96995},
2079        {0.66765, 0.66939, 0.67501, 0.68592, 0.70724, 0.77863,
2080         0.75016, 0.73728, 0.74435, 0.78366, 1.04696},
2081        {0.74464, 0.74636, 0.75200, 0.76293, 0.78405, 0.85561,
2082         0.82716, 0.81427, 0.82133, 0.86064, 1.12396},
2083        {0.82165, 0.82340, 0.82902, 0.83991, 0.86104, 0.93266,
2084         0.90414, 0.89128, 0.89834, 0.93763, 1.20100},
2085        {0.89863, 0.90042, 0.90659, 0.91705, 0.93805, 1.00960,
2086         0.98115, 0.96825, 0.97533, 1.01462, 1.27801},
2087        {0.97563, 0.97740, 0.98310, 0.99391, 1.01504, 1.08659,
2088         1.05814, 1.04526, 1.05233, 1.09163, 1.35503},
2089        {1.05276, 1.05451, 1.06002, 1.07090, 1.09099, 1.16357,
2090         1.13516, 1.12225, 1.12933, 1.16863, 1.43195},
2091        {1.12977, 1.13138, 1.13700, 1.14792, 1.16797, 1.24061,
2092         1.21212, 1.19926, 1.20626, 1.24554, 1.50900},
2093        {1.20664, 1.20839, 1.21400, 1.22490, 1.24606, 1.31772,
2094         1.28914, 1.27382, 1.28329, 1.32262, 1.58550},
2095        {1.28367, 1.28541, 1.29099, 1.30189, 1.32312, 1.39460,
2096         1.36612, 1.34924, 1.36030, 1.39961, 1.66310},
2097        {1.36064, 1.36249, 1.36799, 1.37896, 1.40004, 1.48030,
2098         1.44314, 1.43032, 1.43731, 1.47659, 1.73442},
2099        {1.43762, 1.43937, 1.44497, 1.45618, 1.47704, 1.54932,
2100         1.52012, 1.50725, 1.51430, 1.55357, 1.81708},
2101        {1.51462, 1.51937, 1.52203, 1.53316, 1.55403, 1.62572,
2102         1.59713, 1.58424, 1.59128, 1.63061, 1.89406},
2103        {1.59162, 1.59338, 1.59947, 1.60989, 1.63103, 1.70270,
2104         1.67411, 1.66124, 1.66799, 1.70759, 1.97103},
2105        {1.66874, 1.67037, 1.67597, 1.68687, 1.70814, 1.77969,
2106         1.75112, 1.73806, 1.74530, 1.78457, 2.04794},
2107        {1.74693, 1.74749, 1.75297, 1.76476, 1.78500, 1.85667,
2108         1.82811, 1.81504, 1.82101, 1.86161, 2.12492},
2109        {1.82260, 1.82437, 1.82995, 1.84174, 1.86202, 1.93372,
2110         1.90509, 1.89202, 1.89930, 1.93859, 2.20189},
2111        {1.89964, 1.90135, 1.90693, 1.91789, 1.93952, 2.01080,
2112         1.98207, 1.96921, 1.97628, 2.01384, 2.27887},
2113        {1.97662, 1.97917, 1.98611, 1.99487, 2.01601, 2.08778,
2114         2.05846, 2.04623, 2.05330, 2.09244, 2.35585},
2115        {2.05359, 2.05615, 2.06309, 2.07187, 2.09867, 2.16459,
2116         2.13610, 2.12322, 2.13029, 2.16942, 2.43199},
2117        {2.13063, 2.13233, 2.13795, 2.14886, 2.17008, 2.24199,
2118         2.21310, 2.20020, 2.20727, 2.24659, 2.50983},
2119        {2.20761, 2.20931, 2.21955, 2.22624, 2.24708, 2.32147,
2120         2.29009, 2.27725, 2.28276, 2.32359, 2.58680},
2121        {2.28459, 2.29108, 2.29202, 2.30286, 2.32007, 2.39559,
2122         2.36683, 2.35422, 2.36119, 2.40058, 2.66081},
2123        {2.36153, 2.36806, 2.36889, 2.37985, 2.40092, 2.47828,
2124         2.44381, 2.43099, 2.43819, 2.47750, 2.73779}};
2125
2126   // Garfield simulation at UD = -1800V; vd = 1.450cm/microsec, <driftfield> = 623V/cm
2127   Float_t time1800[ktimebin][kZbin] =
2128       {{0.09166, 0.09199, 0.09300, 0.09470, 0.09709, 0.10026,
2129         0.10434, 0.10983, 0.11831, 0.13987, 0.37802},
2130        {0.06585, 0.06623, 0.06735, 0.06921, 0.07180, 0.07520,
2131         0.07941, 0.08507, 0.09376, 0.11604, 0.35624},
2132        {0.03945, 0.04004, 0.04169, 0.04432, 0.04776, 0.05187,
2133         0.05674, 0.06300, 0.07290, 0.09884, 0.34704},
2134        {0.01108, 0.01287, 0.01717, 0.02280, 0.02880, 0.03476,
2135         0.04095, 0.04909, 0.06419, 0.10589, 0.36843},
2136        {0.01115, 0.01287, 0.01720, 0.02276, 0.02862, 0.03448,
2137         0.04073, 0.04973, 0.06799, 0.11535, 0.37224},
2138        {0.03918, 0.03975, 0.04134, 0.04382, 0.04707, 0.05105,
2139         0.05603, 0.06296, 0.07605, 0.10822, 0.36560},
2140        {0.06498, 0.06532, 0.06635, 0.06809, 0.07058, 0.07395,
2141         0.07855, 0.08495, 0.09632, 0.12488, 0.37730},
2142        {0.09130, 0.09160, 0.09199, 0.09300, 0.09472, 0.00000,
2143         0.11059, 0.11747, 0.13146, 0.16462, 0.42233},
2144        {0.13364, 0.13449, 0.13767, 0.14481, 0.16147, 0.23635,
2145         0.19706, 0.18812, 0.19704, 0.23520, 0.49749},
2146        {0.19698, 0.19844, 0.20311, 0.21236, 0.23082, 0.29920,
2147         0.26936, 0.25927, 0.26732, 0.30601, 0.56871},
2148        {0.26518, 0.26670, 0.27160, 0.28099, 0.29955, 0.36597,
2149         0.33885, 0.32858, 0.33653, 0.37524, 0.63801},
2150        {0.33441, 0.33553, 0.34040, 0.34987, 0.36841, 0.43428,
2151         0.40797, 0.39763, 0.40556, 0.44425, 0.70698},
2152        {0.40296, 0.40447, 0.40934, 0.41881, 0.43737, 0.50306,
2153         0.47695, 0.46662, 0.47455, 0.51329, 0.77600},
2154        {0.47296, 0.47344, 0.47830, 0.48779, 0.50632, 0.57200,
2155         0.54593, 0.53559, 0.54351, 0.58222, 0.84489},
2156        {0.54089, 0.54264, 0.54727, 0.55673, 0.57529, 0.64094,
2157         0.61490, 0.60457, 0.61249, 0.65118, 0.91394},
2158        {0.60987, 0.61138, 0.61624, 0.62573, 0.64428, 0.70989,
2159         0.68397, 0.67354, 0.68147, 0.72015, 0.98291},
2160        {0.67883, 0.68035, 0.68521, 0.69469, 0.71324, 0.77896,
2161         0.75287, 0.74251, 0.75043, 0.78912, 1.04458},
2162        {0.74780, 0.74932, 0.75421, 0.76367, 0.78221, 0.84785,
2163         0.82185, 0.81148, 0.81941, 0.85811, 1.12085},
2164        {0.81690, 0.81830, 0.82316, 0.83263, 0.85120, 0.91683,
2165         0.89077, 0.88045, 0.88837, 0.92707, 1.18976},
2166        {0.88574, 0.88726, 0.89228, 0.90198, 0.92017, 0.98578,
2167         0.95974, 0.94947, 0.95734, 0.99604, 1.25873},
2168        {0.95493, 0.95624, 0.96110, 0.97058, 0.98913, 1.05481,
2169         1.02873, 1.01839, 1.02631, 1.06503, 1.32772},
2170        {1.02392, 1.02524, 1.03008, 1.03955, 1.05810, 1.12378,
2171         1.09757, 1.08605, 1.09530, 1.13399, 1.39669},
2172        {1.09270, 1.09418, 1.09911, 1.10854, 1.12714, 1.19281,
2173         1.16502, 1.15633, 1.16427, 1.20271, 1.46574},
2174        {1.16168, 1.16323, 1.16801, 1.17772, 1.19604, 1.26190,
2175         1.23399, 1.22531, 1.23323, 1.27194, 1.53475},
2176        {1.23061, 1.23214, 1.23698, 1.24669, 1.26503, 1.33073,
2177         1.30461, 1.29428, 1.30220, 1.34091, 1.60372},
2178        {1.29960, 1.30110, 1.30596, 1.31544, 1.33398, 1.39962,
2179         1.37228, 1.36323, 1.37121, 1.40988, 1.67273},
2180        {1.36851, 1.37007, 1.37512, 1.38441, 1.40297, 1.46865,
2181         1.44256, 1.43222, 1.44017, 1.47878, 1.74155},
2182        {1.43752, 1.43904, 1.44773, 1.45338, 1.47220, 1.53759,
2183         1.51136, 1.50119, 1.50914, 1.54775, 1.81050},
2184        {1.50646, 1.50802, 1.51288, 1.52237, 1.54097, 1.60697,
2185         1.58049, 1.57018, 1.57811, 1.61678, 1.87947},
2186        {1.57545, 1.57720, 1.58185, 1.59134, 1.60996, 1.67787,
2187         1.64929, 1.63914, 1.64707, 1.68570, 1.94851},
2188        {1.64442, 1.64617, 1.65081, 1.66035, 1.67893, 1.74684,
2189         1.71826, 1.70745, 1.71604, 1.75310, 2.01748},
2190        {1.71337, 1.71513, 1.71978, 1.72932, 1.74645, 1.81346,
2191         1.78739, 1.77642, 1.78501, 1.82151, 2.08644},
2192        {1.78238, 1.78410, 1.78876, 1.79824, 1.81678, 1.88332,
2193         1.85639, 1.84262, 1.85397, 1.89270, 2.15533},
2194        {1.85135, 1.85306, 1.85778, 1.86728, 1.88580, 1.95615,
2195         1.92536, 1.91171, 1.92283, 1.96165, 2.22428},
2196        {1.92774, 1.92184, 1.92672, 1.93618, 1.95477, 2.02048,
2197         1.99427, 1.98068, 1.99192, 2.03062, 2.29338},
2198        {1.98929, 1.99081, 1.99567, 2.00515, 2.02373, 2.08987,
2199         2.06332, 2.05249, 2.05939, 2.09928, 2.36227},
2200        {2.05829, 2.05978, 2.06464, 2.07414, 2.09272, 2.15850,
2201         2.12928, 2.12194, 2.12987, 2.16825, 2.43083},
2202        {2.12726, 2.12869, 2.13360, 2.14425, 2.16160, 2.22872,
2203         2.20118, 2.19078, 2.19876, 2.23416, 2.50007}};
2204
2205   // Garfield simulation at UD = -1900V; vd = 1.610cm/microsec, <driftfield> = 655V/cm
2206   Float_t time1900[ktimebin][kZbin] =
2207       {{0.09166, 0.09198, 0.09298, 0.09467, 0.09707, 0.10023,
2208         0.10431, 0.10980, 0.11828, 0.13988, 0.37789},
2209        {0.06584, 0.06622, 0.06735, 0.06920, 0.07179, 0.07514,
2210         0.07938, 0.08505, 0.09374, 0.11606, 0.35599},
2211        {0.03945, 0.04002, 0.04169, 0.04432, 0.04775, 0.05185,
2212         0.05672, 0.06298, 0.07290, 0.09888, 0.34695},
2213        {0.01109, 0.01281, 0.01717, 0.02279, 0.02878, 0.03476,
2214         0.04094, 0.04909, 0.06421, 0.10597, 0.36823},
2215        {0.01115, 0.01287, 0.01720, 0.02294, 0.02862, 0.03448,
2216         0.04074, 0.04973, 0.06783, 0.11506, 0.37206},
2217        {0.03940, 0.03975, 0.04135, 0.04386, 0.04708, 0.05106,
2218         0.05604, 0.06293, 0.07592, 0.10787, 0.36484},
2219        {0.06500, 0.06534, 0.06638, 0.06811, 0.07060, 0.07413,
2220         0.07852, 0.08491, 0.09614, 0.12446, 0.37626},
2221        {0.09119, 0.09140, 0.09194, 0.09293, 0.09471, 0.00000,
2222         0.11013, 0.11685, 0.13050, 0.16302, 0.41991},
2223        {0.13111, 0.13190, 0.13466, 0.14091, 0.15554, 0.22409,
2224         0.18846, 0.18167, 0.19113, 0.22854, 0.48995},
2225        {0.18849, 0.18975, 0.19380, 0.20185, 0.21797, 0.28050,
2226         0.25368, 0.24575, 0.25446, 0.29249, 0.55442},
2227        {0.24995, 0.25125, 0.25563, 0.26366, 0.27986, 0.34065,
2228         0.31605, 0.30815, 0.31680, 0.35482, 0.61697},
2229        {0.31187, 0.31324, 0.31745, 0.32580, 0.34205, 0.40217,
2230         0.37825, 0.37031, 0.37897, 0.41696, 0.67890},
2231        {0.37401, 0.37531, 0.37955, 0.38777, 0.40395, 0.46408,
2232         0.44037, 0.43242, 0.44108, 0.47906, 0.74122},
2233        {0.43610, 0.43741, 0.44161, 0.44986, 0.46604, 0.52614,
2234         0.50248, 0.49452, 0.50316, 0.54116, 0.80326},
2235        {0.49820, 0.49988, 0.50372, 0.51196, 0.52814, 0.58822,
2236         0.56459, 0.55661, 0.56527, 0.60326, 0.86526},
2237        {0.56032, 0.56161, 0.56582, 0.57408, 0.59024, 0.65032,
2238         0.62670, 0.61872, 0.62737, 0.66537, 0.92749},
2239        {0.62240, 0.62371, 0.62792, 0.63624, 0.65236, 0.71241,
2240         0.68881, 0.68081, 0.68947, 0.72750, 0.98941},
2241        {0.68449, 0.68581, 0.69002, 0.69828, 0.71444, 0.77452,
2242         0.75089, 0.74295, 0.75157, 0.78957, 1.05157},
2243        {0.74660, 0.74790, 0.75212, 0.76036, 0.77654, 0.83748,
2244         0.81299, 0.80501, 0.81193, 0.85168, 1.11375},
2245        {0.80870, 0.81017, 0.81423, 0.82246, 0.83867, 0.89908,
2246         0.87509, 0.86660, 0.87577, 0.91376, 1.17586},
2247        {0.87080, 0.87233, 0.87632, 0.88458, 0.90074, 0.96083,
2248         0.93718, 0.92922, 0.93787, 0.97588, 1.23794},
2249        {0.93291, 0.93422, 0.93844, 0.94667, 0.96293, 1.02295,
2250         0.99929, 0.99127, 0.99997, 1.03797, 1.29995},
2251        {0.99500, 0.99645, 1.00308, 1.00877, 1.02497, 1.08504,
2252         1.06140, 1.05343, 1.06203, 1.10006, 1.36216},
2253        {1.05712, 1.05926, 1.06262, 1.07092, 1.08706, 1.14754,
2254         1.12350, 1.11550, 1.12417, 1.16218, 1.42427},
2255        {1.11921, 1.12059, 1.12473, 1.13297, 1.14916, 1.21140,
2256         1.18560, 1.17284, 1.18625, 1.22414, 1.48629},
2257        {1.18140, 1.18262, 1.18690, 1.19508, 1.21125, 1.27139,
2258         1.24164, 1.23495, 1.24838, 1.28634, 1.54852},
2259        {1.24340, 1.24473, 1.24901, 1.25732, 1.27336, 1.33358,
2260         1.30793, 1.30179, 1.31047, 1.34848, 1.61066},
2261        {1.30551, 1.30684, 1.31104, 1.32056, 1.33553, 1.39609,
2262         1.37004, 1.36392, 1.37045, 1.41057, 1.67259},
2263        {1.36755, 1.36892, 1.37315, 1.39148, 1.39755, 1.45820,
2264         1.43215, 1.42602, 1.43467, 1.47268, 1.73477},
2265        {1.42966, 1.43101, 1.43549, 1.45359, 1.45976, 1.52031,
2266         1.49601, 1.48811, 1.49677, 1.53477, 1.79691},
2267        {1.49180, 1.49321, 1.49760, 1.51570, 1.52175, 1.58185,
2268         1.55771, 1.55023, 1.55888, 1.59672, 1.85501},
2269        {1.55391, 1.55527, 1.55943, 1.57782, 1.58391, 1.64395,
2270         1.62008, 1.61233, 1.62085, 1.65883, 1.92091},
2271        {1.61599, 1.61732, 1.62154, 1.63993, 1.64612, 1.70608,
2272         1.68237, 1.67108, 1.68301, 1.72110, 1.97931},
2273        {1.67808, 1.67948, 1.68364, 1.70204, 1.70823, 1.76858,
2274         1.74404, 1.73539, 1.74512, 1.78321, 2.04522},
2275        {1.74019, 1.74152, 1.74573, 1.76415, 1.77015, 1.83040,
2276         1.80615, 1.79366, 1.80723, 1.84509, 2.10742},
2277        {1.80235, 1.80362, 1.80783, 1.82626, 1.83227, 1.89246,
2278         1.86795, 1.85405, 1.86938, 1.90720, 2.16953},
2279        {1.86442, 1.86572, 1.86994, 1.88837, 1.89438, 1.95445,
2280         1.93006, 1.92283, 1.93148, 1.96931, 2.23147},
2281        {1.92700, 1.92783, 1.93197, 1.95049, 1.95649, 2.01668,
2282         1.99217, 1.98486, 1.99352, 2.03143, 2.29358}};
2283
2284   // Garfield simulation at UD = -2000V; vd = 1.783cm/microsec, <driftfield> = 688V/cm
2285   Float_t time2000[ktimebin][kZbin] =
2286       {{0.09176, 0.09196, 0.09296, 0.09465, 0.09704, 0.10020,
2287         0.10427, 0.10977, 0.11825, 0.13988, 0.37774},
2288        {0.06583, 0.06620, 0.06735, 0.06918, 0.07177, 0.07513,
2289         0.07936, 0.08503, 0.09372, 0.11606, 0.35586},
2290        {0.03944, 0.04001, 0.04170, 0.04431, 0.04774, 0.05184,
2291         0.05670, 0.06296, 0.07291, 0.09893, 0.34680},
2292        {0.01108, 0.01281, 0.01719, 0.02279, 0.02879, 0.03474,
2293         0.04093, 0.04908, 0.06422, 0.10605, 0.36800},
2294        {0.01114, 0.01287, 0.01720, 0.02276, 0.02863, 0.03449,
2295         0.04073, 0.04970, 0.06774, 0.11478, 0.37179},
2296        {0.03925, 0.03977, 0.04135, 0.04386, 0.04711, 0.05108,
2297         0.05604, 0.06290, 0.07580, 0.10748, 0.36386},
2298        {0.06501, 0.06536, 0.06640, 0.06814, 0.07062, 0.07398,
2299         0.07852, 0.08487, 0.09598, 0.12405, 0.37519},
2300        {0.09109, 0.09127, 0.09188, 0.09292, 0.09472, 0.00000,
2301         0.10964, 0.11630, 0.12960, 0.16150, 0.41765},
2302        {0.12898, 0.12968, 0.13209, 0.13749, 0.15034, 0.21286,
2303         0.18088, 0.17590, 0.18591, 0.22254, 0.48315},
2304        {0.18122, 0.18227, 0.18574, 0.19263, 0.20674, 0.26376,
2305         0.23960, 0.23375, 0.24316, 0.28047, 0.54179},
2306        {0.23674, 0.23784, 0.24142, 0.24847, 0.26264, 0.31810,
2307         0.29602, 0.29008, 0.29944, 0.33675, 0.59795},
2308        {0.29279, 0.29382, 0.29742, 0.30448, 0.31865, 0.37364,
2309         0.35215, 0.34629, 0.35555, 0.39286, 0.65411},
2310        {0.34875, 0.34987, 0.35346, 0.36054, 0.37472, 0.42956,
2311         0.40825, 0.40229, 0.41167, 0.44894, 0.71033},
2312        {0.40484, 0.40594, 0.40954, 0.41660, 0.43077, 0.48560,
2313         0.46433, 0.45840, 0.46772, 0.50500, 0.76632},
2314        {0.46090, 0.46201, 0.46560, 0.47267, 0.48684, 0.54167,
2315         0.52041, 0.51449, 0.52382, 0.56108, 0.82227},
2316        {0.51698, 0.51809, 0.52173, 0.52874, 0.54291, 0.59776,
2317         0.57646, 0.57052, 0.57986, 0.61717, 0.87836},
2318        {0.57306, 0.57418, 0.57782, 0.58513, 0.59899, 0.65380,
2319         0.63255, 0.62661, 0.63594, 0.67325, 0.93460},
2320        {0.62912, 0.63024, 0.63383, 0.64103, 0.65506, 0.70988,
2321         0.68484, 0.68267, 0.69202, 0.72878, 0.99046},
2322        {0.68521, 0.68633, 0.68990, 0.69699, 0.71115, 0.76595,
2323         0.74468, 0.73872, 0.74814, 0.78538, 1.04674},
2324        {0.74127, 0.74239, 0.74605, 0.75303, 0.77022, 0.82204,
2325         0.80078, 0.79484, 0.80416, 0.84147, 1.10261},
2326        {0.79736, 0.79846, 0.80206, 0.80947, 0.82330, 0.87813,
2327         0.85688, 0.85087, 0.86023, 0.89752, 1.15874},
2328        {0.85342, 0.85454, 0.85815, 0.86519, 0.87936, 0.93417,
2329         0.91293, 0.90428, 0.91631, 0.95360, 1.20760},
2330        {0.90949, 0.91061, 0.91423, 0.92128, 0.93544, 0.99026,
2331         0.96807, 0.96305, 0.97239, 1.00967, 1.27078},
2332        {0.96556, 0.96669, 0.97111, 0.97734, 0.99151, 1.04664,
2333         1.02508, 1.01879, 1.02846, 1.06167, 1.32695},
2334        {1.02167, 1.02279, 1.02656, 1.03341, 1.04759, 1.10242,
2335         1.08115, 1.07003, 1.08453, 1.12184, 1.38304},
2336        {1.07776, 1.07883, 1.08242, 1.08950, 1.10384, 1.16422,
2337         1.13725, 1.13133, 1.14061, 1.17793, 1.43910},
2338        {1.13379, 1.13492, 1.13864, 1.14567, 1.15973, 1.21455,
2339         1.19323, 1.18734, 1.19668, 1.23401, 1.49528},
2340        {1.18988, 1.19098, 1.19457, 1.20164, 1.21582, 1.27064,
2341         1.24937, 1.24044, 1.25275, 1.29004, 1.55137},
2342        {1.24592, 1.24706, 1.25087, 1.25773, 1.27188, 1.32670,
2343         1.30544, 1.29953, 1.30883, 1.34613, 1.60743},
2344        {1.30202, 1.30313, 1.30673, 1.31381, 1.32797, 1.38278,
2345         1.36151, 1.35167, 1.36490, 1.40221, 1.66306},
2346        {1.35809, 1.35921, 1.36282, 1.36986, 1.38403, 1.43888,
2347         1.41760, 1.41174, 1.42083, 1.45830, 1.71915},
2348        {1.41419, 1.41528, 1.41890, 1.42595, 1.44011, 1.49496,
2349         1.47368, 1.46769, 1.47706, 1.51436, 1.77523},
2350        {1.47131, 1.47141, 1.47494, 1.48850, 1.49620, 1.55137,
2351         1.52977, 1.51820, 1.53315, 1.57042, 1.83158},
2352        {1.52635, 1.52750, 1.53103, 1.53814, 1.55228, 1.60736,
2353         1.58503, 1.57986, 1.58920, 1.62649, 1.88767},
2354        {1.58418, 1.58355, 1.58711, 1.59526, 1.60833, 1.66316,
2355         1.63345, 1.63261, 1.64556, 1.68204, 1.94359},
2356        {1.64027, 1.63958, 1.64489, 1.65024, 1.66443, 1.71925,
2357         1.69794, 1.69201, 1.70143, 1.73865, 1.99968},
2358        {1.69450, 1.69566, 1.69940, 1.70697, 1.71841, 1.77819,
2359         1.75396, 1.74814, 1.75743, 1.79083, 2.05427},
2360        {1.75054, 1.75221, 1.75527, 1.76306, 1.77662, 1.83428,
2361         1.81006, 1.81173, 1.81345, 1.85076, 2.10289}};
2362
2363   // Garfield simulation at UD = -2100V; vd = 1.959cm/microsec, <driftfield> = 720V/cm
2364   Float_t time2100[ktimebin][kZbin] =
2365       {{0.09160, 0.09194, 0.09294, 0.09462, 0.09701, 0.10017,
2366         0.10424, 0.10974, 0.11823, 0.13988, 0.37762},
2367        {0.06585, 0.06619, 0.06731, 0.06916, 0.07174, 0.07509,
2368         0.07933, 0.08500, 0.09370, 0.11609, 0.35565},
2369        {0.03960, 0.04001, 0.04171, 0.04430, 0.04774, 0.05182,
2370         0.05668, 0.06294, 0.07291, 0.09896, 0.34676},
2371        {0.01109, 0.01280, 0.01716, 0.02279, 0.02876, 0.03474,
2372         0.04096, 0.04908, 0.06424, 0.10612, 0.36790},
2373        {0.01114, 0.01285, 0.01719, 0.02287, 0.02863, 0.03449,
2374         0.04073, 0.04964, 0.06759, 0.11446, 0.37162},
2375        {0.03922, 0.03977, 0.04146, 0.04386, 0.04711, 0.05109,
2376         0.05605, 0.06287, 0.07575, 0.10713, 0.36298},
2377        {0.06504, 0.06538, 0.06641, 0.06818, 0.07064, 0.07426,
2378         0.07852, 0.08483, 0.09581, 0.12363, 0.37424},
2379        {0.09103, 0.09129, 0.09186, 0.09291, 0.09476, 0.00000,
2380         0.10923, 0.11578, 0.12873, 0.16005, 0.41525},
2381        {0.12723, 0.12777, 0.12988, 0.13458, 0.14579, 0.20264,
2382         0.17421, 0.17078, 0.18132, 0.21708, 0.47699},
2383        {0.17508, 0.17601, 0.17897, 0.18487, 0.19698, 0.24881,
2384         0.22737, 0.22337, 0.23348, 0.27000, 0.53032},
2385        {0.22571, 0.22663, 0.22969, 0.23570, 0.24787, 0.29826,
2386         0.27871, 0.27462, 0.28471, 0.32122, 0.58166},
2387        {0.27664, 0.27759, 0.28067, 0.28669, 0.29891, 0.34898,
2388         0.32982, 0.32570, 0.33576, 0.37229, 0.63268},
2389        {0.32766, 0.32862, 0.33170, 0.33778, 0.34988, 0.39973,
2390         0.38088, 0.37675, 0.38680, 0.42333, 0.68159},
2391        {0.37872, 0.37966, 0.38275, 0.38875, 0.40093, 0.45073,
2392         0.43192, 0.42780, 0.43786, 0.47438, 0.73480},
2393        {0.42974, 0.43070, 0.43378, 0.43982, 0.45196, 0.50177,
2394         0.48297, 0.47884, 0.48889, 0.52544, 0.78581},
2395        {0.48081, 0.48175, 0.48482, 0.49084, 0.50302, 0.55290,
2396         0.53398, 0.52988, 0.53994, 0.57647, 0.83687},
2397        {0.53645, 0.53295, 0.53586, 0.54188, 0.55408, 0.60398,
2398         0.58504, 0.58092, 0.59100, 0.62768, 0.88773},
2399        {0.58345, 0.58409, 0.58690, 0.59292, 0.60510, 0.65562,
2400         0.63609, 0.63197, 0.64203, 0.67856, 0.93907},
2401        {0.63397, 0.63490, 0.63795, 0.64403, 0.65613, 0.70612,
2402         0.68714, 0.68301, 0.69294, 0.72955, 0.99000},
2403        {0.68496, 0.68592, 0.68899, 0.69504, 0.70733, 0.75708,
2404         0.73818, 0.73405, 0.74412, 0.78064, 1.04100},
2405        {0.73600, 0.73696, 0.74003, 0.74624, 0.75828, 0.80805,
2406         0.78904, 0.78512, 0.79517, 0.83152, 1.09205},
2407        {0.78709, 0.78801, 0.79108, 0.79709, 0.80931, 0.85906,
2408         0.84027, 0.83614, 0.84621, 0.88269, 1.14058},
2409        {0.83808, 0.83905, 0.84215, 0.84816, 0.86031, 0.91011,
2410         0.89139, 0.88718, 0.89725, 0.93377, 1.19413},
2411        {0.88916, 0.89010, 0.89320, 0.89920, 0.91136, 0.96117,
2412         0.94235, 0.93822, 0.94828, 0.98480, 1.24538},
2413        {0.94036, 0.94113, 0.94422, 0.95023, 0.96241, 1.01220,
2414         0.99310, 0.98927, 0.99933, 1.03585, 1.29629},
2415        {0.99139, 0.99220, 0.99525, 1.00127, 1.01344, 1.06324,
2416         1.04451, 1.04033, 1.04836, 1.08690, 1.34727},
2417        {1.04261, 1.04325, 1.04628, 1.05232, 1.06448, 1.12090,
2418         1.09546, 1.09136, 1.10142, 1.13795, 1.39831},
2419        {1.09331, 1.09429, 1.09742, 1.10336, 1.11557, 1.16547,
2420         1.14658, 1.13642, 1.15247, 1.18898, 1.44936},
2421        {1.14436, 1.14539, 1.14847, 1.15443, 1.16662, 1.21794,
2422         1.19763, 1.19329, 1.20349, 1.23956, 1.50043},
2423        {1.19533, 1.19651, 1.19943, 1.20548, 1.21666, 1.26753,
2424         1.24862, 1.24434, 1.25455, 1.29106, 1.55142},
2425        {1.24638, 1.24756, 1.25046, 1.25648, 1.26764, 1.31858,
2426         1.29967, 1.29538, 1.30499, 1.34211, 1.60250},
2427        {1.29747, 1.29847, 1.30175, 1.30753, 1.31869, 1.36969,
2428         1.35069, 1.34656, 1.35663, 1.39316, 1.64644},
2429        {1.35537, 1.34952, 1.35255, 1.35869, 1.36973, 1.41387,
2430         1.40173, 1.39761, 1.40768, 1.44396, 1.70238},
2431        {1.39956, 1.40056, 1.40380, 1.40961, 1.42178, 1.46492,
2432         1.45278, 1.45423, 1.45872, 1.49522, 1.75557},
2433        {1.45080, 1.45159, 1.45463, 1.46109, 1.47287, 1.52263,
2434         1.50382, 1.50050, 1.50977, 1.54502, 1.80670},
2435        {1.50165, 1.50264, 1.50570, 1.51214, 1.52233, 1.57370,
2436         1.55484, 1.55155, 1.56080, 1.59731, 1.85778},
2437        {1.55269, 1.55364, 1.55675, 1.56274, 1.57491, 1.62598,
2438         1.60590, 1.60259, 1.61185, 1.64836, 1.90883},
2439        {1.60368, 1.60469, 1.60779, 1.61373, 1.62596, 1.67738,
2440         1.65651, 1.65249, 1.66290, 1.69936, 1.95959}};
2441
2442   // Garfield simulation at UD = -2200V; vd = 2.134cm/microsec, <driftfield> = 753V/cm
2443   Float_t time2200[ktimebin][kZbin] =
2444       {{0.09162, 0.09194, 0.09292, 0.09460, 0.09702, 0.10014,
2445         0.10421, 0.10971, 0.11820, 0.13990, 0.37745},
2446        {0.06581, 0.06618, 0.06730, 0.06915, 0.07173, 0.07507,
2447         0.07931, 0.08497, 0.09368, 0.11609, 0.35560},
2448        {0.03947, 0.04001, 0.04167, 0.04429, 0.04772, 0.05183,
2449         0.05667, 0.06293, 0.07292, 0.09900, 0.34673},
2450        {0.01111, 0.01280, 0.01716, 0.02279, 0.02876, 0.03473,
2451         0.04091, 0.04907, 0.06426, 0.10620, 0.36766},
2452        {0.01113, 0.01285, 0.01719, 0.02276, 0.02863, 0.03452,
2453         0.04076, 0.04960, 0.06745, 0.11419, 0.37139},
2454        {0.03923, 0.03978, 0.04137, 0.04387, 0.04713, 0.05110,
2455         0.05605, 0.06284, 0.07551, 0.10677, 0.36210},
2456        {0.06505, 0.06540, 0.06644, 0.06820, 0.07069, 0.07401,
2457         0.07852, 0.08479, 0.09565, 0.12325, 0.37313},
2458        {0.09107, 0.09127, 0.09181, 0.09291, 0.09474, 0.00000,
2459         0.10883, 0.11528, 0.12789, 0.15865, 0.41313},
2460        {0.12559, 0.12622, 0.12800, 0.13206, 0.14166, 0.19331,
2461         0.16832, 0.16632, 0.17724, 0.21218, 0.47098},
2462        {0.16992, 0.17070, 0.17325, 0.17831, 0.18871, 0.23557,
2463         0.21690, 0.21451, 0.22514, 0.26082, 0.52034},
2464        {0.21640, 0.21722, 0.21987, 0.22500, 0.23540, 0.28097,
2465         0.26399, 0.26154, 0.27214, 0.30784, 0.56734},
2466        {0.26318, 0.26400, 0.26679, 0.27181, 0.28220, 0.32739,
2467         0.31090, 0.30842, 0.31902, 0.35474, 0.61415},
2468        {0.31001, 0.31085, 0.31348, 0.31866, 0.32903, 0.37412,
2469         0.35777, 0.35546, 0.36588, 0.40159, 0.66103},
2470        {0.35687, 0.35769, 0.36033, 0.36556, 0.37588, 0.42094,
2471         0.40523, 0.40214, 0.41273, 0.44841, 0.70785},
2472        {0.40372, 0.40457, 0.40723, 0.41234, 0.42273, 0.46778,
2473         0.45148, 0.44903, 0.45961, 0.49526, 0.75486},
2474        {0.45062, 0.45139, 0.45404, 0.45966, 0.46958, 0.51470,
2475         0.49833, 0.49584, 0.50644, 0.54211, 0.80160},
2476        {0.49742, 0.49825, 0.50088, 0.50605, 0.51644, 0.56148,
2477         0.54518, 0.54270, 0.55330, 0.58897, 0.84854},
2478        {0.54427, 0.54510, 0.54774, 0.55290, 0.56329, 0.60846,
2479         0.59203, 0.58955, 0.60014, 0.63578, 0.89528},
2480        {0.59119, 0.59199, 0.59471, 0.59975, 0.61014, 0.65533,
2481         0.63889, 0.63636, 0.64699, 0.68269, 0.94197},
2482        {0.63866, 0.63880, 0.64145, 0.64664, 0.65701, 0.70639,
2483         0.68574, 0.68325, 0.69385, 0.72949, 0.98900},
2484        {0.68483, 0.68566, 0.68831, 0.69347, 0.70386, 0.74890,
2485         0.73260, 0.73010, 0.74069, 0.77638, 1.03320},
2486        {0.73168, 0.73255, 0.73515, 0.74031, 0.75072, 0.79576,
2487         0.77117, 0.77501, 0.78755, 0.82318, 1.08006},
2488        {0.77854, 0.78310, 0.78200, 0.79525, 0.79756, 0.84281,
2489         0.81803, 0.82393, 0.83441, 0.87008, 1.12692},
2490        {0.82541, 0.82642, 0.82916, 0.83528, 0.84442, 0.89086,
2491         0.87315, 0.87079, 0.88125, 0.91694, 1.17648},
2492        {0.87226, 0.87308, 0.87602, 0.88086, 0.89128, 0.93772,
2493         0.92001, 0.91751, 0.92811, 0.95587, 1.22328},
2494        {0.91921, 0.91994, 0.92256, 0.92772, 0.94713, 0.98566,
2495         0.96690, 0.96436, 0.97495, 1.01064, 1.26882},
2496        {0.96790, 0.96679, 0.96941, 0.97463, 0.99399, 1.03001,
2497         1.01376, 1.01112, 1.02181, 1.05749, 1.31568},
2498        {1.01278, 1.01390, 1.01674, 1.02147, 1.03182, 1.07820,
2499         1.06056, 1.05798, 1.06867, 1.10433, 1.36390},
2500        {1.05964, 1.06076, 1.06331, 1.06833, 1.07870, 1.13297,
2501         1.10742, 1.10520, 1.11543, 1.15120, 1.41069},
2502        {1.10664, 1.10762, 1.10997, 1.11519, 1.12556, 1.17531,
2503         1.15427, 1.14620, 1.16229, 1.19805, 1.45758},
2504        {1.15352, 1.15421, 1.15683, 1.16218, 1.17242, 1.21910,
2505         1.20035, 1.19863, 1.20579, 1.24473, 1.50412},
2506        {1.20019, 1.20115, 1.20369, 1.20892, 1.21928, 1.26913,
2507         1.24721, 1.24549, 1.25605, 1.29159, 1.54920},
2508        {1.24707, 1.24846, 1.25052, 1.25602, 1.26608, 1.31558,
2509         1.29448, 1.29232, 1.30293, 1.33675, 1.59798},
2510        {1.29391, 1.29475, 1.29738, 1.30255, 1.31294, 1.36244,
2511         1.34167, 1.33918, 1.34979, 1.38229, 1.64496},
2512        {1.34078, 1.34304, 1.34424, 1.35565, 1.35980, 1.40930,
2513         1.38853, 1.38229, 1.39664, 1.42863, 1.69162},
2514        {1.38762, 1.38847, 1.39110, 1.39627, 1.40666, 1.45183,
2515         1.43539, 1.43289, 1.44348, 1.47549, 1.73876},
2516        {1.43524, 1.43533, 1.43796, 1.44310, 1.45371, 1.49305,
2517         1.48224, 1.47941, 1.49034, 1.52601, 1.78552},
2518        {1.48122, 1.48219, 1.48482, 1.48991, 1.50030, 1.53991,
2519         1.52898, 1.52653, 1.53653, 1.57282, 1.82386}};
2520
2521   if (fTimeStruct1) {
2522     delete [] fTimeStruct1;
2523   }
2524   fTimeStruct1 = new Float_t[ktimebin*kZbin];
2525
2526   if (fTimeStruct2) {
2527     delete [] fTimeStruct2;
2528   }
2529   fTimeStruct2 = new Float_t[ktimebin*kZbin];
2530
2531   for (Int_t ctrt = 0; ctrt < ktimebin; ctrt++) {
2532     for (Int_t ctrz = 0; ctrz <    kZbin; ctrz++) {
2533       if      (vdrift > fVDsmp[6]) {
2534         fTimeStruct1[ctrt+ctrz*ktimebin] = time2100[ctrt][ctrz];
2535         fTimeStruct2[ctrt+ctrz*ktimebin] = time2200[ctrt][ctrz];
2536         fVDlo = fVDsmp[6];
2537         fVDhi = fVDsmp[7];
2538       }
2539       else if (vdrift > fVDsmp[5]) {
2540         fTimeStruct1[ctrt+ctrz*ktimebin] = time2000[ctrt][ctrz];
2541         fTimeStruct2[ctrt+ctrz*ktimebin] = time2100[ctrt][ctrz];
2542         fVDlo = fVDsmp[5];
2543         fVDhi = fVDsmp[6];
2544       }
2545       else if (vdrift > fVDsmp[4]) {
2546         fTimeStruct1[ctrt+ctrz*ktimebin] = time1900[ctrt][ctrz];
2547         fTimeStruct2[ctrt+ctrz*ktimebin] = time2000[ctrt][ctrz];
2548         fVDlo = fVDsmp[4];
2549         fVDhi = fVDsmp[5];
2550       }
2551       else if (vdrift > fVDsmp[3]) {
2552         fTimeStruct1[ctrt+ctrz*ktimebin] = time1800[ctrt][ctrz];
2553         fTimeStruct2[ctrt+ctrz*ktimebin] = time1900[ctrt][ctrz];
2554         fVDlo = fVDsmp[3];
2555         fVDhi = fVDsmp[4];
2556       }
2557       else if (vdrift > fVDsmp[2]) {
2558         fTimeStruct1[ctrt+ctrz*ktimebin] = time1700[ctrt][ctrz];
2559         fTimeStruct2[ctrt+ctrz*ktimebin] = time1800[ctrt][ctrz];
2560         fVDlo = fVDsmp[2];
2561         fVDhi = fVDsmp[3];
2562       }
2563       else if (vdrift > fVDsmp[1]) {
2564         fTimeStruct1[ctrt+ctrz*ktimebin] = time1600[ctrt][ctrz];
2565         fTimeStruct2[ctrt+ctrz*ktimebin] = time1700[ctrt][ctrz];
2566         fVDlo = fVDsmp[1];
2567         fVDhi = fVDsmp[2];
2568       }
2569       else if (vdrift > (fVDsmp[0] - 1.0e-5)) {
2570         fTimeStruct1[ctrt+ctrz*ktimebin] = time1500[ctrt][ctrz];
2571         fTimeStruct2[ctrt+ctrz*ktimebin] = time1600[ctrt][ctrz];
2572         fVDlo = fVDsmp[0];
2573         fVDhi = fVDsmp[1];
2574       }
2575     }
2576   }
2577
2578 }
2579
2580 //_____________________________________________________________________________
2581 void AliTRDdigitizer::RecalcDiffusion(Float_t vdrift)
2582 {
2583   //
2584   // Recalculates the diffusion parameters
2585   //
2586   // The B=0 case is not really included here.
2587   // Should be revisited!
2588   //
2589
2590   if (vdrift == fDiffLastVdrift) {
2591     return;
2592   }
2593
2594   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
2595   if (!simParam) {
2596     AliFatal("Could not get simulation parameters");
2597     return;
2598   }
2599
2600   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
2601   if (!commonParam) {
2602     AliFatal("Could not get common parameters");
2603     return;
2604   }
2605
2606   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
2607   if (!calibration) {
2608     AliFatal("Could not get calibration object");
2609     return;
2610   }
2611
2612   // The magnetic field strength
2613   Double_t x[3] = { 0.0, 0.0, 0.0 };
2614   Double_t b[3];
2615   gAlice->Field(x,b);         // b[] is in kilo Gauss
2616   Float_t field = b[2] * 0.1; // Tesla
2617
2618   fDiffLastVdrift = vdrift;
2619
2620   // DiffusionL
2621   const Int_t kNbL = 5;
2622   Float_t p0L[kNbL] = {  0.007440,  0.007493,  0.007513,  0.007672,  0.007831 };
2623   Float_t p1L[kNbL] = {  0.019252,  0.018912,  0.018636,  0.018012,  0.017343 };
2624   Float_t p2L[kNbL] = { -0.005042, -0.004926, -0.004867, -0.004650, -0.004424 };
2625   Float_t p3L[kNbL] = {  0.000195,  0.000189,  0.000195,  0.000182,  0.000169 };
2626
2627   Int_t ibL = ((Int_t) (10 * (field - 0.15)));
2628   ibL       = TMath::Max(   0,ibL);
2629   ibL       = TMath::Min(kNbL,ibL);
2630
2631   fDiffusionL = p0L[ibL]
2632               + p1L[ibL] * vdrift
2633               + p2L[ibL] * vdrift*vdrift
2634               + p3L[ibL] * vdrift*vdrift*vdrift;
2635
2636   // DiffusionT
2637   const Int_t kNbT = 5;
2638   Float_t p0T[kNbT] = {  0.009550,  0.009599,  0.009674,  0.009757,  0.009850 };
2639   Float_t p1T[kNbT] = {  0.006667,  0.006539,  0.006359,  0.006153,  0.005925 };
2640   Float_t p2T[kNbT] = { -0.000853, -0.000798, -0.000721, -0.000635, -0.000541 };
2641   Float_t p3T[kNbT] = {  0.000131,  0.000122,  0.000111,  0.000098,  0.000085 };
2642
2643   Int_t ibT= ((Int_t) (10 * (field - 0.15)));
2644   ibT      = TMath::Max(   0,ibT);
2645   ibT      = TMath::Min(kNbT,ibT);
2646
2647   fDiffusionT = p0T[ibT]
2648               + p1T[ibT] * vdrift
2649               + p2T[ibT] * vdrift*vdrift
2650               + p3T[ibT] * vdrift*vdrift*vdrift;
2651
2652   // OmegaTau
2653   fOmegaTau = calibration->GetOmegaTau(vdrift,field);
2654
2655   // Lorentzfactor
2656   if (commonParam->ExBOn()) {
2657     fLorentzFactor = 1.0 / (1.0 + fOmegaTau*fOmegaTau);
2658   }
2659   else {
2660     fLorentzFactor = 1.0;
2661   }
2662
2663 }
2664
2665 //_____________________________________________________________________________
2666 Float_t AliTRDdigitizer::GetDiffusionL(Float_t vdrift)
2667 {
2668   //
2669   // Returns the longitudinal diffusion coefficient for a given drift
2670   // velocity <vd> and a B-field <b> for Xe/CO2 (15%).
2671   // The values are according to a GARFIELD simulation.
2672   //
2673
2674   RecalcDiffusion(vdrift);
2675
2676   return fDiffusionL;
2677
2678 }
2679
2680 //_____________________________________________________________________________
2681 Float_t AliTRDdigitizer::GetDiffusionT(Float_t vdrift)
2682 {
2683   //
2684   // Returns the transverse diffusion coefficient for a given drift
2685   // velocity <vd> and a B-field <b> for Xe/CO2 (15%).
2686   // The values are according to a GARFIELD simulation.
2687   //
2688
2689   RecalcDiffusion(vdrift);
2690
2691   return fDiffusionT;
2692
2693 }
2694
2695 //_____________________________________________________________________________
2696 Int_t AliTRDdigitizer::Diffusion(Float_t vdrift, Double_t absdriftlength
2697                                , Double_t &lRow, Double_t &lCol, Double_t &lTime)
2698 {
2699   //
2700   // Applies the diffusion smearing to the position of a single electron.
2701   // Depends on absolute drift length.
2702   //
2703
2704   RecalcDiffusion(vdrift);
2705
2706   Float_t driftSqrt = TMath::Sqrt(absdriftlength);
2707   Float_t sigmaT    = driftSqrt * fDiffusionT;
2708   Float_t sigmaL    = driftSqrt * fDiffusionL;
2709   lRow  = gRandom->Gaus(lRow ,sigmaT);
2710   lCol  = gRandom->Gaus(lCol ,sigmaT * GetLorentzFactor(vdrift));
2711   lTime = gRandom->Gaus(lTime,sigmaL * GetLorentzFactor(vdrift));
2712
2713   return 1;
2714
2715 }
2716
2717 //_____________________________________________________________________________
2718 Float_t AliTRDdigitizer::GetLorentzFactor(Float_t vd)
2719 {
2720   //
2721   // Returns the recalculated Lorentz factor
2722   //
2723
2724   RecalcDiffusion(vd);
2725
2726   return fLorentzFactor;
2727
2728 }
2729
2730 //_____________________________________________________________________________
2731 Int_t AliTRDdigitizer::ExB(Float_t vdrift, Double_t driftlength, Double_t &lCol)
2732 {
2733   //
2734   // Applies E x B effects to the position of a single electron.
2735   // Depends on signed drift length.
2736   //
2737
2738   RecalcDiffusion(vdrift);
2739
2740   lCol = lCol + fOmegaTau * driftlength;
2741
2742   return 1;
2743
2744 }