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