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