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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Transition Radiation Detector version 1 -- slow simulator                //
  21 //                                                                           //
  22 //Begin_Html
  23 /*
  24 <img src="picts/AliTRDfullClass.gif">
  25 */
  26 //End_Html
  27 //                                                                           //
  28 //                                                                           //
  29 ///////////////////////////////////////////////////////////////////////////////
  30
  31 #include <stdlib.h>
  32
  33 #include <TF1.h>
  34 #include <TLorentzVector.h>
  35 #include <TMath.h>
  36 #include <TRandom.h>
  37 #include <TVector.h>
  38 #include <TVirtualMC.h>
  39
  40 #include "AliConst.h"
  41 #include "AliLog.h"
  42 #include "AliMC.h"
  43 #include "AliRun.h"
  44 #include "AliTRDgeometry.h"
  45 #include "AliTRDhit.h"
  46 #include "AliTRDsim.h"
  47 #include "AliTRDv1.h"
  48
  49 ClassImp(AliTRDv1)
  50
  51 //_____________________________________________________________________________
  52 AliTRDv1::AliTRDv1():AliTRD()
  53 {
  54   //
  55   // Default constructor
  56   //
  57
  58   fSensSelect        =  0;
  59   fSensPlane         = -1;
  60   fSensChamber       = -1;
  61   fSensSector        = -1;
  62   fSensSectorRange   =  0;
  63
  64   fDeltaE            = NULL;
  65   fDeltaG            = NULL;
  66   fTR                = NULL;
  67
  68   fStepSize          = 0.1;
  69   fTypeOfStepManager = 2;
  70
  71 }
  72
  73 //_____________________________________________________________________________
  74 AliTRDv1::AliTRDv1(const char *name, const char *title)
  75          :AliTRD(name, title)
  76 {
  77   //
  78   // Standard constructor for Transition Radiation Detector version 1
  79   //
  80
  81   fSensSelect        =  0;
  82   fSensPlane         = -1;
  83   fSensChamber       = -1;
  84   fSensSector        = -1;
  85   fSensSectorRange   =  0;
  86
  87   fDeltaE            = NULL;
  88   fDeltaG            = NULL;
  89   fTR                = NULL;
  90   fStepSize          = 0.1;
  91   fTypeOfStepManager = 2;
  92
  93   SetBufferSize(128000);
  94
  95 }
  96
  97 //_____________________________________________________________________________
  98 AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd)
  99 {
 100   //
 101   // Copy constructor
 102   //
 103
 104   ((AliTRDv1 &) trd).Copy(*this);
 105
 106 }
 107
 108 //_____________________________________________________________________________
 109 AliTRDv1::~AliTRDv1()
 110 {
 111   //
 112   // AliTRDv1 destructor
 113   //
 114
 115   if (fDeltaE) delete fDeltaE;
 116   if (fDeltaG) delete fDeltaG;
 117   if (fTR)     delete fTR;
 118
 119 }
 120
 121 //_____________________________________________________________________________
 122 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
 123 {
 124   //
 125   // Assignment operator
 126   //
 127
 128   if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
 129   return *this;
 130
 131 }
 132
 133 //_____________________________________________________________________________
 134 void AliTRDv1::Copy(TObject &trd) const
 135 {
 136   //
 137   // Copy function
 138   //
 139
 140   ((AliTRDv1 &) trd).fSensSelect        = fSensSelect;
 141   ((AliTRDv1 &) trd).fSensPlane         = fSensPlane;
 142   ((AliTRDv1 &) trd).fSensChamber       = fSensChamber;
 143   ((AliTRDv1 &) trd).fSensSector        = fSensSector;
 144   ((AliTRDv1 &) trd).fSensSectorRange   = fSensSectorRange;
 145
 146   ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager;
 147   ((AliTRDv1 &) trd).fStepSize          = fStepSize;
 148
 149   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
 150   fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
 151   fTR->Copy(*((AliTRDv1 &) trd).fTR);
 152
 153 }
 154
 155 //_____________________________________________________________________________
 156 void AliTRDv1::CreateGeometry()
 157 {
 158   //
 159   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
 160   // This version covers the full azimuth.
 161   //
 162
 163   // Check that FRAME is there otherwise we have no place where to put the TRD
 164   AliModule* frame = gAlice->GetModule("FRAME");
 165   if (!frame) return;
 166
 167   // Define the chambers
 168   AliTRD::CreateGeometry();
 169
 170 }
 171
 172 //_____________________________________________________________________________
 173 void AliTRDv1::CreateMaterials()
 174 {
 175   //
 176   // Create materials for the Transition Radiation Detector version 1
 177   //
 178
 179   AliTRD::CreateMaterials();
 180
 181 }
 182
 183 //_____________________________________________________________________________
 184 void AliTRDv1::CreateTRhit(Int_t det)
 185 {
 186   //
 187   // Creates an electron cluster from a TR photon.
 188   // The photon is assumed to be created a the end of the radiator. The
 189   // distance after which it deposits its energy takes into account the
 190   // absorbtion of the entrance window and of the gas mixture in drift
 191   // volume.
 192   //
 193
 194   // PDG code electron
 195   const Int_t   kPdgElectron = 11;
 196
 197   // Ionization energy
 198   const Float_t kWion        = 22.04;
 199
 200   // Maximum number of TR photons per track
 201   const Int_t   kNTR         = 50;
 202
 203   TLorentzVector mom, pos;
 204
 205   // Create TR at the entrance of the chamber
 206   if (gMC->IsTrackEntering()) {
 207
 208     // Create TR only for electrons
 209     Int_t iPdg = gMC->TrackPid();
 210     if (TMath::Abs(iPdg) != kPdgElectron) return;
 211
 212     Float_t eTR[kNTR];
 213     Int_t   nTR;
 214
 215     // Create TR photons
 216     gMC->TrackMomentum(mom);
 217     Float_t pTot = mom.Rho();
 218     fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
 219     if (nTR > kNTR) {
 220       AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
 221     }
 222
 223     // Loop through the TR photons
 224     for (Int_t iTR = 0; iTR < nTR; iTR++) {
 225
 226       Float_t energyMeV = eTR[iTR] * 0.001;
 227       Float_t energyeV  = eTR[iTR] * 1000.0;
 228       Float_t absLength = 0;
 229       Float_t sigma     = 0;
 230
 231       // Take the absorbtion in the entrance window into account
 232       Double_t muMy = fTR->GetMuMy(energyMeV);
 233       sigma = muMy * fFoilDensity;
 234       if (sigma > 0.0) {
 235         absLength = gRandom->Exp(1.0/sigma);
 236         if (absLength < AliTRDgeometry::MyThick()) continue;
 237       }
 238       else {
 239         continue;
 240       }
 241
 242       // The absorbtion cross sections in the drift gas
 243       // Gas-mixture (Xe/CO2)
 244       Double_t muXe = fTR->GetMuXe(energyMeV);
 245       Double_t muCO = fTR->GetMuCO(energyMeV);
 246       sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
 247
 248       // The distance after which the energy of the TR photon
 249       // is deposited.
 250       if (sigma > 0.0) {
 251         absLength = gRandom->Exp(1.0/sigma);
 252         if (absLength > (AliTRDgeometry::DrThick()
 253                        + AliTRDgeometry::AmThick())) {
 254           continue;
 255         }
 256       }
 257       else {
 258         continue;
 259       }
 260
 261       // The position of the absorbtion
 262       Float_t posHit[3];
 263       gMC->TrackPosition(pos);
 264       posHit[0] = pos[0] + mom[0] / pTot * absLength;
 265       posHit[1] = pos[1] + mom[1] / pTot * absLength;
 266       posHit[2] = pos[2] + mom[2] / pTot * absLength;
 267
 268       // Create the charge
 269       Int_t q = ((Int_t) (energyeV / kWion));
 270
 271       // Add the hit to the array. TR photon hits are marked
 272       // by negative charge
 273       AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE);
 274
 275     }
 276
 277   }
 278
 279 }
 280
 281 //_____________________________________________________________________________
 282 void AliTRDv1::Init()
 283 {
 284   //
 285   // Initialise Transition Radiation Detector after geometry has been built.
 286   //
 287
 288   AliTRD::Init();
 289
 290   AliDebug(1,"Slow simulator\n");
 291   if (fSensSelect) {
 292     if (fSensPlane   >= 0)
 293       AliInfo(Form("Only plane %d is sensitive"));
 294     if (fSensChamber >= 0)
 295       AliInfo(Form("Only chamber %d is sensitive",fSensChamber));
 296     if (fSensSector  >= 0) {
 297       Int_t sens1  = fSensSector;
 298       Int_t sens2  = fSensSector + fSensSectorRange;
 299             sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 300                    * AliTRDgeometry::Nsect();
 301             AliInfo(Form("Only sectors %d - %d are sensitive\n",sens1,sens2-1));
 302     }
 303   }
 304   if (fTR)
 305     AliInfo("TR simulation on")
 306   else
 307     AliInfo("TR simulation off");
 308
 309   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 310   const Float_t kPoti = 12.1;
 311   // Maximum energy (50 keV);
 312   const Float_t kEend = 50000.0;
 313   // Ermilova distribution for the delta-ray spectrum
 314   Float_t poti = TMath::Log(kPoti);
 315   Float_t eEnd = TMath::Log(kEend);
 316
 317   // Ermilova distribution for the delta-ray spectrum
 318   fDeltaE = new TF1("deltae" ,Ermilova    ,poti,eEnd,0);
 319
 320   // Geant3 distribution for the delta-ray spectrum
 321   fDeltaG = new TF1("deltaeg",IntSpecGeant,poti,eEnd,0);
 322
 323   AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
 324
 325 }
 326
 327 //_____________________________________________________________________________
 328 AliTRDsim *AliTRDv1::CreateTR()
 329 {
 330   //
 331   // Enables the simulation of TR
 332   //
 333
 334   fTR = new AliTRDsim();
 335   return fTR;
 336
 337 }
 338
 339 //_____________________________________________________________________________
 340 void AliTRDv1::SetSensPlane(Int_t iplane)
 341 {
 342   //
 343   // Defines the hit-sensitive plane (0-5)
 344   //
 345
 346   if ((iplane < 0) || (iplane > 5)) {
 347     AliWarning(Form("Wrong input value:%d",iplane));
 348     AliWarning("Use standard setting");
 349     fSensPlane  = -1;
 350     fSensSelect =  0;
 351     return;
 352   }
 353
 354   fSensSelect = 1;
 355   fSensPlane  = iplane;
 356
 357 }
 358
 359 //_____________________________________________________________________________
 360 void AliTRDv1::SetSensChamber(Int_t ichamber)
 361 {
 362   //
 363   // Defines the hit-sensitive chamber (0-4)
 364   //
 365
 366   if ((ichamber < 0) || (ichamber > 4)) {
 367     AliWarning(Form("Wrong input value: %d",ichamber));
 368     AliWarning("Use standard setting");
 369     fSensChamber = -1;
 370     fSensSelect  =  0;
 371     return;
 372   }
 373
 374   fSensSelect  = 1;
 375   fSensChamber = ichamber;
 376
 377 }
 378
 379 //_____________________________________________________________________________
 380 void AliTRDv1::SetSensSector(Int_t isector)
 381 {
 382   //
 383   // Defines the hit-sensitive sector (0-17)
 384   //
 385
 386   SetSensSector(isector,1);
 387
 388 }
 389
 390 //_____________________________________________________________________________
 391 void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
 392 {
 393   //
 394   // Defines a range of hit-sensitive sectors. The range is defined by
 395   // <isector> (0-17) as the starting point and <nsector> as the number
 396   // of sectors to be included.
 397   //
 398
 399   if ((isector < 0) || (isector > 17)) {
 400     AliWarning(Form("Wrong input value <isector>: %d",isector));
 401     AliWarning("Use standard setting");
 402     fSensSector      = -1;
 403     fSensSectorRange =  0;
 404     fSensSelect      =  0;
 405     return;
 406   }
 407
 408   if ((nsector < 1) || (nsector > 18)) {
 409     AliWarning(Form("Wrong input value <nsector>: %d",nsector));
 410     AliWarning("Use standard setting");
 411     fSensSector      = -1;
 412     fSensSectorRange =  0;
 413     fSensSelect      =  0;
 414     return;
 415   }
 416
 417   fSensSelect      = 1;
 418   fSensSector      = isector;
 419   fSensSectorRange = nsector;
 420
 421 }
 422
 423 //_____________________________________________________________________________
 424 void AliTRDv1::StepManager()
 425 {
 426   //
 427   // Slow simulator. Every charged track produces electron cluster as hits
 428   // along its path across the drift volume.
 429   //
 430
 431   switch (fTypeOfStepManager) {
 432     case 0  : StepManagerErmilova();  break;  // 0 is Ermilova
 433     case 1  : StepManagerGeant();     break;  // 1 is Geant
 434     case 2  : StepManagerFixedStep(); break;  // 2 is fixed step
 435     default : AliWarning("Not a valid Step Manager.");
 436   }
 437
 438 }
 439
 440 //_____________________________________________________________________________
 441 void AliTRDv1::SelectStepManager(Int_t t)
 442 {
 443   //
 444   // Selects a step manager type:
 445   //   0 - Ermilova
 446   //   1 - Geant3
 447   //   2 - Fixed step size
 448   //
 449
 450   if (t == 1) {
 451     AliWarning("Sorry, Geant parametrization step manager is not implemented yet. Please ask K.Oyama for detail.");
 452   }
 453
 454   fTypeOfStepManager = t;
 455   AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
 456
 457 }
 458
 459 //_____________________________________________________________________________
 460 void AliTRDv1::StepManagerGeant()
 461 {
 462   //
 463   // Slow simulator. Every charged track produces electron cluster as hits
 464   // along its path across the drift volume. The step size is set acording
 465   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 466   // a spectrum taken from Geant3.
 467   //
 468
 469   AliWarning("Not implemented yet.");
 470
 471 }
 472
 473 //_____________________________________________________________________________
 474 void AliTRDv1::StepManagerErmilova()
 475 {
 476   //
 477   // Slow simulator. Every charged track produces electron cluster as hits
 478   // along its path across the drift volume. The step size is set acording
 479   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 480   // a spectrum taken from Ermilova et al.
 481   //
 482
 483   Int_t    pla = 0;
 484   Int_t    cha = 0;
 485   Int_t    sec = 0;
 486   Int_t    det = 0;
 487   Int_t    iPdg;
 488   Int_t    qTot;
 489
 490   Float_t  hits[3];
 491   Double_t random[1];
 492   Float_t  charge;
 493   Float_t  aMass;
 494
 495   Double_t pTot = 0;
 496   Double_t eDelta;
 497   Double_t betaGamma, pp;
 498   Double_t stepSize;
 499
 500   Bool_t   drRegion = kFALSE;
 501   Bool_t   amRegion = kFALSE;
 502
 503   TString  cIdCurrent;
 504   TString  cIdSensDr = "J";
 505   TString  cIdSensAm = "K";
 506   Char_t   cIdChamber[3];
 507            cIdChamber[2] = 0;
 508
 509   TLorentzVector pos, mom;
 510
 511   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 512   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 513   const Int_t    kNdetsec     = kNplan * kNcham;
 514
 515   const Double_t kBig         = 1.0E+12; // Infinitely big
 516   const Float_t  kWion        = 22.04;   // Ionization energy
 517   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
 518
 519   // Minimum energy for the step size adjustment
 520   const Float_t  kEkinMinStep = 1.0e-5;
 521
 522   // Plateau value of the energy-loss for electron in xenon
 523   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 524   //const Double_t kPlateau = 1.70;
 525   // the averaged value (26/3/99)
 526   const Float_t  kPlateau     = 1.55;
 527
 528   const Float_t  kPrim        = 48.0;  // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 529   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 530   const Float_t  kPoti        = 12.1;
 531
 532   const Int_t    kPdgElectron = 11;  // PDG code electron
 533
 534   // Set the maximum step size to a very large number for all
 535   // neutral particles and those outside the driftvolume
 536   gMC->SetMaxStep(kBig);
 537
 538   // Use only charged tracks
 539   if (( gMC->TrackCharge()       ) &&
 540       (!gMC->IsTrackStop()       ) &&
 541       (!gMC->IsTrackDisappeared())) {
 542
 543     // Inside a sensitive volume?
 544     drRegion = kFALSE;
 545     amRegion = kFALSE;
 546     cIdCurrent = gMC->CurrentVolName();
 547     if (cIdSensDr == cIdCurrent[1]) {
 548       drRegion = kTRUE;
 549     }
 550     if (cIdSensAm == cIdCurrent[1]) {
 551       amRegion = kTRUE;
 552     }
 553     if (drRegion || amRegion) {
 554
 555       // The hit coordinates and charge
 556       gMC->TrackPosition(pos);
 557       hits[0] = pos[0];
 558       hits[1] = pos[1];
 559       hits[2] = pos[2];
 560
 561       // The sector number (0 - 17)
 562       // The numbering goes clockwise and starts at y = 0
 563       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 564       if (phi < 90.)
 565         phi = phi + 270.;
 566       else
 567         phi = phi -  90.;
 568       sec = ((Int_t) (phi / 20));
 569
 570       // The plane and chamber number
 571       cIdChamber[0] = cIdCurrent[2];
 572       cIdChamber[1] = cIdCurrent[3];
 573       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 574       cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
 575       pla = ((Int_t) idChamber % kNplan);
 576
 577       // Check on selected volumes
 578       Int_t addthishit = 1;
 579       if (fSensSelect) {
 580         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 581         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 582         if (fSensSector  >= 0) {
 583           Int_t sens1  = fSensSector;
 584           Int_t sens2  = fSensSector + fSensSectorRange;
 585                 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 586                        * AliTRDgeometry::Nsect();
 587           if (sens1 < sens2) {
 588             if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 589           }
 590           else {
 591             if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 592           }
 593         }
 594       }
 595
 596       // Add this hit
 597       if (addthishit) {
 598
 599         // The detector number
 600         det = fGeometry->GetDetector(pla,cha,sec);
 601
 602         // Special hits only in the drift region
 603         if (drRegion) {
 604
 605           // Create a track reference at the entrance and
 606           // exit of each chamber that contain the
 607           // momentum components of the particle
 608           if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
 609             gMC->TrackMomentum(mom);
 610             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 611           }
 612
 613           // Create the hits from TR photons
 614           if (fTR) CreateTRhit(det);
 615
 616         }
 617
 618         // Calculate the energy of the delta-electrons
 619         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 620         eDelta = TMath::Max(eDelta,0.0);
 621
 622         // The number of secondary electrons created
 623         qTot = ((Int_t) (eDelta / kWion) + 1);
 624
 625         // Create a new dEdx hit
 626         if (drRegion) {
 627           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 628                 ,det,hits,qTot,kTRUE);
 629         }
 630         else {
 631           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 632                 ,det,hits,qTot,kFALSE);
 633         }
 634
 635         // Calculate the maximum step size for the next tracking step
 636         // Produce only one hit if Ekin is below cutoff
 637         aMass = gMC->TrackMass();
 638         if ((gMC->Etot() - aMass) > kEkinMinStep) {
 639
 640           // The energy loss according to Bethe Bloch
 641           iPdg  = TMath::Abs(gMC->TrackPid());
 642           if ( (iPdg != kPdgElectron) ||
 643               ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
 644             gMC->TrackMomentum(mom);
 645             pTot      = mom.Rho();
 646             betaGamma = pTot / aMass;
 647             pp        = kPrim * BetheBloch(betaGamma);
 648             // Take charge > 1 into account
 649             charge = gMC->TrackCharge();
 650             if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 651           }
 652           // Electrons above 20 Mev/c are at the plateau
 653           else {
 654             pp = kPrim * kPlateau;
 655           }
 656
 657           if (pp > 0) {
 658             do
 659             gMC->GetRandom()->RndmArray(1, random);
 660             while ((random[0] == 1.) || (random[0] == 0.));
 661             stepSize = - TMath::Log(random[0]) / pp;
 662             gMC->SetMaxStep(stepSize);
 663           }
 664
 665         }
 666
 667       }
 668
 669     }
 670
 671   }
 672
 673 }
 674
 675 //_____________________________________________________________________________
 676 void AliTRDv1::StepManagerFixedStep()
 677 {
 678   //
 679   // Slow simulator. Every charged track produces electron cluster as hits
 680   // along its path across the drift volume. The step size is fixed in
 681   // this version of the step manager.
 682   //
 683
 684   Int_t    pla = 0;
 685   Int_t    cha = 0;
 686   Int_t    sec = 0;
 687   Int_t    det = 0;
 688   Int_t    qTot;
 689
 690   Float_t  hits[3];
 691   Double_t eDep;
 692
 693   Bool_t   drRegion = kFALSE;
 694   Bool_t   amRegion = kFALSE;
 695
 696   TString  cIdCurrent;
 697   TString  cIdSensDr = "J";
 698   TString  cIdSensAm = "K";
 699   Char_t   cIdChamber[3];
 700   cIdChamber[2] = 0;
 701
 702   TLorentzVector pos, mom;
 703
 704   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 705   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 706   const Int_t    kNdetsec     = kNplan * kNcham;
 707
 708   const Double_t kBig         = 1.0E+12;
 709
 710   const Float_t  kWion        = 22.04;   // Ionization energy
 711   const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment
 712
 713   // Set the maximum step size to a very large number for all
 714   // neutral particles and those outside the driftvolume
 715   gMC->SetMaxStep(kBig);
 716
 717   // If not charged track or already stopped or disappeared, just return.
 718   if ((!gMC->TrackCharge()) ||
 719         gMC->IsTrackStop()  ||
 720         gMC->IsTrackDisappeared()) return;
 721
 722   // Inside a sensitive volume?
 723   cIdCurrent = gMC->CurrentVolName();
 724
 725   if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
 726   if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
 727
 728   if ((!drRegion) && (!amRegion)) return;
 729
 730   // The hit coordinates and charge
 731   gMC->TrackPosition(pos);
 732   hits[0] = pos[0];
 733   hits[1] = pos[1];
 734   hits[2] = pos[2];
 735
 736   // The sector number (0 - 17)
 737   // The numbering goes clockwise and starts at y = 0
 738   Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 739   if (phi < 90.) phi += 270.;
 740   else           phi -=  90.;
 741   sec = ((Int_t) (phi / 20.));
 742
 743   // The plane and chamber number
 744   cIdChamber[0] = cIdCurrent[2];
 745   cIdChamber[1] = cIdCurrent[3];
 746   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 747   cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
 748   pla = ((Int_t) idChamber % kNplan);
 749
 750   // Check on selected volumes
 751   Int_t addthishit = 1;
 752   if(fSensSelect) {
 753     if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 754     if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 755     if (fSensSector  >= 0) {
 756       Int_t sens1  = fSensSector;
 757       Int_t sens2  = fSensSector + fSensSectorRange;
 758       sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect();
 759       if (sens1 < sens2) {
 760         if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 761       }
 762       else {
 763         if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 764       }
 765     }
 766   }
 767
 768   if (!addthishit) return;
 769
 770   det = fGeometry->GetDetector(pla,cha,sec);  // The detector number
 771
 772   Int_t trkStat = 0;  // 0: InFlight 1:Entering 2:Exiting
 773
 774   // Special hits only in the drift region
 775   if (drRegion) {
 776
 777     // Create a track reference at the entrance and exit of each
 778     // chamber that contain the momentum components of the particle
 779
 780     if (gMC->IsTrackEntering()) {
 781       gMC->TrackMomentum(mom);
 782       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 783       trkStat = 1;
 784     }
 785     if (gMC->IsTrackExiting()) {
 786       gMC->TrackMomentum(mom);
 787       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 788       trkStat = 2;
 789     }
 790
 791     // Create the hits from TR photons
 792     if (fTR) CreateTRhit(det);
 793
 794   }
 795
 796   // Calculate the charge according to GEANT Edep
 797   // Create a new dEdx hit
 798   eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
 799   qTot = (Int_t) (eDep / kWion);
 800   AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot,drRegion);
 801
 802   // Set Maximum Step Size
 803   // Produce only one hit if Ekin is below cutoff
 804   if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return;
 805   gMC->SetMaxStep(fStepSize);
 806
 807 }
 808
 809 //_____________________________________________________________________________
 810 Double_t AliTRDv1::BetheBloch(Double_t bg)
 811 {
 812   //
 813   // Parametrization of the Bethe-Bloch-curve
 814   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
 815   //
 816
 817   // This parameters have been adjusted to averaged values from GEANT
 818   const Double_t kP1 = 7.17960e-02;
 819   const Double_t kP2 = 8.54196;
 820   const Double_t kP3 = 1.38065e-06;
 821   const Double_t kP4 = 5.30972;
 822   const Double_t kP5 = 2.83798;
 823
 824   // This parameters have been adjusted to Xe-data found in:
 825   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 826   //const Double_t kP1 = 0.76176E-1;
 827   //const Double_t kP2 = 10.632;
 828   //const Double_t kP3 = 3.17983E-6;
 829   //const Double_t kP4 = 1.8631;
 830   //const Double_t kP5 = 1.9479;
 831
 832   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
 833   const Double_t kBgMin = 0.8;
 834   const Double_t kBBMax = 6.83298;
 835   //const Double_t kBgMin = 0.6;
 836   //const Double_t kBBMax = 17.2809;
 837   //const Double_t kBgMin = 0.4;
 838   //const Double_t kBBMax = 82.0;
 839
 840   if (bg > kBgMin) {
 841     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
 842     Double_t aa = TMath::Power(yy,kP4);
 843     Double_t bb = TMath::Power((1./bg),kP5);
 844              bb = TMath::Log(kP3 + bb);
 845     return ((kP2 - aa - bb)*kP1 / aa);
 846   }
 847   else {
 848     return kBBMax;
 849   }
 850
 851 }
 852
 853 //_____________________________________________________________________________
 854 Double_t BetheBlochGeant(Double_t bg)
 855 {
 856   //
 857   // Return dN/dx (number of primary collisions per centimeter)
 858   // for given beta*gamma factor.
 859   //
 860   // Implemented by K.Oyama according to GEANT 3 parametrization shown in
 861   // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
 862   // This must be used as a set with IntSpecGeant.
 863   //
 864
 865   Double_t arr_g[20] = {
 866     1.100000,   1.200000,    1.300000,    1.500000,
 867     1.800000,   2.000000,    2.500000,    3.000000,
 868     4.000000,   7.000000,    10.000000,   20.000000,
 869     40.000000,  70.000000,   100.000000,  300.000000,
 870     600.000000, 1000.000000, 3000.000000, 10000.000000 };
 871
 872   Double_t arr_nc[20] = {
 873     75.009056,   45.508083,   35.299252,   27.116327,
 874     22.734999,   21.411915,   19.934095,   19.449375,
 875     19.344431,   20.185553,   21.027925,   22.912676,
 876     24.933352,   26.504053,   27.387468,   29.566597,
 877     30.353779,   30.787134,   31.129285,   31.157350 };
 878
 879   // betagamma to gamma
 880   Double_t g = TMath::Sqrt( 1. + bg*bg );
 881
 882   // Find the index just before the point we need.
 883   int i;
 884   for( i = 0 ; i < 18 ; i++ )
 885     if( arr_g[i] < g && arr_g[i+1] > g )
 886       break;
 887
 888   // Simple interpolation.
 889   Double_t pp = ((arr_nc[i+1] - arr_nc[i]) /
 890                  (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i];
 891
 892   return pp;
 893
 894 }
 895
 896 //_____________________________________________________________________________
 897 Double_t Ermilova(Double_t *x, Double_t *)
 898 {
 899   //
 900   // Calculates the delta-ray energy distribution according to Ermilova.
 901   // Logarithmic scale !
 902   //
 903
 904   Double_t energy;
 905   Double_t dpos;
 906   Double_t dnde;
 907
 908   Int_t    pos1, pos2;
 909
 910   const Int_t kNv = 31;
 911
 912   Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
 913                      , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
 914                      , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
 915                      , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
 916                      , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
 917                      , 9.4727, 9.9035,10.3735,10.5966,10.8198
 918                      ,11.5129 };
 919
 920   Float_t vye[kNv] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
 921                      , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
 922                      ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
 923                      ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
 924                      ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
 925                      ,  0.04 ,  0.023,  0.015,  0.011,  0.01
 926                      ,  0.004 };
 927
 928   energy = x[0];
 929
 930   // Find the position
 931   pos1 = pos2 = 0;
 932   dpos = 0;
 933   do {
 934     dpos = energy - vxe[pos2++];
 935   }
 936   while (dpos > 0);
 937   pos2--;
 938   if (pos2 > kNv) pos2 = kNv - 1;
 939   pos1 = pos2 - 1;
 940
 941   // Differentiate between the sampling points
 942   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
 943
 944   return dnde;
 945
 946 }
 947
 948 //_____________________________________________________________________________
 949 Double_t IntSpecGeant(Double_t *x, Double_t *)
 950 {
 951   //
 952   // Integrated spectrum from Geant3
 953   //
 954
 955   const Int_t n_pts = 83;
 956   Double_t arr_e[n_pts] = {
 957     2.421257,  2.483278,  2.534301,  2.592230,
 958     2.672067,  2.813299,  3.015059,  3.216819,
 959     3.418579,  3.620338,  3.868209,  3.920198,
 960     3.978284,  4.063923,  4.186264,  4.308605,
 961     4.430946,  4.553288,  4.724261,  4.837736,
 962     4.999842,  5.161949,  5.324056,  5.486163,
 963     5.679688,  5.752998,  5.857728,  5.962457,
 964     6.067185,  6.171914,  6.315653,  6.393674,
 965     6.471694,  6.539689,  6.597658,  6.655627,
 966     6.710957,  6.763648,  6.816338,  6.876198,
 967     6.943227,  7.010257,  7.106285,  7.252151,
 968     7.460531,  7.668911,  7.877290,  8.085670,
 969     8.302979,  8.353585,  8.413120,  8.483500,
 970     8.541030,  8.592857,  8.668865,  8.820485,
 971     9.037086,  9.253686,  9.470286,  9.686887,
 972     9.930838,  9.994655, 10.085822, 10.176990,
 973     10.268158, 10.359325, 10.503614, 10.627565,
 974     10.804637, 10.981709, 11.158781, 11.335854,
 975     11.593397, 11.781165, 12.049404, 12.317644,
 976     12.585884, 12.854123, 14.278421, 16.975889,
 977     20.829416, 24.682943, 28.536469
 978   };
 979   Double_t arr_dndx[n_pts] = {
 980     19.344431, 18.664679, 18.136106, 17.567745,
 981     16.836426, 15.677382, 14.281277, 13.140237,
 982     12.207677, 11.445510, 10.697049, 10.562296,
 983     10.414673, 10.182341,  9.775256,  9.172330,
 984     8.240271,  6.898587,  4.808303,  3.889751,
 985     3.345288,  3.093431,  2.897347,  2.692470,
 986     2.436222,  2.340029,  2.208579,  2.086489,
 987     1.975535,  1.876519,  1.759626,  1.705024,
 988     1.656374,  1.502638,  1.330566,  1.200697,
 989     1.101168,  1.019323,  0.943867,  0.851951,
 990     0.755229,  0.671576,  0.570675,  0.449672,
 991     0.326722,  0.244225,  0.188225,  0.149608,
 992     0.121529,  0.116289,  0.110636,  0.103490,
 993     0.096147,  0.089191,  0.079780,  0.063927,
 994     0.047642,  0.036341,  0.028250,  0.022285,
 995     0.017291,  0.016211,  0.014802,  0.013533,
 996     0.012388,  0.011352,  0.009803,  0.008537,
 997     0.007039,  0.005829,  0.004843,  0.004034,
 998     0.003101,  0.002564,  0.001956,  0.001494,
 999     0.001142,  0.000873,  0.000210,  0.000014,
1000     0.000000,  0.000000,  0.000000
1001   };
1002
1003   Int_t i;
1004   Double_t energy = x[0];
1005   Double_t dnde;
1006
1007   for( i = 0 ; i < n_pts ; i++ )
1008     if( energy < arr_e[i] ) break;
1009
1010   if( i == 0 )
1011     AliErrorGeneral("AliTRDv1","Given energy value is too small or zero");
1012
1013   // Interpolate
1014   dnde = (arr_dndx[i-1] - arr_dndx[i]) / (arr_e[i] - arr_e[i-1]);
1015
1016   return dnde;
1017
1018 }