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