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