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 ////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TLorentzVector.h>
  25 #include <TMath.h>
  26 #include <TRandom.h>
  27 #include <TVirtualMC.h>
  28 #include <TGeoManager.h>
  29 #include <TGeoMatrix.h>
  30 #include <TGeoPhysicalNode.h>
  31
  32 #include "AliTrackReference.h"
  33 #include "AliMC.h"
  34 #include "AliRun.h"
  35 #include "AliGeomManager.h"
  36
  37 #include "AliTRDgeometry.h"
  38 #include "AliTRDCommonParam.h"
  39 #include "AliTRDsimTR.h"
  40 #include "AliTRDv1.h"
  41
  42 ClassImp(AliTRDv1)
  43
  44 //_____________________________________________________________________________
  45 AliTRDv1::AliTRDv1()
  46   :AliTRD()
  47   ,fTRon(kTRUE)
  48   ,fTR(NULL)
  49   ,fStepSize(0)
  50   ,fWion(0)
  51 {
  52   //
  53   // Default constructor
  54   //
  55
  56 }
  57
  58 //_____________________________________________________________________________
  59 AliTRDv1::AliTRDv1(const char *name, const char *title)
  60   :AliTRD(name,title)
  61   ,fTRon(kTRUE)
  62   ,fTR(NULL)
  63   ,fStepSize(0.1)
  64   ,fWion(0)
  65 {
  66   //
  67   // Standard constructor for Transition Radiation Detector version 1
  68   //
  69
  70   SetBufferSize(128000);
  71
  72   if      (AliTRDCommonParam::Instance()->IsXenon()) {
  73     fWion = 23.53; // Ionization energy XeCO2 (85/15)
  74   }
  75   else if (AliTRDCommonParam::Instance()->IsArgon()) {
  76     fWion = 27.21; // Ionization energy ArCO2 (82/18)
  77   }
  78   else {
  79     AliFatal("Wrong gas mixture");
  80     exit(1);
  81   }
  82
  83 }
  84
  85 //_____________________________________________________________________________
  86 AliTRDv1::~AliTRDv1()
  87 {
  88   //
  89   // AliTRDv1 destructor
  90   //
  91
  92   if (fTR) {
  93     delete fTR;
  94     fTR     = 0;
  95   }
  96
  97 }
  98
  99 //_____________________________________________________________________________
 100 void AliTRDv1::AddAlignableVolumes() const
 101 {
 102   //
 103   // Create entries for alignable volumes associating the symbolic volume
 104   // name with the corresponding volume path. Needs to be syncronized with
 105   // eventual changes in the geometry.
 106   //
 107
 108   TString volPath;
 109   TString symName;
 110
 111   TString vpStr   = "ALIC_1/B077_1/BSEGMO";
 112   TString vpApp1  = "_1/BTRD";
 113   TString vpApp2  = "_1";
 114   TString vpApp3a = "/UTR1_1/UTS1_1/UTI1_1/UT";
 115   TString vpApp3b = "/UTR2_1/UTS2_1/UTI2_1/UT";
 116   TString vpApp3c = "/UTR3_1/UTS3_1/UTI3_1/UT";
 117
 118   TString snStr   = "TRD/sm";
 119   TString snApp1  = "/st";
 120   TString snApp2  = "/pl";
 121
 122   //
 123   // The super modules
 124   // The symbolic names are: TRD/sm00
 125   //                           ...
 126   //                         TRD/sm17
 127   //
 128   for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
 129
 130     volPath  = vpStr;
 131     volPath += isector;
 132     volPath += vpApp1;
 133     volPath += isector;
 134     volPath += vpApp2;
 135
 136     symName  = snStr;
 137     symName += Form("%02d",isector);
 138
 139     gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
 140
 141   }
 142
 143   //
 144   // The readout chambers
 145   // The symbolic names are: TRD/sm00/st0/pl0
 146   //                           ...
 147   //                         TRD/sm17/st4/pl5
 148   //
 149   AliGeomManager::ELayerID idTRD1 = AliGeomManager::kTRD1;
 150   Int_t layer, modUID;
 151
 152   for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
 153
 154     if (fGeometry->GetSMstatus(isector) == 0) continue;
 155
 156     for (Int_t istack = 0; istack < AliTRDgeometry::Nstack(); istack++) {
 157       for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) {
 158
 159         layer = idTRD1 + ilayer;
 160         modUID = AliGeomManager::LayerToVolUIDSafe(layer,isector*5+istack);
 161
 162         Int_t idet = AliTRDgeometry::GetDetectorSec(ilayer,istack);
 163
 164         volPath  = vpStr;
 165         volPath += isector;
 166         volPath += vpApp1;
 167         volPath += isector;
 168         volPath += vpApp2;
 169         switch (isector) {
 170         case 13:
 171         case 14:
 172         case 15:
 173           if (istack == 2) {
 174             continue;
 175           }
 176           volPath += vpApp3c;
 177           break;
 178         case 11:
 179         case 12:
 180           volPath += vpApp3b;
 181           break;
 182         default:
 183           volPath += vpApp3a;
 184         };
 185         volPath += Form("%02d",idet);
 186         volPath += vpApp2;
 187
 188         symName  = snStr;
 189         symName += Form("%02d",isector);
 190         symName += snApp1;
 191         symName += istack;
 192         symName += snApp2;
 193         symName += ilayer;
 194
 195         TGeoPNEntry *alignableEntry =
 196           gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data(),modUID);
 197
 198         // Add the tracking to local matrix following the TPC example
 199         if (alignableEntry) {
 200           TGeoHMatrix *globMatrix = alignableEntry->GetGlobalOrig();
 201           Double_t sectorAngle = 20.0 * (isector % 18) + 10.0;
 202           TGeoHMatrix *t2lMatrix  = new TGeoHMatrix();
 203           t2lMatrix->RotateZ(sectorAngle);
 204           t2lMatrix->MultiplyLeft(&(globMatrix->Inverse()));
 205           alignableEntry->SetMatrix(t2lMatrix);
 206         }
 207         else {
 208           AliError(Form("Alignable entry %s is not valid!",symName.Data()));
 209         }
 210
 211       }
 212     }
 213   }
 214
 215 }
 216
 217 //_____________________________________________________________________________
 218 void AliTRDv1::CreateGeometry()
 219 {
 220   //
 221   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
 222   // This version covers the full azimuth.
 223   //
 224
 225   // Check that FRAME is there otherwise we have no place where to put the TRD
 226   AliModule* frame = gAlice->GetModule("FRAME");
 227   if (!frame) {
 228     AliError("TRD needs FRAME to be present\n");
 229     return;
 230   }
 231
 232   // Define the chambers
 233   AliTRD::CreateGeometry();
 234
 235 }
 236
 237 //_____________________________________________________________________________
 238 void AliTRDv1::CreateMaterials()
 239 {
 240   //
 241   // Create materials for the Transition Radiation Detector version 1
 242   //
 243
 244   AliTRD::CreateMaterials();
 245
 246 }
 247
 248 //_____________________________________________________________________________
 249 void AliTRDv1::CreateTRhit(Int_t det)
 250 {
 251   //
 252   // Creates an electron cluster from a TR photon.
 253   // The photon is assumed to be created a the end of the radiator. The
 254   // distance after which it deposits its energy takes into account the
 255   // absorbtion of the entrance window and of the gas mixture in drift
 256   // volume.
 257   //
 258
 259   // Maximum number of TR photons per track
 260   const Int_t   kNTR         = 50;
 261
 262   TLorentzVector mom;
 263   TLorentzVector pos;
 264
 265   Float_t eTR[kNTR];
 266   Int_t   nTR;
 267
 268   // Create TR photons
 269   gMC->TrackMomentum(mom);
 270   Float_t pTot = mom.Rho();
 271   fTR->CreatePhotons(11,pTot,nTR,eTR);
 272   if (nTR > kNTR) {
 273     AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
 274   }
 275
 276   // Loop through the TR photons
 277   for (Int_t iTR = 0; iTR < nTR; iTR++) {
 278
 279     Float_t energyMeV = eTR[iTR] * 0.001;
 280     Float_t energyeV  = eTR[iTR] * 1000.0;
 281     Float_t absLength = 0.0;
 282     Float_t sigma     = 0.0;
 283
 284     // Take the absorbtion in the entrance window into account
 285     Double_t muMy = fTR->GetMuMy(energyMeV);
 286     sigma         = muMy * fFoilDensity;
 287     if (sigma > 0.0) {
 288       absLength = gRandom->Exp(1.0/sigma);
 289       if (absLength < AliTRDgeometry::MyThick()) {
 290         continue;
 291       }
 292     }
 293     else {
 294       continue;
 295     }
 296
 297     // The absorbtion cross sections in the drift gas
 298     // Gas-mixture (Xe/CO2)
 299     Double_t muNo = 0.0;
 300     if      (AliTRDCommonParam::Instance()->IsXenon()) {
 301       muNo = fTR->GetMuXe(energyMeV);
 302     }
 303     else if (AliTRDCommonParam::Instance()->IsArgon()) {
 304       muNo = fTR->GetMuAr(energyMeV);
 305     }
 306     Double_t muCO = fTR->GetMuCO(energyMeV);
 307     sigma = (fGasNobleFraction * muNo + (1.0 - fGasNobleFraction) * muCO)
 308           * fGasDensity
 309           * fTR->GetTemp();
 310
 311     // The distance after which the energy of the TR photon
 312     // is deposited.
 313     if (sigma > 0.0) {
 314       absLength = gRandom->Exp(1.0/sigma);
 315       if (absLength > (AliTRDgeometry::DrThick()
 316                      + AliTRDgeometry::AmThick())) {
 317         continue;
 318       }
 319     }
 320     else {
 321       continue;
 322     }
 323
 324     // The position of the absorbtion
 325     Float_t posHit[3];
 326     gMC->TrackPosition(pos);
 327     posHit[0] = pos[0] + mom[0] / pTot * absLength;
 328     posHit[1] = pos[1] + mom[1] / pTot * absLength;
 329     posHit[2] = pos[2] + mom[2] / pTot * absLength;
 330
 331     // Create the charge
 332     Int_t q = ((Int_t) (energyeV / fWion));
 333
 334     // Add the hit to the array. TR photon hits are marked
 335     // by negative charge
 336     AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 337           ,det
 338           ,posHit
 339           ,-q
 340           ,gMC->TrackTime()*1.0e06
 341           ,kTRUE);
 342
 343   }
 344
 345 }
 346
 347 //_____________________________________________________________________________
 348 void AliTRDv1::Init()
 349 {
 350   //
 351   // Initialise Transition Radiation Detector after geometry has been built.
 352   //
 353
 354   AliTRD::Init();
 355
 356   AliDebug(1,"Slow simulator\n");
 357
 358   // Switch on TR simulation as default
 359   if (!fTRon) {
 360     AliInfo("TR simulation off");
 361   }
 362   else {
 363     fTR = new AliTRDsimTR();
 364   }
 365
 366   AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
 367
 368 }
 369
 370 //_____________________________________________________________________________
 371 void AliTRDv1::StepManager()
 372 {
 373   //
 374   // Slow simulator. Every charged track produces electron cluster as hits
 375   // along its path across the drift volume. The step size is fixed in
 376   // this version of the step manager.
 377   //
 378   // Works for Xe/CO2 as well as Ar/CO2
 379   //
 380
 381   // PDG code electron
 382   const Int_t   kPdgElectron = 11;
 383
 384   Int_t    layer  = 0;
 385   Int_t    stack  = 0;
 386   Int_t    sector = 0;
 387   Int_t    det    = 0;
 388   Int_t    qTot;
 389
 390   Float_t  hits[3];
 391   Double_t eDep;
 392
 393   Bool_t   drRegion = kFALSE;
 394   Bool_t   amRegion = kFALSE;
 395
 396   TString  cIdPath;
 397   Char_t   cIdSector[3];
 398            cIdSector[2]  = 0;
 399
 400   TString  cIdCurrent;
 401   TString  cIdSensDr = "J";
 402   TString  cIdSensAm = "K";
 403   Char_t   cIdChamber[3];
 404            cIdChamber[2] = 0;
 405
 406   TLorentzVector pos;
 407   TLorentzVector mom;
 408
 409   const Int_t    kNlayer      = AliTRDgeometry::Nlayer();
 410   const Int_t    kNstack      = AliTRDgeometry::Nstack();
 411   const Int_t    kNdetsec     = kNlayer * kNstack;
 412
 413   const Double_t kBig         = 1.0e+12;
 414   const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment
 415
 416   // Set the maximum step size to a very large number for all
 417   // neutral particles and those outside the driftvolume
 418   if (!fPrimaryIonisation) gMC->SetMaxStep(kBig);
 419
 420   // If not charged track or already stopped or disappeared, just return.
 421   if ((!gMC->TrackCharge()) ||
 422         gMC->IsTrackDisappeared()) {
 423     return;
 424   }
 425
 426   // Inside a sensitive volume?
 427   cIdCurrent = gMC->CurrentVolName();
 428
 429   if (cIdSensDr == cIdCurrent[1]) {
 430     drRegion = kTRUE;
 431   }
 432   if (cIdSensAm == cIdCurrent[1]) {
 433     amRegion = kTRUE;
 434   }
 435
 436   if ((!drRegion) &&
 437       (!amRegion)) {
 438     return;
 439   }
 440
 441   // The hit coordinates and charge
 442   gMC->TrackPosition(pos);
 443   hits[0] = pos[0];
 444   hits[1] = pos[1];
 445   hits[2] = pos[2];
 446
 447   // The sector number (0 - 17), according to standard coordinate system
 448   cIdPath      = gGeoManager->GetPath();
 449   cIdSector[0] = cIdPath[21];
 450   cIdSector[1] = cIdPath[22];
 451   sector = atoi(cIdSector);
 452
 453   // The plane and chamber number
 454   cIdChamber[0]   = cIdCurrent[2];
 455   cIdChamber[1]   = cIdCurrent[3];
 456   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 457   stack = ((Int_t) idChamber / kNlayer);
 458   layer = ((Int_t) idChamber % kNlayer);
 459
 460   // The detector number
 461   det = fGeometry->GetDetector(layer,stack,sector);
 462
 463   // 0: InFlight 1:Entering 2:Exiting
 464   Int_t trkStat = 0;
 465
 466   // Special hits only in the drift region
 467   if      ((drRegion) &&
 468            (gMC->IsTrackEntering())) {
 469
 470     // Create a track reference at the entrance of each
 471     // chamber that contains the momentum components of the particle
 472     gMC->TrackMomentum(mom);
 473     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 474     trkStat = 1;
 475
 476     // Create the hits from TR photons if electron/positron is
 477     // entering the drift volume
 478     if ((fTR)   &&
 479         (fTRon) &&
 480         (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
 481       CreateTRhit(det);
 482     }
 483
 484   }
 485   else if ((amRegion) &&
 486            (gMC->IsTrackExiting())) {
 487
 488     // Create a track reference at the exit of each
 489     // chamber that contains the momentum components of the particle
 490     gMC->TrackMomentum(mom);
 491     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
 492     trkStat = 2;
 493
 494   }
 495
 496   // Calculate the charge according to GEANT Edep
 497   // Create a new dEdx hit
 498   eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
 499   qTot = (Int_t) (eDep / fWion);
 500   if ((qTot) ||
 501       (trkStat)) {
 502     AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 503           ,det
 504           ,hits
 505           ,qTot
 506           ,gMC->TrackTime()*1.0e06
 507           ,drRegion);
 508   }
 509
 510   // Set Maximum Step Size
 511   // Produce only one hit if Ekin is below cutoff
 512   if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) {
 513     return;
 514   }
 515   if (!fPrimaryIonisation) gMC->SetMaxStep(fStepSize);
 516
 517 }