8caaf1faac6ffc0608d0a364d2e45c90aef2b27e
[u/mrichter/AliRoot.git] / TRD / AliTRDdigitizer.cxx
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 Revision 1.4  2000/05/08 15:53:45  cblume
19 Resolved merge conflict
20
21 Revision 1.3  2000/04/28 14:49:27  cblume
22 Only one declaration of iDict in MakeDigits()
23
24 Revision 1.1.4.1  2000/05/08 14:42:04  cblume
25 Introduced AliTRDdigitsManager
26
27 Revision 1.1  2000/02/28 19:00:13  cblume
28 Add new TRD classes
29
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"
59 #include "AliTRDdataArrayI.h"
60 #include "AliTRDdataArrayF.h"
61 #include "AliTRDdigitsManager.h"
62
63 ClassImp(AliTRDdigitizer)
64
65 //_____________________________________________________________________________
66 AliTRDdigitizer::AliTRDdigitizer():TNamed()
67 {
68   //
69   // AliTRDdigitizer default constructor
70   //
71
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;
92
93 }
94
95 //_____________________________________________________________________________
96 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
97                 :TNamed(name,title)
98 {
99   //
100   // AliTRDdigitizer default constructor
101   //
102
103   fInputFile     = NULL;
104   fDigits        = NULL;
105   fTRD           = NULL;
106   fGeo           = NULL;
107
108   fEvent         = 0;
109
110   Init();
111
112 }
113
114 //_____________________________________________________________________________
115 AliTRDdigitizer::~AliTRDdigitizer()
116 {
117
118   if (fInputFile) {
119     fInputFile->Close();
120     delete fInputFile;
121   }
122
123   if (fDigits) {
124     delete fDigits;
125   }
126
127   if (fPRF) delete fPRF;
128
129 }
130
131 //_____________________________________________________________________________
132 Int_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 //_____________________________________________________________________________
149 Int_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 //_____________________________________________________________________________
164 void 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
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
198 }
199
200 //_____________________________________________________________________________
201 Bool_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
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   }
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 //_____________________________________________________________________________
245 Bool_t AliTRDdigitizer::MakeDigits()
246 {
247   //
248   // Loops through the TRD-hits and creates the digits.
249   //
250
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;
262   Int_t   iDict;
263
264   Int_t   totalSizeDigits = 0;
265   Int_t   totalSizeDict0  = 0;
266   Int_t   totalSizeDict1  = 0;
267   Int_t   totalSizeDict2  = 0;
268
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   }
287
288   // Get the pointer to the hit tree
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");
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
306   Int_t chamBeg = 0;
307   Int_t chamEnd = kNcham;
308   if (fTRD->GetSensChamber() >= 0) {
309     chamBeg = fTRD->GetSensChamber();
310     chamEnd = chamBeg + 1;
311   }
312   Int_t planBeg = 0;
313   Int_t planEnd = kNplan;
314   if (fTRD->GetSensPlane()   >= 0) {
315     planBeg = fTRD->GetSensPlane();
316     planEnd = planBeg + 1;
317   }
318   Int_t sectBeg = 0;
319   Int_t sectEnd = kNsect;
320   if (fTRD->GetSensSector()  >= 0) {
321     sectBeg = fTRD->GetSensSector();
322     sectEnd = sectBeg + 1;
323   }
324
325   Int_t count_hits = 0;
326
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
334         printf("AliTRDdigitizer::MakeDigits -- ");
335         printf("Digitizing chamber %d, plane %d, sector %d.\n"
336               ,iCham,iPlan,iSect);
337
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         }
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
362           Int_t nHit = fTRD->Hits()->GetEntriesFast();
363
364           // Loop through the TRD hits  
365           for (Int_t iHit = 0; iHit < nHit; iHit++) {
366
367             count_hits++;
368
369             AliTRDhit *Hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit);
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;
377             Int_t   plane    = fGeo->GetPlane(detector);
378             Int_t   sector   = fGeo->GetSector(detector);
379             Int_t   chamber  = fGeo->GetChamber(detector);
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];
388             fGeo->Rotate(detector,pos,rot);
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               }
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;
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);
497                   Float_t signalB = signalSum[iRow][iCol][iTime];
498                   if (( rowB < 0) || ( rowB >=  nRowMax)) continue;
499                   if (( colB < 0) || ( colB >=  nColMax)) continue;
500                   if ((timeB < 0) || (timeB >= nTimeMax)) continue;
501                   if (signalB > 0.0) {
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) {
518                       printf("AliTRDdigitizer::MakeDigits -- ");
519                       printf("More than three tracks for one digit!\n");
520                     }
521                   }
522
523                 }
524               }
525             }
526
527           }
528
529         }
530
531         // Add a container for the digits of this detector
532         Digits = fDigits->GetDigits(iDet);        
533         // Allocate memory space for the digits buffer
534         Digits->Allocate(nRowMax,nColMax,nTimeMax);
535
536         // Do the same for the dictionary arrays
537         for (iDict = 0; iDict < kNDict; iDict++) {
538           Dictionary[iDict] = fDigits->GetDictionary(iDet,iDict);
539           Dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeMax);
540         }
541
542         // Create the digits for this chamber
543         for (iRow  = 0; iRow  <  nRowMax; iRow++ ) {
544           for (iCol  = 0; iCol  <  nColMax; iCol++ ) {
545             for (iTime = 0; iTime < nTimeMax; iTime++) {         
546
547               Float_t signalAmp = Signal->GetData(iRow,iCol,iTime);
548
549               // Add the noise
550               signalAmp  = TMath::Max((Float_t) gRandom->Gaus(signalAmp,fNoise)
551                                      ,(Float_t) 0.0);
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
559               if (adc > fADCthreshold) {
560
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                 }
572
573               }
574
575             }
576           }
577         }
578
579         // Compress the arrays
580         Digits->Compress(1,0);
581         for (iDict = 0; iDict < kNDict; iDict++) {
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);
592  
593         // Reset the arrays
594         Signal->Reset();
595         for (iDict = 0; iDict < kNDict; iDict++) {
596           Tracks[iDict]->Reset();
597         }
598
599       }
600     }
601   }
602
603   printf("AliTRDdigitizer::MakeDigits -- ");
604   printf("Total number of analyzed hits = %d\n",count_hits);
605
606   printf("AliTRDdigitizer::MakeDigits -- ");
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 //_____________________________________________________________________________
617 Bool_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"))) { 
625     if (!fDigits->MakeBranch()) return kFALSE;
626   }
627
628   // Store the digits and the dictionary in the tree
629   fDigits->WriteDigits();
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 }