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