[u/mrichter/AliRoot.git] / TRD / AliTRDv1.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
//  Transition Radiation Detector version 1 -- slow simulator             //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include <stdlib.h> 

#include <TLorentzVector.h>
#include <TMath.h>
#include <TRandom.h>
#include <TVector.h>
#include <TVirtualMC.h>
#include <TGeoManager.h>
#include <TGeoMatrix.h>
#include <TGeoPhysicalNode.h>

#include "AliConst.h"
#include "AliLog.h"
#include "AliTrackReference.h"
#include "AliMC.h"
#include "AliRun.h"
#include "AliGeomManager.h"

#include "AliTRDgeometry.h"
#include "AliTRDCommonParam.h"
#include "AliTRDhit.h"
#include "AliTRDsimTR.h"
#include "AliTRDv1.h"

ClassImp(AliTRDv1)
 
//_____________________________________________________________________________
AliTRDv1::AliTRDv1()
  :AliTRD()
  ,fTRon(kFALSE)
  ,fTR(NULL)
  ,fStepSize(0)
  ,fWion(0)
{
  //
  // Default constructor
  //

}

//_____________________________________________________________________________
AliTRDv1::AliTRDv1(const char *name, const char *title) 
  :AliTRD(name,title) 
  ,fTRon(kTRUE)
  ,fTR(NULL)
  ,fStepSize(0.1)
  ,fWion(0)
{
  //
  // Standard constructor for Transition Radiation Detector version 1
  //

  SetBufferSize(128000);

  if      (AliTRDCommonParam::Instance()->IsXenon()) {
    fWion = 23.53; // Ionization energy XeCO2 (85/15)
  }
  else if (AliTRDCommonParam::Instance()->IsArgon()) {
    fWion = 27.21; // Ionization energy ArCO2 (82/18)
  }
  else {
    AliFatal("Wrong gas mixture");
    exit(1);
  }

}

//_____________________________________________________________________________
AliTRDv1::~AliTRDv1()
{
  //
  // AliTRDv1 destructor
  //

  if (fTR) {
    delete fTR;
    fTR     = 0;
  }

}
 
//_____________________________________________________________________________
void AliTRDv1::AddAlignableVolumes() const
{
  //
  // Create entries for alignable volumes associating the symbolic volume
  // name with the corresponding volume path. Needs to be syncronized with
  // eventual changes in the geometry.
  //

  TString volPath;
  TString symName;

  TString vpStr   = "ALIC_1/B077_1/BSEGMO";
  TString vpApp1  = "_1/BTRD";
  TString vpApp2  = "_1";
  TString vpApp3a = "/UTR1_1/UTS1_1/UTI1_1/UT";
  TString vpApp3b = "/UTR2_1/UTS2_1/UTI2_1/UT";
  TString vpApp3c = "/UTR3_1/UTS3_1/UTI3_1/UT";

  TString snStr   = "TRD/sm";
  TString snApp1  = "/st";
  TString snApp2  = "/pl";

  //
  // The super modules
  // The symbolic names are: TRD/sm00
  //                           ...
  //                         TRD/sm17
  //
  for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {

    volPath  = vpStr;
    volPath += isector;
    volPath += vpApp1;
    volPath += isector;
    volPath += vpApp2;

    symName  = snStr;
    symName += Form("%02d",isector);

    gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());

  }

  //
  // The readout chambers
  // The symbolic names are: TRD/sm00/st0/pl0
  //                           ...
  //                         TRD/sm17/st4/pl5
  //
  AliGeomManager::ELayerID idTRD1 = AliGeomManager::kTRD1;
  Int_t layer, modUID;
  
  for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {

    if (fGeometry->GetSMstatus(isector) == 0) continue;

    for (Int_t istack = 0; istack < AliTRDgeometry::Nstack(); istack++) {
      for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) {

	layer = idTRD1 + ilayer;
	modUID = AliGeomManager::LayerToVolUIDSafe(layer,isector*5+istack);

        Int_t idet = AliTRDgeometry::GetDetectorSec(ilayer,istack);

        volPath  = vpStr;
        volPath += isector;
        volPath += vpApp1;
        volPath += isector;
        volPath += vpApp2;
        switch (isector) {
        case 13:
        case 14:
        case 15:
          if (istack == 2) {
            continue;
	  }
          volPath += vpApp3c;
          break;
        case 11:
        case 12:
          volPath += vpApp3b;
          break;
        default:
          volPath += vpApp3a;
	};
        volPath += Form("%02d",idet);
        volPath += vpApp2;

        symName  = snStr;
        symName += Form("%02d",isector);
        symName += snApp1;
        symName += istack;
        symName += snApp2;
        symName += ilayer;

        TGeoPNEntry *alignableEntry = 
	  gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data(),modUID);

	// Add the tracking to local matrix following the TPC example
	if (alignableEntry) {
	  TGeoHMatrix *globMatrix = alignableEntry->GetGlobalOrig();
	  Double_t sectorAngle = 20.0 * (isector % 18) + 10.0;
	  TGeoHMatrix *t2lMatrix  = new TGeoHMatrix();
	  t2lMatrix->RotateZ(sectorAngle);
	  t2lMatrix->MultiplyLeft(&(globMatrix->Inverse()));
	  alignableEntry->SetMatrix(t2lMatrix);
	}
	else {
	  AliError(Form("Alignable entry %s is not valid!",symName.Data()));
	}

      }
    }
  }

}

