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.32  2002/02/12 16:07:21  cblume
  19 Add new constructor
  20
  21 Revision 1.31  2002/02/11 14:27:11  cblume
  22 New pad plane design, new TRF+PRF, tail cancelation, cross talk
  23
  24 Revision 1.30  2001/11/19 08:44:08  cblume
  25 Fix bugs reported by Rene
  26
  27 Revision 1.29  2001/11/14 19:44:25  hristov
  28 Numeric const casted (Alpha)
  29
  30 Revision 1.28  2001/11/14 16:35:58  cblume
  31 Inherits now from AliDetector
  32
  33 Revision 1.27  2001/11/14 10:50:45  cblume
  34 Changes in digits IO. Add merging of summable digits
  35
  36 Revision 1.26  2001/11/06 17:19:41  cblume
  37 Add detailed geometry and simple simulator
  38
  39 Revision 1.25  2001/06/27 09:54:44  cblume
  40 Moved fField initialization to InitDetector()
  41
  42 Revision 1.24  2001/05/21 16:45:47  hristov
  43 Last minute changes (C.Blume)
  44
  45 Revision 1.23  2001/05/07 08:04:48  cblume
  46 New TRF and PRF. Speedup of the code. Digits from amplification region included
  47
  48 Revision 1.22  2001/03/30 14:40:14  cblume
  49 Update of the digitization parameter
  50
  51 Revision 1.21  2001/03/13 09:30:35  cblume
  52 Update of digitization. Moved digit branch definition to AliTRD
  53
  54 Revision 1.20  2001/02/25 20:19:00  hristov
  55 Minor correction: loop variable declared only once for HP, Sun
  56
  57 Revision 1.19  2001/02/14 18:22:26  cblume
  58 Change in the geometry of the padplane
  59
  60 Revision 1.18  2001/01/26 19:56:57  hristov
  61 Major upgrade of AliRoot code
  62
  63 Revision 1.17  2000/12/08 12:53:27  cblume
  64 Change in Copy() function for HP-compiler
  65
  66 Revision 1.16  2000/12/07 12:20:46  cblume
  67 Go back to array compression. Use sampled PRF to speed up digitization
  68
  69 Revision 1.15  2000/11/23 14:34:08  cblume
  70 Fixed bug in expansion routine of arrays (initialize buffers properly)
  71
  72 Revision 1.14  2000/11/20 08:54:44  cblume
  73 Switch off compression as default
  74
  75 Revision 1.13  2000/11/10 14:57:52  cblume
  76 Changes in the geometry constants for the DEC compiler
  77
  78 Revision 1.12  2000/11/01 14:53:20  cblume
  79 Merge with TRD-develop
  80
  81 Revision 1.1.4.9  2000/10/26 17:00:22  cblume
  82 Fixed bug in CheckDetector()
  83
  84 Revision 1.1.4.8  2000/10/23 13:41:35  cblume
  85 Added protection against Log(0) in the gas gain calulation
  86
  87 Revision 1.1.4.7  2000/10/17 02:27:34  cblume
  88 Get rid of global constants
  89
  90 Revision 1.1.4.6  2000/10/16 01:16:53  cblume
  91 Changed timebin 0 to be the one closest to the readout
  92
  93 Revision 1.1.4.5  2000/10/15 23:34:29  cblume
  94 Faster version of the digitizer
  95
  96 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
  97 Made Getters const
  98
  99 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
 100 Replace include files by forward declarations
 101
 102 Revision 1.1.4.2  2000/09/22 14:41:10  cblume
 103 Bug fix in PRF. Included time response. New structure
 104
 105 Revision 1.10  2000/10/05 07:27:53  cblume
 106 Changes in the header-files by FCA
 107
 108 Revision 1.9  2000/10/02 21:28:19  fca
 109 Removal of useless dependecies via forward declarations
 110
 111 Revision 1.8  2000/06/09 11:10:07  cblume
 112 Compiler warnings and coding conventions, next round
 113
 114 Revision 1.7  2000/06/08 18:32:58  cblume
 115 Make code compliant to coding conventions
 116
 117 Revision 1.6  2000/06/07 16:27:32  cblume
 118 Try to remove compiler warnings on Sun and HP
 119
 120 Revision 1.5  2000/05/09 16:38:57  cblume
 121 Removed PadResponse(). Merge problem
 122
 123 Revision 1.4  2000/05/08 15:53:45  cblume
 124 Resolved merge conflict
 125
 126 Revision 1.3  2000/04/28 14:49:27  cblume
 127 Only one declaration of iDict in MakeDigits()
 128
 129 Revision 1.1.4.1  2000/05/08 14:42:04  cblume
 130 Introduced AliTRDdigitsManager
 131
 132 Revision 1.1  2000/02/28 19:00:13  cblume
 133 Add new TRD classes
 134
 135 */
 136
 137 ///////////////////////////////////////////////////////////////////////////////
 138 //                                                                           //
 139 //  Creates and handles digits from TRD hits                                 //
 140 //  Author: C. Blume (C.Blume@gsi.de)                                        //
 141 //                                                                           //
 142 //  The following effects are included:                                      //
 143 //      - Diffusion                                                          //
 144 //      - ExB effects                                                        //
 145 //      - Gas gain including fluctuations                                    //
 146 //      - Pad-response (simple Gaussian approximation)                       //
 147 //      - Time-response                                                      //
 148 //      - Electronics noise                                                  //
 149 //      - Electronics gain                                                   //
 150 //      - Digitization                                                       //
 151 //      - ADC threshold                                                      //
 152 //  The corresponding parameter can be adjusted via the various              //
 153 //  Set-functions. If these parameters are not explicitly set, default       //
 154 //  values are used (see Init-function).                                     //
 155 //  As an example on how to use this class to produce digits from hits       //
 156 //  have a look at the macro hits2digits.C                                   //
 157 //  The production of summable digits is demonstrated in hits2sdigits.C      //
 158 //  and the subsequent conversion of the s-digits into normal digits is      //
 159 //  explained in sdigits2digits.C.                                           //
 160 //                                                                           //
 161 ///////////////////////////////////////////////////////////////////////////////
 162
 163 #include <stdlib.h>
 164
 165 #include <TMath.h>
 166 #include <TVector.h>
 167 #include <TRandom.h>
 168 #include <TROOT.h>
 169 #include <TTree.h>
 170 #include <TFile.h>
 171 #include <TF1.h>
 172 #include <TList.h>
 173 #include <TTask.h>
 174
 175 #include "AliRun.h"
 176 #include "AliMagF.h"
 177 #include "AliRunDigitizer.h"
 178
 179 #include "AliTRD.h"
 180 #include "AliTRDhit.h"
 181 #include "AliTRDdigitizer.h"
 182 #include "AliTRDdataArrayI.h"
 183 #include "AliTRDdataArrayF.h"
 184 #include "AliTRDsegmentArray.h"
 185 #include "AliTRDdigitsManager.h"
 186 #include "AliTRDgeometry.h"
 187
 188 ClassImp(AliTRDdigitizer)
 189
 190 //_____________________________________________________________________________
 191 AliTRDdigitizer::AliTRDdigitizer()
 192 {
 193   //
 194   // AliTRDdigitizer default constructor
 195   //
 196
 197   fInputFile          = NULL;
 198   fDigitsManager      = NULL;
 199   fSDigitsManagerList = NULL;
 200   fSDigitsManager     = NULL;
 201   fTRD                = NULL;
 202   fGeo                = NULL;
 203   fPRFsmp             = NULL;
 204   fTRFsmp             = NULL;
 205   fCTsmp              = NULL;
 206
 207   fMasks              = 0;
 208
 209   fEvent              = 0;
 210   fGasGain            = 0.0;
 211   fNoise              = 0.0;
 212   fChipGain           = 0.0;
 213   fADCoutRange        = 0.0;
 214   fADCinRange         = 0.0;
 215   fADCthreshold       = 0;
 216   fDiffusionOn        = 0;
 217   fDiffusionT         = 0.0;
 218   fDiffusionL         = 0.0;
 219   fElAttachOn         = 0;
 220   fElAttachProp       = 0.0;
 221   fExBOn              = 0;
 222   fOmegaTau           = 0.0;
 223   fPRFOn              = 0;
 224   fTRFOn              = 0;
 225   fCTOn               = 0;
 226   fTCOn               = 0;
 227   fDriftVelocity      = 0.0;
 228   fPadCoupling        = 0.0;
 229   fTimeCoupling       = 0.0;
 230   fTimeBinWidth       = 0.0;
 231   fField              = 0.0;
 232   fTiltingAngle       = 0.0;
 233
 234   fPRFbin             = 0;
 235   fPRFlo              = 0.0;
 236   fPRFhi              = 0.0;
 237   fPRFwid             = 0.0;
 238   fPRFpad             = 0;
 239   fTRFbin             = 0;
 240   fTRFlo              = 0.0;
 241   fTRFhi              = 0.0;
 242   fTRFwid             = 0.0;
 243   fTCnexp             = 0;
 244
 245   fCompress           = kTRUE;
 246   fDebug              = 0;
 247   fSDigits            = kFALSE;
 248   fSDigitsScale       = 0.0;
 249
 250 }
 251
 252 //_____________________________________________________________________________
 253 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
 254                 :AliDigitizer(name,title)
 255 {
 256   //
 257   // AliTRDdigitizer constructor
 258   //
 259
 260   fInputFile          = NULL;
 261
 262   fDigitsManager      = NULL;
 263   fSDigitsManager     = NULL;
 264   fSDigitsManagerList = NULL;
 265
 266   fTRD                = NULL;
 267   fGeo                = NULL;
 268   fPRFsmp             = NULL;
 269   fTRFsmp             = NULL;
 270   fCTsmp              = NULL;
 271
 272   fMasks              = 0;
 273
 274   fEvent              = 0;
 275
 276   fCompress           = kTRUE;
 277   fDebug              = 0;
 278   fSDigits            = kFALSE;
 279
 280   Init();
 281
 282 }
 283
 284 //_____________________________________________________________________________
 285 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager
 286                                 , const Text_t *name, const Text_t *title)
 287                 :AliDigitizer(manager,name,title)
 288 {
 289   //
 290   // AliTRDdigitizer constructor
 291   //
 292
 293   fInputFile          = NULL;
 294
 295   fDigitsManager      = NULL;
 296   fSDigitsManager     = NULL;
 297   fSDigitsManagerList = NULL;
 298
 299   fTRD                = NULL;
 300   fGeo                = NULL;
 301   fPRFsmp             = NULL;
 302   fTRFsmp             = NULL;
 303   fCTsmp              = NULL;
 304
 305   fMasks              = 0;
 306
 307   fEvent              = 0;
 308
 309   fCompress           = kTRUE;
 310   fDebug              = 0;
 311   fSDigits            = kFALSE;
 312
 313   Init();
 314
 315 }
 316
 317 //_____________________________________________________________________________
 318 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager)
 319                 :AliDigitizer(manager,"AliTRDdigitizer","TRD digitizer")
 320 {
 321   //
 322   // AliTRDdigitizer constructor
 323   //
 324
 325   fInputFile          = NULL;
 326
 327   fDigitsManager      = NULL;
 328   fSDigitsManager     = NULL;
 329   fSDigitsManagerList = NULL;
 330
 331   fTRD                = NULL;
 332   fGeo                = NULL;
 333   fPRFsmp             = NULL;
 334   fTRFsmp             = NULL;
 335   fCTsmp              = NULL;
 336
 337   fMasks              = 0;
 338
 339   fEvent              = 0;
 340
 341   fCompress           = kTRUE;
 342   fDebug              = 0;
 343   fSDigits            = kFALSE;
 344
 345   Init();
 346
 347 }
 348
 349 //_____________________________________________________________________________
 350 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
 351 {
 352   //
 353   // AliTRDdigitizer copy constructor
 354   //
 355
 356   ((AliTRDdigitizer &) d).Copy(*this);
 357
 358 }
 359
 360 //_____________________________________________________________________________
 361 AliTRDdigitizer::~AliTRDdigitizer()
 362 {
 363   //
 364   // AliTRDdigitizer destructor
 365   //
 366
 367   if (fInputFile) {
 368     fInputFile->Close();
 369     delete fInputFile;
 370     fInputFile = NULL;
 371   }
 372
 373   if (fDigitsManager) {
 374     delete fDigitsManager;
 375     fDigitsManager = NULL;
 376   }
 377
 378   if (fSDigitsManager) {
 379     delete fSDigitsManager;
 380     fSDigitsManager = NULL;
 381   }
 382
 383   if (fSDigitsManagerList) {
 384     delete fSDigitsManagerList;
 385     fSDigitsManagerList = NULL;
 386   }
 387
 388   if (fTRFsmp) {
 389     delete [] fTRFsmp;
 390     fTRFsmp = NULL;
 391   }
 392
 393   if (fPRFsmp) {
 394     delete [] fPRFsmp;
 395     fPRFsmp = NULL;
 396   }
 397
 398   if (fCTsmp) {
 399     delete [] fCTsmp;
 400     fCTsmp  = NULL;
 401   }
 402
 403   if (fMasks) {
 404     delete [] fMasks;
 405     fMasks = 0;
 406   }
 407
 408 }
 409
 410 //_____________________________________________________________________________
 411 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
 412 {
 413   //
 414   // Assignment operator
 415   //
 416
 417   if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
 418   return *this;
 419
 420 }
 421
 422 //_____________________________________________________________________________
 423 void AliTRDdigitizer::Copy(TObject &d)
 424 {
 425   //
 426   // Copy function
 427   //
 428
 429   Int_t iBin;
 430
 431   ((AliTRDdigitizer &) d).fInputFile          = NULL;
 432   ((AliTRDdigitizer &) d).fSDigitsManagerList = NULL;
 433   ((AliTRDdigitizer &) d).fSDigitsManager     = NULL;
 434   ((AliTRDdigitizer &) d).fDigitsManager      = NULL;
 435   ((AliTRDdigitizer &) d).fTRD                = NULL;
 436   ((AliTRDdigitizer &) d).fGeo                = NULL;
 437
 438   ((AliTRDdigitizer &) d).fMasks              = 0;
 439
 440   ((AliTRDdigitizer &) d).fEvent              = 0;
 441
 442   ((AliTRDdigitizer &) d).fGasGain            = fGasGain;
 443   ((AliTRDdigitizer &) d).fNoise              = fNoise;
 444   ((AliTRDdigitizer &) d).fChipGain           = fChipGain;
 445   ((AliTRDdigitizer &) d).fADCoutRange        = fADCoutRange;
 446   ((AliTRDdigitizer &) d).fADCinRange         = fADCinRange;
 447   ((AliTRDdigitizer &) d).fADCthreshold       = fADCthreshold;
 448   ((AliTRDdigitizer &) d).fDiffusionOn        = fDiffusionOn;
 449   ((AliTRDdigitizer &) d).fDiffusionT         = fDiffusionT;
 450   ((AliTRDdigitizer &) d).fDiffusionL         = fDiffusionL;
 451   ((AliTRDdigitizer &) d).fElAttachOn         = fElAttachOn;
 452   ((AliTRDdigitizer &) d).fElAttachProp       = fElAttachProp;
 453   ((AliTRDdigitizer &) d).fExBOn              = fExBOn;
 454   ((AliTRDdigitizer &) d).fOmegaTau           = fOmegaTau;
 455   ((AliTRDdigitizer &) d).fLorentzFactor      = fLorentzFactor;
 456   ((AliTRDdigitizer &) d).fDriftVelocity      = fDriftVelocity;
 457   ((AliTRDdigitizer &) d).fPadCoupling        = fPadCoupling;
 458   ((AliTRDdigitizer &) d).fTimeCoupling       = fTimeCoupling;
 459   ((AliTRDdigitizer &) d).fTimeBinWidth       = fTimeBinWidth;
 460   ((AliTRDdigitizer &) d).fField              = fField;
 461   ((AliTRDdigitizer &) d).fPRFOn              = fPRFOn;
 462   ((AliTRDdigitizer &) d).fTRFOn              = fTRFOn;
 463   ((AliTRDdigitizer &) d).fCTOn               = fCTOn;
 464   ((AliTRDdigitizer &) d).fTCOn               = fTCOn;
 465   ((AliTRDdigitizer &) d).fTiltingAngle       = fTiltingAngle;
 466
 467   ((AliTRDdigitizer &) d).fCompress           = fCompress;
 468   ((AliTRDdigitizer &) d).fDebug              = fDebug  ;
 469   ((AliTRDdigitizer &) d).fSDigits            = fSDigits;
 470   ((AliTRDdigitizer &) d).fSDigitsScale       = fSDigitsScale;
 471
 472   ((AliTRDdigitizer &) d).fPRFbin             = fPRFbin;
 473   ((AliTRDdigitizer &) d).fPRFlo              = fPRFlo;
 474   ((AliTRDdigitizer &) d).fPRFhi              = fPRFhi;
 475   ((AliTRDdigitizer &) d).fPRFwid             = fPRFwid;
 476   ((AliTRDdigitizer &) d).fPRFpad             = fPRFpad;
 477   if (((AliTRDdigitizer &) d).fPRFsmp) delete [] ((AliTRDdigitizer &) d).fPRFsmp;
 478   ((AliTRDdigitizer &) d).fPRFsmp = new Float_t[fPRFbin];
 479   for (iBin = 0; iBin < fPRFbin; iBin++) {
 480     ((AliTRDdigitizer &) d).fPRFsmp[iBin] = fPRFsmp[iBin];
 481   }
 482   ((AliTRDdigitizer &) d).fTRFbin             = fTRFbin;
 483   ((AliTRDdigitizer &) d).fTRFlo              = fTRFlo;
 484   ((AliTRDdigitizer &) d).fTRFhi              = fTRFhi;
 485   ((AliTRDdigitizer &) d).fTRFwid             = fTRFwid;
 486   if (((AliTRDdigitizer &) d).fTRFsmp) delete [] ((AliTRDdigitizer &) d).fTRFsmp;
 487   ((AliTRDdigitizer &) d).fTRFsmp = new Float_t[fTRFbin];
 488   for (iBin = 0; iBin < fTRFbin; iBin++) {
 489     ((AliTRDdigitizer &) d).fTRFsmp[iBin] = fTRFsmp[iBin];
 490   }
 491
 492   if (((AliTRDdigitizer &) d).fCTsmp)  delete [] ((AliTRDdigitizer &) d).fCTsmp;
 493   ((AliTRDdigitizer &) d).fCTsmp  = new Float_t[fTRFbin];
 494   for (iBin = 0; iBin < fTRFbin; iBin++) {
 495     ((AliTRDdigitizer &) d).fCTsmp[iBin]  = fCTsmp[iBin];
 496   }
 497
 498   ((AliTRDdigitizer &) d).fTCnexp             = fTCnexp;
 499
 500 }
 501
 502 //_____________________________________________________________________________
 503 Float_t AliTRDdigitizer::CrossTalk(Float_t time)
 504 {
 505   //
 506   // Applies the pad-pad capacitive cross talk
 507   //
 508
 509   Int_t iBin = ((Int_t) ((time - fTRFlo) / fTRFwid));
 510   if ((iBin >= 0) && (iBin < fTRFbin)) {
 511     return fCTsmp[iBin];
 512   }
 513   else {
 514     return 0.0;
 515   }
 516
 517 }
 518
 519 //_____________________________________________________________________________
 520 Int_t AliTRDdigitizer::Diffusion(Float_t driftlength, Float_t *xyz)
 521 {
 522   //
 523   // Applies the diffusion smearing to the position of a single electron
 524   //
 525
 526   Float_t driftSqrt = TMath::Sqrt(driftlength);
 527   Float_t sigmaT = driftSqrt * fDiffusionT;
 528   Float_t sigmaL = driftSqrt * fDiffusionL;
 529   xyz[0] = gRandom->Gaus(xyz[0], sigmaL * fLorentzFactor);
 530   xyz[1] = gRandom->Gaus(xyz[1], sigmaT * fLorentzFactor);
 531   xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
 532
 533   return 1;
 534
 535 }
 536
 537 //_____________________________________________________________________________
 538 Int_t AliTRDdigitizer::ExB(Float_t driftlength, Float_t *xyz)
 539 {
 540   //
 541   // Applies E x B effects to the position of a single electron
 542   //
 543
 544   xyz[0] = xyz[0];
 545   xyz[1] = xyz[1] + fOmegaTau * driftlength;
 546   xyz[2] = xyz[2];
 547
 548   return 1;
 549
 550 }
 551
 552 //_____________________________________________________________________________
 553 Int_t AliTRDdigitizer::PadResponse(Float_t signal, Float_t dist
 554                                  , Int_t plane, Float_t *pad)
 555 {
 556   //
 557   // Applies the pad response
 558   //
 559
 560   const Int_t kNplan = AliTRDgeometry::kNplan;
 561
 562   Int_t iBin  = ((Int_t) (( - dist - fPRFlo) / fPRFwid));
 563   Int_t iOff  = plane * fPRFbin;
 564
 565   Int_t iBin0 = iBin - fPRFpad + iOff;
 566   Int_t iBin1 = iBin           + iOff;
 567   Int_t iBin2 = iBin + fPRFpad + iOff;
 568
 569   if ((iBin0 >= 0) && (iBin2 < (fPRFbin*kNplan))) {
 570
 571     pad[0] = signal * fPRFsmp[iBin0];
 572     pad[1] = signal * fPRFsmp[iBin1];
 573     pad[2] = signal * fPRFsmp[iBin2];
 574
 575     return 1;
 576
 577   }
 578   else {
 579
 580     return 0;
 581
 582   }
 583
 584 }
 585
 586 //_____________________________________________________________________________
 587 Float_t AliTRDdigitizer::TimeResponse(Float_t time)
 588 {
 589   //
 590   // Applies the preamp shaper time response
 591   //
 592
 593   Int_t iBin = ((Int_t) ((time - fTRFlo) / fTRFwid));
 594   if ((iBin >= 0) && (iBin < fTRFbin)) {
 595     return fTRFsmp[iBin];
 596   }
 597   else {
 598     return 0.0;
 599   }
 600
 601 }
 602
 603 //_____________________________________________________________________________
 604 Float_t AliTRDdigitizer::Col0Tilted(Float_t col0, Float_t rowOffset
 605                                   , Int_t plane)
 606 {
 607   //
 608   // Calculates col0 for tilted pads
 609   //
 610
 611   Float_t diff = fTiltingAngle * rowOffset;
 612   return (col0 + TMath::Power(-1.0,plane) * diff);
 613
 614 }
 615
 616 //_____________________________________________________________________________
 617 void AliTRDdigitizer::Exec(Option_t* option)
 618 {
 619   //
 620   // Executes the merging
 621   //
 622
 623   Int_t iInput;
 624
 625   AliTRDdigitsManager *sdigitsManager;
 626
 627   TString optionString = option;
 628   if (optionString.Contains("deb")) {
 629     fDebug = 1;
 630     if (optionString.Contains("2")) {
 631       fDebug = 2;
 632     }
 633     printf("<AliTRDdigitizer::Exec> ");
 634     printf("Called with debug option %d\n",fDebug);
 635   }
 636
 637   Int_t nInput = fManager->GetNinputs();
 638   fMasks = new Int_t[nInput];
 639   for (iInput = 0; iInput < nInput; iInput++) {
 640     fMasks[iInput] = fManager->GetMask(iInput);
 641   }
 642
 643   // Set the event number
 644   fEvent = gAlice->GetEvNumber();
 645
 646   // Initialization
 647   InitDetector();
 648
 649   for (iInput = 0; iInput < nInput; iInput++) {
 650
 651     if (fDebug > 0) {
 652       printf("<AliTRDdigitizer::Exec> ");
 653       printf("Add input stream %d\n",iInput);
 654     }
 655
 656     // Read the s-digits via digits manager
 657     sdigitsManager = new AliTRDdigitsManager();
 658     sdigitsManager->SetDebug(fDebug);
 659     sdigitsManager->SetSDigits(kTRUE);
 660     sdigitsManager->ReadDigits(fManager->GetInputTreeTRDS(iInput));
 661
 662     // Add the s-digits to the input list
 663     AddSDigitsManager(sdigitsManager);
 664
 665   }
 666
 667   // Convert the s-digits to normal digits
 668   if (fDebug > 0) {
 669     printf("<AliTRDdigitizer::Exec> ");
 670     printf("Do the conversion\n");
 671   }
 672   SDigits2Digits();
 673
 674   // Store the digits
 675   if (fDebug > 0) {
 676     printf("<AliTRDdigitizer::Exec> ");
 677     printf("Write the digits\n");
 678   }
 679   fDigitsManager->MakeBranch(fManager->GetTreeDTRD());
 680   fDigitsManager->WriteDigits();
 681   if (fDebug > 0) {
 682     printf("<AliTRDdigitizer::Exec> ");
 683     printf("Done\n");
 684   }
 685
 686   DeleteSDigitsManager();
 687
 688 }
 689
 690 //_____________________________________________________________________________
 691 Bool_t AliTRDdigitizer::Init()
 692 {
 693   //
 694   // Initializes the digitization procedure with standard values
 695   //
 696
 697   // The default parameter for the digitization
 698   fGasGain        = 2800.;
 699   fChipGain       = 6.1;
 700   fNoise          = 1000.;
 701   fADCoutRange    = 1023.;          // 10-bit ADC
 702   fADCinRange     = 1000.;          // 1V input range
 703   fADCthreshold   = 1;
 704
 705   // For the summable digits
 706   fSDigitsScale   = 100.;
 707
 708   // The drift velocity (cm / mus)
 709   fDriftVelocity  = 1.5;
 710
 711   // Diffusion on
 712   fDiffusionOn    = 1;
 713
 714   // E x B effects
 715   fExBOn          = 0;
 716
 717   // Propability for electron attachment
 718   fElAttachOn     = 0;
 719   fElAttachProp   = 0.0;
 720
 721   // The pad response function
 722   fPRFOn          = 1;
 723
 724   // The time response function
 725   fTRFOn          = 1;
 726
 727   // The cross talk
 728   fCTOn           = 0;
 729
 730   // The tail cancelation
 731   fTCOn           = 1;
 732
 733   // The number of exponentials
 734   fTCnexp         = 2;
 735
 736   // The pad coupling factor (same number as for the TPC)
 737   fPadCoupling    = 0.5;
 738
 739   // The time coupling factor (same number as for the TPC)
 740   fTimeCoupling   = 0.4;
 741
 742   // The tilting angle for the readout pads
 743   SetTiltingAngle(5.0);
 744
 745   return kTRUE;
 746
 747 }
 748
 749 //_____________________________________________________________________________
 750 Bool_t AliTRDdigitizer::ReInit()
 751 {
 752   //
 753   // Reinitializes the digitization procedure after a change in the parameter
 754   //
 755
 756   if (!fGeo) {
 757     printf("AliTRDdigitizer::ReInit -- ");
 758     printf("No geometry defined. Run InitDetector() first\n");
 759     return kFALSE;
 760   }
 761
 762   // Calculate the time bin width in ns
 763   fTimeBinWidth   = fGeo->GetTimeBinSize() / fDriftVelocity * 1000.0;
 764
 765   // The range and the binwidth for the sampled TRF
 766   fTRFbin = 100;
 767   // Start 0.2 mus before the signal
 768   fTRFlo  = -0.2 * fDriftVelocity;
 769   // End the maximum driftlength after the signal
 770   fTRFhi  = AliTRDgeometry::DrThick()
 771           + fGeo->GetTimeAfter() * fGeo->GetTimeBinSize();
 772   fTRFwid = (fTRFhi - fTRFlo) / ((Float_t) fTRFbin);
 773
 774   // Transverse and longitudinal diffusion coefficients (Xe/CO2)
 775   fDiffusionT     = GetDiffusionT(fDriftVelocity,fField);
 776   fDiffusionL     = GetDiffusionL(fDriftVelocity,fField);
 777
 778   // omega * tau.= tan(Lorentz-angle)
 779   fOmegaTau       = GetOmegaTau(fDriftVelocity,fField);
 780
 781   // The Lorentz factor
 782   if (fExBOn) {
 783     fLorentzFactor = 1.0 / (1.0 + fOmegaTau*fOmegaTau);
 784   }
 785   else {
 786     fLorentzFactor = 1.0;
 787   }
 788
 789   return kTRUE;
 790
 791 }
 792
 793 //_____________________________________________________________________________
 794 void AliTRDdigitizer::SampleTRF()
 795 {
 796   //
 797   // Samples the time response function
 798   //
 799   // New TRF from Venelin Angelov, simulated with CADENCE
 800   // Pad-ground capacitance = 25 pF
 801   // Pad-pad cross talk capacitance = 6 pF
 802   //
 803
 804   Int_t   ipos1;
 805   Int_t   ipos2;
 806   Float_t diff;
 807
 808   const Int_t kNpasa     = 252;
 809
 810   Float_t time[kNpasa]   = { -0.220000, -0.210000, -0.200000, -0.190000
 811                            , -0.180000, -0.170000, -0.160000, -0.150000
 812                            , -0.140000, -0.130000, -0.120000, -0.110000
 813                            , -0.100000, -0.090000, -0.080000, -0.070000
 814                            , -0.060000, -0.050000, -0.040000, -0.030000
 815                            , -0.020000, -0.010000, -0.000000,  0.010000
 816                            ,  0.020000,  0.030000,  0.040000,  0.050000
 817                            ,  0.060000,  0.070000,  0.080000,  0.090000
 818                            ,  0.100000,  0.110000,  0.120000,  0.130000
 819                            ,  0.140000,  0.150000,  0.160000,  0.170000
 820                            ,  0.180000,  0.190000,  0.200000,  0.210000
 821                            ,  0.220000,  0.230000,  0.240000,  0.250000
 822                            ,  0.260000,  0.270000,  0.280000,  0.290000
 823                            ,  0.300000,  0.310000,  0.320000,  0.330000
 824                            ,  0.340000,  0.350000,  0.360000,  0.370000
 825                            ,  0.380000,  0.390000,  0.400000,  0.410000
 826                            ,  0.420000,  0.430000,  0.440000,  0.450000
 827                            ,  0.460000,  0.470000,  0.480000,  0.490000
 828                            ,  0.500000,  0.510000,  0.520000,  0.530000
 829                            ,  0.540000,  0.550000,  0.560000,  0.570000
 830                            ,  0.580000,  0.590000,  0.600000,  0.610000
 831                            ,  0.620000,  0.630000,  0.640000,  0.650000
 832                            ,  0.660000,  0.670000,  0.680000,  0.690000
 833                            ,  0.700000,  0.710000,  0.720000,  0.730000
 834                            ,  0.740000,  0.750000,  0.760000,  0.770000
 835                            ,  0.780000,  0.790000,  0.800000,  0.810000
 836                            ,  0.820000,  0.830000,  0.840000,  0.850000
 837                            ,  0.860000,  0.870000,  0.880000,  0.890000
 838                            ,  0.900000,  0.910000,  0.920000,  0.930000
 839                            ,  0.940000,  0.950000,  0.960000,  0.970000
 840                            ,  0.980000,  0.990000,  1.000000,  1.010000
 841                            ,  1.020000,  1.030000,  1.040000,  1.050000
 842                            ,  1.060000,  1.070000,  1.080000,  1.090000
 843                            ,  1.100000,  1.110000,  1.120000,  1.130000
 844                            ,  1.140000,  1.150000,  1.160000,  1.170000
 845                            ,  1.180000,  1.190000,  1.200000,  1.210000
 846                            ,  1.220000,  1.230000,  1.240000,  1.250000
 847                            ,  1.260000,  1.270000,  1.280000,  1.290000
 848                            ,  1.300000,  1.310000,  1.320000,  1.330000
 849                            ,  1.340000,  1.350000,  1.360000,  1.370000
 850                            ,  1.380000,  1.390000,  1.400000,  1.410000
 851                            ,  1.420000,  1.430000,  1.440000,  1.450000
 852                            ,  1.460000,  1.470000,  1.480000,  1.490000
 853                            ,  1.500000,  1.510000,  1.520000,  1.530000
 854                            ,  1.540000,  1.550000,  1.560000,  1.570000
 855                            ,  1.580000,  1.590000,  1.600000,  1.610000
 856                            ,  1.620000,  1.630000,  1.640000,  1.650000
 857                            ,  1.660000,  1.670000,  1.680000,  1.690000
 858                            ,  1.700000,  1.710000,  1.720000,  1.730000
 859                            ,  1.740000,  1.750000,  1.760000,  1.770000
 860                            ,  1.780000,  1.790000,  1.800000,  1.810000
 861                            ,  1.820000,  1.830000,  1.840000,  1.850000
 862                            ,  1.860000,  1.870000,  1.880000,  1.890000
 863                            ,  1.900000,  1.910000,  1.920000,  1.930000
 864                            ,  1.940000,  1.950000,  1.960000,  1.970000
 865                            ,  1.980000,  1.990000,  2.000000,  2.010000
 866                            ,  2.020000,  2.030000,  2.040000,  2.050000
 867                            ,  2.060000,  2.070000,  2.080000,  2.090000
 868                            ,  2.100000,  2.110000,  2.120000,  2.130000
 869                            ,  2.140000,  2.150000,  2.160000,  2.170000
 870                            ,  2.180000,  2.190000,  2.200000,  2.210000
 871                            ,  2.220000,  2.230000,  2.240000,  2.250000
 872                            ,  2.260000,  2.270000,  2.280000,  2.290000 };
 873
 874   Float_t signal[kNpasa] = {  0.000000,  0.000000,  0.000000,  0.000000
 875                            ,  0.000000,  0.000000,  0.000000,  0.000396
 876                            ,  0.005096,  0.022877,  0.061891,  0.126614
 877                            ,  0.215798,  0.324406,  0.444507,  0.566817
 878                            ,  0.683465,  0.787089,  0.873159,  0.937146
 879                            ,  0.979049,  0.999434,  1.000000,  0.983579
 880                            ,  0.954134,  0.913364,  0.866365,  0.813703
 881                            ,  0.759910,  0.706116,  0.653454,  0.603624
 882                            ,  0.556625,  0.514156,  0.475085,  0.439977
 883                            ,  0.408834,  0.380578,  0.355549,  0.333352
 884                            ,  0.313647,  0.296093,  0.280351,  0.266195
 885                            ,  0.253397,  0.241789,  0.231257,  0.221574
 886                            ,  0.212627,  0.204417,  0.196772,  0.189581
 887                            ,  0.182956,  0.176784,  0.171008,  0.165515
 888                            ,  0.160419,  0.155606,  0.151076,  0.146716
 889                            ,  0.142639,  0.138845,  0.135221,  0.131767
 890                            ,  0.128482,  0.125368,  0.122424,  0.119592
 891                            ,  0.116931,  0.114326,  0.111891,  0.109513
 892                            ,  0.107248,  0.105096,  0.103058,  0.101019
 893                            ,  0.099151,  0.097282,  0.095527,  0.093715
 894                            ,  0.092129,  0.090544,  0.088958,  0.087429
 895                            ,  0.086014,  0.084598,  0.083239,  0.081880
 896                            ,  0.080634,  0.079388,  0.078143,  0.077010
 897                            ,  0.075878,  0.074745,  0.073669,  0.072593
 898                            ,  0.071574,  0.070612,  0.069649,  0.068686
 899                            ,  0.067780,  0.066874,  0.066025,  0.065176
 900                            ,  0.064326,  0.063533,  0.062684,  0.061948
 901                            ,  0.061212,  0.060419,  0.059740,  0.059003
 902                            ,  0.058324,  0.057644,  0.057022,  0.056342
 903                            ,  0.055663,  0.055096,  0.054473,  0.053851
 904                            ,  0.053284,  0.052718,  0.052152,  0.051585
 905                            ,  0.051019,  0.050566,  0.050000,  0.049490
 906                            ,  0.048981,  0.048528,  0.048018,  0.047508
 907                            ,  0.047055,  0.046602,  0.046149,  0.045696
 908                            ,  0.045300,  0.044904,  0.044451,  0.044054
 909                            ,  0.043658,  0.043205,  0.042865,  0.042469
 910                            ,  0.042072,  0.041733,  0.041336,  0.040997
 911                            ,  0.040657,  0.040260,  0.039921,  0.039581
 912                            ,  0.039241,  0.038958,  0.038618,  0.038335
 913                            ,  0.037995,  0.037656,  0.037373,  0.037089
 914                            ,  0.036806,  0.036467,  0.036183,  0.035900
 915                            ,  0.035617,  0.035334,  0.035108,  0.034824
 916                            ,  0.034541,  0.034315,  0.034032,  0.033805
 917                            ,  0.033522,  0.033296,  0.033069,  0.032786
 918                            ,  0.032559,  0.032333,  0.032106,  0.031880
 919                            ,  0.031653,  0.031427,  0.031200,  0.030974
 920                            ,  0.030804,  0.030578,  0.030351,  0.030125
 921                            ,  0.029955,  0.029785,  0.029558,  0.029332
 922                            ,  0.029162,  0.028992,  0.028766,  0.028596
 923                            ,  0.028426,  0.028199,  0.028086,  0.027860
 924                            ,  0.027746,  0.027633,  0.027463,  0.027293
 925                            ,  0.027180,  0.027067,  0.026954,  0.026954
 926                            ,  0.026840,  0.026727,  0.026727,  0.026614
 927                            ,  0.026614,  0.026614,  0.026557,  0.026501
 928                            ,  0.026501,  0.026501,  0.026501,  0.026501
 929                            ,  0.026501,  0.026501,  0.026501,  0.026387
 930                            ,  0.026387,  0.026387,  0.026387,  0.026387
 931                            ,  0.026387,  0.026387,  0.026387,  0.026387
 932                            ,  0.026387,  0.026387,  0.026387,  0.026387
 933                            ,  0.026387,  0.026274,  0.026274,  0.026274
 934                            ,  0.026274,  0.026274,  0.026274,  0.026274
 935                            ,  0.026274,  0.026274,  0.026274,  0.026274
 936                            ,  0.026274,  0.026274,  0.026274,  0.026161 };
 937
 938   Float_t xtalk[kNpasa]  = {  0.000000,  0.000000,  0.000000,  0.000000
 939                            ,  0.000000,  0.000000,  0.000000,  0.000113
 940                            ,  0.000793,  0.003058,  0.007305,  0.013194
 941                            ,  0.019706,  0.025821,  0.030634,  0.033465
 942                            ,  0.034145,  0.032729,  0.029615,  0.025198
 943                            ,  0.019989,  0.014496,  0.009003,  0.003964
 944                            , -0.000510, -0.004190, -0.007191, -0.009400
 945                            , -0.010872, -0.011835, -0.012288, -0.012288
 946                            , -0.012005, -0.011495, -0.010872, -0.010136
 947                            , -0.009343, -0.008607, -0.007871, -0.007191
 948                            , -0.006512, -0.005946, -0.005379, -0.004926
 949                            , -0.004473, -0.004077, -0.003737, -0.003398
 950                            , -0.003114, -0.002831, -0.002605, -0.002378
 951                            , -0.002208, -0.002039, -0.001869, -0.001699
 952                            , -0.001585, -0.001472, -0.001359, -0.001246
 953                            , -0.001132, -0.001019, -0.001019, -0.000906
 954                            , -0.000906, -0.000793, -0.000793, -0.000680
 955                            , -0.000680, -0.000680, -0.000566, -0.000566
 956                            , -0.000566, -0.000566, -0.000453, -0.000453
 957                            , -0.000453, -0.000453, -0.000453, -0.000453
 958                            , -0.000340, -0.000340, -0.000340, -0.000340
 959                            , -0.000340, -0.000340, -0.000340, -0.000340
 960                            , -0.000340, -0.000340, -0.000340, -0.000340
 961                            , -0.000340, -0.000227, -0.000227, -0.000227
 962                            , -0.000227, -0.000227, -0.000227, -0.000227
 963                            , -0.000227, -0.000227, -0.000227, -0.000227
 964                            , -0.000227, -0.000227, -0.000227, -0.000227
 965                            , -0.000227, -0.000227, -0.000227, -0.000227
 966                            , -0.000227, -0.000227, -0.000227, -0.000227
 967                            , -0.000227, -0.000227, -0.000227, -0.000227
 968                            , -0.000227, -0.000227, -0.000227, -0.000227
 969                            , -0.000227, -0.000227, -0.000227, -0.000227
 970                            , -0.000227, -0.000227, -0.000227, -0.000113
 971                            , -0.000113, -0.000113, -0.000113, -0.000113
 972                            , -0.000113, -0.000113, -0.000113, -0.000113
 973                            , -0.000113, -0.000113, -0.000113, -0.000113
 974                            , -0.000113, -0.000113, -0.000113, -0.000113
 975                            , -0.000113, -0.000113, -0.000113, -0.000113
 976                            , -0.000113, -0.000113, -0.000113, -0.000113
 977                            , -0.000113, -0.000113, -0.000113, -0.000113
 978                            , -0.000113, -0.000113, -0.000113, -0.000113
 979                            , -0.000113, -0.000113, -0.000113, -0.000113
 980                            , -0.000113, -0.000113, -0.000113, -0.000113
 981                            , -0.000113, -0.000113, -0.000113, -0.000113
 982                            , -0.000113, -0.000113, -0.000113, -0.000113
 983                            , -0.000113, -0.000113, -0.000113, -0.000113
 984                            , -0.000113, -0.000113, -0.000113, -0.000113
 985                            , -0.000113, -0.000113, -0.000113, -0.000113
 986                            , -0.000113, -0.000113, -0.000113, -0.000113
 987                            , -0.000113, -0.000113, -0.000113, -0.000113
 988                            , -0.000113, -0.000113, -0.000113, -0.000113
 989                            , -0.000113, -0.000113, -0.000113, -0.000113
 990                            , -0.000113, -0.000113, -0.000113,  0.000000
 991                            ,  0.000000,  0.000000,  0.000000,  0.000000
 992                            ,  0.000000,  0.000000,  0.000000,  0.000000
 993                            ,  0.000000,  0.000000,  0.000000,  0.000000
 994                            ,  0.000000,  0.000000,  0.000000,  0.000000
 995                            ,  0.000000,  0.000000,  0.000000,  0.000000
 996                            ,  0.000000,  0.000000,  0.000000,  0.000000
 997                            ,  0.000000,  0.000000,  0.000000,  0.000000
 998                            ,  0.000000,  0.000000,  0.000000,  0.000000
 999                            ,  0.000000,  0.000000,  0.000000,  0.000000
1000                            ,  0.000000,  0.000000,  0.000000,  0.000000 };
1001
1002   // increase CrossTalk to measurements
1003   for (Int_t ipasa = 0; ipasa < kNpasa; ipasa++) {
1004     xtalk[ipasa] *= 1.75;
1005   }
1006
1007   if (fTRFsmp) delete [] fTRFsmp;
1008   fTRFsmp = new Float_t[fTRFbin];
1009   if (fCTsmp)  delete [] fCTsmp;
1010   fCTsmp  = new Float_t[fTRFbin];
1011
1012   Float_t loTRF    = TMath::Max(fTRFlo / fDriftVelocity,time[0]);
1013   Float_t hiTRF    = TMath::Min(fTRFhi / fDriftVelocity,time[kNpasa-1]);
1014   Float_t binWidth = (hiTRF - loTRF) / ((Float_t) fTRFbin);
1015
1016   // Take the linear interpolation
1017   for (Int_t iBin = 0; iBin < fTRFbin; iBin++) {
1018
1019     Float_t bin = (((Float_t) iBin) + 0.5) * binWidth + loTRF;
1020     ipos1 = ipos2 = 0;
1021     diff  = 0;
1022     do {
1023       diff = bin - time[ipos2++];
1024     } while (diff > 0);
1025     ipos2--;
1026     if (ipos2 >= kNpasa) ipos2 = kNpasa - 1;
1027     ipos1 = ipos2 - 1;
1028
1029     fTRFsmp[iBin] = signal[ipos2]
1030                   + diff * (signal[ipos2] - signal[ipos1])
1031                          / (  time[ipos2] -   time[ipos1]);
1032
1033     fCTsmp[iBin]  = xtalk[ipos2]
1034                   + diff * (xtalk[ipos2]  -  xtalk[ipos1])
1035                          / (  time[ipos2] -   time[ipos1]);
1036
1037   }
1038
1039 }
1040
1041 //_____________________________________________________________________________
1042 void AliTRDdigitizer::SamplePRF()
1043 {
1044   //
1045   // Samples the pad response function
1046   //
1047
1048   const Int_t kNplan  = AliTRDgeometry::kNplan;
1049   const Int_t kPRFbin = 61;
1050
1051   Float_t prf[kNplan][kPRFbin] = { { 0.018570, 0.022270, 0.026710, 0.032010
1052                                    , 0.038350, 0.045920, 0.054930, 0.065650
1053                                    , 0.078370, 0.093420, 0.111150, 0.131940
1054                                    , 0.156160, 0.184160, 0.216220, 0.252470
1055                                    , 0.292860, 0.337030, 0.384330, 0.433750
1056                                    , 0.484010, 0.533630, 0.581150, 0.625200
1057                                    , 0.664710, 0.698860, 0.727130, 0.749230
1058                                    , 0.765050, 0.774540, 0.777700, 0.774540
1059                                    , 0.765050, 0.749230, 0.727130, 0.698860
1060                                    , 0.664710, 0.625200, 0.581150, 0.533630
1061                                    , 0.484010, 0.433750, 0.384330, 0.337030
1062                                    , 0.292860, 0.252470, 0.216220, 0.184160
1063                                    , 0.156160, 0.131940, 0.111150, 0.093420
1064                                    , 0.078370, 0.065650, 0.054930, 0.045920
1065                                    , 0.038350, 0.032010, 0.026710, 0.022270
1066                                    , 0.018570                               }
1067                                  , { 0.015730, 0.019040, 0.023030, 0.027840
1068                                    , 0.033650, 0.040650, 0.049060, 0.059160
1069                                    , 0.071260, 0.085710, 0.102910, 0.123270
1070                                    , 0.147240, 0.175220, 0.207590, 0.244540
1071                                    , 0.286090, 0.331920, 0.381350, 0.433290
1072                                    , 0.486290, 0.538710, 0.588870, 0.635280
1073                                    , 0.676760, 0.712460, 0.741890, 0.764810
1074                                    , 0.781150, 0.790930, 0.794180, 0.790930
1075                                    , 0.781150, 0.764810, 0.741890, 0.712460
1076                                    , 0.676760, 0.635280, 0.588870, 0.538710
1077                                    , 0.486290, 0.433290, 0.381350, 0.331920
1078                                    , 0.286090, 0.244540, 0.207590, 0.175220
1079                                    , 0.147240, 0.123270, 0.102910, 0.085710
1080                                    , 0.071260, 0.059160, 0.049060, 0.040650
1081                                    , 0.033650, 0.027840, 0.023030, 0.019040
1082                                    , 0.015730                               }
1083                                  , { 0.013330, 0.016270, 0.019850, 0.024210
1084                                    , 0.029510, 0.035960, 0.043790, 0.053280
1085                                    , 0.064740, 0.078580, 0.095190, 0.115070
1086                                    , 0.138700, 0.166570, 0.199120, 0.236660
1087                                    , 0.279260, 0.326660, 0.378140, 0.432540
1088                                    , 0.488260, 0.543440, 0.596200, 0.644900
1089                                    , 0.688240, 0.725380, 0.755840, 0.779470
1090                                    , 0.796260, 0.806280, 0.809610, 0.806280
1091                                    , 0.796260, 0.779470, 0.755840, 0.725380
1092                                    , 0.688240, 0.644900, 0.596200, 0.543440
1093                                    , 0.488260, 0.432540, 0.378140, 0.326660
1094                                    , 0.279260, 0.236660, 0.199120, 0.166570
1095                                    , 0.138700, 0.115070, 0.095190, 0.078580
1096                                    , 0.064740, 0.053280, 0.043790, 0.035960
1097                                    , 0.029510, 0.024210, 0.019850, 0.016270
1098                                    , 0.013330                               }
1099                                  , { 0.011280, 0.013890, 0.017090, 0.021030
1100                                    , 0.025870, 0.031800, 0.039060, 0.047940
1101                                    , 0.058790, 0.071980, 0.087990, 0.107330
1102                                    , 0.130550, 0.158220, 0.190850, 0.228870
1103                                    , 0.272410, 0.321270, 0.374740, 0.431560
1104                                    , 0.489960, 0.547870, 0.603180, 0.654080
1105                                    , 0.699190, 0.737640, 0.769030, 0.793260
1106                                    , 0.810410, 0.820620, 0.824010, 0.820620
1107                                    , 0.810410, 0.793260, 0.769030, 0.737640
1108                                    , 0.699190, 0.654080, 0.603180, 0.547870
1109                                    , 0.489960, 0.431560, 0.374740, 0.321270
1110                                    , 0.272410, 0.228870, 0.190850, 0.158220
1111                                    , 0.130550, 0.107330, 0.087990, 0.071980
1112                                    , 0.058790, 0.047940, 0.039060, 0.031800
1113                                    , 0.025870, 0.021030, 0.017090, 0.013890
1114                                    , 0.011280                               }
1115                                  , { 0.009550, 0.011860, 0.014720, 0.018270
1116                                    , 0.022660, 0.028100, 0.034820, 0.043120
1117                                    , 0.053340, 0.065900, 0.081280, 0.100040
1118                                    , 0.122800, 0.150180, 0.182800, 0.221170
1119                                    , 0.265550, 0.315790, 0.371180, 0.430370
1120                                    , 0.491430, 0.552030, 0.609840, 0.662860
1121                                    , 0.709630, 0.749290, 0.781490, 0.806220
1122                                    , 0.823650, 0.834000, 0.837430, 0.834000
1123                                    , 0.823650, 0.806220, 0.781490, 0.749290
1124                                    , 0.709630, 0.662860, 0.609840, 0.552030
1125                                    , 0.491430, 0.430370, 0.371180, 0.315790
1126                                    , 0.265550, 0.221170, 0.182800, 0.150180
1127                                    , 0.122800, 0.100040, 0.081280, 0.065900
1128                                    , 0.053340, 0.043120, 0.034820, 0.028100
1129                                    , 0.022660, 0.018270, 0.014720, 0.011860
1130                                    , 0.009550                               }
1131                                  , { 0.008080, 0.010120, 0.012670, 0.015860
1132                                    , 0.019840, 0.024820, 0.031030, 0.038760
1133                                    , 0.048370, 0.060300, 0.075040, 0.093200
1134                                    , 0.115430, 0.142450, 0.174980, 0.213610
1135                                    , 0.258720, 0.310250, 0.367480, 0.429010
1136                                    , 0.492690, 0.555950, 0.616210, 0.671280
1137                                    , 0.719600, 0.760350, 0.793250, 0.818380
1138                                    , 0.836020, 0.846460, 0.849920, 0.846460
1139                                    , 0.836020, 0.818380, 0.793250, 0.760350
1140                                    , 0.719600, 0.671280, 0.616210, 0.555950
1141                                    , 0.492690, 0.429010, 0.367480, 0.310250
1142                                    , 0.258720, 0.213610, 0.174980, 0.142450
1143                                    , 0.115430, 0.093200, 0.075040, 0.060300
1144                                    , 0.048370, 0.038760, 0.031030, 0.024820
1145                                    , 0.019840, 0.015860, 0.012670, 0.010120
1146                                    , 0.008080                               } };
1147
1148   fPRFlo  = -1.5;
1149   fPRFhi  =  1.5;
1150   fPRFbin = kPRFbin;
1151   fPRFwid = (fPRFhi - fPRFlo) / ((Float_t) fPRFbin);
1152   fPRFpad = ((Int_t) (1.0 / fPRFwid));
1153
1154   if (fPRFsmp) delete [] fPRFsmp;
1155   fPRFsmp = new Float_t[kNplan*fPRFbin];
1156   for (Int_t iPla = 0; iPla < kNplan; iPla++) {
1157     for (Int_t iBin = 0; iBin < fPRFbin; iBin++) {
1158       fPRFsmp[iPla*kPRFbin+iBin] = prf[iPla][iBin];
1159     }
1160   }
1161
1162 }
1163
1164 //_____________________________________________________________________________
1165 Bool_t AliTRDdigitizer::Open(const Char_t *file, Int_t nEvent)
1166 {
1167   //
1168   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
1169   //
1170
1171   // Connect the AliRoot file containing Geometry, Kine, and Hits
1172   fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(file);
1173   if (!fInputFile) {
1174     if (fDebug > 0) {
1175       printf("<AliTRDdigitizer::Open> ");
1176       printf("Open the AliROOT-file %s.\n",file);
1177     }
1178     fInputFile = new TFile(file,"UPDATE");
1179   }
1180   else {
1181     if (fDebug > 0) {
1182       printf("<AliTRDdigitizer::Open> ");
1183       printf("%s is already open.\n",file);
1184     }
1185   }
1186
1187   gAlice = (AliRun*) fInputFile->Get("gAlice");
1188   if (gAlice) {
1189     if (fDebug > 0) {
1190       printf("<AliTRDdigitizer::Open> ");
1191       printf("AliRun object found on file.\n");
1192     }
1193   }
1194   else {
1195     printf("<AliTRDdigitizer::Open> ");
1196     printf("Could not find AliRun object.\n");
1197     return kFALSE;
1198   }
1199
1200   fEvent = nEvent;
1201
1202   // Import the Trees for the event nEvent in the file
1203   Int_t nparticles = gAlice->GetEvent(fEvent);
1204   if (nparticles <= 0) {
1205     printf("<AliTRDdigitizer::Open> ");
1206     printf("No entries in the trees for event %d.\n",fEvent);
1207     return kFALSE;
1208   }
1209
1210   if (InitDetector()) {
1211     return MakeBranch();
1212   }
1213   else {
1214     return kFALSE;
1215   }
1216
1217 }
1218
1219 //_____________________________________________________________________________
1220 Bool_t AliTRDdigitizer::InitDetector()
1221 {
1222   //
1223   // Sets the pointer to the TRD detector and the geometry
1224   //
1225
1226   // Get the pointer to the detector class and check for version 1
1227   fTRD = (AliTRD*) gAlice->GetDetector("TRD");
1228   if (fTRD->IsVersion() != 1) {
1229     printf("<AliTRDdigitizer::InitDetector> ");
1230     printf("TRD must be version 1 (slow simulator).\n");
1231     exit(1);
1232   }
1233
1234   // Get the geometry
1235   fGeo = fTRD->GetGeometry();
1236   if (fDebug > 0) {
1237     printf("<AliTRDdigitizer::InitDetector> ");
1238     printf("Geometry version %d\n",fGeo->IsVersion());
1239   }
1240
1241   // The magnetic field strength in Tesla
1242   fField = 0.2 * gAlice->Field()->Factor();
1243
1244   // Create a digits manager
1245   fDigitsManager = new AliTRDdigitsManager();
1246   fDigitsManager->SetSDigits(fSDigits);
1247   fDigitsManager->CreateArrays();
1248   fDigitsManager->SetEvent(fEvent);
1249   fDigitsManager->SetDebug(fDebug);
1250
1251   // The list for the input s-digits manager to be merged
1252   fSDigitsManagerList = new TList();
1253
1254   return ReInit();
1255
1256 }
1257
1258 //_____________________________________________________________________________
1259 Bool_t AliTRDdigitizer::MakeBranch(const Char_t *file)
1260 {
1261   //
1262   // Create the branches for the digits array
1263   //
1264
1265   return fDigitsManager->MakeBranch(file);
1266
1267 }
1268
1269 //_____________________________________________________________________________
1270 Bool_t AliTRDdigitizer::MakeDigits()
1271 {
1272   //
1273   // Creates digits.
1274   //
1275
1276   ///////////////////////////////////////////////////////////////
1277   // Parameter
1278   ///////////////////////////////////////////////////////////////
1279
1280   // Converts number of electrons to fC
1281   const Double_t kEl2fC  = 1.602E-19 * 1.0E15;
1282
1283   ///////////////////////////////////////////////////////////////
1284
1285   // Number of pads included in the pad response
1286   const Int_t kNpad  = 3;
1287
1288   // Number of track dictionary arrays
1289   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1290
1291   // Half the width of the amplification region
1292   const Float_t kAmWidth = AliTRDgeometry::AmThick() / 2.;
1293
1294   Int_t   iRow, iCol, iTime, iPad;
1295   Int_t   iDict  = 0;
1296   Int_t   nBytes = 0;
1297
1298   Int_t   totalSizeDigits = 0;
1299   Int_t   totalSizeDict0  = 0;
1300   Int_t   totalSizeDict1  = 0;
1301   Int_t   totalSizeDict2  = 0;
1302
1303   Int_t   timeTRDbeg = 0;
1304   Int_t   timeTRDend = 1;
1305
1306   Float_t pos[3];
1307   Float_t rot[3];
1308   Float_t xyz[3];
1309   Float_t padSignal[kNpad];
1310   Float_t signalOld[kNpad];
1311
1312   AliTRDdataArrayF *signals = 0;
1313   AliTRDdataArrayI *digits  = 0;
1314   AliTRDdataArrayI *dictionary[kNDict];
1315
1316   // Create a container for the amplitudes
1317   AliTRDsegmentArray *signalsArray
1318                      = new AliTRDsegmentArray("AliTRDdataArrayF"
1319                                              ,AliTRDgeometry::Ndet());
1320
1321   if (fTRFOn) {
1322     timeTRDbeg = ((Int_t) (-fTRFlo / fGeo->GetTimeBinSize())) - 1;
1323     timeTRDend = ((Int_t) ( fTRFhi / fGeo->GetTimeBinSize())) - 1;
1324     if (fDebug > 0) {
1325       printf("<AliTRDdigitizer::MakeDigits> ");
1326       printf("Sample the TRF between -%d and %d\n",timeTRDbeg,timeTRDend);
1327     }
1328   }
1329
1330   Float_t elAttachProp = fElAttachProp / 100.;
1331
1332   // Create the sampled PRF
1333   SamplePRF();
1334
1335   // Create the sampled TRF
1336   SampleTRF();
1337
1338   if (!fGeo) {
1339     printf("<AliTRDdigitizer::MakeDigits> ");
1340     printf("No geometry defined\n");
1341     return kFALSE;
1342   }
1343
1344   if (fDebug > 0) {
1345     printf("<AliTRDdigitizer::MakeDigits> ");
1346     printf("Start creating digits.\n");
1347   }
1348
1349   // Get the pointer to the hit tree
1350   TTree *HitTree = gAlice->TreeH();
1351
1352   // Get the number of entries in the hit tree
1353   // (Number of primary particles creating a hit somewhere)
1354   Int_t nTrack = (Int_t) HitTree->GetEntries();
1355   if (fDebug > 0) {
1356     printf("<AliTRDdigitizer::MakeDigits> ");
1357     printf("Found %d primary particles\n",nTrack);
1358   }
1359
1360   Int_t detectorOld = -1;
1361   Int_t countHits   =  0;
1362
1363   // Loop through all entries in the tree
1364   for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
1365
1366     gAlice->ResetHits();
1367     nBytes += HitTree->GetEvent(iTrack);
1368
1369     // Loop through the TRD hits
1370     Int_t iHit = 0;
1371     AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
1372     while (hit) {
1373
1374       countHits++;
1375       iHit++;
1376
1377               pos[0]      = hit->X();
1378               pos[1]      = hit->Y();
1379               pos[2]      = hit->Z();
1380       Float_t q           = hit->GetCharge();
1381       Int_t   track       = hit->Track();
1382       Int_t   detector    = hit->GetDetector();
1383       Int_t   plane       = fGeo->GetPlane(detector);
1384       Int_t   sector      = fGeo->GetSector(detector);
1385       Int_t   chamber     = fGeo->GetChamber(detector);
1386       Int_t   nRowMax     = fGeo->GetRowMax(plane,chamber,sector);
1387       Int_t   nColMax     = fGeo->GetColMax(plane);
1388       Int_t   nTimeMax    = fGeo->GetTimeMax();
1389       Int_t   nTimeBefore = fGeo->GetTimeBefore();
1390       Int_t   nTimeAfter  = fGeo->GetTimeAfter();
1391       Int_t   nTimeTotal  = fGeo->GetTimeTotal();
1392       Float_t row0        = fGeo->GetRow0(plane,chamber,sector);
1393       Float_t col0        = fGeo->GetCol0(plane);
1394       Float_t time0       = fGeo->GetTime0(plane);
1395       Float_t rowPadSize  = fGeo->GetRowPadSize(plane,chamber,sector);
1396       Float_t colPadSize  = fGeo->GetColPadSize(plane);
1397       Float_t timeBinSize = fGeo->GetTimeBinSize();
1398       Float_t divideRow   = 1.0 / rowPadSize;
1399       Float_t divideCol   = 1.0 / colPadSize;
1400       Float_t divideTime  = 1.0 / timeBinSize;
1401
1402       if (fDebug > 1) {
1403         printf("Analyze hit no. %d ",iHit);
1404         printf("-----------------------------------------------------------\n");
1405         hit->Dump();
1406         printf("plane = %d, sector = %d, chamber = %d\n"
1407               ,plane,sector,chamber);
1408         printf("nRowMax = %d, nColMax = %d, nTimeMax = %d\n"
1409               ,nRowMax,nColMax,nTimeMax);
1410         printf("nTimeBefore = %d, nTimeAfter = %d, nTimeTotal = %d\n"
1411               ,nTimeBefore,nTimeAfter,nTimeTotal);
1412         printf("row0 = %f, col0 = %f, time0 = %f\n"
1413               ,row0,col0,time0);
1414         printf("rowPadSize = %f, colPadSize = %f, timeBinSize = %f\n"
1415                ,rowPadSize,colPadSize,timeBinSize);
1416       }
1417
1418       // Don't analyze test hits and switched off detectors
1419       if ((CheckDetector(plane,chamber,sector)) &&
1420           (((Int_t) q) != 0)) {
1421
1422         if (detector != detectorOld) {
1423
1424           if (fDebug > 1) {
1425             printf("<AliTRDdigitizer::MakeDigits> ");
1426             printf("Get new container. New det = %d, Old det = %d\n"
1427                   ,detector,detectorOld);
1428           }
1429           // Compress the old one if enabled
1430           if ((fCompress) && (detectorOld > -1)) {
1431             if (fDebug > 1) {
1432               printf("<AliTRDdigitizer::MakeDigits> ");
1433               printf("Compress the old container ...");
1434             }
1435             signals->Compress(1,0);
1436             for (iDict = 0; iDict < kNDict; iDict++) {
1437               dictionary[iDict]->Compress(1,0);
1438             }
1439             if (fDebug > 1) printf("done\n");
1440           }
1441           // Get the new container
1442           signals = (AliTRDdataArrayF *) signalsArray->At(detector);
1443           if (signals->GetNtime() == 0) {
1444             // Allocate a new one if not yet existing
1445             if (fDebug > 1) {
1446               printf("<AliTRDdigitizer::MakeDigits> ");
1447               printf("Allocate a new container ... ");
1448             }
1449             signals->Allocate(nRowMax,nColMax,nTimeTotal);
1450           }
1451           else {
1452             // Expand an existing one
1453             if (fCompress) {
1454               if (fDebug > 1) {
1455                 printf("<AliTRDdigitizer::MakeDigits> ");
1456                 printf("Expand an existing container ... ");
1457               }
1458               signals->Expand();
1459             }
1460           }
1461           // The same for the dictionary
1462           for (iDict = 0; iDict < kNDict; iDict++) {
1463             dictionary[iDict] = fDigitsManager->GetDictionary(detector,iDict);
1464             if (dictionary[iDict]->GetNtime() == 0) {
1465               dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1466             }
1467             else {
1468               if (fCompress) dictionary[iDict]->Expand();
1469             }
1470           }
1471           if (fDebug > 1) printf("done\n");
1472           detectorOld = detector;
1473         }
1474
1475         // Rotate the sectors on top of each other
1476         fGeo->Rotate(detector,pos,rot);
1477
1478         // The driftlength. It is negative if the hit is in the
1479         // amplification region.
1480         Float_t driftlength = time0 - rot[0];
1481
1482         // Take also the drift in the amplification region into account
1483         // The drift length is at the moment still the same, regardless of
1484         // the position relativ to the wire. This non-isochronity needs still
1485         // to be implemented.
1486         Float_t driftlengthL = TMath::Abs(driftlength + kAmWidth);
1487         if (fExBOn) driftlengthL /= TMath::Sqrt(fLorentzFactor);
1488
1489         // Loop over all electrons of this hit
1490         // TR photons produce hits with negative charge
1491         Int_t nEl = ((Int_t) TMath::Abs(q));
1492         for (Int_t iEl = 0; iEl < nEl; iEl++) {
1493
1494           xyz[0] = rot[0];
1495           xyz[1] = rot[1];
1496           xyz[2] = rot[2];
1497
1498           // Electron attachment
1499           if (fElAttachOn) {
1500             if (gRandom->Rndm() < (driftlengthL * elAttachProp))
1501               continue;
1502           }
1503
1504           // Apply the diffusion smearing
1505           if (fDiffusionOn) {
1506             if (!(Diffusion(driftlengthL,xyz))) continue;
1507           }
1508
1509           // Apply E x B effects (depends on drift direction)
1510           if (fExBOn) {
1511             if (!(ExB(driftlength+kAmWidth,xyz))) continue;
1512           }
1513
1514           // The electron position after diffusion and ExB in pad coordinates
1515           // The pad row (z-direction)
1516           Float_t rowDist   = xyz[2] - row0;
1517           Int_t   rowE      = ((Int_t) (rowDist * divideRow));
1518           if ((rowE < 0) || (rowE >= nRowMax)) continue;
1519           Float_t rowOffset = ((((Float_t) rowE) + 0.5) * rowPadSize) - rowDist;
1520
1521           // The pad column (rphi-direction)
1522           Float_t col0tilt  =  Col0Tilted(col0,rowOffset,plane);
1523           Float_t colDist   = xyz[1] - col0tilt;
1524           Int_t   colE      = ((Int_t) (colDist * divideCol));
1525           if ((colE < 0) || (colE >= nColMax)) continue;
1526           Float_t colOffset = ((((Float_t) colE) + 0.5) * colPadSize) - colDist;
1527
1528           // The time bin (negative for hits in the amplification region)
1529           // In the amplification region the electrons drift from both sides
1530           // to the middle (anode wire plane)
1531           Float_t timeDist   = time0 - xyz[0];
1532           Float_t timeOffset = 0;
1533           Int_t   timeE      = 0;
1534           if (timeDist > 0) {
1535             // The time bin
1536             timeE      = ((Int_t) (timeDist * divideTime));
1537             // The distance of the position to the middle of the timebin
1538             timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) - timeDist;
1539           }
1540           else {
1541             // Difference between half of the amplification gap width and
1542             // the distance to the anode wire
1543             Float_t anodeDist = kAmWidth - TMath::Abs(timeDist + kAmWidth);
1544             // The time bin
1545             timeE      = -1 * (((Int_t ) (anodeDist * divideTime)) + 1);
1546             // The distance of the position to the middle of the timebin
1547             timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) + anodeDist;
1548           }
1549
1550           // Apply the gas gain including fluctuations
1551           Float_t ggRndm = 0.0;
1552           do {
1553             ggRndm = gRandom->Rndm();
1554           } while (ggRndm <= 0);
1555           Int_t signal = (Int_t) (-fGasGain * TMath::Log(ggRndm));
1556
1557           // Apply the pad response
1558           if (fPRFOn) {
1559             // The distance of the electron to the center of the pad
1560             // in units of pad width
1561             Float_t dist = - colOffset * divideCol;
1562             if (!(PadResponse(signal,dist,plane,padSignal))) continue;
1563           }
1564           else {
1565             padSignal[0] = 0.0;
1566             padSignal[1] = signal;
1567             padSignal[2] = 0.0;
1568           }
1569
1570           // Sample the time response inside the drift region
1571           // + additional time bins before and after.
1572           // The sampling is done always in the middle of the time bin
1573           for (Int_t iTimeBin = TMath::Max(timeE-timeTRDbeg,        -nTimeBefore)
1574                     ;iTimeBin < TMath::Min(timeE+timeTRDend,nTimeMax+nTimeAfter )
1575                     ;iTimeBin++) {
1576
1577             // Apply the time response
1578             Float_t timeResponse = 1.0;
1579             Float_t crossTalk    = 0.0;
1580             Float_t time         = (iTimeBin - timeE) * timeBinSize + timeOffset;
1581             if (fTRFOn) {
1582               timeResponse = TimeResponse(time);
1583             }
1584             if (fCTOn) {
1585               crossTalk    = CrossTalk(time);
1586             }
1587
1588             signalOld[0] = 0.0;
1589             signalOld[1] = 0.0;
1590             signalOld[2] = 0.0;
1591
1592             for (iPad = 0; iPad < kNpad; iPad++) {
1593
1594               Int_t colPos = colE + iPad - 1;
1595               if (colPos <        0) continue;
1596               if (colPos >= nColMax) break;
1597
1598               // Add the signals
1599               // Note: The time bin number is shifted by nTimeBefore to avoid negative
1600               // time bins. This has to be subtracted later.
1601               Int_t iCurrentTimeBin = iTimeBin + nTimeBefore;
1602               signalOld[iPad]  = signals->GetDataUnchecked(rowE,colPos,iCurrentTimeBin);
1603               if( colPos != colE ) {
1604                 signalOld[iPad] += padSignal[iPad] * (timeResponse + crossTalk);
1605               }
1606               else {
1607                 signalOld[iPad] += padSignal[iPad] * timeResponse;
1608               }
1609               signals->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,signalOld[iPad]);
1610
1611               // Store the track index in the dictionary
1612               // Note: We store index+1 in order to allow the array to be compressed
1613               if (signalOld[iPad] > 0) {
1614                 for (iDict = 0; iDict < kNDict; iDict++) {
1615                   Int_t oldTrack = dictionary[iDict]->GetDataUnchecked(rowE
1616                                                                       ,colPos
1617                                                                       ,iCurrentTimeBin);
1618                   if (oldTrack == track+1) break;
1619                   if (oldTrack ==       0) {
1620                     dictionary[iDict]->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,track+1);
1621                     break;
1622                   }
1623                 }
1624               }
1625
1626             } // Loop: pads
1627
1628           } // Loop: time bins
1629
1630         } // Loop: electrons of a single hit
1631
1632       } // If: detector and test hit
1633
1634       hit = (AliTRDhit *) fTRD->NextHit();
1635
1636     } // Loop: hits of one primary track
1637
1638   } // Loop: primary tracks
1639
1640   if (fDebug > 0) {
1641     printf("<AliTRDdigitizer::MakeDigits> ");
1642     printf("Finished analyzing %d hits\n",countHits);
1643   }
1644
1645   // The total conversion factor
1646   Float_t convert = kEl2fC * fPadCoupling * fTimeCoupling * fChipGain;
1647
1648   // Loop through all chambers to finalize the digits
1649   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1650
1651     Int_t plane       = fGeo->GetPlane(iDet);
1652     Int_t sector      = fGeo->GetSector(iDet);
1653     Int_t chamber     = fGeo->GetChamber(iDet);
1654     Int_t nRowMax     = fGeo->GetRowMax(plane,chamber,sector);
1655     Int_t nColMax     = fGeo->GetColMax(plane);
1656     Int_t nTimeMax    = fGeo->GetTimeMax();
1657     Int_t nTimeTotal  = fGeo->GetTimeTotal();
1658
1659     Double_t *inADC  = new Double_t[nTimeTotal];
1660     Double_t *outADC = new Double_t[nTimeTotal];
1661
1662     if (fDebug > 0) {
1663       printf("<AliTRDdigitizer::MakeDigits> ");
1664       printf("Digitization for chamber %d\n",iDet);
1665     }
1666
1667     // Add a container for the digits of this detector
1668     digits = fDigitsManager->GetDigits(iDet);
1669     // Allocate memory space for the digits buffer
1670     digits->Allocate(nRowMax,nColMax,nTimeTotal);
1671
1672     // Get the signal container
1673     signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
1674     if (signals->GetNtime() == 0) {
1675       // Create missing containers
1676       signals->Allocate(nRowMax,nColMax,nTimeTotal);
1677     }
1678     else {
1679       // Expand the container if neccessary
1680       if (fCompress) signals->Expand();
1681     }
1682     // Create the missing dictionary containers
1683     for (iDict = 0; iDict < kNDict; iDict++) {
1684       dictionary[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1685       if (dictionary[iDict]->GetNtime() == 0) {
1686         dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1687       }
1688     }
1689
1690     Int_t nDigits = 0;
1691
1692     // Don't create noise in detectors that are switched off
1693     if (CheckDetector(plane,chamber,sector)) {
1694
1695       // Create the digits for this chamber
1696       for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1697         for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1698
1699           // Create summable digits
1700           if (fSDigits) {
1701
1702             for (iTime = 0; iTime < nTimeTotal; iTime++) {
1703               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1704               signalAmp *= fSDigitsScale;
1705               signalAmp  = TMath::Min(signalAmp,(Float_t) 1.0e9);
1706               Int_t adc  = (Int_t) signalAmp;
1707               nDigits++;
1708               digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1709             }
1710
1711           }
1712           // Create normal digits
1713           else {
1714
1715             for (iTime = 0; iTime < nTimeTotal; iTime++) {
1716               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1717               // Add the noise
1718               signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,fNoise),0.0);
1719               // Convert to mV
1720               signalAmp *= convert;
1721               // Convert to ADC counts. Set the overflow-bit fADCoutRange if the
1722               // signal is larger than fADCinRange
1723               Int_t adc  = 0;
1724               if (signalAmp >= fADCinRange) {
1725                 adc = ((Int_t) fADCoutRange);
1726               }
1727               else {
1728                 adc = ((Int_t) (signalAmp * (fADCoutRange / fADCinRange)));
1729               }
1730               inADC[iTime]  = adc;
1731               outADC[iTime] = adc;
1732             }
1733
1734             // Apply the tail cancelation via the digital filter
1735             if (fTCOn) {
1736               DeConvExp(inADC,outADC,nTimeTotal,fTCnexp);
1737             }
1738
1739             for (iTime = 0; iTime < nTimeTotal; iTime++) {
1740               // Store the amplitude of the digit if above threshold
1741               if (outADC[iTime] > fADCthreshold) {
1742                 if (fDebug > 2) {
1743                   printf("  iRow = %d, iCol = %d, iTime = %d, adc = %f\n"
1744                         ,iRow,iCol,iTime,outADC[iTime]);
1745                 }
1746                 nDigits++;
1747                 digits->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1748               }
1749             }
1750
1751           }
1752
1753         }
1754       }
1755
1756     }
1757
1758     // Compress the arrays
1759     digits->Compress(1,0);
1760     for (iDict = 0; iDict < kNDict; iDict++) {
1761       dictionary[iDict]->Compress(1,0);
1762     }
1763
1764     totalSizeDigits += digits->GetSize();
1765     totalSizeDict0  += dictionary[0]->GetSize();
1766     totalSizeDict1  += dictionary[1]->GetSize();
1767     totalSizeDict2  += dictionary[2]->GetSize();
1768
1769     Float_t nPixel = nRowMax * nColMax * nTimeMax;
1770     if (fDebug > 0) {
1771       printf("<AliTRDdigitizer::MakeDigits> ");
1772       printf("Found %d digits in detector %d (%3.0f).\n"
1773             ,nDigits,iDet
1774             ,100.0 * ((Float_t) nDigits) / nPixel);
1775     }
1776
1777     if (fCompress) signals->Compress(1,0);
1778
1779     delete [] inADC;
1780     delete [] outADC;
1781
1782   }
1783
1784   if (fDebug > 0) {
1785     printf("<AliTRDdigitizer::MakeDigits> ");
1786     printf("Total number of analyzed hits = %d\n",countHits);
1787     printf("<AliTRDdigitizer::MakeDigits> ");
1788     printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
1789                                                       ,totalSizeDict0
1790                                                       ,totalSizeDict1
1791                                                       ,totalSizeDict2);
1792   }
1793
1794   return kTRUE;
1795
1796 }
1797
1798 //_____________________________________________________________________________
1799 void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
1800 {
1801   //
1802   // Add a digits manager for s-digits to the input list.
1803   //
1804
1805   fSDigitsManagerList->Add(man);
1806
1807 }
1808
1809 //_____________________________________________________________________________
1810 void AliTRDdigitizer::DeleteSDigitsManager()
1811 {
1812   //
1813   // Removes digits manager from the input list.
1814   //
1815
1816   fSDigitsManagerList->Delete();
1817
1818 }
1819
1820 //_____________________________________________________________________________
1821 Bool_t AliTRDdigitizer::ConvertSDigits()
1822 {
1823   //
1824   // Converts s-digits to normal digits
1825   //
1826
1827   // Number of track dictionary arrays
1828   const Int_t    kNDict = AliTRDdigitsManager::kNDict;
1829
1830   // Converts number of electrons to fC
1831   const Double_t kEl2fC = 1.602E-19 * 1.0E15;
1832
1833   Int_t iDict = 0;
1834   Int_t iRow;
1835   Int_t iCol;
1836   Int_t iTime;
1837
1838   if (fDebug > 0) {
1839     this->Dump();
1840   }
1841
1842   Double_t sDigitsScale = 1.0 / GetSDigitsScale();
1843   Double_t noise        = GetNoise();
1844   Double_t padCoupling  = GetPadCoupling();
1845   Double_t timeCoupling = GetTimeCoupling();
1846   Double_t chipGain     = GetChipGain();
1847   Double_t convert      = kEl2fC * padCoupling * timeCoupling * chipGain;;
1848   Double_t adcInRange   = GetADCinRange();
1849   Double_t adcOutRange  = GetADCoutRange();
1850   Int_t    adcThreshold = GetADCthreshold();
1851
1852   AliTRDdataArrayI *digitsIn;
1853   AliTRDdataArrayI *digitsOut;
1854   AliTRDdataArrayI *dictionaryIn[kNDict];
1855   AliTRDdataArrayI *dictionaryOut[kNDict];
1856
1857   // Loop through the detectors
1858   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1859
1860     if (fDebug > 0) {
1861       printf("<AliTRDdigitizer::ConvertSDigits> ");
1862       printf("Convert detector %d to digits.\n",iDet);
1863     }
1864
1865     Int_t plane      = fGeo->GetPlane(iDet);
1866     Int_t sector     = fGeo->GetSector(iDet);
1867     Int_t chamber    = fGeo->GetChamber(iDet);
1868     Int_t nRowMax    = fGeo->GetRowMax(plane,chamber,sector);
1869     Int_t nColMax    = fGeo->GetColMax(plane);
1870     Int_t nTimeTotal = fGeo->GetTimeTotal();
1871
1872     Double_t *inADC  = new Double_t[nTimeTotal];
1873     Double_t *outADC = new Double_t[nTimeTotal];
1874
1875     digitsIn  = fSDigitsManager->GetDigits(iDet);
1876     digitsIn->Expand();
1877     digitsOut = fDigitsManager->GetDigits(iDet);
1878     digitsOut->Allocate(nRowMax,nColMax,nTimeTotal);
1879     for (iDict = 0; iDict < kNDict; iDict++) {
1880       dictionaryIn[iDict]  = fSDigitsManager->GetDictionary(iDet,iDict);
1881       dictionaryIn[iDict]->Expand();
1882       dictionaryOut[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1883       dictionaryOut[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1884     }
1885
1886     for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1887       for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1888
1889         for (iTime = 0; iTime < nTimeTotal; iTime++) {
1890           Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1891           signal *= sDigitsScale;
1892           // Add the noise
1893           signal  = TMath::Max((Double_t) gRandom->Gaus(signal,noise),0.0);
1894           // Convert to mV
1895           signal *= convert;
1896           // Convert to ADC counts. Set the overflow-bit adcOutRange if the
1897           // signal is larger than adcInRange
1898           Int_t adc  = 0;
1899           if (signal >= adcInRange) {
1900             adc = ((Int_t) adcOutRange);
1901           }
1902           else {
1903             adc = ((Int_t) (signal * (adcOutRange / adcInRange)));
1904           }
1905           inADC[iTime]  = adc;
1906           outADC[iTime] = adc;
1907         }
1908
1909         // Apply the tail cancelation via the digital filter
1910         if (fTCOn) {
1911           DeConvExp(inADC,outADC,nTimeTotal,fTCnexp);
1912         }
1913
1914         for (iTime = 0; iTime < nTimeTotal; iTime++) {
1915           // Store the amplitude of the digit if above threshold
1916           if (outADC[iTime] > adcThreshold) {
1917             digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1918             // Copy the dictionary
1919             for (iDict = 0; iDict < kNDict; iDict++) {
1920               Int_t track = dictionaryIn[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1921               dictionaryOut[iDict]->SetDataUnchecked(iRow,iCol,iTime,track);
1922             }
1923           }
1924         }
1925
1926       }
1927     }
1928
1929     if (fCompress) {
1930       digitsIn->Compress(1,0);
1931       digitsOut->Compress(1,0);
1932       for (iDict = 0; iDict < kNDict; iDict++) {
1933         dictionaryIn[iDict]->Compress(1,0);
1934         dictionaryOut[iDict]->Compress(1,0);
1935       }
1936     }
1937
1938     delete [] inADC;
1939     delete [] outADC;
1940
1941   }
1942
1943   return kTRUE;
1944
1945 }
1946
1947 //_____________________________________________________________________________
1948 Bool_t AliTRDdigitizer::MergeSDigits()
1949 {
1950   //
1951   // Merges the input s-digits:
1952   //   - The amplitude of the different inputs are summed up.
1953   //   - Of the track IDs from the input dictionaries only one is
1954   //     kept for each input. This works for maximal 3 different merged inputs.
1955   //
1956
1957   // Number of track dictionary arrays
1958   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1959
1960   Int_t iDict = 0;
1961   Int_t jDict = 0;
1962
1963   AliTRDdataArrayI *digitsA;
1964   AliTRDdataArrayI *digitsB;
1965   AliTRDdataArrayI *dictionaryA[kNDict];
1966   AliTRDdataArrayI *dictionaryB[kNDict];
1967
1968   // Get the first s-digits
1969   fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
1970   if (!fSDigitsManager) return kFALSE;
1971
1972   // Loop through the other sets of s-digits
1973   AliTRDdigitsManager *mergeSDigitsManager;
1974   mergeSDigitsManager = (AliTRDdigitsManager *)
1975                         fSDigitsManagerList->After(fSDigitsManager);
1976
1977   if (fDebug > 0) {
1978     if (mergeSDigitsManager) {
1979       printf("<AliTRDdigitizer::MergeSDigits> ");
1980       printf("Merge %d input files.\n",fSDigitsManagerList->GetSize());
1981     }
1982     else {
1983       printf("<AliTRDdigitizer::MergeSDigits> ");
1984       printf("Only one input file.\n");
1985     }
1986   }
1987
1988   Int_t iMerge = 0;
1989   while (mergeSDigitsManager) {
1990
1991     iMerge++;
1992
1993     // Loop through the detectors
1994     for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1995
1996       Int_t plane      = fGeo->GetPlane(iDet);
1997       Int_t sector     = fGeo->GetSector(iDet);
1998       Int_t chamber    = fGeo->GetChamber(iDet);
1999       Int_t nRowMax    = fGeo->GetRowMax(plane,chamber,sector);
2000       Int_t nColMax    = fGeo->GetColMax(plane);
2001       Int_t nTimeTotal = fGeo->GetTimeTotal();
2002
2003       // Loop through the pixels of one detector and add the signals
2004       digitsA = fSDigitsManager->GetDigits(iDet);
2005       digitsB = mergeSDigitsManager->GetDigits(iDet);
2006       digitsA->Expand();
2007       digitsB->Expand();
2008       for (iDict = 0; iDict < kNDict; iDict++) {
2009         dictionaryA[iDict] = fSDigitsManager->GetDictionary(iDet,iDict);
2010         dictionaryB[iDict] = mergeSDigitsManager->GetDictionary(iDet,iDict);
2011         dictionaryA[iDict]->Expand();
2012         dictionaryB[iDict]->Expand();
2013       }
2014
2015       if (fDebug > 0) {
2016         printf("<AliTRDdigitizer::MergeSDigits> ");
2017         printf("Merge detector %d of input no.%d\n",iDet,iMerge+1);
2018       }
2019
2020       for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
2021         for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
2022           for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
2023
2024             // Add the amplitudes of the summable digits
2025             Int_t ampA = digitsA->GetDataUnchecked(iRow,iCol,iTime);
2026             Int_t ampB = digitsB->GetDataUnchecked(iRow,iCol,iTime);
2027             ampA += ampB;
2028             digitsA->SetDataUnchecked(iRow,iCol,iTime,ampA);
2029
2030             // Add the mask to the track id if defined.
2031             for (iDict = 0; iDict < kNDict; iDict++) {
2032               Int_t trackB = dictionaryB[iDict]->GetDataUnchecked(iRow,iCol,iTime);
2033               if ((fMasks) && (trackB > 0)) {
2034                 for (jDict = 0; jDict < kNDict; jDict++) {
2035                   Int_t trackA = dictionaryA[iDict]->GetDataUnchecked(iRow,iCol,iTime);
2036                   if (trackA == 0) {
2037                     trackA = trackB + fMasks[iMerge];
2038                     dictionaryA[iDict]->SetDataUnchecked(iRow,iCol,iTime,trackA);
2039                   }
2040                 }
2041               }
2042             }
2043
2044           }
2045         }
2046       }
2047
2048       if (fCompress) {
2049         digitsA->Compress(1,0);
2050         digitsB->Compress(1,0);
2051         for (iDict = 0; iDict < kNDict; iDict++) {
2052           dictionaryA[iDict]->Compress(1,0);
2053           dictionaryB[iDict]->Compress(1,0);
2054         }
2055       }
2056
2057     }
2058
2059     // The next set of s-digits
2060     mergeSDigitsManager = (AliTRDdigitsManager *)
2061                           fSDigitsManagerList->After(mergeSDigitsManager);
2062
2063   }
2064
2065   return kTRUE;
2066
2067 }
2068
2069 //_____________________________________________________________________________
2070 Bool_t AliTRDdigitizer::SDigits2Digits()
2071 {
2072   //
2073   // Merges the input s-digits and converts them to normal digits
2074   //
2075
2076   if (!MergeSDigits()) return kFALSE;
2077
2078   return ConvertSDigits();
2079
2080 }
2081
2082 //_____________________________________________________________________________
2083 Bool_t AliTRDdigitizer::CheckDetector(Int_t plane, Int_t chamber, Int_t sector)
2084 {
2085   //
2086   // Checks whether a detector is enabled
2087   //
2088
2089   if ((fTRD->GetSensChamber() >=       0) &&
2090       (fTRD->GetSensChamber() != chamber)) return kFALSE;
2091   if ((fTRD->GetSensPlane()   >=       0) &&
2092       (fTRD->GetSensPlane()   !=   plane)) return kFALSE;
2093   if ( fTRD->GetSensSector()  >=       0) {
2094     Int_t sens1 = fTRD->GetSensSector();
2095     Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
2096     sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
2097            * AliTRDgeometry::Nsect();
2098     if (sens1 < sens2) {
2099       if ((sector < sens1) || (sector >= sens2)) return kFALSE;
2100     }
2101     else {
2102       if ((sector < sens1) && (sector >= sens2)) return kFALSE;
2103     }
2104   }
2105
2106   return kTRUE;
2107
2108 }
2109
2110 //_____________________________________________________________________________
2111 Bool_t AliTRDdigitizer::WriteDigits()
2112 {
2113   //
2114   // Writes out the TRD-digits and the dictionaries
2115   //
2116
2117   // Store the digits and the dictionary in the tree
2118   return fDigitsManager->WriteDigits();
2119
2120 }
2121
2122 //_____________________________________________________________________________
2123 void AliTRDdigitizer::SetTiltingAngle(Float_t v)
2124 {
2125   //
2126   // Set the tilting angle for the readout pads
2127   //
2128
2129   fTiltingAngle = TMath::Tan(TMath::Pi()/180.0 * v);
2130
2131 }
2132
2133 //_____________________________________________________________________________
2134 Float_t AliTRDdigitizer::GetTiltingAngle() const
2135 {
2136   //
2137   // Get the tilting angle for the readout pads
2138   //
2139
2140   return TMath::ATan(180.0/TMath::Pi() * fTiltingAngle);
2141
2142 }
2143
2144 //_____________________________________________________________________________
2145 Float_t AliTRDdigitizer::GetDiffusionL(Float_t vd, Float_t b)
2146 {
2147   //
2148   // Returns the longitudinal diffusion coefficient for a given drift
2149   // velocity <vd> and a B-field <b> for Xe/CO2 (15%).
2150   // The values are according to a GARFIELD simulation.
2151   //
2152
2153   const Int_t kNb = 5;
2154   Float_t p0[kNb] = {  0.007440,  0.007493,  0.007513,  0.007672,  0.007831 };
2155   Float_t p1[kNb] = {  0.019252,  0.018912,  0.018636,  0.018012,  0.017343 };
2156   Float_t p2[kNb] = { -0.005042, -0.004926, -0.004867, -0.004650, -0.004424 };
2157   Float_t p3[kNb] = {  0.000195,  0.000189,  0.000195,  0.000182,  0.000169 };
2158
2159   Int_t ib = ((Int_t) (10 * (b - 0.15)));
2160   ib       = TMath::Max(  0,ib);
2161   ib       = TMath::Min(kNb,ib);
2162
2163   Float_t diff = p0[ib]
2164                + p1[ib] * vd
2165                + p2[ib] * vd*vd
2166                + p3[ib] * vd*vd*vd;
2167
2168   return diff;
2169
2170 }
2171
2172 //_____________________________________________________________________________
2173 Float_t AliTRDdigitizer::GetDiffusionT(Float_t vd, Float_t b)
2174 {
2175   //
2176   // Returns the transverse diffusion coefficient for a given drift
2177   // velocity <vd> and a B-field <b> for Xe/CO2 (15%).
2178   // The values are according to a GARFIELD simulation.
2179   //
2180
2181   const Int_t kNb = 5;
2182   Float_t p0[kNb] = {  0.009550,  0.009599,  0.009674,  0.009757,  0.009850 };
2183   Float_t p1[kNb] = {  0.006667,  0.006539,  0.006359,  0.006153,  0.005925 };
2184   Float_t p2[kNb] = { -0.000853, -0.000798, -0.000721, -0.000635, -0.000541 };
2185   Float_t p3[kNb] = {  0.000131,  0.000122,  0.000111,  0.000098,  0.000085 };
2186
2187   Int_t ib = ((Int_t) (10 * (b - 0.15)));
2188   ib       = TMath::Max(  0,ib);
2189   ib       = TMath::Min(kNb,ib);
2190
2191   Float_t diff = p0[ib]
2192                + p1[ib] * vd
2193                + p2[ib] * vd*vd
2194                + p3[ib] * vd*vd*vd;
2195
2196   return diff;
2197
2198 }
2199
2200 //_____________________________________________________________________________
2201 Float_t AliTRDdigitizer::GetOmegaTau(Float_t vd, Float_t b)
2202 {
2203   //
2204   // Returns omega*tau (tan(Lorentz-angle)) for a given drift velocity <vd>
2205   // and a B-field <b> for Xe/CO2 (15%).
2206   // The values are according to a GARFIELD simulation.
2207   //
2208
2209   const Int_t kNb = 5;
2210   Float_t p0[kNb] = {  0.004810,  0.007412,  0.010252,  0.013409,  0.016888 };
2211   Float_t p1[kNb] = {  0.054875,  0.081534,  0.107333,  0.131983,  0.155455 };
2212   Float_t p2[kNb] = { -0.008682, -0.012896, -0.016987, -0.020880, -0.024623 };
2213   Float_t p3[kNb] = {  0.000155,  0.000238,  0.000330,  0.000428,  0.000541 };
2214
2215   Int_t ib = ((Int_t) (10 * (b - 0.15)));
2216   ib       = TMath::Max(  0,ib);
2217   ib       = TMath::Min(kNb,ib);
2218
2219   Float_t alphaL = p0[ib]
2220                  + p1[ib] * vd
2221                  + p2[ib] * vd*vd
2222                  + p3[ib] * vd*vd*vd;
2223
2224   return TMath::Tan(alphaL);
2225
2226 }
2227
2228 //_____________________________________________________________________________
2229 void AliTRDdigitizer::DeConvExp(Double_t *source, Double_t *target
2230                               , Int_t n, Int_t nexp)
2231 {
2232   //
2233   // Does the deconvolution by the digital filter.
2234   //
2235   // Author:        Marcus Gutfleisch, KIP Heidelberg
2236   // Optimized for: New TRF from Venelin Angelov, simulated with CADENCE
2237   //                Pad-ground capacitance = 25 pF
2238   //                Pad-pad cross talk capacitance = 6 pF
2239   //                For 10 MHz digitization, corresponding to 20 time bins
2240   //                in the drift region
2241   //
2242
2243   Double_t rates[2];
2244   Double_t coefficients[2];
2245
2246   /* initialize (coefficient = alpha, rates = lambda) */
2247
2248   if( nexp == 1 ) {
2249     rates[0] = 0.466998;
2250     /* no rescaling */
2251     coefficients[0] = 1.0;
2252   }
2253   if( nexp == 2 ) {
2254     rates[0] = 0.8988162;
2255     coefficients[0] = 0.11392069;
2256     rates[1] = 0.3745688;
2257     coefficients[1] = 0.8860793;
2258     /* no rescaling */
2259     Float_t sumc = coefficients[0]+coefficients[1];
2260     coefficients[0] /= sumc;
2261     coefficients[1] /= sumc;
2262   }
2263
2264   Int_t i, k;
2265   Double_t reminder[2];
2266   Double_t correction, result;
2267
2268   /* attention: computation order is important */
2269   correction=0.0;
2270   for ( k = 0; k < nexp; k++ ) reminder[k]=0.0;
2271
2272   for ( i = 0; i < n; i++ ) {
2273     result = ( source[i] - correction );    /* no rescaling */
2274     target[i] = result;
2275
2276     for ( k = 0; k < nexp; k++ ) reminder[k] = rates[k]
2277                              * ( reminder[k] + coefficients[k] * result);
2278
2279     correction=0.0;
2280     for ( k = 0; k < nexp; k++ ) correction += reminder[k];
2281   }
2282
2283 }
2284