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