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