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