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