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