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