TRD/AliTRDv1.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.27  2001/03/30 14:40:15  cblume
  19 Update of the digitization parameter
  20
  21 Revision 1.26  2000/11/30 17:38:08  cblume
  22 Changes to get in line with new STEER and EVGEN
  23
  24 Revision 1.25  2000/11/15 14:30:16  cblume
  25 Fixed bug in calculating detector no. of extra hit
  26
  27 Revision 1.24  2000/11/10 14:58:36  cblume
  28 Introduce additional hit with amplitude 0 at the chamber borders
  29
  30 Revision 1.23  2000/11/01 14:53:21  cblume
  31 Merge with TRD-develop
  32
  33 Revision 1.17.2.5  2000/10/15 23:40:01  cblume
  34 Remove AliTRDconst
  35
  36 Revision 1.17.2.4  2000/10/06 16:49:46  cblume
  37 Made Getters const
  38
  39 Revision 1.17.2.3  2000/10/04 16:34:58  cblume
  40 Replace include files by forward declarations
  41
  42 Revision 1.17.2.2  2000/09/18 13:50:17  cblume
  43 Include TR photon generation and adapt to new AliTRDhit
  44
  45 Revision 1.22  2000/06/27 13:08:50  cblume
  46 Changed to Copy(TObject &A) to appease the HP-compiler
  47
  48 Revision 1.21  2000/06/09 11:10:07  cblume
  49 Compiler warnings and coding conventions, next round
  50
  51 Revision 1.20  2000/06/08 18:32:58  cblume
  52 Make code compliant to coding conventions
  53
  54 Revision 1.19  2000/06/07 16:27:32  cblume
  55 Try to remove compiler warnings on Sun and HP
  56
  57 Revision 1.18  2000/05/08 16:17:27  cblume
  58 Merge TRD-develop
  59
  60 Revision 1.17.2.1  2000/05/08 14:59:16  cblume
  61 Made inline function non-virtual. Bug fix in setting sensitive chamber
  62
  63 Revision 1.17  2000/02/28 19:10:26  cblume
  64 Include the new TRD classes
  65
  66 Revision 1.16.4.1  2000/02/28 18:04:35  cblume
  67 Change to new hit version, introduce geometry class, and move digitization and clustering to AliTRDdigitizer/AliTRDclusterizerV1
  68
  69 Revision 1.16  1999/11/05 22:50:28  fca
  70 Do not use Atan, removed from ROOT too
  71
  72 Revision 1.15  1999/11/02 17:20:19  fca
  73 initialise nbytes before using it
  74
  75 Revision 1.14  1999/11/02 17:15:54  fca
  76 Correct ansi scoping not accepted by HP compilers
  77
  78 Revision 1.13  1999/11/02 17:14:51  fca
  79 Correct ansi scoping not accepted by HP compilers
  80
  81 Revision 1.12  1999/11/02 16:35:56  fca
  82 New version of TRD introduced
  83
  84 Revision 1.11  1999/11/01 20:41:51  fca
  85 Added protections against using the wrong version of FRAME
  86
  87 Revision 1.10  1999/09/29 09:24:35  fca
  88 Introduction of the Copyright and cvs Log
  89
  90 */
  91
  92 ///////////////////////////////////////////////////////////////////////////////
  93 //                                                                           //
  94 //  Transition Radiation Detector version 1 -- slow simulator                //
  95 //                                                                           //
  96 //Begin_Html
  97 /*
  98 <img src="picts/AliTRDfullClass.gif">
  99 */
 100 //End_Html
 101 //                                                                           //
 102 //                                                                           //
 103 ///////////////////////////////////////////////////////////////////////////////
 104
 105 #include <stdlib.h>
 106
 107 #include <TMath.h>
 108 #include <TVector.h>
 109 #include <TRandom.h>
 110 #include <TF1.h>
 111 #include <TLorentzVector.h>
 112
 113 #include "AliRun.h"
 114 #include "AliMC.h"
 115 #include "AliConst.h"
 116
 117 #include "AliTRDv1.h"
 118 #include "AliTRDhit.h"
 119 #include "AliTRDmatrix.h"
 120 #include "AliTRDgeometry.h"
 121 #include "AliTRDsim.h"
 122
 123 ClassImp(AliTRDv1)
 124
 125 //_____________________________________________________________________________
 126 AliTRDv1::AliTRDv1():AliTRD()
 127 {
 128   //
 129   // Default constructor
 130   //
 131
 132   fIdSensDr        =  0;
 133   fIdSensAm        =  0;
 134
 135   fIdChamber1      =  0;
 136   fIdChamber2      =  0;
 137   fIdChamber3      =  0;
 138
 139   fSensSelect      =  0;
 140   fSensPlane       = -1;
 141   fSensChamber     = -1;
 142   fSensSector      = -1;
 143   fSensSectorRange =  0;
 144
 145   fDeltaE          = NULL;
 146   fTR              = NULL;
 147
 148 }
 149
 150 //_____________________________________________________________________________
 151 AliTRDv1::AliTRDv1(const char *name, const char *title)
 152          :AliTRD(name, title)
 153 {
 154   //
 155   // Standard constructor for Transition Radiation Detector version 1
 156   //
 157
 158   fIdSensDr        =  0;
 159   fIdSensAm        =  0;
 160
 161   fIdChamber1      =  0;
 162   fIdChamber2      =  0;
 163   fIdChamber3      =  0;
 164
 165   fSensSelect      =  0;
 166   fSensPlane       = -1;
 167   fSensChamber     = -1;
 168   fSensSector      = -1;
 169   fSensSectorRange =  0;
 170
 171   fDeltaE          = NULL;
 172   fTR              = NULL;
 173
 174   SetBufferSize(128000);
 175
 176 }
 177
 178 //_____________________________________________________________________________
 179 AliTRDv1::AliTRDv1(const AliTRDv1 &trd)
 180 {
 181   //
 182   // Copy constructor
 183   //
 184
 185   ((AliTRDv1 &) trd).Copy(*this);
 186
 187 }
 188
 189 //_____________________________________________________________________________
 190 AliTRDv1::~AliTRDv1()
 191 {
 192   //
 193   // AliTRDv1 destructor
 194   //
 195
 196   if (fDeltaE) delete fDeltaE;
 197   if (fTR)     delete fTR;
 198
 199 }
 200
 201 //_____________________________________________________________________________
 202 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
 203 {
 204   //
 205   // Assignment operator
 206   //
 207
 208   if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
 209   return *this;
 210
 211 }
 212
 213 //_____________________________________________________________________________
 214 void AliTRDv1::Copy(TObject &trd)
 215 {
 216   //
 217   // Copy function
 218   //
 219
 220   ((AliTRDv1 &) trd).fIdSensDr        = fIdSensDr;
 221   ((AliTRDv1 &) trd).fIdSensAm        = fIdSensAm;
 222
 223   ((AliTRDv1 &) trd).fIdChamber1      = fIdChamber1;
 224   ((AliTRDv1 &) trd).fIdChamber2      = fIdChamber2;
 225   ((AliTRDv1 &) trd).fIdChamber3      = fIdChamber3;
 226
 227   ((AliTRDv1 &) trd).fSensSelect      = fSensSelect;
 228   ((AliTRDv1 &) trd).fSensPlane       = fSensPlane;
 229   ((AliTRDv1 &) trd).fSensChamber     = fSensChamber;
 230   ((AliTRDv1 &) trd).fSensSector      = fSensSector;
 231   ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange;
 232
 233   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
 234   fTR->Copy(*((AliTRDv1 &) trd).fTR);
 235
 236 }
 237
 238 //_____________________________________________________________________________
 239 void AliTRDv1::CreateGeometry()
 240 {
 241   //
 242   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
 243   // This version covers the full azimuth.
 244   //
 245
 246   // Check that FRAME is there otherwise we have no place where to put the TRD
 247   AliModule* frame = gAlice->GetModule("FRAME");
 248   if (!frame) return;
 249
 250   // Define the chambers
 251   AliTRD::CreateGeometry();
 252
 253 }
 254
 255 //_____________________________________________________________________________
 256 void AliTRDv1::CreateMaterials()
 257 {
 258   //
 259   // Create materials for the Transition Radiation Detector version 1
 260   //
 261
 262   AliTRD::CreateMaterials();
 263
 264 }
 265
 266 //_____________________________________________________________________________
 267 void AliTRDv1::CreateTRhit(Int_t det)
 268 {
 269   //
 270   // Creates an electron cluster from a TR photon.
 271   // The photon is assumed to be created a the end of the radiator. The
 272   // distance after which it deposits its energy takes into account the
 273   // absorbtion of the entrance window and of the gas mixture in drift
 274   // volume.
 275   //
 276
 277   // PDG code electron
 278   const Int_t   kPdgElectron = 11;
 279
 280   // Ionization energy
 281   const Float_t kWion        = 22.04;
 282
 283   // Maximum number of TR photons per track
 284   const Int_t   kNTR         = 50;
 285
 286   TLorentzVector mom, pos;
 287   TClonesArray  &lhits = *fHits;
 288
 289   // Create TR at the entrance of the chamber
 290   if (gMC->IsTrackEntering()) {
 291
 292     // Create TR only for electrons
 293     Int_t iPdg = gMC->TrackPid();
 294     if (TMath::Abs(iPdg) != kPdgElectron) return;
 295
 296     Float_t eTR[kNTR];
 297     Int_t   nTR;
 298
 299     // Create TR photons
 300     gMC->TrackMomentum(mom);
 301     Float_t pTot = mom.Rho();
 302     fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
 303     if (nTR > kNTR) {
 304       printf("AliTRDv1::CreateTRhit -- ");
 305       printf("Boundary error: nTR = %d, kNTR = %d\n",nTR,kNTR);
 306       exit(1);
 307     }
 308
 309     // Loop through the TR photons
 310     for (Int_t iTR = 0; iTR < nTR; iTR++) {
 311
 312       Float_t energyMeV = eTR[iTR] * 0.001;
 313       Float_t energyeV  = eTR[iTR] * 1000.0;
 314       Float_t absLength = 0;
 315       Float_t sigma     = 0;
 316
 317       // Take the absorbtion in the entrance window into account
 318       Double_t muMy = fTR->GetMuMy(energyMeV);
 319       sigma = muMy * fFoilDensity;
 320       absLength = gRandom->Exp(sigma);
 321       if (absLength < AliTRDgeometry::MyThick()) continue;
 322
 323       // The absorbtion cross sections in the drift gas
 324       if (fGasMix == 1) {
 325         // Gas-mixture (Xe/CO2)
 326         Double_t muXe = fTR->GetMuXe(energyMeV);
 327         Double_t muCO = fTR->GetMuCO(energyMeV);
 328         sigma = (0.90 * muXe + 0.10 * muCO) * fGasDensity;
 329       }
 330       else {
 331         // Gas-mixture (Xe/Isobutane)
 332         Double_t muXe = fTR->GetMuXe(energyMeV);
 333         Double_t muBu = fTR->GetMuBu(energyMeV);
 334         sigma = (0.97 * muXe + 0.03 * muBu) * fGasDensity;
 335       }
 336
 337       // The distance after which the energy of the TR photon
 338       // is deposited.
 339       absLength = gRandom->Exp(sigma);
 340       if (absLength > AliTRDgeometry::DrThick()) continue;
 341
 342       // The position of the absorbtion
 343       Float_t posHit[3];
 344       gMC->TrackPosition(pos);
 345       posHit[0] = pos[0] + mom[0] / pTot * absLength;
 346       posHit[1] = pos[1] + mom[1] / pTot * absLength;
 347       posHit[2] = pos[2] + mom[2] / pTot * absLength;
 348
 349       // Create the charge
 350       Int_t q = ((Int_t) (energyeV / kWion));
 351
 352       // Add the hit to the array. TR photon hits are marked
 353       // by negative charge
 354       AliTRDhit *hit = new(lhits[fNhits++]) AliTRDhit(fIshunt
 355                                                      ,gAlice->CurrentTrack()
 356                                                      ,det,posHit,-q);
 357       hit->SetTRphoton();
 358
 359     }
 360
 361   }
 362
 363 }
 364
 365 //_____________________________________________________________________________
 366 void AliTRDv1::Init()
 367 {
 368   //
 369   // Initialise Transition Radiation Detector after geometry has been built.
 370   //
 371
 372   AliTRD::Init();
 373
 374   printf("          Slow simulator\n\n");
 375   if (fSensSelect) {
 376     if (fSensPlane   >= 0)
 377       printf("          Only plane %d is sensitive\n",fSensPlane);
 378     if (fSensChamber >= 0)
 379       printf("          Only chamber %d is sensitive\n",fSensChamber);
 380     if (fSensSector  >= 0) {
 381       Int_t sens1  = fSensSector;
 382       Int_t sens2  = fSensSector + fSensSectorRange;
 383             sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 384                    * AliTRDgeometry::Nsect();
 385       printf("          Only sectors %d - %d are sensitive\n",sens1,sens2-1);
 386     }
 387   }
 388   if (fTR)
 389     printf("          TR simulation on\n");
 390   else
 391     printf("          TR simulation off\n");
 392   printf("\n");
 393
 394   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 395   const Float_t kPoti = 12.1;
 396   // Maximum energy (50 keV);
 397   const Float_t kEend = 50000.0;
 398   // Ermilova distribution for the delta-ray spectrum
 399   Float_t poti = TMath::Log(kPoti);
 400   Float_t eEnd = TMath::Log(kEend);
 401   fDeltaE = new TF1("deltae",Ermilova,poti,eEnd,0);
 402
 403   // Identifier of the sensitive volume (drift region)
 404   fIdSensDr   = gMC->VolId("UL05");
 405   // Identifier of the sensitive volume (amplification region)
 406   fIdSensAm   = gMC->VolId("UL06");
 407
 408   // Identifier of the TRD-driftchambers
 409   fIdChamber1 = gMC->VolId("UCIO");
 410   fIdChamber2 = gMC->VolId("UCIM");
 411   fIdChamber3 = gMC->VolId("UCII");
 412
 413   for (Int_t i = 0; i < 80; i++) printf("*");
 414   printf("\n");
 415
 416 }
 417
 418 //_____________________________________________________________________________
 419 AliTRDsim *AliTRDv1::CreateTR()
 420 {
 421   //
 422   // Enables the simulation of TR
 423   //
 424
 425   fTR = new AliTRDsim();
 426   return fTR;
 427
 428 }
 429
 430 //_____________________________________________________________________________
 431 void AliTRDv1::SetSensPlane(Int_t iplane)
 432 {
 433   //
 434   // Defines the hit-sensitive plane (0-5)
 435   //
 436
 437   if ((iplane < 0) || (iplane > 5)) {
 438     printf("Wrong input value: %d\n",iplane);
 439     printf("Use standard setting\n");
 440     fSensPlane  = -1;
 441     fSensSelect =  0;
 442     return;
 443   }
 444
 445   fSensSelect = 1;
 446   fSensPlane  = iplane;
 447
 448 }
 449
 450 //_____________________________________________________________________________
 451 void AliTRDv1::SetSensChamber(Int_t ichamber)
 452 {
 453   //
 454   // Defines the hit-sensitive chamber (0-4)
 455   //
 456
 457   if ((ichamber < 0) || (ichamber > 4)) {
 458     printf("Wrong input value: %d\n",ichamber);
 459     printf("Use standard setting\n");
 460     fSensChamber = -1;
 461     fSensSelect  =  0;
 462     return;
 463   }
 464
 465   fSensSelect  = 1;
 466   fSensChamber = ichamber;
 467
 468 }
 469
 470 //_____________________________________________________________________________
 471 void AliTRDv1::SetSensSector(Int_t isector)
 472 {
 473   //
 474   // Defines the hit-sensitive sector (0-17)
 475   //
 476
 477   SetSensSector(isector,1);
 478
 479 }
 480
 481 //_____________________________________________________________________________
 482 void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
 483 {
 484   //
 485   // Defines a range of hit-sensitive sectors. The range is defined by
 486   // <isector> (0-17) as the starting point and <nsector> as the number
 487   // of sectors to be included.
 488   //
 489
 490   if ((isector < 0) || (isector > 17)) {
 491     printf("Wrong input value <isector>: %d\n",isector);
 492     printf("Use standard setting\n");
 493     fSensSector      = -1;
 494     fSensSectorRange =  0;
 495     fSensSelect      =  0;
 496     return;
 497   }
 498
 499   if ((nsector < 1) || (nsector > 18)) {
 500     printf("Wrong input value <nsector>: %d\n",nsector);
 501     printf("Use standard setting\n");
 502     fSensSector      = -1;
 503     fSensSectorRange =  0;
 504     fSensSelect      =  0;
 505     return;
 506   }
 507
 508   fSensSelect      = 1;
 509   fSensSector      = isector;
 510   fSensSectorRange = nsector;
 511
 512 }
 513
 514 //_____________________________________________________________________________
 515 void AliTRDv1::StepManager()
 516 {
 517   //
 518   // Slow simulator. Every charged track produces electron cluster as hits
 519   // along its path across the drift volume. The step size is set acording
 520   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 521   // a spectrum taken from Ermilova et al.
 522   //
 523
 524   Int_t    iIdSens, icSens;
 525   Int_t    iIdChamber, icChamber;
 526   Int_t    pla = 0;
 527   Int_t    cha = 0;
 528   Int_t    sec = 0;
 529   Int_t    det = 0;
 530   Int_t    iPdg;
 531   Int_t    qTot;
 532
 533   Float_t  hits[3];
 534   Float_t  moms[3];
 535   Float_t  random[1];
 536   Float_t  charge;
 537   Float_t  aMass;
 538
 539   Double_t pTot = 0;
 540   Double_t eDelta;
 541   Double_t betaGamma, pp;
 542   Double_t stepSize;
 543
 544   TLorentzVector pos, mom;
 545   TClonesArray  &lhits = *fHits;
 546
 547   const Double_t kBig     = 1.0E+12;
 548
 549   // Ionization energy
 550   const Float_t  kWion        = 22.04;
 551   // Maximum momentum for e+ e- g
 552   const Float_t  kPTotMaxEl   = 0.002;
 553   // Minimum energy for the step size adjustment
 554   const Float_t  kEkinMinStep = 1.0e-5;
 555   // Plateau value of the energy-loss for electron in xenon
 556   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 557   //const Double_t kPlateau = 1.70;
 558   // the averaged value (26/3/99)
 559   const Float_t  kPlateau     = 1.55;
 560   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 561   const Float_t  kPrim        = 48.0;
 562   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 563   const Float_t  kPoti        = 12.1;
 564
 565   // PDG code electron
 566   const Int_t    kPdgElectron = 11;
 567
 568   // Set the maximum step size to a very large number for all
 569   // neutral particles and those outside the driftvolume
 570   gMC->SetMaxStep(kBig);
 571
 572   // Use only charged tracks
 573   if (( gMC->TrackCharge()       ) &&
 574       (!gMC->IsTrackStop()       ) &&
 575       (!gMC->IsTrackDisappeared())) {
 576
 577     // Inside a sensitive volume?
 578     iIdSens = gMC->CurrentVolID(icSens);
 579     if  ((iIdSens == fIdSensDr) ||
 580          (iIdSens == fIdSensAm)) {
 581
 582       iIdChamber = gMC->CurrentVolOffID(1,icChamber);
 583
 584       // The hit coordinates and charge
 585       gMC->TrackPosition(pos);
 586       hits[0] = pos[0];
 587       hits[1] = pos[1];
 588       hits[2] = pos[2];
 589
 590       // The sector number (0 - 17)
 591       // The numbering goes clockwise and starts at y = 0
 592       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 593       if (phi < 90.)
 594         phi = phi + 270.;
 595       else
 596         phi = phi -  90.;
 597       sec = ((Int_t) (phi / 20));
 598
 599       // The chamber number
 600       //   0: outer left
 601       //   1: middle left
 602       //   2: inner
 603       //   3: middle right
 604       //   4: outer right
 605       if      (iIdChamber == fIdChamber1)
 606         cha = (hits[2] < 0 ? 0 : 4);
 607       else if (iIdChamber == fIdChamber2)
 608         cha = (hits[2] < 0 ? 1 : 3);
 609       else if (iIdChamber == fIdChamber3)
 610         cha = 2;
 611
 612       // The plane number
 613       // The numbering starts at the innermost plane
 614       pla = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6 - 1;
 615
 616       // Check on selected volumes
 617       Int_t addthishit = 1;
 618       if (fSensSelect) {
 619         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 620         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 621         if (fSensSector  >= 0) {
 622           Int_t sens1  = fSensSector;
 623           Int_t sens2  = fSensSector + fSensSectorRange;
 624                 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 625                        * AliTRDgeometry::Nsect();
 626           if (sens1 < sens2) {
 627             if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 628           }
 629           else {
 630             if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 631           }
 632         }
 633       }
 634
 635       // Add this hit
 636       if (addthishit) {
 637
 638         // The detector number
 639         det = fGeometry->GetDetector(pla,cha,sec);
 640
 641         // Special hits and TR photons only in the drift region
 642         if (iIdSens == fIdSensDr) {
 643
 644           // Create some special hits with amplitude 0 at the entrance and
 645           // exit of each chamber that contain the momentum components of the particle
 646           if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
 647             gMC->TrackMomentum(mom);
 648             moms[0] = mom[0];
 649             moms[1] = mom[1];
 650             moms[2] = mom[2];
 651             AliTRDhit *hitTest = new(lhits[fNhits++]) AliTRDhit(fIshunt
 652                                                                ,gAlice->CurrentTrack()
 653                                                                ,det,moms,0);
 654             hitTest->SetTest();
 655           }
 656
 657           // Create the hits from TR photons
 658           if (fTR) CreateTRhit(det);
 659
 660         }
 661
 662         // Calculate the energy of the delta-electrons
 663         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 664         eDelta = TMath::Max(eDelta,0.0);
 665
 666         // The number of secondary electrons created
 667         qTot = ((Int_t) (eDelta / kWion) + 1);
 668
 669         // Create a new dEdx hit
 670         AliTRDhit *hit = new(lhits[fNhits++]) AliTRDhit(fIshunt
 671                                                        ,gAlice->CurrentTrack()
 672                                                        ,det,hits,qTot);
 673         if (iIdSens == fIdSensDr) {
 674           hit->SetDrift();
 675         }
 676         else {
 677           hit->SetAmplification();
 678         }
 679
 680         // Calculate the maximum step size for the next tracking step
 681         // Produce only one hit if Ekin is below cutoff
 682         aMass = gMC->TrackMass();
 683         if ((gMC->Etot() - aMass) > kEkinMinStep) {
 684
 685           // The energy loss according to Bethe Bloch
 686           iPdg  = TMath::Abs(gMC->TrackPid());
 687           if ( (iPdg != kPdgElectron) ||
 688               ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
 689             gMC->TrackMomentum(mom);
 690             pTot      = mom.Rho();
 691             betaGamma = pTot / aMass;
 692             pp        = kPrim * BetheBloch(betaGamma);
 693             // Take charge > 1 into account
 694             charge = gMC->TrackCharge();
 695             if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 696           }
 697           // Electrons above 20 Mev/c are at the plateau
 698           else {
 699             pp = kPrim * kPlateau;
 700           }
 701
 702           if (pp > 0) {
 703             do
 704               gMC->Rndm(random,1);
 705             while ((random[0] == 1.) || (random[0] == 0.));
 706             stepSize = - TMath::Log(random[0]) / pp;
 707             gMC->SetMaxStep(stepSize);
 708           }
 709
 710         }
 711
 712       }
 713
 714     }
 715
 716   }
 717
 718 }
 719
 720 //_____________________________________________________________________________
 721 Double_t AliTRDv1::BetheBloch(Double_t bg)
 722 {
 723   //
 724   // Parametrization of the Bethe-Bloch-curve
 725   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
 726   //
 727
 728   // This parameters have been adjusted to averaged values from GEANT
 729   const Double_t kP1 = 7.17960e-02;
 730   const Double_t kP2 = 8.54196;
 731   const Double_t kP3 = 1.38065e-06;
 732   const Double_t kP4 = 5.30972;
 733   const Double_t kP5 = 2.83798;
 734
 735   // This parameters have been adjusted to Xe-data found in:
 736   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 737   //const Double_t kP1 = 0.76176E-1;
 738   //const Double_t kP2 = 10.632;
 739   //const Double_t kP3 = 3.17983E-6;
 740   //const Double_t kP4 = 1.8631;
 741   //const Double_t kP5 = 1.9479;
 742
 743   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
 744   const Double_t kBgMin = 0.8;
 745   const Double_t kBBMax = 6.83298;
 746   //const Double_t kBgMin = 0.6;
 747   //const Double_t kBBMax = 17.2809;
 748   //const Double_t kBgMin = 0.4;
 749   //const Double_t kBBMax = 82.0;
 750
 751   if (bg > kBgMin) {
 752     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
 753     Double_t aa = TMath::Power(yy,kP4);
 754     Double_t bb = TMath::Power((1./bg),kP5);
 755              bb = TMath::Log(kP3 + bb);
 756     return ((kP2 - aa - bb)*kP1 / aa);
 757   }
 758   else {
 759     return kBBMax;
 760   }
 761
 762 }
 763
 764 //_____________________________________________________________________________
 765 Double_t Ermilova(Double_t *x, Double_t *)
 766 {
 767   //
 768   // Calculates the delta-ray energy distribution according to Ermilova.
 769   // Logarithmic scale !
 770   //
 771
 772   Double_t energy;
 773   Double_t dpos;
 774   Double_t dnde;
 775
 776   Int_t    pos1, pos2;
 777
 778   const Int_t kNv = 31;
 779
 780   Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
 781                      , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
 782                      , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
 783                      , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
 784                      , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
 785                      , 9.4727, 9.9035,10.3735,10.5966,10.8198
 786                      ,11.5129 };
 787
 788   Float_t vye[kNv] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
 789                      , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
 790                      ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
 791                      ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
 792                      ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
 793                      ,  0.04 ,  0.023,  0.015,  0.011,  0.01
 794                      ,  0.004 };
 795
 796   energy = x[0];
 797
 798   // Find the position
 799   pos1 = pos2 = 0;
 800   dpos = 0;
 801   do {
 802     dpos = energy - vxe[pos2++];
 803   }
 804   while (dpos > 0);
 805   pos2--;
 806   if (pos2 > kNv) pos2 = kNv - 1;
 807   pos1 = pos2 - 1;
 808
 809   // Differentiate between the sampling points
 810   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
 811
 812   return dnde;
 813
 814 }