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