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
  34 #include "AliConst.h"
  35 #include "AliLog.h"
  36 #include "AliMC.h"
  37 #include "AliRun.h"
  38
  39 #include "AliTRDgeometry.h"
  40 #include "AliTRDhit.h"
  41 #include "AliTRDsim.h"
  42 #include "AliTRDv1.h"
  43
  44 ClassImp(AliTRDv1)
  45
  46 //_____________________________________________________________________________
  47 AliTRDv1::AliTRDv1()
  48   :AliTRD()
  49   ,fTRon(kFALSE)
  50   ,fTR(NULL)
  51   ,fTypeOfStepManager(0)
  52   ,fStepSize(0)
  53   ,fDeltaE(NULL)
  54   ,fDeltaG(NULL)
  55   ,fTrackLength0(0)
  56   ,fPrimaryTrackPid(0)
  57 {
  58   //
  59   // Default constructor
  60   //
  61
  62 }
  63
  64 //_____________________________________________________________________________
  65 AliTRDv1::AliTRDv1(const char *name, const char *title)
  66   :AliTRD(name,title)
  67   ,fTRon(kTRUE)
  68   ,fTR(NULL)
  69   ,fTypeOfStepManager(2)
  70   ,fStepSize(0.1)
  71   ,fDeltaE(NULL)
  72   ,fDeltaG(NULL)
  73   ,fTrackLength0(0)
  74   ,fPrimaryTrackPid(0)
  75 {
  76   //
  77   // Standard constructor for Transition Radiation Detector version 1
  78   //
  79
  80   SetBufferSize(128000);
  81
  82 }
  83
  84 //_____________________________________________________________________________
  85 AliTRDv1::AliTRDv1(const AliTRDv1 &trd)
  86   :AliTRD(trd)
  87   ,fTRon(trd.fTRon)
  88   ,fTR(NULL)
  89   ,fTypeOfStepManager(trd.fTypeOfStepManager)
  90   ,fStepSize(trd.fStepSize)
  91   ,fDeltaE(NULL)
  92   ,fDeltaG(NULL)
  93   ,fTrackLength0(trd.fTrackLength0)
  94   ,fPrimaryTrackPid(trd.fPrimaryTrackPid)
  95 {
  96   //
  97   // Copy constructor
  98   //
  99
 100   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
 101   fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
 102   fTR->Copy(*((AliTRDv1 &) trd).fTR);
 103
 104 }
 105
 106 //_____________________________________________________________________________
 107 AliTRDv1::~AliTRDv1()
 108 {
 109   //
 110   // AliTRDv1 destructor
 111   //
 112
 113   if (fDeltaE) {
 114     delete fDeltaE;
 115     fDeltaE = 0;
 116   }
 117
 118   if (fDeltaG) {
 119     delete fDeltaG;
 120     fDeltaG = 0;
 121   }
 122
 123   if (fTR) {
 124     delete fTR;
 125     fTR     = 0;
 126   }
 127
 128 }
 129
 130 //_____________________________________________________________________________
 131 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
 132 {
 133   //
 134   // Assignment operator
 135   //
 136
 137   if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
 138
 139   return *this;
 140
 141 }
 142
 143 //_____________________________________________________________________________
 144 void AliTRDv1::Copy(TObject &trd) const
 145 {
 146   //
 147   // Copy function
 148   //
 149
 150   ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager;
 151   ((AliTRDv1 &) trd).fStepSize          = fStepSize;
 152   ((AliTRDv1 &) trd).fTRon              = fTRon;
 153   ((AliTRDv1 &) trd).fTrackLength0      = fTrackLength0;
 154   ((AliTRDv1 &) trd).fPrimaryTrackPid   = fPrimaryTrackPid;
 155
 156   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
 157   fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
 158   fTR->Copy(*((AliTRDv1 &) trd).fTR);
 159
 160 }
 161
 162 //_____________________________________________________________________________
 163 void AliTRDv1::AddAlignableVolumes() const
 164 {
 165   //
 166   // Create entries for alignable volumes associating the symbolic volume
 167   // name with the corresponding volume path. Needs to be syncronized with
 168   // eventual changes in the geometry.
 169   //
 170
 171   TString volPath;
 172   TString symName;
 173
 174   TString vpStr  = "ALIC_1/B077_1/BSEGMO";
 175   TString vpApp1 = "_1/BTRD";
 176   TString vpApp2 = "_1";
 177   TString vpApp3 = "/UTR1_1/UTS1_1/UTI1_1/UT";
 178
 179   TString snStr  = "TRD/sm";
 180   TString snApp1 = "/st";
 181   TString snApp2 = "/pl";
 182
 183   //
 184   // The super modules
 185   // The symbolic names are: TRD/sm00
 186   //                           ...
 187   //                         TRD/sm17
 188   //
 189   for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
 190
 191     volPath  = vpStr;
 192     volPath += isect;
 193     volPath += vpApp1;
 194     volPath += isect;
 195     volPath += vpApp2;
 196
 197     symName  = snStr;
 198     symName += Form("%02d",isect);
 199
 200     gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
 201
 202   }
 203
 204   //
 205   // The readout chambers
 206   // The symbolic names are: TRD/sm00/st0/pl0
 207   //                           ...
 208   //                         TRD/sm17/st4/pl5
 209   //
 210   for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
 211     for (Int_t icham = 0; icham < AliTRDgeometry::Ncham(); icham++) {
 212       for (Int_t iplan = 0; iplan < AliTRDgeometry::Nplan(); iplan++) {
 213
 214         Int_t idet = AliTRDgeometry::GetDetectorSec(iplan,icham);
 215
 216         volPath  = vpStr;
 217         volPath += isect;
 218         volPath += vpApp1;
 219         volPath += isect;
 220         volPath += vpApp2;
 221         volPath += vpApp3;
 222         volPath += Form("%02d",idet);
 223         volPath += vpApp2;
 224
 225         symName  = snStr;
 226         symName += Form("%02d",isect);
 227         symName += snApp1;
 228         symName += icham;
 229         symName += snApp2;
 230         symName += iplan;
 231
 232         printf("volPath=%s\n",volPath.Data());
 233         printf("symName=%s\n",symName.Data());
 234
 235         gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
 236
 237       }
 238     }
 239   }
 240
 241 }
 242
 243 //_____________________________________________________________________________
 244 void AliTRDv1::CreateGeometry()
 245 {
 246   //
 247   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
 248   // This version covers the full azimuth.
 249   //
 250
 251   // Check that FRAME is there otherwise we have no place where to put the TRD
 252   AliModule* frame = gAlice->GetModule("FRAME");
 253   if (!frame) {
 254     AliError("TRD needs FRAME to be present\n");
 255     return;
 256   }
 257
 258   // Define the chambers
 259   AliTRD::CreateGeometry();
 260
 261 }
 262
 263 //_____________________________________________________________________________
 264 void AliTRDv1::CreateMaterials()
 265 {
 266   //
 267   // Create materials for the Transition Radiation Detector version 1
 268   //
 269
 270   AliTRD::CreateMaterials();
 271
 272 }
 273
 274 //_____________________________________________________________________________
 275 void AliTRDv1::CreateTRhit(Int_t det)
 276 {
 277   //
 278   // Creates an electron cluster from a TR photon.
 279   // The photon is assumed to be created a the end of the radiator. The
 280   // distance after which it deposits its energy takes into account the
 281   // absorbtion of the entrance window and of the gas mixture in drift
 282   // volume.
 283   //
 284
 285   // PDG code electron
 286   const Int_t   kPdgElectron = 11;
 287
 288   // Ionization energy
 289   const Float_t kWion        = 23.53;
 290
 291   // Maximum number of TR photons per track
 292   const Int_t   kNTR         = 50;
 293
 294   TLorentzVector mom;
 295   TLorentzVector pos;
 296
 297   // Create TR at the entrance of the chamber
 298   if (gMC->IsTrackEntering()) {
 299
 300     // Create TR only for electrons
 301     Int_t iPdg = gMC->TrackPid();
 302     if (TMath::Abs(iPdg) != kPdgElectron) {
 303       return;
 304     }
 305
 306     Float_t eTR[kNTR];
 307     Int_t   nTR;
 308
 309     // Create TR photons
 310     gMC->TrackMomentum(mom);
 311     Float_t pTot = mom.Rho();
 312     fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
 313     if (nTR > kNTR) {
 314       AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
 315     }
 316
 317     // Loop through the TR photons
 318     for (Int_t iTR = 0; iTR < nTR; iTR++) {
 319
 320       Float_t energyMeV = eTR[iTR] * 0.001;
 321       Float_t energyeV  = eTR[iTR] * 1000.0;
 322       Float_t absLength = 0.0;
 323       Float_t sigma     = 0.0;
 324
 325       // Take the absorbtion in the entrance window into account
 326       Double_t muMy = fTR->GetMuMy(energyMeV);
 327       sigma         = muMy * fFoilDensity;
 328       if (sigma > 0.0) {
 329         absLength = gRandom->Exp(1.0/sigma);
 330         if (absLength < AliTRDgeometry::MyThick()) {
 331           continue;
 332         }
 333       }
 334       else {
 335         continue;
 336       }
 337
 338       // The absorbtion cross sections in the drift gas
 339       // Gas-mixture (Xe/CO2)
 340       Double_t muXe = fTR->GetMuXe(energyMeV);
 341       Double_t muCO = fTR->GetMuCO(energyMeV);
 342       sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
 343
 344       // The distance after which the energy of the TR photon
 345       // is deposited.
 346       if (sigma > 0.0) {
 347         absLength = gRandom->Exp(1.0/sigma);
 348         if (absLength > (AliTRDgeometry::DrThick()
 349                        + AliTRDgeometry::AmThick())) {
 350           continue;
 351         }
 352       }
 353       else {
 354         continue;
 355       }
 356
 357       // The position of the absorbtion
 358       Float_t posHit[3];
 359       gMC->TrackPosition(pos);
 360       posHit[0] = pos[0] + mom[0] / pTot * absLength;
 361       posHit[1] = pos[1] + mom[1] / pTot * absLength;
 362       posHit[2] = pos[2] + mom[2] / pTot * absLength;
 363
 364       // Create the charge
 365       Int_t q = ((Int_t) (energyeV / kWion));
 366
 367       // Add the hit to the array. TR photon hits are marked
 368       // by negative charge
 369       AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 370             ,det
 371             ,posHit
 372             ,-q
 373             ,kTRUE);
 374
 375     }
 376
 377   }
 378
 379 }
 380
 381 //_____________________________________________________________________________
 382 void AliTRDv1::Init()
 383 {
 384   //
 385   // Initialise Transition Radiation Detector after geometry has been built.
 386   //
 387
 388   AliTRD::Init();
 389
 390   AliDebug(1,"Slow simulator\n");
 391
 392   // Switch on TR simulation as default
 393   if (!fTRon) {
 394     AliInfo("TR simulation off");
 395   }
 396   else {
 397     fTR = new AliTRDsim();
 398   }
 399
 400   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 401   const Float_t kPoti = 12.1;
 402   // Maximum energy (50 keV);
 403   const Float_t kEend = 50000.0;
 404   // Ermilova distribution for the delta-ray spectrum
 405   Float_t poti        = TMath::Log(kPoti);
 406   Float_t eEnd        = TMath::Log(kEend);
 407
 408   // Ermilova distribution for the delta-ray spectrum
 409   fDeltaE = new TF1("deltae" ,Ermilova ,poti,eEnd,0);
 410
 411   // Geant3 distribution for the delta-ray spectrum
 412   fDeltaG = new TF1("deltag",IntSpecGeant,2.421257,28.536469,0);
 413
 414   AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
 415
 416 }
 417
 418 //_____________________________________________________________________________
 419 void AliTRDv1::StepManager()
 420 {
 421   //
 422   // Slow simulator. Every charged track produces electron cluster as hits
 423   // along its path across the drift volume.
 424   //
 425
 426   switch (fTypeOfStepManager) {
 427    case 0:
 428     StepManagerErmilova();
 429     break;
 430    case 1:
 431     StepManagerGeant();
 432     break;
 433    case 2:
 434     StepManagerFixedStep();
 435     break;
 436    default:
 437     AliWarning("Not a valid Step Manager.");
 438   }
 439
 440 }
 441
 442 //_____________________________________________________________________________
 443 void AliTRDv1::SelectStepManager(Int_t t)
 444 {
 445   //
 446   // Selects a step manager type:
 447   //   0 - Ermilova
 448   //   1 - Geant3
 449   //   2 - Fixed step size
 450   //
 451
 452   fTypeOfStepManager = t;
 453   AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
 454
 455 }
 456
 457 //_____________________________________________________________________________
 458 void AliTRDv1::StepManagerGeant()
 459 {
 460   //
 461   // Slow simulator. Every charged track produces electron cluster as hits
 462   // along its path across the drift volume. The step size is set acording
 463   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 464   // a spectrum taken from Geant3.
 465   //
 466   // Version by A. Bercuci
 467   //
 468
 469   Int_t    pla = 0;
 470   Int_t    cha = 0;
 471   Int_t    sec = 0;
 472   Int_t    det = 0;
 473   Int_t    iPdg;
 474   Int_t    qTot;
 475
 476   Float_t  hits[3];
 477   Float_t  charge;
 478   Float_t  aMass;
 479
 480   Double_t pTot     = 0;
 481   Double_t eDelta;
 482   Double_t betaGamma;
 483   Double_t pp;
 484   Double_t stepSize = 0;
 485
 486   Bool_t   drRegion = kFALSE;
 487   Bool_t   amRegion = kFALSE;
 488
 489   TString  cIdCurrent;
 490   TString  cIdSensDr = "J";
 491   TString  cIdSensAm = "K";
 492   Char_t   cIdChamber[3];
 493            cIdChamber[2] = 0;
 494
 495   TLorentzVector pos;
 496   TLorentzVector mom;
 497
 498   TArrayI        processes;
 499
 500   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 501   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 502   const Int_t    kNdetsec     = kNplan * kNcham;
 503
 504   const Double_t kBig         = 1.0e+12; // Infinitely big
 505   const Float_t  kWion        = 23.53;   // Ionization energy
 506   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
 507
 508   // Minimum energy for the step size adjustment
 509   const Float_t  kEkinMinStep = 1.0e-5;
 510   // energy threshold for production of delta electrons
 511   const Float_t  kECut        = 1.0e4;
 512   // Parameters entering the parametrized range for delta electrons
 513   const Float_t  kRa          = 5.37e-4;
 514   const Float_t  kRb          = 0.9815;
 515   const Float_t  kRc          = 3.123e-3;
 516   // Gas density -> To be made user adjustable !
 517   // [0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
 518   const Float_t  kRho         = 0.004945 ;
 519
 520   // Plateau value of the energy-loss for electron in xenon
 521   // The averaged value (26/3/99)
 522   const Float_t  kPlateau     = 1.55;
 523   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 524   const Float_t  kPrim        = 19.34;
 525   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 526   const Float_t  kPoti        = 12.1;
 527   // PDG code electron
 528   const Int_t    kPdgElectron = 11;
 529
 530   // Set the maximum step size to a very large number for all
 531   // neutral particles and those outside the driftvolume
 532   gMC->SetMaxStep(kBig);
 533
 534   // Use only charged tracks
 535   if (( gMC->TrackCharge()       ) &&
 536       (!gMC->IsTrackDisappeared())) {
 537
 538     // Inside a sensitive volume?
 539     drRegion = kFALSE;
 540     amRegion = kFALSE;
 541     cIdCurrent = gMC->CurrentVolName();
 542     if (cIdSensDr == cIdCurrent[1]) {
 543       drRegion = kTRUE;
 544     }
 545     if (cIdSensAm == cIdCurrent[1]) {
 546       amRegion = kTRUE;
 547     }
 548     if (drRegion || amRegion) {
 549
 550       // The hit coordinates and charge
 551       gMC->TrackPosition(pos);
 552       hits[0] = pos[0];
 553       hits[1] = pos[1];
 554       hits[2] = pos[2];
 555
 556       // The sector number (0 - 17)
 557       // The numbering goes clockwise and starts at y = 0
 558       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 559       if (phi < 90.0) {
 560         phi = phi + 270.0;
 561       }
 562       else {
 563         phi = phi -  90.0;
 564       }
 565       sec = ((Int_t) (phi / 20.0));
 566
 567       // The plane and chamber number
 568       cIdChamber[0]   = cIdCurrent[2];
 569       cIdChamber[1]   = cIdCurrent[3];
 570       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 571       cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
 572       pla = ((Int_t) idChamber % kNplan);
 573
 574       // Check on selected volumes
 575       Int_t addthishit = 1;
 576
 577       // Add this hit
 578       if (addthishit) {
 579
 580         // The detector number
 581         det = fGeometry->GetDetector(pla,cha,sec);
 582
 583         // Special hits only in the drift region
 584         if (drRegion) {
 585
 586           // Create a track reference at the entrance and
 587           // exit of each chamber that contain the
 588           // momentum components of the particle
 589           if (gMC->IsTrackEntering() ||
 590               gMC->IsTrackExiting()) {
 591             gMC->TrackMomentum(mom);
 592             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 593           }
 594
 595           if (gMC->IsTrackEntering() &&
 596               !gMC->IsNewTrack()) {
 597             // determine if hit belong to primary track
 598             fPrimaryTrackPid = gAlice->GetMCApp()->GetCurrentTrackNumber();
 599             // determine track length when entering the detector
 600             fTrackLength0    = gMC->TrackLength();
 601           }
 602
 603           // Create the hits from TR photons
 604           if (fTR) CreateTRhit(det);
 605
 606         }
 607
 608         // Calculate the energy of the delta-electrons
 609         // modified by Alex Bercuci (A.Bercuci@gsi.de) on 26.01.06
 610         // take into account correlation with the underlying GEANT tracking
 611         // mechanism. see
 612         // http://www-linux.gsi.de/~abercuci/Contributions/TRD/index.html
 613         //
 614         // determine the most significant process (last on the processes list)
 615         // which caused this hit
 616         gMC->StepProcesses(processes);
 617         Int_t nofprocesses = processes.GetSize();
 618         Int_t pid;
 619         if (!nofprocesses) {
 620           pid = 0;
 621         }
 622         else {
 623           pid = processes[nofprocesses-1];
 624         }
 625
 626         // generate Edep according to GEANT parametrisation
 627         eDelta = TMath::Exp(fDeltaG->GetRandom()) - kPoti;
 628         eDelta = TMath::Max(eDelta,0.0);
 629         Float_t prRange = 0.0;
 630         Float_t range   = gMC->TrackLength() - fTrackLength0;
 631         // merge GEANT tracker information with locally cooked one
 632         if (gAlice->GetMCApp()->GetCurrentTrackNumber() == fPrimaryTrackPid) {
 633           if      (pid == 27) {
 634             if (eDelta >= kECut) {
 635               prRange = kRa * eDelta * 0.001
 636                       * (1.0 - kRb / (1.0 + kRc * eDelta * 0.001)) / kRho;
 637               if (prRange >= (3.7 - range)) {
 638                 eDelta *= 0.1;
 639               }
 640             }
 641           }
 642           else if (pid ==  1) {
 643             if (eDelta <  kECut) {
 644               eDelta *= 0.5;
 645             }
 646             else {
 647               prRange = kRa * eDelta * 0.001
 648                       * (1.0 - kRb / (1.0 + kRc * eDelta * 0.001)) / kRho;
 649               if (prRange >= ((AliTRDgeometry::DrThick()
 650                              + AliTRDgeometry::AmThick()) - range)) {
 651                 eDelta *= 0.05;
 652               }
 653               else {
 654                 eDelta *= 0.5;
 655               }
 656             }
 657           }
 658           else {
 659             eDelta = 0.0;
 660           }
 661         }
 662         else {
 663           eDelta = 0.0;
 664         }
 665
 666         // Generate the electron cluster size
 667         if (eDelta == 0.0) {
 668           qTot = 0;
 669         }
 670         else {
 671           qTot = ((Int_t) (eDelta / kWion) + 1);
 672         }
 673
 674         // Create a new dEdx hit
 675         AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 676               ,det
 677               ,hits
 678               ,qTot
 679               ,drRegion);
 680
 681         // Calculate the maximum step size for the next tracking step
 682         // Produce only one hit if Ekin is below cutoff
 683         aMass = gMC->TrackMass();
 684         if ((gMC->Etot() - aMass) > kEkinMinStep) {
 685
 686           // The energy loss according to Bethe Bloch
 687           iPdg = TMath::Abs(gMC->TrackPid());
 688           if ((iPdg != kPdgElectron) ||
 689               ((iPdg == kPdgElectron) &&
 690                (pTot  < kPTotMaxEl))) {
 691             gMC->TrackMomentum(mom);
 692             pTot      = mom.Rho();
 693             betaGamma = pTot / aMass;
 694             pp        = BetheBlochGeant(betaGamma);
 695             // Take charge > 1 into account
 696             charge     = gMC->TrackCharge();
 697             if (TMath::Abs(charge) > 1) {
 698               pp = pp * charge*charge;
 699             }
 700           }
 701           else {
 702             // Electrons above 20 Mev/c are at the plateau
 703             pp = kPrim * kPlateau;
 704           }
 705
 706           Int_t nsteps = 0;
 707           do {
 708             nsteps = gRandom->Poisson(pp);
 709           } while(!nsteps);
 710           stepSize = 1.0 / nsteps;
 711           gMC->SetMaxStep(stepSize);
 712
 713         }
 714
 715       }
 716
 717     }
 718
 719   }
 720
 721 }
 722
 723 //_____________________________________________________________________________
 724 void AliTRDv1::StepManagerErmilova()
 725 {
 726   //
 727   // Slow simulator. Every charged track produces electron cluster as hits
 728   // along its path across the drift volume. The step size is set acording
 729   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 730   // a spectrum taken from Ermilova et al.
 731   //
 732
 733   Int_t    pla = 0;
 734   Int_t    cha = 0;
 735   Int_t    sec = 0;
 736   Int_t    det = 0;
 737   Int_t    iPdg;
 738   Int_t    qTot;
 739
 740   Float_t  hits[3];
 741   Double_t random[1];
 742   Float_t  charge;
 743   Float_t  aMass;
 744
 745   Double_t pTot     = 0.0;
 746   Double_t eDelta;
 747   Double_t betaGamma;
 748   Double_t pp;
 749   Double_t stepSize;
 750
 751   Bool_t   drRegion = kFALSE;
 752   Bool_t   amRegion = kFALSE;
 753
 754   TString  cIdCurrent;
 755   TString  cIdSensDr = "J";
 756   TString  cIdSensAm = "K";
 757   Char_t   cIdChamber[3];
 758            cIdChamber[2] = 0;
 759
 760   TLorentzVector pos;
 761   TLorentzVector mom;
 762
 763   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 764   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 765   const Int_t    kNdetsec     = kNplan * kNcham;
 766
 767   const Double_t kBig         = 1.0e+12; // Infinitely big
 768   const Float_t  kWion        = 23.53;   // Ionization energy
 769   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
 770
 771   // Minimum energy for the step size adjustment
 772   const Float_t  kEkinMinStep = 1.0e-5;
 773
 774   // Plateau value of the energy-loss for electron in xenon
 775   // The averaged value (26/3/99)
 776   const Float_t  kPlateau     = 1.55;
 777   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 778   const Float_t  kPrim        = 48.0;
 779   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 780   const Float_t  kPoti        = 12.1;
 781   // PDG code electron
 782   const Int_t    kPdgElectron = 11;
 783
 784   // Set the maximum step size to a very large number for all
 785   // neutral particles and those outside the driftvolume
 786   gMC->SetMaxStep(kBig);
 787
 788   // Use only charged tracks
 789   if (( gMC->TrackCharge()       ) &&
 790       (!gMC->IsTrackDisappeared())) {
 791
 792     // Inside a sensitive volume?
 793     drRegion = kFALSE;
 794     amRegion = kFALSE;
 795     cIdCurrent = gMC->CurrentVolName();
 796     if (cIdSensDr == cIdCurrent[1]) {
 797       drRegion = kTRUE;
 798     }
 799     if (cIdSensAm == cIdCurrent[1]) {
 800       amRegion = kTRUE;
 801     }
 802     if (drRegion || amRegion) {
 803
 804       // The hit coordinates and charge
 805       gMC->TrackPosition(pos);
 806       hits[0] = pos[0];
 807       hits[1] = pos[1];
 808       hits[2] = pos[2];
 809
 810       // The sector number (0 - 17)
 811       // The numbering goes clockwise and starts at y = 0
 812       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 813       if (phi < 90.0) {
 814         phi = phi + 270.0;
 815       }
 816       else {
 817         phi = phi -  90.0;
 818       }
 819       sec = ((Int_t) (phi / 20.0));
 820
 821       // The plane and chamber number
 822       cIdChamber[0] = cIdCurrent[2];
 823       cIdChamber[1] = cIdCurrent[3];
 824       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 825       cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
 826       pla = ((Int_t) idChamber % kNplan);
 827
 828       // Check on selected volumes
 829       Int_t addthishit = 1;
 830
 831       // Add this hit
 832       if (addthishit) {
 833
 834         // The detector number
 835         det = fGeometry->GetDetector(pla,cha,sec);
 836
 837         // Special hits only in the drift region
 838         if (drRegion) {
 839
 840           // Create a track reference at the entrance and
 841           // exit of each chamber that contain the
 842           // momentum components of the particle
 843           if (gMC->IsTrackEntering() ||
 844               gMC->IsTrackExiting()) {
 845             gMC->TrackMomentum(mom);
 846             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 847           }
 848           // Create the hits from TR photons
 849           if (fTR) {
 850             CreateTRhit(det);
 851           }
 852
 853         }
 854
 855         // Calculate the energy of the delta-electrons
 856         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 857         eDelta = TMath::Max(eDelta,0.0);
 858         // Generate the electron cluster size
 859         if (eDelta == 0.0) {
 860           qTot = 0;
 861         }
 862         else {
 863           qTot = ((Int_t) (eDelta / kWion) + 1);
 864         }
 865
 866         // Create a new dEdx hit
 867         if (drRegion) {
 868           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 869                 ,det
 870                 ,hits
 871                 ,qTot
 872                 ,kTRUE);
 873         }
 874         else {
 875           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 876                 ,det
 877                 ,hits
 878                 ,qTot
 879                 ,kFALSE);
 880         }
 881
 882         // Calculate the maximum step size for the next tracking step
 883         // Produce only one hit if Ekin is below cutoff
 884         aMass = gMC->TrackMass();
 885         if ((gMC->Etot() - aMass) > kEkinMinStep) {
 886
 887           // The energy loss according to Bethe Bloch
 888           iPdg  = TMath::Abs(gMC->TrackPid());
 889           if ((iPdg != kPdgElectron) ||
 890               ((iPdg == kPdgElectron) &&
 891                (pTot  < kPTotMaxEl))) {
 892             gMC->TrackMomentum(mom);
 893             pTot      = mom.Rho();
 894             betaGamma = pTot / aMass;
 895             pp        = kPrim * BetheBloch(betaGamma);
 896             // Take charge > 1 into account
 897             charge = gMC->TrackCharge();
 898             if (TMath::Abs(charge) > 1) {
 899               pp = pp * charge*charge;
 900             }
 901           }
 902           else {
 903             // Electrons above 20 Mev/c are at the plateau
 904             pp = kPrim * kPlateau;
 905           }
 906
 907           if (pp > 0.0) {
 908             do {
 909               gMC->GetRandom()->RndmArray(1,random);
 910             }
 911             while ((random[0] == 1.0) ||
 912                    (random[0] == 0.0));
 913             stepSize = - TMath::Log(random[0]) / pp;
 914             gMC->SetMaxStep(stepSize);
 915           }
 916
 917         }
 918
 919       }
 920
 921     }
 922
 923   }
 924
 925 }
 926
 927 //_____________________________________________________________________________
 928 void AliTRDv1::StepManagerFixedStep()
 929 {
 930   //
 931   // Slow simulator. Every charged track produces electron cluster as hits
 932   // along its path across the drift volume. The step size is fixed in
 933   // this version of the step manager.
 934   //
 935
 936   Int_t    pla = 0;
 937   Int_t    cha = 0;
 938   Int_t    sec = 0;
 939   Int_t    det = 0;
 940   Int_t    qTot;
 941
 942   Float_t  hits[3];
 943   Double_t eDep;
 944
 945   Bool_t   drRegion = kFALSE;
 946   Bool_t   amRegion = kFALSE;
 947
 948   TString  cIdCurrent;
 949   TString  cIdSensDr = "J";
 950   TString  cIdSensAm = "K";
 951   Char_t   cIdChamber[3];
 952   cIdChamber[2] = 0;
 953
 954   TLorentzVector pos;
 955   TLorentzVector mom;
 956
 957   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 958   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 959   const Int_t    kNdetsec     = kNplan * kNcham;
 960
 961   const Double_t kBig         = 1.0e+12;
 962
 963   const Float_t  kWion        = 23.53;   // Ionization energy
 964   const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment
 965
 966   // Set the maximum step size to a very large number for all
 967   // neutral particles and those outside the driftvolume
 968   gMC->SetMaxStep(kBig);
 969
 970   // If not charged track or already stopped or disappeared, just return.
 971   if ((!gMC->TrackCharge()) ||
 972         gMC->IsTrackDisappeared()) return;
 973
 974   // Inside a sensitive volume?
 975   cIdCurrent = gMC->CurrentVolName();
 976
 977   if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
 978   if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
 979
 980   if ((!drRegion) &&
 981       (!amRegion)) {
 982     return;
 983   }
 984
 985   // The hit coordinates and charge
 986   gMC->TrackPosition(pos);
 987   hits[0] = pos[0];
 988   hits[1] = pos[1];
 989   hits[2] = pos[2];
 990
 991   // The sector number (0 - 17)
 992   // The numbering goes clockwise and starts at y = 0
 993   Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 994   if (phi < 90.0) {
 995     phi = phi + 270.0;
 996   }
 997   else {
 998     phi = phi -  90.0;
 999   }
1000   sec = ((Int_t) (phi / 20.0));
1001
1002   // The plane and chamber number
1003   cIdChamber[0]   = cIdCurrent[2];
1004   cIdChamber[1]   = cIdCurrent[3];
1005   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
1006   cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
1007   pla = ((Int_t) idChamber % kNplan);
1008
1009   // Check on selected volumes
1010   Int_t addthishit = 1;
1011
1012   if (!addthishit) {
1013     return;
1014   }
1015
1016   // The detector number
1017   det = fGeometry->GetDetector(pla,cha,sec);
1018
1019   // 0: InFlight 1:Entering 2:Exiting
1020   Int_t trkStat = 0;
1021
1022   // Special hits only in the drift region
1023   if (drRegion) {
1024
1025     // Create a track reference at the entrance and exit of each
1026     // chamber that contain the momentum components of the particle
1027
1028     if (gMC->IsTrackEntering()) {
1029       gMC->TrackMomentum(mom);
1030       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1031       trkStat = 1;
1032     }
1033     if (gMC->IsTrackExiting()) {
1034       gMC->TrackMomentum(mom);
1035       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1036       trkStat = 2;
1037     }
1038
1039     // Create the hits from TR photons
1040     if (fTR) {
1041       CreateTRhit(det);
1042     }
1043
1044   }
1045
1046   // Calculate the charge according to GEANT Edep
1047   // Create a new dEdx hit
1048   eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
1049   qTot = (Int_t) (eDep / kWion);
1050   AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
1051         ,det
1052         ,hits
1053         ,qTot
1054         ,drRegion);
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   for (i = 0; i < npts; i++) {
1254     if (energy < arre[i]) {
1255       break;
1256     }
1257   }
1258
1259   if (i == 0) {
1260     AliErrorGeneral("AliTRDv1::IntSpecGeant","Given energy value is too small or zero");
1261   }
1262
1263   return arrdnde[i];
1264
1265 }