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