TRD/AliTRDdigitizer.cxx

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