TRD/AliTRDv1.cxx

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