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