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