TRD/AliTRDdigitizer.cxx

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