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