TRD/AliTRDdigitizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Creates and handles digits from TRD hits                                 //
  21 //  Authors: C. Blume (blume@ikf.uni-frankfurt.de)                           //
  22 //           C. Lippmann                                                     //
  23 //           B. Vulpescu                                                     //
  24 //                                                                           //
  25 //  The following effects are included:                                      //
  26 //      - Diffusion                                                          //
  27 //      - ExB effects                                                        //
  28 //      - Gas gain including fluctuations                                    //
  29 //      - Pad-response (simple Gaussian approximation)                       //
  30 //      - Time-response                                                      //
  31 //      - Electronics noise                                                  //
  32 //      - Electronics gain                                                   //
  33 //      - Digitization                                                       //
  34 //  The corresponding parameter can be adjusted via the various              //
  35 //  Set-functions. If these parameters are not explicitly set, default       //
  36 //  values are used (see Init-function).                                     //
  37 //                                                                           //
  38 ///////////////////////////////////////////////////////////////////////////////
  39
  40 #include <stdlib.h>
  41
  42 #include <TMath.h>
  43 #include <TVector.h>
  44 #include <TRandom.h>
  45 #include <TROOT.h>
  46 #include <TTree.h>
  47 #include <TFile.h>
  48 #include <TF1.h>
  49 #include <TList.h>
  50 #include <TTask.h>
  51 #include <TGeoManager.h>
  52
  53 #include "AliRun.h"
  54 #include "AliRunLoader.h"
  55 #include "AliLoader.h"
  56 #include "AliConfig.h"
  57 #include "AliMagF.h"
  58 #include "AliRunDigitizer.h"
  59 #include "AliRunLoader.h"
  60 #include "AliLoader.h"
  61 #include "AliLog.h"
  62
  63 #include "AliTRD.h"
  64 #include "AliTRDhit.h"
  65 #include "AliTRDdigitizer.h"
  66 #include "AliTRDdataArrayI.h"
  67 #include "AliTRDdataArrayF.h"
  68 #include "AliTRDsegmentArray.h"
  69 #include "AliTRDdigitsManager.h"
  70 #include "AliTRDgeometry.h"
  71 #include "AliTRDpadPlane.h"
  72 #include "AliTRDcalibDB.h"
  73 #include "AliTRDSimParam.h"
  74 #include "AliTRDCommonParam.h"
  75 #include "AliTRDfeeParam.h"
  76
  77 #include "Cal/AliTRDCalROC.h"
  78 #include "Cal/AliTRDCalDet.h"
  79
  80 ClassImp(AliTRDdigitizer)
  81
  82 //_____________________________________________________________________________
  83 AliTRDdigitizer::AliTRDdigitizer()
  84   :AliDigitizer()
  85   ,fRunLoader(0)
  86   ,fDigitsManager(0)
  87   ,fSDigitsManager(0)
  88   ,fSDigitsManagerList(0)
  89   ,fTRD(0)
  90   ,fGeo(0)
  91   ,fFee(0)
  92   ,fEvent(0)
  93   ,fMasks(0)
  94   ,fCompress(kTRUE)
  95   ,fSDigits(kFALSE)
  96   ,fSDigitsScale(0.0)
  97   ,fMergeSignalOnly(kFALSE)
  98   ,fDiffLastVdrift(0)
  99   ,fDiffusionT(0.0)
 100   ,fDiffusionL(0.0)
 101   ,fOmegaTau(0.0)
 102   ,fLorentzFactor(0.0)
 103   ,fTimeLastVdrift(0)
 104   ,fTimeStruct1(0)
 105   ,fTimeStruct2(0)
 106   ,fVDlo(0)
 107   ,fVDhi(0)
 108 {
 109   //
 110   // AliTRDdigitizer default constructor
 111   //
 112
 113   Init();
 114
 115 }
 116
 117 //_____________________________________________________________________________
 118 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
 119   :AliDigitizer(name,title)
 120   ,fRunLoader(0)
 121   ,fDigitsManager(0)
 122   ,fSDigitsManager(0)
 123   ,fSDigitsManagerList(0)
 124   ,fTRD(0)
 125   ,fGeo(0)
 126   ,fFee(0)
 127   ,fEvent(0)
 128   ,fMasks(0)
 129   ,fCompress(kTRUE)
 130   ,fSDigits(kFALSE)
 131   ,fSDigitsScale(0.0)
 132   ,fMergeSignalOnly(kFALSE)
 133   ,fDiffLastVdrift(0)
 134   ,fDiffusionT(0.0)
 135   ,fDiffusionL(0.0)
 136   ,fOmegaTau(0.0)
 137   ,fLorentzFactor(0.0)
 138   ,fTimeLastVdrift(0)
 139   ,fTimeStruct1(0)
 140   ,fTimeStruct2(0)
 141   ,fVDlo(0)
 142   ,fVDhi(0)
 143 {
 144   //
 145   // AliTRDdigitizer constructor
 146   //
 147
 148   Init();
 149
 150 }
 151
 152 //_____________________________________________________________________________
 153 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager
 154                                , const Text_t *name, const Text_t *title)
 155   :AliDigitizer(manager,name,title)
 156   ,fRunLoader(0)
 157   ,fDigitsManager(0)
 158   ,fSDigitsManager(0)
 159   ,fSDigitsManagerList(0)
 160   ,fTRD(0)
 161   ,fGeo(0)
 162   ,fFee(0)
 163   ,fEvent(0)
 164   ,fMasks(0)
 165   ,fCompress(kTRUE)
 166   ,fSDigits(kFALSE)
 167   ,fSDigitsScale(0.0)
 168   ,fMergeSignalOnly(kFALSE)
 169   ,fDiffLastVdrift(0)
 170   ,fDiffusionT(0.0)
 171   ,fDiffusionL(0.0)
 172   ,fOmegaTau(0.0)
 173   ,fLorentzFactor(0.0)
 174   ,fTimeLastVdrift(0)
 175   ,fTimeStruct1(0)
 176   ,fTimeStruct2(0)
 177   ,fVDlo(0)
 178   ,fVDhi(0)
 179 {
 180   //
 181   // AliTRDdigitizer constructor
 182   //
 183
 184   Init();
 185
 186 }
 187
 188 //_____________________________________________________________________________
 189 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager)
 190   :AliDigitizer(manager,"AliTRDdigitizer","TRD digitizer")
 191   ,fRunLoader(0)
 192   ,fDigitsManager(0)
 193   ,fSDigitsManager(0)
 194   ,fSDigitsManagerList(0)
 195   ,fTRD(0)
 196   ,fGeo(0)
 197   ,fFee(0)
 198   ,fEvent(0)
 199   ,fMasks(0)
 200   ,fCompress(kTRUE)
 201   ,fSDigits(kFALSE)
 202   ,fSDigitsScale(0.0)
 203   ,fMergeSignalOnly(kFALSE)
 204   ,fDiffLastVdrift(0)
 205   ,fDiffusionT(0.0)
 206   ,fDiffusionL(0.0)
 207   ,fOmegaTau(0.0)
 208   ,fLorentzFactor(0.0)
 209   ,fTimeLastVdrift(0)
 210   ,fTimeStruct1(0)
 211   ,fTimeStruct2(0)
 212   ,fVDlo(0)
 213   ,fVDhi(0)
 214 {
 215   //
 216   // AliTRDdigitizer constructor
 217   //
 218
 219   Init();
 220
 221 }
 222
 223 //_____________________________________________________________________________
 224 Bool_t AliTRDdigitizer::Init()
 225 {
 226   //
 227   // Initialize the digitizer with default values
 228   //
 229
 230   fRunLoader          = 0;
 231   fDigitsManager      = 0;
 232   fSDigitsManager     = 0;
 233   fSDigitsManagerList = 0;
 234   fTRD                = 0;
 235   fGeo                = 0;
 236   fFee                = AliTRDfeeParam::Instance();
 237
 238   fEvent              = 0;
 239   fMasks              = 0;
 240   fCompress           = kTRUE;
 241   fSDigits            = kFALSE;
 242   fSDigitsScale       = 100.0;
 243   fMergeSignalOnly    = kFALSE;
 244
 245   fTimeLastVdrift     = -1;
 246   fTimeStruct1        =  0;
 247   fTimeStruct2        =  0;
 248   fVDlo               =  0;
 249   fVDhi               =  0;
 250
 251   fDiffLastVdrift     = -1;
 252   fDiffusionT         =  0.0;
 253   fDiffusionL         =  0.0;
 254   fLorentzFactor      =  0.0;
 255
 256   return AliDigitizer::Init();
 257
 258 }
 259
 260 //_____________________________________________________________________________
 261 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
 262   :AliDigitizer(d)
 263   ,fRunLoader(0)
 264   ,fDigitsManager(0)
 265   ,fSDigitsManager(0)
 266   ,fSDigitsManagerList(0)
 267   ,fTRD(0)
 268   ,fGeo(0)
 269   ,fFee(0)
 270   ,fEvent(0)
 271   ,fMasks(0)
 272   ,fCompress(d.fCompress)
 273   ,fSDigits(d.fSDigits)
 274   ,fSDigitsScale(d.fSDigitsScale)
 275   ,fMergeSignalOnly(d.fMergeSignalOnly)
 276   ,fDiffLastVdrift(-1)
 277   ,fDiffusionT(0.0)
 278   ,fDiffusionL(0.0)
 279   ,fOmegaTau(0.0)
 280   ,fLorentzFactor(0.0)
 281   ,fTimeLastVdrift(-1)
 282   ,fTimeStruct1(0)
 283   ,fTimeStruct2(0)
 284   ,fVDlo(0)
 285   ,fVDhi(0)
 286 {
 287   //
 288   // AliTRDdigitizer copy constructor
 289   //
 290
 291   // Do not copy timestructs, just invalidate lastvdrift.
 292   // Next time they are requested, they get recalculated
 293   if (((AliTRDdigitizer &) d).fTimeStruct1) {
 294     delete [] ((AliTRDdigitizer &) d).fTimeStruct1;
 295     ((AliTRDdigitizer &) d).fTimeStruct1 = 0;
 296   }
 297   if (((AliTRDdigitizer &) d).fTimeStruct2) {
 298     delete [] ((AliTRDdigitizer &) d).fTimeStruct2;
 299     ((AliTRDdigitizer &) d).fTimeStruct2 = 0;
 300   }
 301
 302 }
 303
 304 //_____________________________________________________________________________
 305 AliTRDdigitizer::~AliTRDdigitizer()
 306 {
 307   //
 308   // AliTRDdigitizer destructor
 309   //
 310
 311   if (fDigitsManager) {
 312     delete fDigitsManager;
 313     fDigitsManager      = 0;
 314   }
 315
 316   if (fDigitsManager) {
 317     delete fSDigitsManager;
 318     fSDigitsManager     = 0;
 319   }
 320
 321   if (fSDigitsManagerList) {
 322     fSDigitsManagerList->Delete();
 323     delete fSDigitsManagerList;
 324     fSDigitsManagerList = 0;
 325   }
 326
 327   if (fMasks) {
 328     delete [] fMasks;
 329     fMasks       = 0;
 330   }
 331
 332   if (fTimeStruct1) {
 333     delete [] fTimeStruct1;
 334     fTimeStruct1 = 0;
 335   }
 336
 337   if (fTimeStruct2) {
 338     delete [] fTimeStruct2;
 339     fTimeStruct2 = 0;
 340   }
 341
 342   if (fGeo) {
 343     delete fGeo;
 344     fGeo = 0;
 345   }
 346
 347 }
 348
 349 //_____________________________________________________________________________
 350 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
 351 {
 352   //
 353   // Assignment operator
 354   //
 355
 356   if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
 357
 358   return *this;
 359
 360 }
 361
 362 //_____________________________________________________________________________
 363 void AliTRDdigitizer::Copy(TObject &d) const
 364 {
 365   //
 366   // Copy function
 367   //
 368
 369   ((AliTRDdigitizer &) d).fRunLoader          = 0;
 370   ((AliTRDdigitizer &) d).fDigitsManager      = 0;
 371   ((AliTRDdigitizer &) d).fSDigitsManager     = 0;
 372   ((AliTRDdigitizer &) d).fSDigitsManagerList = 0;
 373   ((AliTRDdigitizer &) d).fTRD                = 0;
 374   ((AliTRDdigitizer &) d).fGeo                = 0;
 375   ((AliTRDdigitizer &) d).fFee                = fFee;
 376   ((AliTRDdigitizer &) d).fEvent              = 0;
 377   ((AliTRDdigitizer &) d).fMasks              = 0;
 378   ((AliTRDdigitizer &) d).fCompress           = fCompress;
 379   ((AliTRDdigitizer &) d).fSDigits            = fSDigits;
 380   ((AliTRDdigitizer &) d).fSDigitsScale       = fSDigitsScale;
 381   ((AliTRDdigitizer &) d).fMergeSignalOnly    = fMergeSignalOnly;
 382
 383   AliTRDdigitizer& target = (AliTRDdigitizer &) d;
 384
 385   // Do not copy timestructs, just invalidate lastvdrift.
 386   // Next time they are requested, they get recalculated
 387   if (target.fTimeStruct1) {
 388     delete [] target.fTimeStruct1;
 389     target.fTimeStruct1 = 0;
 390   }
 391   if (target.fTimeStruct2) {
 392     delete [] target.fTimeStruct2;
 393     target.fTimeStruct2 = 0;
 394   }
 395   target.fTimeLastVdrift = -1;
 396
 397 }
 398
 399 //_____________________________________________________________________________
 400 void AliTRDdigitizer::Exec(Option_t *option)
 401 {
 402   //
 403   // Executes the merging
 404   //
 405
 406   Int_t iInput;
 407
 408   AliTRDdigitsManager *sdigitsManager;
 409
 410   TString optionString = option;
 411   if (optionString.Contains("deb")) {
 412     AliLog::SetClassDebugLevel("AliTRDdigitizer",1);
 413     AliInfo("Called with debug option");
 414   }
 415
 416   // The AliRoot file is already connected by the manager
 417   AliRunLoader *inrl;
 418
 419   if (gAlice) {
 420     AliDebug(1,"AliRun object found on file.");
 421   }
 422   else {
 423     inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
 424     inrl->LoadgAlice();
 425     gAlice = inrl->GetAliRun();
 426     if (!gAlice) {
 427       AliError("Could not find AliRun object.")
 428       return;
 429     }
 430   }
 431
 432   Int_t nInput = fManager->GetNinputs();
 433   fMasks       = new Int_t[nInput];
 434   for (iInput = 0; iInput < nInput; iInput++) {
 435     fMasks[iInput] = fManager->GetMask(iInput);
 436   }
 437
 438   //
 439   // Initialization
 440   //
 441
 442   AliRunLoader *orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
 443
 444   if (InitDetector()) {
 445
 446     AliLoader *ogime = orl->GetLoader("TRDLoader");
 447
 448     TTree *tree = 0;
 449     if (fSDigits) {
 450       // If we produce SDigits
 451       tree = ogime->TreeS();
 452       if (!tree) {
 453         ogime->MakeTree("S");
 454         tree = ogime->TreeS();
 455       }
 456     }
 457     else {
 458       // If we produce Digits
 459       tree = ogime->TreeD();
 460       if (!tree) {
 461         ogime->MakeTree("D");
 462         tree = ogime->TreeD();
 463       }
 464     }
 465
 466     MakeBranch(tree);
 467
 468   }
 469
 470   for (iInput = 0; iInput < nInput; iInput++) {
 471
 472     AliDebug(1,Form("Add input stream %d",iInput));
 473
 474     // Check if the input tree exists
 475     inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
 476     AliLoader *gime = inrl->GetLoader("TRDLoader");
 477
 478     TTree *treees = gime->TreeS();
 479     if (treees == 0x0) {
 480       if (gime->LoadSDigits()) {
 481         AliError(Form("Error Occured while loading S. Digits for input %d.",iInput));
 482         return;
 483       }
 484       treees = gime->TreeS();
 485     }
 486
 487     if (treees == 0x0) {
 488       AliError(Form("Input stream %d does not exist",iInput));
 489       return;
 490     }
 491
 492     // Read the s-digits via digits manager
 493     sdigitsManager = new AliTRDdigitsManager();
 494     sdigitsManager->SetSDigits(kTRUE);
 495
 496     AliRunLoader *rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
 497     AliLoader *gimme = rl->GetLoader("TRDLoader");
 498     if (!gimme->TreeS()) {
 499       gimme->LoadSDigits();
 500     }
 501     sdigitsManager->ReadDigits(gimme->TreeS());
 502
 503     // Add the s-digits to the input list
 504     AddSDigitsManager(sdigitsManager);
 505
 506   }
 507
 508   // Convert the s-digits to normal digits
 509   AliDebug(1,"Do the conversion");
 510   SDigits2Digits();
 511
 512   // Store the digits
 513   AliDebug(1,"Write the digits");
 514   fDigitsManager->WriteDigits();
 515
 516   // Write parameters
 517   orl->CdGAFile();
 518
 519   AliDebug(1,"Done");
 520
 521   DeleteSDigitsManager();
 522
 523 }
 524
 525 //_____________________________________________________________________________
 526 Bool_t AliTRDdigitizer::Open(const Char_t *file, Int_t nEvent)
 527 {
 528   //
 529   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
 530   //
 531   // Connect the AliRoot file containing Geometry, Kine, and Hits
 532   //
 533
 534   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 535
 536   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
 537   if (!fRunLoader) {
 538     fRunLoader = AliRunLoader::Open(file,evfoldname,"UPDATE");
 539   }
 540   if (!fRunLoader) {
 541     AliError(Form("Can not open session for file %s.",file));
 542     return kFALSE;
 543   }
 544
 545   if (!fRunLoader->GetAliRun()) {
 546     fRunLoader->LoadgAlice();
 547   }
 548   gAlice = fRunLoader->GetAliRun();
 549
 550   if (gAlice) {
 551     AliDebug(1,"AliRun object found on file.");
 552   }
 553   else {
 554     AliError("Could not find AliRun object.");
 555     return kFALSE;
 556   }
 557
 558   fEvent = nEvent;
 559
 560   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
 561   if (!loader) {
 562     AliError("Can not get TRD loader from Run Loader");
 563     return kFALSE;
 564   }
 565
 566   if (InitDetector()) {
 567     TTree *tree = 0;
 568     if (fSDigits) {
 569       // If we produce SDigits
 570       tree = loader->TreeS();
 571       if (!tree) {
 572         loader->MakeTree("S");
 573         tree = loader->TreeS();
 574       }
 575     }
 576     else {
 577       // If we produce Digits
 578       tree = loader->TreeD();
 579       if (!tree) {
 580         loader->MakeTree("D");
 581         tree = loader->TreeD();
 582       }
 583     }
 584     return MakeBranch(tree);
 585   }
 586   else {
 587     return kFALSE;
 588   }
 589
 590 }
 591
 592 //_____________________________________________________________________________
 593 Bool_t AliTRDdigitizer::Open(AliRunLoader *runLoader, Int_t nEvent)
 594 {
 595   //
 596   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
 597   //
 598   // Connect the AliRoot file containing Geometry, Kine, and Hits
 599   //
 600
 601   fRunLoader = runLoader;
 602   if (!fRunLoader) {
 603     AliError("RunLoader does not exist");
 604     return kFALSE;
 605   }
 606
 607   if (!fRunLoader->GetAliRun()) {
 608     fRunLoader->LoadgAlice();
 609   }
 610   gAlice = fRunLoader->GetAliRun();
 611
 612   if (gAlice) {
 613     AliDebug(1,"AliRun object found on file.");
 614   }
 615   else {
 616     AliError("Could not find AliRun object.");
 617     return kFALSE;
 618   }
 619
 620   fEvent = nEvent;
 621
 622   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
 623   if (!loader) {
 624     AliError("Can not get TRD loader from Run Loader");
 625     return kFALSE;
 626   }
 627
 628   if (InitDetector()) {
 629     TTree *tree = 0;
 630     if (fSDigits) {
 631       // If we produce SDigits
 632       tree = loader->TreeS();
 633       if (!tree) {
 634         loader->MakeTree("S");
 635         tree = loader->TreeS();
 636       }
 637     }
 638     else {
 639       // If we produce Digits
 640       tree = loader->TreeD();
 641       if (!tree) {
 642         loader->MakeTree("D");
 643         tree = loader->TreeD();
 644       }
 645     }
 646     return MakeBranch(tree);
 647   }
 648   else {
 649     return kFALSE;
 650   }
 651
 652 }
 653
 654 //_____________________________________________________________________________
 655 Bool_t AliTRDdigitizer::InitDetector()
 656 {
 657   //
 658   // Sets the pointer to the TRD detector and the geometry
 659   //
 660
 661   // Get the pointer to the detector class and check for version 1
 662   fTRD = (AliTRD *) gAlice->GetDetector("TRD");
 663   if (!fTRD) {
 664     AliFatal("No TRD module found");
 665     exit(1);
 666   }
 667   if (fTRD->IsVersion() != 1) {
 668     AliFatal("TRD must be version 1 (slow simulator)");
 669     exit(1);
 670   }
 671
 672   // Get the geometry
 673   fGeo = new AliTRDgeometry();
 674
 675   // Create a digits manager
 676   delete fDigitsManager;
 677   fDigitsManager = new AliTRDdigitsManager();
 678   fDigitsManager->SetSDigits(fSDigits);
 679   fDigitsManager->CreateArrays();
 680   fDigitsManager->SetEvent(fEvent);
 681
 682   // The list for the input s-digits manager to be merged
 683   if (fSDigitsManagerList) {
 684     fSDigitsManagerList->Delete();
 685   }
 686   else {
 687     fSDigitsManagerList = new TList();
 688   }
 689
 690   return kTRUE;
 691
 692 }
 693
 694 //_____________________________________________________________________________
 695 Bool_t AliTRDdigitizer::MakeBranch(TTree *tree) const
 696 {
 697   //
 698   // Create the branches for the digits array
 699   //
 700
 701   return fDigitsManager->MakeBranch(tree);
 702
 703 }
 704
 705 //_____________________________________________________________________________
 706 Bool_t AliTRDdigitizer::MakeDigits()
 707 {
 708   //
 709   // Creates digits.
 710   //
 711
 712   ///////////////////////////////////////////////////////////////
 713   // Parameter
 714   ///////////////////////////////////////////////////////////////
 715
 716   // Converts number of electrons to fC
 717   const Double_t kEl2fC  = 1.602e-19 * 1.0e15;
 718
 719   ///////////////////////////////////////////////////////////////
 720
 721   // Number of pads included in the pad response
 722   const Int_t kNpad      = 3;
 723
 724   // Number of track dictionary arrays
 725   const Int_t kNDict     = AliTRDdigitsManager::kNDict;
 726
 727   // Width of the amplification region
 728   const Float_t kAmWidth = AliTRDgeometry::AmThick();
 729   // Width of the drift region
 730   const Float_t kDrWidth = AliTRDgeometry::DrThick();
 731   // Drift + amplification region
 732   const Float_t kDrMin   =          - 0.5 * kAmWidth;
 733   const Float_t kDrMax   = kDrWidth + 0.5 * kAmWidth;
 734
 735   Int_t    iRow;
 736   Int_t    iCol;
 737   Int_t    iTime;
 738   Int_t    iPad;
 739   Int_t    iDict  = 0;
 740   Int_t    nBytes = 0;
 741
 742   Int_t    totalSizeDigits = 0;
 743   Int_t    totalSizeDict0  = 0;
 744   Int_t    totalSizeDict1  = 0;
 745   Int_t    totalSizeDict2  = 0;
 746
 747   Int_t    timeBinTRFend   = 1;
 748
 749   Double_t pos[3];
 750   Double_t loc[3];
 751   Double_t padSignal[kNpad];
 752   Double_t signalOld[kNpad];
 753
 754   AliTRDdataArrayF *signals  = 0;
 755   AliTRDdataArrayI *digits   = 0;
 756   AliTRDdataArrayI *dictionary[kNDict];
 757
 758   AliTRDpadPlane   *padPlane = 0;
 759
 760   AliTRDCalROC     *calVdriftROC          = 0;
 761   Float_t           calVdriftDetValue     = 0.0;
 762   AliTRDCalROC     *calT0ROC              = 0;
 763   Float_t           calT0DetValue         = 0.0;
 764   AliTRDCalROC     *calGainFactorROC      = 0;
 765   Float_t           calGainFactorDetValue = 0.0;
 766
 767   if (!gGeoManager) {
 768     AliFatal("No geometry manager!");
 769   }
 770
 771   if (!fGeo) {
 772     AliError("No geometry defined");
 773     return kFALSE;
 774   }
 775
 776   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 777   if (!simParam) {
 778     AliFatal("Could not get simulation parameters");
 779     return kFALSE;
 780   }
 781
 782   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
 783   if (!commonParam) {
 784     AliFatal("Could not get common parameterss");
 785     return kFALSE;
 786   }
 787
 788   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
 789   if (!calibration) {
 790     AliFatal("Could not get calibration object");
 791     return kFALSE;
 792   }
 793
 794   AliDebug(1,"Start creating digits");
 795
 796   // Create a container for the amplitudes
 797   AliTRDsegmentArray *signalsArray = new AliTRDsegmentArray("AliTRDdataArrayF"
 798                                                            ,AliTRDgeometry::Ndet());
 799
 800   AliLoader *gimme = fRunLoader->GetLoader("TRDLoader");
 801   if (!gimme->TreeH()) {
 802     gimme->LoadHits();
 803   }
 804   TTree *hitTree = gimme->TreeH();
 805   if (hitTree == 0x0) {
 806     AliError("Can not get TreeH");
 807     return kFALSE;
 808   }
 809   fTRD->SetTreeAddress();
 810
 811   // Get the number of entries in the hit tree
 812   // (Number of primary particles creating a hit somewhere)
 813   Int_t nTrack = (Int_t) hitTree->GetEntries();
 814   AliDebug(1,Form("Found %d primary particles",nTrack));
 815   AliDebug(1,Form("Sampling = %.0fMHz"        ,commonParam->GetSamplingFrequency()));
 816   AliDebug(1,Form("Gain     = %d"             ,((Int_t) simParam->GetGasGain())));
 817   AliDebug(1,Form("Noise    = %d"             ,((Int_t) simParam->GetNoise())));
 818   if (simParam->TimeStructOn()) {
 819     AliDebug(1,"Time Structure of drift cells implemented.");
 820   }
 821   else {
 822     AliDebug(1,"Constant drift velocity in drift cells.");
 823   }
 824   if (simParam->TRFOn()) {
 825     timeBinTRFend = ((Int_t) (simParam->GetTRFhi()
 826                   * commonParam->GetSamplingFrequency())) - 1;
 827     AliDebug(1,Form("Sample the TRF up to bin %d",timeBinTRFend));
 828   }
 829
 830   // Get the detector wise calibration objects
 831   const AliTRDCalDet *calVdriftDet     = calibration->GetVdriftDet();
 832   const AliTRDCalDet *calT0Det         = calibration->GetT0Det();
 833   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 834
 835   Int_t   detectorOld  = -1;
 836   Int_t   countHits    =  0;
 837
 838   Int_t   nTimeTotal   = calibration->GetNumberOfTimeBins();
 839   Float_t samplingRate = commonParam->GetSamplingFrequency();
 840   Float_t elAttachProp = simParam->GetElAttachProp() / 100.0;
 841
 842   // Loop through all entries in the tree
 843   for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
 844
 845     gAlice->ResetHits();
 846     nBytes += hitTree->GetEvent(iTrack);
 847
 848     // Loop through the TRD hits
 849     Int_t iHit = 0;
 850     AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
 851     while (hit) {
 852
 853       countHits++;
 854       iHit++;
 855
 856       Float_t q        = hit->GetCharge();
 857       // Don't analyze test hits
 858       if (((Int_t) q) == 0) {
 859         hit = (AliTRDhit *) fTRD->NextHit();
 860         continue;
 861       }
 862
 863       pos[0]           = hit->X();
 864       pos[1]           = hit->Y();
 865       pos[2]           = hit->Z();
 866       Int_t   track    = hit->Track();
 867       Int_t   detector = hit->GetDetector();
 868       Float_t hittime  = hit->GetTime();
 869       Int_t   plane    = fGeo->GetPlane(detector);
 870       Int_t   chamber  = fGeo->GetChamber(detector);
 871       padPlane         = fGeo->GetPadPlane(plane,chamber);
 872       Float_t row0     = padPlane->GetRow0ROC();
 873       Int_t   nRowMax  = padPlane->GetNrows();
 874       Int_t   nColMax  = padPlane->GetNcols();
 875       Int_t   inDrift  = 1;
 876
 877       // Find the current volume with the geo manager
 878       gGeoManager->SetCurrentPoint(pos);
 879       gGeoManager->FindNode();
 880       if (strstr(gGeoManager->GetPath(),"/UK")) {
 881         inDrift = 0;
 882       }
 883
 884       if (detector != detectorOld) {
 885
 886         // Compress the old one if enabled
 887         if ((fCompress) && (detectorOld > -1)) {
 888           signals->Compress(1,0);
 889           for (iDict = 0; iDict < kNDict; iDict++) {
 890             dictionary[iDict]->Compress(1,0);
 891           }
 892         }
 893         // Get the new container
 894         signals = (AliTRDdataArrayF *) signalsArray->At(detector);
 895         if (signals->GetNtime() == 0) {
 896           // Allocate a new one if not yet existing
 897           signals->Allocate(nRowMax,nColMax,nTimeTotal);
 898         }
 899         else {
 900           // Expand an existing one
 901           if (fCompress) {
 902             signals->Expand();
 903           }
 904         }
 905         // The same for the dictionary
 906         for (iDict = 0; iDict < kNDict; iDict++) {
 907           dictionary[iDict] = fDigitsManager->GetDictionary(detector,iDict);
 908           if (dictionary[iDict]->GetNtime() == 0) {
 909             dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
 910           }
 911           else {
 912             if (fCompress) dictionary[iDict]->Expand();
 913           }
 914         }
 915
 916         // Get the calibration objects
 917         calVdriftROC      = calibration->GetVdriftROC(detector);
 918         calVdriftDetValue = calVdriftDet->GetValue(detector);
 919         calT0ROC          = calibration->GetT0ROC(detector);
 920         calT0DetValue     = calT0Det->GetValue(detector);
 921
 922         detectorOld = detector;
 923
 924       }
 925
 926       // Go to the local coordinate system:
 927       // loc[0] - col  direction in amplification or driftvolume
 928       // loc[1] - row  direction in amplification or driftvolume
 929       // loc[2] - time direction in amplification or driftvolume
 930       gGeoManager->MasterToLocal(pos,loc);
 931       if (inDrift) {
 932         // Relative to middle of amplification region
 933         loc[2] = loc[2] - kDrWidth/2.0 - kAmWidth/2.0;
 934       }
 935
 936       // The driftlength [cm] (w/o diffusion yet !).
 937       // It is negative if the hit is between pad plane and anode wires.
 938       Double_t driftlength = -1.0 * loc[2];
 939
 940       // Stupid patch to take care of TR photons that are absorbed
 941       // outside the chamber volume. A real fix would actually need
 942       // a more clever implementation of the TR hit generation
 943       if (q < 0.0) {
 944         if ((loc[1] < padPlane->GetRowEndROC()) ||
 945             (loc[1] > padPlane->GetRow0ROC())) {
 946           hit = (AliTRDhit *) fTRD->NextHit();
 947           continue;
 948         }
 949         if ((driftlength < kDrMin) ||
 950             (driftlength > kDrMax)) {
 951           hit = (AliTRDhit *) fTRD->NextHit();
 952           continue;
 953         }
 954       }
 955
 956       // Get row and col of unsmeared electron to retrieve drift velocity
 957       // The pad row (z-direction)
 958       Int_t    rowE         = padPlane->GetPadRowNumberROC(loc[1]);
 959       if (rowE < 0) {
 960         hit = (AliTRDhit *) fTRD->NextHit();
 961         continue;
 962       }
 963       Double_t rowOffset    = padPlane->GetPadRowOffsetROC(rowE,loc[1]);
 964
 965       // The pad column (rphi-direction)
 966       Double_t offsetTilt   = padPlane->GetTiltOffset(rowOffset);
 967       Int_t    colE         = padPlane->GetPadColNumber(loc[0]+offsetTilt);
 968       if (colE < 0) {
 969         hit = (AliTRDhit *) fTRD->NextHit();
 970         continue;
 971       }
 972       Double_t colOffset    = padPlane->GetPadColOffset(colE,loc[0]+offsetTilt);
 973
 974       Float_t driftvelocity = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
 975
 976       // Normalized drift length
 977       Double_t absdriftlength = TMath::Abs(driftlength);
 978       if (commonParam->ExBOn()) {
 979         absdriftlength /= TMath::Sqrt(GetLorentzFactor(driftvelocity));
 980       }
 981
 982       // Loop over all electrons of this hit
 983       // TR photons produce hits with negative charge
 984       Int_t nEl = ((Int_t) TMath::Abs(q));
 985       for (Int_t iEl = 0; iEl < nEl; iEl++) {
 986
 987         // Now the real local coordinate system of the ROC
 988         // column direction: locC
 989         // row direction:    locR
 990         // time direction:   locT
 991         // locR and locC are identical to the coordinates of the corresponding
 992         // volumina of the drift or amplification region.
 993         // locT is defined relative to the wire plane (i.e. middle of amplification
 994         // region), meaming locT = 0, and is negative for hits coming from the
 995         // drift region.
 996         Double_t locC = loc[0];
 997         Double_t locR = loc[1];
 998         Double_t locT = loc[2];
 999
1000         // Electron attachment
1001         if (simParam->ElAttachOn()) {
1002           if (gRandom->Rndm() < (absdriftlength * elAttachProp)) {
1003             continue;
1004           }
1005         }
1006
1007         // Apply the diffusion smearing
1008         if (simParam->DiffusionOn()) {
1009           if (!(Diffusion(driftvelocity,absdriftlength,locR,locC,locT))) {
1010             continue;
1011           }
1012         }
1013
1014         // Apply E x B effects (depends on drift direction)
1015         if (commonParam->ExBOn()) {
1016           if (!(ExB(driftvelocity,driftlength,locC))) {
1017             continue;
1018           }
1019         }
1020
1021         // The electron position after diffusion and ExB in pad coordinates.
1022         // The pad row (z-direction)
1023         rowE       = padPlane->GetPadRowNumberROC(locR);
1024         if (rowE < 0) continue;
1025         rowOffset  = padPlane->GetPadRowOffsetROC(rowE,locR);
1026
1027         // The pad column (rphi-direction)
1028         offsetTilt = padPlane->GetTiltOffset(rowOffset);
1029         colE       = padPlane->GetPadColNumber(locC+offsetTilt);
1030         if (colE < 0) continue;
1031         colOffset  = padPlane->GetPadColOffset(colE,locC+offsetTilt);
1032
1033         // Also re-retrieve drift velocity because col and row may have changed
1034         driftvelocity = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
1035         Float_t t0    = calT0DetValue     + calT0ROC->GetValue(colE,rowE);
1036
1037         // Convert the position to drift time [mus], using either constant drift velocity or
1038         // time structure of drift cells (non-isochronity, GARFIELD calculation).
1039         // Also add absolute time of hits to take pile-up events into account properly
1040         Double_t drifttime;
1041         if (simParam->TimeStructOn()) {
1042           // Get z-position with respect to anode wire
1043           Double_t zz  =  row0 - locR + simParam->GetAnodeWireOffset();
1044           zz -= ((Int_t)(2 * zz)) / 2.0;
1045           if (zz > 0.25) {
1046             zz  = 0.5 - zz;
1047           }
1048           // Use drift time map (GARFIELD)
1049           drifttime = TimeStruct(driftvelocity,0.5*kAmWidth-1.0*locT,zz)
1050                     + hittime;
1051         }
1052         else {
1053           // Use constant drift velocity
1054           drifttime = -1.0 * locT / driftvelocity
1055                     + hittime;
1056         }
1057
1058         // Apply the gas gain including fluctuations
1059         Double_t ggRndm = 0.0;
1060         do {
1061           ggRndm = gRandom->Rndm();
1062         } while (ggRndm <= 0);
1063         Int_t signal = (Int_t) (-(simParam->GetGasGain()) * TMath::Log(ggRndm));
1064
1065         // Apply the pad response
1066         if (simParam->PRFOn()) {
1067           // The distance of the electron to the center of the pad
1068           // in units of pad width
1069           Double_t dist = (colOffset - 0.5*padPlane->GetColSize(colE))
1070                         / padPlane->GetColSize(colE);
1071           // This is still the fixed parametrization, i.e. not dependent on
1072           // calibration values !!!!
1073           if (!(calibration->PadResponse(signal,dist,plane,padSignal))) continue;
1074         }
1075         else {
1076           padSignal[0] = 0.0;
1077           padSignal[1] = signal;
1078           padSignal[2] = 0.0;
1079         }
1080
1081         // The time bin (always positive), with t0 distortion
1082         Double_t timeBinIdeal = drifttime * samplingRate + t0;
1083         // Protection
1084         if (TMath::Abs(timeBinIdeal) > 2*nTimeTotal) {
1085           timeBinIdeal = 2 * nTimeTotal;
1086         }
1087         Int_t    timeBinTruncated = (Int_t) timeBinIdeal;
1088         // The distance of the position to the middle of the timebin
1089         Double_t timeOffset       = ((Float_t) timeBinTruncated
1090                                   + 0.5 - timeBinIdeal) / samplingRate;
1091
1092         // Sample the time response inside the drift region
1093         // + additional time bins before and after.
1094         // The sampling is done always in the middle of the time bin
1095         for (Int_t iTimeBin = TMath::Max(timeBinTruncated,0)
1096             ;iTimeBin < TMath::Min(timeBinTruncated+timeBinTRFend,nTimeTotal)
1097             ;iTimeBin++) {
1098
1099           // Apply the time response
1100           Double_t timeResponse = 1.0;
1101           Double_t crossTalk    = 0.0;
1102           Double_t time         = (iTimeBin - timeBinTruncated) / samplingRate + timeOffset;
1103           if (simParam->TRFOn()) {
1104             timeResponse = simParam->TimeResponse(time);
1105           }
1106           if (simParam->CTOn()) {
1107             crossTalk    = simParam->CrossTalk(time);
1108           }
1109
1110           signalOld[0] = 0.0;
1111           signalOld[1] = 0.0;
1112           signalOld[2] = 0.0;
1113
1114           for (iPad = 0; iPad < kNpad; iPad++) {
1115
1116             Int_t colPos = colE + iPad - 1;
1117             if (colPos <        0) continue;
1118             if (colPos >= nColMax) break;
1119
1120             // Add the signals
1121             Int_t iCurrentTimeBin = iTimeBin;
1122             signalOld[iPad]  = signals->GetDataUnchecked(rowE,colPos,iCurrentTimeBin);
1123             if (colPos != colE) {
1124               signalOld[iPad] += padSignal[iPad] * (timeResponse + crossTalk);
1125             }
1126             else {
1127               signalOld[iPad] += padSignal[iPad] * timeResponse;
1128             }
1129             signals->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,signalOld[iPad]);
1130
1131             // Store the track index in the dictionary
1132             // Note: We store index+1 in order to allow the array to be compressed
1133             if (signalOld[iPad] > 0) {
1134               for (iDict = 0; iDict < kNDict; iDict++) {
1135                 Int_t oldTrack = dictionary[iDict]->GetDataUnchecked(rowE
1136                                                                     ,colPos
1137                                                                     ,iCurrentTimeBin);
1138                 if (oldTrack == track+1) break;
1139                 if (oldTrack ==       0) {
1140                   dictionary[iDict]->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,track+1);
1141                   break;
1142                 }
1143               }
1144             }
1145
1146           } // Loop: pads
1147
1148         } // Loop: time bins
1149
1150       } // Loop: electrons of a single hit
1151
1152       hit = (AliTRDhit *) fTRD->NextHit();
1153
1154     } // Loop: hits of one primary track
1155
1156   } // Loop: primary tracks
1157
1158   AliDebug(1,Form("Finished analyzing %d hits",countHits));
1159
1160   //____________________________________________________________________
1161   //
1162   // Create (s)digits
1163   //____________________________________________________________________
1164   //
1165
1166   // The coupling factor
1167   Double_t coupling = simParam->GetPadCoupling()
1168                     * simParam->GetTimeCoupling();
1169
1170   // The conversion factor
1171   Double_t convert  = kEl2fC
1172                     * simParam->GetChipGain();
1173
1174   // Loop through all chambers to finalize the digits
1175   Int_t iDetBeg = 0;
1176   Int_t iDetEnd = AliTRDgeometry::Ndet();
1177   for (Int_t iDet = iDetBeg; iDet < iDetEnd; iDet++) {
1178
1179     Int_t plane      = fGeo->GetPlane(iDet);
1180     Int_t sector     = fGeo->GetSector(iDet);
1181     Int_t chamber    = fGeo->GetChamber(iDet);
1182     Int_t nRowMax    = fGeo->GetRowMax(plane,chamber,sector);
1183     Int_t nColMax    = fGeo->GetColMax(plane);
1184
1185     Double_t *inADC  = new Double_t[nTimeTotal];
1186     Double_t *outADC = new Double_t[nTimeTotal];
1187
1188     // Add a container for the digits of this detector
1189     digits = fDigitsManager->GetDigits(iDet);
1190     // Allocate memory space for the digits buffer
1191     if (digits->GetNtime() == 0) {
1192       digits->Allocate(nRowMax,nColMax,nTimeTotal);
1193     }
1194
1195     // Get the signal container
1196     signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
1197     if (signals->GetNtime() == 0) {
1198       // Create missing containers
1199       signals->Allocate(nRowMax,nColMax,nTimeTotal);
1200     }
1201     else {
1202       // Expand the container if neccessary
1203       if (fCompress) signals->Expand();
1204     }
1205     // Create the missing dictionary containers
1206     for (iDict = 0; iDict < kNDict; iDict++) {
1207       dictionary[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1208       if (dictionary[iDict]->GetNtime() == 0) {
1209         dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1210       }
1211     }
1212
1213     Int_t nDigits = 0;
1214
1215     // Don't create noise in detectors that are switched off / not installed, etc.
1216     if (( calibration->IsChamberInstalled(iDet)) &&
1217         (!calibration->IsChamberMasked(iDet))    &&
1218         ( fGeo->ChamberInGeometry(iDet))) {
1219
1220       // Get the calibration objects
1221       calGainFactorROC      = calibration->GetGainFactorROC(iDet);
1222       calGainFactorDetValue = calGainFactorDet->GetValue(iDet);
1223
1224       // Create the digits for this chamber
1225       for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1226         for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1227
1228           // Check whether pad is masked
1229           // Bridged pads are not considered yet!!!
1230           if (calibration->IsPadMasked(iDet,iCol,iRow)) {
1231             continue;
1232           }
1233
1234           // Create summable digits
1235           if (fSDigits) {
1236
1237             for (iTime = 0; iTime < nTimeTotal; iTime++) {
1238               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1239               signalAmp *= fSDigitsScale;
1240               signalAmp  = TMath::Min(signalAmp,(Float_t) 1.0e9);
1241               Int_t adc  = (Int_t) signalAmp;
1242               if (adc > 0) nDigits++;
1243               digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1244             }
1245
1246           }
1247           // Create normal digits
1248           else {
1249
1250             Float_t padgain = calGainFactorDetValue
1251                             * calGainFactorROC->GetValue(iCol,iRow);
1252             if (padgain <= 0) {
1253               AliError(Form("Not a valid gain %f, %d %d %d",padgain,iDet,iCol,iRow));
1254             }
1255
1256             for (iTime = 0; iTime < nTimeTotal; iTime++) {
1257
1258               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1259
1260               // Pad and time coupling
1261               signalAmp *= coupling;
1262               // Gain factors
1263               signalAmp *= padgain;
1264
1265               // Add the noise, starting from minus ADC baseline in electrons
1266               Double_t baselineEl = simParam->GetADCbaseline()
1267                                   * (simParam->GetADCinRange() / simParam->GetADCoutRange())
1268                                   / convert;
1269               signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,simParam->GetNoise())
1270                                      ,-baselineEl);
1271               // Convert to mV
1272               signalAmp *= convert;
1273               // Add ADC baseline in mV
1274               signalAmp += simParam->GetADCbaseline()
1275                          * (simParam->GetADCinRange() / simParam->GetADCoutRange());
1276               // Convert to ADC counts. Set the overflow-bit fADCoutRange if the
1277               // signal is larger than fADCinRange
1278               Int_t adc  = 0;
1279               if (signalAmp >= simParam->GetADCinRange()) {
1280                 adc = ((Int_t) simParam->GetADCoutRange());
1281               }
1282               else {
1283                 adc = TMath::Nint(signalAmp * (simParam->GetADCoutRange()
1284                                             / simParam->GetADCinRange()));
1285               }
1286
1287               inADC[iTime]  = adc;
1288               outADC[iTime] = adc;
1289
1290             }
1291
1292             for (iTime = 0; iTime < nTimeTotal; iTime++) {
1293               //if (outADC[iTime] != 0 ) {   // Now this is enough because there is ZS in raw simulator
1294               // Reintroduce threshold as a temporary fix for the digits file size problem
1295               if (outADC[iTime] > (simParam->GetADCbaseline() + 3)) {
1296                 nDigits++;
1297                 digits->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
1298               }
1299             }
1300
1301           }
1302
1303         }
1304       }
1305
1306     }
1307
1308     // Compress the arrays
1309     digits->Compress(1,0);
1310     for (iDict = 0; iDict < kNDict; iDict++) {
1311       dictionary[iDict]->Compress(1,0);
1312     }
1313
1314     totalSizeDigits += digits->GetSize();
1315     totalSizeDict0  += dictionary[0]->GetSize();
1316     totalSizeDict1  += dictionary[1]->GetSize();
1317     totalSizeDict2  += dictionary[2]->GetSize();
1318
1319     if (nDigits > 0) {
1320       Float_t nPixel = nRowMax * nColMax * nTimeTotal;
1321       AliDebug(1,Form("Found %d digits in detector %d (%3.0f)."
1322                      ,nDigits,iDet
1323                      ,100.0 * ((Float_t) nDigits) / nPixel));
1324     }
1325
1326     if (fCompress) {
1327       signals->Compress(1,0);
1328     }
1329
1330     delete [] inADC;
1331     delete [] outADC;
1332
1333   }
1334
1335   if (signalsArray) {
1336     delete signalsArray;
1337     signalsArray = 0;
1338   }
1339
1340   AliDebug(1,Form("Total number of analyzed hits = %d",countHits));
1341   AliDebug(1,Form("Total digits data size = %d, %d, %d, %d",totalSizeDigits
1342                                                            ,totalSizeDict0
1343                                                            ,totalSizeDict1
1344                                                            ,totalSizeDict2));
1345
1346   return kTRUE;
1347
1348 }
1349
1350 //_____________________________________________________________________________
1351 void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
1352 {
1353   //
1354   // Add a digits manager for s-digits to the input list.
1355   //
1356
1357   fSDigitsManagerList->Add(man);
1358
1359 }
1360
1361 //_____________________________________________________________________________
1362 void AliTRDdigitizer::DeleteSDigitsManager()
1363 {
1364   //
1365   // Removes digits manager from the input list.
1366   //
1367
1368   fSDigitsManagerList->Delete();
1369
1370 }
1371
1372 //_____________________________________________________________________________
1373 Bool_t AliTRDdigitizer::ConvertSDigits()
1374 {
1375   //
1376   // Converts s-digits to normal digits
1377   //
1378
1379   // Number of track dictionary arrays
1380   const Int_t    kNDict = AliTRDdigitsManager::kNDict;
1381
1382   // Converts number of electrons to fC
1383   const Double_t kEl2fC = 1.602e-19 * 1.0e15;
1384
1385   Int_t iDict = 0;
1386   Int_t iRow;
1387   Int_t iCol;
1388   Int_t iTime;
1389
1390   AliTRDCalROC *calGainFactorROC      = 0;
1391   Float_t       calGainFactorDetValue = 0.0;
1392
1393   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
1394   if (!simParam) {
1395     AliFatal("Could not get simulation parameters");
1396     return kFALSE;
1397   }
1398
1399   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1400   if (!calibration) {
1401     AliFatal("Could not get calibration object");
1402     return kFALSE;
1403   }
1404
1405   Double_t sDigitsScale = 1.0 / GetSDigitsScale();
1406   Double_t noise        = simParam->GetNoise();
1407   Double_t padCoupling  = simParam->GetPadCoupling();
1408   Double_t timeCoupling = simParam->GetTimeCoupling();
1409   Double_t chipGain     = simParam->GetChipGain();
1410   Double_t coupling     = padCoupling * timeCoupling;
1411   Double_t convert      = kEl2fC * chipGain;
1412   Double_t adcInRange   = simParam->GetADCinRange();
1413   Double_t adcOutRange  = simParam->GetADCoutRange();
1414   Int_t    adcBaseline  = simParam->GetADCbaseline();
1415   Int_t    nTimeTotal   = calibration->GetNumberOfTimeBins();
1416
1417   AliTRDdataArrayI *digitsIn;
1418   AliTRDdataArrayI *digitsOut;
1419   AliTRDdataArrayI *dictionaryIn[kNDict];
1420   AliTRDdataArrayI *dictionaryOut[kNDict];
1421
1422   // Get the detector wise calibration objects
1423   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
1424
1425   // Loop through the detectors
1426   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1427
1428     Int_t plane      = fGeo->GetPlane(iDet);
1429     Int_t sector     = fGeo->GetSector(iDet);
1430     Int_t chamber    = fGeo->GetChamber(iDet);
1431     Int_t nRowMax    = fGeo->GetRowMax(plane,chamber,sector);
1432     Int_t nColMax    = fGeo->GetColMax(plane);
1433
1434     Double_t *inADC  = new Double_t[nTimeTotal];
1435     Double_t *outADC = new Double_t[nTimeTotal];
1436
1437     digitsIn  = fSDigitsManager->GetDigits(iDet);
1438     digitsIn->Expand();
1439     digitsOut = fDigitsManager->GetDigits(iDet);
1440     digitsOut->Allocate(nRowMax,nColMax,nTimeTotal);
1441     for (iDict = 0; iDict < kNDict; iDict++) {
1442       dictionaryIn[iDict]  = fSDigitsManager->GetDictionary(iDet,iDict);
1443       dictionaryIn[iDict]->Expand();
1444       dictionaryOut[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1445       dictionaryOut[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1446     }
1447
1448     // Don't create noise in detectors that are switched off / not installed, etc.
1449     if (( calibration->IsChamberInstalled(iDet)) &&
1450         (!calibration->IsChamberMasked(iDet))    &&
1451         ( fGeo->ChamberInGeometry(iDet))) {
1452
1453       // Get the calibration objects
1454       calGainFactorROC      = calibration->GetGainFactorROC(iDet);
1455       calGainFactorDetValue = calGainFactorDet->GetValue(iDet);
1456
1457       for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1458         for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1459
1460           // Check whether pad is masked
1461           // Bridged pads are not considered yet!!!
1462           if (calibration->IsPadMasked(iDet,iCol,iRow)) {
1463             continue;
1464           }
1465
1466           Float_t padgain = calGainFactorDetValue
1467                           * calGainFactorROC->GetValue(iCol,iRow);
1468           if (padgain <= 0) {
1469             AliError(Form("Not a valid gain %f, %d %d %d",padgain,iDet,iCol,iRow));
1470           }
1471
1472           for (iTime = 0; iTime < nTimeTotal; iTime++) {
1473
1474             // Scale s-digits to normal digits
1475             Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1476             signal         *= sDigitsScale;
1477             // Pad and time coupling
1478             signal *= coupling;
1479             // Gain factors
1480             signal *= padgain;
1481             // Add the noise, starting from minus ADC baseline in electrons
1482             Double_t baselineEl = adcBaseline * (adcInRange / adcOutRange) / convert;
1483             signal  = TMath::Max((Double_t) gRandom->Gaus(signal,noise),-baselineEl);
1484             // Convert to mV
1485             signal *= convert;
1486             // add ADC baseline in mV
1487             signal += adcBaseline * (adcInRange / adcOutRange);
1488             // Convert to ADC counts. Set the overflow-bit adcOutRange if the
1489             // signal is larger than adcInRange
1490             Int_t adc  = 0;
1491             if (signal >= adcInRange) {
1492               adc = ((Int_t) adcOutRange);
1493             }
1494             else {
1495               adc = TMath::Nint(signal * (adcOutRange / adcInRange));
1496             }
1497             inADC[iTime]  = adc;
1498             outADC[iTime] = adc;
1499
1500           }
1501
1502           for (iTime = 0; iTime < nTimeTotal; iTime++) {
1503             // Store the amplitude of the digit if above threshold
1504             //if (outADC[iTime] != 0) {  // now this is ok because there is ZS in raw simulation
1505             // Reintroduce the threshold as a temporary fix for the digits file size
1506             if (outADC[iTime] > (adcBaseline + 3)) {
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 }