LCOV - code coverage report
Current view: top level - TRD/TRDsim - AliTRDv1.cxx (source / functions) Hit Total Coverage
Test: coverage.info Lines: 199 219 90.9 %
Date: 2016-06-14 17:26:59 Functions: 14 14 100.0 %

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

Generated by: LCOV version 1.11