Changed timebin 0 to be the one closest to the readout
[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.9  2000/10/02 21:28:19  fca
19 Removal of useless dependecies via forward declarations
20
21 Revision 1.8  2000/06/09 11:10:07  cblume
22 Compiler warnings and coding conventions, next round
23
24 Revision 1.7  2000/06/08 18:32:58  cblume
25 Make code compliant to coding conventions
26
27 Revision 1.6  2000/06/07 16:27:32  cblume
28 Try to remove compiler warnings on Sun and HP
29
30 Revision 1.5  2000/05/09 16:38:57  cblume
31 Removed PadResponse(). Merge problem
32
33 Revision 1.4  2000/05/08 15:53:45  cblume
34 Resolved merge conflict
35
36 Revision 1.3  2000/04/28 14:49:27  cblume
37 Only one declaration of iDict in MakeDigits()
38
39 Revision 1.1.4.1  2000/05/08 14:42:04  cblume
40 Introduced AliTRDdigitsManager
41
42 Revision 1.1  2000/02/28 19:00:13  cblume
43 Add new TRD classes
44
45 */
46
47 ///////////////////////////////////////////////////////////////////////////////
48 //                                                                           //
49 //  Creates and handles digits from TRD hits                                 //
50 //                                                                           //
51 //  The following effects are included:                                      //
52 //      - Diffusion                                                          //
53 //      - ExB effects                                                        //
54 //      - Gas gain including fluctuations                                    //
55 //      - Pad-response (simple Gaussian approximation)                       //
56 //      - Electronics noise                                                  //
57 //      - Electronics gain                                                   //
58 //      - Digitization                                                       //
59 //      - ADC threshold                                                      //
60 //  The corresponding parameter can be adjusted via the various              //
61 //  Set-functions. If these parameters are not explicitly set, default       //
62 //  values are used (see Init-function).                                     //
63 //  To produce digits from a root-file with TRD-hits use the                 //
64 //  slowDigitsCreate.C macro.                                                //
65 //                                                                           //
66 ///////////////////////////////////////////////////////////////////////////////
67
68 #include <stdlib.h>
69
70 #include <TMath.h>
71 #include <TVector.h>
72 #include <TRandom.h>
73 #include <TROOT.h>
74 #include <TTree.h>
75
76 #include "AliTRD.h"
77 #include "AliTRDdigitizer.h"
78 #include "AliTRDdataArrayI.h"
79 #include "AliTRDdataArrayF.h"
80 #include "AliTRDdigitsManager.h"
81
82 ClassImp(AliTRDdigitizer)
83
84 //_____________________________________________________________________________
85 AliTRDdigitizer::AliTRDdigitizer():TNamed()
86 {
87   //
88   // AliTRDdigitizer default constructor
89   //
90
91   fInputFile     = NULL;
92   fDigits        = NULL;
93   fTRD           = NULL;
94   fGeo           = NULL;
95   fPRF           = NULL;
96
97   fEvent         = 0;
98   fGasGain       = 0.0;
99   fNoise         = 0.0;
100   fChipGain      = 0.0;
101   fADCoutRange   = 0.0;
102   fADCinRange    = 0.0;
103   fADCthreshold  = 0;
104   fDiffusionOn   = 0;
105   fDiffusionT    = 0.0;
106   fDiffusionL    = 0.0;
107   fElAttachOn    = 0;
108   fElAttachProp  = 0.0;
109   fExBOn         = 0;
110   fLorentzAngle  = 0.0;
111
112 }
113
114 //_____________________________________________________________________________
115 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
116                 :TNamed(name,title)
117 {
118   //
119   // AliTRDdigitizer default constructor
120   //
121
122   fInputFile     = NULL;
123   fDigits        = NULL;
124   fTRD           = NULL;
125   fGeo           = NULL;
126
127   fEvent         = 0;
128
129   Init();
130
131 }
132
133 //_____________________________________________________________________________
134 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
135 {
136   //
137   // AliTRDdigitizer copy constructor
138   //
139
140   ((AliTRDdigitizer &) d).Copy(*this);
141
142 }
143
144 //_____________________________________________________________________________
145 AliTRDdigitizer::~AliTRDdigitizer()
146 {
147   //
148   // AliTRDdigitizer destructor
149   //
150
151   if (fInputFile) {
152     fInputFile->Close();
153     delete fInputFile;
154   }
155
156   if (fDigits) {
157     delete fDigits;
158   }
159
160   if (fPRF) delete fPRF;
161
162 }
163
164 //_____________________________________________________________________________
165 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
166 {
167   //
168   // Assignment operator
169   //
170
171   if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
172   return *this;
173
174 }
175
176 //_____________________________________________________________________________
177 void AliTRDdigitizer::Copy(TObject &d)
178 {
179   //
180   // Copy function
181   //
182
183   ((AliTRDdigitizer &) d).fInputFile     = NULL;
184   ((AliTRDdigitizer &) d).fDigits        = NULL;
185   ((AliTRDdigitizer &) d).fTRD           = NULL;
186   ((AliTRDdigitizer &) d).fGeo           = NULL;
187
188   ((AliTRDdigitizer &) d).fEvent         = 0;
189
190   ((AliTRDdigitizer &) d).fGasGain       = fGasGain;
191   ((AliTRDdigitizer &) d).fNoise         = fNoise;
192   ((AliTRDdigitizer &) d).fChipGain      = fChipGain;
193   ((AliTRDdigitizer &) d).fADCoutRange   = fADCoutRange;
194   ((AliTRDdigitizer &) d).fADCinRange    = fADCinRange;
195   ((AliTRDdigitizer &) d).fADCthreshold  = fADCthreshold;
196   ((AliTRDdigitizer &) d).fDiffusionOn   = fDiffusionOn; 
197   ((AliTRDdigitizer &) d).fDiffusionT    = fDiffusionT;
198   ((AliTRDdigitizer &) d).fDiffusionL    = fDiffusionL;
199   ((AliTRDdigitizer &) d).fElAttachOn    = fElAttachOn;
200   ((AliTRDdigitizer &) d).fElAttachProp  = fElAttachProp;
201   ((AliTRDdigitizer &) d).fExBOn         = fExBOn;
202   ((AliTRDdigitizer &) d).fLorentzAngle  = fLorentzAngle;
203   ((AliTRDdigitizer &) d).fLorentzFactor = fLorentzFactor;
204
205   fPRF->Copy(*((AliTRDdigitizer &) d).fPRF);
206
207 }
208
209 //_____________________________________________________________________________
210 Int_t AliTRDdigitizer::Diffusion(Float_t driftlength, Float_t *xyz)
211 {
212   //
213   // Applies the diffusion smearing to the position of a single electron
214   //
215
216   Float_t driftSqrt = TMath::Sqrt(driftlength);
217   Float_t sigmaT = driftSqrt * fDiffusionT;
218   Float_t sigmaL = driftSqrt * fDiffusionL;
219   xyz[0] = gRandom->Gaus(xyz[0], sigmaL * fLorentzFactor);
220   xyz[1] = gRandom->Gaus(xyz[1], sigmaT * fLorentzFactor);
221   xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
222   return 1;
223
224 }
225
226 //_____________________________________________________________________________
227 Int_t AliTRDdigitizer::ExB(Float_t driftlength, Float_t *xyz)
228 {
229   //
230   // Applies E x B effects to the position of a single electron
231   //
232
233   xyz[0] = xyz[0];
234   xyz[1] = xyz[1] + fLorentzAngle * driftlength;
235   xyz[2] = xyz[2];
236
237   return 1;
238
239 }
240
241 //_____________________________________________________________________________
242 void AliTRDdigitizer::Init()
243 {
244   //
245   // Initializes the digitization procedure with standard values
246   //
247
248   // The default parameter for the digitization
249   fGasGain       = 2.0E3;
250   fNoise         = 3000.;
251   fChipGain      = 10.;
252   fADCoutRange   = 255.;
253   fADCinRange    = 2000.;
254   fADCthreshold  = 1;
255
256   // Transverse and longitudinal diffusion coefficients (Xe/Isobutane)
257   fDiffusionOn   = 1;
258   fDiffusionT    = 0.060;
259   fDiffusionL    = 0.017;
260
261   // Propability for electron attachment
262   fElAttachOn    = 0;
263   fElAttachProp  = 0.0;
264
265   // E x B effects
266   fExBOn         = 0;
267   // omega * tau. (tau ~ 12 * 10^-12, B = 0.2T)
268   fLorentzAngle  = 17.6 * 12.0 * 0.2 * 0.01;
269
270   // The pad response function
271   fPRF           = new TF1("PRF","[0]*([1]+exp(-x*x/(2.0*[2])))",-2,2);
272   fPRF->SetParameter(0, 0.8872);
273   fPRF->SetParameter(1,-0.00573);
274   fPRF->SetParameter(2, 0.454 * 0.454);
275
276 }
277
278 //_____________________________________________________________________________
279 Bool_t AliTRDdigitizer::Open(const Char_t *name, Int_t nEvent)
280 {
281   //
282   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
283   //
284
285   // Connect the AliRoot file containing Geometry, Kine, and Hits
286   fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
287   if (!fInputFile) {
288     printf("AliTRDdigitizer::Open -- ");
289     printf("Open the ALIROOT-file %s.\n",name);
290     fInputFile = new TFile(name,"UPDATE");
291   }
292   else {
293     printf("AliTRDdigitizer::Open -- ");
294     printf("%s is already open.\n",name);
295   }
296
297   gAlice = (AliRun*) fInputFile->Get("gAlice");
298   if (gAlice) {
299     printf("AliTRDdigitizer::Open -- ");
300     printf("AliRun object found on file.\n");
301   }
302   else {
303     printf("AliTRDdigitizer::Open -- ");
304     printf("Could not find AliRun object.\n");
305     return kFALSE;
306   }
307
308   fEvent = nEvent;
309
310   // Import the Trees for the event nEvent in the file
311   Int_t nparticles = gAlice->GetEvent(fEvent);
312   if (nparticles <= 0) {
313     printf("AliTRDdigitizer::Open -- ");
314     printf("No entries in the trees for event %d.\n",fEvent);
315     return kFALSE;
316   }
317
318   // Get the pointer to the detector class and check for version 1
319   fTRD = (AliTRD*) gAlice->GetDetector("TRD");
320   if (fTRD->IsVersion() != 1) {
321     printf("AliTRDdigitizer::Open -- ");
322     printf("TRD must be version 1 (slow simulator).\n");
323     exit(1);
324   }
325
326   // Get the geometry
327   fGeo = fTRD->GetGeometry();
328   printf("AliTRDdigitizer::Open -- ");
329   printf("Geometry version %d\n",fGeo->IsVersion());
330
331   return kTRUE;
332
333 }
334
335 //_____________________________________________________________________________
336 Bool_t AliTRDdigitizer::MakeDigits()
337 {
338   //
339   // Loops through the TRD-hits and creates the digits.
340   //
341
342   ///////////////////////////////////////////////////////////////
343   // Parameter 
344   ///////////////////////////////////////////////////////////////
345
346   // Converts number of electrons to fC
347   const Float_t kEl2fC  = 1.602E-19 * 1.0E15; 
348
349   ///////////////////////////////////////////////////////////////
350
351   Int_t   iRow, iCol, iTime;
352   Int_t   nBytes = 0;
353   Int_t   iDict;
354
355   Int_t   totalSizeDigits = 0;
356   Int_t   totalSizeDict0  = 0;
357   Int_t   totalSizeDict1  = 0;
358   Int_t   totalSizeDict2  = 0;
359
360   AliTRDdataArrayI *digits;
361   AliTRDdataArrayI *dictionary[kNDict];
362
363   if (!fGeo) {
364     printf("AliTRDdigitizer::MakeDigits -- ");
365     printf("No geometry defined\n");
366     return kFALSE;
367   }
368
369   // Create a digits manager
370   fDigits = new AliTRDdigitsManager();
371
372   // Create detector arrays to keep the signal and track numbers
373   AliTRDdataArrayF *signal = new AliTRDdataArrayF();
374   AliTRDdataArrayI *tracks[kNDict];
375   for (iDict = 0; iDict < kNDict; iDict++) {
376     tracks[iDict] = new AliTRDdataArrayI();
377   }
378
379   // Get the pointer to the hit tree
380   TTree *hitTree = gAlice->TreeH();
381
382   // Get the number of entries in the hit tree
383   // (Number of primary particles creating a hit somewhere)
384   Int_t nTrack = (Int_t) hitTree->GetEntries();
385
386   printf("AliTRDdigitizer::MakeDigits -- ");
387   printf("Start creating digits.\n");
388
389   // The Lorentz factor
390   if (fExBOn) {
391     fLorentzFactor = 1.0 / (1.0 + fLorentzAngle*fLorentzAngle);
392   }
393   else {
394     fLorentzFactor = 1.0;
395   }
396
397   Int_t chamBeg = 0;
398   Int_t chamEnd = kNcham;
399   if (fTRD->GetSensChamber()  >= 0) {
400     chamBeg = fTRD->GetSensChamber();
401     chamEnd = chamBeg + 1;
402   }
403   Int_t planBeg = 0;
404   Int_t planEnd = kNplan;
405   if (fTRD->GetSensPlane()    >= 0) {
406     planBeg = fTRD->GetSensPlane();
407     planEnd = planBeg + 1;
408   }
409   Int_t sectBeg = 0;
410   Int_t sectEnd = kNsect;
411
412   Int_t countHits = 0;
413
414   // Loop through all the chambers
415   for (Int_t iCham = chamBeg; iCham < chamEnd; iCham++) {
416     for (Int_t iPlan = planBeg; iPlan < planEnd; iPlan++) {
417       for (Int_t iSect = sectBeg; iSect < sectEnd; iSect++) {
418
419         if (fTRD->GetSensSector() >= 0) {
420           Int_t sens1 = fTRD->GetSensSector();
421           Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
422           sens2 -= ((Int_t) (sens2 / kNsect)) * kNsect;
423           if (sens1 < sens2) {
424             if ((iSect < sens1) || (iSect >= sens2)) continue;
425           }
426           else {
427             if ((iSect < sens1) && (iSect >= sens2)) continue;
428           }
429         }
430
431         Int_t nDigits = 0;
432
433         printf("AliTRDdigitizer::MakeDigits -- ");
434         printf("Digitizing chamber %d, plane %d, sector %d.\n"
435               ,iCham,iPlan,iSect);
436
437         Int_t   iDet        = fGeo->GetDetector(iPlan,iCham,iSect);
438         Int_t   nRowMax     = fGeo->GetRowMax(iPlan,iCham,iSect);
439         Int_t   nColMax     = fGeo->GetColMax(iPlan);
440         Int_t   nTimeMax    = fGeo->GetTimeMax();
441         Float_t row0        = fGeo->GetRow0(iPlan,iCham,iSect);
442         Float_t col0        = fGeo->GetCol0(iPlan);
443         Float_t time0       = fGeo->GetTime0(iPlan);
444         Float_t rowPadSize  = fGeo->GetRowPadSize();
445         Float_t colPadSize  = fGeo->GetColPadSize();
446         Float_t timeBinSize = fGeo->GetTimeBinSize();
447
448         // Adjust the size of the detector arrays
449         signal->Allocate(nRowMax,nColMax,nTimeMax);
450         for (iDict = 0; iDict < kNDict; iDict++) {
451           tracks[iDict]->Allocate(nRowMax,nColMax,nTimeMax);
452         }
453
454         // Loop through all entries in the tree
455         for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
456
457           gAlice->ResetHits();
458           nBytes += hitTree->GetEvent(iTrack);
459
460           // Get the number of hits in the TRD created by this particle
461           Int_t nHit = fTRD->Hits()->GetEntriesFast();
462
463           // Loop through the TRD hits  
464           for (Int_t iHit = 0; iHit < nHit; iHit++) {
465
466             countHits++;
467
468             AliTRDhit *hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit);
469             Float_t pos[3];
470                     pos[0]   = hit->X();
471                     pos[1]   = hit->Y();
472                     pos[2]   = hit->Z();
473             Float_t q        = hit->GetCharge();
474             Int_t   track    = hit->Track();
475             Int_t   detector = hit->GetDetector();
476             Int_t   plane    = fGeo->GetPlane(detector);
477             Int_t   sector   = fGeo->GetSector(detector);
478             Int_t   chamber  = fGeo->GetChamber(detector);
479
480             if ((sector  != iSect) ||
481                 (plane   != iPlan) ||
482                 (chamber != iCham)) 
483               continue;
484
485             // Rotate the sectors on top of each other
486             Float_t rot[3];
487             fGeo->Rotate(detector,pos,rot);
488
489             // The hit position in pad coordinates (center pad)
490             // The pad row (z-direction)
491             Int_t  rowH = (Int_t) ((rot[2] -  row0) /  rowPadSize);
492             // The pad column (rphi-direction)  
493             Int_t  colH = (Int_t) ((rot[1] -  col0) /  colPadSize);
494             // The time bucket
495             Int_t timeH = (Int_t) ((rot[0] - time0) / timeBinSize);
496
497             // Array to sum up the signal in a box surrounding the
498             // hit postition
499             const Int_t kTimeBox = 7;
500             const Int_t  kColBox = 9;
501             const Int_t  kRowBox = 7;
502             Float_t signalSum[kRowBox][kColBox][kTimeBox];
503             for (iRow  = 0;  iRow <  kRowBox; iRow++ ) {
504               for (iCol  = 0;  iCol <  kColBox; iCol++ ) {
505                 for (iTime = 0; iTime < kTimeBox; iTime++) {
506                   signalSum[iRow][iCol][iTime] = 0;
507                 }
508               }
509             }
510
511             // Loop over all electrons of this hit
512             Int_t nEl = (Int_t) q;
513             for (Int_t iEl = 0; iEl < nEl; iEl++) {
514
515               // The driftlength
516               Float_t driftlength = rot[0] - time0;
517               if ((driftlength <        0) || 
518                   (driftlength > kDrThick)) break;
519               Float_t driftlengthL = driftlength;
520               if (fExBOn) driftlengthL /= TMath::Sqrt(fLorentzFactor);
521               Float_t xyz[3];
522               xyz[0] = rot[0];
523               xyz[1] = rot[1];
524               xyz[2] = rot[2];
525
526               // Electron attachment
527               if (fElAttachOn) {
528                 if (gRandom->Rndm() < (driftlengthL * fElAttachProp / 100.)) continue;
529               }
530
531               // Apply the diffusion smearing
532               if (fDiffusionOn) {
533                 if (!(Diffusion(driftlengthL,xyz))) continue;
534               }
535
536               // Apply E x B effects
537               if (fExBOn) { 
538                 if (!(ExB(driftlength,xyz))) continue;   
539               }
540
541               // The electron position and the distance to the hit position
542               // in pad units
543               // The pad row (z-direction)
544               Int_t  rowE = (Int_t) ((xyz[2] -  row0) /  rowPadSize);
545               Int_t  rowD =  rowH -  rowE;
546               // The pad column (rphi-direction)
547               Int_t  colE = (Int_t) ((xyz[1] -  col0) /  colPadSize);
548               Int_t  colD =  colH -  colE;
549               // The time bucket
550               Int_t timeE = (Int_t) ((xyz[0] - time0) / timeBinSize);
551               Int_t timeD = timeH - timeE;
552
553               // Apply the gas gain including fluctuations
554               Int_t signal = (Int_t) (-fGasGain * TMath::Log(gRandom->Rndm()));
555
556               // The distance of the electron to the center of the pad 
557               // in units of pad width
558               Float_t dist = (xyz[1] - col0 - (colE + 0.5) * colPadSize) 
559                            / colPadSize;
560
561               // Sum up the signal in the different pixels
562               // and apply the pad response
563               Int_t  rowIdx =  rowD + (Int_t) ( kRowBox / 2);
564               Int_t  colIdx =  colD + (Int_t) ( kColBox / 2);
565               Int_t timeIdx = timeD + (Int_t) (kTimeBox / 2);
566
567               if (( rowIdx < 0) || ( rowIdx >  kRowBox)) {
568                 printf("AliTRDdigitizer::MakeDigits -- ");
569                 printf("Boundary error. rowIdx = %d (%d)\n", rowIdx, kRowBox);
570                 continue;
571               }
572               if (( colIdx < 0) || ( colIdx >  kColBox)) {
573                 printf("AliTRDdigitizer::MakeDigits -- ");
574                 printf("Boundary error. colIdx = %d (%d)\n", colIdx, kColBox);
575                 continue;
576               }
577               if ((timeIdx < 0) || (timeIdx > kTimeBox)) {
578                 printf("AliTRDdigitizer::MakeDigits -- ");
579                 printf("Boundary error. timeIdx = %d (%d)\n",timeIdx,kTimeBox);
580                 continue;
581               }
582               signalSum[rowIdx][colIdx-1][timeIdx] += fPRF->Eval(dist-1.0,0,0) * signal;
583               signalSum[rowIdx][colIdx  ][timeIdx] += fPRF->Eval(dist    ,0,0) * signal;
584               signalSum[rowIdx][colIdx+1][timeIdx] += fPRF->Eval(dist+1.0,0,0) * signal;
585
586             }
587
588             // Add the padcluster to the detector matrix
589             for (iRow  = 0;  iRow <  kRowBox; iRow++ ) {
590               for (iCol  = 0;  iCol <  kColBox; iCol++ ) {
591                 for (iTime = 0; iTime < kTimeBox; iTime++) {
592
593                   Int_t  rowB =  rowH + iRow  - (Int_t) ( kRowBox / 2); 
594                   Int_t  colB =  colH + iCol  - (Int_t) ( kColBox / 2);
595                   Int_t timeB = timeH + iTime - (Int_t) (kTimeBox / 2);
596                   Float_t signalB = signalSum[iRow][iCol][iTime];
597                   if (( rowB < 0) || ( rowB >=  nRowMax)) continue;
598                   if (( colB < 0) || ( colB >=  nColMax)) continue;
599                   if ((timeB < 0) || (timeB >= nTimeMax)) continue;
600                   if (signalB > 0.0) {
601
602                     // Add the signal sum  
603                     signalB += signal->GetData(rowB,colB,timeB);
604                     signal->SetData(rowB,colB,timeB,signalB);  
605                     // Store the track index in the dictionary
606                     // Note: We store index+1 in order to allow the array to be compressed
607                     for (iDict = 0; iDict < kNDict; iDict++) {
608                       Int_t oldTrack = tracks[iDict]->GetData(rowB,colB,timeB);
609                       if (oldTrack == track+1) break;
610                       if (oldTrack ==      -1) break;
611                       if (oldTrack ==       0) {
612                         tracks[iDict]->SetData(rowB,colB,timeB,track+1);
613                         break;
614                       }
615                     }
616                     if (iDict == kNDict) {
617                       printf("AliTRDdigitizer::MakeDigits -- ");
618                       printf("More than three tracks for one digit!\n");
619                     }
620                   }
621
622                 }
623               }
624             }
625
626           }
627
628         }
629
630         // Add a container for the digits of this detector
631         digits = fDigits->GetDigits(iDet);        
632         // Allocate memory space for the digits buffer
633         digits->Allocate(nRowMax,nColMax,nTimeMax);
634
635         // Do the same for the dictionary arrays
636         for (iDict = 0; iDict < kNDict; iDict++) {
637           dictionary[iDict] = fDigits->GetDictionary(iDet,iDict);
638           dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeMax);
639         }
640
641         // Create the digits for this chamber
642         for (iRow  = 0; iRow  <  nRowMax; iRow++ ) {
643           for (iCol  = 0; iCol  <  nColMax; iCol++ ) {
644             for (iTime = 0; iTime < nTimeMax; iTime++) {         
645
646               Float_t signalAmp = signal->GetData(iRow,iCol,iTime);
647
648               // Add the noise
649               signalAmp  = TMath::Max((Float_t) gRandom->Gaus(signalAmp,fNoise)
650                                      ,(Float_t) 0.0);
651               // Convert to fC
652               signalAmp *= kEl2fC;
653               // Convert to mV
654               signalAmp *= fChipGain;
655               // Convert to ADC counts
656               Int_t adc  = (Int_t) (signalAmp * (fADCoutRange / fADCinRange));
657
658               if (adc > fADCthreshold) {
659
660                 nDigits++;
661
662                 // Store the amplitude of the digit
663                 digits->SetData(iRow,iCol,iTime,adc);
664
665                 // Store the track index in the dictionary
666                 // Note: We store index+1 in order to allow the array to be compressed
667                 for (iDict = 0; iDict < kNDict; iDict++) {
668                   dictionary[iDict]->SetData(iRow,iCol,iTime
669                                             ,tracks[iDict]->GetData(iRow,iCol,iTime));
670                 }
671
672               }
673
674             }
675           }
676         }
677
678         // Compress the arrays
679         digits->Compress(1,0);
680         for (iDict = 0; iDict < kNDict; iDict++) {
681           dictionary[iDict]->Compress(1,0);
682         }
683
684         totalSizeDigits += digits->GetSize();
685         totalSizeDict0  += dictionary[0]->GetSize();
686         totalSizeDict1  += dictionary[1]->GetSize();
687         totalSizeDict2  += dictionary[2]->GetSize();
688
689         printf("AliTRDdigitizer::MakeDigits -- ");
690         printf("Number of digits found: %d.\n",nDigits);
691  
692         // Reset the arrays
693         signal->Reset();
694         for (iDict = 0; iDict < kNDict; iDict++) {
695           tracks[iDict]->Reset();
696         }
697
698       }
699     }
700   }
701
702   printf("AliTRDdigitizer::MakeDigits -- ");
703   printf("Total number of analyzed hits = %d\n",countHits);
704
705   printf("AliTRDdigitizer::MakeDigits -- ");
706   printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
707                                                     ,totalSizeDict0
708                                                     ,totalSizeDict1
709                                                     ,totalSizeDict2);        
710
711   return kTRUE;
712
713 }
714
715 //_____________________________________________________________________________
716 Bool_t AliTRDdigitizer::WriteDigits()
717 {
718   //
719   // Writes out the TRD-digits and the dictionaries
720   //
721
722   // Create the branches
723   if (!(gAlice->TreeD()->GetBranch("TRDdigits"))) { 
724     if (!fDigits->MakeBranch()) return kFALSE;
725   }
726
727   // Store the digits and the dictionary in the tree
728   fDigits->WriteDigits();
729
730   // Write the new tree into the input file (use overwrite option)
731   Char_t treeName[7];
732   sprintf(treeName,"TreeD%d",fEvent);
733   printf("AliTRDdigitizer::WriteDigits -- ");
734   printf("Write the digits tree %s for event %d.\n"
735         ,treeName,fEvent);
736   gAlice->TreeD()->Write(treeName,2);
737  
738   return kTRUE;
739
740 }