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