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