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