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