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