Dimitri just makes it work
[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$
eda4336d 18Revision 1.10 2000/10/05 07:27:53 cblume
19Changes in the header-files by FCA
20
6798b56e 21Revision 1.9 2000/10/02 21:28:19 fca
22Removal of useless dependecies via forward declarations
23
94de3818 24Revision 1.8 2000/06/09 11:10:07 cblume
25Compiler warnings and coding conventions, next round
26
dd9a6ee3 27Revision 1.7 2000/06/08 18:32:58 cblume
28Make code compliant to coding conventions
29
8230f242 30Revision 1.6 2000/06/07 16:27:32 cblume
31Try to remove compiler warnings on Sun and HP
32
9d0b222b 33Revision 1.5 2000/05/09 16:38:57 cblume
34Removed PadResponse(). Merge problem
35
f0a7bf65 36Revision 1.4 2000/05/08 15:53:45 cblume
37Resolved merge conflict
38
da581aea 39Revision 1.3 2000/04/28 14:49:27 cblume
40Only one declaration of iDict in MakeDigits()
41
42Revision 1.1.4.1 2000/05/08 14:42:04 cblume
43Introduced AliTRDdigitsManager
28329a48 44
1befd3b2 45Revision 1.1 2000/02/28 19:00:13 cblume
46Add new TRD classes
47
f7336fa3 48*/
49
50///////////////////////////////////////////////////////////////////////////////
51// //
52// Creates and handles digits from TRD hits //
53// //
54// The following effects are included: //
55// - Diffusion //
56// - ExB effects //
57// - Gas gain including fluctuations //
58// - Pad-response (simple Gaussian approximation) //
59// - Electronics noise //
60// - Electronics gain //
61// - Digitization //
62// - ADC threshold //
63// The corresponding parameter can be adjusted via the various //
64// Set-functions. If these parameters are not explicitly set, default //
65// values are used (see Init-function). //
66// To produce digits from a root-file with TRD-hits use the //
67// slowDigitsCreate.C macro. //
68// //
69///////////////////////////////////////////////////////////////////////////////
70
6798b56e 71#include <stdlib.h>
72
f7336fa3 73#include <TMath.h>
74#include <TVector.h>
75#include <TRandom.h>
94de3818 76#include <TROOT.h>
77#include <TTree.h>
f7336fa3 78
79#include "AliTRD.h"
80#include "AliTRDdigitizer.h"
da581aea 81#include "AliTRDdataArrayI.h"
82#include "AliTRDdataArrayF.h"
83#include "AliTRDdigitsManager.h"
f7336fa3 84
85ClassImp(AliTRDdigitizer)
86
87//_____________________________________________________________________________
88AliTRDdigitizer::AliTRDdigitizer():TNamed()
89{
90 //
91 // AliTRDdigitizer default constructor
92 //
93
da581aea 94 fInputFile = NULL;
95 fDigits = NULL;
96 fTRD = NULL;
97 fGeo = NULL;
98 fPRF = NULL;
99
100 fEvent = 0;
101 fGasGain = 0.0;
102 fNoise = 0.0;
103 fChipGain = 0.0;
104 fADCoutRange = 0.0;
105 fADCinRange = 0.0;
106 fADCthreshold = 0;
107 fDiffusionOn = 0;
108 fDiffusionT = 0.0;
109 fDiffusionL = 0.0;
110 fElAttachOn = 0;
111 fElAttachProp = 0.0;
112 fExBOn = 0;
113 fLorentzAngle = 0.0;
f7336fa3 114
115}
116
117//_____________________________________________________________________________
118AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
119 :TNamed(name,title)
120{
121 //
122 // AliTRDdigitizer default constructor
123 //
124
da581aea 125 fInputFile = NULL;
126 fDigits = NULL;
127 fTRD = NULL;
128 fGeo = NULL;
f7336fa3 129
da581aea 130 fEvent = 0;
f7336fa3 131
132 Init();
133
134}
135
136//_____________________________________________________________________________
dd9a6ee3 137AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
8230f242 138{
139 //
140 // AliTRDdigitizer copy constructor
141 //
142
dd9a6ee3 143 ((AliTRDdigitizer &) d).Copy(*this);
8230f242 144
145}
146
147//_____________________________________________________________________________
f7336fa3 148AliTRDdigitizer::~AliTRDdigitizer()
149{
8230f242 150 //
151 // AliTRDdigitizer destructor
152 //
f7336fa3 153
154 if (fInputFile) {
155 fInputFile->Close();
156 delete fInputFile;
157 }
158
da581aea 159 if (fDigits) {
160 delete fDigits;
f7336fa3 161 }
162
da581aea 163 if (fPRF) delete fPRF;
f7336fa3 164
165}
166
167//_____________________________________________________________________________
dd9a6ee3 168AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
169{
170 //
171 // Assignment operator
172 //
173
174 if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
175 return *this;
176
177}
178
179//_____________________________________________________________________________
180void AliTRDdigitizer::Copy(TObject &d)
8230f242 181{
182 //
183 // Copy function
184 //
185
dd9a6ee3 186 ((AliTRDdigitizer &) d).fInputFile = NULL;
187 ((AliTRDdigitizer &) d).fDigits = NULL;
188 ((AliTRDdigitizer &) d).fTRD = NULL;
189 ((AliTRDdigitizer &) d).fGeo = NULL;
190
191 ((AliTRDdigitizer &) d).fEvent = 0;
192
193 ((AliTRDdigitizer &) d).fGasGain = fGasGain;
194 ((AliTRDdigitizer &) d).fNoise = fNoise;
195 ((AliTRDdigitizer &) d).fChipGain = fChipGain;
196 ((AliTRDdigitizer &) d).fADCoutRange = fADCoutRange;
197 ((AliTRDdigitizer &) d).fADCinRange = fADCinRange;
198 ((AliTRDdigitizer &) d).fADCthreshold = fADCthreshold;
199 ((AliTRDdigitizer &) d).fDiffusionOn = fDiffusionOn;
200 ((AliTRDdigitizer &) d).fDiffusionT = fDiffusionT;
201 ((AliTRDdigitizer &) d).fDiffusionL = fDiffusionL;
202 ((AliTRDdigitizer &) d).fElAttachOn = fElAttachOn;
203 ((AliTRDdigitizer &) d).fElAttachProp = fElAttachProp;
204 ((AliTRDdigitizer &) d).fExBOn = fExBOn;
205 ((AliTRDdigitizer &) d).fLorentzAngle = fLorentzAngle;
206 ((AliTRDdigitizer &) d).fLorentzFactor = fLorentzFactor;
207
208 fPRF->Copy(*((AliTRDdigitizer &) d).fPRF);
8230f242 209
210}
211
212//_____________________________________________________________________________
f7336fa3 213Int_t AliTRDdigitizer::Diffusion(Float_t driftlength, Float_t *xyz)
214{
215 //
216 // Applies the diffusion smearing to the position of a single electron
217 //
218
219 Float_t driftSqrt = TMath::Sqrt(driftlength);
220 Float_t sigmaT = driftSqrt * fDiffusionT;
221 Float_t sigmaL = driftSqrt * fDiffusionL;
222 xyz[0] = gRandom->Gaus(xyz[0], sigmaL * fLorentzFactor);
223 xyz[1] = gRandom->Gaus(xyz[1], sigmaT * fLorentzFactor);
224 xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
225 return 1;
226
227}
228
229//_____________________________________________________________________________
230Int_t AliTRDdigitizer::ExB(Float_t driftlength, Float_t *xyz)
231{
232 //
233 // Applies E x B effects to the position of a single electron
234 //
235
236 xyz[0] = xyz[0];
237 xyz[1] = xyz[1] + fLorentzAngle * driftlength;
238 xyz[2] = xyz[2];
239
240 return 1;
241
242}
243
244//_____________________________________________________________________________
245void AliTRDdigitizer::Init()
246{
247 //
248 // Initializes the digitization procedure with standard values
249 //
250
251 // The default parameter for the digitization
252 fGasGain = 2.0E3;
253 fNoise = 3000.;
254 fChipGain = 10.;
255 fADCoutRange = 255.;
256 fADCinRange = 2000.;
257 fADCthreshold = 1;
258
259 // Transverse and longitudinal diffusion coefficients (Xe/Isobutane)
260 fDiffusionOn = 1;
261 fDiffusionT = 0.060;
262 fDiffusionL = 0.017;
263
264 // Propability for electron attachment
265 fElAttachOn = 0;
266 fElAttachProp = 0.0;
267
268 // E x B effects
269 fExBOn = 0;
270 // omega * tau. (tau ~ 12 * 10^-12, B = 0.2T)
271 fLorentzAngle = 17.6 * 12.0 * 0.2 * 0.01;
272
da581aea 273 // The pad response function
274 fPRF = new TF1("PRF","[0]*([1]+exp(-x*x/(2.0*[2])))",-2,2);
275 fPRF->SetParameter(0, 0.8872);
276 fPRF->SetParameter(1,-0.00573);
277 fPRF->SetParameter(2, 0.454 * 0.454);
278
f7336fa3 279}
280
281//_____________________________________________________________________________
282Bool_t AliTRDdigitizer::Open(const Char_t *name, Int_t nEvent)
283{
284 //
285 // Opens a ROOT-file with TRD-hits and reads in the hit-tree
286 //
287
288 // Connect the AliRoot file containing Geometry, Kine, and Hits
289 fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
290 if (!fInputFile) {
291 printf("AliTRDdigitizer::Open -- ");
292 printf("Open the ALIROOT-file %s.\n",name);
293 fInputFile = new TFile(name,"UPDATE");
294 }
295 else {
296 printf("AliTRDdigitizer::Open -- ");
297 printf("%s is already open.\n",name);
298 }
299
da581aea 300 gAlice = (AliRun*) fInputFile->Get("gAlice");
301 if (gAlice) {
302 printf("AliTRDdigitizer::Open -- ");
303 printf("AliRun object found on file.\n");
304 }
305 else {
306 printf("AliTRDdigitizer::Open -- ");
307 printf("Could not find AliRun object.\n");
308 return kFALSE;
309 }
f7336fa3 310
311 fEvent = nEvent;
312
313 // Import the Trees for the event nEvent in the file
314 Int_t nparticles = gAlice->GetEvent(fEvent);
315 if (nparticles <= 0) {
316 printf("AliTRDdigitizer::Open -- ");
317 printf("No entries in the trees for event %d.\n",fEvent);
318 return kFALSE;
319 }
320
dd9a6ee3 321 // Get the pointer to the detector class and check for version 1
322 fTRD = (AliTRD*) gAlice->GetDetector("TRD");
323 if (fTRD->IsVersion() != 1) {
324 printf("AliTRDdigitizer::Open -- ");
325 printf("TRD must be version 1 (slow simulator).\n");
326 exit(1);
327 }
328
329 // Get the geometry
330 fGeo = fTRD->GetGeometry();
331 printf("AliTRDdigitizer::Open -- ");
332 printf("Geometry version %d\n",fGeo->IsVersion());
333
f7336fa3 334 return kTRUE;
335
336}
337
338//_____________________________________________________________________________
f7336fa3 339Bool_t AliTRDdigitizer::MakeDigits()
340{
341 //
342 // Loops through the TRD-hits and creates the digits.
343 //
344
f7336fa3 345 ///////////////////////////////////////////////////////////////
346 // Parameter
347 ///////////////////////////////////////////////////////////////
348
349 // Converts number of electrons to fC
8230f242 350 const Float_t kEl2fC = 1.602E-19 * 1.0E15;
f7336fa3 351
352 ///////////////////////////////////////////////////////////////
353
354 Int_t iRow, iCol, iTime;
355 Int_t nBytes = 0;
28329a48 356 Int_t iDict;
f7336fa3 357
358 Int_t totalSizeDigits = 0;
359 Int_t totalSizeDict0 = 0;
360 Int_t totalSizeDict1 = 0;
361 Int_t totalSizeDict2 = 0;
362
8230f242 363 AliTRDdataArrayI *digits;
364 AliTRDdataArrayI *dictionary[kNDict];
da581aea 365
366 if (!fGeo) {
367 printf("AliTRDdigitizer::MakeDigits -- ");
368 printf("No geometry defined\n");
369 return kFALSE;
370 }
371
372 // Create a digits manager
373 fDigits = new AliTRDdigitsManager();
374
375 // Create detector arrays to keep the signal and track numbers
8230f242 376 AliTRDdataArrayF *signal = new AliTRDdataArrayF();
377 AliTRDdataArrayI *tracks[kNDict];
da581aea 378 for (iDict = 0; iDict < kNDict; iDict++) {
8230f242 379 tracks[iDict] = new AliTRDdataArrayI();
da581aea 380 }
f7336fa3 381
382 // Get the pointer to the hit tree
8230f242 383 TTree *hitTree = gAlice->TreeH();
da581aea 384
385 // Get the number of entries in the hit tree
386 // (Number of primary particles creating a hit somewhere)
8230f242 387 Int_t nTrack = (Int_t) hitTree->GetEntries();
da581aea 388
389 printf("AliTRDdigitizer::MakeDigits -- ");
390 printf("Start creating digits.\n");
f7336fa3 391
392 // The Lorentz factor
393 if (fExBOn) {
394 fLorentzFactor = 1.0 / (1.0 + fLorentzAngle*fLorentzAngle);
395 }
396 else {
397 fLorentzFactor = 1.0;
398 }
399
f7336fa3 400 Int_t chamBeg = 0;
401 Int_t chamEnd = kNcham;
9d0b222b 402 if (fTRD->GetSensChamber() >= 0) {
da581aea 403 chamBeg = fTRD->GetSensChamber();
404 chamEnd = chamBeg + 1;
f7336fa3 405 }
406 Int_t planBeg = 0;
407 Int_t planEnd = kNplan;
9d0b222b 408 if (fTRD->GetSensPlane() >= 0) {
da581aea 409 planBeg = fTRD->GetSensPlane();
f7336fa3 410 planEnd = planBeg + 1;
411 }
412 Int_t sectBeg = 0;
413 Int_t sectEnd = kNsect;
f7336fa3 414
8230f242 415 Int_t countHits = 0;
da581aea 416
f7336fa3 417 // Loop through all the chambers
418 for (Int_t iCham = chamBeg; iCham < chamEnd; iCham++) {
419 for (Int_t iPlan = planBeg; iPlan < planEnd; iPlan++) {
420 for (Int_t iSect = sectBeg; iSect < sectEnd; iSect++) {
421
9d0b222b 422 if (fTRD->GetSensSector() >= 0) {
423 Int_t sens1 = fTRD->GetSensSector();
424 Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
425 sens2 -= ((Int_t) (sens2 / kNsect)) * kNsect;
dd9a6ee3 426 if (sens1 < sens2) {
9d0b222b 427 if ((iSect < sens1) || (iSect >= sens2)) continue;
dd9a6ee3 428 }
429 else {
9d0b222b 430 if ((iSect < sens1) && (iSect >= sens2)) continue;
dd9a6ee3 431 }
9d0b222b 432 }
433
f7336fa3 434 Int_t nDigits = 0;
435
f7336fa3 436 printf("AliTRDdigitizer::MakeDigits -- ");
437 printf("Digitizing chamber %d, plane %d, sector %d.\n"
438 ,iCham,iPlan,iSect);
439
da581aea 440 Int_t iDet = fGeo->GetDetector(iPlan,iCham,iSect);
441 Int_t nRowMax = fGeo->GetRowMax(iPlan,iCham,iSect);
442 Int_t nColMax = fGeo->GetColMax(iPlan);
443 Int_t nTimeMax = fGeo->GetTimeMax();
444 Float_t row0 = fGeo->GetRow0(iPlan,iCham,iSect);
445 Float_t col0 = fGeo->GetCol0(iPlan);
446 Float_t time0 = fGeo->GetTime0(iPlan);
447 Float_t rowPadSize = fGeo->GetRowPadSize();
448 Float_t colPadSize = fGeo->GetColPadSize();
449 Float_t timeBinSize = fGeo->GetTimeBinSize();
450
451 // Adjust the size of the detector arrays
8230f242 452 signal->Allocate(nRowMax,nColMax,nTimeMax);
da581aea 453 for (iDict = 0; iDict < kNDict; iDict++) {
8230f242 454 tracks[iDict]->Allocate(nRowMax,nColMax,nTimeMax);
da581aea 455 }
f7336fa3 456
457 // Loop through all entries in the tree
458 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
459
460 gAlice->ResetHits();
8230f242 461 nBytes += hitTree->GetEvent(iTrack);
f7336fa3 462
463 // Get the number of hits in the TRD created by this particle
da581aea 464 Int_t nHit = fTRD->Hits()->GetEntriesFast();
f7336fa3 465
466 // Loop through the TRD hits
467 for (Int_t iHit = 0; iHit < nHit; iHit++) {
468
8230f242 469 countHits++;
f7336fa3 470
8230f242 471 AliTRDhit *hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit);
f7336fa3 472 Float_t pos[3];
94de3818 473 pos[0] = hit->X();
474 pos[1] = hit->Y();
475 pos[2] = hit->Z();
8230f242 476 Float_t q = hit->GetCharge();
94de3818 477 Int_t track = hit->Track();
8230f242 478 Int_t detector = hit->GetDetector();
da581aea 479 Int_t plane = fGeo->GetPlane(detector);
480 Int_t sector = fGeo->GetSector(detector);
481 Int_t chamber = fGeo->GetChamber(detector);
f7336fa3 482
483 if ((sector != iSect) ||
484 (plane != iPlan) ||
485 (chamber != iCham))
486 continue;
487
488 // Rotate the sectors on top of each other
489 Float_t rot[3];
da581aea 490 fGeo->Rotate(detector,pos,rot);
f7336fa3 491
492 // The hit position in pad coordinates (center pad)
493 // The pad row (z-direction)
494 Int_t rowH = (Int_t) ((rot[2] - row0) / rowPadSize);
495 // The pad column (rphi-direction)
496 Int_t colH = (Int_t) ((rot[1] - col0) / colPadSize);
497 // The time bucket
498 Int_t timeH = (Int_t) ((rot[0] - time0) / timeBinSize);
499
500 // Array to sum up the signal in a box surrounding the
501 // hit postition
8230f242 502 const Int_t kTimeBox = 7;
503 const Int_t kColBox = 9;
504 const Int_t kRowBox = 7;
505 Float_t signalSum[kRowBox][kColBox][kTimeBox];
506 for (iRow = 0; iRow < kRowBox; iRow++ ) {
507 for (iCol = 0; iCol < kColBox; iCol++ ) {
508 for (iTime = 0; iTime < kTimeBox; iTime++) {
f7336fa3 509 signalSum[iRow][iCol][iTime] = 0;
510 }
511 }
512 }
513
514 // Loop over all electrons of this hit
515 Int_t nEl = (Int_t) q;
516 for (Int_t iEl = 0; iEl < nEl; iEl++) {
517
518 // The driftlength
519 Float_t driftlength = rot[0] - time0;
520 if ((driftlength < 0) ||
521 (driftlength > kDrThick)) break;
522 Float_t driftlengthL = driftlength;
523 if (fExBOn) driftlengthL /= TMath::Sqrt(fLorentzFactor);
524 Float_t xyz[3];
525 xyz[0] = rot[0];
526 xyz[1] = rot[1];
527 xyz[2] = rot[2];
528
529 // Electron attachment
530 if (fElAttachOn) {
531 if (gRandom->Rndm() < (driftlengthL * fElAttachProp / 100.)) continue;
532 }
533
534 // Apply the diffusion smearing
535 if (fDiffusionOn) {
536 if (!(Diffusion(driftlengthL,xyz))) continue;
537 }
538
539 // Apply E x B effects
540 if (fExBOn) {
541 if (!(ExB(driftlength,xyz))) continue;
542 }
543
544 // The electron position and the distance to the hit position
545 // in pad units
546 // The pad row (z-direction)
547 Int_t rowE = (Int_t) ((xyz[2] - row0) / rowPadSize);
548 Int_t rowD = rowH - rowE;
549 // The pad column (rphi-direction)
550 Int_t colE = (Int_t) ((xyz[1] - col0) / colPadSize);
551 Int_t colD = colH - colE;
552 // The time bucket
553 Int_t timeE = (Int_t) ((xyz[0] - time0) / timeBinSize);
554 Int_t timeD = timeH - timeE;
555
556 // Apply the gas gain including fluctuations
eda4336d 557 Float_t ggRndm = 0.0;
558 do {
559 ggRndm = gRandom->Rndm();
560 } while (ggRndm <= 0);
561 Int_t signal = (Int_t) (-fGasGain * TMath::Log(ggRndm));
f7336fa3 562
563 // The distance of the electron to the center of the pad
564 // in units of pad width
565 Float_t dist = (xyz[1] - col0 - (colE + 0.5) * colPadSize)
566 / colPadSize;
567
568 // Sum up the signal in the different pixels
569 // and apply the pad response
8230f242 570 Int_t rowIdx = rowD + (Int_t) ( kRowBox / 2);
571 Int_t colIdx = colD + (Int_t) ( kColBox / 2);
572 Int_t timeIdx = timeD + (Int_t) (kTimeBox / 2);
f7336fa3 573
8230f242 574 if (( rowIdx < 0) || ( rowIdx > kRowBox)) {
f7336fa3 575 printf("AliTRDdigitizer::MakeDigits -- ");
8230f242 576 printf("Boundary error. rowIdx = %d (%d)\n", rowIdx, kRowBox);
f7336fa3 577 continue;
578 }
8230f242 579 if (( colIdx < 0) || ( colIdx > kColBox)) {
f7336fa3 580 printf("AliTRDdigitizer::MakeDigits -- ");
8230f242 581 printf("Boundary error. colIdx = %d (%d)\n", colIdx, kColBox);
f7336fa3 582 continue;
583 }
8230f242 584 if ((timeIdx < 0) || (timeIdx > kTimeBox)) {
f7336fa3 585 printf("AliTRDdigitizer::MakeDigits -- ");
8230f242 586 printf("Boundary error. timeIdx = %d (%d)\n",timeIdx,kTimeBox);
f7336fa3 587 continue;
588 }
da581aea 589 signalSum[rowIdx][colIdx-1][timeIdx] += fPRF->Eval(dist-1.0,0,0) * signal;
590 signalSum[rowIdx][colIdx ][timeIdx] += fPRF->Eval(dist ,0,0) * signal;
591 signalSum[rowIdx][colIdx+1][timeIdx] += fPRF->Eval(dist+1.0,0,0) * signal;
f7336fa3 592
593 }
594
595 // Add the padcluster to the detector matrix
8230f242 596 for (iRow = 0; iRow < kRowBox; iRow++ ) {
597 for (iCol = 0; iCol < kColBox; iCol++ ) {
598 for (iTime = 0; iTime < kTimeBox; iTime++) {
f7336fa3 599
8230f242 600 Int_t rowB = rowH + iRow - (Int_t) ( kRowBox / 2);
601 Int_t colB = colH + iCol - (Int_t) ( kColBox / 2);
602 Int_t timeB = timeH + iTime - (Int_t) (kTimeBox / 2);
f7336fa3 603 Float_t signalB = signalSum[iRow][iCol][iTime];
da581aea 604 if (( rowB < 0) || ( rowB >= nRowMax)) continue;
605 if (( colB < 0) || ( colB >= nColMax)) continue;
606 if ((timeB < 0) || (timeB >= nTimeMax)) continue;
f7336fa3 607 if (signalB > 0.0) {
da581aea 608
609 // Add the signal sum
8230f242 610 signalB += signal->GetData(rowB,colB,timeB);
611 signal->SetData(rowB,colB,timeB,signalB);
da581aea 612 // Store the track index in the dictionary
613 // Note: We store index+1 in order to allow the array to be compressed
614 for (iDict = 0; iDict < kNDict; iDict++) {
8230f242 615 Int_t oldTrack = tracks[iDict]->GetData(rowB,colB,timeB);
da581aea 616 if (oldTrack == track+1) break;
617 if (oldTrack == -1) break;
618 if (oldTrack == 0) {
8230f242 619 tracks[iDict]->SetData(rowB,colB,timeB,track+1);
da581aea 620 break;
621 }
622 }
623 if (iDict == kNDict) {
f7336fa3 624 printf("AliTRDdigitizer::MakeDigits -- ");
da581aea 625 printf("More than three tracks for one digit!\n");
626 }
627 }
f7336fa3 628
629 }
630 }
631 }
632
633 }
634
635 }
636
637 // Add a container for the digits of this detector
8230f242 638 digits = fDigits->GetDigits(iDet);
f7336fa3 639 // Allocate memory space for the digits buffer
8230f242 640 digits->Allocate(nRowMax,nColMax,nTimeMax);
f7336fa3 641
da581aea 642 // Do the same for the dictionary arrays
28329a48 643 for (iDict = 0; iDict < kNDict; iDict++) {
8230f242 644 dictionary[iDict] = fDigits->GetDictionary(iDet,iDict);
645 dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeMax);
f7336fa3 646 }
647
da581aea 648 // Create the digits for this chamber
f7336fa3 649 for (iRow = 0; iRow < nRowMax; iRow++ ) {
650 for (iCol = 0; iCol < nColMax; iCol++ ) {
651 for (iTime = 0; iTime < nTimeMax; iTime++) {
652
8230f242 653 Float_t signalAmp = signal->GetData(iRow,iCol,iTime);
f7336fa3 654
655 // Add the noise
1befd3b2 656 signalAmp = TMath::Max((Float_t) gRandom->Gaus(signalAmp,fNoise)
657 ,(Float_t) 0.0);
f7336fa3 658 // Convert to fC
8230f242 659 signalAmp *= kEl2fC;
f7336fa3 660 // Convert to mV
661 signalAmp *= fChipGain;
662 // Convert to ADC counts
663 Int_t adc = (Int_t) (signalAmp * (fADCoutRange / fADCinRange));
664
da581aea 665 if (adc > fADCthreshold) {
f7336fa3 666
da581aea 667 nDigits++;
668
669 // Store the amplitude of the digit
8230f242 670 digits->SetData(iRow,iCol,iTime,adc);
da581aea 671
672 // Store the track index in the dictionary
673 // Note: We store index+1 in order to allow the array to be compressed
674 for (iDict = 0; iDict < kNDict; iDict++) {
8230f242 675 dictionary[iDict]->SetData(iRow,iCol,iTime
676 ,tracks[iDict]->GetData(iRow,iCol,iTime));
da581aea 677 }
f7336fa3 678
da581aea 679 }
f7336fa3 680
681 }
682 }
683 }
684
685 // Compress the arrays
8230f242 686 digits->Compress(1,0);
28329a48 687 for (iDict = 0; iDict < kNDict; iDict++) {
8230f242 688 dictionary[iDict]->Compress(1,0);
f7336fa3 689 }
690
8230f242 691 totalSizeDigits += digits->GetSize();
692 totalSizeDict0 += dictionary[0]->GetSize();
693 totalSizeDict1 += dictionary[1]->GetSize();
694 totalSizeDict2 += dictionary[2]->GetSize();
f7336fa3 695
696 printf("AliTRDdigitizer::MakeDigits -- ");
697 printf("Number of digits found: %d.\n",nDigits);
da581aea 698
699 // Reset the arrays
8230f242 700 signal->Reset();
da581aea 701 for (iDict = 0; iDict < kNDict; iDict++) {
8230f242 702 tracks[iDict]->Reset();
da581aea 703 }
f7336fa3 704
705 }
706 }
707 }
708
709 printf("AliTRDdigitizer::MakeDigits -- ");
8230f242 710 printf("Total number of analyzed hits = %d\n",countHits);
da581aea 711
712 printf("AliTRDdigitizer::MakeDigits -- ");
f7336fa3 713 printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
714 ,totalSizeDict0
715 ,totalSizeDict1
716 ,totalSizeDict2);
717
718 return kTRUE;
719
720}
721
722//_____________________________________________________________________________
f7336fa3 723Bool_t AliTRDdigitizer::WriteDigits()
724{
725 //
726 // Writes out the TRD-digits and the dictionaries
727 //
728
729 // Create the branches
730 if (!(gAlice->TreeD()->GetBranch("TRDdigits"))) {
da581aea 731 if (!fDigits->MakeBranch()) return kFALSE;
f7336fa3 732 }
733
da581aea 734 // Store the digits and the dictionary in the tree
735 fDigits->WriteDigits();
f7336fa3 736
737 // Write the new tree into the input file (use overwrite option)
738 Char_t treeName[7];
739 sprintf(treeName,"TreeD%d",fEvent);
740 printf("AliTRDdigitizer::WriteDigits -- ");
741 printf("Write the digits tree %s for event %d.\n"
742 ,treeName,fEvent);
743 gAlice->TreeD()->Write(treeName,2);
744
745 return kTRUE;
746
747}