Fixed bug in expansion routine of arrays (initialize buffers properly)
[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.14  2000/11/20 08:54:44  cblume
19 Switch off compression as default
20
21 Revision 1.13  2000/11/10 14:57:52  cblume
22 Changes in the geometry constants for the DEC compiler
23
24 Revision 1.12  2000/11/01 14:53:20  cblume
25 Merge with TRD-develop
26
27 Revision 1.1.4.9  2000/10/26 17:00:22  cblume
28 Fixed bug in CheckDetector()
29
30 Revision 1.1.4.8  2000/10/23 13:41:35  cblume
31 Added protection against Log(0) in the gas gain calulation
32
33 Revision 1.1.4.7  2000/10/17 02:27:34  cblume
34 Get rid of global constants
35
36 Revision 1.1.4.6  2000/10/16 01:16:53  cblume
37 Changed timebin 0 to be the one closest to the readout
38
39 Revision 1.1.4.5  2000/10/15 23:34:29  cblume
40 Faster version of the digitizer
41
42 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
43 Made Getters const
44
45 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
46 Replace include files by forward declarations
47
48 Revision 1.1.4.2  2000/09/22 14:41:10  cblume
49 Bug fix in PRF. Included time response. New structure
50
51 Revision 1.10  2000/10/05 07:27:53  cblume
52 Changes in the header-files by FCA
53
54 Revision 1.9  2000/10/02 21:28:19  fca
55 Removal of useless dependecies via forward declarations
56
57 Revision 1.8  2000/06/09 11:10:07  cblume
58 Compiler warnings and coding conventions, next round
59
60 Revision 1.7  2000/06/08 18:32:58  cblume
61 Make code compliant to coding conventions
62
63 Revision 1.6  2000/06/07 16:27:32  cblume
64 Try to remove compiler warnings on Sun and HP
65
66 Revision 1.5  2000/05/09 16:38:57  cblume
67 Removed PadResponse(). Merge problem
68
69 Revision 1.4  2000/05/08 15:53:45  cblume
70 Resolved merge conflict
71
72 Revision 1.3  2000/04/28 14:49:27  cblume
73 Only one declaration of iDict in MakeDigits()
74
75 Revision 1.1.4.1  2000/05/08 14:42:04  cblume
76 Introduced AliTRDdigitsManager
77
78 Revision 1.1  2000/02/28 19:00:13  cblume
79 Add new TRD classes
80
81 */
82
83 ///////////////////////////////////////////////////////////////////////////////
84 //                                                                           //
85 //  Creates and handles digits from TRD hits                                 //
86 //                                                                           //
87 //  The following effects are included:                                      //
88 //      - Diffusion                                                          //
89 //      - ExB effects                                                        //
90 //      - Gas gain including fluctuations                                    //
91 //      - Pad-response (simple Gaussian approximation)                       //
92 //      - Electronics noise                                                  //
93 //      - Electronics gain                                                   //
94 //      - Digitization                                                       //
95 //      - ADC threshold                                                      //
96 //  The corresponding parameter can be adjusted via the various              //
97 //  Set-functions. If these parameters are not explicitly set, default       //
98 //  values are used (see Init-function).                                     //
99 //  To produce digits from a root-file with TRD-hits use the                 //
100 //  slowDigitsCreate.C macro.                                                //
101 //                                                                           //
102 ///////////////////////////////////////////////////////////////////////////////
103
104 #include <stdlib.h>
105
106 #include <TMath.h>
107 #include <TVector.h>
108 #include <TRandom.h>
109 #include <TROOT.h>
110 #include <TTree.h>
111 #include <TFile.h>
112 #include <TF1.h>
113
114 #include "AliRun.h"
115
116 #include "AliTRD.h"
117 #include "AliTRDhit.h"
118 #include "AliTRDdigitizer.h"
119 #include "AliTRDdataArrayI.h"
120 #include "AliTRDdataArrayF.h"
121 #include "AliTRDsegmentArray.h"
122 #include "AliTRDdigitsManager.h"
123 #include "AliTRDgeometry.h"
124
125 ClassImp(AliTRDdigitizer)
126
127 //_____________________________________________________________________________
128 AliTRDdigitizer::AliTRDdigitizer():TNamed()
129 {
130   //
131   // AliTRDdigitizer default constructor
132   //
133
134   fInputFile     = NULL;
135   fDigits        = NULL;
136   fTRD           = NULL;
137   fGeo           = NULL;
138   fPRF           = NULL;
139   fTRF           = NULL;
140   fTRFint        = NULL;
141
142   fEvent         = 0;
143   fGasGain       = 0.0;
144   fNoise         = 0.0;
145   fChipGain      = 0.0;
146   fADCoutRange   = 0.0;
147   fADCinRange    = 0.0;
148   fADCthreshold  = 0;
149   fDiffusionOn   = 0;
150   fDiffusionT    = 0.0;
151   fDiffusionL    = 0.0;
152   fElAttachOn    = 0;
153   fElAttachProp  = 0.0;
154   fExBOn         = 0;
155   fOmegaTau      = 0.0;
156   fPRFOn         = 0;
157   fTRFOn         = 0;
158   fDriftVelocity = 0.0;
159
160   fTRFbin        = 0;
161   fTRFlo         = 0.0;
162   fTRFhi         = 0.0;
163   fTRFwid        = 0.0;
164   fCompress      = kTRUE;
165   fVerbose       = 1;
166
167 }
168
169 //_____________________________________________________________________________
170 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
171                 :TNamed(name,title)
172 {
173   //
174   // AliTRDdigitizer default constructor
175   //
176
177   fInputFile     = NULL;
178   fDigits        = NULL;
179   fTRD           = NULL;
180   fGeo           = NULL;
181
182   fEvent         = 0;
183
184   fCompress      = kFALSE;
185   fVerbose       = 1;
186
187   Init();
188
189 }
190
191 //_____________________________________________________________________________
192 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
193 {
194   //
195   // AliTRDdigitizer copy constructor
196   //
197
198   ((AliTRDdigitizer &) d).Copy(*this);
199
200 }
201
202 //_____________________________________________________________________________
203 AliTRDdigitizer::~AliTRDdigitizer()
204 {
205   //
206   // AliTRDdigitizer destructor
207   //
208
209   if (fInputFile) {
210     fInputFile->Close();
211     delete fInputFile;
212   }
213
214   if (fDigits) {
215     delete fDigits;
216   }
217
218   if (fPRF) delete fPRF;
219   if (fTRF) delete fTRF;
220
221 }
222
223 //_____________________________________________________________________________
224 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
225 {
226   //
227   // Assignment operator
228   //
229
230   if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
231   return *this;
232
233 }
234
235 //_____________________________________________________________________________
236 void AliTRDdigitizer::Copy(TObject &d)
237 {
238   //
239   // Copy function
240   //
241
242   ((AliTRDdigitizer &) d).fInputFile     = NULL;
243   ((AliTRDdigitizer &) d).fDigits        = NULL;
244   ((AliTRDdigitizer &) d).fTRD           = NULL;
245   ((AliTRDdigitizer &) d).fGeo           = NULL;
246
247   ((AliTRDdigitizer &) d).fEvent         = 0;
248
249   ((AliTRDdigitizer &) d).fGasGain       = fGasGain;
250   ((AliTRDdigitizer &) d).fNoise         = fNoise;
251   ((AliTRDdigitizer &) d).fChipGain      = fChipGain;
252   ((AliTRDdigitizer &) d).fADCoutRange   = fADCoutRange;
253   ((AliTRDdigitizer &) d).fADCinRange    = fADCinRange;
254   ((AliTRDdigitizer &) d).fADCthreshold  = fADCthreshold;
255   ((AliTRDdigitizer &) d).fDiffusionOn   = fDiffusionOn; 
256   ((AliTRDdigitizer &) d).fDiffusionT    = fDiffusionT;
257   ((AliTRDdigitizer &) d).fDiffusionL    = fDiffusionL;
258   ((AliTRDdigitizer &) d).fElAttachOn    = fElAttachOn;
259   ((AliTRDdigitizer &) d).fElAttachProp  = fElAttachProp;
260   ((AliTRDdigitizer &) d).fExBOn         = fExBOn;
261   ((AliTRDdigitizer &) d).fOmegaTau      = fOmegaTau;
262   ((AliTRDdigitizer &) d).fLorentzFactor = fLorentzFactor;
263   ((AliTRDdigitizer &) d).fPRFOn         = fPRFOn;
264   ((AliTRDdigitizer &) d).fTRFOn         = fTRFOn;
265
266   ((AliTRDdigitizer &) d).fCompress      = fCompress;
267   ((AliTRDdigitizer &) d).fVerbose       = fVerbose;
268
269   fPRF->Copy(*((AliTRDdigitizer &) d).fPRF);
270   fTRF->Copy(*((AliTRDdigitizer &) d).fTRF);
271
272   ((AliTRDdigitizer &) d).fTRFbin        = fTRFbin;
273   ((AliTRDdigitizer &) d).fTRFlo         = fTRFlo;
274   ((AliTRDdigitizer &) d).fTRFhi         = fTRFhi;
275   ((AliTRDdigitizer &) d).fTRFwid        = fTRFwid;
276   if (((AliTRDdigitizer &) d).fTRFint) delete ((AliTRDdigitizer &) d).fTRFint;
277   ((AliTRDdigitizer &) d).fTRFint = new Float_t[fTRFbin];
278   for (Int_t iBin = 0; iBin < fTRFbin; iBin++) {
279     ((AliTRDdigitizer &) d).fTRFint[iBin] = fTRFint[iBin];
280   }                                                                             
281
282 }
283
284 //_____________________________________________________________________________
285 Int_t AliTRDdigitizer::Diffusion(Float_t driftlength, Float_t *xyz)
286 {
287   //
288   // Applies the diffusion smearing to the position of a single electron
289   //
290
291   Float_t driftSqrt = TMath::Sqrt(driftlength);
292   Float_t sigmaT = driftSqrt * fDiffusionT;
293   Float_t sigmaL = driftSqrt * fDiffusionL;
294   xyz[0] = gRandom->Gaus(xyz[0], sigmaL * fLorentzFactor);
295   xyz[1] = gRandom->Gaus(xyz[1], sigmaT * fLorentzFactor);
296   xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
297
298   return 1;
299
300 }
301
302 //_____________________________________________________________________________
303 Int_t AliTRDdigitizer::ExB(Float_t driftlength, Float_t *xyz)
304 {
305   //
306   // Applies E x B effects to the position of a single electron
307   //
308
309   xyz[0] = xyz[0];
310   xyz[1] = xyz[1] + fOmegaTau * driftlength;
311   xyz[2] = xyz[2];
312
313   return 1;
314
315 }
316
317 //_____________________________________________________________________________
318 Int_t AliTRDdigitizer::PadResponse(Float_t signal, Float_t dist, Float_t *pad)
319 {
320   //
321   // Applies the pad response
322   //
323
324   if (fPRF) {
325     pad[0] = TMath::Max(fPRF->Eval(-1.0 - dist,0,0) * signal,0.0);
326     pad[1] = TMath::Max(fPRF->Eval(     - dist,0,0) * signal,0.0);
327     pad[2] = TMath::Max(fPRF->Eval( 1.0 - dist,0,0) * signal,0.0);
328     return 1;
329   }
330   else {
331     return 0;
332   }
333
334 }
335
336 //_____________________________________________________________________________
337 Float_t AliTRDdigitizer::TimeResponse(Float_t time)
338 {
339   //
340   // Applies the preamp shaper time response
341   //
342
343   Int_t iBin = ((Int_t) ((time - fTRFlo) / fTRFwid)); 
344   if ((iBin >= 0) && (iBin < fTRFbin)) {
345     return fTRFint[iBin];
346   }
347   else {
348     return 0.0;
349   }    
350
351 }
352
353 //_____________________________________________________________________________
354 void AliTRDdigitizer::Init()
355 {
356   //
357   // Initializes the digitization procedure with standard values
358   //
359
360   // The default parameter for the digitization
361   fGasGain       = 8.0E3;
362   fNoise         = 3000.;
363   fChipGain      = 10.;
364   fADCoutRange   = 255.;
365   fADCinRange    = 2000.;
366   fADCthreshold  = 1;
367
368   // Transverse and longitudinal diffusion coefficients (Xe/Isobutane)
369   fDiffusionOn   = 1;
370   fDiffusionT    = 0.060;
371   fDiffusionL    = 0.017;
372
373   // Propability for electron attachment
374   fElAttachOn    = 0;
375   fElAttachProp  = 0.0;
376
377   // E x B effects
378   fExBOn         = 0;
379   // omega * tau. (tau ~ 12 * 10^-12, B = 0.2T)
380   fOmegaTau      = 17.6 * 12.0 * 0.2 * 0.01;
381
382   // The pad response function
383   fPRFOn         = 1;
384   fPRF           = new TF1("PRF","[0]*([1]+exp(-x*x/(2.0*[2])))",-3,3);
385   fPRF->SetParameter(0, 0.8872);
386   fPRF->SetParameter(1,-0.00573);
387   fPRF->SetParameter(2, 0.454 * 0.454);
388
389   // The drift velocity (cm / mus)
390   fDriftVelocity = 1.0;
391
392   // The time response function
393   fTRFOn         = 1;
394   Float_t loTRF  = -200.0;
395   Float_t hiTRF  = 1000.0;
396   fTRF           = new TF1("TRF",TRFlandau,loTRF,hiTRF,3);
397   fTRF->SetParameter(0,  1.0 / 24.24249);
398   fTRF->SetParameter(1,  0.0);
399   fTRF->SetParameter(2, 25.0);
400   fTRFbin        = 1200;
401   fTRFlo         = loTRF * fDriftVelocity / 1000.0;
402   fTRFhi         = hiTRF * fDriftVelocity / 1000.0;
403   fTRFwid        = (fTRFhi - fTRFlo) / ((Float_t) fTRFbin);
404
405 }
406
407 //_____________________________________________________________________________
408 void AliTRDdigitizer::IntegrateTRF()
409 {
410   //
411   // Integrates the time response function over the time bin size
412   //
413
414   if (fTRFint) delete fTRFint;
415   fTRFint = new Float_t[fTRFbin];
416   Float_t hiTRF   = fTRFhi                 / fDriftVelocity * 1000.0;
417   Float_t loTRF   = fTRFlo                 / fDriftVelocity * 1000.0;
418   Float_t timeBin = fGeo->GetTimeBinSize() / fDriftVelocity * 1000.0;
419   Float_t binWidth = (hiTRF - loTRF) / ((Float_t) fTRFbin);
420   for (Int_t iBin = 0; iBin < fTRFbin; iBin++) {
421     Float_t bin = iBin * binWidth + loTRF - 0.5 * timeBin;
422     fTRFint[iBin] = fTRF->Integral(bin,bin + timeBin);
423   }
424
425 }
426
427 //_____________________________________________________________________________
428 Bool_t AliTRDdigitizer::Open(const Char_t *name, Int_t nEvent)
429 {
430   //
431   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
432   //
433
434   // Connect the AliRoot file containing Geometry, Kine, and Hits
435   fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
436   if (!fInputFile) {
437     printf("AliTRDdigitizer::Open -- ");
438     printf("Open the ALIROOT-file %s.\n",name);
439     fInputFile = new TFile(name,"UPDATE");
440   }
441   else {
442     printf("AliTRDdigitizer::Open -- ");
443     printf("%s is already open.\n",name);
444   }
445
446   gAlice = (AliRun*) fInputFile->Get("gAlice");
447   if (gAlice) {
448     printf("AliTRDdigitizer::Open -- ");
449     printf("AliRun object found on file.\n");
450   }
451   else {
452     printf("AliTRDdigitizer::Open -- ");
453     printf("Could not find AliRun object.\n");
454     return kFALSE;
455   }
456
457   fEvent = nEvent;
458
459   // Import the Trees for the event nEvent in the file
460   Int_t nparticles = gAlice->GetEvent(fEvent);
461   if (nparticles <= 0) {
462     printf("AliTRDdigitizer::Open -- ");
463     printf("No entries in the trees for event %d.\n",fEvent);
464     return kFALSE;
465   }
466
467   return InitDetector();
468
469 }
470
471 //_____________________________________________________________________________
472 Bool_t AliTRDdigitizer::InitDetector()
473 {
474   //
475   // Sets the pointer to the TRD detector and the geometry
476   //
477
478   // Get the pointer to the detector class and check for version 1
479   fTRD = (AliTRD*) gAlice->GetDetector("TRD");
480   if (fTRD->IsVersion() != 1) {
481     printf("AliTRDdigitizer::Open -- ");
482     printf("TRD must be version 1 (slow simulator).\n");
483     exit(1);
484   }
485
486   // Get the geometry
487   fGeo = fTRD->GetGeometry();
488   printf("AliTRDdigitizer::Open -- ");
489   printf("Geometry version %d\n",fGeo->IsVersion());
490
491   return kTRUE;
492
493 }
494
495 //_____________________________________________________________________________
496 Bool_t AliTRDdigitizer::MakeDigits()
497 {
498   //
499   // Loops through the TRD-hits and creates the digits.
500   //
501
502   ///////////////////////////////////////////////////////////////
503   // Parameter 
504   ///////////////////////////////////////////////////////////////
505
506   // Converts number of electrons to fC
507   const Float_t kEl2fC  = 1.602E-19 * 1.0E15; 
508
509   ///////////////////////////////////////////////////////////////
510
511   // Number of pads included in the pad response
512   const Int_t kNpad  = 3;
513
514   // Number of track dictionary arrays
515   const Int_t kNDict = AliTRDdigitsManager::kNDict;
516
517   Int_t   iRow, iCol, iTime;
518   Int_t   iDict;
519   Int_t   nBytes = 0;
520
521   Int_t   totalSizeDigits = 0;
522   Int_t   totalSizeDict0  = 0;
523   Int_t   totalSizeDict1  = 0;
524   Int_t   totalSizeDict2  = 0;
525
526   AliTRDdataArrayF *signals = 0;
527   AliTRDdataArrayI *digits  = 0;
528   AliTRDdataArrayI *dictionary[kNDict];
529
530   // Create a digits manager
531   fDigits = new AliTRDdigitsManager();
532
533   // Create a container for the amplitudes
534   AliTRDsegmentArray *signalsArray 
535                      = new AliTRDsegmentArray("AliTRDdataArrayF",AliTRDgeometry::Ndet());
536
537   if (!fGeo) {
538     printf("AliTRDdigitizer::MakeDigits -- ");
539     printf("No geometry defined\n");
540     return kFALSE;
541   }
542
543   printf("AliTRDdigitizer::MakeDigits -- ");
544   printf("Start creating digits.\n");
545   if (fVerbose > 0) this->Dump();
546
547   // The Lorentz factor
548   if (fExBOn) {
549     fLorentzFactor = 1.0 / (1.0 + fOmegaTau*fOmegaTau);
550   }
551   else {
552     fLorentzFactor = 1.0;
553   }
554
555   // Create the integrated TRF
556   IntegrateTRF();
557
558   // Get the pointer to the hit tree
559   TTree *HitTree = gAlice->TreeH();
560
561   // Get the number of entries in the hit tree
562   // (Number of primary particles creating a hit somewhere)
563   Int_t nTrack = (Int_t) HitTree->GetEntries();
564   if (fVerbose > 0) {
565     printf("AliTRDdigitizer::MakeDigits -- ");
566     printf("Found %d primary particles\n",nTrack);
567   } 
568
569   Int_t detectorOld = -1;
570   Int_t countHits   =  0; 
571
572   // Loop through all entries in the tree
573   for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
574
575     gAlice->ResetHits();
576     nBytes += HitTree->GetEvent(iTrack);
577
578     // Get the number of hits in the TRD created by this particle
579     Int_t nHit = fTRD->Hits()->GetEntriesFast();
580     if (fVerbose > 0) {
581       printf("AliTRDdigitizer::MakeDigits -- ");
582       printf("Found %d hits for primary particle %d\n",nHit,iTrack);
583     }
584
585     // Loop through the TRD hits  
586     for (Int_t iHit = 0; iHit < nHit; iHit++) {
587
588       countHits++;
589
590       AliTRDhit *hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit);
591       Float_t pos[3];
592               pos[0]   = hit->X();
593               pos[1]   = hit->Y();
594               pos[2]   = hit->Z();
595       Float_t q        = hit->GetCharge();
596       Int_t   track    = hit->Track();
597       Int_t   detector = hit->GetDetector();
598       Int_t   plane    = fGeo->GetPlane(detector);
599       Int_t   sector   = fGeo->GetSector(detector);
600       Int_t   chamber  = fGeo->GetChamber(detector);
601
602       if (!(CheckDetector(plane,chamber,sector))) continue;
603
604       Int_t   nRowMax     = fGeo->GetRowMax(plane,chamber,sector);
605       Int_t   nColMax     = fGeo->GetColMax(plane);
606       Int_t   nTimeMax    = fGeo->GetTimeMax();
607       Float_t row0        = fGeo->GetRow0(plane,chamber,sector);
608       Float_t col0        = fGeo->GetCol0(plane);
609       Float_t time0       = fGeo->GetTime0(plane);
610       Float_t rowPadSize  = fGeo->GetRowPadSize();
611       Float_t colPadSize  = fGeo->GetColPadSize();
612       Float_t timeBinSize = fGeo->GetTimeBinSize();
613
614       if (fVerbose > 1) {
615         printf("Analyze hit no. %d ",iHit);
616         printf("-----------------------------------------------------------\n");
617         hit->Dump();
618         printf("plane = %d, sector = %d, chamber = %d\n"
619               ,plane,sector,chamber);
620         printf("nRowMax = %d, nColMax = %d, nTimeMax = %d\n" 
621               ,nRowMax,nColMax,nTimeMax);
622         printf("row0 = %f, col0 = %f, time0 = %f\n"
623               ,row0,col0,time0);
624         printf("rowPadSize = %f, colPadSize = %f, timeBinSize = %f\n"
625                ,rowPadSize,colPadSize,timeBinSize); 
626       }
627        
628       // Don't analyze test hits with amplitude 0.
629       if (((Int_t) q) == 0) continue;
630
631       if (detector != detectorOld) {
632         if (fVerbose > 1) {
633           printf("AliTRDdigitizer::MakeDigits -- ");
634           printf("Get new container. New det = %d, Old det = %d\n"
635                 ,detector,detectorOld);
636         }
637         // Compress the old one if enabled
638         if ((fCompress) && (detectorOld > -1)) {
639           if (fVerbose > 1) {
640             printf("AliTRDdigitizer::MakeDigits -- ");
641             printf("Compress the old container ... ");
642           }
643           signals->Compress(1,0);
644           for (iDict = 0; iDict < kNDict; iDict++) {
645             dictionary[iDict]->Compress(1,0);
646           }
647           if (fVerbose > 1) printf("done\n");
648         }
649         // Get the new container
650         signals = (AliTRDdataArrayF *) signalsArray->At(detector);
651         if (signals->GetNtime() == 0) {
652           // Allocate a new one if not yet existing
653           if (fVerbose > 1) {
654             printf("AliTRDdigitizer::MakeDigits -- ");
655             printf("Allocate a new container ... ");
656           }
657           signals->Allocate(nRowMax,nColMax,nTimeMax);
658         }
659         else {
660           // Expand an existing one
661           if (fCompress) {
662             if (fVerbose > 1) {
663               printf("AliTRDdigitizer::MakeDigits -- ");
664               printf("Expand an existing container ... ");
665             }
666             signals->Expand();
667           }
668         }
669         // The same for the dictionary
670         for (iDict = 0; iDict < kNDict; iDict++) {       
671           dictionary[iDict] = fDigits->GetDictionary(detector,iDict);
672           if (dictionary[iDict]->GetNtime() == 0) {
673             dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeMax);
674           }
675           else {
676             if (fCompress) dictionary[iDict]->Expand();
677           }
678         }      
679         if (fVerbose > 1) printf("done\n");
680         detectorOld = detector;
681       }
682
683       // Rotate the sectors on top of each other
684       Float_t rot[3];
685       fGeo->Rotate(detector,pos,rot);
686
687       // The driftlength
688       Float_t driftlength = time0 - rot[0];
689       if ((driftlength < 0) || 
690           (driftlength > AliTRDgeometry::DrThick())) continue;
691       Float_t driftlengthL = driftlength;
692       if (fExBOn) driftlengthL /= TMath::Sqrt(fLorentzFactor);
693
694       // The hit position in pad coordinates (center pad)
695       // The pad row (z-direction)
696       Int_t  rowH = (Int_t) ((rot[2] -  row0) /  rowPadSize);
697       // The pad column (rphi-direction)  
698       Int_t  colH = (Int_t) ((rot[1] -  col0) /  colPadSize);
699       // The time bucket
700       Int_t timeH = (Int_t) (driftlength      / timeBinSize);
701       if (fVerbose > 1) {
702         printf("rowH = %d, colH = %d, timeH = %d\n"
703               ,rowH,colH,timeH);
704       }
705
706       // Loop over all electrons of this hit
707       // TR photons produce hits with negative charge
708       Int_t nEl = ((Int_t) TMath::Abs(q));
709       for (Int_t iEl = 0; iEl < nEl; iEl++) {
710
711         Float_t xyz[3];
712         xyz[0] = rot[0];
713         xyz[1] = rot[1];
714         xyz[2] = rot[2];
715
716         // Electron attachment
717         if (fElAttachOn) {
718           if (gRandom->Rndm() < (driftlengthL * fElAttachProp / 100.)) 
719             continue;
720         }
721
722         // Apply the diffusion smearing
723         if (fDiffusionOn) {
724           if (!(Diffusion(driftlengthL,xyz))) continue;
725         }
726
727         // Apply E x B effects
728         if (fExBOn) { 
729           if (!(ExB(driftlength,xyz))) continue;   
730         }
731
732         // The electron position 
733         // The pad row (z-direction)
734         Int_t  rowE = (Int_t) ((xyz[2] -  row0) /  rowPadSize);
735         // The pad column (rphi-direction)
736         Int_t  colE = (Int_t) ((xyz[1] -  col0) /  colPadSize);
737         // The time bucket
738         Int_t timeE = (Int_t) ((time0 - xyz[0]) / timeBinSize);
739
740         if (( rowE < 0) || ( rowE >=  nRowMax)) continue;
741         if (( colE < 0) || ( colE >=  nColMax)) continue;
742         if ((timeE < 0) || (timeE >= nTimeMax)) continue;
743
744         // Apply the gas gain including fluctuations
745         Float_t ggRndm = 0.0;
746         do {
747           ggRndm = gRandom->Rndm();
748         } while (ggRndm <= 0);
749         Int_t signal = (Int_t) (-fGasGain * TMath::Log(ggRndm));
750
751         if (fVerbose > 2) {
752           printf("  electron no. %d, signal = %d\n",iEl,signal);
753           printf("  rowE = %d, colE = %d, timeE = %d\n"
754                 ,rowE,colE,timeE);
755         }
756
757         // Apply the pad response 
758         Float_t padSignal[kNpad];
759         if (fPRFOn) {
760           // The distance of the electron to the center of the pad 
761           // in units of pad width
762           Float_t dist = (xyz[1] - col0 - (colE + 0.5) * colPadSize) 
763                        / colPadSize;
764           if (!(PadResponse(signal,dist,padSignal))) continue;
765         }
766         else {
767           padSignal[0] = 0.0;
768           padSignal[1] = signal;
769           padSignal[2] = 0.0;
770         }
771
772         // The distance of the position to the beginning of the timebin
773         Float_t timeOffset = (time0 - timeE * timeBinSize) - xyz[0];
774         Int_t timeTRDbeg = 0;
775         Int_t timeTRDend = 1;
776         if (fTRFOn) {
777           timeTRDbeg =  2;
778           timeTRDend = 11;
779         }
780         for (Int_t iTimeBin = TMath::Max(timeE - timeTRDbeg,       0)
781                  ; iTimeBin < TMath::Min(timeE + timeTRDend,nTimeMax) 
782                  ; iTimeBin++) {
783
784           // Apply the time response
785           Float_t timeResponse = 1.0;
786           if (fTRFOn) {
787             Float_t time = (iTimeBin - timeE) * timeBinSize + timeOffset;
788             timeResponse = TimeResponse(time);
789           }
790
791           // Add the signals
792           Float_t signalOld[kNpad] = { 0.0, 0.0, 0.0 };
793           for (Int_t iPad = 0; iPad < kNpad; iPad++) {
794             Int_t colPos = colE + iPad - 1;
795             if (colPos <        0) continue;
796             if (colPos >= nColMax) break;
797             signalOld[iPad]  = signals->GetData(rowE,colPos,iTimeBin);
798             signalOld[iPad] += padSignal[iPad] * timeResponse;
799             signals->SetData(rowE,colPos,iTimeBin,signalOld[iPad]);
800           }
801           if (fVerbose > 3) {
802             printf("    iTimeBin = %d, timeResponse = %f\n"
803                   ,iTimeBin,timeResponse);
804             printf("    pad-signal = %f, %f, %f\n"
805                    ,signalOld[0],signalOld[1],signalOld[2]);
806           }
807
808           // Store the track index in the dictionary
809           // Note: We store index+1 in order to allow the array to be compressed
810           for (iDict = 0; iDict < kNDict; iDict++) {
811             Int_t oldTrack = dictionary[iDict]->GetData(rowE,colE,timeE);
812             if (oldTrack == track+1) break;
813             //if (oldTrack ==      -1) break;
814             if (oldTrack ==       0) {
815               dictionary[iDict]->SetData(rowE,colE,timeE,track+1);
816               if (fVerbose > 3) {
817                 printf("    track index = %d\n",track); 
818               }
819               break;
820             }
821           }
822           if ((fVerbose > 1) && (iDict == kNDict)) {
823             printf("AliTRDdigitizer::MakeDigits -- ");
824             printf("More than three tracks for one digit!\n");
825           }
826
827         }
828
829       }
830
831     }
832
833   } // All hits finished
834
835   printf("AliTRDdigitizer::MakeDigits -- ");
836   printf("Finished analyzing %d hits\n",countHits);
837
838   // Loop through all chambers to finalize the digits
839   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
840
841     Int_t plane    = fGeo->GetPlane(iDet);
842     Int_t sector   = fGeo->GetSector(iDet);
843     Int_t chamber  = fGeo->GetChamber(iDet);
844     Int_t nRowMax  = fGeo->GetRowMax(plane,chamber,sector);
845     Int_t nColMax  = fGeo->GetColMax(plane);
846     Int_t nTimeMax = fGeo->GetTimeMax();
847
848     //if (!(CheckDetector(plane,chamber,sector))) continue;
849     if (fVerbose > 0) {
850       printf("AliTRDdigitizer::MakeDigits -- ");
851       printf("Digitization for chamber %d\n",iDet);
852     }
853
854     // Add a container for the digits of this detector
855     digits = fDigits->GetDigits(iDet);        
856     // Allocate memory space for the digits buffer
857     digits->Allocate(nRowMax,nColMax,nTimeMax);
858
859     // Get the signal container
860     signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
861     if (signals->GetNtime() == 0) {
862       // Create missing containers
863       signals->Allocate(nRowMax,nColMax,nTimeMax);      
864     }
865     else {
866       // Expand the container if neccessary
867       if (fCompress) signals->Expand();
868     }
869     // Create the missing dictionary containers
870     for (iDict = 0; iDict < kNDict; iDict++) {       
871       dictionary[iDict] = fDigits->GetDictionary(iDet,iDict);
872       if (dictionary[iDict]->GetNtime() == 0) {
873         dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeMax);
874       }
875     }
876
877     Int_t nDigits = 0;
878
879     // Create the digits for this chamber
880     for (iRow  = 0; iRow  <  nRowMax; iRow++ ) {
881       for (iCol  = 0; iCol  <  nColMax; iCol++ ) {
882         for (iTime = 0; iTime < nTimeMax; iTime++) {         
883
884           Float_t signalAmp = signals->GetData(iRow,iCol,iTime);
885
886           // Add the noise
887           signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,fNoise),0.0);
888           // Convert to fC
889           signalAmp *= kEl2fC;
890           // Convert to mV
891           signalAmp *= fChipGain;
892           // Convert to ADC counts. Set the overflow-bit fADCoutRange if the 
893           // signal is larger than fADCinRange
894           Int_t adc  = 0;
895           if (signalAmp >= fADCinRange) {
896             adc = ((Int_t) fADCoutRange);
897           }
898           else {
899             adc = ((Int_t) (signalAmp * (fADCoutRange / fADCinRange)));
900           }
901
902           // Store the amplitude of the digit if above threshold
903           if (adc > fADCthreshold) {
904             if (fVerbose > 2) {
905               printf("  iRow = %d, iCol = %d, iTime = %d\n"
906                     ,iRow,iCol,iTime);
907               printf("  signal = %f, adc = %d\n",signalAmp,adc);
908             }
909             nDigits++;
910             digits->SetData(iRow,iCol,iTime,adc);
911           }
912
913         }
914       }
915     }
916
917     // Compress the arrays
918     digits->Compress(1,0);
919     for (iDict = 0; iDict < kNDict; iDict++) {
920       dictionary[iDict]->Compress(1,0);
921     }
922
923     totalSizeDigits += digits->GetSize();
924     totalSizeDict0  += dictionary[0]->GetSize();
925     totalSizeDict1  += dictionary[1]->GetSize();
926     totalSizeDict2  += dictionary[2]->GetSize();
927
928     Float_t nPixel = nRowMax * nColMax * nTimeMax;
929     printf("AliTRDdigitizer::MakeDigits -- ");
930     printf("Found %d digits in detector %d (%3.0f).\n"
931           ,nDigits,iDet
932           ,100.0 * ((Float_t) nDigits) / nPixel);
933  
934     if (fCompress) signals->Compress(1,0);
935
936   }
937
938   printf("AliTRDdigitizer::MakeDigits -- ");
939   printf("Total number of analyzed hits = %d\n",countHits);
940
941   printf("AliTRDdigitizer::MakeDigits -- ");
942   printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
943                                                     ,totalSizeDict0
944                                                     ,totalSizeDict1
945                                                     ,totalSizeDict2);        
946
947   return kTRUE;
948
949 }
950
951 //_____________________________________________________________________________
952 Bool_t AliTRDdigitizer::CheckDetector(Int_t plane, Int_t chamber, Int_t sector)
953 {
954   //
955   // Checks whether a detector is enabled
956   //
957
958   if ((fTRD->GetSensChamber() >=       0) &&
959       (fTRD->GetSensChamber() != chamber)) return kFALSE;
960   if ((fTRD->GetSensPlane()   >=       0) &&
961       (fTRD->GetSensPlane()   !=   plane)) return kFALSE;
962   if ( fTRD->GetSensSector()  >=       0) {
963     Int_t sens1 = fTRD->GetSensSector();
964     Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
965     sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
966            * AliTRDgeometry::Nsect();
967     if (sens1 < sens2) {
968       if ((sector < sens1) || (sector >= sens2)) return kFALSE;
969     }
970     else {
971       if ((sector < sens1) && (sector >= sens2)) return kFALSE;
972     }
973   }
974
975   return kTRUE;
976
977 }
978
979 //_____________________________________________________________________________
980 Bool_t AliTRDdigitizer::WriteDigits()
981 {
982   //
983   // Writes out the TRD-digits and the dictionaries
984   //
985
986   // Create the branches
987   if (!(gAlice->TreeD()->GetBranch("TRDdigits"))) { 
988     if (!fDigits->MakeBranch()) return kFALSE;
989   }
990
991   // Store the digits and the dictionary in the tree
992   fDigits->WriteDigits();
993
994   // Write the new tree into the input file (use overwrite option)
995   Char_t treeName[7];
996   sprintf(treeName,"TreeD%d",fEvent);
997   printf("AliTRDdigitizer::WriteDigits -- ");
998   printf("Write the digits tree %s for event %d.\n"
999         ,treeName,fEvent);
1000   gAlice->TreeD()->Write(treeName,2);
1001  
1002   return kTRUE;
1003
1004 }
1005
1006 //_____________________________________________________________________________
1007 void AliTRDdigitizer::SetPRF(TF1 *prf)
1008 {
1009   //
1010   // Defines a new pad response function
1011   //
1012
1013   if (fPRF) delete fPRF;
1014   fPRF = prf;     
1015
1016 }
1017
1018 //_____________________________________________________________________________
1019 void AliTRDdigitizer::SetTRF(TF1 *trf)
1020 {
1021   //
1022   // Defines a new time response function
1023   //
1024
1025   if (fTRF) delete fTRF;
1026   fTRF = trf;     
1027
1028 }
1029
1030 //_____________________________________________________________________________
1031 Double_t TRFlandau(Double_t *x, Double_t *par)
1032 {
1033
1034   Double_t xx = x[0];
1035   Double_t landau = par[0] * TMath::Landau(xx,par[1],par[2]);
1036
1037   return landau;
1038
1039 }