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 ////////////////////////////////////////////////////////////////////////////
  23
  24 #include <stdlib.h>
  25
  26 #include <TF1.h>
  27 #include <TLorentzVector.h>
  28 #include <TMath.h>
  29 #include <TRandom.h>
  30 #include <TVector.h>
  31 #include <TVirtualMC.h>
  32 #include <TGeoManager.h>
  33 #include <TGeoPhysicalNode.h>
  34
  35 #include "AliConst.h"
  36 #include "AliLog.h"
  37 #include "AliTrackReference.h"
  38 #include "AliMC.h"
  39 #include "AliRun.h"
  40
  41 #include "AliTRDgeometry.h"
  42 #include "AliTRDhit.h"
  43 #include "AliTRDsimTR.h"
  44 #include "AliTRDv1.h"
  45
  46 ClassImp(AliTRDv1)
  47
  48 //_____________________________________________________________________________
  49 AliTRDv1::AliTRDv1()
  50   :AliTRD()
  51   ,fTRon(kFALSE)
  52   ,fTR(NULL)
  53   ,fTypeOfStepManager(0)
  54   ,fStepSize(0)
  55   ,fDeltaE(NULL)
  56   ,fDeltaG(NULL)
  57   ,fTrackLength0(0)
  58   ,fPrimaryTrackPid(0)
  59 {
  60   //
  61   // Default constructor
  62   //
  63
  64 }
  65
  66 //_____________________________________________________________________________
  67 AliTRDv1::AliTRDv1(const char *name, const char *title)
  68   :AliTRD(name,title)
  69   ,fTRon(kTRUE)
  70   ,fTR(NULL)
  71   ,fTypeOfStepManager(2)
  72   ,fStepSize(0.1)
  73   ,fDeltaE(NULL)
  74   ,fDeltaG(NULL)
  75   ,fTrackLength0(0)
  76   ,fPrimaryTrackPid(0)
  77 {
  78   //
  79   // Standard constructor for Transition Radiation Detector version 1
  80   //
  81
  82   SetBufferSize(128000);
  83
  84 }
  85
  86 //_____________________________________________________________________________
  87 AliTRDv1::~AliTRDv1()
  88 {
  89   //
  90   // AliTRDv1 destructor
  91   //
  92
  93   if (fDeltaE) {
  94     delete fDeltaE;
  95     fDeltaE = 0;
  96   }
  97
  98   if (fDeltaG) {
  99     delete fDeltaG;
 100     fDeltaG = 0;
 101   }
 102
 103   if (fTR) {
 104     delete fTR;
 105     fTR     = 0;
 106   }
 107
 108 }
 109
 110 //_____________________________________________________________________________
 111 void AliTRDv1::AddAlignableVolumes() const
 112 {
 113   //
 114   // Create entries for alignable volumes associating the symbolic volume
 115   // name with the corresponding volume path. Needs to be syncronized with
 116   // eventual changes in the geometry.
 117   //
 118
 119   TString volPath;
 120   TString symName;
 121
 122   TString vpStr  = "ALIC_1/B077_1/BSEGMO";
 123   TString vpApp1 = "_1/BTRD";
 124   TString vpApp2 = "_1";
 125   TString vpApp3 = "/UTR1_1/UTS1_1/UTI1_1/UT";
 126
 127   TString snStr  = "TRD/sm";
 128   TString snApp1 = "/st";
 129   TString snApp2 = "/pl";
 130
 131   //
 132   // The super modules
 133   // The symbolic names are: TRD/sm00
 134   //                           ...
 135   //                         TRD/sm17
 136   //
 137   for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
 138
 139     volPath  = vpStr;
 140     volPath += isect;
 141     volPath += vpApp1;
 142     volPath += isect;
 143     volPath += vpApp2;
 144
 145     symName  = snStr;
 146     symName += Form("%02d",isect);
 147
 148     gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
 149
 150   }
 151
 152   //
 153   // The readout chambers
 154   // The symbolic names are: TRD/sm00/st0/pl0
 155   //                           ...
 156   //                         TRD/sm17/st4/pl5
 157   //
 158   for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
 159
 160     if (fGeometry->GetSMstatus(isect) == 0) continue;
 161
 162     for (Int_t icham = 0; icham < AliTRDgeometry::Ncham(); icham++) {
 163       for (Int_t iplan = 0; iplan < AliTRDgeometry::Nplan(); iplan++) {
 164
 165         Int_t idet = AliTRDgeometry::GetDetectorSec(iplan,icham);
 166
 167         volPath  = vpStr;
 168         volPath += isect;
 169         volPath += vpApp1;
 170         volPath += isect;
 171         volPath += vpApp2;
 172         volPath += vpApp3;
 173         volPath += Form("%02d",idet);
 174         volPath += vpApp2;
 175
 176         symName  = snStr;
 177         symName += Form("%02d",isect);
 178         symName += snApp1;
 179         symName += icham;
 180         symName += snApp2;
 181         symName += iplan;
 182
 183         TGeoPNEntry *alignableEntry =
 184           gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
 185
 186         // Add the tracking to local matrix following the TPC example
 187         if (alignableEntry) {
 188           const char *path = alignableEntry->GetTitle();
 189           if (!gGeoManager->cd(path)) {
 190             AliFatal(Form("Volume path %s not valid!",path));
 191           }
 192           TGeoHMatrix *globMatrix = gGeoManager->GetCurrentMatrix();
 193           Double_t sectorAngle = 20.0 * (isect % 18) + 10.0;
 194           TGeoHMatrix *t2lMatrix  = new TGeoHMatrix();
 195           t2lMatrix->RotateZ(sectorAngle);
 196           t2lMatrix->MultiplyLeft(&(globMatrix->Inverse()));
 197           alignableEntry->SetMatrix(t2lMatrix);
 198         }
 199         else {
 200           AliError(Form("Alignable entry %s is not valid!",symName.Data()));
 201         }
 202
 203       }
 204     }
 205   }
 206
 207 }
 208
 209 //_____________________________________________________________________________
 210 void AliTRDv1::CreateGeometry()
 211 {
 212   //
 213   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
 214   // This version covers the full azimuth.
 215   //
 216
 217   // Check that FRAME is there otherwise we have no place where to put the TRD
 218   AliModule* frame = gAlice->GetModule("FRAME");
 219   if (!frame) {
 220     AliError("TRD needs FRAME to be present\n");
 221     return;
 222   }
 223
 224   // Define the chambers
 225   AliTRD::CreateGeometry();
 226
 227 }
 228
 229 //_____________________________________________________________________________
 230 void AliTRDv1::CreateMaterials()
 231 {
 232   //
 233   // Create materials for the Transition Radiation Detector version 1
 234   //
 235
 236   AliTRD::CreateMaterials();
 237
 238 }
 239
 240 //_____________________________________________________________________________
 241 void AliTRDv1::CreateTRhit(Int_t det)
 242 {
 243   //
 244   // Creates an electron cluster from a TR photon.
 245   // The photon is assumed to be created a the end of the radiator. The
 246   // distance after which it deposits its energy takes into account the
 247   // absorbtion of the entrance window and of the gas mixture in drift
 248   // volume.
 249   //
 250
 251   // Ionization energy
 252   const Float_t kWion        = 23.53;
 253
 254   // Maximum number of TR photons per track
 255   const Int_t   kNTR         = 50;
 256
 257   TLorentzVector mom;
 258   TLorentzVector pos;
 259
 260   Float_t eTR[kNTR];
 261   Int_t   nTR;
 262
 263   // Create TR photons
 264   gMC->TrackMomentum(mom);
 265   Float_t pTot = mom.Rho();
 266   fTR->CreatePhotons(11,pTot,nTR,eTR);
 267   if (nTR > kNTR) {
 268     AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
 269   }
 270
 271   // Loop through the TR photons
 272   for (Int_t iTR = 0; iTR < nTR; iTR++) {
 273
 274     Float_t energyMeV = eTR[iTR] * 0.001;
 275     Float_t energyeV  = eTR[iTR] * 1000.0;
 276     Float_t absLength = 0.0;
 277     Float_t sigma     = 0.0;
 278
 279     // Take the absorbtion in the entrance window into account
 280     Double_t muMy = fTR->GetMuMy(energyMeV);
 281     sigma         = muMy * fFoilDensity;
 282     if (sigma > 0.0) {
 283       absLength = gRandom->Exp(1.0/sigma);
 284       if (absLength < AliTRDgeometry::MyThick()) {
 285         continue;
 286       }
 287     }
 288     else {
 289       continue;
 290     }
 291
 292     // The absorbtion cross sections in the drift gas
 293     // Gas-mixture (Xe/CO2)
 294     Double_t muXe = fTR->GetMuXe(energyMeV);
 295     Double_t muCO = fTR->GetMuCO(energyMeV);
 296     sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
 297
 298     // The distance after which the energy of the TR photon
 299     // is deposited.
 300     if (sigma > 0.0) {
 301       absLength = gRandom->Exp(1.0/sigma);
 302       if (absLength > (AliTRDgeometry::DrThick()
 303                      + AliTRDgeometry::AmThick())) {
 304         continue;
 305       }
 306     }
 307     else {
 308       continue;
 309     }
 310
 311     // The position of the absorbtion
 312     Float_t posHit[3];
 313     gMC->TrackPosition(pos);
 314     posHit[0] = pos[0] + mom[0] / pTot * absLength;
 315     posHit[1] = pos[1] + mom[1] / pTot * absLength;
 316     posHit[2] = pos[2] + mom[2] / pTot * absLength;
 317
 318     // Create the charge
 319     Int_t q = ((Int_t) (energyeV / kWion));
 320
 321     // Add the hit to the array. TR photon hits are marked
 322     // by negative charge
 323     AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 324           ,det
 325           ,posHit
 326           ,-q
 327           ,gMC->TrackTime()*1.0e06
 328           ,kTRUE);
 329
 330   }
 331
 332 }
 333
 334 //_____________________________________________________________________________
 335 void AliTRDv1::Init()
 336 {
 337   //
 338   // Initialise Transition Radiation Detector after geometry has been built.
 339   //
 340
 341   AliTRD::Init();
 342
 343   AliDebug(1,"Slow simulator\n");
 344
 345   // Switch on TR simulation as default
 346   if (!fTRon) {
 347     AliInfo("TR simulation off");
 348   }
 349   else {
 350     fTR = new AliTRDsimTR();
 351   }
 352
 353   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 354   const Float_t kPoti = 12.1;
 355   // Maximum energy (50 keV);
 356   const Float_t kEend = 50000.0;
 357   // Ermilova distribution for the delta-ray spectrum
 358   Float_t poti        = TMath::Log(kPoti);
 359   Float_t eEnd        = TMath::Log(kEend);
 360
 361   // Ermilova distribution for the delta-ray spectrum
 362   fDeltaE = new TF1("deltae" ,Ermilova ,poti,eEnd,0);
 363
 364   // Geant3 distribution for the delta-ray spectrum
 365   fDeltaG = new TF1("deltag",IntSpecGeant,2.421257,28.536469,0);
 366
 367   AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
 368
 369 }
 370
 371 //_____________________________________________________________________________
 372 void AliTRDv1::StepManager()
 373 {
 374   //
 375   // Slow simulator. Every charged track produces electron cluster as hits
 376   // along its path across the drift volume.
 377   //
 378
 379   switch (fTypeOfStepManager) {
 380    case 0:
 381     StepManagerErmilova();
 382     break;
 383    case 1:
 384     StepManagerGeant();
 385     break;
 386    case 2:
 387     StepManagerFixedStep();
 388     break;
 389    default:
 390     AliWarning("Not a valid Step Manager.");
 391   }
 392
 393 }
 394
 395 //_____________________________________________________________________________
 396 void AliTRDv1::SelectStepManager(Int_t t)
 397 {
 398   //
 399   // Selects a step manager type:
 400   //   0 - Ermilova
 401   //   1 - Geant3
 402   //   2 - Fixed step size
 403   //
 404
 405   fTypeOfStepManager = t;
 406   AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
 407
 408 }
 409
 410 //_____________________________________________________________________________
 411 void AliTRDv1::StepManagerGeant()
 412 {
 413   //
 414   // Slow simulator. Every charged track produces electron cluster as hits
 415   // along its path across the drift volume. The step size is set acording
 416   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 417   // a spectrum taken from Geant3.
 418   //
 419   // Version by A. Bercuci
 420   //
 421
 422   Int_t    pla = 0;
 423   Int_t    cha = 0;
 424   Int_t    sec = 0;
 425   Int_t    det = 0;
 426   Int_t    iPdg;
 427   Int_t    qTot;
 428
 429   Float_t  hits[3];
 430   Float_t  charge;
 431   Float_t  aMass;
 432
 433   Double_t pTot     = 0;
 434   Double_t eDelta;
 435   Double_t betaGamma;
 436   Double_t pp;
 437   Double_t stepSize = 0;
 438
 439   Bool_t   drRegion = kFALSE;
 440   Bool_t   amRegion = kFALSE;
 441
 442   TString  cIdPath;
 443   Char_t   cIdSector[3];
 444            cIdSector[2]  = 0;
 445
 446   TString  cIdCurrent;
 447   TString  cIdSensDr = "J";
 448   TString  cIdSensAm = "K";
 449   Char_t   cIdChamber[3];
 450            cIdChamber[2] = 0;
 451
 452   TLorentzVector pos;
 453   TLorentzVector mom;
 454
 455   TArrayI        processes;
 456
 457   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 458   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 459   const Int_t    kNdetsec     = kNplan * kNcham;
 460
 461   const Double_t kBig         = 1.0e+12; // Infinitely big
 462   const Float_t  kWion        = 23.53;   // Ionization energy
 463   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
 464
 465   // Minimum energy for the step size adjustment
 466   const Float_t  kEkinMinStep = 1.0e-5;
 467   // energy threshold for production of delta electrons
 468   const Float_t  kECut        = 1.0e4;
 469   // Parameters entering the parametrized range for delta electrons
 470   const Float_t  kRa          = 5.37e-4;
 471   const Float_t  kRb          = 0.9815;
 472   const Float_t  kRc          = 3.123e-3;
 473   // Gas density -> To be made user adjustable !
 474   // [0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
 475   const Float_t  kRho         = 0.004945 ;
 476
 477   // Plateau value of the energy-loss for electron in xenon
 478   // The averaged value (26/3/99)
 479   const Float_t  kPlateau     = 1.55;
 480   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 481   const Float_t  kPrim        = 19.34;
 482   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 483   const Float_t  kPoti        = 12.1;
 484   // PDG code electron
 485   const Int_t    kPdgElectron = 11;
 486
 487   // Set the maximum step size to a very large number for all
 488   // neutral particles and those outside the driftvolume
 489   gMC->SetMaxStep(kBig);
 490
 491   // Use only charged tracks
 492   if (( gMC->TrackCharge()       ) &&
 493       (!gMC->IsTrackDisappeared())) {
 494
 495     // Inside a sensitive volume?
 496     drRegion = kFALSE;
 497     amRegion = kFALSE;
 498     cIdCurrent = gMC->CurrentVolName();
 499     if (cIdSensDr == cIdCurrent[1]) {
 500       drRegion = kTRUE;
 501     }
 502     if (cIdSensAm == cIdCurrent[1]) {
 503       amRegion = kTRUE;
 504     }
 505     if (drRegion || amRegion) {
 506
 507       // The hit coordinates and charge
 508       gMC->TrackPosition(pos);
 509       hits[0] = pos[0];
 510       hits[1] = pos[1];
 511       hits[2] = pos[2];
 512
 513       // The sector number (0 - 17), according to standard coordinate system
 514       cIdPath      = gGeoManager->GetPath();
 515       cIdSector[0] = cIdPath[21];
 516       cIdSector[1] = cIdPath[22];
 517       sec = atoi(cIdSector);
 518
 519       // The plane and chamber number
 520       cIdChamber[0]   = cIdCurrent[2];
 521       cIdChamber[1]   = cIdCurrent[3];
 522       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 523       cha = ((Int_t) idChamber / kNplan);
 524       pla = ((Int_t) idChamber % kNplan);
 525
 526       // The detector number
 527       det = fGeometry->GetDetector(pla,cha,sec);
 528
 529       // Special hits only in the drift region
 530       if      ((drRegion) &&
 531                (gMC->IsTrackEntering())) {
 532
 533         // Create a track reference at the entrance of each
 534         // chamber that contains the momentum components of the particle
 535         gMC->TrackMomentum(mom);
 536         AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 537
 538         // Create the hits from TR photons if electron/positron is
 539         // entering the drift volume
 540         if ((fTR) &&
 541             (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
 542           CreateTRhit(det);
 543         }
 544
 545       }
 546       else if ((amRegion) &&
 547                (gMC->IsTrackExiting())) {
 548
 549         // Create a track reference at the exit of each
 550         // chamber that contains the momentum components of the particle
 551         gMC->TrackMomentum(mom);
 552         AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 553
 554       }
 555
 556       // Calculate the energy of the delta-electrons
 557       // modified by Alex Bercuci (A.Bercuci@gsi.de) on 26.01.06
 558       // take into account correlation with the underlying GEANT tracking
 559       // mechanism. see
 560       // http://www-linux.gsi.de/~abercuci/Contributions/TRD/index.html
 561       //
 562       // determine the most significant process (last on the processes list)
 563       // which caused this hit
 564       gMC->StepProcesses(processes);
 565       Int_t nofprocesses = processes.GetSize();
 566       Int_t pid;
 567       if (!nofprocesses) {
 568         pid = 0;
 569       }
 570       else {
 571         pid =   processes[nofprocesses-1];
 572       }
 573
 574       // Generate Edep according to GEANT parametrisation
 575       eDelta = TMath::Exp(fDeltaG->GetRandom()) - kPoti;
 576       eDelta = TMath::Max(eDelta,0.0);
 577       Float_t prRange = 0.0;
 578       Float_t range   = gMC->TrackLength() - fTrackLength0;
 579       // merge GEANT tracker information with locally cooked one
 580       if (gAlice->GetMCApp()->GetCurrentTrackNumber() == fPrimaryTrackPid) {
 581         if      (pid == 27) {
 582           if (eDelta >= kECut) {
 583             prRange = kRa * eDelta * 0.001
 584                     * (1.0 - kRb / (1.0 + kRc * eDelta * 0.001)) / kRho;
 585             if (prRange >= (3.7 - range)) {
 586               eDelta *= 0.1;
 587             }
 588           }
 589         }
 590         else if (pid ==  1) {
 591           if (eDelta <  kECut) {
 592             eDelta *= 0.5;
 593           }
 594           else {
 595             prRange = kRa * eDelta * 0.001
 596                     * (1.0 - kRb / (1.0 + kRc * eDelta * 0.001)) / kRho;
 597             if (prRange >= ((AliTRDgeometry::DrThick()
 598                            + AliTRDgeometry::AmThick()) - range)) {
 599               eDelta *= 0.05;
 600             }
 601             else {
 602               eDelta *= 0.5;
 603             }
 604           }
 605         }
 606         else {
 607           eDelta = 0.0;
 608         }
 609       }
 610       else {
 611         eDelta = 0.0;
 612       }
 613
 614       // Generate the electron cluster size
 615       if (eDelta > 0.0) {
 616
 617         qTot = ((Int_t) (eDelta / kWion) + 1);
 618
 619         // Create a new dEdx hit
 620         AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 621               ,det
 622               ,hits
 623               ,qTot
 624               ,gMC->TrackTime()*1.0e06
 625               ,drRegion);
 626
 627       }
 628
 629       // Calculate the maximum step size for the next tracking step
 630       // Produce only one hit if Ekin is below cutoff
 631       aMass = gMC->TrackMass();
 632       if ((gMC->Etot() - aMass) > kEkinMinStep) {
 633
 634         // The energy loss according to Bethe Bloch
 635         iPdg = TMath::Abs(gMC->TrackPid());
 636         if ((iPdg != kPdgElectron) ||
 637             ((iPdg == kPdgElectron) &&
 638              (pTot  < kPTotMaxEl))) {
 639           gMC->TrackMomentum(mom);
 640           pTot      = mom.Rho();
 641           betaGamma = pTot / aMass;
 642           pp        = BetheBlochGeant(betaGamma);
 643           // Take charge > 1 into account
 644           charge     = gMC->TrackCharge();
 645           if (TMath::Abs(charge) > 1) {
 646             pp = pp * charge*charge;
 647           }
 648         }
 649         else {
 650           // Electrons above 20 Mev/c are at the plateau
 651           pp = kPrim * kPlateau;
 652         }
 653
 654         Int_t nsteps = 0;
 655         do {
 656           nsteps = gRandom->Poisson(pp);
 657         } while(!nsteps);
 658         stepSize = 1.0 / nsteps;
 659         gMC->SetMaxStep(stepSize);
 660
 661       }
 662
 663     }
 664
 665   }
 666
 667 }
 668
 669 //_____________________________________________________________________________
 670 void AliTRDv1::StepManagerErmilova()
 671 {
 672   //
 673   // Slow simulator. Every charged track produces electron cluster as hits
 674   // along its path across the drift volume. The step size is set acording
 675   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 676   // a spectrum taken from Ermilova et al.
 677   //
 678
 679   Int_t    pla = 0;
 680   Int_t    cha = 0;
 681   Int_t    sec = 0;
 682   Int_t    det = 0;
 683   Int_t    iPdg;
 684   Int_t    qTot;
 685
 686   Float_t  hits[3];
 687   Double_t random[1];
 688   Float_t  charge;
 689   Float_t  aMass;
 690
 691   Double_t pTot     = 0.0;
 692   Double_t eDelta;
 693   Double_t betaGamma;
 694   Double_t pp;
 695   Double_t stepSize;
 696
 697   Bool_t   drRegion = kFALSE;
 698   Bool_t   amRegion = kFALSE;
 699
 700   TString  cIdPath;
 701   Char_t   cIdSector[3];
 702            cIdSector[2]  = 0;
 703
 704   TString  cIdCurrent;
 705   TString  cIdSensDr = "J";
 706   TString  cIdSensAm = "K";
 707   Char_t   cIdChamber[3];
 708            cIdChamber[2] = 0;
 709
 710   TLorentzVector pos;
 711   TLorentzVector mom;
 712
 713   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 714   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 715   const Int_t    kNdetsec     = kNplan * kNcham;
 716
 717   const Double_t kBig         = 1.0e+12; // Infinitely big
 718   const Float_t  kWion        = 23.53;   // Ionization energy
 719   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
 720
 721   // Minimum energy for the step size adjustment
 722   const Float_t  kEkinMinStep = 1.0e-5;
 723
 724   // Plateau value of the energy-loss for electron in xenon
 725   // The averaged value (26/3/99)
 726   const Float_t  kPlateau     = 1.55;
 727   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 728   const Float_t  kPrim        = 48.0;
 729   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 730   const Float_t  kPoti        = 12.1;
 731   // PDG code electron
 732   const Int_t    kPdgElectron = 11;
 733
 734   // Set the maximum step size to a very large number for all
 735   // neutral particles and those outside the driftvolume
 736   gMC->SetMaxStep(kBig);
 737
 738   // Use only charged tracks
 739   if (( gMC->TrackCharge()       ) &&
 740       (!gMC->IsTrackDisappeared())) {
 741
 742     // Inside a sensitive volume?
 743     drRegion = kFALSE;
 744     amRegion = kFALSE;
 745     cIdCurrent = gMC->CurrentVolName();
 746     if (cIdSensDr == cIdCurrent[1]) {
 747       drRegion = kTRUE;
 748     }
 749     if (cIdSensAm == cIdCurrent[1]) {
 750       amRegion = kTRUE;
 751     }
 752     if (drRegion || amRegion) {
 753
 754       // The hit coordinates and charge
 755       gMC->TrackPosition(pos);
 756       hits[0] = pos[0];
 757       hits[1] = pos[1];
 758       hits[2] = pos[2];
 759
 760       // The sector number (0 - 17), according to standard coordinate system
 761       cIdPath      = gGeoManager->GetPath();
 762       cIdSector[0] = cIdPath[21];
 763       cIdSector[1] = cIdPath[22];
 764       sec = atoi(cIdSector);
 765
 766       // The plane and chamber number
 767       cIdChamber[0] = cIdCurrent[2];
 768       cIdChamber[1] = cIdCurrent[3];
 769       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 770       cha = ((Int_t) idChamber / kNplan);
 771       pla = ((Int_t) idChamber % kNplan);
 772
 773       // The detector number
 774       det = fGeometry->GetDetector(pla,cha,sec);
 775
 776       // Special hits only in the drift region
 777       if      ((drRegion) &&
 778                (gMC->IsTrackEntering())) {
 779
 780         // Create a track reference at the entrance of each
 781         // chamber that contains the momentum components of the particle
 782         gMC->TrackMomentum(mom);
 783         AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 784
 785         // Create the hits from TR photons if electron/positron is
 786         // entering the drift volume
 787         if ((fTR) &&
 788             (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
 789           CreateTRhit(det);
 790         }
 791
 792       }
 793       else if ((amRegion) &&
 794                (gMC->IsTrackExiting())) {
 795
 796         // Create a track reference at the exit of each
 797         // chamber that contains the momentum components of the particle
 798         gMC->TrackMomentum(mom);
 799         AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 800
 801       }
 802
 803       // Calculate the energy of the delta-electrons
 804       eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 805       eDelta = TMath::Max(eDelta,0.0);
 806
 807       // Generate the electron cluster size
 808       if (eDelta > 0.0) {
 809
 810         qTot = ((Int_t) (eDelta / kWion) + 1);
 811
 812         // Create a new dEdx hit
 813         if (drRegion) {
 814           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 815                 ,det
 816                 ,hits
 817                 ,qTot
 818                 ,gMC->TrackTime()*1.0e06
 819                 ,kTRUE);
 820         }
 821         else {
 822           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 823                 ,det
 824                 ,hits
 825                 ,qTot
 826                 ,gMC->TrackTime()*1.0e06
 827                 ,kFALSE);
 828         }
 829
 830       }
 831
 832       // Calculate the maximum step size for the next tracking step
 833       // Produce only one hit if Ekin is below cutoff
 834       aMass = gMC->TrackMass();
 835       if ((gMC->Etot() - aMass) > kEkinMinStep) {
 836
 837         // The energy loss according to Bethe Bloch
 838         iPdg  = TMath::Abs(gMC->TrackPid());
 839         if ((iPdg != kPdgElectron) ||
 840             ((iPdg == kPdgElectron) &&
 841              (pTot  < kPTotMaxEl))) {
 842           gMC->TrackMomentum(mom);
 843           pTot      = mom.Rho();
 844           betaGamma = pTot / aMass;
 845           pp        = kPrim * BetheBloch(betaGamma);
 846           // Take charge > 1 into account
 847           charge = gMC->TrackCharge();
 848           if (TMath::Abs(charge) > 1) {
 849             pp = pp * charge*charge;
 850           }
 851         }
 852         else {
 853           // Electrons above 20 Mev/c are at the plateau
 854           pp = kPrim * kPlateau;
 855         }
 856
 857         if (pp > 0.0) {
 858           do {
 859             gMC->GetRandom()->RndmArray(1,random);
 860           }
 861           while ((random[0] == 1.0) ||
 862                  (random[0] == 0.0));
 863           stepSize = - TMath::Log(random[0]) / pp;
 864           gMC->SetMaxStep(stepSize);
 865         }
 866
 867       }
 868
 869     }
 870
 871   }
 872
 873 }
 874
 875 //_____________________________________________________________________________
 876 void AliTRDv1::StepManagerFixedStep()
 877 {
 878   //
 879   // Slow simulator. Every charged track produces electron cluster as hits
 880   // along its path across the drift volume. The step size is fixed in
 881   // this version of the step manager.
 882   //
 883
 884   // PDG code electron
 885   const Int_t   kPdgElectron = 11;
 886
 887   Int_t    pla = 0;
 888   Int_t    cha = 0;
 889   Int_t    sec = 0;
 890   Int_t    det = 0;
 891   Int_t    qTot;
 892
 893   Float_t  hits[3];
 894   Double_t eDep;
 895
 896   Bool_t   drRegion = kFALSE;
 897   Bool_t   amRegion = kFALSE;
 898
 899   TString  cIdPath;
 900   Char_t   cIdSector[3];
 901            cIdSector[2]  = 0;
 902
 903   TString  cIdCurrent;
 904   TString  cIdSensDr = "J";
 905   TString  cIdSensAm = "K";
 906   Char_t   cIdChamber[3];
 907            cIdChamber[2] = 0;
 908
 909   TLorentzVector pos;
 910   TLorentzVector mom;
 911
 912   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 913   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 914   const Int_t    kNdetsec     = kNplan * kNcham;
 915
 916   const Double_t kBig         = 1.0e+12;
 917
 918   const Float_t  kWion        = 23.53;   // Ionization energy
 919   const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment
 920
 921   // Set the maximum step size to a very large number for all
 922   // neutral particles and those outside the driftvolume
 923   gMC->SetMaxStep(kBig);
 924
 925   // If not charged track or already stopped or disappeared, just return.
 926   if ((!gMC->TrackCharge()) ||
 927         gMC->IsTrackDisappeared()) {
 928     return;
 929   }
 930
 931   // Inside a sensitive volume?
 932   cIdCurrent = gMC->CurrentVolName();
 933
 934   if (cIdSensDr == cIdCurrent[1]) {
 935     drRegion = kTRUE;
 936   }
 937   if (cIdSensAm == cIdCurrent[1]) {
 938     amRegion = kTRUE;
 939   }
 940
 941   if ((!drRegion) &&
 942       (!amRegion)) {
 943     return;
 944   }
 945
 946   // The hit coordinates and charge
 947   gMC->TrackPosition(pos);
 948   hits[0] = pos[0];
 949   hits[1] = pos[1];
 950   hits[2] = pos[2];
 951
 952   // The sector number (0 - 17), according to standard coordinate system
 953   cIdPath      = gGeoManager->GetPath();
 954   cIdSector[0] = cIdPath[21];
 955   cIdSector[1] = cIdPath[22];
 956   sec = atoi(cIdSector);
 957
 958   // The plane and chamber number
 959   cIdChamber[0]   = cIdCurrent[2];
 960   cIdChamber[1]   = cIdCurrent[3];
 961   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 962   cha = ((Int_t) idChamber / kNplan);
 963   pla = ((Int_t) idChamber % kNplan);
 964
 965   // The detector number
 966   det = fGeometry->GetDetector(pla,cha,sec);
 967
 968   // 0: InFlight 1:Entering 2:Exiting
 969   Int_t trkStat = 0;
 970
 971   // Special hits only in the drift region
 972   if      ((drRegion) &&
 973            (gMC->IsTrackEntering())) {
 974
 975     // Create a track reference at the entrance of each
 976     // chamber that contains the momentum components of the particle
 977     gMC->TrackMomentum(mom);
 978     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 979     trkStat = 1;
 980
 981     // Create the hits from TR photons if electron/positron is
 982     // entering the drift volume
 983     if ((fTR) &&
 984         (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
 985       CreateTRhit(det);
 986     }
 987
 988   }
 989   else if ((amRegion) &&
 990            (gMC->IsTrackExiting())) {
 991
 992     // Create a track reference at the exit of each
 993     // chamber that contains the momentum components of the particle
 994     gMC->TrackMomentum(mom);
 995     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 996     trkStat = 2;
 997
 998   }
 999
1000   // Calculate the charge according to GEANT Edep
1001   // Create a new dEdx hit
1002   eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
1003   qTot = (Int_t) (eDep / kWion);
1004   if ((qTot) ||
1005       (trkStat)) {
1006     AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
1007           ,det
1008           ,hits
1009           ,qTot
1010           ,gMC->TrackTime()*1.0e06
1011           ,drRegion);
1012   }
1013
1014   // Set Maximum Step Size
1015   // Produce only one hit if Ekin is below cutoff
1016   if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) {
1017     return;
1018   }
1019   gMC->SetMaxStep(fStepSize);
1020
1021 }
1022
1023 //_____________________________________________________________________________
1024 Double_t AliTRDv1::BetheBloch(Double_t bg)
1025 {
1026   //
1027   // Parametrization of the Bethe-Bloch-curve
1028   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
1029   //
1030
1031   // This parameters have been adjusted to averaged values from GEANT
1032   const Double_t kP1    = 7.17960e-02;
1033   const Double_t kP2    = 8.54196;
1034   const Double_t kP3    = 1.38065e-06;
1035   const Double_t kP4    = 5.30972;
1036   const Double_t kP5    = 2.83798;
1037
1038   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
1039   const Double_t kBgMin = 0.8;
1040   const Double_t kBBMax = 6.83298;
1041
1042   if (bg > kBgMin) {
1043     Double_t yy = bg / TMath::Sqrt(1.0 + bg*bg);
1044     Double_t aa = TMath::Power(yy,kP4);
1045     Double_t bb = TMath::Power((1.0/bg),kP5);
1046              bb = TMath::Log(kP3 + bb);
1047     return ((kP2 - aa - bb) * kP1 / aa);
1048   }
1049   else {
1050     return kBBMax;
1051   }
1052
1053 }
1054
1055 //_____________________________________________________________________________
1056 Double_t AliTRDv1::BetheBlochGeant(Double_t bg)
1057 {
1058   //
1059   // Return dN/dx (number of primary collisions per centimeter)
1060   // for given beta*gamma factor.
1061   //
1062   // Implemented by K.Oyama according to GEANT 3 parametrization shown in
1063   // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
1064   // This must be used as a set with IntSpecGeant.
1065   //
1066
1067   Int_t    i = 0;
1068
1069   Double_t arrG[20]  = {     1.100000,     1.200000,     1.300000,     1.500000
1070                        ,     1.800000,     2.000000,     2.500000,     3.000000
1071                        ,     4.000000,     7.000000,    10.000000,    20.000000
1072                        ,    40.000000,    70.000000,   100.000000,   300.000000
1073                        ,   600.000000,  1000.000000,  3000.000000, 10000.000000 };
1074
1075   Double_t arrNC[20] = {    75.009056,    45.508083,    35.299252,    27.116327
1076                        ,    22.734999,    21.411915,    19.934095,    19.449375
1077                        ,    19.344431,    20.185553,    21.027925,    22.912676
1078                        ,    24.933352,    26.504053,    27.387468,    29.566597
1079                        ,    30.353779,    30.787134,    31.129285,    31.157350 };
1080
1081   // Betagamma to gamma
1082   Double_t g = TMath::Sqrt(1.0 + bg*bg);
1083
1084   // Find the index just before the point we need.
1085   for (i = 0; i < 18; i++) {
1086     if ((arrG[i]   < g) &&
1087         (arrG[i+1] > g)) {
1088       break;
1089     }
1090   }
1091
1092   // Simple interpolation.
1093   Double_t pp = ((arrNC[i+1] - arrNC[i]) / (arrG[i+1]  - arrG[i]))
1094               * (g - arrG[i]) + arrNC[i];
1095
1096   return pp;
1097
1098 }
1099
1100 //_____________________________________________________________________________
1101 Double_t Ermilova(Double_t *x, Double_t *)
1102 {
1103   //
1104   // Calculates the delta-ray energy distribution according to Ermilova.
1105   // Logarithmic scale !
1106   //
1107
1108   Double_t energy;
1109   Double_t dpos;
1110   Double_t dnde;
1111
1112   Int_t    pos1;
1113   Int_t    pos2;
1114
1115   const Int_t kNv = 31;
1116
1117   Float_t vxe[kNv] = {  2.3026,  2.9957,  3.4012,  3.6889,  3.9120
1118                      ,  4.0943,  4.2485,  4.3820,  4.4998,  4.6052
1119                      ,  4.7005,  5.0752,  5.2983,  5.7038,  5.9915
1120                      ,  6.2146,  6.5221,  6.9078,  7.3132,  7.6009
1121                      ,  8.0064,  8.5172,  8.6995,  8.9872,  9.2103
1122                      ,  9.4727,  9.9035, 10.3735, 10.5966, 10.8198
1123                      , 11.5129 };
1124
1125   Float_t vye[kNv] = { 80.0,    31.0,    23.3,    21.1,    21.0
1126                      , 20.9,    20.8,    20.0,    16.0,    11.0
1127                      ,  8.0,     6.0,     5.2,     4.6,     4.0
1128                      ,  3.5,     3.0,     1.4,     0.67,    0.44
1129                      ,  0.3,     0.18,    0.12,    0.08,    0.056
1130                      ,  0.04,    0.023,   0.015,   0.011,    0.01
1131                      ,  0.004  };
1132
1133   energy = x[0];
1134
1135   // Find the position
1136   pos1 = 0;
1137   pos2 = 0;
1138   dpos = 0;
1139   do {
1140     dpos = energy - vxe[pos2++];
1141   }
1142   while (dpos > 0);
1143   pos2--;
1144   if (pos2 > kNv) {
1145     pos2 = kNv - 1;
1146   }
1147   pos1 = pos2 - 1;
1148
1149   // Differentiate between the sampling points
1150   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
1151
1152   return dnde;
1153
1154 }
1155
1156 //_____________________________________________________________________________
1157 Double_t IntSpecGeant(Double_t *x, Double_t *)
1158 {
1159   //
1160   // Integrated spectrum from Geant3
1161   //
1162
1163   const Int_t npts = 83;
1164   Double_t arre[npts]    = {  2.421257,     2.483278,    2.534301,     2.592230
1165                            ,  2.672067,     2.813299,    3.015059,     3.216819
1166                            ,  3.418579,     3.620338,    3.868209,     3.920198
1167                            ,  3.978284,     4.063923,    4.186264,     4.308605
1168                            ,  4.430946,     4.553288,    4.724261,     4.837736
1169                            ,  4.999842,     5.161949,    5.324056,     5.486163
1170                            ,  5.679688,     5.752998,    5.857728,     5.962457
1171                            ,  6.067185,     6.171914,    6.315653,     6.393674
1172                            ,  6.471694,     6.539689,    6.597658,     6.655627
1173                            ,  6.710957,     6.763648,    6.816338,     6.876198
1174                            ,  6.943227,     7.010257,    7.106285,     7.252151
1175                            ,  7.460531,     7.668911,    7.877290,     8.085670
1176                            ,  8.302979,     8.353585,    8.413120,     8.483500
1177                            ,  8.541030,     8.592857,    8.668865,     8.820485
1178                            ,  9.037086,     9.253686,    9.470286,     9.686887
1179                            ,  9.930838,     9.994655,   10.085822,    10.176990
1180                            , 10.268158,    10.359325,   10.503614,    10.627565
1181                            , 10.804637,    10.981709,   11.158781,    11.335854
1182                            , 11.593397,    11.781165,   12.049404,    12.317644
1183                            , 12.585884,    12.854123,   14.278421,    16.975889
1184                            , 20.829416,    24.682943,   28.536469 };
1185
1186   Double_t arrdnde[npts] = { 10.960000,    10.960000,   10.359500,     9.811340
1187                            ,  9.1601500,    8.206670,    6.919630,     5.655430
1188                            ,  4.6221300,    3.777610,    3.019560,     2.591950
1189                            ,  2.5414600,    2.712920,    3.327460,     4.928240
1190                            ,  7.6185300,   10.966700,   12.225800,     8.094750
1191                            ,  3.3586900,    1.553650,    1.209600,     1.263840
1192                            ,  1.3241100,    1.312140,    1.255130,     1.165770
1193                            ,  1.0594500,    0.945450,    0.813231,     0.699837
1194                            ,  0.6235580,    2.260990,    2.968350,     2.240320
1195                            ,  1.7988300,    1.553300,    1.432070,     1.535520
1196                            ,  1.4429900,    1.247990,    1.050750,     0.829549
1197                            ,  0.5900280,    0.395897,    0.268741,     0.185320
1198                            ,  0.1292120,    0.103545,    0.0949525,    0.101535
1199                            ,  0.1276380,    0.134216,    0.123816,     0.104557
1200                            ,  0.0751843,    0.0521745,   0.0373546,    0.0275391
1201                            ,  0.0204713,    0.0169234,   0.0154552,    0.0139194
1202                            ,  0.0125592,    0.0113638,   0.0107354,    0.0102137
1203                            ,  0.00845984,   0.00683338,  0.00556836,   0.00456874
1204                            ,  0.0036227,    0.00285991,  0.00226664,   0.00172234
1205                            ,  0.00131226,   0.00100284,  0.000465492,  7.26607e-05
1206                            ,  3.63304e-06,  0.0000000,   0.0000000   };
1207
1208   Int_t    i;
1209   Double_t energy = x[0];
1210
1211   if (energy >= arre[npts-1]) {
1212     return 0.0;
1213   }
1214
1215   for (i = 0; i < npts; i++) {
1216     if (energy < arre[i]) {
1217       break;
1218     }
1219   }
1220
1221   if (i == 0) {
1222     AliErrorGeneral("AliTRDv1::IntSpecGeant","Given energy value is too small or zero");
1223   }
1224
1225   return arrdnde[i];
1226
1227 }