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