TRD/AliTRDv2.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 /*
  17 $Log$
  18 Revision 1.13  1999/09/29 09:24:35  fca
  19 Introduction of the Copyright and cvs Log
  20
  21 */
  22
  23 ///////////////////////////////////////////////////////////////////////////////
  24 //                                                                           //
  25 //  Transition Radiation Detector version 2 -- detailed simulation           //
  26 //                                                                           //
  27 //Begin_Html
  28 /*
  29 <img src="picts/AliTRDv2Class.gif">
  30 */
  31 //End_Html
  32 //                                                                           //
  33 //                                                                           //
  34 ///////////////////////////////////////////////////////////////////////////////
  35
  36 #include <TMath.h>
  37 #include <TVector.h>
  38 #include <TRandom.h>
  39
  40 #include "AliTRDv2.h"
  41 #include "AliTRDmatrix.h"
  42 #include "AliRun.h"
  43 #include "AliMC.h"
  44 #include "AliConst.h"
  45
  46 ClassImp(AliTRDv2)
  47
  48 //_____________________________________________________________________________
  49 AliTRDv2::AliTRDv2(const char *name, const char *title)
  50          :AliTRD(name, title)
  51 {
  52   //
  53   // Standard constructor for Transition Radiation Detector version 2
  54   //
  55
  56   fIdSens       = 0;
  57
  58   fIdSpace1     = 0;
  59   fIdSpace2     = 0;
  60   fIdSpace3     = 0;
  61
  62   fIdChamber1   = 0;
  63   fIdChamber2   = 0;
  64   fIdChamber3   = 0;
  65
  66   fSensSelect   = 0;
  67   fSensPlane    = 0;
  68   fSensChamber  = 0;
  69   fSensSector   = 0;
  70
  71   Int_t iplan;
  72
  73   for (iplan = 0; iplan < kNplan; iplan++) {
  74     for (Int_t icham = 0; icham < kNcham; icham++) {
  75       fRowMax[iplan][icham] = 0;
  76     }
  77     fColMax[iplan] = 0;
  78   }
  79   fTimeMax      = 0;
  80
  81   fRowPadSize   = 0;
  82   fColPadSize   = 0;
  83   fTimeBinSize  = 0;
  84
  85   fGasGain      = 0;
  86   fNoise        = 0;
  87   fChipGain     = 0;
  88   fADCoutRange  = 0;
  89   fADCinRange   = 0;
  90   fADCthreshold = 0;
  91
  92   fDiffusionT   = 0;
  93   fDiffusionL   = 0;
  94
  95   fDeltaE       = NULL;
  96
  97   SetBufferSize(128000);
  98
  99 }
 100
 101 //_____________________________________________________________________________
 102 AliTRDv2::~AliTRDv2()
 103 {
 104
 105   if (fDeltaE) delete fDeltaE;
 106
 107 }
 108
 109 //_____________________________________________________________________________
 110 void AliTRDv2::CreateGeometry()
 111 {
 112   //
 113   // Create the GEANT geometry for the Transition Radiation Detector - Version 2
 114   // This version covers the full azimuth.
 115   //
 116   // Author:  Christoph Blume (C.Blume@gsi.de) 20/07/99
 117   //
 118
 119   Float_t xpos, ypos, zpos;
 120
 121   // Check that FRAME is there otherwise we have no place where to put the TRD
 122   AliModule* FRAME = gAlice->GetModule("FRAME");
 123   if (!FRAME) return;
 124
 125   // Define the chambers
 126   AliTRD::CreateGeometry();
 127
 128   // Position the the TRD-sectors in all TRD-volumes in the spaceframe
 129   xpos     = 0.;
 130   ypos     = 0.;
 131   zpos     = 0.;
 132   gMC->Gspos("TRD ",1,"BTR1",xpos,ypos,zpos,0,"ONLY");
 133   gMC->Gspos("TRD ",2,"BTR2",xpos,ypos,zpos,0,"ONLY");
 134   gMC->Gspos("TRD ",3,"BTR3",xpos,ypos,zpos,0,"ONLY");
 135
 136 }
 137
 138 //_____________________________________________________________________________
 139 void AliTRDv2::CreateMaterials()
 140 {
 141   //
 142   // Create materials for the Transition Radiation Detector version 2
 143   //
 144
 145   AliTRD::CreateMaterials();
 146
 147 }
 148
 149 //_____________________________________________________________________________
 150 void AliTRDv2::Diffusion(Float_t driftlength, Float_t *xyz)
 151 {
 152   //
 153   // Applies the diffusion smearing to the position of a single electron
 154   //
 155
 156   if ((driftlength >        0) &&
 157       (driftlength < kDrThick)) {
 158     Float_t driftSqrt = TMath::Sqrt(driftlength);
 159     Float_t sigmaT = driftSqrt * fDiffusionT;
 160     Float_t sigmaL = driftSqrt * fDiffusionL;
 161     xyz[0] = gRandom->Gaus(xyz[0], sigmaL);
 162     xyz[1] = gRandom->Gaus(xyz[1], sigmaT);
 163     xyz[2] = gRandom->Gaus(xyz[2], sigmaT);
 164   }
 165   else {
 166     xyz[0] = 0.0;
 167     xyz[1] = 0.0;
 168     xyz[2] = 0.0;
 169   }
 170
 171 }
 172
 173 //_____________________________________________________________________________
 174 void AliTRDv2::Hits2Digits()
 175 {
 176   //
 177   // Creates TRD digits from hits. This procedure includes the following:
 178   //      - Diffusion
 179   //      - Gas gain including fluctuations
 180   //      - Pad-response (simple Gaussian approximation)
 181   //      - Electronics noise
 182   //      - Electronics gain
 183   //      - Digitization
 184   //      - ADC threshold
 185   // The corresponding parameter can be adjusted via the various Set-functions.
 186   // If these parameters are not explicitly set, default values are used (see
 187   // Init-function).
 188   // To produce digits from a galice.root file with TRD-hits use the
 189   // digitsCreate.C macro.
 190   //
 191
 192   printf(" Start creating digits\n");
 193
 194   ///////////////////////////////////////////////////////////////
 195   // Parameter
 196   ///////////////////////////////////////////////////////////////
 197
 198   // Converts number of electrons to fC
 199   const Float_t el2fC = 1.602E-19 * 1.0E15;
 200
 201   ///////////////////////////////////////////////////////////////
 202
 203   Int_t nBytes = 0;
 204
 205   AliTRDhit *TRDhit;
 206
 207   Int_t iplan;
 208   Int_t iRow;
 209
 210   // Position of pad 0,0,0
 211   //
 212   // chambers seen from the top:
 213   //     +----------------------------+
 214   //     |                            |
 215   //     |                            |     ^
 216   //     |                            | rphi|
 217   //     |                            |     |
 218   //     |0                           |     |
 219   //     +----------------------------+     +------>
 220   //                                             z
 221   // chambers seen from the side:           ^
 222   //     +----------------------------+ time|
 223   //     |                            |     |
 224   //     |0                           |     |
 225   //     +----------------------------+     +------>
 226   //                                             z
 227   //
 228   // The pad row (z-direction)
 229   Float_t row0[kNplan][kNcham];
 230   for (iplan = 0; iplan < kNplan; iplan++) {
 231     row0[iplan][0] = -fClengthI[iplan]/2. - fClengthM[iplan] - fClengthO[iplan]
 232                    + kCcthick;
 233     row0[iplan][1] = -fClengthI[iplan]/2. - fClengthM[iplan]
 234                    + kCcthick;
 235     row0[iplan][2] = -fClengthI[iplan]/2.
 236                    + kCcthick;
 237     row0[iplan][3] =  fClengthI[iplan]/2.
 238                    + kCcthick;
 239     row0[iplan][4] =  fClengthI[iplan]/2. + fClengthM[iplan]
 240                    + kCcthick;
 241   }
 242   // The pad column (rphi-direction)
 243   Float_t col0[kNplan];
 244   for (iplan = 0; iplan < kNplan; iplan++) {
 245     col0[iplan]    = -fCwidth[iplan]/2. + kCcthick;
 246   }
 247   // The time bucket
 248   Float_t time0[kNplan];
 249   for (iplan = 0; iplan < kNplan; iplan++) {
 250     time0[iplan]   = kRmin + kCcframe/2. + kDrZpos - 0.5 * kDrThick
 251                            + iplan * (kCheight + kCspace);
 252   }
 253
 254   // Get the pointer to the hit tree
 255   TTree *HitTree    = gAlice->TreeH();
 256   // Get the pointer to the digits tree
 257   TTree *DigitsTree = gAlice->TreeD();
 258
 259   // Get the number of entries in the hit tree
 260   // (Number of primary particles creating a hit somewhere)
 261   Int_t nTrack = (Int_t) HitTree->GetEntries();
 262
 263   Int_t chamBeg = 0;
 264   Int_t chamEnd = kNcham;
 265   if (fSensChamber) chamEnd = chamBeg = fSensChamber;
 266   Int_t planBeg = 0;
 267   Int_t planEnd = kNplan;
 268   if (fSensPlane)   planEnd = planBeg = fSensPlane;
 269   Int_t sectBeg = 0;
 270   Int_t sectEnd = kNsect;
 271   if (fSensSector)  sectEnd = sectBeg = fSensSector;
 272
 273   // Loop through all the chambers
 274   for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
 275     for (iplan = planBeg; iplan < planEnd; iplan++) {
 276       for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
 277
 278         printf(" Digitizing chamber %d, plane %d, sector %d\n"
 279               ,icham+1,iplan+1,isect+1);
 280
 281         // Create a detector matrix to keep the signal and track numbers
 282         AliTRDmatrix *matrix = new AliTRDmatrix(fRowMax[iplan][icham]
 283                                                ,fColMax[iplan]
 284                                                ,fTimeMax
 285                                                ,isect+1,icham+1,iplan+1);
 286
 287         // Loop through all entries in the tree
 288         for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
 289
 290           gAlice->ResetHits();
 291           nBytes += HitTree->GetEvent(iTrack);
 292
 293           // Get the number of hits in the TRD created by this particle
 294           Int_t nHit = fHits->GetEntriesFast();
 295
 296           // Loop through the TRD hits
 297           for (Int_t iHit = 0; iHit < nHit; iHit++) {
 298
 299             if (!(TRDhit = (AliTRDhit *) fHits->UncheckedAt(iHit)))
 300               continue;
 301
 302             Float_t x       = TRDhit->fX;
 303             Float_t y       = TRDhit->fY;
 304             Float_t z       = TRDhit->fZ;
 305             Float_t q       = TRDhit->fQ;
 306             Int_t   track   = TRDhit->fTrack;
 307             Int_t   plane   = TRDhit->fPlane;
 308             Int_t   sector  = TRDhit->fSector;
 309             Int_t   chamber = TRDhit->fChamber;
 310
 311             if ((sector  != isect+1) ||
 312                 (plane   != iplan+1) ||
 313                 (chamber != icham+1))
 314               continue;
 315
 316             // Rotate the sectors on top of each other
 317             Float_t phi  = 2.0 * kPI /  (Float_t) kNsect
 318                                * ((Float_t) sector - 0.5);
 319             Float_t xRot = -x * TMath::Cos(phi) + y * TMath::Sin(phi);
 320             Float_t yRot =  x * TMath::Sin(phi) + y * TMath::Cos(phi);
 321             Float_t zRot =  z;
 322
 323             // The hit position in pad coordinates (center pad)
 324             // The pad row (z-direction)
 325             Int_t rowH  = (Int_t) ((zRot -  row0[iplan][icham]) / fRowPadSize);
 326             // The pad column (rphi-direction)
 327             Int_t colH  = (Int_t) ((yRot -  col0[iplan]       ) / fColPadSize);
 328             // The time bucket
 329             Int_t timeH = (Int_t) ((xRot - time0[iplan]       ) / fTimeBinSize);
 330
 331             // Array to sum up the signal in a box surrounding the
 332             // hit postition
 333             const Int_t timeBox = 5;
 334             const Int_t  colBox = 7;
 335             const Int_t  rowBox = 5;
 336             Float_t signalSum[rowBox][colBox][timeBox];
 337             for (iRow  = 0;  iRow <  rowBox; iRow++ ) {
 338               for (Int_t iCol  = 0;  iCol <  colBox; iCol++ ) {
 339                 for (Int_t iTime = 0; iTime < timeBox; iTime++) {
 340                   signalSum[iRow][iCol][iTime] = 0;
 341                 }
 342               }
 343             }
 344
 345             // Loop over all electrons of this hit
 346             Int_t nEl = (Int_t) q;
 347             for (Int_t iEl = 0; iEl < nEl; iEl++) {
 348
 349               // Apply the diffusion smearing
 350               Float_t driftlength = xRot - time0[iplan];
 351               Float_t xyz[3];
 352               xyz[0] = xRot;
 353               xyz[1] = yRot;
 354               xyz[2] = zRot;
 355               Diffusion(driftlength,xyz);
 356
 357               // At this point absorption effects that depend on the
 358               // driftlength could be taken into account.
 359
 360               // The electron position and the distance to the hit position
 361               // in pad units
 362               // The pad row (z-direction)
 363               Int_t  rowE = (Int_t) ((xyz[2] -  row0[iplan][icham]) / fRowPadSize);
 364               Int_t  rowD =  rowH -  rowE;
 365               // The pad column (rphi-direction)
 366               Int_t  colE = (Int_t) ((xyz[1] -  col0[iplan]       ) / fColPadSize);
 367               Int_t  colD =  colH -  colE;
 368               // The time bucket
 369               Int_t timeE = (Int_t) ((xyz[0] - time0[iplan]       ) / fTimeBinSize);
 370               Int_t timeD = timeH - timeE;
 371
 372               // Apply the gas gain including fluctuations
 373               Int_t signal = (Int_t) (-fGasGain * TMath::Log(gRandom->Rndm()));
 374
 375               // The distance of the electron to the center of the pad
 376               // in units of pad width
 377               Float_t dist = (xyz[1] - col0[iplan] - (colE + 0.5) * fColPadSize)
 378                            / fColPadSize;
 379
 380               // Sum up the signal in the different pixels
 381               // and apply the pad response
 382               Int_t  rowIdx =  rowD + (Int_t) ( rowBox / 2);
 383               Int_t  colIdx =  colD + (Int_t) ( colBox / 2);
 384               Int_t timeIdx = timeD + (Int_t) (timeBox / 2);
 385               signalSum[rowIdx][colIdx-1][timeIdx] += PadResponse(dist-1.) * signal;
 386               signalSum[rowIdx][colIdx  ][timeIdx] += PadResponse(dist   ) * signal;
 387               signalSum[rowIdx][colIdx+1][timeIdx] += PadResponse(dist+1.) * signal;
 388
 389             }
 390
 391             // Add the padcluster to the detector matrix
 392             for (iRow  = 0;  iRow <  rowBox; iRow++ ) {
 393               for (Int_t iCol  = 0;  iCol <  colBox; iCol++ ) {
 394                 for (Int_t iTime = 0; iTime < timeBox; iTime++) {
 395
 396                   Int_t  rowB =  rowH + iRow  - (Int_t) ( rowBox / 2);
 397                   Int_t  colB =  colH + iCol  - (Int_t) ( colBox / 2);
 398                   Int_t timeB = timeH + iTime - (Int_t) (timeBox / 2);
 399
 400                   Float_t signalB = signalSum[iRow][iCol][iTime];
 401                   if (signalB > 0.0) {
 402                     matrix->AddSignal(rowB,colB,timeB,signalB);
 403                     if (!(matrix->AddTrack(rowB,colB,timeB,track)))
 404                       printf("More than three tracks in a pixel!\n");
 405                   }
 406
 407                 }
 408               }
 409             }
 410
 411           }
 412
 413         }
 414
 415         // Create the hits for this chamber
 416         for (Int_t iRow  = 0; iRow  <  fRowMax[iplan][icham]; iRow++ ) {
 417           for (Int_t iCol  = 0; iCol  <  fColMax[iplan]         ; iCol++ ) {
 418             for (Int_t iTime = 0; iTime < fTimeMax                ; iTime++) {
 419
 420               Float_t signalAmp = matrix->GetSignal(iRow,iCol,iTime);
 421
 422               // Add the noise
 423               signalAmp  = TMath::Max(gRandom->Gaus(signalAmp,fNoise),(Float_t) 0.0);
 424               // Convert to fC
 425               signalAmp *= el2fC;
 426               // Convert to mV
 427               signalAmp *= fChipGain;
 428               // Convert to ADC counts
 429               Int_t adc  = (Int_t) (signalAmp * (fADCoutRange / fADCinRange));
 430
 431               // Apply threshold on ADC value
 432               if (adc > fADCthreshold) {
 433
 434                 Int_t trackSave[3];
 435                 for (Int_t ii = 0; ii < 3; ii++) {
 436                   trackSave[ii] = matrix->GetTrack(iRow,iCol,iTime,ii);
 437                 }
 438
 439                 Int_t digits[7];
 440                 digits[0] = matrix->GetSector();
 441                 digits[1] = matrix->GetChamber();
 442                 digits[2] = matrix->GetPlane();
 443                 digits[3] = iRow;
 444                 digits[4] = iCol;
 445                 digits[5] = iTime;
 446                 digits[6] = adc;
 447
 448                 // Add this digit to the TClonesArray
 449                 AddDigit(trackSave,digits);
 450
 451               }
 452
 453             }
 454           }
 455         }
 456
 457         // Clean up
 458         delete matrix;
 459
 460       }
 461     }
 462   }
 463
 464   // Fill the digits-tree
 465   DigitsTree->Fill();
 466
 467 }
 468
 469 //_____________________________________________________________________________
 470 void AliTRDv2::Init()
 471 {
 472   //
 473   // Initialise Transition Radiation Detector after geometry has been built.
 474   // Includes the default settings of all parameter for the digitization.
 475   //
 476
 477   AliTRD::Init();
 478
 479   if (fSensPlane)
 480     printf("          Only plane %d is sensitive\n",fSensPlane);
 481   if (fSensChamber)
 482     printf("          Only chamber %d is sensitive\n",fSensChamber);
 483   if (fSensSector)
 484     printf("          Only sector %d is sensitive\n",fSensSector);
 485
 486   for (Int_t i = 0; i < 80; i++) printf("*");
 487   printf("\n");
 488
 489   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 490   const Float_t kPoti = 12.1;
 491   // Maximum energy (50 keV);
 492   const Float_t kEend = 50000.0;
 493   // Ermilova distribution for the delta-ray spectrum
 494   Float_t Poti = TMath::Log(kPoti);
 495   Float_t Eend = TMath::Log(kEend);
 496   fDeltaE  = new TF1("deltae",Ermilova,Poti,Eend,0);
 497
 498   // Identifier of the sensitive volume (drift region)
 499   fIdSens     = gMC->VolId("UL05");
 500
 501   // Identifier of the TRD-spaceframe volumina
 502   fIdSpace1   = gMC->VolId("B028");
 503   fIdSpace2   = gMC->VolId("B029");
 504   fIdSpace3   = gMC->VolId("B030");
 505
 506   // Identifier of the TRD-driftchambers
 507   fIdChamber1 = gMC->VolId("UCIO");
 508   fIdChamber2 = gMC->VolId("UCIM");
 509   fIdChamber3 = gMC->VolId("UCII");
 510
 511   // The default pad dimensions
 512   if (!(fRowPadSize))  fRowPadSize  = 4.5;
 513   if (!(fColPadSize))  fColPadSize  = 1.0;
 514   if (!(fTimeBinSize)) fTimeBinSize = 0.1;
 515
 516   // The maximum number of pads
 517   for (Int_t iplan = 0; iplan < kNplan; iplan++) {
 518     // Rows
 519     fRowMax[iplan][0] = 1 + TMath::Nint((fClengthO[iplan] - 2. * kCcthick)
 520                                                           / fRowPadSize - 0.5);
 521     fRowMax[iplan][1] = 1 + TMath::Nint((fClengthM[iplan] - 2. * kCcthick)
 522                                                           / fRowPadSize - 0.5);
 523     fRowMax[iplan][2] = 1 + TMath::Nint((fClengthI[iplan] - 2. * kCcthick)
 524                                                           / fRowPadSize - 0.5);
 525     fRowMax[iplan][3] = 1 + TMath::Nint((fClengthM[iplan] - 2. * kCcthick)
 526                                                           / fRowPadSize - 0.5);
 527     fRowMax[iplan][4] = 1 + TMath::Nint((fClengthO[iplan] - 2. * kCcthick)
 528                                                           / fRowPadSize - 0.5);
 529     // Columns
 530     fColMax[iplan]    = 1 + TMath::Nint((fCwidth[iplan]   - 2. * kCcthick)
 531                                                           / fColPadSize - 0.5);
 532   }
 533   // Time buckets
 534   fTimeMax = 1 + TMath::Nint(kDrThick / fTimeBinSize - 0.5);
 535
 536   // The default parameter for the digitization
 537   if (!(fGasGain))      fGasGain      = 2.0E3;
 538   if (!(fNoise))        fNoise        = 3000.;
 539   if (!(fChipGain))     fChipGain     = 10.;
 540   if (!(fADCoutRange))  fADCoutRange  = 255.;
 541   if (!(fADCinRange))   fADCinRange   = 2000.;
 542   if (!(fADCthreshold)) fADCthreshold = 0;
 543
 544   // Transverse and longitudinal diffusion coefficients (Xe/Isobutane)
 545   if (!(fDiffusionT))   fDiffusionT   = 0.060;
 546   if (!(fDiffusionL))   fDiffusionL   = 0.017;
 547
 548 }
 549
 550 //_____________________________________________________________________________
 551 void AliTRDv2::MakeBranch(Option_t* option)
 552 {
 553   //
 554   // Create Tree branches for the TRD digits.
 555   //
 556
 557   Int_t  buffersize = 4000;
 558   Char_t branchname[10];
 559
 560   sprintf(branchname,"%s",GetName());
 561
 562   AliDetector::MakeBranch(option);
 563
 564   Char_t *D = strstr(option,"D");
 565   if (fDigits && gAlice->TreeD() && D) {
 566     gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
 567     printf("Making Branch %s for digits\n",branchname);
 568   }
 569
 570 }
 571
 572 //_____________________________________________________________________________
 573 Float_t AliTRDv2::PadResponse(Float_t x)
 574 {
 575   //
 576   // The pad response for the chevron pads.
 577   // We use a simple Gaussian approximation which should be good
 578   // enough for our purpose.
 579   //
 580
 581   // The parameters for the response function
 582   const Float_t aa  =  0.8872;
 583   const Float_t bb  = -0.00573;
 584   const Float_t cc  =  0.454;
 585   const Float_t cc2 =  cc*cc;
 586
 587   Float_t pr = aa * (bb + TMath::Exp(-x*x / (2. * cc2)));
 588
 589   //TF1 *funPR = new TF1("funPR","[0]*([1]+exp(-x*x /(2.*[2])))",-3,3);
 590   //funPR->SetParameter(0,aa );
 591   //funPR->SetParameter(1,bb );
 592   //funPR->SetParameter(2,cc2);
 593   //
 594   //Float_t pr = funPR->Eval(distance);
 595   //
 596   //delete funPR;
 597
 598   return (pr);
 599
 600 }
 601
 602 //_____________________________________________________________________________
 603 void AliTRDv2::SetSensPlane(Int_t iplane)
 604 {
 605   //
 606   // Defines the hit-sensitive plane (1-6)
 607   //
 608
 609   if ((iplane < 0) || (iplane > 6)) {
 610     printf("Wrong input value: %d\n",iplane);
 611     printf("Use standard setting\n");
 612     fSensPlane  = 0;
 613     fSensSelect = 0;
 614     return;
 615   }
 616
 617   fSensSelect = 1;
 618   fSensPlane  = iplane;
 619
 620 }
 621
 622 //_____________________________________________________________________________
 623 void AliTRDv2::SetSensChamber(Int_t ichamber)
 624 {
 625   //
 626   // Defines the hit-sensitive chamber (1-5)
 627   //
 628
 629   if ((ichamber < 0) || (ichamber > 5)) {
 630     printf("Wrong input value: %d\n",ichamber);
 631     printf("Use standard setting\n");
 632     fSensChamber = 0;
 633     fSensSelect  = 0;
 634     return;
 635   }
 636
 637   fSensSelect  = 1;
 638   fSensChamber = ichamber;
 639
 640 }
 641
 642 //_____________________________________________________________________________
 643 void AliTRDv2::SetSensSector(Int_t isector)
 644 {
 645   //
 646   // Defines the hit-sensitive sector (1-18)
 647   //
 648
 649   if ((isector < 0) || (isector > 18)) {
 650     printf("Wrong input value: %d\n",isector);
 651     printf("Use standard setting\n");
 652     fSensSector = 0;
 653     fSensSelect = 0;
 654     return;
 655   }
 656
 657   fSensSelect = 1;
 658   fSensSector = isector;
 659
 660 }
 661
 662 //_____________________________________________________________________________
 663 void AliTRDv2::StepManager()
 664 {
 665   //
 666   // Called at every step in the Transition Radiation Detector version 2.
 667   // Slow simulator. Every charged track produces electron cluster as hits
 668   // along its path across the drift volume. The step size is set acording
 669   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 670   // a spectrum taken from Ermilova et al.
 671   //
 672
 673   Int_t    iIdSens, icSens;
 674   Int_t    iIdSpace, icSpace;
 675   Int_t    iIdChamber, icChamber;
 676   Int_t    vol[3];
 677   Int_t    iPid;
 678
 679   Int_t    secMap1[10] = {  3,  7,  8,  9, 10, 11,  2,  1, 18, 17 };
 680   Int_t    secMap2[ 5] = { 16, 15, 14, 13, 12 };
 681   Int_t    secMap3[ 3] = {  5,  6,  4 };
 682
 683   Float_t  hits[4];
 684   Float_t  random[1];
 685   Float_t  charge;
 686   Float_t  aMass;
 687
 688   Double_t pTot;
 689   Double_t qTot;
 690   Double_t eDelta;
 691   Double_t betaGamma, pp;
 692
 693   TLorentzVector pos, mom;
 694   TClonesArray  &lhits = *fHits;
 695
 696   const Double_t kBig = 1.0E+12;
 697
 698   // Ionization energy
 699   const Float_t kWion    = 22.04;
 700   // Maximum energy for e+ e- g for the step-size calculation
 701   const Float_t kPTotMax = 0.002;
 702   // Plateau value of the energy-loss for electron in xenon
 703   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 704   //const Double_t kPlateau = 1.70;
 705   // the averaged value (26/3/99)
 706   const Float_t kPlateau = 1.55;
 707   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 708   const Float_t kPrim    = 48.0;
 709   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 710   const Float_t kPoti    = 12.1;
 711
 712   // Set the maximum step size to a very large number for all
 713   // neutral particles and those outside the driftvolume
 714   gMC->SetMaxStep(kBig);
 715
 716   // Use only charged tracks
 717   if (( gMC->TrackCharge()       ) &&
 718       (!gMC->IsTrackStop()       ) &&
 719       (!gMC->IsTrackDisappeared())) {
 720
 721     // Inside a sensitive volume?
 722     iIdSens = gMC->CurrentVolID(icSens);
 723     if (iIdSens == fIdSens) {
 724
 725       iIdSpace   = gMC->CurrentVolOffID(4,icSpace  );
 726       iIdChamber = gMC->CurrentVolOffID(1,icChamber);
 727
 728       // Calculate the energy of the delta-electrons
 729       eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 730       eDelta = TMath::Max(eDelta,0.0);
 731
 732       // The number of secondary electrons created
 733       qTot = (Double_t) ((Int_t) (eDelta / kWion) + 1);
 734
 735       // The hit coordinates and charge
 736       gMC->TrackPosition(pos);
 737       hits[0] = pos[0];
 738       hits[1] = pos[1];
 739       hits[2] = pos[2];
 740       hits[3] = qTot;
 741
 742       // The sector number
 743       if      (iIdSpace == fIdSpace1)
 744         vol[0] = secMap1[icSpace-1];
 745       else if (iIdSpace == fIdSpace2)
 746         vol[0] = secMap2[icSpace-1];
 747       else if (iIdSpace == fIdSpace3)
 748         vol[0] = secMap3[icSpace-1];
 749
 750       // The chamber number
 751       //   1: outer left
 752       //   2: middle left
 753       //   3: inner
 754       //   4: middle right
 755       //   5: outer right
 756       if      (iIdChamber == fIdChamber1)
 757         vol[1] = (hits[2] < 0 ? 1 : 5);
 758       else if (iIdChamber == fIdChamber2)
 759         vol[1] = (hits[2] < 0 ? 2 : 4);
 760       else if (iIdChamber == fIdChamber3)
 761         vol[1] = 3;
 762
 763       // The plane number
 764       vol[2] = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6;
 765
 766       // Check on selected volumes
 767       Int_t addthishit = 1;
 768       if (fSensSelect) {
 769         if ((fSensPlane)   && (vol[2] != fSensPlane  )) addthishit = 0;
 770         if ((fSensChamber) && (vol[1] != fSensChamber)) addthishit = 0;
 771         if ((fSensSector)  && (vol[0] != fSensSector )) addthishit = 0;
 772       }
 773
 774       // Add this hit
 775       if (addthishit) {
 776
 777         new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
 778
 779         // The energy loss according to Bethe Bloch
 780         gMC->TrackMomentum(mom);
 781         pTot = mom.Rho();
 782         iPid = gMC->TrackPid();
 783         if ( (iPid >  3) ||
 784             ((iPid <= 3) && (pTot < kPTotMax))) {
 785           aMass     = gMC->TrackMass();
 786           betaGamma = pTot / aMass;
 787           pp        = kPrim * BetheBloch(betaGamma);
 788           // Take charge > 1 into account
 789           charge = gMC->TrackCharge();
 790           if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 791         }
 792         // Electrons above 20 Mev/c are at the plateau
 793         else {
 794           pp = kPrim * kPlateau;
 795         }
 796
 797         // Calculate the maximum step size for the next tracking step
 798         if (pp > 0) {
 799           do
 800             gMC->Rndm(random,1);
 801           while ((random[0] == 1.) || (random[0] == 0.));
 802           gMC->SetMaxStep( - TMath::Log(random[0]) / pp);
 803         }
 804
 805       }
 806       else {
 807         // set step size to maximal value
 808         gMC->SetMaxStep(kBig);
 809       }
 810
 811     }
 812
 813   }
 814
 815 }
 816
 817 //_____________________________________________________________________________
 818 Double_t AliTRDv2::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   if (bg > 0) {
 841     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
 842     Double_t aa = TMath::Power(yy,kP4);
 843     Double_t bb = TMath::Power((1./bg),kP5);
 844              bb = TMath::Log(kP3 + bb);
 845     return ((kP2 - aa - bb)*kP1 / aa);
 846   }
 847   else
 848     return 0;
 849
 850 }
 851
 852 //_____________________________________________________________________________
 853 Double_t Ermilova(Double_t *x, Double_t *)
 854 {
 855   //
 856   // Calculates the delta-ray energy distribution according to Ermilova.
 857   // Logarithmic scale !
 858   //
 859
 860   Double_t energy;
 861   Double_t dpos;
 862   Double_t dnde;
 863
 864   Int_t    pos1, pos2;
 865
 866   const Int_t nV = 31;
 867
 868   Float_t vxe[nV] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
 869                     , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
 870                     , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
 871                     , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
 872                     , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
 873                     , 9.4727, 9.9035,10.3735,10.5966,10.8198
 874                     ,11.5129 };
 875
 876   Float_t vye[nV] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
 877                     , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
 878                     ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
 879                     ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
 880                     ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
 881                     ,  0.04 ,  0.023,  0.015,  0.011,  0.01
 882                     ,  0.004 };
 883
 884   energy = x[0];
 885
 886   // Find the position
 887   pos1 = pos2 = 0;
 888   dpos = 0;
 889   do {
 890     dpos = energy - vxe[pos2++];
 891   }
 892   while (dpos > 0);
 893   pos2--;
 894   if (pos2 > nV) pos2 = nV;
 895   pos1 = pos2 - 1;
 896
 897   // Differentiate between the sampling points
 898   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
 899
 900   return dnde;
 901
 902 }