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