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