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