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