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