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 /*
  17 $Log$
  18 Revision 1.17.2.1  2000/05/08 14:59:16  cblume
  19 Made inline function non-virtual. Bug fix in setting sensitive chamber
  20
  21 Revision 1.17  2000/02/28 19:10:26  cblume
  22 Include the new TRD classes
  23
  24 Revision 1.16.4.1  2000/02/28 18:04:35  cblume
  25 Change to new hit version, introduce geometry class, and move digitization and clustering to AliTRDdigitizer/AliTRDclusterizerV1
  26
  27 Revision 1.16  1999/11/05 22:50:28  fca
  28 Do not use Atan, removed from ROOT too
  29
  30 Revision 1.15  1999/11/02 17:20:19  fca
  31 initialise nbytes before using it
  32
  33 Revision 1.14  1999/11/02 17:15:54  fca
  34 Correct ansi scoping not accepted by HP compilers
  35
  36 Revision 1.13  1999/11/02 17:14:51  fca
  37 Correct ansi scoping not accepted by HP compilers
  38
  39 Revision 1.12  1999/11/02 16:35:56  fca
  40 New version of TRD introduced
  41
  42 Revision 1.11  1999/11/01 20:41:51  fca
  43 Added protections against using the wrong version of FRAME
  44
  45 Revision 1.10  1999/09/29 09:24:35  fca
  46 Introduction of the Copyright and cvs Log
  47
  48 */
  49
  50 ///////////////////////////////////////////////////////////////////////////////
  51 //                                                                           //
  52 //  Transition Radiation Detector version 2 -- slow simulator                //
  53 //                                                                           //
  54 //Begin_Html
  55 /*
  56 <img src="picts/AliTRDfullClass.gif">
  57 */
  58 //End_Html
  59 //                                                                           //
  60 //                                                                           //
  61 ///////////////////////////////////////////////////////////////////////////////
  62
  63 #include <TMath.h>
  64 #include <TVector.h>
  65 #include <TRandom.h>
  66
  67 #include "AliRun.h"
  68 #include "AliMC.h"
  69 #include "AliConst.h"
  70
  71 #include "AliTRDv1.h"
  72 #include "AliTRDmatrix.h"
  73 #include "AliTRDgeometry.h"
  74
  75 ClassImp(AliTRDv1)
  76
  77 //_____________________________________________________________________________
  78 AliTRDv1::AliTRDv1(const char *name, const char *title)
  79          :AliTRD(name, title)
  80 {
  81   //
  82   // Standard constructor for Transition Radiation Detector version 1
  83   //
  84
  85   fIdSens        =  0;
  86
  87   fIdChamber1    =  0;
  88   fIdChamber2    =  0;
  89   fIdChamber3    =  0;
  90
  91   fSensSelect    =  0;
  92   fSensPlane     = -1;
  93   fSensChamber   = -1;
  94   fSensSector    = -1;
  95
  96   fDeltaE        = NULL;
  97
  98   SetBufferSize(128000);
  99
 100 }
 101
 102 //_____________________________________________________________________________
 103 AliTRDv1::~AliTRDv1()
 104 {
 105
 106   if (fDeltaE) delete fDeltaE;
 107
 108 }
 109
 110 //_____________________________________________________________________________
 111 void AliTRDv1::CreateGeometry()
 112 {
 113   //
 114   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
 115   // This version covers the full azimuth.
 116   //
 117
 118   // Check that FRAME is there otherwise we have no place where to put the TRD
 119   AliModule* FRAME = gAlice->GetModule("FRAME");
 120   if (!FRAME) return;
 121
 122   // Define the chambers
 123   AliTRD::CreateGeometry();
 124
 125 }
 126
 127 //_____________________________________________________________________________
 128 void AliTRDv1::CreateMaterials()
 129 {
 130   //
 131   // Create materials for the Transition Radiation Detector version 1
 132   //
 133
 134   AliTRD::CreateMaterials();
 135
 136 }
 137
 138 //_____________________________________________________________________________
 139 void AliTRDv1::Init()
 140 {
 141   //
 142   // Initialise Transition Radiation Detector after geometry has been built.
 143   //
 144
 145   AliTRD::Init();
 146
 147   printf("          Slow simulator\n\n");
 148   if (fSensSelect) {
 149     if (fSensPlane   >= 0)
 150       printf("          Only plane %d is sensitive\n",fSensPlane);
 151     if (fSensChamber >= 0)
 152       printf("          Only chamber %d is sensitive\n",fSensChamber);
 153     if (fSensSector  >= 0)
 154       printf("          Only sector %d is sensitive\n",fSensSector);
 155   }
 156   printf("\n");
 157
 158   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 159   const Float_t kPoti = 12.1;
 160   // Maximum energy (50 keV);
 161   const Float_t kEend = 50000.0;
 162   // Ermilova distribution for the delta-ray spectrum
 163   Float_t Poti = TMath::Log(kPoti);
 164   Float_t Eend = TMath::Log(kEend);
 165   fDeltaE  = new TF1("deltae",Ermilova,Poti,Eend,0);
 166
 167   // Identifier of the sensitive volume (drift region)
 168   fIdSens     = gMC->VolId("UL05");
 169
 170   // Identifier of the TRD-driftchambers
 171   fIdChamber1 = gMC->VolId("UCIO");
 172   fIdChamber2 = gMC->VolId("UCIM");
 173   fIdChamber3 = gMC->VolId("UCII");
 174
 175   for (Int_t i = 0; i < 80; i++) printf("*");
 176   printf("\n");
 177
 178 }
 179
 180 //_____________________________________________________________________________
 181 void AliTRDv1::SetSensPlane(Int_t iplane)
 182 {
 183   //
 184   // Defines the hit-sensitive plane (0-5)
 185   //
 186
 187   if ((iplane < 0) || (iplane > 5)) {
 188     printf("Wrong input value: %d\n",iplane);
 189     printf("Use standard setting\n");
 190     fSensPlane  = -1;
 191     fSensSelect =  0;
 192     return;
 193   }
 194
 195   fSensSelect = 1;
 196   fSensPlane  = iplane;
 197
 198 }
 199
 200 //_____________________________________________________________________________
 201 void AliTRDv1::SetSensChamber(Int_t ichamber)
 202 {
 203   //
 204   // Defines the hit-sensitive chamber (0-4)
 205   //
 206
 207   if ((ichamber < 0) || (ichamber > 4)) {
 208     printf("Wrong input value: %d\n",ichamber);
 209     printf("Use standard setting\n");
 210     fSensChamber = -1;
 211     fSensSelect  =  0;
 212     return;
 213   }
 214
 215   fSensSelect  = 1;
 216   fSensChamber = ichamber;
 217
 218 }
 219
 220 //_____________________________________________________________________________
 221 void AliTRDv1::SetSensSector(Int_t isector)
 222 {
 223   //
 224   // Defines the hit-sensitive sector (0-17)
 225   //
 226
 227   if ((isector < 0) || (isector > 17)) {
 228     printf("Wrong input value: %d\n",isector);
 229     printf("Use standard setting\n");
 230     fSensSector = -1;
 231     fSensSelect =  0;
 232     return;
 233   }
 234
 235   fSensSelect = 1;
 236   fSensSector = isector;
 237
 238 }
 239
 240 //_____________________________________________________________________________
 241 void AliTRDv1::StepManager()
 242 {
 243   //
 244   // Slow simulator. Every charged track produces electron cluster as hits
 245   // along its path across the drift volume. The step size is set acording
 246   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 247   // a spectrum taken from Ermilova et al.
 248   //
 249
 250   Int_t    iIdSens, icSens;
 251   Int_t    iIdSpace, icSpace;
 252   Int_t    iIdChamber, icChamber;
 253   Int_t    pla = 0;
 254   Int_t    cha = 0;
 255   Int_t    sec = 0;
 256   Int_t    iPdg;
 257
 258   Float_t  hits[4];
 259   Float_t  random[1];
 260   Float_t  charge;
 261   Float_t  aMass;
 262
 263   Double_t pTot;
 264   Double_t qTot;
 265   Double_t eDelta;
 266   Double_t betaGamma, pp;
 267
 268   TLorentzVector pos, mom;
 269   TClonesArray  &lhits = *fHits;
 270
 271   const Double_t kBig     = 1.0E+12;
 272
 273   // Ionization energy
 274   const Float_t  kWion    = 22.04;
 275   // Maximum energy for e+ e- g for the step-size calculation
 276   const Float_t  kPTotMax = 0.002;
 277   // Plateau value of the energy-loss for electron in xenon
 278   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 279   //const Double_t kPlateau = 1.70;
 280   // the averaged value (26/3/99)
 281   const Float_t  kPlateau = 1.55;
 282   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 283   const Float_t  kPrim    = 48.0;
 284   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 285   const Float_t  kPoti    = 12.1;
 286
 287   // PDG code electron
 288   const Int_t    pdgElectron = 11;
 289
 290   // Set the maximum step size to a very large number for all
 291   // neutral particles and those outside the driftvolume
 292   gMC->SetMaxStep(kBig);
 293
 294   // Use only charged tracks
 295   if (( gMC->TrackCharge()       ) &&
 296       (!gMC->IsTrackStop()       ) &&
 297       (!gMC->IsTrackDisappeared())) {
 298
 299     // Inside a sensitive volume?
 300     iIdSens = gMC->CurrentVolID(icSens);
 301     if (iIdSens == fIdSens) {
 302
 303       iIdSpace   = gMC->CurrentVolOffID(4,icSpace  );
 304       iIdChamber = gMC->CurrentVolOffID(1,icChamber);
 305
 306       // Calculate the energy of the delta-electrons
 307       eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 308       eDelta = TMath::Max(eDelta,0.0);
 309
 310       // The number of secondary electrons created
 311       qTot = (Double_t) ((Int_t) (eDelta / kWion) + 1);
 312
 313       // The hit coordinates and charge
 314       gMC->TrackPosition(pos);
 315       hits[0] = pos[0];
 316       hits[1] = pos[1];
 317       hits[2] = pos[2];
 318       hits[3] = qTot;
 319
 320       // The sector number (0 - 17)
 321       // The numbering goes clockwise and starts at y = 0
 322       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 323       if (phi < 90.)
 324         phi = phi + 270.;
 325       else
 326         phi = phi -  90.;
 327       sec = ((Int_t) (phi / 20));
 328
 329       // The chamber number
 330       //   0: outer left
 331       //   1: middle left
 332       //   2: inner
 333       //   3: middle right
 334       //   4: outer right
 335       if      (iIdChamber == fIdChamber1)
 336         cha = (hits[2] < 0 ? 0 : 4);
 337       else if (iIdChamber == fIdChamber2)
 338         cha = (hits[2] < 0 ? 1 : 3);
 339       else if (iIdChamber == fIdChamber3)
 340         cha = 2;
 341
 342       // The plane number
 343       // The numbering starts at the innermost plane
 344       pla = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6 - 1;
 345
 346       // Check on selected volumes
 347       Int_t addthishit = 1;
 348       if (fSensSelect) {
 349         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 350         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 351         if ((fSensSector  >= 0) && (sec != fSensSector )) addthishit = 0;
 352       }
 353
 354       // Add this hit
 355       if (addthishit) {
 356
 357         new(lhits[fNhits++]) AliTRDhit(fIshunt
 358                                       ,gAlice->CurrentTrack()
 359                                       ,fGeometry->GetDetector(pla,cha,sec)
 360                                       ,hits);
 361
 362         // The energy loss according to Bethe Bloch
 363         gMC->TrackMomentum(mom);
 364         pTot = mom.Rho();
 365         iPdg = TMath::Abs(gMC->TrackPid());
 366         if ( (iPdg != pdgElectron) ||
 367             ((iPdg == pdgElectron) && (pTot < kPTotMax))) {
 368           aMass     = gMC->TrackMass();
 369           betaGamma = pTot / aMass;
 370           pp        = kPrim * BetheBloch(betaGamma);
 371           // Take charge > 1 into account
 372           charge = gMC->TrackCharge();
 373           if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 374         }
 375         // Electrons above 20 Mev/c are at the plateau
 376         else {
 377           pp = kPrim * kPlateau;
 378         }
 379
 380         // Calculate the maximum step size for the next tracking step
 381         if (pp > 0) {
 382           do
 383             gMC->Rndm(random,1);
 384           while ((random[0] == 1.) || (random[0] == 0.));
 385           gMC->SetMaxStep( - TMath::Log(random[0]) / pp);
 386         }
 387
 388       }
 389       else {
 390         // set step size to maximal value
 391         gMC->SetMaxStep(kBig);
 392       }
 393
 394     }
 395
 396   }
 397
 398 }
 399
 400 //_____________________________________________________________________________
 401 Double_t AliTRDv1::BetheBloch(Double_t bg)
 402 {
 403   //
 404   // Parametrization of the Bethe-Bloch-curve
 405   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
 406   //
 407
 408   // This parameters have been adjusted to averaged values from GEANT
 409   const Double_t kP1 = 7.17960e-02;
 410   const Double_t kP2 = 8.54196;
 411   const Double_t kP3 = 1.38065e-06;
 412   const Double_t kP4 = 5.30972;
 413   const Double_t kP5 = 2.83798;
 414
 415   // This parameters have been adjusted to Xe-data found in:
 416   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 417   //const Double_t kP1 = 0.76176E-1;
 418   //const Double_t kP2 = 10.632;
 419   //const Double_t kP3 = 3.17983E-6;
 420   //const Double_t kP4 = 1.8631;
 421   //const Double_t kP5 = 1.9479;
 422
 423   if (bg > 0) {
 424     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
 425     Double_t aa = TMath::Power(yy,kP4);
 426     Double_t bb = TMath::Power((1./bg),kP5);
 427              bb = TMath::Log(kP3 + bb);
 428     return ((kP2 - aa - bb)*kP1 / aa);
 429   }
 430   else
 431     return 0;
 432
 433 }
 434
 435 //_____________________________________________________________________________
 436 Double_t Ermilova(Double_t *x, Double_t *)
 437 {
 438   //
 439   // Calculates the delta-ray energy distribution according to Ermilova.
 440   // Logarithmic scale !
 441   //
 442
 443   Double_t energy;
 444   Double_t dpos;
 445   Double_t dnde;
 446
 447   Int_t    pos1, pos2;
 448
 449   const Int_t nV = 31;
 450
 451   Float_t vxe[nV] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
 452                     , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
 453                     , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
 454                     , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
 455                     , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
 456                     , 9.4727, 9.9035,10.3735,10.5966,10.8198
 457                     ,11.5129 };
 458
 459   Float_t vye[nV] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
 460                     , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
 461                     ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
 462                     ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
 463                     ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
 464                     ,  0.04 ,  0.023,  0.015,  0.011,  0.01
 465                     ,  0.004 };
 466
 467   energy = x[0];
 468
 469   // Find the position
 470   pos1 = pos2 = 0;
 471   dpos = 0;
 472   do {
 473     dpos = energy - vxe[pos2++];
 474   }
 475   while (dpos > 0);
 476   pos2--;
 477   if (pos2 > nV) pos2 = nV;
 478   pos1 = pos2 - 1;
 479
 480   // Differentiate between the sampling points
 481   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
 482
 483   return dnde;
 484
 485 }