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