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