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