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