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