TRD/AliTRDv2.cxx

   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 <TVector.h>
  16
  17 #include "AliTRDv2.h"
  18 #include "AliRun.h"
  19 #include "AliMC.h"
  20 #include "AliConst.h"
  21
  22 ClassImp(AliTRDv2)
  23
  24 //_____________________________________________________________________________
  25 AliTRDv2::AliTRDv2(const char *name, const char *title)
  26          :AliTRD(name, title)
  27 {
  28   //
  29   // Standard constructor for Transition Radiation Detector version 2
  30   //
  31   for (Int_t icham = 0; icham < ncham; ++icham) {
  32     fIdSensI[icham] = 0;
  33     fIdSensN[icham] = 0;
  34     fIdSensO[icham] = 0;
  35   }
  36   fDeltaE = NULL;
  37
  38   SetBufferSize(128000);
  39 }
  40
  41 AliTRDv2::~AliTRDv2()
  42 {
  43    if (fDeltaE)  delete fDeltaE;
  44 }
  45
  46 //_____________________________________________________________________________
  47 void AliTRDv2::CreateGeometry()
  48 {
  49   //
  50   // Create geometry for the Transition Radiation Detector version 2
  51   // This version covers the full azimuth.
  52   // --- Author :  Christoph Blume (GSI) 20/5/99
  53   //
  54   // --- Volume names :
  55   //       TRD         --> Mother TRD volume                              (Al)
  56   //       UTRS        --> Sectors of the sub-detector                    (Al)
  57   //       UTRI        --> Inner part of the detector frame               (Air)
  58   //     The chambers
  59   //       UCI1-6      --> The frame of the inner chambers                (C)
  60   //       UCN1-6      --> The frame of the neighbouring chambers         (C)
  61   //       UCO1-6      --> The frame of the outer chambers                (C)
  62   //       UII1-6      --> The inner part of the inner chambers           (Air)
  63   //       UIN1-6      --> The inner part of the neighbouring chambers    (Air)
  64   //       UIO1-6      --> The inner part of the outer chambers           (Air)
  65   //     The layers inside a chamber
  66   //       UT0I(N,O)   --> Radiator seal                                  (G10)
  67   //       UT1I(N,O)   --> Radiator                                       (CO2)
  68   //       UT2I(N,O)   --> Polyethylene of radiator                       (PE)
  69   //       UT3I(N,O)   --> Entrance window                                (Mylar)
  70   //       UXI(N,O)1-6 --> Gas volume (sensitive)                         (Xe/Isobutane)
  71   //       UT5I(N,O)   --> Pad plane                                      (Cu)
  72   //       UT6I(N,O)   --> Support structure                              (G10)
  73   //       UT7I(N,O)   --> FEE + signal lines                             (Cu)
  74   //       UT8I(N,O)   --> Polyethylene of cooling device                 (PE)
  75   //       UT9I(N,O)   --> Cooling water                                  (Water)
  76   //
  77   //Begin_Html
  78   /*
  79     <img src="gif/AliTRDv2.gif">
  80   */
  81   //End_Html
  82   //Begin_Html
  83   /*
  84     <img src="gif/AliTRDv2Tree.gif">
  85   */
  86   //End_Html
  87
  88   Float_t xpos, ypos, zpos;
  89   Int_t   idmat[2];
  90
  91   const Int_t nparmo = 10;
  92   const Int_t nparfr =  4;
  93   const Int_t nparch =  3;
  94   const Int_t nparic =  4;
  95   const Int_t nparnc =  4;
  96   const Int_t nparoc = 11;
  97
  98   Float_t par_mo[nparmo];
  99   Float_t par_fr[nparfr];
 100   Float_t par_ch[nparch];
 101   Float_t par_ic[nparic];
 102   Float_t par_nc[nparnc];
 103   Float_t par_oc[nparoc];
 104
 105   Int_t *idtmed = gAlice->Idtmed();
 106
 107   AliMC* pMC = AliMC::GetMC();
 108
 109   //////////////////////////////////////////////////////////////////////////
 110   //     Definition of Volumes
 111   //////////////////////////////////////////////////////////////////////////
 112
 113   // Definition of the mother volume for the TRD (Al)
 114   par_mo[0] =   0.;
 115   par_mo[1] = 360.;
 116   par_mo[2] = nsect;
 117   par_mo[3] = 2.;
 118   par_mo[4] = -zmax1;
 119   par_mo[5] = rmin;
 120   par_mo[6] = rmax;
 121   par_mo[7] =  zmax1;
 122   par_mo[8] = rmin;
 123   par_mo[9] = rmax;
 124   pMC->Gsvolu("TRD ", "PGON", idtmed[1301-1], par_mo, nparmo);
 125   pMC->Gsdvn("UTRS", "TRD ", nsect, 2);
 126
 127   // The minimal width of a sector in rphi-direction
 128   Float_t widmi = rmin * TMath::Sin(kPI/nsect);
 129   // The maximal width of a sector in rphi-direction
 130   Float_t widma = rmax * TMath::Sin(kPI/nsect);
 131   // The total thickness of the spaceframe (Al + Air)
 132   Float_t frame = widmi - (widpl1 / 2);
 133
 134   // Definition of the inner part of the detector frame (Air)
 135   par_fr[0] = widmi - alframe / 2.;
 136   par_fr[1] = widma - alframe / 2.;
 137   par_fr[2] = zmax1;
 138   par_fr[3] = (rmax - rmin) / 2;
 139   pMC->Gsvolu("UTRI", "TRD1", idtmed[1302-1], par_fr, nparfr);
 140
 141   // Some parameter for the chambers
 142   Float_t lendifc = (zmax1 - zmax2) / nmodul;
 143   Float_t heightc = (rmax  - rmin ) / nmodul;
 144   Float_t widdifc = (widma - widmi) / nmodul;
 145
 146   // Definition of the chambers
 147   Char_t ctagc[5], ctagi[5];
 148   for (Int_t icham = 1; icham <= ncham; ++icham) {
 149
 150     // Carbon frame of the inner chambers (C)
 151     par_ch[0] = widmi + (icham-1) * widdifc - frame;
 152     par_ch[1] = zleni   / 2.;
 153     par_ch[2] = heightc / 2.;
 154     sprintf(ctagc,"UCI%1d",icham);
 155     pMC->Gsvolu(ctagc, "BOX ", idtmed[1307-1], par_ch, nparch);
 156     // Inner part of the inner chambers (Air)
 157     par_ch[0] -= ccframe;
 158     par_ch[1] -= ccframe;
 159     sprintf(ctagc,"UII%1d",icham);
 160     pMC->Gsvolu(ctagc, "BOX ", idtmed[1302-1], par_ch, nparch);
 161
 162     // Carbon frame of the neighbouring chambers (C)
 163     par_ch[0] = widmi + (icham-1) * widdifc - frame;
 164     par_ch[1] = zlenn   / 2.;
 165     par_ch[2] = heightc / 2.;
 166     sprintf(ctagc,"UCN%1d",icham);
 167     pMC->Gsvolu(ctagc, "BOX ", idtmed[1307-1], par_ch, nparch);
 168     // Inner part of the neighbouring chambers (Air)
 169     par_ch[0] -= ccframe;
 170     par_ch[1] -= ccframe;
 171     sprintf(ctagc,"UIN%1d",icham);
 172     pMC->Gsvolu(ctagc, "BOX ", idtmed[1302-1], par_ch, nparch);
 173
 174     // Carbon frame of the outer chambers (C)
 175     par_ch[0] = widmi + (icham-1) * widdifc - frame;
 176     par_ch[1] = (icham - 6) * lendifc / 2. + zleno   / 2.;
 177     par_ch[2] = heightc / 2.;
 178     sprintf(ctagc,"UCO%1d",icham);
 179     pMC->Gsvolu(ctagc, "BOX ", idtmed[1307-1], par_ch, nparch);
 180     // Inner part of the outer chambers (Air)
 181     par_ch[0] -= ccframe;
 182     par_ch[1] -= ccframe;
 183     sprintf(ctagc,"UIO%1d",icham);
 184     pMC->Gsvolu(ctagc, "BOX ", idtmed[1302-1], par_ch, nparch);
 185
 186   }
 187
 188   // Definition of the layers in each inner chamber
 189   par_ic[0] = -1.;
 190   par_ic[1] = -1.;
 191   // G10 layer (radiator layer)
 192   par_ic[2] = sethick / 2;
 193   pMC->Gsvolu("UT0I", "BOX ", idtmed[1313-1], par_ic, nparic);
 194   // CO2 layer (radiator)
 195   par_ic[2] = rathick / 2;
 196   pMC->Gsvolu("UT1I", "BOX ", idtmed[1312-1], par_ic, nparic);
 197   // PE layer (radiator)
 198   par_ic[2] = pethick / 2;
 199   pMC->Gsvolu("UT2I", "BOX ", idtmed[1303-1], par_ic, nparic);
 200   // Mylar layer (entrance window + HV cathode)
 201   par_ic[2] = mythick / 2;
 202   pMC->Gsvolu("UT3I", "BOX ", idtmed[1308-1], par_ic, nparic);
 203   // Xe/Isobutane layer (gasvolume)
 204   par_ic[2] = xethick / 2.;
 205   for (Int_t icham = 1; icham <= 6; ++icham) {
 206     sprintf(ctagc,"UXI%1d",icham);
 207     pMC->Gsvolu(ctagc, "BOX ", idtmed[1309-1], par_ic, nparic);
 208   }
 209   // Cu layer (pad plane)
 210   par_ic[2] = cuthick / 2;
 211   pMC->Gsvolu("UT5I", "BOX ", idtmed[1305-1], par_ic, nparic);
 212   // G10 layer (support structure)
 213   par_ic[2] = suthick / 2;
 214   pMC->Gsvolu("UT6I", "BOX ", idtmed[1313-1], par_ic, nparic);
 215   // Cu layer (FEE + signal lines)
 216   par_ic[2] = fethick / 2;
 217   pMC->Gsvolu("UT7I", "BOX ", idtmed[1305-1], par_ic, nparic);
 218   // PE layer (cooling devices)
 219   par_ic[2] = cothick / 2;
 220   pMC->Gsvolu("UT8I", "BOX ", idtmed[1303-1], par_ic, nparic);
 221   // Water layer (cooling)
 222   par_ic[2] = wathick / 2;
 223   pMC->Gsvolu("UT9I", "BOX ", idtmed[1314-1], par_ic, nparic);
 224
 225   // Definition of the layers in each neighbouring chamber
 226   par_nc[0] = -1.;
 227   par_nc[1] = -1.;
 228   // G10 layer (radiator layer)
 229   par_nc[2] = sethick / 2;
 230   pMC->Gsvolu("UT0N", "BOX ", idtmed[1313-1], par_nc, nparnc);
 231   // CO2 layer (radiator)
 232   par_nc[2] = rathick / 2;
 233   pMC->Gsvolu("UT1N", "BOX ", idtmed[1312-1], par_nc, nparnc);
 234   // PE layer (radiator)
 235   par_nc[2] = pethick / 2;
 236   pMC->Gsvolu("UT2N", "BOX ", idtmed[1303-1], par_nc, nparnc);
 237   // Mylar layer (entrance window + HV cathode)
 238   par_nc[2] = mythick / 2;
 239   pMC->Gsvolu("UT3N", "BOX ", idtmed[1308-1], par_nc, nparnc);
 240   // Xe/Isobutane layer (gasvolume)
 241   par_nc[2] = xethick / 2.;
 242   for (Int_t icham = 1; icham <= 6; ++icham) {
 243     sprintf(ctagc,"UXN%1d",icham);
 244     pMC->Gsvolu(ctagc, "BOX ", idtmed[1309-1], par_nc, nparnc);
 245   }
 246   // Cu layer (pad plane)
 247   par_nc[2] = cuthick / 2;
 248   pMC->Gsvolu("UT5N", "BOX ", idtmed[1305-1], par_nc, nparnc);
 249   // G10 layer (support structure)
 250   par_nc[2] = suthick / 2;
 251   pMC->Gsvolu("UT6N", "BOX ", idtmed[1313-1], par_nc, nparnc);
 252   // Cu layer (FEE + signal lines)
 253   par_nc[2] = fethick / 2;
 254   pMC->Gsvolu("UT7N", "BOX ", idtmed[1305-1], par_nc, nparnc);
 255   // PE layer (cooling devices)
 256   par_nc[2] = cothick / 2;
 257   pMC->Gsvolu("UT8N", "BOX ", idtmed[1303-1], par_nc, nparnc);
 258   // Water layer (cooling)
 259   par_nc[2] = wathick / 2;
 260   pMC->Gsvolu("UT9N", "BOX ", idtmed[1314-1], par_nc, nparnc);
 261
 262   // Definition of the layers in each outer chamber
 263   par_oc[0] = -1.;
 264   par_oc[1] = -1.;
 265   // G10 layer (radiator layer)
 266   par_oc[2] = sethick / 2;
 267   pMC->Gsvolu("UT0O", "BOX ", idtmed[1313-1], par_oc, nparoc);
 268   // CO2 layer (radiator)
 269   par_oc[2] = rathick / 2;
 270   pMC->Gsvolu("UT1O", "BOX ", idtmed[1312-1], par_oc, nparoc);
 271   // PE layer (radiator)
 272   par_oc[2] = pethick / 2;
 273   pMC->Gsvolu("UT2O", "BOX ", idtmed[1303-1], par_oc, nparoc);
 274   // Mylar layer (entrance window + HV cathode)
 275   par_oc[2] = mythick / 2;
 276   pMC->Gsvolu("UT3O", "BOX ", idtmed[1308-1], par_oc, nparoc);
 277   // Xe/Isobutane layer (gasvolume)
 278   par_oc[2] = xethick / 2.;
 279   for (Int_t icham = 1; icham <= 6; ++icham) {
 280     sprintf(ctagc,"UXO%1d",icham);
 281     pMC->Gsvolu(ctagc, "BOX ", idtmed[1309-1], par_oc, nparoc);
 282   }
 283   // Cu layer (pad plane)
 284   par_oc[2] = cuthick / 2;
 285   pMC->Gsvolu("UT5O", "BOX ", idtmed[1305-1], par_oc, nparoc);
 286   // G10 layer (support structure)
 287   par_oc[2] = suthick / 2;
 288   pMC->Gsvolu("UT6O", "BOX ", idtmed[1313-1], par_oc, nparoc);
 289   // Cu layer (FEE + signal lines)
 290   par_oc[2] = fethick / 2;
 291   pMC->Gsvolu("UT7O", "BOX ", idtmed[1305-1], par_oc, nparoc);
 292   // PE layer (cooling devices)
 293   par_oc[2] = cothick / 2;
 294   pMC->Gsvolu("UT8O", "BOX ", idtmed[1303-1], par_oc, nparoc);
 295   // Water layer (cooling)
 296   par_oc[2] = wathick / 2;
 297   pMC->Gsvolu("UT9O", "BOX ", idtmed[1314-1], par_oc, nparoc);
 298
 299   //////////////////////////////////////////////////////////////////////////
 300   //     Positioning of Volumes
 301   //////////////////////////////////////////////////////////////////////////
 302
 303   // The rotation matrices
 304   AliMatrix(idmat[0],  90.,  90., 180.,   0.,  90.,   0.);
 305   AliMatrix(idmat[1],  90.,  90.,   0.,   0.,  90.,   0.);
 306
 307   // Position of the layers in a chamber
 308   pMC->Gspos("UT2I", 1, "UT1I", 0., 0., pezpos, 0, "ONLY");
 309   pMC->Gspos("UT2N", 1, "UT1N", 0., 0., pezpos, 0, "ONLY");
 310   pMC->Gspos("UT2O", 1, "UT1O", 0., 0., pezpos, 0, "ONLY");
 311   for (Int_t icham = 1; icham <= ncham; ++icham) {
 312     // The inner chambers
 313     sprintf(ctagi,"UII%1d",icham);
 314     sprintf(ctagc,"UXI%1d",icham);
 315     pMC->Gspos("UT9I", icham, ctagi, 0., 0., wazpos, 0, "ONLY");
 316     pMC->Gspos("UT8I", icham, ctagi, 0., 0., cozpos, 0, "ONLY");
 317     pMC->Gspos("UT7I", icham, ctagi, 0., 0., fezpos, 0, "ONLY");
 318     pMC->Gspos("UT6I", icham, ctagi, 0., 0., suzpos, 0, "ONLY");
 319     pMC->Gspos("UT5I", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
 320     pMC->Gspos(ctagc ,     1, ctagi, 0., 0., xezpos, 0, "ONLY");
 321     pMC->Gspos("UT3I", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
 322     pMC->Gspos("UT1I", icham, ctagi, 0., 0., razpos, 0, "ONLY");
 323     pMC->Gspos("UT0I", icham, ctagi, 0., 0., sezpos, 0, "ONLY");
 324     // The neighbouring chambers
 325     sprintf(ctagi,"UIN%1d",icham);
 326     sprintf(ctagc,"UXN%1d",icham);
 327     pMC->Gspos("UT9N", icham, ctagi, 0., 0., wazpos, 0, "ONLY");
 328     pMC->Gspos("UT8N", icham, ctagi, 0., 0., cozpos, 0, "ONLY");
 329     pMC->Gspos("UT7N", icham, ctagi, 0., 0., fezpos, 0, "ONLY");
 330     pMC->Gspos("UT6N", icham, ctagi, 0., 0., suzpos, 0, "ONLY");
 331     pMC->Gspos("UT5N", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
 332     pMC->Gspos(ctagc ,     1, ctagi, 0., 0., xezpos, 0, "ONLY");
 333     pMC->Gspos("UT3N", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
 334     pMC->Gspos("UT1N", icham, ctagi, 0., 0., razpos, 0, "ONLY");
 335     pMC->Gspos("UT0N", icham, ctagi, 0., 0., sezpos, 0, "ONLY");
 336     // The outer chambers
 337     sprintf(ctagi,"UIO%1d",icham);
 338     sprintf(ctagc,"UXO%1d",icham);
 339     pMC->Gspos("UT9O", icham, ctagi, 0., 0., wazpos, 0, "ONLY");
 340     pMC->Gspos("UT8O", icham, ctagi, 0., 0., cozpos, 0, "ONLY");
 341     pMC->Gspos("UT7O", icham, ctagi, 0., 0., fezpos, 0, "ONLY");
 342     pMC->Gspos("UT6O", icham, ctagi, 0., 0., suzpos, 0, "ONLY");
 343     pMC->Gspos("UT5O", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
 344     pMC->Gspos(ctagc ,     1, ctagi, 0., 0., xezpos, 0, "ONLY");
 345     pMC->Gspos("UT3O", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
 346     pMC->Gspos("UT1O", icham, ctagi, 0., 0., razpos, 0, "ONLY");
 347     pMC->Gspos("UT0O", icham, ctagi, 0., 0., sezpos, 0, "ONLY");
 348   }
 349
 350   // Position of the inner part of the chambers in the carbon-frames
 351   for (Int_t icham = 1; icham <= ncham; ++icham) {
 352     xpos = 0.;
 353     ypos = 0.;
 354     zpos = 0.;
 355     // The inner chambers
 356     sprintf(ctagi,"UII%1d",icham);
 357     sprintf(ctagc,"UCI%1d",icham);
 358     pMC->Gspos(ctagi, 1, ctagc, xpos, ypos, zpos, 0, "ONLY");
 359     // The neighbouring chambers
 360     sprintf(ctagi,"UIN%1d",icham);
 361     sprintf(ctagc,"UCN%1d",icham);
 362     pMC->Gspos(ctagi, 1, ctagc, xpos, ypos, zpos, 0, "ONLY");
 363     // The outer chambers
 364     sprintf(ctagi,"UIO%1d",icham);
 365     sprintf(ctagc,"UCO%1d",icham);
 366     pMC->Gspos(ctagi, 1, ctagc, xpos, ypos, zpos, 0, "ONLY");
 367   }
 368
 369   // Position of the chambers in the full TRD-setup
 370   for (Int_t icham = 1; icham <= ncham; ++icham) {
 371     // The inner chambers
 372     xpos = 0.;
 373     ypos = 0.;
 374     zpos = (icham-0.5) * heightc - (rmax - rmin) / 2;
 375     sprintf(ctagc,"UCI%1d",icham);
 376     pMC->Gspos(ctagc, 1, "UTRI", xpos, ypos, zpos, 0, "ONLY");
 377     // The neighbouring chambers
 378     xpos = 0.;
 379     ypos = (zleni + zlenn) / 2.;
 380     zpos = (icham-0.5) * heightc - (rmax - rmin) / 2;
 381     sprintf(ctagc,"UCN%1d",icham);
 382     pMC->Gspos(ctagc, 1, "UTRI", xpos, ypos, zpos, 0, "ONLY");
 383     ypos = -ypos;
 384     sprintf(ctagc,"UCN%1d",icham);
 385     pMC->Gspos(ctagc, 2, "UTRI", xpos, ypos, zpos, 0, "ONLY");
 386     // The outer chambers
 387     xpos = 0.;
 388     ypos = (zleni / 2. + zlenn + zmax2 + (icham-1) * lendifc) / 2.;
 389     zpos = (icham-0.5) * heightc - (rmax-rmin)/2;
 390     sprintf(ctagc,"UCO%1d",icham);
 391     pMC->Gspos(ctagc, 1, "UTRI", xpos, ypos, zpos, 0, "ONLY");
 392     ypos = -ypos;
 393     sprintf(ctagc,"UCO%1d",icham);
 394     pMC->Gspos(ctagc, 2, "UTRI", xpos, ypos, zpos, 0, "ONLY");
 395   }
 396
 397   // Position of the inner part of the detector frame
 398   xpos = (rmax + rmin) / 2;
 399   ypos = 0.;
 400   zpos = 0.;
 401   pMC->Gspos("UTRI", 1, "UTRS", xpos, ypos, zpos, idmat[0], "ONLY");
 402
 403   // Position of the TRD mother volume in the ALICE experiment
 404   xpos = 0.;
 405   ypos = 0.;
 406   zpos = 0.;
 407   pMC->Gspos("TRD ", 1, "ALIC", xpos, ypos, zpos,        0, "ONLY");
 408
 409 }
 410
 411 //_____________________________________________________________________________
 412 void AliTRDv2::DrawModule()
 413 {
 414   //
 415   // Draw a shaded view of the Transition Radiation Detector version 1
 416   //
 417
 418   AliMC* pMC = AliMC::GetMC();
 419
 420   // Set everything unseen
 421   pMC->Gsatt("*", "seen", -1);
 422
 423   // Set ALIC mother transparent
 424   pMC->Gsatt("ALIC","SEEN",0);
 425
 426   // Set the volumes visible
 427   pMC->Gsatt("TRD ","SEEN",0);
 428   pMC->Gsatt("UTRS","SEEN",0);
 429   pMC->Gsatt("UTRI","SEEN",0);
 430   Char_t ctag[5];
 431   for (Int_t icham = 0; icham < ncham; ++icham) {
 432     sprintf(ctag,"UCI%1d",icham+1);
 433     pMC->Gsatt(ctag,"SEEN",0);
 434     sprintf(ctag,"UCN%1d",icham+1);
 435     pMC->Gsatt(ctag,"SEEN",0);
 436     sprintf(ctag,"UCO%1d",icham+1);
 437     pMC->Gsatt(ctag,"SEEN",0);
 438     sprintf(ctag,"UII%1d",icham+1);
 439     pMC->Gsatt(ctag,"SEEN",0);
 440     sprintf(ctag,"UIN%1d",icham+1);
 441     pMC->Gsatt(ctag,"SEEN",0);
 442     sprintf(ctag,"UIO%1d",icham+1);
 443     pMC->Gsatt(ctag,"SEEN",0);
 444     sprintf(ctag,"UXI%1d",icham+1);
 445     pMC->Gsatt(ctag,"SEEN",1);
 446     sprintf(ctag,"UXN%1d",icham+1);
 447     pMC->Gsatt(ctag,"SEEN",1);
 448     sprintf(ctag,"UXO%1d",icham+1);
 449     pMC->Gsatt(ctag,"SEEN",1);
 450   }
 451   pMC->Gsatt("UT1I","SEEN",1);
 452   pMC->Gsatt("UT1N","SEEN",1);
 453   pMC->Gsatt("UT1O","SEEN",1);
 454
 455   pMC->Gdopt("hide", "on");
 456   pMC->Gdopt("shad", "on");
 457   pMC->Gsatt("*", "fill", 7);
 458   pMC->SetClipBox(".");
 459   pMC->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
 460   pMC->DefaultRange();
 461   pMC->Gdraw("alic", 40, 30, 0, 12, 9.4, .021, .021);
 462   pMC->Gdhead(1111, "Transition Radiation Detector Version 2");
 463   pMC->Gdman(18, 4, "MAN");
 464   pMC->Gdopt("hide", "off");
 465 }
 466
 467 //_____________________________________________________________________________
 468 void AliTRDv2::CreateMaterials()
 469 {
 470   //
 471   // Create materials for the Transition Radiation Detector version 2
 472   //
 473   AliTRD::CreateMaterials();
 474 }
 475
 476 //_____________________________________________________________________________
 477 void AliTRDv2::Init()
 478 {
 479   //
 480   // Initialise Transition Radiation Detector after geometry has been built
 481   //
 482
 483   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 484   const Float_t kPoti = 12.1;
 485   // Maximum energy (50 keV);
 486   const Float_t kEend = 50000.0;
 487
 488   AliTRD::Init();
 489
 490   AliMC* pMC = AliMC::GetMC();
 491
 492   // Get the sensitive volumes
 493   Char_t ctag[5];
 494   for (Int_t icham = 0; icham < ncham; ++icham) {
 495     sprintf(ctag,"UXI%1d",icham+1);
 496     fIdSensI[icham] = pMC->VolId(ctag);
 497     sprintf(ctag,"UXN%1d",icham+1);
 498     fIdSensN[icham] = pMC->VolId(ctag);
 499     sprintf(ctag,"UXO%1d",icham+1);
 500     fIdSensO[icham] = pMC->VolId(ctag);
 501   }
 502
 503   Float_t Poti = TMath::Log(kPoti);
 504   Float_t Eend = TMath::Log(kEend);
 505
 506   // Ermilova distribution for the delta-ray spectrum
 507   fDeltaE  = new TF1("deltae",Ermilova,Poti,Eend,0);
 508
 509 }
 510
 511 //_____________________________________________________________________________
 512 void AliTRDv2::StepManager()
 513 {
 514   //
 515   // Called at every step in the Transition Radiation Detector version 2
 516   //
 517
 518   Int_t          idSens, icSens, id;
 519   Int_t          iPla, iCha, iSec;
 520   Int_t          iOut;
 521   Int_t          vol[3];
 522   Int_t          iPid;
 523
 524   const Double_t kBig = 1.0E+12;
 525
 526   Float_t        hits[4];
 527   Float_t        mom[4];
 528   Float_t        random[1];
 529   Float_t        charge;
 530   Float_t        aMass;
 531
 532   Double_t       pTot;
 533   Double_t       qTot;
 534   Double_t       eDelta;
 535   Double_t       betaGamma, pp;
 536
 537   TClonesArray  &lhits = *fHits;
 538
 539   AliMC* pMC = AliMC::GetMC();
 540
 541   // Ionization energy
 542   const Float_t kWion    = 22.04;
 543   // Maximum energy for e+ e- g for the step-size calculation
 544   const Float_t kPTotMax = 0.002;
 545   // Plateau value of the energy-loss for electron in xenon
 546   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 547   //const Double_t kPlateau = 1.70;
 548   // the averaged value (26/3/99)
 549   const Float_t kPlateau = 1.55;
 550   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 551   const Float_t kPrim    = 48.0;
 552   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 553   const Float_t kPoti    = 12.1;
 554
 555   // Set the maximum step size to a very large number for all
 556   // neutral particles and those outside the driftvolume
 557   pMC->SetMaxStep(kBig);
 558
 559   // Use only charged tracks
 560   if (( pMC->TrackCharge()   ) &&
 561       (!pMC->TrackStop()     ) &&
 562       (!pMC->TrackDisappear())) {
 563
 564     // Find the sensitive volume
 565     idSens = pMC->CurrentVol(0,icSens);
 566     iPla   = 0;
 567     iOut   = 0;
 568     for (Int_t icham = 0; icham < ncham; ++icham) {
 569       if (idSens == fIdSensI[icham]) {
 570         iOut = 0;
 571         iPla = icham + 1;
 572       }
 573       if (idSens == fIdSensN[icham]) {
 574         iOut = 1;
 575         iPla = icham + 1;
 576       }
 577       if (idSens == fIdSensO[icham]) {
 578         iOut = 2;
 579         iPla = icham + 1;
 580       }
 581     }
 582
 583     // Inside a sensitive volume?
 584     if (iPla) {
 585
 586       // Calculate the energy of the delta-electrons
 587       eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 588       eDelta = TMath::Max(eDelta,0.0);
 589
 590       // The number of secondary electrons created
 591       qTot = (Double_t) ((Int_t) (eDelta / kWion) + 1);
 592
 593       // The sector number
 594       id = pMC->CurrentVolOff(4,0,iSec);
 595
 596       // The chamber number
 597       //   1: outer left
 598       //   2: neighbouring left
 599       //   3: inner
 600       //   4: neighbouring right
 601       //   5: outer right
 602       id = pMC->CurrentVolOff(2,0,iCha);
 603       if (iCha == 1)
 604         iCha = 3 + iOut;
 605       else
 606         iCha = 3 - iOut;
 607
 608       vol[0]  = iSec;
 609       vol[1]  = iCha;
 610       vol[2]  = iPla;
 611
 612       // Check on selected volumes
 613       Int_t addthishit = 1;
 614       if (fSensSelect) {
 615         if ((fSensPlane)   && (vol[2] != fSensPlane  )) addthishit = 0;
 616         if ((fSensChamber) && (vol[1] != fSensChamber)) addthishit = 0;
 617         if ((fSensSector)  && (vol[0] != fSensSector )) addthishit = 0;
 618       }
 619
 620       if (addthishit) {
 621
 622         // Add this hit
 623         pMC->TrackPosition(hits);
 624         hits[3] = qTot;
 625         new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
 626
 627         // The energy loss according to Bethe Bloch
 628         pMC->TrackMomentum(mom);
 629         pTot = mom[3];
 630         iPid = pMC->TrackPid();
 631         if ( (iPid >  3) ||
 632             ((iPid <= 3) && (pTot < kPTotMax))) {
 633           aMass     = pMC->TrackMass();
 634           betaGamma = pTot / aMass;
 635           pp        = kPrim * BetheBloch(betaGamma);
 636           // Take charge > 1 into account
 637           charge = pMC->TrackCharge();
 638           if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 639         }
 640         // Electrons above 20 Mev/c are at the plateau
 641         else {
 642           pp = kPrim * kPlateau;
 643         }
 644
 645         // Calculate the maximum step size for the next tracking step
 646         if (pp > 0) {
 647           do
 648             pMC->Rndm(random,1);
 649           while ((random[0] == 1.) || (random[0] == 0.));
 650           pMC->SetMaxStep( - TMath::Log(random[0]) / pp);
 651         }
 652
 653       }
 654       else {
 655         // set step size to maximal value
 656         pMC->SetMaxStep(kBig);
 657       }
 658
 659     }
 660
 661   }
 662
 663 }
 664
 665 //_____________________________________________________________________________
 666 Double_t AliTRDv2::BetheBloch(Double_t bg)
 667 {
 668   //
 669   // Parametrization of the Bethe-Bloch-curve
 670   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
 671   //
 672
 673   // The parameters have been adjusted to Xe-data found in:
 674   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 675   //const Double_t kP1 = 0.76176E-1;
 676   //const Double_t kP2 = 10.632;
 677   //const Double_t kP3 = 3.17983E-6;
 678   //const Double_t kP4 = 1.8631;
 679   //const Double_t kP5 = 1.9479;
 680
 681   // This parameters have been adjusted to averaged values from GEANT
 682   const Double_t kP1 = 7.17960e-02;
 683   const Double_t kP2 = 8.54196;
 684   const Double_t kP3 = 1.38065e-06;
 685   const Double_t kP4 = 5.30972;
 686   const Double_t kP5 = 2.83798;
 687
 688   if (bg > 0) {
 689     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
 690     Double_t aa = TMath::Power(yy,kP4);
 691     Double_t bb = TMath::Power((1./bg),kP5);
 692              bb = TMath::Log(kP3 + bb);
 693     return ((kP2 - aa - bb)*kP1 / aa);
 694   }
 695   else
 696     return 0;
 697
 698 }
 699
 700 //_____________________________________________________________________________
 701 Double_t Ermilova(Double_t *x, Double_t *par)
 702 {
 703   //
 704   // Calculates the delta-ray energy distribution according to Ermilova
 705   // Logarithmic scale !
 706   //
 707
 708   Double_t energy;
 709   Double_t dpos;
 710   Double_t dnde;
 711
 712   Int_t    pos1, pos2;
 713
 714   const Int_t nV = 31;
 715
 716   Float_t vxe[nV] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
 717                     , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
 718                     , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
 719                     , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
 720                     , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
 721                     , 9.4727, 9.9035,10.3735,10.5966,10.8198
 722                     ,11.5129 };
 723
 724   Float_t vye[nV] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
 725                     , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
 726                     ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
 727                     ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
 728                     ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
 729                     ,  0.04 ,  0.023,  0.015,  0.011,  0.01
 730                     ,  0.004 };
 731
 732   energy = x[0];
 733
 734   // Find the position
 735   pos1 = pos2 = 0;
 736   dpos = 0;
 737   do {
 738     dpos = energy - vxe[pos2++];
 739   }
 740   while (dpos > 0);
 741   pos2--;
 742   if (pos2 > nV) pos2 = nV;
 743   pos1 = pos2 - 1;
 744
 745   // Differentiate between the sampling points
 746   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
 747
 748   return dnde;
 749
 750 }