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

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  Transition Radiation Detector version 2 -- detailed simulation           //
//                                                                           //
//Begin_Html
/*
<img src="gif/AliTRDv2Class.gif">
*/
//End_Html
//                                                                           //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include <TMath.h>
#include <TRandom.h>
#include <TVector.h>
#include <TGeometry.h>
#include <TNode.h>
#include <TPGON.h> 

#include "GParticle.h"
#include "AliTRDv2.h"
#include "AliRun.h"
#include "AliConst.h"
#include "AliMC.h"

ClassImp(AliTRDv2)

//_____________________________________________________________________________
AliTRDv2::AliTRDv2(const char *name, const char *title) 
         :AliTRD(name, title) 
{
  //
  // Standard constructor for Transition Radiation Detector version 2
  //
  fIdSenO1 = fIdSenO2 = fIdSenO3 = fIdSenO4 = fIdSenO5 = fIdSenO6 = 0;
  fIdSenI1 = fIdSenI2 = fIdSenI3 = fIdSenI4 = fIdSenI5 = fIdSenI6 = 0;
  SetBufferSize(128000);
}
 
//_____________________________________________________________________________
void AliTRDv2::CreateGeometry()
{
  //
  // Create geometry for the Transition Radiation Detector version 2
  // This version covers the full azimuth. 
  //
  //Begin_Html
  /*
    <img src="gif/AliTRDv2.gif">
  */
  //End_Html
  //Begin_Html
  /*
    <img src="gif/AliTRDv2Tree.gif">
  */
  //End_Html


  // --- Name Conventions : 
  
  //     The mother volume and the support structure 
  //        TRD     --> Mother TRD volume                          (Air) 
  //        UTRS    --> The sectors of the detector                (Air) 
  //        UTSP    --> The aluminum part of the support structure (Al) 
  //        UTII(O) --> The inner parts of the support structure   (Air) 
  
  //     The chambers 
  //        UCI1-6  --> The frame of the inner chambers            (C) 
  //        UCO1-6  --> The frame of the outer chambers            (C) 
  //        UII1-6  --> The inner part of the inner chambers       (Air) 
  //        UIO1-6  --> The inner part of the outer chambers       (Air) 
  
  //     The layers inside a chamber 
  //        UT1I(O) --> Radiator layer                             (CO2) 
  //        UT2I(O) --> Polyethylene layer                         (PE) 
  //        UT3I(O) --> Mylar layer                                (Mylar) 
  //        UXI1-6  --> Xe/C02 layer in the inner chambers         (Xe/C02) 
  //        UXO1-6  --> Xe/C02 layer in the outer chambers         (Xe/C02) 
  //        UT5I(O) --> Cu layer (pads/sensitive)                  (Cu) 
  //        UT6I(O) --> Kapton layer                               (Kapton) 
  //        UT7I(O) --> NOMEX layer                                (C) 
  //        UT8I(O) --> Readout layer                              (Al) 
  
  
  // --- Contains geometry information 
  
  // --- Number of sectors in the full detector 
  // --- Number of modules in each sector 
  // --- z-Coordinates of the TRD-frame 
  // --- r-Coordinates of the TRD-frame 
  // --- Thickness of the aluminium of the support frame 
  // --- Thickness of the interior of the support frame 
  // --- Thickness of the carbon chamber frame 
  // --- Thickness and z-position of the PE-layer in the radiator 
  // --- Thickness and z-position of the radiator 
  // --- Thickness and z-position of the mylar-layer 
  // --- Thickness and z-position of the Xe/C02-layer 
  // --- Thickness and z-position of the Cu-layer (Pads) 
  // --- Thickness and z-position of the kapton-layer 
  // --- Thickness and z-position of the NOMEX-layer 
  //     Simple C-layer for the time being 
  // --- Thickness and z-position of the readout-layer 
  // --- Parameter for the arrays 
  // --- Number of different chambers 
  // --- Number of rotation matrices 

  AliMC* pMC = AliMC::GetMC();
  
  Float_t xpos, ypos, zpos;
  Int_t icham;
  Int_t idmat[2];
  Float_t widma, widmi;
  Float_t lendifc, widdifc, heightc;
  Float_t par_ic[3], par_oc[3], par_mo[10], par_sp[4], par_ch[3];
  
  Int_t *idtmed = gAlice->Idtmed();
  
  //************************************************************************
  
  //     Definition of Volumes 
  
  //************************************************************************
  
  const Int_t ncham = 6; //Number of different chambers
  
  widmi  = rmin*TMath::Sin(kPI/nsect);
  widma  = rmax*TMath::Sin(kPI/nsect);
  // --- Some parameter for the chambers 
  lendifc = (zmax1-zmax2)/nmodul;
  heightc = (rmax-rmin)/nmodul;
  widdifc = (widma - widmi)/nmodul;
  // --- Definition of the Mother volume for the TRD (Air) 
  par_mo[0] = 0.;
  par_mo[1] = 360.;
  par_mo[2] = nsect;
  par_mo[3] = 2.;
  par_mo[4] = -zmax1;
  par_mo[5] = rmin;
  par_mo[6] = rmax;
  par_mo[7] = zmax1;
  par_mo[8] = rmin;
  par_mo[9] = rmax;
  pMC->Gsvolu("TRD ", "PGON", idtmed[1301], par_mo, 10);
  // --- Divide the mother volume into sectors 
  pMC->Gsdvn("UTRS", "TRD ", nsect, 2);
  // --- Definition of the aluminum part of the support structure (Al) 
  par_sp[0] = widmi - inframe/2;
  par_sp[1] = widma - inframe/2;
  par_sp[2] = zmax1/2;
  par_sp[3] = (rmax-rmin)/2;
  pMC->Gsvolu("UTSP", "TRD1", idtmed[1300], par_sp, 4);
  // --- Definition of the inner part of the support structure (Air) 
  par_sp[0] = widmi - inframe/2 - alfram1/2;
  par_sp[1] = widma - inframe/2 - alfram1/2;
  par_sp[2] = zmax1/4 -alfram2/2;
  par_sp[3] = (rmax-rmin)/2;
  pMC->Gsvolu("UTII", "TRD1", idtmed[1301], par_sp, 4);
  pMC->Gsvolu("UTIO", "TRD1", idtmed[1301], par_sp, 4);
  // --- Definition of the chambers 
  char ctagc[5], ctagi[5], ctagx[5];
  for (icham = 1; icham <= ncham; ++icham) {
    // --- Carbon frame of the inner chambers (C) 
    par_ch[0] = widmi + (icham-1) * widdifc - (inframe+alfram1)/2;
    par_ch[1] = zmax1/4 -alfram2/2;
    par_ch[2] = heightc/2.;
    sprintf(ctagc,"UCI%1d",icham);
    pMC->Gsvolu(ctagc, "BOX ", idtmed[1306], par_ch, 3);
    // --- Inner part of the inner chambers (Air) 
    par_ch[0] -= ccframe;
    par_ch[1] -= ccframe;
    sprintf(ctagi,"UII%1d",icham);
    pMC->Gsvolu(ctagi, "BOX ", idtmed[1301], par_ch, 3);
    // --- Carbon frame of the outer chambers (C) 
    par_ch[0] = widmi + (icham - 1) * widdifc - 2.;
    par_ch[1] = (icham - 6) * lendifc / 2. + zmax1/4 -alfram2/2;
    par_ch[2] = heightc / 2.;
    sprintf(ctagc,"UCO%1d",icham);
    pMC->Gsvolu(ctagc, "BOX ", idtmed[1306], par_ch, 3);
    // --- Inner part of the outer chambers (Air) 
    par_ch[0] -= ccframe;
    par_ch[1] -= ccframe;
    sprintf(ctagi,"UIO%1d",icham);
    pMC->Gsvolu(ctagi, "BOX ", idtmed[1301], par_ch, 3);
  }
  // --- Definition of the layers in each inner chamber 
  par_ic[0] = -1.;
  par_ic[1] = -1.;
  // --- Radiator layer 
  par_ic[2] = rathick/2;;
  pMC->Gsvolu("UT1I", "BOX ", idtmed[1311], par_ic, 3);
  // --- Polyethylene layer 
  par_ic[2] = pethick/2;
  pMC->Gsvolu("UT2I", "BOX ", idtmed[1302], par_ic, 3);
  // --- Mylar layer 
  par_ic[2] = mythick/2;
  pMC->Gsvolu("UT3I", "BOX ", idtmed[1307], par_ic, 3);
  // --- Xe/CO2 layer 
  par_ic[2] = 1.8;
  for (icham = 1; icham <= 6; ++icham) {
    sprintf(ctagx,"UXI%1d",icham);
    pMC->Gsvolu(ctagx, "BOX ", idtmed[1308], par_ic, 3);
  }
  // --- Cu layer 
  par_ic[2] = cuthick/2;
  pMC->Gsvolu("UT5I", "BOX ", idtmed[1304], par_ic, 3);
  // --- Kapton layer 
  par_ic[2] = kathick/2;
  pMC->Gsvolu("UT6I", "BOX ", idtmed[1310], par_ic, 3);
  // --- NOMEX layer 
  par_ic[2] = nothick/2;
  pMC->Gsvolu("UT7I", "BOX ", idtmed[1309], par_ic, 3);
  // --- Read out layer 
  par_ic[2] = rothick/2;
  pMC->Gsvolu("UT8I", "BOX ", idtmed[1305], par_ic, 3);
  // --- Definition of the layers in each outer chamber 
  par_oc[0] = -1.;
  par_oc[1] = -1.;
  // --- Radiator layer 
  par_oc[2] = rathick/2;;
  pMC->Gsvolu("UT1O", "BOX ", idtmed[1311], par_oc, 3);
  // --- Polyethylene layer 
  par_oc[2] = pethick/2;
  pMC->Gsvolu("UT2O", "BOX ", idtmed[1302], par_oc, 3);
  // --- Mylar layer 
  par_oc[2] = mythick/2;
  pMC->Gsvolu("UT3O", "BOX ", idtmed[1307], par_oc, 3);
  // --- Xe/CO2 layer 
  par_oc[2] = xethick/2;
  for (icham = 1; icham <= ncham; ++icham) {
    sprintf(ctagx,"UXO%1d",icham);
    pMC->Gsvolu(ctagx, "BOX ", idtmed[1308], par_oc, 3);
  }
  // --- Cu layer 
  par_oc[2] = cuthick/2;
  pMC->Gsvolu("UT5O", "BOX ", idtmed[1304], par_oc, 3);
  // --- Kapton layer 
  par_oc[2] = kathick/2;
  pMC->Gsvolu("UT6O", "BOX ", idtmed[1310], par_oc, 3);
  // --- NOMEX layer 
  par_oc[2] = nothick/2;
  pMC->Gsvolu("UT7O", "BOX ", idtmed[1309], par_oc, 3);
  // --- Read out layer 
  par_oc[2] = rothick/2;
  pMC->Gsvolu("UT8O", "BOX ", idtmed[1305], par_oc, 3);
  //*************************************************************************
  
  //     Positioning of Volumes 
  
  //************************************************************************
  // --- The rotation matrices 
  AliMatrix(idmat[0], 90., 90., 180., 0., 90., 0.);
  AliMatrix(idmat[1], 90., 90.,   0., 0., 90., 0.);
  // --- Position of the layers in a chamber 
  for (icham = 1; icham <= ncham; ++icham) {
    // --- The inner chambers 
    sprintf(ctagi,"UII%1d",icham);
    sprintf(ctagx,"UXI%1d",icham);
    pMC->Gspos("UT8I", icham, ctagi, 0., 0., rozpos, 0, "ONLY");
    pMC->Gspos("UT7I", icham, ctagi, 0., 0., nozpos, 0, "ONLY");
    pMC->Gspos("UT6I", icham, ctagi, 0., 0., kazpos, 0, "ONLY");
    pMC->Gspos("UT5I", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
    pMC->Gspos(ctagx,  1,     ctagi, 0., 0., xezpos, 0, "ONLY");
    pMC->Gspos("UT3I", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
    pMC->Gspos("UT1I", icham, ctagi, 0., 0., razpos, 0, "ONLY");
    // --- The outer chambers 
    sprintf(ctagi,"UIO%d",icham);
    sprintf(ctagx,"UXO%d",icham);
    pMC->Gspos("UT8O", icham, ctagi, 0., 0., rozpos, 0, "ONLY");
    pMC->Gspos("UT7O", icham, ctagi, 0., 0., nozpos, 0, "ONLY");
    pMC->Gspos("UT6O", icham, ctagi, 0., 0., kazpos, 0, "ONLY");
    pMC->Gspos("UT5O", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
    pMC->Gspos(ctagx,  1,     ctagi, 0., 0., xezpos, 0, "ONLY");
    pMC->Gspos("UT3O", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
    pMC->Gspos("UT1O", icham, ctagi, 0., 0., razpos, 0, "ONLY");
  }
  pMC->Gspos("UT2I", 1, "UT1I", 0., 0., 0., 0, "ONLY");
  pMC->Gspos("UT2O", 1, "UT1O", 0., 0., 0., 0, "ONLY");
  // --- Position of the inner part of the chambers in the carbon-frames 
  for (icham = 1; icham <= ncham; ++icham) {
    // --- The inner chambers 
    sprintf(ctagi,"UII%1d",icham);
    sprintf(ctagc,"UCI%1d",icham);
    pMC->Gspos(ctagi, 1, ctagc, 0., 0., 0., 0, "ONLY");
    // --- The outer chambers 
    sprintf(ctagi,"UIO%1d",icham);
    sprintf(ctagc,"UCO%1d",icham);
    pMC->Gspos(ctagi, 1, ctagc, 0., 0., 0., 0, "ONLY");
  }
  // --- Position of the chambers in the full TRD-setup 
  for (icham = 1; icham <= ncham; ++icham) {
    // --- The inner chambers 
    xpos = 0.;
    ypos = 0.;
    zpos = (icham - .5) * heightc - (rmax-rmin)/2;
    sprintf(ctagc,"UCI%1d",icham);
    pMC->Gspos(ctagc, 1, "UTII", xpos, ypos, zpos, 0, "ONLY");
    // --- The outer chambers 
    xpos = 0.;
    ypos = 0. - (icham - 6) * lendifc / 2.;
    zpos = (icham - .5) * heightc - (rmax-rmin)/2;
    sprintf(ctagc,"UCO%1d",icham);
    pMC->Gspos(ctagc, 1, "UTIO", xpos, ypos, zpos, 0, "ONLY");
  }
  // --- Position of the inner parts of the support frame 
  xpos = 0.;
  ypos = zmax1/4;
  zpos = 0.;
  pMC->Gspos("UTIO", 1, "UTSP", xpos,-ypos, zpos, 0, "ONLY");
  pMC->Gspos("UTII", 1, "UTSP", xpos, ypos, zpos, 0, "ONLY");
  // --- Position of the support frame in the TRD-sectors 
  xpos = (rmax+rmin)/2;
  ypos = 0.;
  zpos = zmax1/2;
  pMC->Gspos("UTSP", 1, "UTRS", xpos, ypos, zpos, idmat[0], "ONLY");
  pMC->Gspos("UTSP", 2, "UTRS", xpos, ypos,-zpos, idmat[1], "ONLY");
  // --- Position of TRD mother volume in ALICE experiment 
  pMC->Gspos("TRD ", 1, "ALIC", 0., 0., 0., 0, "ONLY");
}

//_____________________________________________________________________________
void AliTRDv2::DrawModule()
{
  //
  // Draw a shaded view of the Transition Radiation Detector version 1
  //

  AliMC* pMC = AliMC::GetMC();
  
  // Set everything unseen
  pMC->Gsatt("*", "seen", -1);
  // 
  // Set ALIC mother transparent
  pMC->Gsatt("ALIC","SEEN",0);
  //
  // Set the volumes visible
  pMC->Gsatt("TRD","SEEN",0);
  pMC->Gsatt("UTRS","SEEN",0);
  pMC->Gsatt("UTSP","SEEN",0);
  pMC->Gsatt("UTII","SEEN",0);
  pMC->Gsatt("UTIO","SEEN",0);
  pMC->Gsatt("UCI1","SEEN",0);
  pMC->Gsatt("UII1","SEEN",0);
  pMC->Gsatt("UCO1","SEEN",0);
  pMC->Gsatt("UIO1","SEEN",0);
  pMC->Gsatt("UCI2","SEEN",0);
  pMC->Gsatt("UII2","SEEN",0);
  pMC->Gsatt("UCO2","SEEN",0);
  pMC->Gsatt("UIO2","SEEN",0);
  pMC->Gsatt("UCI3","SEEN",0);
  pMC->Gsatt("UII3","SEEN",0);
  pMC->Gsatt("UCO3","SEEN",0);
  pMC->Gsatt("UIO3","SEEN",0);
  pMC->Gsatt("UCI4","SEEN",0);
  pMC->Gsatt("UII4","SEEN",0);
  pMC->Gsatt("UCO4","SEEN",0);
  pMC->Gsatt("UIO4","SEEN",0);
  pMC->Gsatt("UCI5","SEEN",0);
  pMC->Gsatt("UII5","SEEN",0);
  pMC->Gsatt("UCO5","SEEN",0);
  pMC->Gsatt("UIO5","SEEN",0);
  pMC->Gsatt("UCI6","SEEN",0);
  pMC->Gsatt("UII6","SEEN",0);
  pMC->Gsatt("UCO6","SEEN",0);
  pMC->Gsatt("UIO6","SEEN",0);
  pMC->Gsatt("UT1I","SEEN",1);
  pMC->Gsatt("UXI1","SEEN",1);
  pMC->Gsatt("UXI2","SEEN",1);
  pMC->Gsatt("UXI3","SEEN",1);
  pMC->Gsatt("UXI4","SEEN",1);
  pMC->Gsatt("UXI5","SEEN",1);
  pMC->Gsatt("UXI6","SEEN",1);
  pMC->Gsatt("UT1O","SEEN",1);
  pMC->Gsatt("UXO1","SEEN",1);
  pMC->Gsatt("UXO2","SEEN",1);
  pMC->Gsatt("UXO3","SEEN",1);
  pMC->Gsatt("UXO4","SEEN",1);
  pMC->Gsatt("UXO5","SEEN",1);
  pMC->Gsatt("UXO6","SEEN",1);
  //
  pMC->Gdopt("hide", "on");
  pMC->Gdopt("shad", "on");
  pMC->Gsatt("*", "fill", 7);
  pMC->SetClipBox(".");
  pMC->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
  pMC->DefaultRange();
  pMC->Gdraw("alic", 40, 30, 0, 12, 9.4, .021, .021);
  pMC->Gdhead(1111, "Transition Radiation Detector Version 2");
  pMC->Gdman(18, 4, "MAN");
  pMC->Gdopt("hide", "off");
}

//_____________________________________________________________________________
void AliTRDv2::CreateMaterials()
{
  //
  // Create materials for the Transition Radiation Detector version 2
  //
  printf("TRD: Slow simulation with fine geometry\n");
  AliTRD::CreateMaterials();
}

//_____________________________________________________________________________
void AliTRDv2::Init() 
{
  //
  // Initialise Transition Radiation Detector after geometry has been built
  //
  AliTRD::Init();
  AliMC* pMC = AliMC::GetMC();
  fIdSenI1 = pMC->VolId("UXI1");
  fIdSenI2 = pMC->VolId("UXI2");
  fIdSenI3 = pMC->VolId("UXI3");
  fIdSenI4 = pMC->VolId("UXI4");
  fIdSenI5 = pMC->VolId("UXI5");
  fIdSenI6 = pMC->VolId("UXI6");
  
  fIdSenO1 = pMC->VolId("UXO1");
  fIdSenO2 = pMC->VolId("UXO2");
  fIdSenO3 = pMC->VolId("UXO3");
  fIdSenO4 = pMC->VolId("UXO4");
  fIdSenO5 = pMC->VolId("UXO5");
  fIdSenO6 = pMC->VolId("UXO6");
}

//_____________________________________________________________________________
void AliTRDv2::StepManager()
{
  //
  // Called at every step in the Transition Radiation Detector version 2
  //
  Int_t          idSens, icSens;
  Int_t          iPla, iCha, iSec;
  Int_t          iOut;
  Int_t          vol[3]; 
  Int_t          iPid;
  
  const Float_t  kBig   = 1.0E+12;
  
  Float_t        random[1];
  Float_t        charge;
  Float_t        betaGamma, pp;
  Float_t        aMass;
  Float_t        eLos, qTot;
  Float_t        hits[4];
  Float_t        mom[4];
  Float_t        pTot;
  
  TClonesArray  &lhits = *fHits;
  
  AliMC* pMC = AliMC::GetMC();
  
  // Ionization energy
  // taken from: Ionization Measurements in High Energy Physics, Springer
  const Double_t kWion    = 23.0E-9;
  // Maximum energy for e+ e- g for the step-size calculation
  const Float_t  kPTotMax = 0.002;
  // Plateau value of the energy-loss for electron in xenon
  // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
  const Float_t  kPlateau = 1.70;
  // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
  // taken from: Ionization Measurements in High Energy Physics, Springer
  const Float_t  kPrim    = 43.68;
  
  // Set the maximum step size to a very large number for all 
  // neutral particles and those outside the driftvolume
  pMC->SetMaxStep(kBig); 
  
  // Use only charged tracks 
  if (( pMC->TrackCharge()   ) &&
      (!pMC->TrackStop()     ) && 
      (!pMC->TrackDisappear())) {
    
    // Find the sensitive volume
    idSens = pMC->CurrentVol(0,icSens);
    iPla   = 0;
    iOut   = 0;
    if      (idSens == fIdSenI1) iPla = 1; 
    else if (idSens == fIdSenI2) iPla = 2;
    else if (idSens == fIdSenI3) iPla = 3;
    else if (idSens == fIdSenI4) iPla = 4;
    else if (idSens == fIdSenI5) iPla = 5;
    else if (idSens == fIdSenI6) iPla = 6;
    else if (idSens == fIdSenO1) {
      iPla = 1;
      iOut = 1; }
    else if (idSens == fIdSenO2) {
      iPla = 2;
      iOut = 1; }
    else if (idSens == fIdSenO3) {
      iPla = 3;
      iOut = 1; }
    else if (idSens == fIdSenO4) {
      iPla = 4;
      iOut = 1; }
    else if (idSens == fIdSenO5) {
      iPla = 5;
      iOut = 1; }
    else if (idSens == fIdSenO6) {
      iPla = 6;
      iOut = 1; }
    
    // Inside a sensitive volume?
    if (iPla) {
      
      // Calculate the energy loss 
      // ((1/E^2.0) distribution for the fluctuations)
      pMC->Rndm(random,1);
      eLos = Eloss(random[0]);
      
      // The amount of charge created
      qTot = (Float_t) ((Int_t) (eLos / kWion) + 1);
      
      // The sector number
      pMC->CurrentVolOff(5,0,iSec);
      
      // The chamber number
      pMC->CurrentVolOff(4,0,iCha);
      if (iOut) 
        iCha = (iCha - 1) * 3 + 1;
      else      
        iCha =  iCha + 1;
      
      vol[0]  = iSec;
      vol[1]  = iCha;
      vol[2]  = iPla;

      pMC->TrackPosition(hits);
      hits[3] = qTot;
      
      // Add the hit
      new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
      
      pMC->TrackMomentum(mom);
      pTot = mom[3];
      
      // New step size for electrons only if momentum is small enough
      iPid = pMC->TrackPid();
      if ( (iPid >  3) ||
	   ((iPid <= 3) && (pTot < kPTotMax))) {
        aMass     = pMC->TrackMass();
        betaGamma = pTot / aMass;
        pp        = kPrim * BetheBloch(betaGamma);
	// Take charge > 1 into account
        charge = pMC->TrackCharge();
        if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
      }
      // Electrons above 20 Mev/c are at the plateau
      else {
        pp = kPrim * kPlateau;
      }
      
      // Calculate the maximum step size for the next tracking step
      do 
        pMC->Rndm(random,1);
      while ((random[0] == 1.) || (random[0] == 0.));
      pMC->SetMaxStep( - TMath::Log(random[0]) / pp);
      
    }
  }
}
Commit	Line	Data
fe4da5cc	1	///////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// Transition Radiation Detector version 2 -- detailed simulation //
	4	// //
	5	//Begin_Html
	6	/*
	7	<img src="gif/AliTRDv2Class.gif">
	8	*/
	9	//End_Html
	10	// //
	11	// //
	12	///////////////////////////////////////////////////////////////////////////////
	13
	14	#include <TMath.h>
	15	#include <TRandom.h>
	16	#include <TVector.h>
	17	#include <TGeometry.h>
	18	#include <TNode.h>
	19	#include <TPGON.h>
	20
	21	#include "GParticle.h"
	22	#include "AliTRDv2.h"
	23	#include "AliRun.h"
	24	#include "AliConst.h"
	25	#include "AliMC.h"
	26
	27	ClassImp(AliTRDv2)
	28
	29	//_____________________________________________________________________________
	30	AliTRDv2::AliTRDv2(const char name, const char title)
	31	:AliTRD(name, title)
	32	{
	33	//
	34	// Standard constructor for Transition Radiation Detector version 2
	35	//
	36	fIdSenO1 = fIdSenO2 = fIdSenO3 = fIdSenO4 = fIdSenO5 = fIdSenO6 = 0;
	37	fIdSenI1 = fIdSenI2 = fIdSenI3 = fIdSenI4 = fIdSenI5 = fIdSenI6 = 0;
	38	SetBufferSize(128000);
	39	}
	40
	41	//_____________________________________________________________________________
	42	void AliTRDv2::CreateGeometry()
	43	{
	44	//
	45	// Create geometry for the Transition Radiation Detector version 2
	46	// This version covers the full azimuth.
	47	//
	48	//Begin_Html
	49	/*
	50	<img src="gif/AliTRDv2.gif">
	51	*/
	52	//End_Html
	53	//Begin_Html
	54	/*
	55	<img src="gif/AliTRDv2Tree.gif">
	56	*/
	57	//End_Html
	58
	59
	60	// --- Name Conventions :
	61
	62	// The mother volume and the support structure
	63	// TRD --> Mother TRD volume (Air)
	64	// UTRS --> The sectors of the detector (Air)
65	// UTSP --> The aluminum part of the support structure (Al)
66	// UTII(O) --> The inner parts of the support structure (Air)
67
68	// The chambers
69	// UCI1-6 --> The frame of the inner chambers (C)
70	// UCO1-6 --> The frame of the outer chambers (C)
71	// UII1-6 --> The inner part of the inner chambers (Air)
72	// UIO1-6 --> The inner part of the outer chambers (Air)
73
74	// The layers inside a chamber
75	// UT1I(O) --> Radiator layer (CO2)
76	// UT2I(O) --> Polyethylene layer (PE)
77	// UT3I(O) --> Mylar layer (Mylar)
78	// UXI1-6 --> Xe/C02 layer in the inner chambers (Xe/C02)
79	// UXO1-6 --> Xe/C02 layer in the outer chambers (Xe/C02)
80	// UT5I(O) --> Cu layer (pads/sensitive) (Cu)
81	// UT6I(O) --> Kapton layer (Kapton)
82	// UT7I(O) --> NOMEX layer (C)
83	// UT8I(O) --> Readout layer (Al)
84
85
86	// --- Contains geometry information
87
88	// --- Number of sectors in the full detector
89	// --- Number of modules in each sector
90	// --- z-Coordinates of the TRD-frame
91	// --- r-Coordinates of the TRD-frame
92	// --- Thickness of the aluminium of the support frame
93	// --- Thickness of the interior of the support frame
94	// --- Thickness of the carbon chamber frame
95	// --- Thickness and z-position of the PE-layer in the radiator
96	// --- Thickness and z-position of the radiator
97	// --- Thickness and z-position of the mylar-layer
98	// --- Thickness and z-position of the Xe/C02-layer
99	// --- Thickness and z-position of the Cu-layer (Pads)
100	// --- Thickness and z-position of the kapton-layer
101	// --- Thickness and z-position of the NOMEX-layer
102	// Simple C-layer for the time being
103	// --- Thickness and z-position of the readout-layer
104	// --- Parameter for the arrays
105	// --- Number of different chambers
106	// --- Number of rotation matrices
107
108	AliMC* pMC = AliMC::GetMC();
109
110	Float_t xpos, ypos, zpos;
111	Int_t icham;
112	Int_t idmat[2];
113	Float_t widma, widmi;
114	Float_t lendifc, widdifc, heightc;
115	Float_t par_ic[3], par_oc[3], par_mo[10], par_sp[4], par_ch[3];
116
117	Int_t *idtmed = gAlice->Idtmed();
118
119	//************************************************************************
120
121	// Definition of Volumes
122
123	//************************************************************************
124
125	const Int_t ncham = 6; //Number of different chambers
126
127	widmi = rmin*TMath::Sin(kPI/nsect);
128	widma = rmax*TMath::Sin(kPI/nsect);
129	// --- Some parameter for the chambers
130	lendifc = (zmax1-zmax2)/nmodul;
131	heightc = (rmax-rmin)/nmodul;
132	widdifc = (widma - widmi)/nmodul;
133	// --- Definition of the Mother volume for the TRD (Air)
134	par_mo[0] = 0.;
135	par_mo[1] = 360.;
136	par_mo[2] = nsect;
137	par_mo[3] = 2.;
138	par_mo[4] = -zmax1;
139	par_mo[5] = rmin;
140	par_mo[6] = rmax;
141	par_mo[7] = zmax1;
142	par_mo[8] = rmin;
143	par_mo[9] = rmax;
144	pMC->Gsvolu("TRD ", "PGON", idtmed[1301], par_mo, 10);
145	// --- Divide the mother volume into sectors
146	pMC->Gsdvn("UTRS", "TRD ", nsect, 2);
147	// --- Definition of the aluminum part of the support structure (Al)
148	par_sp[0] = widmi - inframe/2;
149	par_sp[1] = widma - inframe/2;
150	par_sp[2] = zmax1/2;
151	par_sp[3] = (rmax-rmin)/2;
152	pMC->Gsvolu("UTSP", "TRD1", idtmed[1300], par_sp, 4);
153	// --- Definition of the inner part of the support structure (Air)
154	par_sp[0] = widmi - inframe/2 - alfram1/2;
155	par_sp[1] = widma - inframe/2 - alfram1/2;
156	par_sp[2] = zmax1/4 -alfram2/2;
157	par_sp[3] = (rmax-rmin)/2;
158	pMC->Gsvolu("UTII", "TRD1", idtmed[1301], par_sp, 4);
159	pMC->Gsvolu("UTIO", "TRD1", idtmed[1301], par_sp, 4);
160	// --- Definition of the chambers
161	char ctagc[5], ctagi[5], ctagx[5];
162	for (icham = 1; icham <= ncham; ++icham) {
163	// --- Carbon frame of the inner chambers (C)
164	par_ch[0] = widmi + (icham-1) * widdifc - (inframe+alfram1)/2;
165	par_ch[1] = zmax1/4 -alfram2/2;
166	par_ch[2] = heightc/2.;
167	sprintf(ctagc,"UCI%1d",icham);
168	pMC->Gsvolu(ctagc, "BOX ", idtmed[1306], par_ch, 3);
169	// --- Inner part of the inner chambers (Air)
170	par_ch[0] -= ccframe;
171	par_ch[1] -= ccframe;
172	sprintf(ctagi,"UII%1d",icham);
173	pMC->Gsvolu(ctagi, "BOX ", idtmed[1301], par_ch, 3);
174	// --- Carbon frame of the outer chambers (C)
175	par_ch[0] = widmi + (icham - 1) * widdifc - 2.;
176	par_ch[1] = (icham - 6) * lendifc / 2. + zmax1/4 -alfram2/2;
177	par_ch[2] = heightc / 2.;
178	sprintf(ctagc,"UCO%1d",icham);
179	pMC->Gsvolu(ctagc, "BOX ", idtmed[1306], par_ch, 3);
180	// --- Inner part of the outer chambers (Air)
181	par_ch[0] -= ccframe;
182	par_ch[1] -= ccframe;
183	sprintf(ctagi,"UIO%1d",icham);
184	pMC->Gsvolu(ctagi, "BOX ", idtmed[1301], par_ch, 3);
185	}
186	// --- Definition of the layers in each inner chamber
187	par_ic[0] = -1.;
188	par_ic[1] = -1.;
189	// --- Radiator layer
190	par_ic[2] = rathick/2;;
191	pMC->Gsvolu("UT1I", "BOX ", idtmed[1311], par_ic, 3);
192	// --- Polyethylene layer
193	par_ic[2] = pethick/2;
194	pMC->Gsvolu("UT2I", "BOX ", idtmed[1302], par_ic, 3);
195	// --- Mylar layer
196	par_ic[2] = mythick/2;
197	pMC->Gsvolu("UT3I", "BOX ", idtmed[1307], par_ic, 3);
198	// --- Xe/CO2 layer
199	par_ic[2] = 1.8;
200	for (icham = 1; icham <= 6; ++icham) {
201	sprintf(ctagx,"UXI%1d",icham);
202	pMC->Gsvolu(ctagx, "BOX ", idtmed[1308], par_ic, 3);
203	}
204	// --- Cu layer
205	par_ic[2] = cuthick/2;
206	pMC->Gsvolu("UT5I", "BOX ", idtmed[1304], par_ic, 3);
207	// --- Kapton layer
208	par_ic[2] = kathick/2;
209	pMC->Gsvolu("UT6I", "BOX ", idtmed[1310], par_ic, 3);
210	// --- NOMEX layer
211	par_ic[2] = nothick/2;
212	pMC->Gsvolu("UT7I", "BOX ", idtmed[1309], par_ic, 3);
213	// --- Read out layer
214	par_ic[2] = rothick/2;
215	pMC->Gsvolu("UT8I", "BOX ", idtmed[1305], par_ic, 3);
216	// --- Definition of the layers in each outer chamber
217	par_oc[0] = -1.;
218	par_oc[1] = -1.;
219	// --- Radiator layer
220	par_oc[2] = rathick/2;;
221	pMC->Gsvolu("UT1O", "BOX ", idtmed[1311], par_oc, 3);
222	// --- Polyethylene layer
223	par_oc[2] = pethick/2;
224	pMC->Gsvolu("UT2O", "BOX ", idtmed[1302], par_oc, 3);
225	// --- Mylar layer
226	par_oc[2] = mythick/2;
227	pMC->Gsvolu("UT3O", "BOX ", idtmed[1307], par_oc, 3);
228	// --- Xe/CO2 layer
229	par_oc[2] = xethick/2;
230	for (icham = 1; icham <= ncham; ++icham) {
231	sprintf(ctagx,"UXO%1d",icham);
232	pMC->Gsvolu(ctagx, "BOX ", idtmed[1308], par_oc, 3);
233	}
234	// --- Cu layer
235	par_oc[2] = cuthick/2;
236	pMC->Gsvolu("UT5O", "BOX ", idtmed[1304], par_oc, 3);
237	// --- Kapton layer
238	par_oc[2] = kathick/2;
239	pMC->Gsvolu("UT6O", "BOX ", idtmed[1310], par_oc, 3);
240	// --- NOMEX layer
241	par_oc[2] = nothick/2;
242	pMC->Gsvolu("UT7O", "BOX ", idtmed[1309], par_oc, 3);
243	// --- Read out layer
244	par_oc[2] = rothick/2;
245	pMC->Gsvolu("UT8O", "BOX ", idtmed[1305], par_oc, 3);
246	//*************************************************************************
247
248	// Positioning of Volumes
249
250	//************************************************************************
251	// --- The rotation matrices
252	AliMatrix(idmat[0], 90., 90., 180., 0., 90., 0.);
253	AliMatrix(idmat[1], 90., 90., 0., 0., 90., 0.);
254	// --- Position of the layers in a chamber
255	for (icham = 1; icham <= ncham; ++icham) {
256	// --- The inner chambers
257	sprintf(ctagi,"UII%1d",icham);
258	sprintf(ctagx,"UXI%1d",icham);
259	pMC->Gspos("UT8I", icham, ctagi, 0., 0., rozpos, 0, "ONLY");
260	pMC->Gspos("UT7I", icham, ctagi, 0., 0., nozpos, 0, "ONLY");
261	pMC->Gspos("UT6I", icham, ctagi, 0., 0., kazpos, 0, "ONLY");
262	pMC->Gspos("UT5I", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
263	pMC->Gspos(ctagx, 1, ctagi, 0., 0., xezpos, 0, "ONLY");
264	pMC->Gspos("UT3I", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
265	pMC->Gspos("UT1I", icham, ctagi, 0., 0., razpos, 0, "ONLY");
266	// --- The outer chambers
267	sprintf(ctagi,"UIO%d",icham);
268	sprintf(ctagx,"UXO%d",icham);
269	pMC->Gspos("UT8O", icham, ctagi, 0., 0., rozpos, 0, "ONLY");
270	pMC->Gspos("UT7O", icham, ctagi, 0., 0., nozpos, 0, "ONLY");
271	pMC->Gspos("UT6O", icham, ctagi, 0., 0., kazpos, 0, "ONLY");
272	pMC->Gspos("UT5O", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
273	pMC->Gspos(ctagx, 1, ctagi, 0., 0., xezpos, 0, "ONLY");
274	pMC->Gspos("UT3O", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
275	pMC->Gspos("UT1O", icham, ctagi, 0., 0., razpos, 0, "ONLY");
276	}
277	pMC->Gspos("UT2I", 1, "UT1I", 0., 0., 0., 0, "ONLY");
278	pMC->Gspos("UT2O", 1, "UT1O", 0., 0., 0., 0, "ONLY");
279	// --- Position of the inner part of the chambers in the carbon-frames
280	for (icham = 1; icham <= ncham; ++icham) {
281	// --- The inner chambers
282	sprintf(ctagi,"UII%1d",icham);
283	sprintf(ctagc,"UCI%1d",icham);
284	pMC->Gspos(ctagi, 1, ctagc, 0., 0., 0., 0, "ONLY");
285	// --- The outer chambers
286	sprintf(ctagi,"UIO%1d",icham);
287	sprintf(ctagc,"UCO%1d",icham);
288	pMC->Gspos(ctagi, 1, ctagc, 0., 0., 0., 0, "ONLY");
289	}
290	// --- Position of the chambers in the full TRD-setup
291	for (icham = 1; icham <= ncham; ++icham) {
292	// --- The inner chambers
293	xpos = 0.;
294	ypos = 0.;
295	zpos = (icham - .5) * heightc - (rmax-rmin)/2;
296	sprintf(ctagc,"UCI%1d",icham);
297	pMC->Gspos(ctagc, 1, "UTII", xpos, ypos, zpos, 0, "ONLY");
298	// --- The outer chambers
299	xpos = 0.;
300	ypos = 0. - (icham - 6) * lendifc / 2.;
301	zpos = (icham - .5) * heightc - (rmax-rmin)/2;
302	sprintf(ctagc,"UCO%1d",icham);
303	pMC->Gspos(ctagc, 1, "UTIO", xpos, ypos, zpos, 0, "ONLY");
304	}
305	// --- Position of the inner parts of the support frame
306	xpos = 0.;
307	ypos = zmax1/4;
308	zpos = 0.;
309	pMC->Gspos("UTIO", 1, "UTSP", xpos,-ypos, zpos, 0, "ONLY");
310	pMC->Gspos("UTII", 1, "UTSP", xpos, ypos, zpos, 0, "ONLY");
311	// --- Position of the support frame in the TRD-sectors
312	xpos = (rmax+rmin)/2;
313	ypos = 0.;
314	zpos = zmax1/2;
315	pMC->Gspos("UTSP", 1, "UTRS", xpos, ypos, zpos, idmat[0], "ONLY");
316	pMC->Gspos("UTSP", 2, "UTRS", xpos, ypos,-zpos, idmat[1], "ONLY");
317	// --- Position of TRD mother volume in ALICE experiment
318	pMC->Gspos("TRD ", 1, "ALIC", 0., 0., 0., 0, "ONLY");
319	}
320
321	//_____________________________________________________________________________
05e51f55	322	void AliTRDv2::DrawModule()
fe4da5cc	323	{
	324	//
	325	// Draw a shaded view of the Transition Radiation Detector version 1
	326	//
	327
	328	AliMC* pMC = AliMC::GetMC();
	329
	330	// Set everything unseen
	331	pMC->Gsatt("*", "seen", -1);
	332	//
	333	// Set ALIC mother transparent
	334	pMC->Gsatt("ALIC","SEEN",0);
	335	//
	336	// Set the volumes visible
	337	pMC->Gsatt("TRD","SEEN",0);
	338	pMC->Gsatt("UTRS","SEEN",0);
	339	pMC->Gsatt("UTSP","SEEN",0);
	340	pMC->Gsatt("UTII","SEEN",0);
	341	pMC->Gsatt("UTIO","SEEN",0);
	342	pMC->Gsatt("UCI1","SEEN",0);
	343	pMC->Gsatt("UII1","SEEN",0);
	344	pMC->Gsatt("UCO1","SEEN",0);
	345	pMC->Gsatt("UIO1","SEEN",0);
	346	pMC->Gsatt("UCI2","SEEN",0);
	347	pMC->Gsatt("UII2","SEEN",0);
	348	pMC->Gsatt("UCO2","SEEN",0);
	349	pMC->Gsatt("UIO2","SEEN",0);
	350	pMC->Gsatt("UCI3","SEEN",0);
	351	pMC->Gsatt("UII3","SEEN",0);
	352	pMC->Gsatt("UCO3","SEEN",0);
	353	pMC->Gsatt("UIO3","SEEN",0);
	354	pMC->Gsatt("UCI4","SEEN",0);
	355	pMC->Gsatt("UII4","SEEN",0);
	356	pMC->Gsatt("UCO4","SEEN",0);
	357	pMC->Gsatt("UIO4","SEEN",0);
	358	pMC->Gsatt("UCI5","SEEN",0);
	359	pMC->Gsatt("UII5","SEEN",0);
	360	pMC->Gsatt("UCO5","SEEN",0);
	361	pMC->Gsatt("UIO5","SEEN",0);
	362	pMC->Gsatt("UCI6","SEEN",0);
	363	pMC->Gsatt("UII6","SEEN",0);
	364	pMC->Gsatt("UCO6","SEEN",0);
	365	pMC->Gsatt("UIO6","SEEN",0);
	366	pMC->Gsatt("UT1I","SEEN",1);
	367	pMC->Gsatt("UXI1","SEEN",1);
	368	pMC->Gsatt("UXI2","SEEN",1);
	369	pMC->Gsatt("UXI3","SEEN",1);
	370	pMC->Gsatt("UXI4","SEEN",1);
	371	pMC->Gsatt("UXI5","SEEN",1);
	372	pMC->Gsatt("UXI6","SEEN",1);
	373	pMC->Gsatt("UT1O","SEEN",1);
	374	pMC->Gsatt("UXO1","SEEN",1);
	375	pMC->Gsatt("UXO2","SEEN",1);
	376	pMC->Gsatt("UXO3","SEEN",1);
	377	pMC->Gsatt("UXO4","SEEN",1);
	378	pMC->Gsatt("UXO5","SEEN",1);
	379	pMC->Gsatt("UXO6","SEEN",1);
	380	//
	381	pMC->Gdopt("hide", "on");
	382	pMC->Gdopt("shad", "on");
	383	pMC->Gsatt("*", "fill", 7);
	384	pMC->SetClipBox(".");
	385	pMC->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
	386	pMC->DefaultRange();
387	pMC->Gdraw("alic", 40, 30, 0, 12, 9.4, .021, .021);
388	pMC->Gdhead(1111, "Transition Radiation Detector Version 2");
389	pMC->Gdman(18, 4, "MAN");
390	pMC->Gdopt("hide", "off");
391	}
392
393	//_____________________________________________________________________________
394	void AliTRDv2::CreateMaterials()
395	{
396	//
397	// Create materials for the Transition Radiation Detector version 2
398	//
399	printf("TRD: Slow simulation with fine geometry\n");
400	AliTRD::CreateMaterials();
401	}
402
403	//_____________________________________________________________________________
404	void AliTRDv2::Init()
405	{
406	//
407	// Initialise Transition Radiation Detector after geometry has been built
408	//
409	AliTRD::Init();
410	AliMC* pMC = AliMC::GetMC();
411	fIdSenI1 = pMC->VolId("UXI1");
412	fIdSenI2 = pMC->VolId("UXI2");
413	fIdSenI3 = pMC->VolId("UXI3");
414	fIdSenI4 = pMC->VolId("UXI4");
415	fIdSenI5 = pMC->VolId("UXI5");
416	fIdSenI6 = pMC->VolId("UXI6");
417
418	fIdSenO1 = pMC->VolId("UXO1");
419	fIdSenO2 = pMC->VolId("UXO2");
420	fIdSenO3 = pMC->VolId("UXO3");
421	fIdSenO4 = pMC->VolId("UXO4");
422	fIdSenO5 = pMC->VolId("UXO5");
423	fIdSenO6 = pMC->VolId("UXO6");
424	}
425
426	//_____________________________________________________________________________
427	void AliTRDv2::StepManager()
428	{
429	//
430	// Called at every step in the Transition Radiation Detector version 2
431	//
432	Int_t idSens, icSens;
433	Int_t iPla, iCha, iSec;
434	Int_t iOut;
435	Int_t vol[3];
436	Int_t iPid;
437
438	const Float_t kBig = 1.0E+12;
439
440	Float_t random[1];
441	Float_t charge;
442	Float_t betaGamma, pp;
443	Float_t aMass;
444	Float_t eLos, qTot;
445	Float_t hits[4];
446	Float_t mom[4];
447	Float_t pTot;
448
449	TClonesArray &lhits = *fHits;
450
451	AliMC* pMC = AliMC::GetMC();
452
453	// Ionization energy
454	// taken from: Ionization Measurements in High Energy Physics, Springer
455	const Double_t kWion = 23.0E-9;
456	// Maximum energy for e+ e- g for the step-size calculation
457	const Float_t kPTotMax = 0.002;
458	// Plateau value of the energy-loss for electron in xenon
459	// taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
460	const Float_t kPlateau = 1.70;
461	// dN1/dx\|min for the gas mixture (90% Xe + 10% CO2)
462	// taken from: Ionization Measurements in High Energy Physics, Springer
463	const Float_t kPrim = 43.68;
464
465	// Set the maximum step size to a very large number for all
466	// neutral particles and those outside the driftvolume
467	pMC->SetMaxStep(kBig);
468
469	// Use only charged tracks
470	if (( pMC->TrackCharge() ) &&
471	(!pMC->TrackStop() ) &&
472	(!pMC->TrackDisappear())) {
473
474	// Find the sensitive volume
475	idSens = pMC->CurrentVol(0,icSens);
476	iPla = 0;
477	iOut = 0;
478	if (idSens == fIdSenI1) iPla = 1;
479	else if (idSens == fIdSenI2) iPla = 2;
480	else if (idSens == fIdSenI3) iPla = 3;
481	else if (idSens == fIdSenI4) iPla = 4;
482	else if (idSens == fIdSenI5) iPla = 5;
483	else if (idSens == fIdSenI6) iPla = 6;
484	else if (idSens == fIdSenO1) {
485	iPla = 1;
486	iOut = 1; }
487	else if (idSens == fIdSenO2) {
488	iPla = 2;
489	iOut = 1; }
490	else if (idSens == fIdSenO3) {
491	iPla = 3;
492	iOut = 1; }
493	else if (idSens == fIdSenO4) {
494	iPla = 4;
495	iOut = 1; }
496	else if (idSens == fIdSenO5) {
497	iPla = 5;
498	iOut = 1; }
499	else if (idSens == fIdSenO6) {
500	iPla = 6;
501	iOut = 1; }
502
503	// Inside a sensitive volume?
504	if (iPla) {
505
506	// Calculate the energy loss
507	// ((1/E^2.0) distribution for the fluctuations)
508	pMC->Rndm(random,1);
509	eLos = Eloss(random[0]);
510
511	// The amount of charge created
512	qTot = (Float_t) ((Int_t) (eLos / kWion) + 1);
513
514	// The sector number
515	pMC->CurrentVolOff(5,0,iSec);
516
517	// The chamber number
518	pMC->CurrentVolOff(4,0,iCha);
519	if (iOut)
520	iCha = (iCha - 1) * 3 + 1;
521	else
522	iCha = iCha + 1;
523
524	vol[0] = iSec;
525	vol[1] = iCha;
526	vol[2] = iPla;
527
528	pMC->TrackPosition(hits);
529	hits[3] = qTot;
530
531	// Add the hit
532	new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
533
534	pMC->TrackMomentum(mom);
535	pTot = mom[3];
536
537	// New step size for electrons only if momentum is small enough
538	iPid = pMC->TrackPid();
539	if ( (iPid > 3) \|\|
540	((iPid <= 3) && (pTot < kPTotMax))) {
541	aMass = pMC->TrackMass();
542	betaGamma = pTot / aMass;
543	pp = kPrim * BetheBloch(betaGamma);
544	// Take charge > 1 into account
545	charge = pMC->TrackCharge();
546	if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
547	}
548	// Electrons above 20 Mev/c are at the plateau
549	else {
550	pp = kPrim * kPlateau;
551	}
552
553	// Calculate the maximum step size for the next tracking step
554	do
555	pMC->Rndm(random,1);
556	while ((random[0] == 1.) \|\| (random[0] == 0.));
557	pMC->SetMaxStep( - TMath::Log(random[0]) / pp);
558
559	}
560	}
561	}