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