8a195b2f95c56792715d25e4575db31020660143
[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.22  2001/03/30 14:40:14  cblume
19 Update of the digitization parameter
20
21 Revision 1.21  2001/03/13 09:30:35  cblume
22 Update of digitization. Moved digit branch definition to AliTRD
23
24 Revision 1.20  2001/02/25 20:19:00  hristov
25 Minor correction: loop variable declared only once for HP, Sun
26
27 Revision 1.19  2001/02/14 18:22:26  cblume
28 Change in the geometry of the padplane
29
30 Revision 1.18  2001/01/26 19:56:57  hristov
31 Major upgrade of AliRoot code
32
33 Revision 1.17  2000/12/08 12:53:27  cblume
34 Change in Copy() function for HP-compiler
35
36 Revision 1.16  2000/12/07 12:20:46  cblume
37 Go back to array compression. Use sampled PRF to speed up digitization
38
39 Revision 1.15  2000/11/23 14:34:08  cblume
40 Fixed bug in expansion routine of arrays (initialize buffers properly)
41
42 Revision 1.14  2000/11/20 08:54:44  cblume
43 Switch off compression as default
44
45 Revision 1.13  2000/11/10 14:57:52  cblume
46 Changes in the geometry constants for the DEC compiler
47
48 Revision 1.12  2000/11/01 14:53:20  cblume
49 Merge with TRD-develop
50
51 Revision 1.1.4.9  2000/10/26 17:00:22  cblume
52 Fixed bug in CheckDetector()
53
54 Revision 1.1.4.8  2000/10/23 13:41:35  cblume
55 Added protection against Log(0) in the gas gain calulation
56
57 Revision 1.1.4.7  2000/10/17 02:27:34  cblume
58 Get rid of global constants
59
60 Revision 1.1.4.6  2000/10/16 01:16:53  cblume
61 Changed timebin 0 to be the one closest to the readout
62
63 Revision 1.1.4.5  2000/10/15 23:34:29  cblume
64 Faster version of the digitizer
65
66 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
67 Made Getters const
68
69 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
70 Replace include files by forward declarations
71
72 Revision 1.1.4.2  2000/09/22 14:41:10  cblume
73 Bug fix in PRF. Included time response. New structure
74
75 Revision 1.10  2000/10/05 07:27:53  cblume
76 Changes in the header-files by FCA
77
78 Revision 1.9  2000/10/02 21:28:19  fca
79 Removal of useless dependecies via forward declarations
80
81 Revision 1.8  2000/06/09 11:10:07  cblume
82 Compiler warnings and coding conventions, next round
83
84 Revision 1.7  2000/06/08 18:32:58  cblume
85 Make code compliant to coding conventions
86
87 Revision 1.6  2000/06/07 16:27:32  cblume
88 Try to remove compiler warnings on Sun and HP
89
90 Revision 1.5  2000/05/09 16:38:57  cblume
91 Removed PadResponse(). Merge problem
92
93 Revision 1.4  2000/05/08 15:53:45  cblume
94 Resolved merge conflict
95
96 Revision 1.3  2000/04/28 14:49:27  cblume
97 Only one declaration of iDict in MakeDigits()
98
99 Revision 1.1.4.1  2000/05/08 14:42:04  cblume
100 Introduced AliTRDdigitsManager
101
102 Revision 1.1  2000/02/28 19:00:13  cblume
103 Add new TRD classes
104
105 */
106
107 ///////////////////////////////////////////////////////////////////////////////
108 //                                                                           //
109 //  Creates and handles digits from TRD hits                                 //
110 //                                                                           //
111 //  The following effects are included:                                      //
112 //      - Diffusion                                                          //
113 //      - ExB effects                                                        //
114 //      - Gas gain including fluctuations                                    //
115 //      - Pad-response (simple Gaussian approximation)                       //
116 //      - Electronics noise                                                  //
117 //      - Electronics gain                                                   //
118 //      - Digitization                                                       //
119 //      - ADC threshold                                                      //
120 //  The corresponding parameter can be adjusted via the various              //
121 //  Set-functions. If these parameters are not explicitly set, default       //
122 //  values are used (see Init-function).                                     //
123 //  To produce digits from a root-file with TRD-hits use the                 //
124 //  slowDigitsCreate.C macro.                                                //
125 //                                                                           //
126 ///////////////////////////////////////////////////////////////////////////////
127
128 #include <stdlib.h>
129
130 #include <TMath.h>
131 #include <TVector.h>
132 #include <TRandom.h>
133 #include <TROOT.h>
134 #include <TTree.h>
135 #include <TFile.h>
136 #include <TF1.h>
137
138 #include "AliRun.h"
139
140 #include "AliTRD.h"
141 #include "AliTRDhit.h"
142 #include "AliTRDdigitizer.h"
143 #include "AliTRDdataArrayI.h"
144 #include "AliTRDdataArrayF.h"
145 #include "AliTRDsegmentArray.h"
146 #include "AliTRDdigitsManager.h"
147 #include "AliTRDgeometry.h"
148
149 ClassImp(AliTRDdigitizer)
150
151 //_____________________________________________________________________________
152 AliTRDdigitizer::AliTRDdigitizer():TNamed()
153 {
154   //
155   // AliTRDdigitizer default constructor
156   //
157
158   fInputFile      = NULL;
159   fDigits         = NULL;
160   fTRD            = NULL;
161   fGeo            = NULL;
162   fPRF            = NULL;
163   fPRFsmp         = NULL;
164   fTRFsmp         = NULL;
165
166   fEvent          = 0;
167   fGasGain        = 0.0;
168   fNoise          = 0.0;
169   fChipGain       = 0.0;
170   fSinRange       = 0.0;
171   fSoutRange      = 0.0;
172   fADCoutRange    = 0.0;
173   fADCinRange     = 0.0;
174   fADCthreshold   = 0;
175   fDiffusionOn    = 0;
176   fDiffusionT     = 0.0;
177   fDiffusionL     = 0.0;
178   fElAttachOn     = 0;
179   fElAttachProp   = 0.0;
180   fExBOn          = 0;
181   fOmegaTau       = 0.0;
182   fPRFOn          = 0;
183   fTRFOn          = 0;
184   fDriftVelocity  = 0.0;
185   fPadCoupling    = 0.0;
186   fTimeCoupling   = 0.0;
187   fTimeBinWidth   = 0.0;
188
189   fPRFbin         = 0;
190   fPRFlo          = 0.0;
191   fPRFhi          = 0.0;
192   fPRFwid         = 0.0;
193   fPRFpad         = 0;
194   fTRFbin         = 0;
195   fTRFlo          = 0.0;
196   fTRFhi          = 0.0;
197   fTRFwid         = 0.0;
198
199   fCompress       = kTRUE;
200   fVerbose        = 1;
201   fSDigits        = kFALSE;
202
203 }
204
205 //_____________________________________________________________________________
206 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
207                 :TNamed(name,title)
208 {
209   //
210   // AliTRDdigitizer default constructor
211   //
212
213   fInputFile     = NULL;
214   fDigits        = NULL;
215   fTRD           = NULL;
216   fGeo           = NULL;
217   fPRF           = NULL;
218   fPRFsmp        = NULL;
219   fTRFsmp        = NULL;
220
221   fEvent         = 0;
222
223   fCompress      = kTRUE;
224   fVerbose       = 1;
225   fSDigits       = kFALSE;
226
227   Init();
228
229 }
230
231 //_____________________________________________________________________________
232 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
233 {
234   //
235   // AliTRDdigitizer copy constructor
236   //
237
238   ((AliTRDdigitizer &) d).Copy(*this);
239
240 }
241
242 //_____________________________________________________________________________
243 AliTRDdigitizer::~AliTRDdigitizer()
244 {
245   //
246   // AliTRDdigitizer destructor
247   //
248
249   if (fInputFile) {
250     fInputFile->Close();
251     delete fInputFile;
252   }
253
254   if (fDigits) {
255     delete fDigits;
256   }
257
258   if (fPRF) delete fPRF;
259
260 }
261
262 //_____________________________________________________________________________
263 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
264 {
265   //
266   // Assignment operator
267   //
268
269   if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
270   return *this;
271
272 }
273
274 //_____________________________________________________________________________
275 void AliTRDdigitizer::Copy(TObject &d)
276 {
277   //
278   // Copy function
279   //
280
281   Int_t iBin;
282
283   ((AliTRDdigitizer &) d).fInputFile      = NULL;
284   ((AliTRDdigitizer &) d).fDigits         = NULL;
285   ((AliTRDdigitizer &) d).fTRD            = NULL;
286   ((AliTRDdigitizer &) d).fGeo            = NULL;
287
288   ((AliTRDdigitizer &) d).fEvent          = 0;
289
290   ((AliTRDdigitizer &) d).fGasGain        = fGasGain;
291   ((AliTRDdigitizer &) d).fNoise          = fNoise;
292   ((AliTRDdigitizer &) d).fChipGain       = fChipGain;
293   ((AliTRDdigitizer &) d).fSoutRange      = fSoutRange;
294   ((AliTRDdigitizer &) d).fSinRange       = fSinRange;
295   ((AliTRDdigitizer &) d).fADCoutRange    = fADCoutRange;
296   ((AliTRDdigitizer &) d).fADCinRange     = fADCinRange;
297   ((AliTRDdigitizer &) d).fADCthreshold   = fADCthreshold;
298   ((AliTRDdigitizer &) d).fDiffusionOn    = fDiffusionOn; 
299   ((AliTRDdigitizer &) d).fDiffusionT     = fDiffusionT;
300   ((AliTRDdigitizer &) d).fDiffusionL     = fDiffusionL;
301   ((AliTRDdigitizer &) d).fElAttachOn     = fElAttachOn;
302   ((AliTRDdigitizer &) d).fElAttachProp   = fElAttachProp;
303   ((AliTRDdigitizer &) d).fExBOn          = fExBOn;
304   ((AliTRDdigitizer &) d).fOmegaTau       = fOmegaTau;
305   ((AliTRDdigitizer &) d).fLorentzFactor  = fLorentzFactor;
306   ((AliTRDdigitizer &) d).fDriftVelocity  = fDriftVelocity;
307   ((AliTRDdigitizer &) d).fPadCoupling    = fPadCoupling;
308   ((AliTRDdigitizer &) d).fTimeCoupling   = fTimeCoupling;
309   ((AliTRDdigitizer &) d).fTimeBinWidth   = fTimeBinWidth;
310   ((AliTRDdigitizer &) d).fPRFOn          = fPRFOn;
311   ((AliTRDdigitizer &) d).fTRFOn          = fTRFOn;
312
313   ((AliTRDdigitizer &) d).fCompress       = fCompress;
314   ((AliTRDdigitizer &) d).fVerbose        = fVerbose;
315   ((AliTRDdigitizer &) d).fSDigits        = fSDigits;
316
317   fPRF->Copy(*((AliTRDdigitizer &) d).fPRF);
318
319   ((AliTRDdigitizer &) d).fPRFbin         = fPRFbin;
320   ((AliTRDdigitizer &) d).fPRFlo          = fPRFlo;
321   ((AliTRDdigitizer &) d).fPRFhi          = fPRFhi;
322   ((AliTRDdigitizer &) d).fPRFwid         = fPRFwid;
323   ((AliTRDdigitizer &) d).fPRFpad         = fPRFpad;
324   if (((AliTRDdigitizer &) d).fPRFsmp) delete ((AliTRDdigitizer &) d).fPRFsmp;
325   ((AliTRDdigitizer &) d).fPRFsmp = new Float_t[fPRFbin];
326   for (iBin = 0; iBin < fPRFbin; iBin++) {
327     ((AliTRDdigitizer &) d).fPRFsmp[iBin] = fPRFsmp[iBin];
328   }                                                                             
329   ((AliTRDdigitizer &) d).fTRFbin         = fTRFbin;
330   ((AliTRDdigitizer &) d).fTRFlo          = fTRFlo;
331   ((AliTRDdigitizer &) d).fTRFhi          = fTRFhi;
332   ((AliTRDdigitizer &) d).fTRFwid         = fTRFwid;
333   if (((AliTRDdigitizer &) d).fTRFsmp) delete ((AliTRDdigitizer &) d).fTRFsmp;
334   ((AliTRDdigitizer &) d).fTRFsmp = new Float_t[fTRFbin];
335   for (iBin = 0; iBin < fTRFbin; iBin++) {
336     ((AliTRDdigitizer &) d).fTRFsmp[iBin] = fTRFsmp[iBin];
337   }                                      
338                                        
339 }
340
341 //_____________________________________________________________________________
342 Int_t AliTRDdigitizer::Diffusion(Float_t driftlength, Float_t *xyz)
343 {
344   //
345   // Applies the diffusion smearing to the position of a single electron
346   //
347
348   Float_t driftSqrt = TMath::Sqrt(driftlength);
349   Float_t sigmaT = driftSqrt * fDiffusionT;
350   Float_t sigmaL = driftSqrt * fDiffusionL;
351   xyz[0] = gRandom->Gaus(xyz[0], sigmaL * fLorentzFactor);
352   xyz[1] = gRandom->Gaus(xyz[1], sigmaT * fLorentzFactor);
353   xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
354
355   return 1;
356
357 }
358
359 //_____________________________________________________________________________
360 Int_t AliTRDdigitizer::ExB(Float_t driftlength, Float_t *xyz)
361 {
362   //
363   // Applies E x B effects to the position of a single electron
364   //
365
366   xyz[0] = xyz[0];
367   xyz[1] = xyz[1] + fOmegaTau * driftlength;
368   xyz[2] = xyz[2];
369
370   return 1;
371
372 }
373
374 //_____________________________________________________________________________
375 Int_t AliTRDdigitizer::PadResponse(Float_t signal, Float_t dist, Float_t *pad)
376 {
377   //
378   // Applies the pad response
379   //
380
381   Int_t iBin =  ((Int_t) (( - dist - fPRFlo) / fPRFwid));
382
383   Int_t iBin0 = iBin - fPRFpad;
384   Int_t iBin1 = iBin;
385   Int_t iBin2 = iBin + fPRFpad;
386
387   if ((iBin0 >= 0) && (iBin2 < fPRFbin)) {
388
389     pad[0] = signal * fPRFsmp[iBin0];
390     pad[1] = signal * fPRFsmp[iBin1];
391     pad[2] = signal * fPRFsmp[iBin2];
392
393     return 1;
394
395   }
396   else {
397
398     return 0;
399
400   }
401
402 }
403
404 //_____________________________________________________________________________
405 Float_t AliTRDdigitizer::TimeResponse(Float_t time)
406 {
407   //
408   // Applies the preamp shaper time response
409   //
410
411   Int_t iBin = ((Int_t) ((time - fTRFlo) / fTRFwid)); 
412   if ((iBin >= 0) && (iBin < fTRFbin)) {
413     return fTRFsmp[iBin];
414   }
415   else {
416     return 0.0;
417   }    
418
419 }
420
421 //_____________________________________________________________________________
422 void AliTRDdigitizer::Init()
423 {
424   //
425   // Initializes the digitization procedure with standard values
426   //
427
428   // The default parameter for the digitization
429   // Updated to new TRF 200 ns   
430   fGasGain        = 1600.;
431   fChipGain       = 8.0;
432   fNoise          = 1000.;
433   fADCoutRange    = 1023.;          // 10-bit ADC
434   fADCinRange     = 1000.;          // 1V input range
435   fADCthreshold   = 1;
436
437   // For the summable digits
438   fSinRange       = 1000000.;
439   fSoutRange      = 1000000.;
440
441   // Transverse and longitudinal diffusion coefficients (Xe/Isobutane)
442   fDiffusionOn    = 1;
443   fDiffusionT     = 0.060;
444   fDiffusionL     = 0.017;
445
446   // Propability for electron attachment
447   fElAttachOn     = 0;
448   fElAttachProp   = 0.0;
449
450   // E x B effects
451   fExBOn          = 0;
452   // omega * tau.= arctan(Lorentz-angle)
453   fOmegaTau       = 0.19438031;
454
455   // The pad response function
456   fPRFOn          =  1;
457   fPRFlo          = -3.0;
458   fPRFhi          =  3.0;
459   fPRFbin         = 120;
460   fPRFwid         = (fPRFhi - fPRFlo) / ((Float_t) fPRFbin);
461   fPRFpad         = ((Int_t) (1.0 / fPRFwid));
462   // New PRF from Bogdan  25/04/01
463   fPRF            = new TF1("PRF"
464                            ,"[0]*([1]+exp(-pow(sqrt(x*x),[3])/(2.0*[2])))"
465                            ,fPRFlo,fPRFhi);
466   fPRF->SetParameter(0, 0.8303); 
467   fPRF->SetParameter(1,-0.00392); 
468   fPRF->SetParameter(2, 0.472 * 0.472); 
469   fPRF->SetParameter(3, 2.19); 
470
471   // The time response function
472   fTRFOn          =  1;
473
474   // The drift velocity (cm / mus)
475   fDriftVelocity  = 2.0;
476
477   // The pad coupling factor (same number as for the TPC)
478   fPadCoupling    = 0.5;
479
480   // The time coupling factor (same number as for the TPC)
481   fTimeCoupling   = 0.4;
482
483 }
484
485 //_____________________________________________________________________________
486 void AliTRDdigitizer::ReInit()
487 {
488   //
489   // Reinitializes the digitization procedure after a change in the parameter
490   //
491
492   if (!fGeo) {
493     printf("AliTRDdigitizer::ReInit -- ");
494     printf("No geometry defined. Run InitDetector() first\n");
495     exit(1);
496   }
497
498   // Calculate the time bin width in ns
499   fTimeBinWidth   = fGeo->GetTimeBinSize() / fDriftVelocity * 1000.0;
500
501   // The range and the binwidth for the sampled TRF 
502   fTRFbin = 100;
503   // Start 0.8 mus before the signal
504   fTRFlo  = -0.8 * fDriftVelocity;
505   // End the maximum driftlength after the signal 
506   fTRFhi  = AliTRDgeometry::DrThick() 
507           + fGeo->GetTimeAfter() * fGeo->GetTimeBinSize();
508   fTRFwid = (fTRFhi - fTRFlo) / ((Float_t) fTRFbin);
509
510   // The Lorentz factor
511   if (fExBOn) {
512     fLorentzFactor = 1.0 / (1.0 + fOmegaTau*fOmegaTau);
513   }
514   else {
515     fLorentzFactor = 1.0;
516   }
517
518 }
519
520 //_____________________________________________________________________________
521 void AliTRDdigitizer::SampleTRF()
522 {
523   //
524   // Samples the time response function
525   // It is defined according to Vasiles simulation of the preamp shaper
526   // output and includes the effect of the ion tail (based on Tariqs 
527   // Garfield simulation) and a shaping time of 125 ns FWHM
528   //
529
530   Int_t   ipos1;
531   Int_t   ipos2;
532   Float_t diff;
533
534   const Float_t kShift = 0.0;
535   const Float_t kScale = 0.5;
536   //const Float_t kScale = 1.0;
537
538   const Int_t kNpasa = 36;
539   Float_t time[kNpasa]   = { -2.80,     -2.40,     -2.00,     -1.60
540                            , -1.20,     -0.80,     -0.60,     -0.40
541                            , -0.30,     -0.20,     -0.10,      0.00
542                            ,  0.10,      0.20,      0.30,      0.40
543                            ,  0.60,      0.80,      1.20,      1.60
544                            ,  2.00,      2.40,      2.80,      3.20
545                            ,  3.60,      4.00,      4.40,      4.80
546                            ,  5.20,      5.60,      7.20,      9.20
547                            , 11.20,     13.20,     15.20,     17.20     };
548   Float_t signal[kNpasa] = {  0.000000,  0.000000,  0.000000,  0.000000
549                            ,  0.000000,  0.000000,  0.015385,  0.086154
550                            ,  0.236923,  0.452308,  0.726154,  1.003077
551                            ,  0.953846,  0.652307,  0.332308,  0.181539
552                            ,  0.120000,  0.083077,  0.049231,  0.024615
553                            ,  0.015385,  0.009231,  0.003077,  0.000000
554                            , -0.003077, -0.006154, -0.009231, -0.012308
555                            , -0.015385, -0.018462, -0.018462, -0.018462
556                            , -0.015385, -0.012308, -0.009231, -0.006154 };
557   for (Int_t ipasa = 0; ipasa < kNpasa; ipasa++) {
558     time[ipasa] = kScale * time[ipasa] + kShift; 
559   }
560
561   if (fTRFsmp) delete fTRFsmp;
562   fTRFsmp = new Float_t[fTRFbin];
563
564   Float_t loTRF    = TMath::Max(fTRFlo / fDriftVelocity,time[0]);
565   Float_t hiTRF    = TMath::Min(fTRFhi / fDriftVelocity,time[kNpasa-1]);
566   Float_t binWidth = (hiTRF - loTRF) / ((Float_t) fTRFbin);
567
568   // Take the linear interpolation
569   for (Int_t iBin = 0; iBin < fTRFbin; iBin++) {
570
571     Float_t bin = (((Float_t) iBin) + 0.5) * binWidth + loTRF;
572     ipos1 = ipos2 = 0;
573     diff  = 0;
574     do {
575       diff = bin - time[ipos2++];
576     } while (diff > 0);
577     ipos2--;
578     if (ipos2 > kNpasa) ipos2 = kNpasa - 1;
579     ipos1 = ipos2 - 1;
580
581     fTRFsmp[iBin] = signal[ipos2] 
582                   + diff * (signal[ipos2] - signal[ipos1]) 
583                          / (  time[ipos2] -   time[ipos1]);
584
585   }
586
587 }
588
589 //_____________________________________________________________________________
590 void AliTRDdigitizer::SamplePRF()
591 {
592   //
593   // Samples the pad response function
594   //
595
596   if (fPRFsmp) delete fPRFsmp;
597   fPRFsmp = new Float_t[fPRFbin];
598   for (Int_t iBin = 0; iBin < fPRFbin; iBin++) {
599     Float_t bin = (((Float_t ) iBin) + 0.5) * fPRFwid + fPRFlo;
600     fPRFsmp[iBin] = TMath::Max(fPRF->Eval(bin),0.0);
601   }
602
603 }
604
605 //_____________________________________________________________________________
606 Bool_t AliTRDdigitizer::Open(const Char_t *name, Int_t nEvent)
607 {
608   //
609   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
610   //
611
612   // Connect the AliRoot file containing Geometry, Kine, and Hits
613   fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
614   if (!fInputFile) {
615     printf("AliTRDdigitizer::Open -- ");
616     printf("Open the ALIROOT-file %s.\n",name);
617     fInputFile = new TFile(name,"UPDATE");
618   }
619   else {
620     printf("AliTRDdigitizer::Open -- ");
621     printf("%s is already open.\n",name);
622   }
623
624   gAlice = (AliRun*) fInputFile->Get("gAlice");
625   if (gAlice) {
626     printf("AliTRDdigitizer::Open -- ");
627     printf("AliRun object found on file.\n");
628   }
629   else {
630     printf("AliTRDdigitizer::Open -- ");
631     printf("Could not find AliRun object.\n");
632     return kFALSE;
633   }
634
635   fEvent = nEvent;
636
637   // Import the Trees for the event nEvent in the file
638   Int_t nparticles = gAlice->GetEvent(fEvent);
639   if (nparticles <= 0) {
640     printf("AliTRDdigitizer::Open -- ");
641     printf("No entries in the trees for event %d.\n",fEvent);
642     return kFALSE;
643   }
644
645   return InitDetector();
646
647 }
648
649 //_____________________________________________________________________________
650 Bool_t AliTRDdigitizer::InitDetector()
651 {
652   //
653   // Sets the pointer to the TRD detector and the geometry
654   //
655
656   // Get the pointer to the detector class and check for version 1
657   fTRD = (AliTRD*) gAlice->GetDetector("TRD");
658   if (fTRD->IsVersion() != 1) {
659     printf("AliTRDdigitizer::InitDetector -- ");
660     printf("TRD must be version 1 (slow simulator).\n");
661     exit(1);
662   }
663
664   // Get the geometry
665   fGeo = fTRD->GetGeometry();
666   printf("AliTRDdigitizer::InitDetector -- ");
667   printf("Geometry version %d\n",fGeo->IsVersion());
668
669   ReInit();
670
671   return kTRUE;
672
673 }
674
675 //_____________________________________________________________________________
676 Bool_t AliTRDdigitizer::SumSDigits()
677 {
678   //
679   // Sums up the summable digits and creates final digits
680   // Not yet implemented
681   //
682
683   return kFALSE;
684
685 }
686
687 //_____________________________________________________________________________
688 Bool_t AliTRDdigitizer::MakeDigits()
689 {
690   //
691   // Creates digits.
692   //
693
694   ///////////////////////////////////////////////////////////////
695   // Parameter 
696   ///////////////////////////////////////////////////////////////
697
698   // Converts number of electrons to fC
699   const Double_t kEl2fC  = 1.602E-19 * 1.0E15; 
700
701   ///////////////////////////////////////////////////////////////
702
703   // Number of pads included in the pad response
704   const Int_t kNpad  = 3;
705
706   // Number of track dictionary arrays
707   const Int_t kNDict = AliTRDdigitsManager::kNDict;
708
709   // Half the width of the amplification region
710   const Float_t kAmWidth = AliTRDgeometry::AmThick() / 2.;
711
712   Int_t   iRow, iCol, iTime, iPad;
713   Int_t   iDict  = 0;
714   Int_t   nBytes = 0;
715
716   Int_t   totalSizeDigits = 0;
717   Int_t   totalSizeDict0  = 0;
718   Int_t   totalSizeDict1  = 0;
719   Int_t   totalSizeDict2  = 0;
720
721   Int_t   timeTRDbeg = 0;
722   Int_t   timeTRDend = 1;
723
724   Float_t pos[3];
725   Float_t rot[3];
726   Float_t xyz[3];
727   Float_t padSignal[kNpad];
728   Float_t signalOld[kNpad];
729
730   AliTRDdataArrayF *signals = 0;
731   AliTRDdataArrayI *digits  = 0;
732   AliTRDdataArrayI *dictionary[kNDict];
733
734   // Create a digits manager
735   fDigits = new AliTRDdigitsManager();
736
737   // Create a container for the amplitudes
738   AliTRDsegmentArray *signalsArray 
739                      = new AliTRDsegmentArray("AliTRDdataArrayF",AliTRDgeometry::Ndet());
740
741   if (fTRFOn) {
742     timeTRDbeg = ((Int_t) (-fTRFlo / fGeo->GetTimeBinSize())) - 1;
743     timeTRDend = ((Int_t) ( fTRFhi / fGeo->GetTimeBinSize())) - 1;
744     printf("AliTRDdigitizer::MakeDigits -- ");
745     printf("Sample the TRF between -%d and %d\n",timeTRDbeg,timeTRDend);
746   }
747
748   Float_t elAttachProp = fElAttachProp / 100.; 
749
750   // Create the sampled PRF
751   SamplePRF();
752
753   // Create the sampled TRF
754   SampleTRF();
755
756   if (!fGeo) {
757     printf("AliTRDdigitizer::MakeDigits -- ");
758     printf("No geometry defined\n");
759     return kFALSE;
760   }
761
762   printf("AliTRDdigitizer::MakeDigits -- ");
763   printf("Start creating digits.\n");
764   if (fVerbose > 0) this->Dump();
765
766   // Get the pointer to the hit tree
767   TTree *HitTree = gAlice->TreeH();
768
769   // Get the number of entries in the hit tree
770   // (Number of primary particles creating a hit somewhere)
771   Int_t nTrack = (Int_t) HitTree->GetEntries();
772   if (fVerbose > 0) {
773     printf("AliTRDdigitizer::MakeDigits -- ");
774     printf("Found %d primary particles\n",nTrack);
775   } 
776
777   Int_t detectorOld = -1;
778   Int_t countHits   =  0; 
779
780   // Loop through all entries in the tree
781   for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
782
783     gAlice->ResetHits();
784     nBytes += HitTree->GetEvent(iTrack);
785
786     // Get the number of hits in the TRD created by this particle
787     Int_t nHit = fTRD->Hits()->GetEntriesFast();
788     if (fVerbose > 0) {
789       printf("AliTRDdigitizer::MakeDigits -- ");
790       printf("Found %d hits for primary particle %d\n",nHit,iTrack);
791     }
792
793     // Loop through the TRD hits  
794     for (Int_t iHit = 0; iHit < nHit; iHit++) {
795
796       countHits++;
797
798       AliTRDhit *hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit);
799               pos[0]   = hit->X();
800               pos[1]   = hit->Y();
801               pos[2]   = hit->Z();
802       Float_t q        = hit->GetCharge();
803       Int_t   track    = hit->Track();
804       Int_t   detector = hit->GetDetector();
805       Int_t   plane    = fGeo->GetPlane(detector);
806       Int_t   sector   = fGeo->GetSector(detector);
807       Int_t   chamber  = fGeo->GetChamber(detector);
808
809       if (!(CheckDetector(plane,chamber,sector))) continue;
810
811       Int_t   nRowMax     = fGeo->GetRowMax(plane,chamber,sector);
812       Int_t   nColMax     = fGeo->GetColMax(plane);
813       Int_t   nTimeMax    = fGeo->GetTimeMax();
814       Int_t   nTimeBefore = fGeo->GetTimeBefore();
815       Int_t   nTimeAfter  = fGeo->GetTimeAfter();
816       Int_t   nTimeTotal  = fGeo->GetTimeTotal();
817       Float_t row0        = fGeo->GetRow0(plane,chamber,sector);
818       Float_t col0        = fGeo->GetCol0(plane);
819       Float_t time0       = fGeo->GetTime0(plane);
820       Float_t rowPadSize  = fGeo->GetRowPadSize(plane,chamber,sector);
821       Float_t colPadSize  = fGeo->GetColPadSize(plane);
822       Float_t timeBinSize = fGeo->GetTimeBinSize();
823       Float_t divideRow   = 1.0 / rowPadSize;
824       Float_t divideCol   = 1.0 / colPadSize;
825       Float_t divideTime  = 1.0 / timeBinSize;
826
827       if (fVerbose > 1) {
828         printf("Analyze hit no. %d ",iHit);
829         printf("-----------------------------------------------------------\n");
830         hit->Dump();
831         printf("plane = %d, sector = %d, chamber = %d\n"
832               ,plane,sector,chamber);
833         printf("nRowMax = %d, nColMax = %d, nTimeMax = %d\n" 
834               ,nRowMax,nColMax,nTimeMax);
835         printf("nTimeBefore = %d, nTimeAfter = %d, nTimeTotal = %d\n"
836               ,nTimeBefore,nTimeAfter,nTimeTotal);
837         printf("row0 = %f, col0 = %f, time0 = %f\n"
838               ,row0,col0,time0);
839         printf("rowPadSize = %f, colPadSize = %f, timeBinSize = %f\n"
840                ,rowPadSize,colPadSize,timeBinSize); 
841       }
842        
843       // Don't analyze test hits
844       if (hit->FromTest()) continue;
845
846       if (detector != detectorOld) {
847
848         if (fVerbose > 1) {
849           printf("AliTRDdigitizer::MakeDigits -- ");
850           printf("Get new container. New det = %d, Old det = %d\n"
851                 ,detector,detectorOld);
852         }
853         // Compress the old one if enabled
854         if ((fCompress) && (detectorOld > -1)) {
855           if (fVerbose > 1) {
856             printf("AliTRDdigitizer::MakeDigits -- ");
857             printf("Compress the old container ...");
858           }
859           signals->Compress(1,0);
860           for (iDict = 0; iDict < kNDict; iDict++) {
861             dictionary[iDict]->Compress(1,0);
862           }
863           if (fVerbose > 1) printf("done\n");
864         }
865         // Get the new container
866         signals = (AliTRDdataArrayF *) signalsArray->At(detector);
867         if (signals->GetNtime() == 0) {
868           // Allocate a new one if not yet existing
869           if (fVerbose > 1) {
870             printf("AliTRDdigitizer::MakeDigits -- ");
871             printf("Allocate a new container ... ");
872           }
873           signals->Allocate(nRowMax,nColMax,nTimeTotal);
874         }
875         else {
876           // Expand an existing one
877           if (fCompress) {
878             if (fVerbose > 1) {
879               printf("AliTRDdigitizer::MakeDigits -- ");
880               printf("Expand an existing container ... ");
881             }
882             signals->Expand();
883           }
884         }
885         // The same for the dictionary
886         for (iDict = 0; iDict < kNDict; iDict++) {       
887           dictionary[iDict] = fDigits->GetDictionary(detector,iDict);
888           if (dictionary[iDict]->GetNtime() == 0) {
889             dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
890           }
891           else {
892             if (fCompress) dictionary[iDict]->Expand();
893           }
894         }      
895         if (fVerbose > 1) printf("done\n");
896         detectorOld = detector;
897       }
898
899       // Rotate the sectors on top of each other
900       fGeo->Rotate(detector,pos,rot);
901
902       // The driftlength. It is negative if the hit is in the 
903       // amplification region.
904       Float_t driftlength = time0 - rot[0];
905
906       // Take also the drift in the amplification region into account
907       // The drift length is at the moment still the same, regardless of
908       // the position relativ to the wire. This non-isochronity needs still
909       // to be implemented.
910       Float_t driftlengthL = TMath::Abs(driftlength + kAmWidth);
911       if (fExBOn) driftlengthL /= TMath::Sqrt(fLorentzFactor);
912
913       // Loop over all electrons of this hit
914       // TR photons produce hits with negative charge
915       Int_t nEl = ((Int_t) TMath::Abs(q));
916       for (Int_t iEl = 0; iEl < nEl; iEl++) {
917
918         xyz[0] = rot[0];
919         xyz[1] = rot[1];
920         xyz[2] = rot[2];
921
922         // Electron attachment
923         if (fElAttachOn) {
924           if (gRandom->Rndm() < (driftlengthL * elAttachProp)) 
925             continue;
926         }
927
928         // Apply the diffusion smearing
929         if (fDiffusionOn) {
930           if (!(Diffusion(driftlengthL,xyz))) continue;
931         }
932
933         // Apply E x B effects (depends on drift direction)
934         if (fExBOn) { 
935           if (!(ExB(driftlength+kAmWidth,xyz))) continue;   
936         }
937
938         // The electron position after diffusion and ExB in pad coordinates 
939         // The pad row (z-direction)
940         Int_t  rowE = ((Int_t) ((xyz[2] -  row0) * divideRow));
941         if ((rowE < 0) || (rowE >= nRowMax)) continue;   
942
943         // The pad column (rphi-direction)
944         Int_t  colE = ((Int_t) ((xyz[1] -  col0) * divideCol));
945         if ((colE < 0) || (colE >= nColMax)) continue;   
946   
947         // The time bin (negative for hits in the amplification region)
948         // In the amplification region the electrons drift from both sides
949         // to the middle (anode wire plane)
950         Float_t timeDist   = time0 - xyz[0];
951         Float_t timeOffset = 0;
952         Int_t   timeE      = 0;
953         if (timeDist > 0) {
954           // The time bin
955           timeE      = ((Int_t) (timeDist * divideTime));
956           // The distance of the position to the middle of the timebin
957           timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) - timeDist;
958         }
959         else {
960           // Difference between half of the amplification gap width and
961           // the distance to the anode wire
962           Float_t anodeDist = kAmWidth - TMath::Abs(timeDist + kAmWidth);
963           // The time bin
964           timeE      = -1 * (((Int_t ) (anodeDist * divideTime)) + 1);
965           // The distance of the position to the middle of the timebin
966           timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) + anodeDist;
967         }
968  
969         // Apply the gas gain including fluctuations
970         Float_t ggRndm = 0.0;
971         do {
972           ggRndm = gRandom->Rndm();
973         } while (ggRndm <= 0);
974         Int_t signal = (Int_t) (-fGasGain * TMath::Log(ggRndm));
975
976         // Apply the pad response 
977         if (fPRFOn) {
978           // The distance of the electron to the center of the pad 
979           // in units of pad width
980           Float_t dist = (xyz[1] - col0 - (colE + 0.5) * colPadSize) 
981                        * divideCol;
982           if (!(PadResponse(signal,dist,padSignal))) continue;
983         }
984         else {
985           padSignal[0] = 0.0;
986           padSignal[1] = signal;
987           padSignal[2] = 0.0;
988         }
989
990         // Sample the time response inside the drift region
991         // + additional time bins before and after.
992         // The sampling is done always in the middle of the time bin
993         for (Int_t iTimeBin = TMath::Max(timeE-timeTRDbeg,        -nTimeBefore) 
994                   ;iTimeBin < TMath::Min(timeE+timeTRDend,nTimeMax+nTimeAfter ) 
995                   ;iTimeBin++) {
996
997           // Apply the time response
998           Float_t timeResponse = 1.0;
999           if (fTRFOn) {
1000             Float_t time = (iTimeBin - timeE) * timeBinSize + timeOffset;
1001             timeResponse = TimeResponse(time);
1002           }
1003
1004           signalOld[0] = 0.0;
1005           signalOld[1] = 0.0;
1006           signalOld[2] = 0.0;
1007
1008           for (iPad = 0; iPad < kNpad; iPad++) {
1009
1010             Int_t colPos = colE + iPad - 1;
1011             if (colPos <        0) continue;
1012             if (colPos >= nColMax) break;
1013
1014             // Add the signals
1015             // Note: The time bin number is shifted by nTimeBefore to avoid negative
1016             // time bins. This has to be subtracted lateron.
1017             Int_t iCurrentTimeBin = iTimeBin + nTimeBefore;
1018             signalOld[iPad]  = signals->GetDataUnchecked(rowE,colPos,iCurrentTimeBin);
1019             signalOld[iPad] += padSignal[iPad] * timeResponse;
1020             signals->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,signalOld[iPad]);
1021
1022             // Store the track index in the dictionary
1023             // Note: We store index+1 in order to allow the array to be compressed
1024             if (signalOld[iPad] > 0) {
1025               for (iDict = 0; iDict < kNDict; iDict++) {
1026                 Int_t oldTrack = dictionary[iDict]->GetDataUnchecked(rowE
1027                                                                     ,colPos
1028                                                                     ,iCurrentTimeBin);
1029                 if (oldTrack == track+1) break;
1030                 if (oldTrack ==       0) {
1031                   dictionary[iDict]->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,track+1);
1032                   break;
1033                 }
1034               }
1035             }
1036
1037           }
1038
1039         }
1040
1041       }
1042
1043     }
1044
1045   } // All hits finished
1046
1047   printf("AliTRDdigitizer::MakeDigits -- ");
1048   printf("Finished analyzing %d hits\n",countHits);
1049
1050   // The total conversion factor
1051   Float_t convert = kEl2fC * fPadCoupling * fTimeCoupling * fChipGain;
1052
1053   // Loop through all chambers to finalize the digits
1054   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1055
1056     Int_t plane       = fGeo->GetPlane(iDet);
1057     Int_t sector      = fGeo->GetSector(iDet);
1058     Int_t chamber     = fGeo->GetChamber(iDet);
1059     Int_t nRowMax     = fGeo->GetRowMax(plane,chamber,sector);
1060     Int_t nColMax     = fGeo->GetColMax(plane);
1061     Int_t nTimeMax    = fGeo->GetTimeMax();
1062     Int_t nTimeTotal  = fGeo->GetTimeTotal();
1063
1064     if (fVerbose > 0) {
1065       printf("AliTRDdigitizer::MakeDigits -- ");
1066       printf("Digitization for chamber %d\n",iDet);
1067     }
1068
1069     // Add a container for the digits of this detector
1070     digits = fDigits->GetDigits(iDet);        
1071     // Allocate memory space for the digits buffer
1072     digits->Allocate(nRowMax,nColMax,nTimeTotal);
1073
1074     // Get the signal container
1075     signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
1076     if (signals->GetNtime() == 0) {
1077       // Create missing containers
1078       signals->Allocate(nRowMax,nColMax,nTimeTotal);      
1079     }
1080     else {
1081       // Expand the container if neccessary
1082       if (fCompress) signals->Expand();
1083     }
1084     // Create the missing dictionary containers
1085     for (iDict = 0; iDict < kNDict; iDict++) {       
1086       dictionary[iDict] = fDigits->GetDictionary(iDet,iDict);
1087       if (dictionary[iDict]->GetNtime() == 0) {
1088         dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1089       }
1090     }
1091
1092     Int_t nDigits = 0;
1093
1094     // Don't create noise in detectors that are switched off
1095     if (CheckDetector(plane,chamber,sector)) {
1096
1097       // Create the digits for this chamber
1098       for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1099         for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1100           for (iTime = 0; iTime < nTimeTotal; iTime++) {         
1101
1102             // Create summable digits
1103             if (fSDigits) {
1104
1105               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1106               Int_t adc  = 0;
1107               if (signalAmp >= fSinRange) {
1108                 adc = ((Int_t) fSoutRange);
1109               }
1110               else {
1111                 adc = ((Int_t) (signalAmp * (fSoutRange / fSinRange)));
1112               }
1113               nDigits++;
1114               digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1115
1116             }
1117             // Create normal digits
1118             else {
1119
1120               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1121
1122               // Add the noise
1123               signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,fNoise),0.0);
1124               // Convert to mV
1125               signalAmp *= convert;
1126               // Convert to ADC counts. Set the overflow-bit fADCoutRange if the 
1127               // signal is larger than fADCinRange
1128               Int_t adc  = 0;
1129               if (signalAmp >= fADCinRange) {
1130                 adc = ((Int_t) fADCoutRange);
1131               }
1132               else {
1133                 adc = ((Int_t) (signalAmp * (fADCoutRange / fADCinRange)));
1134               }
1135
1136               // Store the amplitude of the digit if above threshold
1137               if (adc > fADCthreshold) {
1138                 if (fVerbose > 2) {
1139                   printf("  iRow = %d, iCol = %d, iTime = %d\n"
1140                         ,iRow,iCol,iTime);
1141                   printf("  signal = %f, adc = %d\n",signalAmp,adc);
1142                 }
1143                 nDigits++;
1144                 digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1145               }
1146
1147             }
1148
1149           }
1150         }
1151       }
1152
1153     }
1154
1155     // Compress the arrays
1156     digits->Compress(1,0);
1157     for (iDict = 0; iDict < kNDict; iDict++) {
1158       dictionary[iDict]->Compress(1,0);
1159     }
1160
1161     totalSizeDigits += digits->GetSize();
1162     totalSizeDict0  += dictionary[0]->GetSize();
1163     totalSizeDict1  += dictionary[1]->GetSize();
1164     totalSizeDict2  += dictionary[2]->GetSize();
1165
1166     Float_t nPixel = nRowMax * nColMax * nTimeMax;
1167     printf("AliTRDdigitizer::MakeDigits -- ");
1168     printf("Found %d digits in detector %d (%3.0f).\n"
1169           ,nDigits,iDet
1170           ,100.0 * ((Float_t) nDigits) / nPixel);
1171  
1172     if (fCompress) signals->Compress(1,0);
1173
1174   }
1175
1176   printf("AliTRDdigitizer::MakeDigits -- ");
1177   printf("Total number of analyzed hits = %d\n",countHits);
1178
1179   printf("AliTRDdigitizer::MakeDigits -- ");
1180   printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
1181                                                     ,totalSizeDict0
1182                                                     ,totalSizeDict1
1183                                                     ,totalSizeDict2);        
1184
1185   return kTRUE;
1186
1187 }
1188
1189 //_____________________________________________________________________________
1190 Bool_t AliTRDdigitizer::CheckDetector(Int_t plane, Int_t chamber, Int_t sector)
1191 {
1192   //
1193   // Checks whether a detector is enabled
1194   //
1195
1196   if ((fTRD->GetSensChamber() >=       0) &&
1197       (fTRD->GetSensChamber() != chamber)) return kFALSE;
1198   if ((fTRD->GetSensPlane()   >=       0) &&
1199       (fTRD->GetSensPlane()   !=   plane)) return kFALSE;
1200   if ( fTRD->GetSensSector()  >=       0) {
1201     Int_t sens1 = fTRD->GetSensSector();
1202     Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
1203     sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
1204            * AliTRDgeometry::Nsect();
1205     if (sens1 < sens2) {
1206       if ((sector < sens1) || (sector >= sens2)) return kFALSE;
1207     }
1208     else {
1209       if ((sector < sens1) && (sector >= sens2)) return kFALSE;
1210     }
1211   }
1212
1213   return kTRUE;
1214
1215 }
1216
1217 //_____________________________________________________________________________
1218 Bool_t AliTRDdigitizer::WriteDigits()
1219 {
1220   //
1221   // Writes out the TRD-digits and the dictionaries
1222   //
1223
1224   // Create the branches
1225   if (!(gAlice->TreeD()->GetBranch("TRDdigits"))) { 
1226     return kFALSE;
1227   }
1228
1229   // Store the digits and the dictionary in the tree
1230   fDigits->WriteDigits();
1231
1232   // Write the new tree into the input file (use overwrite option)
1233   Char_t treeName[15];
1234   sprintf(treeName,"TreeD%d",fEvent);
1235   printf("AliTRDdigitizer::WriteDigits -- ");
1236   printf("Write the digits tree %s for event %d.\n"
1237         ,treeName,fEvent);
1238   gAlice->TreeD()->Write(treeName,TObject::kOverwrite);
1239  
1240   return kTRUE;
1241
1242 }
1243
1244 //_____________________________________________________________________________
1245 void AliTRDdigitizer::SetPRF(TF1 *prf)
1246 {
1247   //
1248   // Defines a new pad response function
1249   //
1250
1251   if (fPRF) delete fPRF;
1252   fPRF = prf;     
1253
1254 }
1255