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