Inherits now from AliDetector
[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 /*
17 $Log$
18 Revision 1.27  2001/11/14 10:50:45  cblume
19 Changes in digits IO. Add merging of summable digits
20
21 Revision 1.26  2001/11/06 17:19:41  cblume
22 Add detailed geometry and simple simulator
23
24 Revision 1.25  2001/06/27 09:54:44  cblume
25 Moved fField initialization to InitDetector()
26
27 Revision 1.24  2001/05/21 16:45:47  hristov
28 Last minute changes (C.Blume)
29
30 Revision 1.23  2001/05/07 08:04:48  cblume
31 New TRF and PRF. Speedup of the code. Digits from amplification region included
32
33 Revision 1.22  2001/03/30 14:40:14  cblume
34 Update of the digitization parameter
35
36 Revision 1.21  2001/03/13 09:30:35  cblume
37 Update of digitization. Moved digit branch definition to AliTRD
38
39 Revision 1.20  2001/02/25 20:19:00  hristov
40 Minor correction: loop variable declared only once for HP, Sun
41
42 Revision 1.19  2001/02/14 18:22:26  cblume
43 Change in the geometry of the padplane
44
45 Revision 1.18  2001/01/26 19:56:57  hristov
46 Major upgrade of AliRoot code
47
48 Revision 1.17  2000/12/08 12:53:27  cblume
49 Change in Copy() function for HP-compiler
50
51 Revision 1.16  2000/12/07 12:20:46  cblume
52 Go back to array compression. Use sampled PRF to speed up digitization
53
54 Revision 1.15  2000/11/23 14:34:08  cblume
55 Fixed bug in expansion routine of arrays (initialize buffers properly)
56
57 Revision 1.14  2000/11/20 08:54:44  cblume
58 Switch off compression as default
59
60 Revision 1.13  2000/11/10 14:57:52  cblume
61 Changes in the geometry constants for the DEC compiler
62
63 Revision 1.12  2000/11/01 14:53:20  cblume
64 Merge with TRD-develop
65
66 Revision 1.1.4.9  2000/10/26 17:00:22  cblume
67 Fixed bug in CheckDetector()
68
69 Revision 1.1.4.8  2000/10/23 13:41:35  cblume
70 Added protection against Log(0) in the gas gain calulation
71
72 Revision 1.1.4.7  2000/10/17 02:27:34  cblume
73 Get rid of global constants
74
75 Revision 1.1.4.6  2000/10/16 01:16:53  cblume
76 Changed timebin 0 to be the one closest to the readout
77
78 Revision 1.1.4.5  2000/10/15 23:34:29  cblume
79 Faster version of the digitizer
80
81 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
82 Made Getters const
83
84 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
85 Replace include files by forward declarations
86
87 Revision 1.1.4.2  2000/09/22 14:41:10  cblume
88 Bug fix in PRF. Included time response. New structure
89
90 Revision 1.10  2000/10/05 07:27:53  cblume
91 Changes in the header-files by FCA
92
93 Revision 1.9  2000/10/02 21:28:19  fca
94 Removal of useless dependecies via forward declarations
95
96 Revision 1.8  2000/06/09 11:10:07  cblume
97 Compiler warnings and coding conventions, next round
98
99 Revision 1.7  2000/06/08 18:32:58  cblume
100 Make code compliant to coding conventions
101
102 Revision 1.6  2000/06/07 16:27:32  cblume
103 Try to remove compiler warnings on Sun and HP
104
105 Revision 1.5  2000/05/09 16:38:57  cblume
106 Removed PadResponse(). Merge problem
107
108 Revision 1.4  2000/05/08 15:53:45  cblume
109 Resolved merge conflict
110
111 Revision 1.3  2000/04/28 14:49:27  cblume
112 Only one declaration of iDict in MakeDigits()
113
114 Revision 1.1.4.1  2000/05/08 14:42:04  cblume
115 Introduced AliTRDdigitsManager
116
117 Revision 1.1  2000/02/28 19:00:13  cblume
118 Add new TRD classes
119
120 */
121
122 ///////////////////////////////////////////////////////////////////////////////
123 //                                                                           //
124 //  Creates and handles digits from TRD hits                                 //
125 //  Author: C. Blume (C.Blume@gsi.de)                                        //
126 //                                                                           //
127 //  The following effects are included:                                      //
128 //      - Diffusion                                                          //
129 //      - ExB effects                                                        //
130 //      - Gas gain including fluctuations                                    //
131 //      - Pad-response (simple Gaussian approximation)                       //
132 //      - Time-response                                                      //
133 //      - Electronics noise                                                  //
134 //      - Electronics gain                                                   //
135 //      - Digitization                                                       //
136 //      - ADC threshold                                                      //
137 //  The corresponding parameter can be adjusted via the various              //
138 //  Set-functions. If these parameters are not explicitly set, default       //
139 //  values are used (see Init-function).                                     //
140 //  As an example on how to use this class to produce digits from hits       //
141 //  have a look at the macro hits2digits.C                                   //
142 //  The production of summable digits is demonstrated in hits2sdigits.C      //
143 //  and the subsequent conversion of the s-digits into normal digits is      //
144 //  explained in sdigits2digits.C.                                           //
145 //                                                                           //
146 ///////////////////////////////////////////////////////////////////////////////
147
148 #include <stdlib.h>
149
150 #include <TMath.h>
151 #include <TVector.h>
152 #include <TRandom.h>
153 #include <TROOT.h>
154 #include <TTree.h>
155 #include <TFile.h>
156 #include <TF1.h>
157 #include <TList.h>
158 #include <TTask.h>
159
160 #include "AliRun.h"
161 #include "AliMagF.h"
162 #include "AliRunDigitizer.h"
163
164 #include "AliTRD.h"
165 #include "AliTRDhit.h"
166 #include "AliTRDdigitizer.h"
167 #include "AliTRDdataArrayI.h"
168 #include "AliTRDdataArrayF.h"
169 #include "AliTRDsegmentArray.h"
170 #include "AliTRDdigitsManager.h"
171 #include "AliTRDgeometry.h"
172
173 ClassImp(AliTRDdigitizer)
174
175 //_____________________________________________________________________________
176 AliTRDdigitizer::AliTRDdigitizer()
177 {
178   //
179   // AliTRDdigitizer default constructor
180   //
181
182   fInputFile          = NULL;
183   fDigitsManager      = NULL;
184   fSDigitsManagerList = NULL;
185   fSDigitsManager     = NULL;
186   fTRD                = NULL;
187   fGeo                = NULL;
188   fPRFsmp             = NULL;
189   fTRFsmp             = NULL;
190
191   fEvent              = 0;
192   fGasGain            = 0.0;
193   fNoise              = 0.0;
194   fChipGain           = 0.0;
195   fADCoutRange        = 0.0;
196   fADCinRange         = 0.0;
197   fADCthreshold       = 0;
198   fDiffusionOn        = 0;
199   fDiffusionT         = 0.0;
200   fDiffusionL         = 0.0;
201   fElAttachOn         = 0;
202   fElAttachProp       = 0.0;
203   fExBOn              = 0;
204   fOmegaTau           = 0.0;
205   fPRFOn              = 0;
206   fTRFOn              = 0;
207   fDriftVelocity      = 0.0;
208   fPadCoupling        = 0.0;
209   fTimeCoupling       = 0.0;
210   fTimeBinWidth       = 0.0;
211   fField              = 0.0;
212
213   fPRFbin             = 0;
214   fPRFlo              = 0.0;
215   fPRFhi              = 0.0;
216   fPRFwid             = 0.0;
217   fPRFpad             = 0;
218   fTRFbin             = 0;
219   fTRFlo              = 0.0;
220   fTRFhi              = 0.0;
221   fTRFwid             = 0.0;
222
223   fCompress           = kTRUE;
224   fVerbose            = 0;
225   fSDigits            = kFALSE;
226   fSDigitsScale       = 0.0;
227
228 }
229
230 //_____________________________________________________________________________
231 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
232                 :AliDigitizer(name,title)
233 {
234   //
235   // AliTRDdigitizer constructor
236   //
237
238   fInputFile          = NULL;
239
240   fDigitsManager      = NULL;
241   fSDigitsManager     = NULL;
242   fSDigitsManagerList = NULL;
243
244   fTRD                = NULL;
245   fGeo                = NULL;
246   fPRFsmp             = NULL;
247   fTRFsmp             = NULL;
248
249   fEvent              = 0;
250
251   fCompress           = kTRUE;
252   fVerbose            = 0;
253   fSDigits            = kFALSE;
254
255   Init();
256
257 }
258
259 //_____________________________________________________________________________
260 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager
261                                 , const Text_t *name, const Text_t *title)
262                 :AliDigitizer(manager,name,title)
263 {
264   //
265   // AliTRDdigitizer constructor
266   //
267
268   fInputFile          = NULL;
269
270   fDigitsManager      = NULL;
271   fSDigitsManager     = NULL;
272   fSDigitsManagerList = NULL;
273
274   fTRD                = NULL;
275   fGeo                = NULL;
276   fPRFsmp             = NULL;
277   fTRFsmp             = NULL;
278
279   fEvent              = 0;
280
281   fCompress           = kTRUE;
282   fVerbose            = 0;
283   fSDigits            = kFALSE;
284
285   Init();
286
287 }
288
289 //_____________________________________________________________________________
290 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
291 {
292   //
293   // AliTRDdigitizer copy constructor
294   //
295
296   ((AliTRDdigitizer &) d).Copy(*this);
297
298 }
299
300 //_____________________________________________________________________________
301 AliTRDdigitizer::~AliTRDdigitizer()
302 {
303   //
304   // AliTRDdigitizer destructor
305   //
306
307   if (fInputFile) {
308     fInputFile->Close();
309     delete fInputFile;
310     fInputFile = NULL;
311   }
312
313   if (fDigitsManager) {
314     delete fDigitsManager;
315     fDigitsManager = NULL;
316   }
317
318   if (fSDigitsManager) {
319     delete fSDigitsManager;
320     fSDigitsManager = NULL;
321   }
322
323   if (fSDigitsManagerList) {
324     fSDigitsManagerList->Delete();
325     delete fSDigitsManagerList;
326     fSDigitsManagerList = NULL;
327   }
328 }
329
330 //_____________________________________________________________________________
331 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
332 {
333   //
334   // Assignment operator
335   //
336
337   if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
338   return *this;
339
340 }
341
342 //_____________________________________________________________________________
343 void AliTRDdigitizer::Copy(TObject &d)
344 {
345   //
346   // Copy function
347   //
348
349   Int_t iBin;
350
351   ((AliTRDdigitizer &) d).fInputFile          = NULL;
352   ((AliTRDdigitizer &) d).fSDigitsManagerList = NULL;
353   ((AliTRDdigitizer &) d).fSDigitsManager     = NULL;
354   ((AliTRDdigitizer &) d).fDigitsManager      = NULL;
355   ((AliTRDdigitizer &) d).fTRD                = NULL;
356   ((AliTRDdigitizer &) d).fGeo                = NULL;
357
358   ((AliTRDdigitizer &) d).fEvent              = 0;
359
360   ((AliTRDdigitizer &) d).fGasGain            = fGasGain;
361   ((AliTRDdigitizer &) d).fNoise              = fNoise;
362   ((AliTRDdigitizer &) d).fChipGain           = fChipGain;
363   ((AliTRDdigitizer &) d).fADCoutRange        = fADCoutRange;
364   ((AliTRDdigitizer &) d).fADCinRange         = fADCinRange;
365   ((AliTRDdigitizer &) d).fADCthreshold       = fADCthreshold;
366   ((AliTRDdigitizer &) d).fDiffusionOn        = fDiffusionOn; 
367   ((AliTRDdigitizer &) d).fDiffusionT         = fDiffusionT;
368   ((AliTRDdigitizer &) d).fDiffusionL         = fDiffusionL;
369   ((AliTRDdigitizer &) d).fElAttachOn         = fElAttachOn;
370   ((AliTRDdigitizer &) d).fElAttachProp       = fElAttachProp;
371   ((AliTRDdigitizer &) d).fExBOn              = fExBOn;
372   ((AliTRDdigitizer &) d).fOmegaTau           = fOmegaTau;
373   ((AliTRDdigitizer &) d).fLorentzFactor      = fLorentzFactor;
374   ((AliTRDdigitizer &) d).fDriftVelocity      = fDriftVelocity;
375   ((AliTRDdigitizer &) d).fPadCoupling        = fPadCoupling;
376   ((AliTRDdigitizer &) d).fTimeCoupling       = fTimeCoupling;
377   ((AliTRDdigitizer &) d).fTimeBinWidth       = fTimeBinWidth;
378   ((AliTRDdigitizer &) d).fField              = fField;
379   ((AliTRDdigitizer &) d).fPRFOn              = fPRFOn;
380   ((AliTRDdigitizer &) d).fTRFOn              = fTRFOn;
381
382   ((AliTRDdigitizer &) d).fCompress           = fCompress;
383   ((AliTRDdigitizer &) d).fVerbose            = fVerbose;
384   ((AliTRDdigitizer &) d).fSDigits            = fSDigits;
385   ((AliTRDdigitizer &) d).fSDigitsScale       = fSDigitsScale;
386
387   ((AliTRDdigitizer &) d).fPRFbin             = fPRFbin;
388   ((AliTRDdigitizer &) d).fPRFlo              = fPRFlo;
389   ((AliTRDdigitizer &) d).fPRFhi              = fPRFhi;
390   ((AliTRDdigitizer &) d).fPRFwid             = fPRFwid;
391   ((AliTRDdigitizer &) d).fPRFpad             = fPRFpad;
392   if (((AliTRDdigitizer &) d).fPRFsmp) delete ((AliTRDdigitizer &) d).fPRFsmp;
393   ((AliTRDdigitizer &) d).fPRFsmp = new Float_t[fPRFbin];
394   for (iBin = 0; iBin < fPRFbin; iBin++) {
395     ((AliTRDdigitizer &) d).fPRFsmp[iBin] = fPRFsmp[iBin];
396   }                                                                             
397   ((AliTRDdigitizer &) d).fTRFbin             = fTRFbin;
398   ((AliTRDdigitizer &) d).fTRFlo              = fTRFlo;
399   ((AliTRDdigitizer &) d).fTRFhi              = fTRFhi;
400   ((AliTRDdigitizer &) d).fTRFwid             = fTRFwid;
401   if (((AliTRDdigitizer &) d).fTRFsmp) delete ((AliTRDdigitizer &) d).fTRFsmp;
402   ((AliTRDdigitizer &) d).fTRFsmp = new Float_t[fTRFbin];
403   for (iBin = 0; iBin < fTRFbin; iBin++) {
404     ((AliTRDdigitizer &) d).fTRFsmp[iBin] = fTRFsmp[iBin];
405   }                                      
406                                        
407 }
408
409 //_____________________________________________________________________________
410 Int_t AliTRDdigitizer::Diffusion(Float_t driftlength, Float_t *xyz)
411 {
412   //
413   // Applies the diffusion smearing to the position of a single electron
414   //
415
416   Float_t driftSqrt = TMath::Sqrt(driftlength);
417   Float_t sigmaT = driftSqrt * fDiffusionT;
418   Float_t sigmaL = driftSqrt * fDiffusionL;
419   xyz[0] = gRandom->Gaus(xyz[0], sigmaL * fLorentzFactor);
420   xyz[1] = gRandom->Gaus(xyz[1], sigmaT * fLorentzFactor);
421   xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
422
423   return 1;
424
425 }
426
427 //_____________________________________________________________________________
428 Int_t AliTRDdigitizer::ExB(Float_t driftlength, Float_t *xyz)
429 {
430   //
431   // Applies E x B effects to the position of a single electron
432   //
433
434   xyz[0] = xyz[0];
435   xyz[1] = xyz[1] + fOmegaTau * driftlength;
436   xyz[2] = xyz[2];
437
438   return 1;
439
440 }
441
442 //_____________________________________________________________________________
443 Int_t AliTRDdigitizer::PadResponse(Float_t signal, Float_t dist, Float_t *pad)
444 {
445   //
446   // Applies the pad response
447   //
448
449   Int_t iBin =  ((Int_t) (( - dist - fPRFlo) / fPRFwid));
450
451   Int_t iBin0 = iBin - fPRFpad;
452   Int_t iBin1 = iBin;
453   Int_t iBin2 = iBin + fPRFpad;
454
455   if ((iBin0 >= 0) && (iBin2 < fPRFbin)) {
456
457     pad[0] = signal * fPRFsmp[iBin0];
458     pad[1] = signal * fPRFsmp[iBin1];
459     pad[2] = signal * fPRFsmp[iBin2];
460
461     return 1;
462
463   }
464   else {
465
466     return 0;
467
468   }
469
470 }
471
472 //_____________________________________________________________________________
473 Float_t AliTRDdigitizer::TimeResponse(Float_t time)
474 {
475   //
476   // Applies the preamp shaper time response
477   //
478
479   Int_t iBin = ((Int_t) ((time - fTRFlo) / fTRFwid)); 
480   if ((iBin >= 0) && (iBin < fTRFbin)) {
481     return fTRFsmp[iBin];
482   }
483   else {
484     return 0.0;
485   }    
486
487 }
488
489 //_____________________________________________________________________________
490 Bool_t AliTRDdigitizer::Init()
491 {
492   //
493   // Initializes the digitization procedure with standard values
494   //
495
496   // The default parameter for the digitization
497   fGasGain        = 2800.;
498   fChipGain       = 6.1;
499   fNoise          = 1000.;
500   fADCoutRange    = 1023.;          // 10-bit ADC
501   fADCinRange     = 1000.;          // 1V input range
502   fADCthreshold   = 1;
503
504   // For the summable digits
505   fSDigitsScale   = 100.;
506
507   // The drift velocity (cm / mus)
508   fDriftVelocity  = 1.5;
509
510   // Diffusion on
511   fDiffusionOn    = 1;
512
513   // E x B effects
514   fExBOn          = 0;
515
516   // Propability for electron attachment
517   fElAttachOn     = 0;
518   fElAttachProp   = 0.0;
519
520   // The pad response function
521   fPRFOn          =  1;
522
523   // The time response function
524   fTRFOn          =  1;
525
526   // The pad coupling factor (same number as for the TPC)
527   fPadCoupling    = 0.5;
528
529   // The time coupling factor (same number as for the TPC)
530   fTimeCoupling   = 0.4;
531
532   return kTRUE;
533
534 }
535
536 //_____________________________________________________________________________
537 Bool_t AliTRDdigitizer::ReInit()
538 {
539   //
540   // Reinitializes the digitization procedure after a change in the parameter
541   //
542
543   if (!fGeo) {
544     printf("AliTRDdigitizer::ReInit -- ");
545     printf("No geometry defined. Run InitDetector() first\n");
546     return kFALSE;
547   }
548
549   // Calculate the time bin width in ns
550   fTimeBinWidth   = fGeo->GetTimeBinSize() / fDriftVelocity * 1000.0;
551
552   // The range and the binwidth for the sampled TRF 
553   fTRFbin = 100;
554   // Start 0.2 mus before the signal
555   fTRFlo  = -0.2 * fDriftVelocity;
556   // End the maximum driftlength after the signal 
557   fTRFhi  = AliTRDgeometry::DrThick() 
558           + fGeo->GetTimeAfter() * fGeo->GetTimeBinSize();
559   fTRFwid = (fTRFhi - fTRFlo) / ((Float_t) fTRFbin);
560
561   // Transverse and longitudinal diffusion coefficients (Xe/CO2)
562   fDiffusionT     = GetDiffusionT(fDriftVelocity,fField);
563   fDiffusionL     = GetDiffusionL(fDriftVelocity,fField);
564
565   // omega * tau.= tan(Lorentz-angle)
566   fOmegaTau       = GetOmegaTau(fDriftVelocity,fField);
567
568   // The Lorentz factor
569   if (fExBOn) {
570     fLorentzFactor = 1.0 / (1.0 + fOmegaTau*fOmegaTau);
571   }
572   else {
573     fLorentzFactor = 1.0;
574   }
575
576   return kTRUE;
577
578 }
579
580 //_____________________________________________________________________________
581 void AliTRDdigitizer::SampleTRF()
582 {
583   //
584   // Samples the time response function
585   // It is defined according to Vasiles simulation of the preamp shaper
586   // output and includes the effect of the ion tail (based on Tariqs 
587   // Garfield simulation) and a shaping time of 125 ns FWHM
588   //
589
590   Int_t   ipos1;
591   Int_t   ipos2;
592   Float_t diff;
593
594   const Int_t kNpasa = 200;
595   Float_t time[kNpasa]   = {  -0.280000, -0.270000, -0.260000, -0.250000
596                             , -0.240000, -0.230000, -0.220000, -0.210000
597                             , -0.200000, -0.190000, -0.180000, -0.170000
598                             , -0.160000, -0.150000, -0.140000, -0.130000
599                             , -0.120000, -0.110000, -0.100000, -0.090000
600                             , -0.080000, -0.070000, -0.060000, -0.050000
601                             , -0.040000, -0.030000, -0.020000, -0.010000
602                             , -0.000000,  0.010000,  0.020000,  0.030000
603                             ,  0.040000,  0.050000,  0.060000,  0.070000
604                             ,  0.080000,  0.090000,  0.100000,  0.110000
605                             ,  0.120000,  0.130000,  0.140000,  0.150000
606                             ,  0.160000,  0.170000,  0.180000,  0.190000
607                             ,  0.200000,  0.210000,  0.220000,  0.230000
608                             ,  0.240000,  0.250000,  0.260000,  0.270000
609                             ,  0.280000,  0.290000,  0.300000,  0.310000
610                             ,  0.320000,  0.330000,  0.340000,  0.350000
611                             ,  0.360000,  0.370000,  0.380000,  0.390000
612                             ,  0.400000,  0.410000,  0.420000,  0.430000
613                             ,  0.440000,  0.450000,  0.460000,  0.470000
614                             ,  0.480000,  0.490000,  0.500000,  0.510000
615                             ,  0.520000,  0.530000,  0.540000,  0.550000
616                             ,  0.560000,  0.570000,  0.580000,  0.590000
617                             ,  0.600000,  0.610000,  0.620000,  0.630000
618                             ,  0.640000,  0.650000,  0.660000,  0.670000
619                             ,  0.680000,  0.690000,  0.700000,  0.710000
620                             ,  0.720000,  0.730000,  0.740000,  0.750000
621                             ,  0.760000,  0.770000,  0.780000,  0.790000
622                             ,  0.800000,  0.810000,  0.820000,  0.830000
623                             ,  0.840000,  0.850000,  0.860000,  0.870000
624                             ,  0.880000,  0.890000,  0.900000,  0.910000
625                             ,  0.920000,  0.930000,  0.940000,  0.950000
626                             ,  0.960000,  0.970000,  0.980000,  0.990000
627                             ,  1.000000,  1.010000,  1.020000,  1.030000
628                             ,  1.040000,  1.050000,  1.060000,  1.070000
629                             ,  1.080000,  1.090000,  1.100000,  1.110000
630                             ,  1.120000,  1.130000,  1.140000,  1.150000
631                             ,  1.160000,  1.170000,  1.180000,  1.190000
632                             ,  1.200000,  1.210000,  1.220000,  1.230000
633                             ,  1.240000,  1.250000,  1.260000,  1.270000
634                             ,  1.280000,  1.290000,  1.300000,  1.310000
635                             ,  1.320000,  1.330000,  1.340000,  1.350000
636                             ,  1.360000,  1.370000,  1.380000,  1.390000
637                             ,  1.400000,  1.410000,  1.420000,  1.430000
638                             ,  1.440000,  1.450000,  1.460000,  1.470000
639                             ,  1.480000,  1.490000,  1.500000,  1.510000
640                             ,  1.520000,  1.530000,  1.540000,  1.550000
641                             ,  1.560000,  1.570000,  1.580000,  1.590000
642                             ,  1.600000,  1.610000,  1.620000,  1.630000
643                             ,  1.640000,  1.650000,  1.660000,  1.670000
644                             ,  1.680000,  1.690000,  1.700000,  1.710000 };
645
646   Float_t signal[kNpasa] = {   0.000000,  0.000000,  0.000000,  0.000000 
647                             ,  0.000000,  0.000000,  0.000000,  0.000000
648                             ,  0.000000,  0.000000,  0.000000,  0.000000
649                             ,  0.000000,  0.000000,  0.000000,  0.000098
650                             ,  0.003071,  0.020056,  0.066053,  0.148346
651                             ,  0.263120,  0.398496,  0.540226,  0.674436
652                             ,  0.790977,  0.883083,  0.947744,  0.985714
653                             ,  0.999248,  0.992105,  0.967669,  0.930827
654                             ,  0.884586,  0.833083,  0.778571,  0.723684
655                             ,  0.669173,  0.617293,  0.567669,  0.521805
656                             ,  0.479699,  0.440977,  0.405639,  0.373985
657                             ,  0.345526,  0.320038,  0.297256,  0.276917
658                             ,  0.258797,  0.242632,  0.228195,  0.215301
659                             ,  0.203759,  0.193383,  0.184023,  0.175564
660                             ,  0.167895,  0.160940,  0.154549,  0.148722
661                             ,  0.143308,  0.138346,  0.133722,  0.129398
662                             ,  0.125376,  0.121617,  0.118045,  0.114699
663                             ,  0.111541,  0.108571,  0.105714,  0.103008
664                             ,  0.100414,  0.097970,  0.095602,  0.093346
665                             ,  0.091165,  0.089060,  0.087068,  0.085150
666                             ,  0.083308,  0.081541,  0.079812,  0.078158
667                             ,  0.076541,  0.075000,  0.073496,  0.072068
668                             ,  0.070677,  0.069286,  0.068008,  0.066729
669                             ,  0.065489,  0.064286,  0.063120,  0.061992
670                             ,  0.060902,  0.059850,  0.058797,  0.057820
671                             ,  0.056842,  0.055902,  0.054962,  0.054060
672                             ,  0.053158,  0.052293,  0.051466,  0.050639
673                             ,  0.049850,  0.049060,  0.048308,  0.047556
674                             ,  0.046842,  0.046128,  0.045451,  0.044774
675                             ,  0.044098,  0.043459,  0.042820,  0.042218
676                             ,  0.041617,  0.041015,  0.040451,  0.039887
677                             ,  0.039323,  0.038797,  0.038271,  0.037744
678                             ,  0.037237,  0.036744,  0.036259,  0.035786
679                             ,  0.035323,  0.034872,  0.034429,  0.033996
680                             ,  0.033575,  0.033162,  0.032756,  0.032361
681                             ,  0.031974,  0.031594,  0.031222,  0.030857
682                             ,  0.030496,  0.030143,  0.029793,  0.029451
683                             ,  0.029109,  0.028774,  0.028444,  0.028113
684                             ,  0.027793,  0.027477,  0.027165,  0.026861
685                             ,  0.026564,  0.026271,  0.025981,  0.025699
686                             ,  0.025421,  0.025147,  0.024880,  0.024613
687                             ,  0.024353,  0.024094,  0.023842,  0.023590
688                             ,  0.023346,  0.023102,  0.022865,  0.022628
689                             ,  0.022398,  0.022173,  0.021951,  0.021733
690                             ,  0.021519,  0.021308,  0.021098,  0.020891
691                             ,  0.020688,  0.020485,  0.020286,  0.020090
692                             ,  0.019895,  0.019707,  0.019519,  0.019335
693                             ,  0.019150,  0.018974,  0.018797,  0.018624
694                             ,  0.018451,  0.018282,  0.018113,  0.017947
695                             ,  0.017782,  0.017617,  0.017455,  0.017297 };
696
697   if (fTRFsmp) delete fTRFsmp;
698   fTRFsmp = new Float_t[fTRFbin];
699
700   Float_t loTRF    = TMath::Max(fTRFlo / fDriftVelocity,time[0]);
701   Float_t hiTRF    = TMath::Min(fTRFhi / fDriftVelocity,time[kNpasa-1]);
702   Float_t binWidth = (hiTRF - loTRF) / ((Float_t) fTRFbin);
703
704   // Take the linear interpolation
705   for (Int_t iBin = 0; iBin < fTRFbin; iBin++) {
706
707     Float_t bin = (((Float_t) iBin) + 0.5) * binWidth + loTRF;
708     ipos1 = ipos2 = 0;
709     diff  = 0;
710     do {
711       diff = bin - time[ipos2++];
712     } while (diff > 0);
713     ipos2--;
714     if (ipos2 > kNpasa) ipos2 = kNpasa - 1;
715     ipos1 = ipos2 - 1;
716
717     fTRFsmp[iBin] = signal[ipos2] 
718                   + diff * (signal[ipos2] - signal[ipos1]) 
719                          / (  time[ipos2] -   time[ipos1]);
720
721   }
722
723 }
724
725 //_____________________________________________________________________________
726 void AliTRDdigitizer::SamplePRF()
727 {
728   //
729   // Samples the pad response function
730   //
731
732   const Int_t kPRFbin = 61;
733   Float_t prf[kPRFbin] = { 0.002340, 0.003380, 0.004900, 0.007080, 0.010220
734                          , 0.014740, 0.021160, 0.030230, 0.042800, 0.059830
735                          , 0.082030, 0.109700, 0.142550, 0.179840, 0.220610
736                          , 0.263980, 0.309180, 0.355610, 0.402790, 0.450350
737                          , 0.497930, 0.545190, 0.591740, 0.637100, 0.680610
738                          , 0.721430, 0.758400, 0.790090, 0.814720, 0.830480
739                          , 0.835930, 0.830480, 0.814710, 0.790070, 0.758390
740                          , 0.721410, 0.680590, 0.637080, 0.591730, 0.545180
741                          , 0.497920, 0.450340, 0.402790, 0.355610, 0.309190
742                          , 0.263990, 0.220630, 0.179850, 0.142570, 0.109720
743                          , 0.082040, 0.059830, 0.042820, 0.030230, 0.021170
744                          , 0.014740, 0.010230, 0.007080, 0.004900, 0.003380
745                          , 0.002340 };
746
747   fPRFlo  = -1.5;
748   fPRFhi  =  1.5;
749   fPRFbin = kPRFbin;
750   fPRFwid = (fPRFhi - fPRFlo) / ((Float_t) fPRFbin);
751   fPRFpad = ((Int_t) (1.0 / fPRFwid));
752
753   if (fPRFsmp) delete fPRFsmp;
754   fPRFsmp = new Float_t[fPRFbin];
755   for (Int_t iBin = 0; iBin < fPRFbin; iBin++) {
756     fPRFsmp[iBin] = prf[iBin];
757   }
758
759 }
760
761 //_____________________________________________________________________________
762 Bool_t AliTRDdigitizer::Open(const Char_t *name, Int_t nEvent)
763 {
764   //
765   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
766   //
767
768   // Connect the AliRoot file containing Geometry, Kine, and Hits
769   fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
770   if (!fInputFile) {
771     if (fVerbose > 0) {
772       printf("AliTRDdigitizer::Open -- ");
773       printf("Open the AliROOT-file %s.\n",name);
774     }
775     fInputFile = new TFile(name,"UPDATE");
776   }
777   else {
778     if (fVerbose > 0) {
779       printf("AliTRDdigitizer::Open -- ");
780       printf("%s is already open.\n",name);
781     }
782   }
783
784   gAlice = (AliRun*) fInputFile->Get("gAlice");
785   if (gAlice) {
786     if (fVerbose > 0) {
787       printf("AliTRDdigitizer::Open -- ");
788       printf("AliRun object found on file.\n");
789     }
790   }
791   else {
792     printf("AliTRDdigitizer::Open -- ");
793     printf("Could not find AliRun object.\n");
794     return kFALSE;
795   }
796
797   fEvent = nEvent;
798
799   // Import the Trees for the event nEvent in the file
800   Int_t nparticles = gAlice->GetEvent(fEvent);
801   if (nparticles <= 0) {
802     printf("AliTRDdigitizer::Open -- ");
803     printf("No entries in the trees for event %d.\n",fEvent);
804     return kFALSE;
805   }
806
807   if (InitDetector()) {
808     return MakeBranch();
809   }
810   else {
811     return kFALSE;
812   }
813
814 }
815
816 //_____________________________________________________________________________
817 Bool_t AliTRDdigitizer::InitDetector()
818 {
819   //
820   // Sets the pointer to the TRD detector and the geometry
821   //
822
823   // Get the pointer to the detector class and check for version 1
824   fTRD = (AliTRD*) gAlice->GetDetector("TRD");
825   if (fTRD->IsVersion() != 1) {
826     printf("AliTRDdigitizer::InitDetector -- ");
827     printf("TRD must be version 1 (slow simulator).\n");
828     exit(1);
829   }
830
831   // Get the geometry
832   fGeo = fTRD->GetGeometry();
833   if (fVerbose > 0) {
834     printf("AliTRDdigitizer::InitDetector -- ");
835     printf("Geometry version %d\n",fGeo->IsVersion());
836   }
837
838   // The magnetic field strength in Tesla
839   fField = 0.2 * gAlice->Field()->Factor();
840
841   // Create a digits manager
842   fDigitsManager = new AliTRDdigitsManager();
843   fDigitsManager->SetSDigits(fSDigits);
844   fDigitsManager->CreateArrays();
845   fDigitsManager->SetEvent(fEvent);
846   fDigitsManager->SetVerbose(fVerbose);
847
848   // The list for the input s-digits manager to be merged
849   fSDigitsManagerList = new TList();
850
851   return ReInit();
852
853 }
854
855 //_____________________________________________________________________________
856 Bool_t AliTRDdigitizer::MakeBranch(const Char_t *file)
857 {
858   // 
859   // Create the branches for the digits array
860   //
861
862   return fDigitsManager->MakeBranch(file);
863
864 }
865
866 //_____________________________________________________________________________
867 Bool_t AliTRDdigitizer::MakeDigits()
868 {
869   //
870   // Creates digits.
871   //
872
873   ///////////////////////////////////////////////////////////////
874   // Parameter 
875   ///////////////////////////////////////////////////////////////
876
877   // Converts number of electrons to fC
878   const Double_t kEl2fC  = 1.602E-19 * 1.0E15; 
879
880   ///////////////////////////////////////////////////////////////
881
882   // Number of pads included in the pad response
883   const Int_t kNpad  = 3;
884
885   // Number of track dictionary arrays
886   const Int_t kNDict = AliTRDdigitsManager::kNDict;
887
888   // Half the width of the amplification region
889   const Float_t kAmWidth = AliTRDgeometry::AmThick() / 2.;
890
891   Int_t   iRow, iCol, iTime, iPad;
892   Int_t   iDict  = 0;
893   Int_t   nBytes = 0;
894
895   Int_t   totalSizeDigits = 0;
896   Int_t   totalSizeDict0  = 0;
897   Int_t   totalSizeDict1  = 0;
898   Int_t   totalSizeDict2  = 0;
899
900   Int_t   timeTRDbeg = 0;
901   Int_t   timeTRDend = 1;
902
903   Float_t pos[3];
904   Float_t rot[3];
905   Float_t xyz[3];
906   Float_t padSignal[kNpad];
907   Float_t signalOld[kNpad];
908
909   AliTRDdataArrayF *signals = 0;
910   AliTRDdataArrayI *digits  = 0;
911   AliTRDdataArrayI *dictionary[kNDict];
912
913   // Create a container for the amplitudes
914   AliTRDsegmentArray *signalsArray 
915                      = new AliTRDsegmentArray("AliTRDdataArrayF",AliTRDgeometry::Ndet());
916
917   if (fTRFOn) {
918     timeTRDbeg = ((Int_t) (-fTRFlo / fGeo->GetTimeBinSize())) - 1;
919     timeTRDend = ((Int_t) ( fTRFhi / fGeo->GetTimeBinSize())) - 1;
920     if (fVerbose > 0) {
921       printf("AliTRDdigitizer::MakeDigits -- ");
922       printf("Sample the TRF between -%d and %d\n",timeTRDbeg,timeTRDend);
923     }
924   }
925
926   Float_t elAttachProp = fElAttachProp / 100.; 
927
928   // Create the sampled PRF
929   SamplePRF();
930
931   // Create the sampled TRF
932   SampleTRF();
933
934   if (!fGeo) {
935     printf("AliTRDdigitizer::MakeDigits -- ");
936     printf("No geometry defined\n");
937     return kFALSE;
938   }
939
940   if (fVerbose > 0) {
941     printf("AliTRDdigitizer::MakeDigits -- ");
942     printf("Start creating digits.\n");
943   }
944
945   // Get the pointer to the hit tree
946   TTree *HitTree = gAlice->TreeH();
947
948   // Get the number of entries in the hit tree
949   // (Number of primary particles creating a hit somewhere)
950   Int_t nTrack = (Int_t) HitTree->GetEntries();
951   if (fVerbose > 0) {
952     printf("AliTRDdigitizer::MakeDigits -- ");
953     printf("Found %d primary particles\n",nTrack);
954   } 
955
956   Int_t detectorOld = -1;
957   Int_t countHits   =  0; 
958
959   // Loop through all entries in the tree
960   for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
961
962     gAlice->ResetHits();
963     nBytes += HitTree->GetEvent(iTrack);
964
965     // Get the number of hits in the TRD created by this particle
966     Int_t nHit = fTRD->Hits()->GetEntriesFast();
967     if (fVerbose > 0) {
968       printf("AliTRDdigitizer::MakeDigits -- ");
969       printf("Found %d hits for primary particle %d\n",nHit,iTrack);
970     }
971
972     // Loop through the TRD hits  
973     for (Int_t iHit = 0; iHit < nHit; iHit++) {
974
975       countHits++;
976
977       AliTRDhit *hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit);
978               pos[0]   = hit->X();
979               pos[1]   = hit->Y();
980               pos[2]   = hit->Z();
981       Float_t q        = hit->GetCharge();
982       Int_t   track    = hit->Track();
983       Int_t   detector = hit->GetDetector();
984       Int_t   plane    = fGeo->GetPlane(detector);
985       Int_t   sector   = fGeo->GetSector(detector);
986       Int_t   chamber  = fGeo->GetChamber(detector);
987
988       if (!(CheckDetector(plane,chamber,sector))) continue;
989
990       Int_t   nRowMax     = fGeo->GetRowMax(plane,chamber,sector);
991       Int_t   nColMax     = fGeo->GetColMax(plane);
992       Int_t   nTimeMax    = fGeo->GetTimeMax();
993       Int_t   nTimeBefore = fGeo->GetTimeBefore();
994       Int_t   nTimeAfter  = fGeo->GetTimeAfter();
995       Int_t   nTimeTotal  = fGeo->GetTimeTotal();
996       Float_t row0        = fGeo->GetRow0(plane,chamber,sector);
997       Float_t col0        = fGeo->GetCol0(plane);
998       Float_t time0       = fGeo->GetTime0(plane);
999       Float_t rowPadSize  = fGeo->GetRowPadSize(plane,chamber,sector);
1000       Float_t colPadSize  = fGeo->GetColPadSize(plane);
1001       Float_t timeBinSize = fGeo->GetTimeBinSize();
1002       Float_t divideRow   = 1.0 / rowPadSize;
1003       Float_t divideCol   = 1.0 / colPadSize;
1004       Float_t divideTime  = 1.0 / timeBinSize;
1005
1006       if (fVerbose > 1) {
1007         printf("Analyze hit no. %d ",iHit);
1008         printf("-----------------------------------------------------------\n");
1009         hit->Dump();
1010         printf("plane = %d, sector = %d, chamber = %d\n"
1011               ,plane,sector,chamber);
1012         printf("nRowMax = %d, nColMax = %d, nTimeMax = %d\n" 
1013               ,nRowMax,nColMax,nTimeMax);
1014         printf("nTimeBefore = %d, nTimeAfter = %d, nTimeTotal = %d\n"
1015               ,nTimeBefore,nTimeAfter,nTimeTotal);
1016         printf("row0 = %f, col0 = %f, time0 = %f\n"
1017               ,row0,col0,time0);
1018         printf("rowPadSize = %f, colPadSize = %f, timeBinSize = %f\n"
1019                ,rowPadSize,colPadSize,timeBinSize); 
1020       }
1021        
1022       // Don't analyze test hits
1023       if (hit->FromTest()) continue;
1024
1025       if (detector != detectorOld) {
1026
1027         if (fVerbose > 1) {
1028           printf("AliTRDdigitizer::MakeDigits -- ");
1029           printf("Get new container. New det = %d, Old det = %d\n"
1030                 ,detector,detectorOld);
1031         }
1032         // Compress the old one if enabled
1033         if ((fCompress) && (detectorOld > -1)) {
1034           if (fVerbose > 1) {
1035             printf("AliTRDdigitizer::MakeDigits -- ");
1036             printf("Compress the old container ...");
1037           }
1038           signals->Compress(1,0);
1039           for (iDict = 0; iDict < kNDict; iDict++) {
1040             dictionary[iDict]->Compress(1,0);
1041           }
1042           if (fVerbose > 1) printf("done\n");
1043         }
1044         // Get the new container
1045         signals = (AliTRDdataArrayF *) signalsArray->At(detector);
1046         if (signals->GetNtime() == 0) {
1047           // Allocate a new one if not yet existing
1048           if (fVerbose > 1) {
1049             printf("AliTRDdigitizer::MakeDigits -- ");
1050             printf("Allocate a new container ... ");
1051           }
1052           signals->Allocate(nRowMax,nColMax,nTimeTotal);
1053         }
1054         else {
1055           // Expand an existing one
1056           if (fCompress) {
1057             if (fVerbose > 1) {
1058               printf("AliTRDdigitizer::MakeDigits -- ");
1059               printf("Expand an existing container ... ");
1060             }
1061             signals->Expand();
1062           }
1063         }
1064         // The same for the dictionary
1065         for (iDict = 0; iDict < kNDict; iDict++) {       
1066           dictionary[iDict] = fDigitsManager->GetDictionary(detector,iDict);
1067           if (dictionary[iDict]->GetNtime() == 0) {
1068             dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1069           }
1070           else {
1071             if (fCompress) dictionary[iDict]->Expand();
1072           }
1073         }      
1074         if (fVerbose > 1) printf("done\n");
1075         detectorOld = detector;
1076       }
1077
1078       // Rotate the sectors on top of each other
1079       fGeo->Rotate(detector,pos,rot);
1080
1081       // The driftlength. It is negative if the hit is in the 
1082       // amplification region.
1083       Float_t driftlength = time0 - rot[0];
1084
1085       // Take also the drift in the amplification region into account
1086       // The drift length is at the moment still the same, regardless of
1087       // the position relativ to the wire. This non-isochronity needs still
1088       // to be implemented.
1089       Float_t driftlengthL = TMath::Abs(driftlength + kAmWidth);
1090       if (fExBOn) driftlengthL /= TMath::Sqrt(fLorentzFactor);
1091
1092       // Loop over all electrons of this hit
1093       // TR photons produce hits with negative charge
1094       Int_t nEl = ((Int_t) TMath::Abs(q));
1095       for (Int_t iEl = 0; iEl < nEl; iEl++) {
1096
1097         xyz[0] = rot[0];
1098         xyz[1] = rot[1];
1099         xyz[2] = rot[2];
1100
1101         // Electron attachment
1102         if (fElAttachOn) {
1103           if (gRandom->Rndm() < (driftlengthL * elAttachProp)) 
1104             continue;
1105         }
1106
1107         // Apply the diffusion smearing
1108         if (fDiffusionOn) {
1109           if (!(Diffusion(driftlengthL,xyz))) continue;
1110         }
1111
1112         // Apply E x B effects (depends on drift direction)
1113         if (fExBOn) { 
1114           if (!(ExB(driftlength+kAmWidth,xyz))) continue;   
1115         }
1116
1117         // The electron position after diffusion and ExB in pad coordinates 
1118         // The pad row (z-direction)
1119         Int_t  rowE = ((Int_t) ((xyz[2] -  row0) * divideRow));
1120         if ((rowE < 0) || (rowE >= nRowMax)) continue;   
1121
1122         // The pad column (rphi-direction)
1123         Int_t  colE = ((Int_t) ((xyz[1] -  col0) * divideCol));
1124         if ((colE < 0) || (colE >= nColMax)) continue;   
1125   
1126         // The time bin (negative for hits in the amplification region)
1127         // In the amplification region the electrons drift from both sides
1128         // to the middle (anode wire plane)
1129         Float_t timeDist   = time0 - xyz[0];
1130         Float_t timeOffset = 0;
1131         Int_t   timeE      = 0;
1132         if (timeDist > 0) {
1133           // The time bin
1134           timeE      = ((Int_t) (timeDist * divideTime));
1135           // The distance of the position to the middle of the timebin
1136           timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) - timeDist;
1137         }
1138         else {
1139           // Difference between half of the amplification gap width and
1140           // the distance to the anode wire
1141           Float_t anodeDist = kAmWidth - TMath::Abs(timeDist + kAmWidth);
1142           // The time bin
1143           timeE      = -1 * (((Int_t ) (anodeDist * divideTime)) + 1);
1144           // The distance of the position to the middle of the timebin
1145           timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) + anodeDist;
1146         }
1147  
1148         // Apply the gas gain including fluctuations
1149         Float_t ggRndm = 0.0;
1150         do {
1151           ggRndm = gRandom->Rndm();
1152         } while (ggRndm <= 0);
1153         Int_t signal = (Int_t) (-fGasGain * TMath::Log(ggRndm));
1154
1155         // Apply the pad response 
1156         if (fPRFOn) {
1157           // The distance of the electron to the center of the pad 
1158           // in units of pad width
1159           Float_t dist = (xyz[1] - col0 - (colE + 0.5) * colPadSize) 
1160                        * divideCol;
1161           if (!(PadResponse(signal,dist,padSignal))) continue;
1162         }
1163         else {
1164           padSignal[0] = 0.0;
1165           padSignal[1] = signal;
1166           padSignal[2] = 0.0;
1167         }
1168
1169         // Sample the time response inside the drift region
1170         // + additional time bins before and after.
1171         // The sampling is done always in the middle of the time bin
1172         for (Int_t iTimeBin = TMath::Max(timeE-timeTRDbeg,        -nTimeBefore) 
1173                   ;iTimeBin < TMath::Min(timeE+timeTRDend,nTimeMax+nTimeAfter ) 
1174                   ;iTimeBin++) {
1175
1176           // Apply the time response
1177           Float_t timeResponse = 1.0;
1178           if (fTRFOn) {
1179             Float_t time = (iTimeBin - timeE) * timeBinSize + timeOffset;
1180             timeResponse = TimeResponse(time);
1181           }
1182
1183           signalOld[0] = 0.0;
1184           signalOld[1] = 0.0;
1185           signalOld[2] = 0.0;
1186
1187           for (iPad = 0; iPad < kNpad; iPad++) {
1188
1189             Int_t colPos = colE + iPad - 1;
1190             if (colPos <        0) continue;
1191             if (colPos >= nColMax) break;
1192
1193             // Add the signals
1194             // Note: The time bin number is shifted by nTimeBefore to avoid negative
1195             // time bins. This has to be subtracted later.
1196             Int_t iCurrentTimeBin = iTimeBin + nTimeBefore;
1197             signalOld[iPad]  = signals->GetDataUnchecked(rowE,colPos,iCurrentTimeBin);
1198             signalOld[iPad] += padSignal[iPad] * timeResponse;
1199             signals->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,signalOld[iPad]);
1200
1201             // Store the track index in the dictionary
1202             // Note: We store index+1 in order to allow the array to be compressed
1203             if (signalOld[iPad] > 0) {
1204               for (iDict = 0; iDict < kNDict; iDict++) {
1205                 Int_t oldTrack = dictionary[iDict]->GetDataUnchecked(rowE
1206                                                                     ,colPos
1207                                                                     ,iCurrentTimeBin);
1208                 if (oldTrack == track+1) break;
1209                 if (oldTrack ==       0) {
1210                   dictionary[iDict]->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,track+1);
1211                   break;
1212                 }
1213               }
1214             }
1215
1216           }
1217
1218         }
1219
1220       }
1221
1222     }
1223
1224   } // All hits finished
1225
1226   if (fVerbose > 0) {
1227     printf("AliTRDdigitizer::MakeDigits -- ");
1228     printf("Finished analyzing %d hits\n",countHits);
1229   }
1230
1231   // The total conversion factor
1232   Float_t convert = kEl2fC * fPadCoupling * fTimeCoupling * fChipGain;
1233
1234   // Loop through all chambers to finalize the digits
1235   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1236
1237     Int_t plane       = fGeo->GetPlane(iDet);
1238     Int_t sector      = fGeo->GetSector(iDet);
1239     Int_t chamber     = fGeo->GetChamber(iDet);
1240     Int_t nRowMax     = fGeo->GetRowMax(plane,chamber,sector);
1241     Int_t nColMax     = fGeo->GetColMax(plane);
1242     Int_t nTimeMax    = fGeo->GetTimeMax();
1243     Int_t nTimeTotal  = fGeo->GetTimeTotal();
1244
1245     if (fVerbose > 0) {
1246       printf("AliTRDdigitizer::MakeDigits -- ");
1247       printf("Digitization for chamber %d\n",iDet);
1248     }
1249
1250     // Add a container for the digits of this detector
1251     digits = fDigitsManager->GetDigits(iDet);        
1252     // Allocate memory space for the digits buffer
1253     digits->Allocate(nRowMax,nColMax,nTimeTotal);
1254
1255     // Get the signal container
1256     signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
1257     if (signals->GetNtime() == 0) {
1258       // Create missing containers
1259       signals->Allocate(nRowMax,nColMax,nTimeTotal);      
1260     }
1261     else {
1262       // Expand the container if neccessary
1263       if (fCompress) signals->Expand();
1264     }
1265     // Create the missing dictionary containers
1266     for (iDict = 0; iDict < kNDict; iDict++) {       
1267       dictionary[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1268       if (dictionary[iDict]->GetNtime() == 0) {
1269         dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1270       }
1271     }
1272
1273     Int_t nDigits = 0;
1274
1275     // Don't create noise in detectors that are switched off
1276     if (CheckDetector(plane,chamber,sector)) {
1277
1278       // Create the digits for this chamber
1279       for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1280         for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1281           for (iTime = 0; iTime < nTimeTotal; iTime++) {         
1282
1283             // Create summable digits
1284             if (fSDigits) {
1285
1286               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1287               signalAmp *= fSDigitsScale;
1288               signalAmp  = TMath::Min(signalAmp,1.0e9);
1289               Int_t adc  = (Int_t) signalAmp;
1290               nDigits++;
1291               digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1292
1293             }
1294             // Create normal digits
1295             else {
1296
1297               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1298
1299               // Add the noise
1300               signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,fNoise),0.0);
1301               // Convert to mV
1302               signalAmp *= convert;
1303               // Convert to ADC counts. Set the overflow-bit fADCoutRange if the 
1304               // signal is larger than fADCinRange
1305               Int_t adc  = 0;
1306               if (signalAmp >= fADCinRange) {
1307                 adc = ((Int_t) fADCoutRange);
1308               }
1309               else {
1310                 adc = ((Int_t) (signalAmp * (fADCoutRange / fADCinRange)));
1311               }
1312
1313               // Store the amplitude of the digit if above threshold
1314               if (adc > fADCthreshold) {
1315                 if (fVerbose > 2) {
1316                   printf("  iRow = %d, iCol = %d, iTime = %d\n"
1317                         ,iRow,iCol,iTime);
1318                   printf("  signal = %f, adc = %d\n",signalAmp,adc);
1319                 }
1320                 nDigits++;
1321                 digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1322               }
1323
1324             }
1325
1326           }
1327         }
1328       }
1329
1330     }
1331
1332     // Compress the arrays
1333     digits->Compress(1,0);
1334     for (iDict = 0; iDict < kNDict; iDict++) {
1335       dictionary[iDict]->Compress(1,0);
1336     }
1337
1338     totalSizeDigits += digits->GetSize();
1339     totalSizeDict0  += dictionary[0]->GetSize();
1340     totalSizeDict1  += dictionary[1]->GetSize();
1341     totalSizeDict2  += dictionary[2]->GetSize();
1342
1343     Float_t nPixel = nRowMax * nColMax * nTimeMax;
1344     if (fVerbose > 0) {
1345       printf("AliTRDdigitizer::MakeDigits -- ");
1346       printf("Found %d digits in detector %d (%3.0f).\n"
1347             ,nDigits,iDet
1348             ,100.0 * ((Float_t) nDigits) / nPixel);
1349     } 
1350
1351     if (fCompress) signals->Compress(1,0);
1352
1353   }
1354
1355   if (fVerbose > 0) {
1356     printf("AliTRDdigitizer::MakeDigits -- ");
1357     printf("Total number of analyzed hits = %d\n",countHits);
1358     printf("AliTRDdigitizer::MakeDigits -- ");
1359     printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
1360                                                       ,totalSizeDict0
1361                                                       ,totalSizeDict1
1362                                                       ,totalSizeDict2);        
1363   }
1364
1365   return kTRUE;
1366
1367 }
1368
1369 //_____________________________________________________________________________
1370 void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
1371 {
1372   //
1373   // Add a digits manager for s-digits to the input list.
1374   //
1375
1376   fSDigitsManagerList->Add(man);
1377
1378 }
1379
1380 //_____________________________________________________________________________
1381 Bool_t AliTRDdigitizer::ConvertSDigits()
1382 {
1383   //
1384   // Converts s-digits to normal digits
1385   //
1386
1387   // Number of track dictionary arrays
1388   const Int_t    kNDict = AliTRDdigitsManager::kNDict;
1389
1390   // Converts number of electrons to fC
1391   const Double_t kEl2fC = 1.602E-19 * 1.0E15; 
1392
1393   Int_t iDict = 0;
1394
1395   if (fVerbose > 0) {
1396     this->Dump();
1397   }
1398
1399   Double_t sDigitsScale = 1.0 / GetSDigitsScale();
1400   Double_t noise        = GetNoise();
1401   Double_t padCoupling  = GetPadCoupling();
1402   Double_t timeCoupling = GetTimeCoupling();
1403   Double_t chipGain     = GetChipGain();
1404   Double_t convert      = kEl2fC * padCoupling * timeCoupling * chipGain;;
1405   Double_t adcInRange   = GetADCinRange();
1406   Double_t adcOutRange  = GetADCoutRange();
1407   Int_t    adcThreshold = GetADCthreshold();
1408
1409   AliTRDdataArrayI *digitsIn;
1410   AliTRDdataArrayI *digitsOut;
1411   AliTRDdataArrayI *dictionaryIn[kNDict];
1412   AliTRDdataArrayI *dictionaryOut[kNDict];
1413
1414   // Loop through the detectors
1415   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1416
1417     if (fVerbose > 0) {
1418       printf("AliTRDdigitizer::ConvertSDigits -- ");
1419       printf("Convert detector %d to digits.\n",iDet);
1420     }
1421
1422     Int_t plane      = fGeo->GetPlane(iDet);
1423     Int_t sector     = fGeo->GetSector(iDet);
1424     Int_t chamber    = fGeo->GetChamber(iDet);
1425     Int_t nRowMax    = fGeo->GetRowMax(plane,chamber,sector);
1426     Int_t nColMax    = fGeo->GetColMax(plane);
1427     Int_t nTimeTotal = fGeo->GetTimeTotal();
1428
1429     digitsIn  = fSDigitsManager->GetDigits(iDet);
1430     digitsIn->Expand();
1431     digitsOut = fDigitsManager->GetDigits(iDet);
1432     digitsOut->Allocate(nRowMax,nColMax,nTimeTotal);
1433     for (iDict = 0; iDict < kNDict; iDict++) {
1434       dictionaryIn[iDict]  = fSDigitsManager->GetDictionary(iDet,iDict);
1435       dictionaryIn[iDict]->Expand();
1436       dictionaryOut[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1437       dictionaryOut[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1438     }
1439
1440     for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1441       for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1442         for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {         
1443
1444           Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1445           signal *= sDigitsScale;
1446           // Add the noise
1447           signal  = TMath::Max((Double_t) gRandom->Gaus(signal,noise),0.0);
1448           // Convert to mV
1449           signal *= convert;
1450           // Convert to ADC counts. Set the overflow-bit adcOutRange if the 
1451           // signal is larger than adcInRange
1452           Int_t adc  = 0;
1453           if (signal >= adcInRange) {
1454             adc = ((Int_t) adcOutRange);
1455           }
1456           else {
1457             adc = ((Int_t) (signal * (adcOutRange / adcInRange)));
1458           }
1459           // Store the amplitude of the digit if above threshold
1460           if (adc > adcThreshold) {
1461             digitsOut->SetDataUnchecked(iRow,iCol,iTime,adc);
1462           }
1463           // Copy the dictionary
1464           for (iDict = 0; iDict < kNDict; iDict++) { 
1465             Int_t track = dictionaryIn[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1466             dictionaryOut[iDict]->SetDataUnchecked(iRow,iCol,iTime,track);
1467           }
1468
1469         }
1470       }
1471     }
1472
1473     if (fCompress) {
1474       digitsIn->Compress(1,0);
1475       digitsOut->Compress(1,0);
1476       for (iDict = 0; iDict < kNDict; iDict++) {
1477         dictionaryIn[iDict]->Compress(1,0);
1478         dictionaryOut[iDict]->Compress(1,0);
1479       }
1480     }
1481
1482   }    
1483
1484   return kTRUE;
1485
1486 }
1487
1488 //_____________________________________________________________________________
1489 Bool_t AliTRDdigitizer::MergeSDigits()
1490 {
1491   //
1492   // Merges the input s-digits:
1493   //   - The amplitude of the different inputs are summed up.
1494   //   - Of the track IDs from the input dictionaries only one is
1495   //     kept for each input. This works for maximal 3 different merged inputs.
1496   //
1497
1498   // Number of track dictionary arrays
1499   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1500
1501   Int_t iDict = 0;
1502
1503   AliTRDdataArrayI *digitsA;
1504   AliTRDdataArrayI *digitsB;
1505   AliTRDdataArrayI *dictionaryA[kNDict];
1506   AliTRDdataArrayI *dictionaryB[kNDict];
1507
1508   // Get the first s-digits
1509   fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
1510   if (!fSDigitsManager) return kFALSE;
1511
1512   // Loop through the other sets of s-digits
1513   AliTRDdigitsManager *mergeSDigitsManager;
1514   mergeSDigitsManager = (AliTRDdigitsManager *) 
1515                         fSDigitsManagerList->After(fSDigitsManager);
1516
1517   if (fVerbose > 0) {
1518     if (mergeSDigitsManager) {
1519       printf("AliTRDdigitizer::MergeSDigits -- ");
1520       printf("Merge serveral input files.\n");
1521     }
1522     else {
1523       printf("AliTRDdigitizer::MergeSDigits -- ");
1524       printf("Only one input file.\n");
1525     }
1526   }
1527
1528   Int_t iMerge = 0;
1529   while (mergeSDigitsManager) {
1530
1531     iMerge++;
1532
1533     // Loop through the detectors
1534     for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1535
1536       Int_t plane      = fGeo->GetPlane(iDet);
1537       Int_t sector     = fGeo->GetSector(iDet);
1538       Int_t chamber    = fGeo->GetChamber(iDet);
1539       Int_t nRowMax    = fGeo->GetRowMax(plane,chamber,sector);
1540       Int_t nColMax    = fGeo->GetColMax(plane);
1541       Int_t nTimeTotal = fGeo->GetTimeTotal();
1542
1543       // Loop through the pixels of one detector and add the signals
1544       digitsA = fSDigitsManager->GetDigits(iDet);
1545       digitsB = mergeSDigitsManager->GetDigits(iDet);
1546       digitsA->Expand();
1547       digitsB->Expand();
1548       for (iDict = 0; iDict < kNDict; iDict++) {
1549         dictionaryA[iDict] = fSDigitsManager->GetDictionary(iDet,iDict);
1550         dictionaryB[iDict] = mergeSDigitsManager->GetDictionary(iDet,iDict);
1551         dictionaryA[iDict]->Expand();
1552         dictionaryB[iDict]->Expand();
1553       }
1554
1555       if (fVerbose > 0) {
1556         printf("AliTRDdigitizer::MergeSDigits -- ");
1557         printf("Merge detector %d of input no.%d.\n",iDet,iMerge);
1558       }
1559
1560       for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1561         for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1562           for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {         
1563
1564             // Add the amplitudes of the summable digits 
1565             Int_t ampA = digitsA->GetDataUnchecked(iRow,iCol,iTime);
1566             Int_t ampB = digitsB->GetDataUnchecked(iRow,iCol,iTime);
1567             ampA += ampB;
1568             digitsA->SetDataUnchecked(iRow,iCol,iTime,ampA);
1569
1570             // Take only one track from each input
1571             Int_t track = dictionaryB[0]->GetDataUnchecked(iRow,iCol,iTime);
1572             if (iMerge < kNDict) {
1573               dictionaryA[iMerge]->SetDataUnchecked(iRow,iCol,iTime,track);
1574             }
1575
1576           }
1577         }
1578       }
1579
1580       if (fCompress) {
1581         digitsA->Compress(1,0);
1582         digitsB->Compress(1,0);
1583         for (iDict = 0; iDict < kNDict; iDict++) {
1584           dictionaryA[iDict]->Compress(1,0);
1585           dictionaryB[iDict]->Compress(1,0);
1586         }
1587       }
1588
1589     }    
1590
1591     // The next set of s-digits
1592     mergeSDigitsManager = (AliTRDdigitsManager *) 
1593                           fSDigitsManagerList->After(mergeSDigitsManager);
1594
1595   }
1596
1597   return kTRUE;
1598
1599 }
1600
1601 //_____________________________________________________________________________
1602 Bool_t AliTRDdigitizer::SDigits2Digits()
1603 {
1604   //
1605   // Merges the input s-digits and converts them to normal digits
1606   //
1607
1608   if (!MergeSDigits()) return kFALSE;
1609
1610   return ConvertSDigits();
1611
1612 }
1613
1614 //_____________________________________________________________________________
1615 Bool_t AliTRDdigitizer::CheckDetector(Int_t plane, Int_t chamber, Int_t sector)
1616 {
1617   //
1618   // Checks whether a detector is enabled
1619   //
1620
1621   if ((fTRD->GetSensChamber() >=       0) &&
1622       (fTRD->GetSensChamber() != chamber)) return kFALSE;
1623   if ((fTRD->GetSensPlane()   >=       0) &&
1624       (fTRD->GetSensPlane()   !=   plane)) return kFALSE;
1625   if ( fTRD->GetSensSector()  >=       0) {
1626     Int_t sens1 = fTRD->GetSensSector();
1627     Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
1628     sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
1629            * AliTRDgeometry::Nsect();
1630     if (sens1 < sens2) {
1631       if ((sector < sens1) || (sector >= sens2)) return kFALSE;
1632     }
1633     else {
1634       if ((sector < sens1) && (sector >= sens2)) return kFALSE;
1635     }
1636   }
1637
1638   return kTRUE;
1639
1640 }
1641
1642 //_____________________________________________________________________________
1643 Bool_t AliTRDdigitizer::WriteDigits()
1644 {
1645   //
1646   // Writes out the TRD-digits and the dictionaries
1647   //
1648
1649   // Store the digits and the dictionary in the tree
1650   return fDigitsManager->WriteDigits();
1651
1652 }
1653
1654 //_____________________________________________________________________________
1655 Float_t AliTRDdigitizer::GetDiffusionL(Float_t vd, Float_t b)
1656 {
1657   //
1658   // Returns the longitudinal diffusion coefficient for a given drift 
1659   // velocity <vd> and a B-field <b> for Xe/CO2 (15%).
1660   // The values are according to a GARFIELD simulation.
1661   //
1662
1663   const Int_t kNb = 5;
1664   Float_t p0[kNb] = {  0.007440,  0.007493,  0.007513,  0.007672,  0.007831 };
1665   Float_t p1[kNb] = {  0.019252,  0.018912,  0.018636,  0.018012,  0.017343 };
1666   Float_t p2[kNb] = { -0.005042, -0.004926, -0.004867, -0.004650, -0.004424 };
1667   Float_t p3[kNb] = {  0.000195,  0.000189,  0.000195,  0.000182,  0.000169 };
1668
1669   Int_t ib = ((Int_t) (10 * (b - 0.15)));
1670   ib       = TMath::Max(  0,ib);
1671   ib       = TMath::Min(kNb,ib);
1672
1673   Float_t diff = p0[ib] 
1674                + p1[ib] * vd
1675                + p2[ib] * vd*vd
1676                + p3[ib] * vd*vd*vd;
1677
1678   return diff;
1679
1680 }
1681
1682 //_____________________________________________________________________________
1683 Float_t AliTRDdigitizer::GetDiffusionT(Float_t vd, Float_t b)
1684 {
1685   //
1686   // Returns the transverse diffusion coefficient for a given drift 
1687   // velocity <vd> and a B-field <b> for Xe/CO2 (15%).
1688   // The values are according to a GARFIELD simulation.
1689   //
1690
1691   const Int_t kNb = 5;
1692   Float_t p0[kNb] = {  0.009550,  0.009599,  0.009674,  0.009757,  0.009850 };
1693   Float_t p1[kNb] = {  0.006667,  0.006539,  0.006359,  0.006153,  0.005925 };
1694   Float_t p2[kNb] = { -0.000853, -0.000798, -0.000721, -0.000635, -0.000541 };
1695   Float_t p3[kNb] = {  0.000131,  0.000122,  0.000111,  0.000098,  0.000085 };
1696
1697   Int_t ib = ((Int_t) (10 * (b - 0.15)));
1698   ib       = TMath::Max(  0,ib);
1699   ib       = TMath::Min(kNb,ib);
1700
1701   Float_t diff = p0[ib] 
1702                + p1[ib] * vd
1703                + p2[ib] * vd*vd
1704                + p3[ib] * vd*vd*vd;
1705
1706   return diff;
1707
1708 }
1709
1710 //_____________________________________________________________________________
1711 Float_t AliTRDdigitizer::GetOmegaTau(Float_t vd, Float_t b)
1712 {
1713   //
1714   // Returns omega*tau (tan(Lorentz-angle)) for a given drift velocity <vd> 
1715   // and a B-field <b> for Xe/CO2 (15%).
1716   // The values are according to a GARFIELD simulation.
1717   //
1718
1719   const Int_t kNb = 5;
1720   Float_t p0[kNb] = {  0.004810,  0.007412,  0.010252,  0.013409,  0.016888 };
1721   Float_t p1[kNb] = {  0.054875,  0.081534,  0.107333,  0.131983,  0.155455 };
1722   Float_t p2[kNb] = { -0.008682, -0.012896, -0.016987, -0.020880, -0.024623 };
1723   Float_t p3[kNb] = {  0.000155,  0.000238,  0.000330,  0.000428,  0.000541 };
1724
1725   Int_t ib = ((Int_t) (10 * (b - 0.15)));
1726   ib       = TMath::Max(  0,ib);
1727   ib       = TMath::Min(kNb,ib);
1728
1729   Float_t alphaL = p0[ib] 
1730                  + p1[ib] * vd
1731                  + p2[ib] * vd*vd
1732                  + p3[ib] * vd*vd*vd;
1733
1734   return TMath::Tan(alphaL);
1735
1736 }
1737