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