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