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