TRD/AliTRDv1.cxx

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