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