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