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