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 }