Update TRD code from C.Blume
[u/mrichter/AliRoot.git] / TRD / AliTRDv2.cxx
1 ///////////////////////////////////////////////////////////////////////////////
2 //                                                                           //
3 //  Transition Radiation Detector version 2 -- detailed simulation           //
4 //                                                                           //
5 //Begin_Html
6 /*
7 <img src="picts/AliTRDv2Class.gif">
8 */
9 //End_Html
10 //                                                                           //
11 //                                                                           //
12 ///////////////////////////////////////////////////////////////////////////////
13
14 #include <TMath.h>
15 #include <TVector.h>
16 #include <TRandom.h>
17
18 #include "AliTRDv2.h"
19 #include "AliTRDmatrix.h"
20 #include "AliRun.h"
21 #include "AliMC.h"
22 #include "AliConst.h"
23
24 ClassImp(AliTRDv2)
25
26 //_____________________________________________________________________________
27 AliTRDv2::AliTRDv2(const char *name, const char *title) 
28          :AliTRD(name, title) 
29 {
30   //
31   // Standard constructor for Transition Radiation Detector version 2
32   //
33
34   fIdSens       = 0;
35
36   fIdSpace1     = 0;
37   fIdSpace2     = 0;
38   fIdSpace3     = 0;
39
40   fIdChamber1   = 0;
41   fIdChamber2   = 0;
42   fIdChamber3   = 0;
43
44   fSensSelect   = 0;
45   fSensPlane    = 0;
46   fSensChamber  = 0;
47   fSensSector   = 0;
48
49   for (Int_t iplan = 0; iplan < kNplan; iplan++) {
50     for (Int_t icham = 0; icham < kNcham; icham++) {
51       fRowMax[iplan][icham] = 0;
52     }
53     fColMax[iplan] = 0;
54   }
55   fTimeMax      = 0;
56
57   fRowPadSize   = 0;
58   fColPadSize   = 0;
59   fTimeBinSize  = 0;
60
61   fGasGain      = 0;
62   fNoise        = 0;
63   fChipGain     = 0;
64   fADCoutRange  = 0;
65   fADCinRange   = 0;
66   fADCthreshold = 0;
67
68   fDiffusionT   = 0;
69   fDiffusionL   = 0;
70
71   fDeltaE       = NULL;
72
73   SetBufferSize(128000);
74
75 }
76
77 //_____________________________________________________________________________
78 AliTRDv2::~AliTRDv2()
79 {
80
81   if (fDeltaE) delete fDeltaE;
82
83 }
84  
85 //_____________________________________________________________________________
86 void AliTRDv2::CreateGeometry()
87 {
88   //
89   // Create the GEANT geometry for the Transition Radiation Detector - Version 2
90   // This version covers the full azimuth. 
91   //
92   // Author:  Christoph Blume (C.Blume@gsi.de) 20/07/99 
93   //
94
95   Float_t xpos, ypos, zpos;
96
97   // Check that FRAME is there otherwise we have no place where to put the TRD
98   AliModule* FRAME = gAlice->GetModule("FRAME");
99   if (!FRAME) return;
100
101   // Define the chambers
102   AliTRD::CreateGeometry();
103
104   // Position the the TRD-sectors in all TRD-volumes in the spaceframe
105   xpos     = 0.;
106   ypos     = 0.;
107   zpos     = 0.;
108   gMC->Gspos("TRD ",1,"BTR1",xpos,ypos,zpos,0,"ONLY");
109   gMC->Gspos("TRD ",2,"BTR2",xpos,ypos,zpos,0,"ONLY");
110   gMC->Gspos("TRD ",3,"BTR3",xpos,ypos,zpos,0,"ONLY");
111
112 }
113
114 //_____________________________________________________________________________
115 void AliTRDv2::CreateMaterials()
116 {
117   //
118   // Create materials for the Transition Radiation Detector version 2
119   //
120
121   AliTRD::CreateMaterials();
122
123 }
124
125 //_____________________________________________________________________________
126 void AliTRDv2::Diffusion(Float_t driftlength, Float_t *xyz)
127 {
128   //
129   // Applies the diffusion smearing to the position of a single electron
130   //
131
132   if ((driftlength >        0) && 
133       (driftlength < kDrThick)) {
134     Float_t driftSqrt = TMath::Sqrt(driftlength);
135     Float_t sigmaT = driftSqrt * fDiffusionT;
136     Float_t sigmaL = driftSqrt * fDiffusionL;
137     xyz[0] = gRandom->Gaus(xyz[0], sigmaL);
138     xyz[1] = gRandom->Gaus(xyz[1], sigmaT);
139     xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
140   }
141   else {
142     xyz[0] = 0.0;
143     xyz[1] = 0.0;
144     xyz[2] = 0.0;
145   }
146
147 }
148
149 //_____________________________________________________________________________
150 void AliTRDv2::Hits2Digits()
151 {
152   //
153   // Creates TRD digits from hits. This procedure includes the following:
154   //      - Diffusion
155   //      - Gas gain including fluctuations
156   //      - Pad-response (simple Gaussian approximation)
157   //      - Electronics noise
158   //      - Electronics gain
159   //      - Digitization
160   //      - ADC threshold
161   // The corresponding parameter can be adjusted via the various Set-functions.
162   // If these parameters are not explicitly set, default values are used (see
163   // Init-function).
164   // To produce digits from a galice.root file with TRD-hits use the
165   // digitsCreate.C macro.
166   //
167
168   printf(" Start creating digits\n");
169
170   ///////////////////////////////////////////////////////////////
171   // Parameter 
172   ///////////////////////////////////////////////////////////////
173
174   // Converts number of electrons to fC
175   const Float_t el2fC = 1.602E-19 * 1.0E15; 
176
177   ///////////////////////////////////////////////////////////////
178
179   Int_t nBytes = 0;
180
181   AliTRDhit *TRDhit;
182
183   // Position of pad 0,0,0 
184   // 
185   // chambers seen from the top:
186   //     +----------------------------+
187   //     |                            |
188   //     |                            |     ^
189   //     |                            | rphi|
190   //     |                            |     |
191   //     |0                           |     | 
192   //     +----------------------------+     +------>
193   //                                             z 
194   // chambers seen from the side:           ^
195   //     +----------------------------+ time|
196   //     |                            |     |
197   //     |0                           |     |
198   //     +----------------------------+     +------>
199   //                                             z
200   //                                             
201   // The pad row (z-direction)
202   Float_t row0[kNplan][kNcham];
203   for (Int_t iplan = 0; iplan < kNplan; iplan++) {
204     row0[iplan][0] = -fClengthI[iplan]/2. - fClengthM[iplan] - fClengthO[iplan] 
205                    + kCcthick; 
206     row0[iplan][1] = -fClengthI[iplan]/2. - fClengthM[iplan]                    
207                    + kCcthick;
208     row0[iplan][2] = -fClengthI[iplan]/2.                                       
209                    + kCcthick;
210     row0[iplan][3] =  fClengthI[iplan]/2.                                       
211                    + kCcthick; 
212     row0[iplan][4] =  fClengthI[iplan]/2. + fClengthM[iplan]                    
213                    + kCcthick; 
214   }
215   // The pad column (rphi-direction)  
216   Float_t col0[kNplan];
217   for (Int_t iplan = 0; iplan < kNplan; iplan++) {
218     col0[iplan]    = -fCwidth[iplan]/2. + kCcthick;
219   }
220   // The time bucket
221   Float_t time0[kNplan];
222   for (Int_t iplan = 0; iplan < kNplan; iplan++) {
223     time0[iplan]   = kRmin + kCcframe/2. + kDrZpos - 0.5 * kDrThick
224                            + iplan * (kCheight + kCspace);
225   } 
226
227   // Get the pointer to the hit tree
228   TTree *HitTree    = gAlice->TreeH();
229   // Get the pointer to the digits tree
230   TTree *DigitsTree = gAlice->TreeD();
231
232   // Get the number of entries in the hit tree
233   // (Number of primary particles creating a hit somewhere)
234   Int_t nTrack = (Int_t) HitTree->GetEntries();
235
236   Int_t chamBeg = 0;
237   Int_t chamEnd = kNcham;
238   if (fSensChamber) chamEnd = chamBeg = fSensChamber;
239   Int_t planBeg = 0;
240   Int_t planEnd = kNplan;
241   if (fSensPlane)   planEnd = planBeg = fSensPlane;
242   Int_t sectBeg = 0;
243   Int_t sectEnd = kNsect;
244   if (fSensSector)  sectEnd = sectBeg = fSensSector;
245
246   // Loop through all the chambers
247   for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
248     for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
249       for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
250
251         printf(" Digitizing chamber %d, plane %d, sector %d\n"
252               ,icham+1,iplan+1,isect+1);
253
254         // Create a detector matrix to keep the signal and track numbers
255         AliTRDmatrix *matrix = new AliTRDmatrix(fRowMax[iplan][icham]
256                                                ,fColMax[iplan]
257                                                ,fTimeMax
258                                                ,isect+1,icham+1,iplan+1);
259
260         // Loop through all entries in the tree
261         for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
262
263           gAlice->ResetHits();
264           nBytes += HitTree->GetEvent(iTrack);
265
266           // Get the number of hits in the TRD created by this particle
267           Int_t nHit = fHits->GetEntriesFast();
268
269           // Loop through the TRD hits  
270           for (Int_t iHit = 0; iHit < nHit; iHit++) {
271
272             if (!(TRDhit = (AliTRDhit *) fHits->UncheckedAt(iHit))) 
273               continue;
274
275             Float_t x       = TRDhit->fX;
276             Float_t y       = TRDhit->fY;
277             Float_t z       = TRDhit->fZ;
278             Float_t q       = TRDhit->fQ;
279             Int_t   track   = TRDhit->fTrack;
280             Int_t   plane   = TRDhit->fPlane;
281             Int_t   sector  = TRDhit->fSector;
282             Int_t   chamber = TRDhit->fChamber;        
283
284             if ((sector  != isect+1) ||
285                 (plane   != iplan+1) ||
286                 (chamber != icham+1)) 
287               continue;
288
289             // Rotate the sectors on top of each other
290             Float_t phi  = 2.0 * kPI /  (Float_t) kNsect 
291                                * ((Float_t) sector - 0.5);
292             Float_t xRot = -x * TMath::Cos(phi) + y * TMath::Sin(phi);
293             Float_t yRot =  x * TMath::Sin(phi) + y * TMath::Cos(phi);
294             Float_t zRot =  z;
295
296             // The hit position in pad coordinates (center pad)
297             // The pad row (z-direction)
298             Int_t rowH  = (Int_t) ((zRot -  row0[iplan][icham]) / fRowPadSize);
299             // The pad column (rphi-direction)  
300             Int_t colH  = (Int_t) ((yRot -  col0[iplan]       ) / fColPadSize);
301             // The time bucket
302             Int_t timeH = (Int_t) ((xRot - time0[iplan]       ) / fTimeBinSize);
303
304             // Array to sum up the signal in a box surrounding the
305             // hit postition
306             const Int_t timeBox = 5;
307             const Int_t  colBox = 7;
308             const Int_t  rowBox = 5;
309             Float_t signalSum[rowBox][colBox][timeBox];
310             for (Int_t iRow  = 0;  iRow <  rowBox; iRow++ ) {
311               for (Int_t iCol  = 0;  iCol <  colBox; iCol++ ) {
312                 for (Int_t iTime = 0; iTime < timeBox; iTime++) {
313                   signalSum[iRow][iCol][iTime] = 0;
314                 }
315               }
316             }
317
318             // Loop over all electrons of this hit
319             Int_t nEl = (Int_t) q;
320             for (Int_t iEl = 0; iEl < nEl; iEl++) {
321
322               // Apply the diffusion smearing
323               Float_t driftlength = xRot - time0[iplan];
324               Float_t xyz[3];
325               xyz[0] = xRot;
326               xyz[1] = yRot;
327               xyz[2] = zRot;
328               Diffusion(driftlength,xyz);
329
330               // At this point absorption effects that depend on the 
331               // driftlength could be taken into account.              
332
333               // The electron position and the distance to the hit position
334               // in pad units
335               // The pad row (z-direction)
336               Int_t  rowE = (Int_t) ((xyz[2] -  row0[iplan][icham]) / fRowPadSize);
337               Int_t  rowD =  rowH -  rowE;
338               // The pad column (rphi-direction)
339               Int_t  colE = (Int_t) ((xyz[1] -  col0[iplan]       ) / fColPadSize);
340               Int_t  colD =  colH -  colE;
341               // The time bucket
342               Int_t timeE = (Int_t) ((xyz[0] - time0[iplan]       ) / fTimeBinSize);
343               Int_t timeD = timeH - timeE;
344
345               // Apply the gas gain including fluctuations
346               Int_t signal = (Int_t) (-fGasGain * TMath::Log(gRandom->Rndm()));
347
348               // The distance of the electron to the center of the pad 
349               // in units of pad width
350               Float_t dist = (xyz[1] - col0[iplan] - (colE + 0.5) * fColPadSize) 
351                            / fColPadSize;
352
353               // Sum up the signal in the different pixels
354               // and apply the pad response
355               Int_t  rowIdx =  rowD + (Int_t) ( rowBox / 2);
356               Int_t  colIdx =  colD + (Int_t) ( colBox / 2);
357               Int_t timeIdx = timeD + (Int_t) (timeBox / 2);
358               signalSum[rowIdx][colIdx-1][timeIdx] += PadResponse(dist-1.) * signal;
359               signalSum[rowIdx][colIdx  ][timeIdx] += PadResponse(dist   ) * signal;
360               signalSum[rowIdx][colIdx+1][timeIdx] += PadResponse(dist+1.) * signal;
361
362             }
363             
364             // Add the padcluster to the detector matrix
365             for (Int_t iRow  = 0;  iRow <  rowBox; iRow++ ) {
366               for (Int_t iCol  = 0;  iCol <  colBox; iCol++ ) {
367                 for (Int_t iTime = 0; iTime < timeBox; iTime++) {
368
369                   Int_t  rowB =  rowH + iRow  - (Int_t) ( rowBox / 2); 
370                   Int_t  colB =  colH + iCol  - (Int_t) ( colBox / 2);
371                   Int_t timeB = timeH + iTime - (Int_t) (timeBox / 2);
372
373                   Float_t signalB = signalSum[iRow][iCol][iTime];
374                   if (signalB > 0.0) {
375                     matrix->AddSignal(rowB,colB,timeB,signalB);
376                     if (!(matrix->AddTrack(rowB,colB,timeB,track))) 
377                       printf("More than three tracks in a pixel!\n");
378                   }
379
380                 }
381               }
382             }
383
384           }
385
386         }
387
388         // Create the hits for this chamber
389         for (Int_t iRow  = 0; iRow  <  fRowMax[iplan][icham]; iRow++ ) {
390           for (Int_t iCol  = 0; iCol  <  fColMax[iplan]         ; iCol++ ) {
391             for (Int_t iTime = 0; iTime < fTimeMax                ; iTime++) {         
392
393               Float_t signalAmp = matrix->GetSignal(iRow,iCol,iTime);
394
395               // Add the noise
396               signalAmp  = TMath::Max(gRandom->Gaus(signalAmp,fNoise),0.0);
397               // Convert to fC
398               signalAmp *= el2fC;
399               // Convert to mV
400               signalAmp *= fChipGain;
401               // Convert to ADC counts
402               Int_t adc  = (Int_t) (signalAmp * (fADCoutRange / fADCinRange));
403
404               // Apply threshold on ADC value
405               if (adc > fADCthreshold) {
406
407                 Int_t trackSave[3];
408                 for (Int_t ii = 0; ii < 3; ii++) {
409                   trackSave[ii] = matrix->GetTrack(iRow,iCol,iTime,ii);
410                 }
411
412                 Int_t digits[7];
413                 digits[0] = matrix->GetSector();
414                 digits[1] = matrix->GetChamber();
415                 digits[2] = matrix->GetPlane();
416                 digits[3] = iRow;
417                 digits[4] = iCol;
418                 digits[5] = iTime;
419                 digits[6] = adc;
420
421                 // Add this digit to the TClonesArray
422                 AddDigit(trackSave,digits);
423
424               }
425
426             }
427           }
428         }
429
430         // Clean up
431         delete matrix;
432
433       }
434     }
435   }
436
437   // Fill the digits-tree
438   DigitsTree->Fill();
439
440 }
441
442 //_____________________________________________________________________________
443 void AliTRDv2::Init() 
444 {
445   //
446   // Initialise Transition Radiation Detector after geometry has been built.
447   // Includes the default settings of all parameter for the digitization.
448   //
449
450   AliTRD::Init();
451
452   if (fSensPlane)
453     printf("          Only plane %d is sensitive\n",fSensPlane);
454   if (fSensChamber)   
455     printf("          Only chamber %d is sensitive\n",fSensChamber);
456   if (fSensSector)
457     printf("          Only sector %d is sensitive\n",fSensSector);
458
459   for (Int_t i = 0; i < 80; i++) printf("*");
460   printf("\n");
461
462   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
463   const Float_t kPoti = 12.1;
464   // Maximum energy (50 keV);
465   const Float_t kEend = 50000.0;
466   // Ermilova distribution for the delta-ray spectrum
467   Float_t Poti = TMath::Log(kPoti);
468   Float_t Eend = TMath::Log(kEend);
469   fDeltaE  = new TF1("deltae",Ermilova,Poti,Eend,0);
470
471   // Identifier of the sensitive volume (drift region)
472   fIdSens     = gMC->VolId("UL05");
473
474   // Identifier of the TRD-spaceframe volumina
475   fIdSpace1   = gMC->VolId("B028");
476   fIdSpace2   = gMC->VolId("B029");
477   fIdSpace3   = gMC->VolId("B030");
478
479   // Identifier of the TRD-driftchambers
480   fIdChamber1 = gMC->VolId("UCIO");
481   fIdChamber2 = gMC->VolId("UCIM");
482   fIdChamber3 = gMC->VolId("UCII");
483
484   // The default pad dimensions
485   if (!(fRowPadSize))  fRowPadSize  = 4.5;
486   if (!(fColPadSize))  fColPadSize  = 1.0;
487   if (!(fTimeBinSize)) fTimeBinSize = 0.1;
488
489   // The maximum number of pads
490   for (Int_t iplan = 0; iplan < kNplan; iplan++) {
491     // Rows 
492     fRowMax[iplan][0] = 1 + TMath::Nint((fClengthO[iplan] - 2. * kCcthick) 
493                                                           / fRowPadSize - 0.5);
494     fRowMax[iplan][1] = 1 + TMath::Nint((fClengthM[iplan] - 2. * kCcthick) 
495                                                           / fRowPadSize - 0.5);
496     fRowMax[iplan][2] = 1 + TMath::Nint((fClengthI[iplan] - 2. * kCcthick) 
497                                                           / fRowPadSize - 0.5);
498     fRowMax[iplan][3] = 1 + TMath::Nint((fClengthM[iplan] - 2. * kCcthick) 
499                                                           / fRowPadSize - 0.5);
500     fRowMax[iplan][4] = 1 + TMath::Nint((fClengthO[iplan] - 2. * kCcthick) 
501                                                           / fRowPadSize - 0.5);
502     // Columns
503     fColMax[iplan]    = 1 + TMath::Nint((fCwidth[iplan]   - 2. * kCcthick) 
504                                                           / fColPadSize - 0.5);
505   }
506   // Time buckets
507   fTimeMax = 1 + TMath::Nint(kDrThick / fTimeBinSize - 0.5);
508
509   // The default parameter for the digitization
510   if (!(fGasGain))      fGasGain      = 2.0E3;
511   if (!(fNoise))        fNoise        = 3000.;
512   if (!(fChipGain))     fChipGain     = 10.;
513   if (!(fADCoutRange))  fADCoutRange  = 255.;
514   if (!(fADCinRange))   fADCinRange   = 2000.;
515   if (!(fADCthreshold)) fADCthreshold = 0;
516
517   // Transverse and longitudinal diffusion coefficients (Xe/Isobutane)
518   if (!(fDiffusionT))   fDiffusionT   = 0.060;
519   if (!(fDiffusionL))   fDiffusionL   = 0.017;
520
521 }
522
523 //_____________________________________________________________________________
524 void AliTRDv2::MakeBranch(Option_t* option)
525 {
526   //
527   // Create Tree branches for the TRD digits.
528   //
529
530   Int_t  buffersize = 4000;
531   Char_t branchname[10];
532
533   sprintf(branchname,"%s",GetName());
534
535   AliDetector::MakeBranch(option); 
536
537   Char_t *D = strstr(option,"D");
538   if (fDigits && gAlice->TreeD() && D) {
539     gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
540     printf("Making Branch %s for digits\n",branchname);
541   }
542
543 }
544
545 //_____________________________________________________________________________
546 Float_t AliTRDv2::PadResponse(Float_t x)
547 {
548   //
549   // The pad response for the chevron pads. 
550   // We use a simple Gaussian approximation which should be good
551   // enough for our purpose.
552   //
553
554   // The parameters for the response function
555   const Float_t aa  =  0.8872;
556   const Float_t bb  = -0.00573;
557   const Float_t cc  =  0.454;
558   const Float_t cc2 =  cc*cc;
559
560   Float_t pr = aa * (bb + TMath::Exp(-x*x / (2. * cc2)));
561
562   //TF1 *funPR = new TF1("funPR","[0]*([1]+exp(-x*x /(2.*[2])))",-3,3);
563   //funPR->SetParameter(0,aa );
564   //funPR->SetParameter(1,bb );
565   //funPR->SetParameter(2,cc2);
566   //
567   //Float_t pr = funPR->Eval(distance);
568   //
569   //delete funPR;
570
571   return (pr);
572
573 }
574
575 //_____________________________________________________________________________
576 void AliTRDv2::SetSensPlane(Int_t iplane)
577 {
578   //
579   // Defines the hit-sensitive plane (1-6)
580   //
581
582   if ((iplane < 0) || (iplane > 6)) {
583     printf("Wrong input value: %d\n",iplane);
584     printf("Use standard setting\n");
585     fSensPlane  = 0;
586     fSensSelect = 0;
587     return;
588   }
589
590   fSensSelect = 1;
591   fSensPlane  = iplane;
592
593 }
594
595 //_____________________________________________________________________________
596 void AliTRDv2::SetSensChamber(Int_t ichamber)
597 {
598   //
599   // Defines the hit-sensitive chamber (1-5)
600   //
601
602   if ((ichamber < 0) || (ichamber > 5)) {
603     printf("Wrong input value: %d\n",ichamber);
604     printf("Use standard setting\n");
605     fSensChamber = 0;
606     fSensSelect  = 0;
607     return;
608   }
609
610   fSensSelect  = 1;
611   fSensChamber = ichamber;
612
613 }
614
615 //_____________________________________________________________________________
616 void AliTRDv2::SetSensSector(Int_t isector)
617 {
618   //
619   // Defines the hit-sensitive sector (1-18)
620   //
621
622   if ((isector < 0) || (isector > 18)) {
623     printf("Wrong input value: %d\n",isector);
624     printf("Use standard setting\n");
625     fSensSector = 0;
626     fSensSelect = 0;
627     return;
628   }
629
630   fSensSelect = 1;
631   fSensSector = isector;
632
633 }
634
635 //_____________________________________________________________________________
636 void AliTRDv2::StepManager()
637 {
638   //
639   // Called at every step in the Transition Radiation Detector version 2.
640   // Slow simulator. Every charged track produces electron cluster as hits 
641   // along its path across the drift volume. The step size is set acording
642   // to Bethe-Bloch. The energy distribution of the delta electrons follows
643   // a spectrum taken from Ermilova et al.
644   //
645
646   Int_t    iIdSens, icSens;
647   Int_t    iIdSpace, icSpace;
648   Int_t    iIdChamber, icChamber;
649   Int_t    vol[3]; 
650   Int_t    iPid;
651
652   Int_t    secMap1[10] = {  3,  7,  8,  9, 10, 11,  2,  1, 18, 17 };
653   Int_t    secMap2[ 5] = { 16, 15, 14, 13, 12 };
654   Int_t    secMap3[ 3] = {  5,  6,  4 };
655
656   Float_t  hits[4];
657   Float_t  random[1];
658   Float_t  charge;
659   Float_t  aMass;
660
661   Double_t pTot;
662   Double_t qTot;
663   Double_t eDelta;
664   Double_t betaGamma, pp;
665
666   TLorentzVector pos, mom;
667   TClonesArray  &lhits = *fHits;
668
669   const Double_t kBig = 1.0E+12;
670
671   // Ionization energy
672   const Float_t kWion    = 22.04;
673   // Maximum energy for e+ e- g for the step-size calculation
674   const Float_t kPTotMax = 0.002;
675   // Plateau value of the energy-loss for electron in xenon
676   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
677   //const Double_t kPlateau = 1.70;
678   // the averaged value (26/3/99)
679   const Float_t kPlateau = 1.55;
680   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
681   const Float_t kPrim    = 48.0;
682   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
683   const Float_t kPoti    = 12.1;
684
685   // Set the maximum step size to a very large number for all 
686   // neutral particles and those outside the driftvolume
687   gMC->SetMaxStep(kBig); 
688
689   // Use only charged tracks 
690   if (( gMC->TrackCharge()       ) &&
691       (!gMC->IsTrackStop()       ) && 
692       (!gMC->IsTrackDisappeared())) {
693
694     // Inside a sensitive volume?
695     iIdSens = gMC->CurrentVolID(icSens);
696     if (iIdSens == fIdSens) { 
697
698       iIdSpace   = gMC->CurrentVolOffID(4,icSpace  );
699       iIdChamber = gMC->CurrentVolOffID(1,icChamber);
700
701       // Calculate the energy of the delta-electrons
702       eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
703       eDelta = TMath::Max(eDelta,0.0);
704
705       // The number of secondary electrons created
706       qTot = (Double_t) ((Int_t) (eDelta / kWion) + 1);
707
708       // The hit coordinates and charge
709       gMC->TrackPosition(pos);
710       hits[0] = pos[0];
711       hits[1] = pos[1];
712       hits[2] = pos[2];
713       hits[3] = qTot;
714
715       // The sector number
716       if      (iIdSpace == fIdSpace1) 
717         vol[0] = secMap1[icSpace-1];
718       else if (iIdSpace == fIdSpace2) 
719         vol[0] = secMap2[icSpace-1];
720       else if (iIdSpace == fIdSpace3) 
721         vol[0] = secMap3[icSpace-1];
722
723       // The chamber number 
724       //   1: outer left
725       //   2: middle left
726       //   3: inner
727       //   4: middle right
728       //   5: outer right
729       if      (iIdChamber == fIdChamber1)
730         vol[1] = (hits[2] < 0 ? 1 : 5);
731       else if (iIdChamber == fIdChamber2)       
732         vol[1] = (hits[2] < 0 ? 2 : 4);
733       else if (iIdChamber == fIdChamber3)       
734         vol[1] = 3;
735
736       // The plane number
737       vol[2] = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6;
738
739       // Check on selected volumes
740       Int_t addthishit = 1;
741       if (fSensSelect) {
742         if ((fSensPlane)   && (vol[2] != fSensPlane  )) addthishit = 0;
743         if ((fSensChamber) && (vol[1] != fSensChamber)) addthishit = 0;
744         if ((fSensSector)  && (vol[0] != fSensSector )) addthishit = 0;
745       }
746
747       // Add this hit
748       if (addthishit) {
749
750         new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
751
752         // The energy loss according to Bethe Bloch
753         gMC->TrackMomentum(mom);
754         pTot = mom.Rho();
755         iPid = gMC->TrackPid();
756         if ( (iPid >  3) ||
757             ((iPid <= 3) && (pTot < kPTotMax))) {
758           aMass     = gMC->TrackMass();
759           betaGamma = pTot / aMass;
760           pp        = kPrim * BetheBloch(betaGamma);
761           // Take charge > 1 into account
762           charge = gMC->TrackCharge();
763           if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
764         }
765         // Electrons above 20 Mev/c are at the plateau
766         else {
767           pp = kPrim * kPlateau;
768         }
769       
770         // Calculate the maximum step size for the next tracking step
771         if (pp > 0) {
772           do 
773             gMC->Rndm(random,1);
774           while ((random[0] == 1.) || (random[0] == 0.));
775           gMC->SetMaxStep( - TMath::Log(random[0]) / pp);
776         }
777
778       }
779       else {
780         // set step size to maximal value
781         gMC->SetMaxStep(kBig); 
782       }
783
784     }
785
786   }
787
788 }
789
790 //_____________________________________________________________________________
791 Double_t AliTRDv2::BetheBloch(Double_t bg) 
792 {
793   //
794   // Parametrization of the Bethe-Bloch-curve
795   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
796   //
797
798   // This parameters have been adjusted to averaged values from GEANT
799   const Double_t kP1 = 7.17960e-02;
800   const Double_t kP2 = 8.54196;
801   const Double_t kP3 = 1.38065e-06;
802   const Double_t kP4 = 5.30972;
803   const Double_t kP5 = 2.83798;
804
805   // This parameters have been adjusted to Xe-data found in:
806   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
807   //const Double_t kP1 = 0.76176E-1;
808   //const Double_t kP2 = 10.632;
809   //const Double_t kP3 = 3.17983E-6;
810   //const Double_t kP4 = 1.8631;
811   //const Double_t kP5 = 1.9479;
812
813   if (bg > 0) {
814     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
815     Double_t aa = TMath::Power(yy,kP4);
816     Double_t bb = TMath::Power((1./bg),kP5);
817              bb = TMath::Log(kP3 + bb);
818     return ((kP2 - aa - bb)*kP1 / aa);
819   }
820   else
821     return 0;
822
823 }
824
825 //_____________________________________________________________________________
826 Double_t Ermilova(Double_t *x, Double_t *)
827 {
828   //
829   // Calculates the delta-ray energy distribution according to Ermilova.
830   // Logarithmic scale !
831   //
832
833   Double_t energy;
834   Double_t dpos;
835   Double_t dnde;
836
837   Int_t    pos1, pos2;
838
839   const Int_t nV = 31;
840
841   Float_t vxe[nV] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120  
842                     , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
843                     , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
844                     , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
845                     , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
846                     , 9.4727, 9.9035,10.3735,10.5966,10.8198
847                     ,11.5129 };
848
849   Float_t vye[nV] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
850                     , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
851                     ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
852                     ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
853                     ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
854                     ,  0.04 ,  0.023,  0.015,  0.011,  0.01
855                     ,  0.004 };
856
857   energy = x[0];
858
859   // Find the position 
860   pos1 = pos2 = 0;
861   dpos = 0;
862   do {
863     dpos = energy - vxe[pos2++];
864   } 
865   while (dpos > 0);
866   pos2--; 
867   if (pos2 > nV) pos2 = nV;
868   pos1 = pos2 - 1;
869
870   // Differentiate between the sampling points
871   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
872
873   return dnde;
874
875 }