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