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