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