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.27  2001/11/14 10:50:45  cblume
  19 Changes in digits IO. Add merging of summable digits
  20
  21 Revision 1.26  2001/11/06 17:19:41  cblume
  22 Add detailed geometry and simple simulator
  23
  24 Revision 1.25  2001/06/27 09:54:44  cblume
  25 Moved fField initialization to InitDetector()
  26
  27 Revision 1.24  2001/05/21 16:45:47  hristov
  28 Last minute changes (C.Blume)
  29
  30 Revision 1.23  2001/05/07 08:04:48  cblume
  31 New TRF and PRF. Speedup of the code. Digits from amplification region included
  32
  33 Revision 1.22  2001/03/30 14:40:14  cblume
  34 Update of the digitization parameter
  35
  36 Revision 1.21  2001/03/13 09:30:35  cblume
  37 Update of digitization. Moved digit branch definition to AliTRD
  38
  39 Revision 1.20  2001/02/25 20:19:00  hristov
  40 Minor correction: loop variable declared only once for HP, Sun
  41
  42 Revision 1.19  2001/02/14 18:22:26  cblume
  43 Change in the geometry of the padplane
  44
  45 Revision 1.18  2001/01/26 19:56:57  hristov
  46 Major upgrade of AliRoot code
  47
  48 Revision 1.17  2000/12/08 12:53:27  cblume
  49 Change in Copy() function for HP-compiler
  50
  51 Revision 1.16  2000/12/07 12:20:46  cblume
  52 Go back to array compression. Use sampled PRF to speed up digitization
  53
  54 Revision 1.15  2000/11/23 14:34:08  cblume
  55 Fixed bug in expansion routine of arrays (initialize buffers properly)
  56
  57 Revision 1.14  2000/11/20 08:54:44  cblume
  58 Switch off compression as default
  59
  60 Revision 1.13  2000/11/10 14:57:52  cblume
  61 Changes in the geometry constants for the DEC compiler
  62
  63 Revision 1.12  2000/11/01 14:53:20  cblume
  64 Merge with TRD-develop
  65
  66 Revision 1.1.4.9  2000/10/26 17:00:22  cblume
  67 Fixed bug in CheckDetector()
  68
  69 Revision 1.1.4.8  2000/10/23 13:41:35  cblume
  70 Added protection against Log(0) in the gas gain calulation
  71
  72 Revision 1.1.4.7  2000/10/17 02:27:34  cblume
  73 Get rid of global constants
  74
  75 Revision 1.1.4.6  2000/10/16 01:16:53  cblume
  76 Changed timebin 0 to be the one closest to the readout
  77
  78 Revision 1.1.4.5  2000/10/15 23:34:29  cblume
  79 Faster version of the digitizer
  80
  81 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
  82 Made Getters const
  83
  84 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
  85 Replace include files by forward declarations
  86
  87 Revision 1.1.4.2  2000/09/22 14:41:10  cblume
  88 Bug fix in PRF. Included time response. New structure
  89
  90 Revision 1.10  2000/10/05 07:27:53  cblume
  91 Changes in the header-files by FCA
  92
  93 Revision 1.9  2000/10/02 21:28:19  fca
  94 Removal of useless dependecies via forward declarations
  95
  96 Revision 1.8  2000/06/09 11:10:07  cblume
  97 Compiler warnings and coding conventions, next round
  98
  99 Revision 1.7  2000/06/08 18:32:58  cblume
 100 Make code compliant to coding conventions
 101
 102 Revision 1.6  2000/06/07 16:27:32  cblume
 103 Try to remove compiler warnings on Sun and HP
 104
 105 Revision 1.5  2000/05/09 16:38:57  cblume
 106 Removed PadResponse(). Merge problem
 107
 108 Revision 1.4  2000/05/08 15:53:45  cblume
 109 Resolved merge conflict
 110
 111 Revision 1.3  2000/04/28 14:49:27  cblume
 112 Only one declaration of iDict in MakeDigits()
 113
 114 Revision 1.1.4.1  2000/05/08 14:42:04  cblume
 115 Introduced AliTRDdigitsManager
 116
 117 Revision 1.1  2000/02/28 19:00:13  cblume
 118 Add new TRD classes
 119
 120 */
 121
 122 ///////////////////////////////////////////////////////////////////////////////
 123 //                                                                           //
 124 //  Creates and handles digits from TRD hits                                 //
 125 //  Author: C. Blume (C.Blume@gsi.de)                                        //
 126 //                                                                           //
 127 //  The following effects are included:                                      //
 128 //      - Diffusion                                                          //
 129 //      - ExB effects                                                        //
 130 //      - Gas gain including fluctuations                                    //
 131 //      - Pad-response (simple Gaussian approximation)                       //
 132 //      - Time-response                                                      //
 133 //      - Electronics noise                                                  //
 134 //      - Electronics gain                                                   //
 135 //      - Digitization                                                       //
 136 //      - ADC threshold                                                      //
 137 //  The corresponding parameter can be adjusted via the various              //
 138 //  Set-functions. If these parameters are not explicitly set, default       //
 139 //  values are used (see Init-function).                                     //
 140 //  As an example on how to use this class to produce digits from hits       //
 141 //  have a look at the macro hits2digits.C                                   //
 142 //  The production of summable digits is demonstrated in hits2sdigits.C      //
 143 //  and the subsequent conversion of the s-digits into normal digits is      //
 144 //  explained in sdigits2digits.C.                                           //
 145 //                                                                           //
 146 ///////////////////////////////////////////////////////////////////////////////
 147
 148 #include <stdlib.h>
 149
 150 #include <TMath.h>
 151 #include <TVector.h>
 152 #include <TRandom.h>
 153 #include <TROOT.h>
 154 #include <TTree.h>
 155 #include <TFile.h>
 156 #include <TF1.h>
 157 #include <TList.h>
 158 #include <TTask.h>
 159
 160 #include "AliRun.h"
 161 #include "AliMagF.h"
 162 #include "AliRunDigitizer.h"
 163
 164 #include "AliTRD.h"
 165 #include "AliTRDhit.h"
 166 #include "AliTRDdigitizer.h"
 167 #include "AliTRDdataArrayI.h"
 168 #include "AliTRDdataArrayF.h"
 169 #include "AliTRDsegmentArray.h"
 170 #include "AliTRDdigitsManager.h"
 171 #include "AliTRDgeometry.h"
 172
 173 ClassImp(AliTRDdigitizer)
 174
 175 //_____________________________________________________________________________
 176 AliTRDdigitizer::AliTRDdigitizer()
 177 {
 178   //
 179   // AliTRDdigitizer default constructor
 180   //
 181
 182   fInputFile          = NULL;
 183   fDigitsManager      = NULL;
 184   fSDigitsManagerList = NULL;
 185   fSDigitsManager     = NULL;
 186   fTRD                = NULL;
 187   fGeo                = NULL;
 188   fPRFsmp             = NULL;
 189   fTRFsmp             = NULL;
 190
 191   fEvent              = 0;
 192   fGasGain            = 0.0;
 193   fNoise              = 0.0;
 194   fChipGain           = 0.0;
 195   fADCoutRange        = 0.0;
 196   fADCinRange         = 0.0;
 197   fADCthreshold       = 0;
 198   fDiffusionOn        = 0;
 199   fDiffusionT         = 0.0;
 200   fDiffusionL         = 0.0;
 201   fElAttachOn         = 0;
 202   fElAttachProp       = 0.0;
 203   fExBOn              = 0;
 204   fOmegaTau           = 0.0;
 205   fPRFOn              = 0;
 206   fTRFOn              = 0;
 207   fDriftVelocity      = 0.0;
 208   fPadCoupling        = 0.0;
 209   fTimeCoupling       = 0.0;
 210   fTimeBinWidth       = 0.0;
 211   fField              = 0.0;
 212
 213   fPRFbin             = 0;
 214   fPRFlo              = 0.0;
 215   fPRFhi              = 0.0;
 216   fPRFwid             = 0.0;
 217   fPRFpad             = 0;
 218   fTRFbin             = 0;
 219   fTRFlo              = 0.0;
 220   fTRFhi              = 0.0;
 221   fTRFwid             = 0.0;
 222
 223   fCompress           = kTRUE;
 224   fVerbose            = 0;
 225   fSDigits            = kFALSE;
 226   fSDigitsScale       = 0.0;
 227
 228 }
 229
 230 //_____________________________________________________________________________
 231 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
 232                 :AliDigitizer(name,title)
 233 {
 234   //
 235   // AliTRDdigitizer constructor
 236   //
 237
 238   fInputFile          = NULL;
 239
 240   fDigitsManager      = NULL;
 241   fSDigitsManager     = NULL;
 242   fSDigitsManagerList = NULL;
 243
 244   fTRD                = NULL;
 245   fGeo                = NULL;
 246   fPRFsmp             = NULL;
 247   fTRFsmp             = NULL;
 248
 249   fEvent              = 0;
 250
 251   fCompress           = kTRUE;
 252   fVerbose            = 0;
 253   fSDigits            = kFALSE;
 254
 255   Init();
 256
 257 }
 258
 259 //_____________________________________________________________________________
 260 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager
 261                                 , const Text_t *name, const Text_t *title)
 262                 :AliDigitizer(manager,name,title)
 263 {
 264   //
 265   // AliTRDdigitizer constructor
 266   //
 267
 268   fInputFile          = NULL;
 269
 270   fDigitsManager      = NULL;
 271   fSDigitsManager     = NULL;
 272   fSDigitsManagerList = NULL;
 273
 274   fTRD                = NULL;
 275   fGeo                = NULL;
 276   fPRFsmp             = NULL;
 277   fTRFsmp             = NULL;
 278
 279   fEvent              = 0;
 280
 281   fCompress           = kTRUE;
 282   fVerbose            = 0;
 283   fSDigits            = kFALSE;
 284
 285   Init();
 286
 287 }
 288
 289 //_____________________________________________________________________________
 290 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
 291 {
 292   //
 293   // AliTRDdigitizer copy constructor
 294   //
 295
 296   ((AliTRDdigitizer &) d).Copy(*this);
 297
 298 }
 299
 300 //_____________________________________________________________________________
 301 AliTRDdigitizer::~AliTRDdigitizer()
 302 {
 303   //
 304   // AliTRDdigitizer destructor
 305   //
 306
 307   if (fInputFile) {
 308     fInputFile->Close();
 309     delete fInputFile;
 310     fInputFile = NULL;
 311   }
 312
 313   if (fDigitsManager) {
 314     delete fDigitsManager;
 315     fDigitsManager = NULL;
 316   }
 317
 318   if (fSDigitsManager) {
 319     delete fSDigitsManager;
 320     fSDigitsManager = NULL;
 321   }
 322
 323   if (fSDigitsManagerList) {
 324     fSDigitsManagerList->Delete();
 325     delete fSDigitsManagerList;
 326     fSDigitsManagerList = NULL;
 327   }
 328 }
 329
 330 //_____________________________________________________________________________
 331 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
 332 {
 333   //
 334   // Assignment operator
 335   //
 336
 337   if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
 338   return *this;
 339
 340 }
 341
 342 //_____________________________________________________________________________
 343 void AliTRDdigitizer::Copy(TObject &d)
 344 {
 345   //
 346   // Copy function
 347   //
 348
 349   Int_t iBin;
 350
 351   ((AliTRDdigitizer &) d).fInputFile          = NULL;
 352   ((AliTRDdigitizer &) d).fSDigitsManagerList = NULL;
 353   ((AliTRDdigitizer &) d).fSDigitsManager     = NULL;
 354   ((AliTRDdigitizer &) d).fDigitsManager      = NULL;
 355   ((AliTRDdigitizer &) d).fTRD                = NULL;
 356   ((AliTRDdigitizer &) d).fGeo                = NULL;
 357
 358   ((AliTRDdigitizer &) d).fEvent              = 0;
 359
 360   ((AliTRDdigitizer &) d).fGasGain            = fGasGain;
 361   ((AliTRDdigitizer &) d).fNoise              = fNoise;
 362   ((AliTRDdigitizer &) d).fChipGain           = fChipGain;
 363   ((AliTRDdigitizer &) d).fADCoutRange        = fADCoutRange;
 364   ((AliTRDdigitizer &) d).fADCinRange         = fADCinRange;
 365   ((AliTRDdigitizer &) d).fADCthreshold       = fADCthreshold;
 366   ((AliTRDdigitizer &) d).fDiffusionOn        = fDiffusionOn;
 367   ((AliTRDdigitizer &) d).fDiffusionT         = fDiffusionT;
 368   ((AliTRDdigitizer &) d).fDiffusionL         = fDiffusionL;
 369   ((AliTRDdigitizer &) d).fElAttachOn         = fElAttachOn;
 370   ((AliTRDdigitizer &) d).fElAttachProp       = fElAttachProp;
 371   ((AliTRDdigitizer &) d).fExBOn              = fExBOn;
 372   ((AliTRDdigitizer &) d).fOmegaTau           = fOmegaTau;
 373   ((AliTRDdigitizer &) d).fLorentzFactor      = fLorentzFactor;
 374   ((AliTRDdigitizer &) d).fDriftVelocity      = fDriftVelocity;
 375   ((AliTRDdigitizer &) d).fPadCoupling        = fPadCoupling;
 376   ((AliTRDdigitizer &) d).fTimeCoupling       = fTimeCoupling;
 377   ((AliTRDdigitizer &) d).fTimeBinWidth       = fTimeBinWidth;
 378   ((AliTRDdigitizer &) d).fField              = fField;
 379   ((AliTRDdigitizer &) d).fPRFOn              = fPRFOn;
 380   ((AliTRDdigitizer &) d).fTRFOn              = fTRFOn;
 381
 382   ((AliTRDdigitizer &) d).fCompress           = fCompress;
 383   ((AliTRDdigitizer &) d).fVerbose            = fVerbose;
 384   ((AliTRDdigitizer &) d).fSDigits            = fSDigits;
 385   ((AliTRDdigitizer &) d).fSDigitsScale       = fSDigitsScale;
 386
 387   ((AliTRDdigitizer &) d).fPRFbin             = fPRFbin;
 388   ((AliTRDdigitizer &) d).fPRFlo              = fPRFlo;
 389   ((AliTRDdigitizer &) d).fPRFhi              = fPRFhi;
 390   ((AliTRDdigitizer &) d).fPRFwid             = fPRFwid;
 391   ((AliTRDdigitizer &) d).fPRFpad             = fPRFpad;
 392   if (((AliTRDdigitizer &) d).fPRFsmp) delete ((AliTRDdigitizer &) d).fPRFsmp;
 393   ((AliTRDdigitizer &) d).fPRFsmp = new Float_t[fPRFbin];
 394   for (iBin = 0; iBin < fPRFbin; iBin++) {
 395     ((AliTRDdigitizer &) d).fPRFsmp[iBin] = fPRFsmp[iBin];
 396   }
 397   ((AliTRDdigitizer &) d).fTRFbin             = fTRFbin;
 398   ((AliTRDdigitizer &) d).fTRFlo              = fTRFlo;
 399   ((AliTRDdigitizer &) d).fTRFhi              = fTRFhi;
 400   ((AliTRDdigitizer &) d).fTRFwid             = fTRFwid;
 401   if (((AliTRDdigitizer &) d).fTRFsmp) delete ((AliTRDdigitizer &) d).fTRFsmp;
 402   ((AliTRDdigitizer &) d).fTRFsmp = new Float_t[fTRFbin];
 403   for (iBin = 0; iBin < fTRFbin; iBin++) {
 404     ((AliTRDdigitizer &) d).fTRFsmp[iBin] = fTRFsmp[iBin];
 405   }
 406
 407 }
 408
 409 //_____________________________________________________________________________
 410 Int_t AliTRDdigitizer::Diffusion(Float_t driftlength, Float_t *xyz)
 411 {
 412   //
 413   // Applies the diffusion smearing to the position of a single electron
 414   //
 415
 416   Float_t driftSqrt = TMath::Sqrt(driftlength);
 417   Float_t sigmaT = driftSqrt * fDiffusionT;
 418   Float_t sigmaL = driftSqrt * fDiffusionL;
 419   xyz[0] = gRandom->Gaus(xyz[0], sigmaL * fLorentzFactor);
 420   xyz[1] = gRandom->Gaus(xyz[1], sigmaT * fLorentzFactor);
 421   xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
 422
 423   return 1;
 424
 425 }
 426
 427 //_____________________________________________________________________________
 428 Int_t AliTRDdigitizer::ExB(Float_t driftlength, Float_t *xyz)
 429 {
 430   //
 431   // Applies E x B effects to the position of a single electron
 432   //
 433
 434   xyz[0] = xyz[0];
 435   xyz[1] = xyz[1] + fOmegaTau * driftlength;
 436   xyz[2] = xyz[2];
 437
 438   return 1;
 439
 440 }
 441
 442 //_____________________________________________________________________________
 443 Int_t AliTRDdigitizer::PadResponse(Float_t signal, Float_t dist, Float_t *pad)
 444 {
 445   //
 446   // Applies the pad response
 447   //
 448
 449   Int_t iBin =  ((Int_t) (( - dist - fPRFlo) / fPRFwid));
 450
 451   Int_t iBin0 = iBin - fPRFpad;
 452   Int_t iBin1 = iBin;
 453   Int_t iBin2 = iBin + fPRFpad;
 454
 455   if ((iBin0 >= 0) && (iBin2 < fPRFbin)) {
 456
 457     pad[0] = signal * fPRFsmp[iBin0];
 458     pad[1] = signal * fPRFsmp[iBin1];
 459     pad[2] = signal * fPRFsmp[iBin2];
 460
 461     return 1;
 462
 463   }
 464   else {
 465
 466     return 0;
 467
 468   }
 469
 470 }
 471
 472 //_____________________________________________________________________________
 473 Float_t AliTRDdigitizer::TimeResponse(Float_t time)
 474 {
 475   //
 476   // Applies the preamp shaper time response
 477   //
 478
 479   Int_t iBin = ((Int_t) ((time - fTRFlo) / fTRFwid));
 480   if ((iBin >= 0) && (iBin < fTRFbin)) {
 481     return fTRFsmp[iBin];
 482   }
 483   else {
 484     return 0.0;
 485   }
 486
 487 }
 488
 489 //_____________________________________________________________________________
 490 Bool_t AliTRDdigitizer::Init()
 491 {
 492   //
 493   // Initializes the digitization procedure with standard values
 494   //
 495
 496   // The default parameter for the digitization
 497   fGasGain        = 2800.;
 498   fChipGain       = 6.1;
 499   fNoise          = 1000.;
 500   fADCoutRange    = 1023.;          // 10-bit ADC
 501   fADCinRange     = 1000.;          // 1V input range
 502   fADCthreshold   = 1;
 503
 504   // For the summable digits
 505   fSDigitsScale   = 100.;
 506
 507   // The drift velocity (cm / mus)
 508   fDriftVelocity  = 1.5;
 509
 510   // Diffusion on
 511   fDiffusionOn    = 1;
 512
 513   // E x B effects
 514   fExBOn          = 0;
 515
 516   // Propability for electron attachment
 517   fElAttachOn     = 0;
 518   fElAttachProp   = 0.0;
 519
 520   // The pad response function
 521   fPRFOn          =  1;
 522
 523   // The time response function
 524   fTRFOn          =  1;
 525
 526   // The pad coupling factor (same number as for the TPC)
 527   fPadCoupling    = 0.5;
 528
 529   // The time coupling factor (same number as for the TPC)
 530   fTimeCoupling   = 0.4;
 531
 532   return kTRUE;
 533
 534 }
 535
 536 //_____________________________________________________________________________
 537 Bool_t AliTRDdigitizer::ReInit()
 538 {
 539   //
 540   // Reinitializes the digitization procedure after a change in the parameter
 541   //
 542
 543   if (!fGeo) {
 544     printf("AliTRDdigitizer::ReInit -- ");
 545     printf("No geometry defined. Run InitDetector() first\n");
 546     return kFALSE;
 547   }
 548
 549   // Calculate the time bin width in ns
 550   fTimeBinWidth   = fGeo->GetTimeBinSize() / fDriftVelocity * 1000.0;
 551
 552   // The range and the binwidth for the sampled TRF
 553   fTRFbin = 100;
 554   // Start 0.2 mus before the signal
 555   fTRFlo  = -0.2 * fDriftVelocity;
 556   // End the maximum driftlength after the signal
 557   fTRFhi  = AliTRDgeometry::DrThick()
 558           + fGeo->GetTimeAfter() * fGeo->GetTimeBinSize();
 559   fTRFwid = (fTRFhi - fTRFlo) / ((Float_t) fTRFbin);
 560
 561   // Transverse and longitudinal diffusion coefficients (Xe/CO2)
 562   fDiffusionT     = GetDiffusionT(fDriftVelocity,fField);
 563   fDiffusionL     = GetDiffusionL(fDriftVelocity,fField);
 564
 565   // omega * tau.= tan(Lorentz-angle)
 566   fOmegaTau       = GetOmegaTau(fDriftVelocity,fField);
 567
 568   // The Lorentz factor
 569   if (fExBOn) {
 570     fLorentzFactor = 1.0 / (1.0 + fOmegaTau*fOmegaTau);
 571   }
 572   else {
 573     fLorentzFactor = 1.0;
 574   }
 575
 576   return kTRUE;
 577
 578 }
 579
 580 //_____________________________________________________________________________
 581 void AliTRDdigitizer::SampleTRF()
 582 {
 583   //
 584   // Samples the time response function
 585   // It is defined according to Vasiles simulation of the preamp shaper
 586   // output and includes the effect of the ion tail (based on Tariqs
 587   // Garfield simulation) and a shaping time of 125 ns FWHM
 588   //
 589
 590   Int_t   ipos1;
 591   Int_t   ipos2;
 592   Float_t diff;
 593
 594   const Int_t kNpasa = 200;
 595   Float_t time[kNpasa]   = {  -0.280000, -0.270000, -0.260000, -0.250000
 596                             , -0.240000, -0.230000, -0.220000, -0.210000
 597                             , -0.200000, -0.190000, -0.180000, -0.170000
 598                             , -0.160000, -0.150000, -0.140000, -0.130000
 599                             , -0.120000, -0.110000, -0.100000, -0.090000
 600                             , -0.080000, -0.070000, -0.060000, -0.050000
 601                             , -0.040000, -0.030000, -0.020000, -0.010000
 602                             , -0.000000,  0.010000,  0.020000,  0.030000
 603                             ,  0.040000,  0.050000,  0.060000,  0.070000
 604                             ,  0.080000,  0.090000,  0.100000,  0.110000
 605                             ,  0.120000,  0.130000,  0.140000,  0.150000
 606                             ,  0.160000,  0.170000,  0.180000,  0.190000
 607                             ,  0.200000,  0.210000,  0.220000,  0.230000
 608                             ,  0.240000,  0.250000,  0.260000,  0.270000
 609                             ,  0.280000,  0.290000,  0.300000,  0.310000
 610                             ,  0.320000,  0.330000,  0.340000,  0.350000
 611                             ,  0.360000,  0.370000,  0.380000,  0.390000
 612                             ,  0.400000,  0.410000,  0.420000,  0.430000
 613                             ,  0.440000,  0.450000,  0.460000,  0.470000
 614                             ,  0.480000,  0.490000,  0.500000,  0.510000
 615                             ,  0.520000,  0.530000,  0.540000,  0.550000
 616                             ,  0.560000,  0.570000,  0.580000,  0.590000
 617                             ,  0.600000,  0.610000,  0.620000,  0.630000
 618                             ,  0.640000,  0.650000,  0.660000,  0.670000
 619                             ,  0.680000,  0.690000,  0.700000,  0.710000
 620                             ,  0.720000,  0.730000,  0.740000,  0.750000
 621                             ,  0.760000,  0.770000,  0.780000,  0.790000
 622                             ,  0.800000,  0.810000,  0.820000,  0.830000
 623                             ,  0.840000,  0.850000,  0.860000,  0.870000
 624                             ,  0.880000,  0.890000,  0.900000,  0.910000
 625                             ,  0.920000,  0.930000,  0.940000,  0.950000
 626                             ,  0.960000,  0.970000,  0.980000,  0.990000
 627                             ,  1.000000,  1.010000,  1.020000,  1.030000
 628                             ,  1.040000,  1.050000,  1.060000,  1.070000
 629                             ,  1.080000,  1.090000,  1.100000,  1.110000
 630                             ,  1.120000,  1.130000,  1.140000,  1.150000
 631                             ,  1.160000,  1.170000,  1.180000,  1.190000
 632                             ,  1.200000,  1.210000,  1.220000,  1.230000
 633                             ,  1.240000,  1.250000,  1.260000,  1.270000
 634                             ,  1.280000,  1.290000,  1.300000,  1.310000
 635                             ,  1.320000,  1.330000,  1.340000,  1.350000
 636                             ,  1.360000,  1.370000,  1.380000,  1.390000
 637                             ,  1.400000,  1.410000,  1.420000,  1.430000
 638                             ,  1.440000,  1.450000,  1.460000,  1.470000
 639                             ,  1.480000,  1.490000,  1.500000,  1.510000
 640                             ,  1.520000,  1.530000,  1.540000,  1.550000
 641                             ,  1.560000,  1.570000,  1.580000,  1.590000
 642                             ,  1.600000,  1.610000,  1.620000,  1.630000
 643                             ,  1.640000,  1.650000,  1.660000,  1.670000
 644                             ,  1.680000,  1.690000,  1.700000,  1.710000 };
 645
 646   Float_t signal[kNpasa] = {   0.000000,  0.000000,  0.000000,  0.000000
 647                             ,  0.000000,  0.000000,  0.000000,  0.000000
 648                             ,  0.000000,  0.000000,  0.000000,  0.000000
 649                             ,  0.000000,  0.000000,  0.000000,  0.000098
 650                             ,  0.003071,  0.020056,  0.066053,  0.148346
 651                             ,  0.263120,  0.398496,  0.540226,  0.674436
 652                             ,  0.790977,  0.883083,  0.947744,  0.985714
 653                             ,  0.999248,  0.992105,  0.967669,  0.930827
 654                             ,  0.884586,  0.833083,  0.778571,  0.723684
 655                             ,  0.669173,  0.617293,  0.567669,  0.521805
 656                             ,  0.479699,  0.440977,  0.405639,  0.373985
 657                             ,  0.345526,  0.320038,  0.297256,  0.276917
 658                             ,  0.258797,  0.242632,  0.228195,  0.215301
 659                             ,  0.203759,  0.193383,  0.184023,  0.175564
 660                             ,  0.167895,  0.160940,  0.154549,  0.148722
 661                             ,  0.143308,  0.138346,  0.133722,  0.129398
 662                             ,  0.125376,  0.121617,  0.118045,  0.114699
 663                             ,  0.111541,  0.108571,  0.105714,  0.103008
 664                             ,  0.100414,  0.097970,  0.095602,  0.093346
 665                             ,  0.091165,  0.089060,  0.087068,  0.085150
 666                             ,  0.083308,  0.081541,  0.079812,  0.078158
 667                             ,  0.076541,  0.075000,  0.073496,  0.072068
 668                             ,  0.070677,  0.069286,  0.068008,  0.066729
 669                             ,  0.065489,  0.064286,  0.063120,  0.061992
 670                             ,  0.060902,  0.059850,  0.058797,  0.057820
 671                             ,  0.056842,  0.055902,  0.054962,  0.054060
 672                             ,  0.053158,  0.052293,  0.051466,  0.050639
 673                             ,  0.049850,  0.049060,  0.048308,  0.047556
 674                             ,  0.046842,  0.046128,  0.045451,  0.044774
 675                             ,  0.044098,  0.043459,  0.042820,  0.042218
 676                             ,  0.041617,  0.041015,  0.040451,  0.039887
 677                             ,  0.039323,  0.038797,  0.038271,  0.037744
 678                             ,  0.037237,  0.036744,  0.036259,  0.035786
 679                             ,  0.035323,  0.034872,  0.034429,  0.033996
 680                             ,  0.033575,  0.033162,  0.032756,  0.032361
 681                             ,  0.031974,  0.031594,  0.031222,  0.030857
 682                             ,  0.030496,  0.030143,  0.029793,  0.029451
 683                             ,  0.029109,  0.028774,  0.028444,  0.028113
 684                             ,  0.027793,  0.027477,  0.027165,  0.026861
 685                             ,  0.026564,  0.026271,  0.025981,  0.025699
 686                             ,  0.025421,  0.025147,  0.024880,  0.024613
 687                             ,  0.024353,  0.024094,  0.023842,  0.023590
 688                             ,  0.023346,  0.023102,  0.022865,  0.022628
 689                             ,  0.022398,  0.022173,  0.021951,  0.021733
 690                             ,  0.021519,  0.021308,  0.021098,  0.020891
 691                             ,  0.020688,  0.020485,  0.020286,  0.020090
 692                             ,  0.019895,  0.019707,  0.019519,  0.019335
 693                             ,  0.019150,  0.018974,  0.018797,  0.018624
 694                             ,  0.018451,  0.018282,  0.018113,  0.017947
 695                             ,  0.017782,  0.017617,  0.017455,  0.017297 };
 696
 697   if (fTRFsmp) delete fTRFsmp;
 698   fTRFsmp = new Float_t[fTRFbin];
 699
 700   Float_t loTRF    = TMath::Max(fTRFlo / fDriftVelocity,time[0]);
 701   Float_t hiTRF    = TMath::Min(fTRFhi / fDriftVelocity,time[kNpasa-1]);
 702   Float_t binWidth = (hiTRF - loTRF) / ((Float_t) fTRFbin);
 703
 704   // Take the linear interpolation
 705   for (Int_t iBin = 0; iBin < fTRFbin; iBin++) {
 706
 707     Float_t bin = (((Float_t) iBin) + 0.5) * binWidth + loTRF;
 708     ipos1 = ipos2 = 0;
 709     diff  = 0;
 710     do {
 711       diff = bin - time[ipos2++];
 712     } while (diff > 0);
 713     ipos2--;
 714     if (ipos2 > kNpasa) ipos2 = kNpasa - 1;
 715     ipos1 = ipos2 - 1;
 716
 717     fTRFsmp[iBin] = signal[ipos2]
 718                   + diff * (signal[ipos2] - signal[ipos1])
 719                          / (  time[ipos2] -   time[ipos1]);
 720
 721   }
 722
 723 }
 724
 725 //_____________________________________________________________________________
 726 void AliTRDdigitizer::SamplePRF()
 727 {
 728   //
 729   // Samples the pad response function
 730   //
 731
 732   const Int_t kPRFbin = 61;
 733   Float_t prf[kPRFbin] = { 0.002340, 0.003380, 0.004900, 0.007080, 0.010220
 734                          , 0.014740, 0.021160, 0.030230, 0.042800, 0.059830
 735                          , 0.082030, 0.109700, 0.142550, 0.179840, 0.220610
 736                          , 0.263980, 0.309180, 0.355610, 0.402790, 0.450350
 737                          , 0.497930, 0.545190, 0.591740, 0.637100, 0.680610
 738                          , 0.721430, 0.758400, 0.790090, 0.814720, 0.830480
 739                          , 0.835930, 0.830480, 0.814710, 0.790070, 0.758390
 740                          , 0.721410, 0.680590, 0.637080, 0.591730, 0.545180
 741                          , 0.497920, 0.450340, 0.402790, 0.355610, 0.309190
 742                          , 0.263990, 0.220630, 0.179850, 0.142570, 0.109720
 743                          , 0.082040, 0.059830, 0.042820, 0.030230, 0.021170
 744                          , 0.014740, 0.010230, 0.007080, 0.004900, 0.003380
 745                          , 0.002340 };
 746
 747   fPRFlo  = -1.5;
 748   fPRFhi  =  1.5;
 749   fPRFbin = kPRFbin;
 750   fPRFwid = (fPRFhi - fPRFlo) / ((Float_t) fPRFbin);
 751   fPRFpad = ((Int_t) (1.0 / fPRFwid));
 752
 753   if (fPRFsmp) delete fPRFsmp;
 754   fPRFsmp = new Float_t[fPRFbin];
 755   for (Int_t iBin = 0; iBin < fPRFbin; iBin++) {
 756     fPRFsmp[iBin] = prf[iBin];
 757   }
 758
 759 }
 760
 761 //_____________________________________________________________________________
 762 Bool_t AliTRDdigitizer::Open(const Char_t *name, Int_t nEvent)
 763 {
 764   //
 765   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
 766   //
 767
 768   // Connect the AliRoot file containing Geometry, Kine, and Hits
 769   fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
 770   if (!fInputFile) {
 771     if (fVerbose > 0) {
 772       printf("AliTRDdigitizer::Open -- ");
 773       printf("Open the AliROOT-file %s.\n",name);
 774     }
 775     fInputFile = new TFile(name,"UPDATE");
 776   }
 777   else {
 778     if (fVerbose > 0) {
 779       printf("AliTRDdigitizer::Open -- ");
 780       printf("%s is already open.\n",name);
 781     }
 782   }
 783
 784   gAlice = (AliRun*) fInputFile->Get("gAlice");
 785   if (gAlice) {
 786     if (fVerbose > 0) {
 787       printf("AliTRDdigitizer::Open -- ");
 788       printf("AliRun object found on file.\n");
 789     }
 790   }
 791   else {
 792     printf("AliTRDdigitizer::Open -- ");
 793     printf("Could not find AliRun object.\n");
 794     return kFALSE;
 795   }
 796
 797   fEvent = nEvent;
 798
 799   // Import the Trees for the event nEvent in the file
 800   Int_t nparticles = gAlice->GetEvent(fEvent);
 801   if (nparticles <= 0) {
 802     printf("AliTRDdigitizer::Open -- ");
 803     printf("No entries in the trees for event %d.\n",fEvent);
 804     return kFALSE;
 805   }
 806
 807   if (InitDetector()) {
 808     return MakeBranch();
 809   }
 810   else {
 811     return kFALSE;
 812   }
 813
 814 }
 815
 816 //_____________________________________________________________________________
 817 Bool_t AliTRDdigitizer::InitDetector()
 818 {
 819   //
 820   // Sets the pointer to the TRD detector and the geometry
 821   //
 822
 823   // Get the pointer to the detector class and check for version 1
 824   fTRD = (AliTRD*) gAlice->GetDetector("TRD");
 825   if (fTRD->IsVersion() != 1) {
 826     printf("AliTRDdigitizer::InitDetector -- ");
 827     printf("TRD must be version 1 (slow simulator).\n");
 828     exit(1);
 829   }
 830
 831   // Get the geometry
 832   fGeo = fTRD->GetGeometry();
 833   if (fVerbose > 0) {
 834     printf("AliTRDdigitizer::InitDetector -- ");
 835     printf("Geometry version %d\n",fGeo->IsVersion());
 836   }
 837
 838   // The magnetic field strength in Tesla
 839   fField = 0.2 * gAlice->Field()->Factor();
 840
 841   // Create a digits manager
 842   fDigitsManager = new AliTRDdigitsManager();
 843   fDigitsManager->SetSDigits(fSDigits);
 844   fDigitsManager->CreateArrays();
 845   fDigitsManager->SetEvent(fEvent);
 846   fDigitsManager->SetVerbose(fVerbose);
 847
 848   // The list for the input s-digits manager to be merged
 849   fSDigitsManagerList = new TList();
 850
 851   return ReInit();
 852
 853 }
 854
 855 //_____________________________________________________________________________
 856 Bool_t AliTRDdigitizer::MakeBranch(const Char_t *file)
 857 {
 858   //
 859   // Create the branches for the digits array
 860   //
 861
 862   return fDigitsManager->MakeBranch(file);
 863
 864 }
 865
 866 //_____________________________________________________________________________
 867 Bool_t AliTRDdigitizer::MakeDigits()
 868 {
 869   //
 870   // Creates digits.
 871   //
 872
 873   ///////////////////////////////////////////////////////////////
 874   // Parameter
 875   ///////////////////////////////////////////////////////////////
 876
 877   // Converts number of electrons to fC
 878   const Double_t kEl2fC  = 1.602E-19 * 1.0E15;
 879
 880   ///////////////////////////////////////////////////////////////
 881
 882   // Number of pads included in the pad response
 883   const Int_t kNpad  = 3;
 884
 885   // Number of track dictionary arrays
 886   const Int_t kNDict = AliTRDdigitsManager::kNDict;
 887
 888   // Half the width of the amplification region
 889   const Float_t kAmWidth = AliTRDgeometry::AmThick() / 2.;
 890
 891   Int_t   iRow, iCol, iTime, iPad;
 892   Int_t   iDict  = 0;
 893   Int_t   nBytes = 0;
 894
 895   Int_t   totalSizeDigits = 0;
 896   Int_t   totalSizeDict0  = 0;
 897   Int_t   totalSizeDict1  = 0;
 898   Int_t   totalSizeDict2  = 0;
 899
 900   Int_t   timeTRDbeg = 0;
 901   Int_t   timeTRDend = 1;
 902
 903   Float_t pos[3];
 904   Float_t rot[3];
 905   Float_t xyz[3];
 906   Float_t padSignal[kNpad];
 907   Float_t signalOld[kNpad];
 908
 909   AliTRDdataArrayF *signals = 0;
 910   AliTRDdataArrayI *digits  = 0;
 911   AliTRDdataArrayI *dictionary[kNDict];
 912
 913   // Create a container for the amplitudes
 914   AliTRDsegmentArray *signalsArray
 915                      = new AliTRDsegmentArray("AliTRDdataArrayF",AliTRDgeometry::Ndet());
 916
 917   if (fTRFOn) {
 918     timeTRDbeg = ((Int_t) (-fTRFlo / fGeo->GetTimeBinSize())) - 1;
 919     timeTRDend = ((Int_t) ( fTRFhi / fGeo->GetTimeBinSize())) - 1;
 920     if (fVerbose > 0) {
 921       printf("AliTRDdigitizer::MakeDigits -- ");
 922       printf("Sample the TRF between -%d and %d\n",timeTRDbeg,timeTRDend);
 923     }
 924   }
 925
 926   Float_t elAttachProp = fElAttachProp / 100.;
 927
 928   // Create the sampled PRF
 929   SamplePRF();
 930
 931   // Create the sampled TRF
 932   SampleTRF();
 933
 934   if (!fGeo) {
 935     printf("AliTRDdigitizer::MakeDigits -- ");
 936     printf("No geometry defined\n");
 937     return kFALSE;
 938   }
 939
 940   if (fVerbose > 0) {
 941     printf("AliTRDdigitizer::MakeDigits -- ");
 942     printf("Start creating digits.\n");
 943   }
 944
 945   // Get the pointer to the hit tree
 946   TTree *HitTree = gAlice->TreeH();
 947
 948   // Get the number of entries in the hit tree
 949   // (Number of primary particles creating a hit somewhere)
 950   Int_t nTrack = (Int_t) HitTree->GetEntries();
 951   if (fVerbose > 0) {
 952     printf("AliTRDdigitizer::MakeDigits -- ");
 953     printf("Found %d primary particles\n",nTrack);
 954   }
 955
 956   Int_t detectorOld = -1;
 957   Int_t countHits   =  0;
 958
 959   // Loop through all entries in the tree
 960   for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
 961
 962     gAlice->ResetHits();
 963     nBytes += HitTree->GetEvent(iTrack);
 964
 965     // Get the number of hits in the TRD created by this particle
 966     Int_t nHit = fTRD->Hits()->GetEntriesFast();
 967     if (fVerbose > 0) {
 968       printf("AliTRDdigitizer::MakeDigits -- ");
 969       printf("Found %d hits for primary particle %d\n",nHit,iTrack);
 970     }
 971
 972     // Loop through the TRD hits
 973     for (Int_t iHit = 0; iHit < nHit; iHit++) {
 974
 975       countHits++;
 976
 977       AliTRDhit *hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit);
 978               pos[0]   = hit->X();
 979               pos[1]   = hit->Y();
 980               pos[2]   = hit->Z();
 981       Float_t q        = hit->GetCharge();
 982       Int_t   track    = hit->Track();
 983       Int_t   detector = hit->GetDetector();
 984       Int_t   plane    = fGeo->GetPlane(detector);
 985       Int_t   sector   = fGeo->GetSector(detector);
 986       Int_t   chamber  = fGeo->GetChamber(detector);
 987
 988       if (!(CheckDetector(plane,chamber,sector))) continue;
 989
 990       Int_t   nRowMax     = fGeo->GetRowMax(plane,chamber,sector);
 991       Int_t   nColMax     = fGeo->GetColMax(plane);
 992       Int_t   nTimeMax    = fGeo->GetTimeMax();
 993       Int_t   nTimeBefore = fGeo->GetTimeBefore();
 994       Int_t   nTimeAfter  = fGeo->GetTimeAfter();
 995       Int_t   nTimeTotal  = fGeo->GetTimeTotal();
 996       Float_t row0        = fGeo->GetRow0(plane,chamber,sector);
 997       Float_t col0        = fGeo->GetCol0(plane);
 998       Float_t time0       = fGeo->GetTime0(plane);
 999       Float_t rowPadSize  = fGeo->GetRowPadSize(plane,chamber,sector);
1000       Float_t colPadSize  = fGeo->GetColPadSize(plane);
1001       Float_t timeBinSize = fGeo->GetTimeBinSize();
1002       Float_t divideRow   = 1.0 / rowPadSize;
1003       Float_t divideCol   = 1.0 / colPadSize;
1004       Float_t divideTime  = 1.0 / timeBinSize;
1005
1006       if (fVerbose > 1) {
1007         printf("Analyze hit no. %d ",iHit);
1008         printf("-----------------------------------------------------------\n");
1009         hit->Dump();
1010         printf("plane = %d, sector = %d, chamber = %d\n"
1011               ,plane,sector,chamber);
1012         printf("nRowMax = %d, nColMax = %d, nTimeMax = %d\n"
1013               ,nRowMax,nColMax,nTimeMax);
1014         printf("nTimeBefore = %d, nTimeAfter = %d, nTimeTotal = %d\n"
1015               ,nTimeBefore,nTimeAfter,nTimeTotal);
1016         printf("row0 = %f, col0 = %f, time0 = %f\n"
1017               ,row0,col0,time0);
1018         printf("rowPadSize = %f, colPadSize = %f, timeBinSize = %f\n"
1019                ,rowPadSize,colPadSize,timeBinSize);
1020       }
1021
1022       // Don't analyze test hits
1023       if (hit->FromTest()) continue;
1024
1025       if (detector != detectorOld) {
1026
1027         if (fVerbose > 1) {
1028           printf("AliTRDdigitizer::MakeDigits -- ");
1029           printf("Get new container. New det = %d, Old det = %d\n"
1030                 ,detector,detectorOld);
1031         }
1032         // Compress the old one if enabled
1033         if ((fCompress) && (detectorOld > -1)) {
1034           if (fVerbose > 1) {
1035             printf("AliTRDdigitizer::MakeDigits -- ");
1036             printf("Compress the old container ...");
1037           }
1038           signals->Compress(1,0);
1039           for (iDict = 0; iDict < kNDict; iDict++) {
1040             dictionary[iDict]->Compress(1,0);
1041           }
1042           if (fVerbose > 1) printf("done\n");
1043         }
1044         // Get the new container
1045         signals = (AliTRDdataArrayF *) signalsArray->At(detector);
1046         if (signals->GetNtime() == 0) {
1047           // Allocate a new one if not yet existing
1048           if (fVerbose > 1) {
1049             printf("AliTRDdigitizer::MakeDigits -- ");
1050             printf("Allocate a new container ... ");
1051           }
1052           signals->Allocate(nRowMax,nColMax,nTimeTotal);
1053         }
1054         else {
1055           // Expand an existing one
1056           if (fCompress) {
1057             if (fVerbose > 1) {
1058               printf("AliTRDdigitizer::MakeDigits -- ");
1059               printf("Expand an existing container ... ");
1060             }
1061             signals->Expand();
1062           }
1063         }
1064         // The same for the dictionary
1065         for (iDict = 0; iDict < kNDict; iDict++) {
1066           dictionary[iDict] = fDigitsManager->GetDictionary(detector,iDict);
1067           if (dictionary[iDict]->GetNtime() == 0) {
1068             dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1069           }
1070           else {
1071             if (fCompress) dictionary[iDict]->Expand();
1072           }
1073         }
1074         if (fVerbose > 1) printf("done\n");
1075         detectorOld = detector;
1076       }
1077
1078       // Rotate the sectors on top of each other
1079       fGeo->Rotate(detector,pos,rot);
1080
1081       // The driftlength. It is negative if the hit is in the
1082       // amplification region.
1083       Float_t driftlength = time0 - rot[0];
1084
1085       // Take also the drift in the amplification region into account
1086       // The drift length is at the moment still the same, regardless of
1087       // the position relativ to the wire. This non-isochronity needs still
1088       // to be implemented.
1089       Float_t driftlengthL = TMath::Abs(driftlength + kAmWidth);
1090       if (fExBOn) driftlengthL /= TMath::Sqrt(fLorentzFactor);
1091
1092       // Loop over all electrons of this hit
1093       // TR photons produce hits with negative charge
1094       Int_t nEl = ((Int_t) TMath::Abs(q));
1095       for (Int_t iEl = 0; iEl < nEl; iEl++) {
1096
1097         xyz[0] = rot[0];
1098         xyz[1] = rot[1];
1099         xyz[2] = rot[2];
1100
1101         // Electron attachment
1102         if (fElAttachOn) {
1103           if (gRandom->Rndm() < (driftlengthL * elAttachProp))
1104             continue;
1105         }
1106
1107         // Apply the diffusion smearing
1108         if (fDiffusionOn) {
1109           if (!(Diffusion(driftlengthL,xyz))) continue;
1110         }
1111
1112         // Apply E x B effects (depends on drift direction)
1113         if (fExBOn) {
1114           if (!(ExB(driftlength+kAmWidth,xyz))) continue;
1115         }
1116
1117         // The electron position after diffusion and ExB in pad coordinates
1118         // The pad row (z-direction)
1119         Int_t  rowE = ((Int_t) ((xyz[2] -  row0) * divideRow));
1120         if ((rowE < 0) || (rowE >= nRowMax)) continue;
1121
1122         // The pad column (rphi-direction)
1123         Int_t  colE = ((Int_t) ((xyz[1] -  col0) * divideCol));
1124         if ((colE < 0) || (colE >= nColMax)) continue;
1125
1126         // The time bin (negative for hits in the amplification region)
1127         // In the amplification region the electrons drift from both sides
1128         // to the middle (anode wire plane)
1129         Float_t timeDist   = time0 - xyz[0];
1130         Float_t timeOffset = 0;
1131         Int_t   timeE      = 0;
1132         if (timeDist > 0) {
1133           // The time bin
1134           timeE      = ((Int_t) (timeDist * divideTime));
1135           // The distance of the position to the middle of the timebin
1136           timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) - timeDist;
1137         }
1138         else {
1139           // Difference between half of the amplification gap width and
1140           // the distance to the anode wire
1141           Float_t anodeDist = kAmWidth - TMath::Abs(timeDist + kAmWidth);
1142           // The time bin
1143           timeE      = -1 * (((Int_t ) (anodeDist * divideTime)) + 1);
1144           // The distance of the position to the middle of the timebin
1145           timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) + anodeDist;
1146         }
1147
1148         // Apply the gas gain including fluctuations
1149         Float_t ggRndm = 0.0;
1150         do {
1151           ggRndm = gRandom->Rndm();
1152         } while (ggRndm <= 0);
1153         Int_t signal = (Int_t) (-fGasGain * TMath::Log(ggRndm));
1154
1155         // Apply the pad response
1156         if (fPRFOn) {
1157           // The distance of the electron to the center of the pad
1158           // in units of pad width
1159           Float_t dist = (xyz[1] - col0 - (colE + 0.5) * colPadSize)
1160                        * divideCol;
1161           if (!(PadResponse(signal,dist,padSignal))) continue;
1162         }
1163         else {
1164           padSignal[0] = 0.0;
1165           padSignal[1] = signal;
1166           padSignal[2] = 0.0;
1167         }
1168
1169         // Sample the time response inside the drift region
1170         // + additional time bins before and after.
1171         // The sampling is done always in the middle of the time bin
1172         for (Int_t iTimeBin = TMath::Max(timeE-timeTRDbeg,        -nTimeBefore)
1173                   ;iTimeBin < TMath::Min(timeE+timeTRDend,nTimeMax+nTimeAfter )
1174                   ;iTimeBin++) {
1175
1176           // Apply the time response
1177           Float_t timeResponse = 1.0;
1178           if (fTRFOn) {
1179             Float_t time = (iTimeBin - timeE) * timeBinSize + timeOffset;
1180             timeResponse = TimeResponse(time);
1181           }
1182
1183           signalOld[0] = 0.0;
1184           signalOld[1] = 0.0;
1185           signalOld[2] = 0.0;
1186
1187           for (iPad = 0; iPad < kNpad; iPad++) {
1188
1189             Int_t colPos = colE + iPad - 1;
1190             if (colPos <        0) continue;
1191             if (colPos >= nColMax) break;
1192
1193             // Add the signals
1194             // Note: The time bin number is shifted by nTimeBefore to avoid negative
1195             // time bins. This has to be subtracted later.
1196             Int_t iCurrentTimeBin = iTimeBin + nTimeBefore;
1197             signalOld[iPad]  = signals->GetDataUnchecked(rowE,colPos,iCurrentTimeBin);
1198             signalOld[iPad] += padSignal[iPad] * timeResponse;
1199             signals->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,signalOld[iPad]);
1200
1201             // Store the track index in the dictionary
1202             // Note: We store index+1 in order to allow the array to be compressed
1203             if (signalOld[iPad] > 0) {
1204               for (iDict = 0; iDict < kNDict; iDict++) {
1205                 Int_t oldTrack = dictionary[iDict]->GetDataUnchecked(rowE
1206                                                                     ,colPos
1207                                                                     ,iCurrentTimeBin);
1208                 if (oldTrack == track+1) break;
1209                 if (oldTrack ==       0) {
1210                   dictionary[iDict]->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,track+1);
1211                   break;
1212                 }
1213               }
1214             }
1215
1216           }
1217
1218         }
1219
1220       }
1221
1222     }
1223
1224   } // All hits finished
1225
1226   if (fVerbose > 0) {
1227     printf("AliTRDdigitizer::MakeDigits -- ");
1228     printf("Finished analyzing %d hits\n",countHits);
1229   }
1230
1231   // The total conversion factor
1232   Float_t convert = kEl2fC * fPadCoupling * fTimeCoupling * fChipGain;
1233
1234   // Loop through all chambers to finalize the digits
1235   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1236
1237     Int_t plane       = fGeo->GetPlane(iDet);
1238     Int_t sector      = fGeo->GetSector(iDet);
1239     Int_t chamber     = fGeo->GetChamber(iDet);
1240     Int_t nRowMax     = fGeo->GetRowMax(plane,chamber,sector);
1241     Int_t nColMax     = fGeo->GetColMax(plane);
1242     Int_t nTimeMax    = fGeo->GetTimeMax();
1243     Int_t nTimeTotal  = fGeo->GetTimeTotal();
1244
1245     if (fVerbose > 0) {
1246       printf("AliTRDdigitizer::MakeDigits -- ");
1247       printf("Digitization for chamber %d\n",iDet);
1248     }
1249
1250     // Add a container for the digits of this detector
1251     digits = fDigitsManager->GetDigits(iDet);
1252     // Allocate memory space for the digits buffer
1253     digits->Allocate(nRowMax,nColMax,nTimeTotal);
1254
1255     // Get the signal container
1256     signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
1257     if (signals->GetNtime() == 0) {
1258       // Create missing containers
1259       signals->Allocate(nRowMax,nColMax,nTimeTotal);
1260     }
1261     else {
1262       // Expand the container if neccessary
1263       if (fCompress) signals->Expand();
1264     }
1265     // Create the missing dictionary containers
1266     for (iDict = 0; iDict < kNDict; iDict++) {
1267       dictionary[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1268       if (dictionary[iDict]->GetNtime() == 0) {
1269         dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1270       }
1271     }
1272
1273     Int_t nDigits = 0;
1274
1275     // Don't create noise in detectors that are switched off
1276     if (CheckDetector(plane,chamber,sector)) {
1277
1278       // Create the digits for this chamber
1279       for (iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1280         for (iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1281           for (iTime = 0; iTime < nTimeTotal; iTime++) {
1282
1283             // Create summable digits
1284             if (fSDigits) {
1285
1286               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1287               signalAmp *= fSDigitsScale;
1288               signalAmp  = TMath::Min(signalAmp,1.0e9);
1289               Int_t adc  = (Int_t) signalAmp;
1290               nDigits++;
1291               digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1292
1293             }
1294             // Create normal digits
1295             else {
1296
1297               Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1298
1299               // Add the noise
1300               signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,fNoise),0.0);
1301               // Convert to mV
1302               signalAmp *= convert;
1303               // Convert to ADC counts. Set the overflow-bit fADCoutRange if the
1304               // signal is larger than fADCinRange
1305               Int_t adc  = 0;
1306               if (signalAmp >= fADCinRange) {
1307                 adc = ((Int_t) fADCoutRange);
1308               }
1309               else {
1310                 adc = ((Int_t) (signalAmp * (fADCoutRange / fADCinRange)));
1311               }
1312
1313               // Store the amplitude of the digit if above threshold
1314               if (adc > fADCthreshold) {
1315                 if (fVerbose > 2) {
1316                   printf("  iRow = %d, iCol = %d, iTime = %d\n"
1317                         ,iRow,iCol,iTime);
1318                   printf("  signal = %f, adc = %d\n",signalAmp,adc);
1319                 }
1320                 nDigits++;
1321                 digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1322               }
1323
1324             }
1325
1326           }
1327         }
1328       }
1329
1330     }
1331
1332     // Compress the arrays
1333     digits->Compress(1,0);
1334     for (iDict = 0; iDict < kNDict; iDict++) {
1335       dictionary[iDict]->Compress(1,0);
1336     }
1337
1338     totalSizeDigits += digits->GetSize();
1339     totalSizeDict0  += dictionary[0]->GetSize();
1340     totalSizeDict1  += dictionary[1]->GetSize();
1341     totalSizeDict2  += dictionary[2]->GetSize();
1342
1343     Float_t nPixel = nRowMax * nColMax * nTimeMax;
1344     if (fVerbose > 0) {
1345       printf("AliTRDdigitizer::MakeDigits -- ");
1346       printf("Found %d digits in detector %d (%3.0f).\n"
1347             ,nDigits,iDet
1348             ,100.0 * ((Float_t) nDigits) / nPixel);
1349     }
1350
1351     if (fCompress) signals->Compress(1,0);
1352
1353   }
1354
1355   if (fVerbose > 0) {
1356     printf("AliTRDdigitizer::MakeDigits -- ");
1357     printf("Total number of analyzed hits = %d\n",countHits);
1358     printf("AliTRDdigitizer::MakeDigits -- ");
1359     printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
1360                                                       ,totalSizeDict0
1361                                                       ,totalSizeDict1
1362                                                       ,totalSizeDict2);
1363   }
1364
1365   return kTRUE;
1366
1367 }
1368
1369 //_____________________________________________________________________________
1370 void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
1371 {
1372   //
1373   // Add a digits manager for s-digits to the input list.
1374   //
1375
1376   fSDigitsManagerList->Add(man);
1377
1378 }
1379
1380 //_____________________________________________________________________________
1381 Bool_t AliTRDdigitizer::ConvertSDigits()
1382 {
1383   //
1384   // Converts s-digits to normal digits
1385   //
1386
1387   // Number of track dictionary arrays
1388   const Int_t    kNDict = AliTRDdigitsManager::kNDict;
1389
1390   // Converts number of electrons to fC
1391   const Double_t kEl2fC = 1.602E-19 * 1.0E15;
1392
1393   Int_t iDict = 0;
1394
1395   if (fVerbose > 0) {
1396     this->Dump();
1397   }
1398
1399   Double_t sDigitsScale = 1.0 / GetSDigitsScale();
1400   Double_t noise        = GetNoise();
1401   Double_t padCoupling  = GetPadCoupling();
1402   Double_t timeCoupling = GetTimeCoupling();
1403   Double_t chipGain     = GetChipGain();
1404   Double_t convert      = kEl2fC * padCoupling * timeCoupling * chipGain;;
1405   Double_t adcInRange   = GetADCinRange();
1406   Double_t adcOutRange  = GetADCoutRange();
1407   Int_t    adcThreshold = GetADCthreshold();
1408
1409   AliTRDdataArrayI *digitsIn;
1410   AliTRDdataArrayI *digitsOut;
1411   AliTRDdataArrayI *dictionaryIn[kNDict];
1412   AliTRDdataArrayI *dictionaryOut[kNDict];
1413
1414   // Loop through the detectors
1415   for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1416
1417     if (fVerbose > 0) {
1418       printf("AliTRDdigitizer::ConvertSDigits -- ");
1419       printf("Convert detector %d to digits.\n",iDet);
1420     }
1421
1422     Int_t plane      = fGeo->GetPlane(iDet);
1423     Int_t sector     = fGeo->GetSector(iDet);
1424     Int_t chamber    = fGeo->GetChamber(iDet);
1425     Int_t nRowMax    = fGeo->GetRowMax(plane,chamber,sector);
1426     Int_t nColMax    = fGeo->GetColMax(plane);
1427     Int_t nTimeTotal = fGeo->GetTimeTotal();
1428
1429     digitsIn  = fSDigitsManager->GetDigits(iDet);
1430     digitsIn->Expand();
1431     digitsOut = fDigitsManager->GetDigits(iDet);
1432     digitsOut->Allocate(nRowMax,nColMax,nTimeTotal);
1433     for (iDict = 0; iDict < kNDict; iDict++) {
1434       dictionaryIn[iDict]  = fSDigitsManager->GetDictionary(iDet,iDict);
1435       dictionaryIn[iDict]->Expand();
1436       dictionaryOut[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1437       dictionaryOut[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1438     }
1439
1440     for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1441       for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1442         for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
1443
1444           Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1445           signal *= sDigitsScale;
1446           // Add the noise
1447           signal  = TMath::Max((Double_t) gRandom->Gaus(signal,noise),0.0);
1448           // Convert to mV
1449           signal *= convert;
1450           // Convert to ADC counts. Set the overflow-bit adcOutRange if the
1451           // signal is larger than adcInRange
1452           Int_t adc  = 0;
1453           if (signal >= adcInRange) {
1454             adc = ((Int_t) adcOutRange);
1455           }
1456           else {
1457             adc = ((Int_t) (signal * (adcOutRange / adcInRange)));
1458           }
1459           // Store the amplitude of the digit if above threshold
1460           if (adc > adcThreshold) {
1461             digitsOut->SetDataUnchecked(iRow,iCol,iTime,adc);
1462           }
1463           // Copy the dictionary
1464           for (iDict = 0; iDict < kNDict; iDict++) {
1465             Int_t track = dictionaryIn[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1466             dictionaryOut[iDict]->SetDataUnchecked(iRow,iCol,iTime,track);
1467           }
1468
1469         }
1470       }
1471     }
1472
1473     if (fCompress) {
1474       digitsIn->Compress(1,0);
1475       digitsOut->Compress(1,0);
1476       for (iDict = 0; iDict < kNDict; iDict++) {
1477         dictionaryIn[iDict]->Compress(1,0);
1478         dictionaryOut[iDict]->Compress(1,0);
1479       }
1480     }
1481
1482   }
1483
1484   return kTRUE;
1485
1486 }
1487
1488 //_____________________________________________________________________________
1489 Bool_t AliTRDdigitizer::MergeSDigits()
1490 {
1491   //
1492   // Merges the input s-digits:
1493   //   - The amplitude of the different inputs are summed up.
1494   //   - Of the track IDs from the input dictionaries only one is
1495   //     kept for each input. This works for maximal 3 different merged inputs.
1496   //
1497
1498   // Number of track dictionary arrays
1499   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1500
1501   Int_t iDict = 0;
1502
1503   AliTRDdataArrayI *digitsA;
1504   AliTRDdataArrayI *digitsB;
1505   AliTRDdataArrayI *dictionaryA[kNDict];
1506   AliTRDdataArrayI *dictionaryB[kNDict];
1507
1508   // Get the first s-digits
1509   fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
1510   if (!fSDigitsManager) return kFALSE;
1511
1512   // Loop through the other sets of s-digits
1513   AliTRDdigitsManager *mergeSDigitsManager;
1514   mergeSDigitsManager = (AliTRDdigitsManager *)
1515                         fSDigitsManagerList->After(fSDigitsManager);
1516
1517   if (fVerbose > 0) {
1518     if (mergeSDigitsManager) {
1519       printf("AliTRDdigitizer::MergeSDigits -- ");
1520       printf("Merge serveral input files.\n");
1521     }
1522     else {
1523       printf("AliTRDdigitizer::MergeSDigits -- ");
1524       printf("Only one input file.\n");
1525     }
1526   }
1527
1528   Int_t iMerge = 0;
1529   while (mergeSDigitsManager) {
1530
1531     iMerge++;
1532
1533     // Loop through the detectors
1534     for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1535
1536       Int_t plane      = fGeo->GetPlane(iDet);
1537       Int_t sector     = fGeo->GetSector(iDet);
1538       Int_t chamber    = fGeo->GetChamber(iDet);
1539       Int_t nRowMax    = fGeo->GetRowMax(plane,chamber,sector);
1540       Int_t nColMax    = fGeo->GetColMax(plane);
1541       Int_t nTimeTotal = fGeo->GetTimeTotal();
1542
1543       // Loop through the pixels of one detector and add the signals
1544       digitsA = fSDigitsManager->GetDigits(iDet);
1545       digitsB = mergeSDigitsManager->GetDigits(iDet);
1546       digitsA->Expand();
1547       digitsB->Expand();
1548       for (iDict = 0; iDict < kNDict; iDict++) {
1549         dictionaryA[iDict] = fSDigitsManager->GetDictionary(iDet,iDict);
1550         dictionaryB[iDict] = mergeSDigitsManager->GetDictionary(iDet,iDict);
1551         dictionaryA[iDict]->Expand();
1552         dictionaryB[iDict]->Expand();
1553       }
1554
1555       if (fVerbose > 0) {
1556         printf("AliTRDdigitizer::MergeSDigits -- ");
1557         printf("Merge detector %d of input no.%d.\n",iDet,iMerge);
1558       }
1559
1560       for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1561         for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1562           for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
1563
1564             // Add the amplitudes of the summable digits
1565             Int_t ampA = digitsA->GetDataUnchecked(iRow,iCol,iTime);
1566             Int_t ampB = digitsB->GetDataUnchecked(iRow,iCol,iTime);
1567             ampA += ampB;
1568             digitsA->SetDataUnchecked(iRow,iCol,iTime,ampA);
1569
1570             // Take only one track from each input
1571             Int_t track = dictionaryB[0]->GetDataUnchecked(iRow,iCol,iTime);
1572             if (iMerge < kNDict) {
1573               dictionaryA[iMerge]->SetDataUnchecked(iRow,iCol,iTime,track);
1574             }
1575
1576           }
1577         }
1578       }
1579
1580       if (fCompress) {
1581         digitsA->Compress(1,0);
1582         digitsB->Compress(1,0);
1583         for (iDict = 0; iDict < kNDict; iDict++) {
1584           dictionaryA[iDict]->Compress(1,0);
1585           dictionaryB[iDict]->Compress(1,0);
1586         }
1587       }
1588
1589     }
1590
1591     // The next set of s-digits
1592     mergeSDigitsManager = (AliTRDdigitsManager *)
1593                           fSDigitsManagerList->After(mergeSDigitsManager);
1594
1595   }
1596
1597   return kTRUE;
1598
1599 }
1600
1601 //_____________________________________________________________________________
1602 Bool_t AliTRDdigitizer::SDigits2Digits()
1603 {
1604   //
1605   // Merges the input s-digits and converts them to normal digits
1606   //
1607
1608   if (!MergeSDigits()) return kFALSE;
1609
1610   return ConvertSDigits();
1611
1612 }
1613
1614 //_____________________________________________________________________________
1615 Bool_t AliTRDdigitizer::CheckDetector(Int_t plane, Int_t chamber, Int_t sector)
1616 {
1617   //
1618   // Checks whether a detector is enabled
1619   //
1620
1621   if ((fTRD->GetSensChamber() >=       0) &&
1622       (fTRD->GetSensChamber() != chamber)) return kFALSE;
1623   if ((fTRD->GetSensPlane()   >=       0) &&
1624       (fTRD->GetSensPlane()   !=   plane)) return kFALSE;
1625   if ( fTRD->GetSensSector()  >=       0) {
1626     Int_t sens1 = fTRD->GetSensSector();
1627     Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
1628     sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
1629            * AliTRDgeometry::Nsect();
1630     if (sens1 < sens2) {
1631       if ((sector < sens1) || (sector >= sens2)) return kFALSE;
1632     }
1633     else {
1634       if ((sector < sens1) && (sector >= sens2)) return kFALSE;
1635     }
1636   }
1637
1638   return kTRUE;
1639
1640 }
1641
1642 //_____________________________________________________________________________
1643 Bool_t AliTRDdigitizer::WriteDigits()
1644 {
1645   //
1646   // Writes out the TRD-digits and the dictionaries
1647   //
1648
1649   // Store the digits and the dictionary in the tree
1650   return fDigitsManager->WriteDigits();
1651
1652 }
1653
1654 //_____________________________________________________________________________
1655 Float_t AliTRDdigitizer::GetDiffusionL(Float_t vd, Float_t b)
1656 {
1657   //
1658   // Returns the longitudinal diffusion coefficient for a given drift
1659   // velocity <vd> and a B-field <b> for Xe/CO2 (15%).
1660   // The values are according to a GARFIELD simulation.
1661   //
1662
1663   const Int_t kNb = 5;
1664   Float_t p0[kNb] = {  0.007440,  0.007493,  0.007513,  0.007672,  0.007831 };
1665   Float_t p1[kNb] = {  0.019252,  0.018912,  0.018636,  0.018012,  0.017343 };
1666   Float_t p2[kNb] = { -0.005042, -0.004926, -0.004867, -0.004650, -0.004424 };
1667   Float_t p3[kNb] = {  0.000195,  0.000189,  0.000195,  0.000182,  0.000169 };
1668
1669   Int_t ib = ((Int_t) (10 * (b - 0.15)));
1670   ib       = TMath::Max(  0,ib);
1671   ib       = TMath::Min(kNb,ib);
1672
1673   Float_t diff = p0[ib]
1674                + p1[ib] * vd
1675                + p2[ib] * vd*vd
1676                + p3[ib] * vd*vd*vd;
1677
1678   return diff;
1679
1680 }
1681
1682 //_____________________________________________________________________________
1683 Float_t AliTRDdigitizer::GetDiffusionT(Float_t vd, Float_t b)
1684 {
1685   //
1686   // Returns the transverse diffusion coefficient for a given drift
1687   // velocity <vd> and a B-field <b> for Xe/CO2 (15%).
1688   // The values are according to a GARFIELD simulation.
1689   //
1690
1691   const Int_t kNb = 5;
1692   Float_t p0[kNb] = {  0.009550,  0.009599,  0.009674,  0.009757,  0.009850 };
1693   Float_t p1[kNb] = {  0.006667,  0.006539,  0.006359,  0.006153,  0.005925 };
1694   Float_t p2[kNb] = { -0.000853, -0.000798, -0.000721, -0.000635, -0.000541 };
1695   Float_t p3[kNb] = {  0.000131,  0.000122,  0.000111,  0.000098,  0.000085 };
1696
1697   Int_t ib = ((Int_t) (10 * (b - 0.15)));
1698   ib       = TMath::Max(  0,ib);
1699   ib       = TMath::Min(kNb,ib);
1700
1701   Float_t diff = p0[ib]
1702                + p1[ib] * vd
1703                + p2[ib] * vd*vd
1704                + p3[ib] * vd*vd*vd;
1705
1706   return diff;
1707
1708 }
1709
1710 //_____________________________________________________________________________
1711 Float_t AliTRDdigitizer::GetOmegaTau(Float_t vd, Float_t b)
1712 {
1713   //
1714   // Returns omega*tau (tan(Lorentz-angle)) for a given drift velocity <vd>
1715   // and a B-field <b> for Xe/CO2 (15%).
1716   // The values are according to a GARFIELD simulation.
1717   //
1718
1719   const Int_t kNb = 5;
1720   Float_t p0[kNb] = {  0.004810,  0.007412,  0.010252,  0.013409,  0.016888 };
1721   Float_t p1[kNb] = {  0.054875,  0.081534,  0.107333,  0.131983,  0.155455 };
1722   Float_t p2[kNb] = { -0.008682, -0.012896, -0.016987, -0.020880, -0.024623 };
1723   Float_t p3[kNb] = {  0.000155,  0.000238,  0.000330,  0.000428,  0.000541 };
1724
1725   Int_t ib = ((Int_t) (10 * (b - 0.15)));
1726   ib       = TMath::Max(  0,ib);
1727   ib       = TMath::Min(kNb,ib);
1728
1729   Float_t alphaL = p0[ib]
1730                  + p1[ib] * vd
1731                  + p2[ib] * vd*vd
1732                  + p3[ib] * vd*vd*vd;
1733
1734   return TMath::Tan(alphaL);
1735
1736 }
1737