//_____________________________________________________________________________
void AliTRDv1::CreateGeometry()
{
  //
  // Create the GEANT geometry for the Transition Radiation Detector - Version 1
  // This version covers the full azimuth. 
  //

  // Check that FRAME is there otherwise we have no place where to put the TRD
  AliModule* frame = gAlice->GetModule("FRAME");
  if (!frame) {
    AliError("TRD needs FRAME to be present\n");
    return;
  }

  // Define the chambers
  AliTRD::CreateGeometry();

}

//_____________________________________________________________________________
void AliTRDv1::CreateMaterials()
{
  //
  // Create materials for the Transition Radiation Detector version 1
  //

  AliTRD::CreateMaterials();

}

//_____________________________________________________________________________
void AliTRDv1::CreateTRhit(Int_t det)
{
  //
  // Creates an electron cluster from a TR photon.
  // The photon is assumed to be created a the end of the radiator. The 
  // distance after which it deposits its energy takes into account the 
  // absorbtion of the entrance window and of the gas mixture in drift
  // volume.
  //

  // Maximum number of TR photons per track
  const Int_t   kNTR         = 50;

  TLorentzVector mom;
  TLorentzVector pos;

  Float_t eTR[kNTR];
  Int_t   nTR;

  // Create TR photons
  gMC->TrackMomentum(mom);
  Float_t pTot = mom.Rho();
  fTR->CreatePhotons(11,pTot,nTR,eTR);
  if (nTR > kNTR) {
    AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
  }

  // Loop through the TR photons
  for (Int_t iTR = 0; iTR < nTR; iTR++) {

    Float_t energyMeV = eTR[iTR] * 0.001;
    Float_t energyeV  = eTR[iTR] * 1000.0;
    Float_t absLength = 0.0;
    Float_t sigma     = 0.0;

    // Take the absorbtion in the entrance window into account
    Double_t muMy = fTR->GetMuMy(energyMeV);
    sigma         = muMy * fFoilDensity;
    if (sigma > 0.0) {
      absLength = gRandom->Exp(1.0/sigma);
      if (absLength < AliTRDgeometry::MyThick()) {
        continue;
      }
    }
    else {
      continue;
    }

    // The absorbtion cross sections in the drift gas
    // Gas-mixture (Xe/CO2)
    Double_t muNo = 0.0;
    if      (AliTRDCommonParam::Instance()->IsXenon()) {
      muNo = fTR->GetMuXe(energyMeV);
    }
    else if (AliTRDCommonParam::Instance()->IsArgon()) {
      muNo = fTR->GetMuAr(energyMeV);
    }
    Double_t muCO = fTR->GetMuCO(energyMeV);
    sigma = (fGasNobleFraction * muNo + (1.0 - fGasNobleFraction) * muCO) 
          * fGasDensity 
          * fTR->GetTemp();

    // The distance after which the energy of the TR photon
    // is deposited.
    if (sigma > 0.0) {
      absLength = gRandom->Exp(1.0/sigma);
      if (absLength > (AliTRDgeometry::DrThick()
                     + AliTRDgeometry::AmThick())) {
        continue;
      }
    }
    else {
      continue;
    }

    // The position of the absorbtion
    Float_t posHit[3];
    gMC->TrackPosition(pos);
    posHit[0] = pos[0] + mom[0] / pTot * absLength;
    posHit[1] = pos[1] + mom[1] / pTot * absLength;
    posHit[2] = pos[2] + mom[2] / pTot * absLength;

    // Create the charge 
    Int_t q = ((Int_t) (energyeV / fWion));

    // Add the hit to the array. TR photon hits are marked 
    // by negative charge
    AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
          ,det
          ,posHit
          ,-q
          ,gMC->TrackTime()*1.0e06
          ,kTRUE);

  }

}

