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