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