//_____________________________________________________________________________
void AliTRDv1::Init() 
{
  //
  // Initialise Transition Radiation Detector after geometry has been built.
  //

  AliTRD::Init();

  AliDebug(1,"Slow simulator\n");

  // Switch on TR simulation as default
  if (!fTRon) {
    AliInfo("TR simulation off");
  }
  else {
    fTR = new AliTRDsimTR();
  }

  AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");

}

//_____________________________________________________________________________
void AliTRDv1::StepManager()
{
  //
  // Slow simulator. Every charged track produces electron cluster as hits 
  // along its path across the drift volume. The step size is fixed in
  // this version of the step manager.
  //
  // Works for Xe/CO2 as well as Ar/CO2
  //

  // PDG code electron
  const Int_t   kPdgElectron = 11;

  Int_t    layer  = 0;
  Int_t    stack  = 0;
  Int_t    sector = 0;
  Int_t    det    = 0;
  Int_t    qTot;

  Float_t  hits[3];
  Double_t eDep;

  Bool_t   drRegion = kFALSE;
  Bool_t   amRegion = kFALSE;

  TString  cIdPath;
  Char_t   cIdSector[3];
           cIdSector[2]  = 0;

  TString  cIdCurrent;
  TString  cIdSensDr = "J";
  TString  cIdSensAm = "K";
  Char_t   cIdChamber[3];
           cIdChamber[2] = 0;

  TLorentzVector pos;
  TLorentzVector mom;

  const Int_t    kNlayer      = AliTRDgeometry::Nlayer();
  const Int_t    kNstack      = AliTRDgeometry::Nstack();
  const Int_t    kNdetsec     = kNlayer * kNstack;

  const Double_t kBig         = 1.0e+12;
  const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment

  // Set the maximum step size to a very large number for all 
  // neutral particles and those outside the driftvolume
  gMC->SetMaxStep(kBig); 

  // If not charged track or already stopped or disappeared, just return.
  if ((!gMC->TrackCharge()) || 
        gMC->IsTrackDisappeared()) {
    return;
  }

  // Inside a sensitive volume?
  cIdCurrent = gMC->CurrentVolName();

  if (cIdSensDr == cIdCurrent[1]) {
    drRegion = kTRUE;
  }
  if (cIdSensAm == cIdCurrent[1]) {
    amRegion = kTRUE;
  }

  if ((!drRegion) && 
      (!amRegion)) {
    return;
  }

  // The hit coordinates and charge
  gMC->TrackPosition(pos);
  hits[0] = pos[0];
  hits[1] = pos[1];
  hits[2] = pos[2];

  // The sector number (0 - 17), according to standard coordinate system
  cIdPath      = gGeoManager->GetPath();
  cIdSector[0] = cIdPath[21];
  cIdSector[1] = cIdPath[22];
  sector = atoi(cIdSector);

  // The plane and chamber number
  cIdChamber[0]   = cIdCurrent[2];
  cIdChamber[1]   = cIdCurrent[3];
  Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
  stack = ((Int_t) idChamber / kNlayer);
  layer = ((Int_t) idChamber % kNlayer);

  // The detector number
  det = fGeometry->GetDetector(layer,stack,sector);

  // 0: InFlight 1:Entering 2:Exiting
  Int_t trkStat = 0;

  // Special hits only in the drift region
  if      ((drRegion) &&
           (gMC->IsTrackEntering())) {

    // Create a track reference at the entrance of each
    // chamber that contains the momentum components of the particle
    gMC->TrackMomentum(mom);
    AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
    trkStat = 1;

    // Create the hits from TR photons if electron/positron is
    // entering the drift volume
    if ((fTR) && 
        (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
      CreateTRhit(det);
    }

  }
  else if ((amRegion) && 
           (gMC->IsTrackExiting())) {

    // Create a track reference at the exit of each
    // chamber that contains the momentum components of the particle
    gMC->TrackMomentum(mom);
    AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
    trkStat = 2;

  }
  
  // Calculate the charge according to GEANT Edep
  // Create a new dEdx hit
  eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
  qTot = (Int_t) (eDep / fWion);
  if ((qTot) ||
      (trkStat)) {
    AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
          ,det
          ,hits
          ,qTot
          ,gMC->TrackTime()*1.0e06
          ,drRegion);
  }

  // Set Maximum Step Size
  // Produce only one hit if Ekin is below cutoff
  if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) {
    return;
  }
  gMC->SetMaxStep(fStepSize);

}
Commit	Line	Data
	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	}