TRD/AliTRD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Transition Radiation Detector                                            //
  21 //  This class contains the basic functions for the Transition Radiation     //
  22 //  Detector.                                                                //
  23 //                                                                           //
  24 ///////////////////////////////////////////////////////////////////////////////
  25
  26 #include <stdlib.h>
  27 #include <Riostream.h>
  28
  29 #include <TFile.h>
  30 #include <TGeometry.h>
  31 #include <TLorentzVector.h>
  32 #include <TMath.h>
  33 #include <TNode.h>
  34 #include <TPGON.h>
  35 #include <TParticle.h>
  36 #include <TROOT.h>
  37 #include <TTree.h>
  38 #include <TVirtualMC.h>
  39
  40 #include "AliConst.h"
  41 #include "AliDigit.h"
  42 #include "AliLoader.h"
  43 #include "AliMagF.h"
  44 #include "AliRun.h"
  45 #include "AliTRD.h"
  46 #include "AliTRDdigit.h"
  47 #include "AliTRDdigitizer.h"
  48 #include "AliTRDdigitsManager.h"
  49 #include "AliTRDgeometryFull.h"
  50 #include "AliTRDgeometryHole.h"
  51 #include "AliTRDhit.h"
  52 #include "AliTRDpoints.h"
  53 #include "AliTRDtrackHits.h"
  54 #include "AliTRDrawData.h"
  55 #include "AliTrackReference.h"
  56 #include "AliMC.h"
  57
  58 ClassImp(AliTRD)
  59
  60 //_____________________________________________________________________________
  61 AliTRD::AliTRD()
  62 {
  63   //
  64   // Default constructor
  65   //
  66
  67   fIshunt        = 0;
  68   fGasMix        = 0;
  69   fHits          = 0;
  70   fDigits        = 0;
  71
  72   fGeometry      = 0;
  73
  74   fGasDensity    = 0;
  75   fFoilDensity   = 0;
  76
  77   fDrawTR        = 0;
  78   fDisplayType   = 0;
  79
  80   fTrackHits     = 0;
  81   fHitType       = 0;
  82
  83 }
  84
  85 //_____________________________________________________________________________
  86 AliTRD::AliTRD(const char *name, const char *title)
  87        : AliDetector(name,title)
  88 {
  89   //
  90   // Standard constructor for the TRD
  91   //
  92
  93   // Check that FRAME is there otherwise we have no place where to
  94   // put TRD
  95   AliModule* frame = gAlice->GetModule("FRAME");
  96   if (!frame) {
  97     Error("Ctor","TRD needs FRAME to be present\n");
  98     exit(1);
  99   }
 100
 101   // Define the TRD geometry
 102   if ((frame->IsVersion() == 0) ||
 103       (frame->IsVersion() == 1)) {
 104     fGeometry = new AliTRDgeometryFull();
 105   }
 106   else {
 107     Error("Ctor","Could not find valid FRAME version\n");
 108     exit(1);
 109   }
 110
 111   // Save the geometry
 112   TDirectory* saveDir = gDirectory;
 113   gAlice->GetRunLoader()->CdGAFile();
 114   fGeometry->Write("TRDgeometry");
 115   saveDir->cd();
 116
 117   // Allocate the hit array
 118   fHits           = new TClonesArray("AliTRDhit"     ,405);
 119   gAlice->GetMCApp()->AddHitList(fHits);
 120
 121   // Allocate the digits array
 122   fDigits         = 0;
 123
 124   fIshunt        = 0;
 125   fGasMix        = 1;
 126
 127   fGasDensity    = 0;
 128   fFoilDensity   = 0;
 129
 130   fDrawTR        = 0;
 131   fDisplayType   = 0;
 132
 133   fTrackHits     = 0;
 134   fHitType       = 2;
 135
 136   SetMarkerColor(kWhite);
 137
 138 }
 139
 140 //_____________________________________________________________________________
 141 AliTRD::AliTRD(const AliTRD &trd):AliDetector(trd)
 142 {
 143   //
 144   // Copy constructor
 145   //
 146
 147   ((AliTRD &) trd).Copy(*this);
 148
 149 }
 150
 151 //_____________________________________________________________________________
 152 AliTRD::~AliTRD()
 153 {
 154   //
 155   // TRD destructor
 156   //
 157
 158   fIshunt = 0;
 159
 160   if (fGeometry) {
 161     delete fGeometry;
 162     fGeometry  = 0;
 163   }
 164   if (fHits) {
 165     delete fHits;
 166     fHits      = 0;
 167   }
 168 //  if (fRecPoints) {
 169 //    delete fRecPoints;
 170 //    fRecPoints = 0;
 171 //  }
 172   if (fTrackHits) {
 173     delete fTrackHits;
 174     fTrackHits = 0;
 175   }
 176
 177 }
 178
 179 //_____________________________________________________________________________
 180 void AliTRD::Hits2Digits()
 181 {
 182   //
 183   // Create digits
 184   //
 185   AliTRDdigitizer digitizer("TRDdigitizer","TRD digitizer class");
 186   digitizer.SetDebug(GetDebug());
 187
 188   // Initialization
 189   digitizer.InitDetector();
 190
 191   if (!fLoader->TreeH()) fLoader->LoadHits("read");
 192   fLoader->LoadDigits("recreate");
 193   AliRunLoader* runLoader = fLoader->GetRunLoader();
 194
 195   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
 196     digitizer.Open(runLoader->GetFileName().Data(), iEvent);
 197     digitizer.MakeDigits();
 198     digitizer.WriteDigits();
 199   }
 200
 201   fLoader->UnloadHits();
 202   fLoader->UnloadDigits();
 203 }
 204
 205 //_____________________________________________________________________________
 206 void AliTRD::Hits2SDigits()
 207 {
 208   //
 209   // Create summable digits
 210   //
 211   AliTRDdigitizer digitizer("TRDdigitizer","TRD digitizer class");
 212   // For the summable digits
 213   digitizer.SetSDigits(kTRUE);
 214   digitizer.SetDebug(GetDebug());
 215
 216   // Initialization
 217   digitizer.InitDetector();
 218
 219   if (!fLoader->TreeH()) fLoader->LoadHits("read");
 220   fLoader->LoadSDigits("recreate");
 221   AliRunLoader* runLoader = fLoader->GetRunLoader();
 222
 223   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
 224     digitizer.Open(runLoader->GetFileName().Data(), iEvent);
 225     digitizer.MakeDigits();
 226     digitizer.WriteDigits();
 227   }
 228
 229   fLoader->UnloadHits();
 230   fLoader->UnloadSDigits();
 231   }
 232
 233 //_____________________________________________________________________________
 234 AliDigitizer* AliTRD::CreateDigitizer(AliRunDigitizer* manager) const
 235 {
 236   return new AliTRDdigitizer(manager);
 237 }
 238
 239 //_____________________________________________________________________________
 240 void AliTRD::SDigits2Digits()
 241 {
 242   //
 243   // Create final digits from summable digits
 244   //
 245
 246    // Create the TRD digitizer
 247   AliTRDdigitizer digitizer("TRDdigitizer","TRD digitizer class");
 248   digitizer.SetDebug(GetDebug());
 249
 250   // Set the parameter
 251   digitizer.SetEvent(gAlice->GetEvNumber());
 252
 253   // Initialization
 254   digitizer.InitDetector();
 255
 256   // Read the s-digits via digits manager
 257   AliTRDdigitsManager sdigitsManager;
 258
 259   sdigitsManager.SetDebug(GetDebug());
 260   sdigitsManager.SetSDigits(kTRUE);
 261   sdigitsManager.CreateArrays();
 262
 263   if (!fLoader->TreeS())
 264     if (fLoader->LoadSDigits("read"))
 265      {
 266        Error("SDigits2Digits","Error while reading SDigits for event %d",gAlice->GetEvNumber());
 267        return;
 268      }
 269   if (!fLoader->TreeS()) return;
 270
 271   sdigitsManager.ReadDigits(fLoader->TreeS());
 272
 273   // Add the s-digits to the input list
 274   digitizer.AddSDigitsManager(&sdigitsManager);
 275
 276   // Convert the s-digits to normal digits
 277   digitizer.SDigits2Digits();
 278
 279   // Store the digits
 280   if (!fLoader->TreeD()) fLoader->MakeTree("D");
 281   if (digitizer.MakeBranch(fLoader->TreeD())){
 282     digitizer.WriteDigits();
 283   }
 284
 285 }
 286
 287 //_____________________________________________________________________________
 288 void AliTRD::Digits2Raw()
 289 {
 290   //
 291   // convert digits of the current event to raw data
 292   //
 293
 294   fLoader->LoadDigits();
 295   TTree* digits = fLoader->TreeD();
 296   if (!digits) {
 297     Error("Digits2Raw", "no digits tree");
 298     return;
 299   }
 300
 301   AliTRDrawData rawWriter;
 302 //  rawWriter.SetDebug(2);
 303   if (!rawWriter.Digits2Raw(digits)) {
 304     Error("AliTRD::Digits2Raw","The raw writer could not load the digits tree");
 305   }
 306
 307   fLoader->UnloadDigits();
 308
 309 }
 310
 311 //_____________________________________________________________________________
 312 void AliTRD::AddHit(Int_t track, Int_t det, Float_t *hits, Int_t q
 313                   , Bool_t inDrift)
 314 {
 315   //
 316   // Add a hit for the TRD
 317   //
 318   // The data structure is set according to fHitType:
 319   //   bit0: standard TClonesArray
 320   //   bit1: compressed trackHits structure
 321   //
 322
 323   if (fHitType & 1) {
 324     TClonesArray &lhits = *fHits;
 325     new(lhits[fNhits++]) AliTRDhit(fIshunt,track,det,hits,q);
 326   }
 327
 328   if (fHitType > 1) {
 329     AddHit2(track,det,hits,q,inDrift);
 330   }
 331
 332 }
 333
 334 //_____________________________________________________________________________
 335 void AliTRD::BuildGeometry()
 336 {
 337   //
 338   // Create the ROOT TNode geometry for the TRD
 339   //
 340
 341   TNode *node, *top;
 342   TPGON *pgon;
 343
 344   Float_t rmin, rmax;
 345   Float_t zmax1, zmax2;
 346
 347   Int_t   iPlan;
 348
 349   const Int_t kColorTRD = 46;
 350
 351   // Find the top node alice
 352   top = gAlice->GetGeometry()->GetNode("alice");
 353
 354   if      (fDisplayType == 0) {
 355
 356     pgon = new TPGON("S_TRD","TRD","void",0,360,AliTRDgeometry::Nsect(),4);
 357     rmin = AliTRDgeometry::Rmin();
 358     rmax = AliTRDgeometry::Rmax();
 359     pgon->DefineSection(0,-AliTRDgeometry::Zmax1(),rmax,rmax);
 360     pgon->DefineSection(1,-AliTRDgeometry::Zmax2(),rmin,rmax);
 361     pgon->DefineSection(2, AliTRDgeometry::Zmax2(),rmin,rmax);
 362     pgon->DefineSection(3, AliTRDgeometry::Zmax1(),rmax,rmax);
 363     top->cd();
 364     node = new TNode("TRD","TRD","S_TRD",0,0,0,"");
 365     node->SetLineColor(kColorTRD);
 366     fNodes->Add(node);
 367
 368   }
 369   else if (fDisplayType == 1) {
 370
 371     Char_t name[7];
 372
 373     Float_t slope = (AliTRDgeometry::Zmax1() - AliTRDgeometry::Zmax2())
 374                   / (AliTRDgeometry::Rmax()  - AliTRDgeometry::Rmin());
 375
 376     rmin  = AliTRDgeometry::Rmin() + AliTRDgeometry::CraHght();
 377     rmax  = rmin                   + AliTRDgeometry::CdrHght();
 378
 379     Float_t thickness = rmin - AliTRDgeometry::Rmin();
 380     zmax2 = AliTRDgeometry::Zmax2() + slope * thickness;
 381     zmax1 = zmax2 + slope * AliTRDgeometry::DrThick();
 382
 383     for (iPlan = 0; iPlan < AliTRDgeometry::Nplan(); iPlan++) {
 384
 385       sprintf(name,"S_TR1%d",iPlan);
 386       pgon  = new TPGON(name,"TRD","void",0,360,AliTRDgeometry::Nsect(),4);
 387       pgon->DefineSection(0,-zmax1,rmax,rmax);
 388       pgon->DefineSection(1,-zmax2,rmin,rmax);
 389       pgon->DefineSection(2, zmax2,rmin,rmax);
 390       pgon->DefineSection(3, zmax1,rmax,rmax);
 391       top->cd();
 392       node = new TNode("TRD","TRD",name,0,0,0,"");
 393       node->SetLineColor(kColorTRD);
 394       fNodes->Add(node);
 395
 396       Float_t height = AliTRDgeometry::Cheight() + AliTRDgeometry::Cspace();
 397       rmin  = rmin  + height;
 398       rmax  = rmax  + height;
 399       zmax1 = zmax1 + slope * height;
 400       zmax2 = zmax2 + slope * height;
 401
 402     }
 403
 404     thickness += AliTRDgeometry::DrThick();
 405     rmin  = AliTRDgeometry::Rmin() + thickness;
 406     rmax  = rmin + AliTRDgeometry::AmThick();
 407     zmax2 = AliTRDgeometry::Zmax2() + slope * thickness;
 408     zmax1 = zmax2 + slope * AliTRDgeometry::AmThick();
 409
 410     for (iPlan = 0; iPlan < AliTRDgeometry::Nplan(); iPlan++) {
 411
 412       sprintf(name,"S_TR2%d",iPlan);
 413       pgon  = new TPGON(name,"TRD","void",0,360,AliTRDgeometry::Nsect(),4);
 414       pgon->DefineSection(0,-zmax1,rmax,rmax);
 415       pgon->DefineSection(1,-zmax2,rmin,rmax);
 416       pgon->DefineSection(2, zmax2,rmin,rmax);
 417       pgon->DefineSection(3, zmax1,rmax,rmax);
 418       top->cd();
 419       node = new TNode("TRD","TRD",name,0,0,0,"");
 420       node->SetLineColor(kColorTRD);
 421       fNodes->Add(node);
 422
 423       Float_t height = AliTRDgeometry::Cheight() + AliTRDgeometry::Cspace();
 424       rmin  = rmin  + height;
 425       rmax  = rmax  + height;
 426       zmax1 = zmax1 + slope * height;
 427       zmax2 = zmax2 + slope * height;
 428
 429     }
 430
 431   }
 432
 433 }
 434
 435 //_____________________________________________________________________________
 436 void AliTRD::Copy(TObject &trd)
 437 {
 438   //
 439   // Copy function
 440   //
 441
 442   ((AliTRD &) trd).fGasMix      = fGasMix;
 443   ((AliTRD &) trd).fGeometry    = fGeometry;
 444   ((AliTRD &) trd).fGasDensity  = fGasDensity;
 445   ((AliTRD &) trd).fFoilDensity = fFoilDensity;
 446   ((AliTRD &) trd).fDrawTR      = fDrawTR;
 447   ((AliTRD &) trd).fDisplayType = fDisplayType;
 448   ((AliTRD &) trd).fHitType     = fHitType;
 449
 450   //AliDetector::Copy(trd);
 451
 452 }
 453
 454 //_____________________________________________________________________________
 455 void AliTRD::CreateGeometry()
 456 {
 457   //
 458   // Creates the volumes for the TRD chambers
 459   //
 460
 461   // Check that FRAME is there otherwise we have no place where to put the TRD
 462   AliModule* frame = gAlice->GetModule("FRAME");
 463   if (!frame) {
 464     printf(" The TRD needs the FRAME to be defined first\n");
 465     return;
 466   }
 467
 468   fGeometry->CreateGeometry(fIdtmed->GetArray() - 1299);
 469
 470 }
 471
 472 //_____________________________________________________________________________
 473 void AliTRD::CreateMaterials()
 474 {
 475   //
 476   // Create the materials for the TRD
 477   // Origin Y.Foka
 478   //
 479
 480   Int_t   isxfld = gAlice->Field()->Integ();
 481   Float_t sxmgmx = gAlice->Field()->Max();
 482
 483   // For polyethilene (CH2)
 484   Float_t ape[2] = { 12., 1. };
 485   Float_t zpe[2] = {  6., 1. };
 486   Float_t wpe[2] = {  1., 2. };
 487   Float_t dpe    = 0.95;
 488
 489   // For mylar (C5H4O2)
 490   Float_t amy[3] = { 12., 1., 16. };
 491   Float_t zmy[3] = {  6., 1.,  8. };
 492   Float_t wmy[3] = {  5., 4.,  2. };
 493   Float_t dmy    = 1.39;
 494
 495   // For CO2
 496   Float_t aco[2] = { 12., 16. };
 497   Float_t zco[2] = {  6.,  8. };
 498   Float_t wco[2] = {  1.,  2. };
 499   Float_t dco    = 0.001977;
 500
 501   // For water
 502   Float_t awa[2] = {  1., 16. };
 503   Float_t zwa[2] = {  1.,  8. };
 504   Float_t wwa[2] = {  2.,  1. };
 505   Float_t dwa    = 1.0;
 506
 507   // For isobutane (C4H10)
 508   Float_t ais[2] = { 12.,  1. };
 509   Float_t zis[2] = {  6.,  1. };
 510   Float_t wis[2] = {  4., 10. };
 511   Float_t dis    = 0.00267;
 512
 513   // For plexiglas (C5H8O2)
 514   Float_t apg[3] = { 12.011 ,  1.0    , 15.9994 };
 515   Float_t zpg[3] = {  6.0   ,  1.0    ,  8.0    };
 516   Float_t wpg[3] = {  5.0   ,  8.0    ,  2.0    };
 517   Float_t dpg    = 1.18;
 518
 519   // For Xe/CO2-gas-mixture
 520   // Xe-content of the Xe/CO2-mixture (85% / 15%)
 521   Float_t fxc    = .85;
 522   // Xe-content of the Xe/Isobutane-mixture (97% / 3%)
 523   Float_t fxi    = .97;
 524   Float_t dxe    = .005858;
 525
 526   // General tracking parameter
 527   Float_t tmaxfd = -10.;
 528   Float_t stemax = -1e10;
 529   Float_t deemax = -0.1;
 530   Float_t epsil  =  1e-4;
 531   Float_t stmin  = -0.001;
 532
 533   Float_t absl, radl, d, buf[1];
 534   Float_t agm[2], zgm[2], wgm[2];
 535   Float_t dgm1, dgm2;
 536   Int_t   nbuf;
 537
 538   //////////////////////////////////////////////////////////////////////////
 539   //     Define Materials
 540   //////////////////////////////////////////////////////////////////////////
 541
 542   AliMaterial( 1, "Al"   ,  26.98, 13.0, 2.7     ,     8.9 ,    37.2);
 543   AliMaterial( 2, "Air"  ,  14.61,  7.3, 0.001205, 30420.0 , 67500.0);
 544   AliMaterial( 4, "Xe"   , 131.29, 54.0, dxe     ,  1447.59,     0.0);
 545   AliMaterial( 5, "Cu"   ,  63.54, 29.0, 8.96    ,     1.43,    14.8);
 546   AliMaterial( 6, "C"    ,  12.01,  6.0, 2.265   ,    18.8 ,    74.4);
 547   AliMaterial(12, "G10"  ,  20.00, 10.0, 1.7     ,    19.4 ,   999.0);
 548   AliMaterial(15, "Sn"   , 118.71, 50.0, 7.31    ,     1.21,    14.8);
 549   AliMaterial(16, "Si"   ,  28.09, 14.0, 2.33    ,     9.36,    37.2);
 550   AliMaterial(17, "Epoxy",  17.75,  8.9, 1.8     ,    21.82,   999.0);
 551
 552   // Mixtures
 553   AliMixture(3, "Polyethilene",   ape, zpe, dpe, -2, wpe);
 554   AliMixture(7, "Mylar",          amy, zmy, dmy, -3, wmy);
 555   AliMixture(8, "CO2",            aco, zco, dco, -2, wco);
 556   AliMixture(9, "Isobutane",      ais, zis, dis, -2, wis);
 557   AliMixture(13,"Water",          awa, zwa, dwa, -2, wwa);
 558   AliMixture(14,"Plexiglas",      apg, zpg, dpg, -3, wpg);
 559
 560   // Gas mixtures
 561   Char_t namate[21]="";
 562   // Xe/CO2-mixture
 563   // Get properties of Xe
 564   gMC->Gfmate((*fIdmate)[4], namate, agm[0], zgm[0], d, radl, absl, buf, nbuf);
 565   // Get properties of CO2
 566   gMC->Gfmate((*fIdmate)[8], namate, agm[1], zgm[1], d, radl, absl, buf, nbuf);
 567   // Create gas mixture
 568   wgm[0] = fxc;
 569   wgm[1] = 1. - fxc;
 570   dgm1   = wgm[0] * dxe + wgm[1] * dco;
 571   AliMixture(10, "Gas mixture 1", agm, zgm, dgm1,  2, wgm);
 572   // Xe/Isobutane-mixture
 573   // Get properties of Xe
 574   gMC->Gfmate((*fIdmate)[4], namate, agm[0], zgm[0], d, radl, absl, buf, nbuf);
 575   // Get properties of Isobutane
 576   gMC->Gfmate((*fIdmate)[9], namate, agm[1], zgm[1], d, radl, absl, buf, nbuf);
 577   // Create gas mixture
 578   wgm[0] = fxi;
 579   wgm[1] = 1. - fxi;
 580   dgm2   = wgm[0] * dxe + wgm[1] * dis;
 581   AliMixture(11, "Gas mixture 2", agm, zgm, dgm2,  2, wgm);
 582
 583   //////////////////////////////////////////////////////////////////////////
 584   //     Tracking Media Parameters
 585   //////////////////////////////////////////////////////////////////////////
 586
 587   // Al Frame
 588   AliMedium(1, "Al Frame",   1, 0, isxfld, sxmgmx
 589                 , tmaxfd, stemax, deemax, epsil, stmin);
 590   // Air
 591   AliMedium(2, "Air",        2, 0, isxfld, sxmgmx
 592                 , tmaxfd, stemax, deemax, epsil, stmin);
 593   // Polyethilene
 594   AliMedium(3, "Radiator",   3, 0, isxfld, sxmgmx
 595                 , tmaxfd, stemax, deemax, epsil, stmin);
 596   // Xe
 597   AliMedium(4, "Xe",         4, 1, isxfld, sxmgmx
 598                 , tmaxfd, stemax, deemax, epsil, stmin);
 599   // Cu pads
 600   AliMedium(5, "Padplane",   5, 1, isxfld, sxmgmx
 601                 , tmaxfd, stemax, deemax, epsil, stmin);
 602   // Fee + cables
 603   AliMedium(6, "Readout",    1, 0, isxfld, sxmgmx
 604                 , tmaxfd, stemax, deemax, epsil, stmin);
 605   // C frame
 606   AliMedium(7, "C Frame",    6, 0, isxfld, sxmgmx
 607                 , tmaxfd, stemax, deemax, epsil, stmin);
 608   // Mylar foils
 609   AliMedium(8, "Mylar",      7, 0, isxfld, sxmgmx
 610                 , tmaxfd, stemax, deemax, epsil, stmin);
 611   if (fGasMix == 1) {
 612     // Gas-mixture (Xe/CO2)
 613     AliMedium(9, "Gas-mix",   10, 1, isxfld, sxmgmx
 614                   , tmaxfd, stemax, deemax, epsil, stmin);
 615   }
 616   else {
 617     // Gas-mixture (Xe/Isobutane)
 618     AliMedium(9, "Gas-mix",   11, 1, isxfld, sxmgmx
 619                   , tmaxfd, stemax, deemax, epsil, stmin);
 620   }
 621   // Nomex-honeycomb (use carbon for the time being)
 622   AliMedium(10, "Nomex",      6, 0, isxfld, sxmgmx
 623                 , tmaxfd, stemax, deemax, epsil, stmin);
 624   // Kapton foils (use Mylar for the time being)
 625   AliMedium(11, "Kapton",     7, 0, isxfld, sxmgmx
 626                 , tmaxfd, stemax, deemax, epsil, stmin);
 627   // Gas-filling of the radiator
 628   AliMedium(12, "CO2",        8, 0, isxfld, sxmgmx
 629                 , tmaxfd, stemax, deemax, epsil, stmin);
 630   // G10-plates
 631   AliMedium(13, "G10-plates",12, 0, isxfld, sxmgmx
 632                 , tmaxfd, stemax, deemax, epsil, stmin);
 633   // Cooling water
 634   AliMedium(14, "Water",     13, 0, isxfld, sxmgmx
 635                 , tmaxfd, stemax, deemax, epsil, stmin);
 636   // Rohacell (plexiglas) for the radiator
 637   AliMedium(15, "Rohacell",  14, 0, isxfld, sxmgmx
 638                 , tmaxfd, stemax, deemax, epsil, stmin);
 639   // Al layer in MCMs
 640   AliMedium(16, "MCM-Al"  ,   1, 0, isxfld, sxmgmx
 641                 , tmaxfd, stemax, deemax, epsil, stmin);
 642   // Sn layer in MCMs
 643   AliMedium(17, "MCM-Sn"  ,  15, 0, isxfld, sxmgmx
 644                 , tmaxfd, stemax, deemax, epsil, stmin);
 645   // Cu layer in MCMs
 646   AliMedium(18, "MCM-Cu"  ,   5, 0, isxfld, sxmgmx
 647                 , tmaxfd, stemax, deemax, epsil, stmin);
 648   // G10 layer in MCMs
 649   AliMedium(19, "MCM-G10" ,  12, 0, isxfld, sxmgmx
 650                 , tmaxfd, stemax, deemax, epsil, stmin);
 651   // Si in readout chips
 652   AliMedium(20, "Chip-Si" ,  16, 0, isxfld, sxmgmx
 653                 , tmaxfd, stemax, deemax, epsil, stmin);
 654   // Epoxy in readout chips
 655   AliMedium(21, "Chip-Ep" ,  17, 0, isxfld, sxmgmx
 656                 , tmaxfd, stemax, deemax, epsil, stmin);
 657   // PE in connectors
 658   AliMedium(22, "Conn-PE" ,   3, 0, isxfld, sxmgmx
 659                 , tmaxfd, stemax, deemax, epsil, stmin);
 660   // Cu in connectors
 661   AliMedium(23, "Chip-Cu" ,   5, 0, isxfld, sxmgmx
 662                 , tmaxfd, stemax, deemax, epsil, stmin);
 663   // Al of cooling pipes
 664   AliMedium(24, "Cooling" ,   1, 0, isxfld, sxmgmx
 665                 , tmaxfd, stemax, deemax, epsil, stmin);
 666   // Cu in services
 667   AliMedium(25, "Serv-Cu" ,   5, 0, isxfld, sxmgmx
 668                 , tmaxfd, stemax, deemax, epsil, stmin);
 669
 670   // Save the density values for the TRD absorbtion
 671   fFoilDensity = dmy;
 672   if (fGasMix == 1)
 673     fGasDensity = dgm1;
 674   else
 675     fGasDensity = dgm2;
 676
 677 }
 678
 679 //_____________________________________________________________________________
 680 void AliTRD::DrawModule() const
 681 {
 682   //
 683   // Draw a shaded view of the Transition Radiation Detector version 0
 684   //
 685
 686   // Set everything unseen
 687   gMC->Gsatt("*"   ,"SEEN",-1);
 688
 689   // Set ALIC mother transparent
 690   gMC->Gsatt("ALIC","SEEN", 0);
 691
 692   // Set the volumes visible
 693   if (fGeometry->IsVersion() == 0) {
 694     gMC->Gsatt("B071","SEEN", 0);
 695     gMC->Gsatt("B074","SEEN", 0);
 696     gMC->Gsatt("B075","SEEN", 0);
 697     gMC->Gsatt("B077","SEEN", 0);
 698     gMC->Gsatt("BTR1","SEEN", 0);
 699     gMC->Gsatt("BTR2","SEEN", 0);
 700     gMC->Gsatt("BTR3","SEEN", 0);
 701     gMC->Gsatt("UTR1","SEEN", 0);
 702     gMC->Gsatt("UTR2","SEEN", 0);
 703     gMC->Gsatt("UTR3","SEEN", 0);
 704   }
 705   else {
 706     gMC->Gsatt("B071","SEEN", 0);
 707     gMC->Gsatt("B074","SEEN", 0);
 708     gMC->Gsatt("B075","SEEN", 0);
 709     gMC->Gsatt("B077","SEEN", 0);
 710     gMC->Gsatt("BTR1","SEEN", 0);
 711     gMC->Gsatt("BTR2","SEEN", 0);
 712     gMC->Gsatt("BTR3","SEEN", 0);
 713     gMC->Gsatt("UTR1","SEEN", 0);
 714     if (fGeometry->GetPHOShole())
 715       gMC->Gsatt("UTR2","SEEN", 0);
 716     if (fGeometry->GetRICHhole())
 717       gMC->Gsatt("UTR3","SEEN", 0);
 718   }
 719 //   gMC->Gsatt("UCII","SEEN", 0);
 720 //   gMC->Gsatt("UCIM","SEEN", 0);
 721 //   gMC->Gsatt("UCIO","SEEN", 0);
 722 //   gMC->Gsatt("UL02","SEEN", 1);
 723 //   gMC->Gsatt("UL05","SEEN", 1);
 724 //   gMC->Gsatt("UL06","SEEN", 1);
 725
 726   gMC->Gdopt("hide", "on");
 727   gMC->Gdopt("shad", "on");
 728   gMC->Gsatt("*", "fill", 7);
 729   gMC->SetClipBox(".");
 730   gMC->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
 731   gMC->DefaultRange();
 732   gMC->Gdraw("alic", 40, 30, 0, 12, 9.4, .021, .021);
 733   gMC->Gdhead(1111, "Transition Radiation Detector");
 734   gMC->Gdman(18, 4, "MAN");
 735
 736 }
 737
 738 //_____________________________________________________________________________
 739 Int_t AliTRD::DistancetoPrimitive(Int_t , Int_t ) const
 740 {
 741   //
 742   // Distance between the mouse and the TRD detector on the screen
 743   // Dummy routine
 744
 745   return 9999;
 746
 747 }
 748
 749 //_____________________________________________________________________________
 750 void AliTRD::Init()
 751 {
 752   //
 753   // Initialize the TRD detector after the geometry has been created
 754   //
 755
 756   Int_t i;
 757
 758   if (fDebug) {
 759     printf("\n%s: ",ClassName());
 760     for (i = 0; i < 35; i++) printf("*");
 761     printf(" TRD_INIT ");
 762     for (i = 0; i < 35; i++) printf("*");
 763     printf("\n");
 764   }
 765
 766   if (fGeometry->IsVersion() == 1) {
 767     printf("%s: Full geometry version initialized\n",ClassName());
 768     if (fGeometry->GetPHOShole())
 769       printf("%s: Leave space in front of PHOS free\n",ClassName());
 770     if (fGeometry->GetRICHhole())
 771       printf("%s: Leave space in front of RICH free\n",ClassName());
 772   }
 773   else {
 774     printf("%s: Not a valid geometry\n",ClassName());
 775   }
 776
 777   if (fGasMix == 1) {
 778     printf("%s: Gas Mixture: 85%% Xe + 15%% CO2\n",ClassName());
 779   }
 780   else {
 781     printf("%s: Gas Mixture: 97%% Xe + 3%% Isobutane\n",ClassName());
 782   }
 783
 784 }
 785
 786 //_____________________________________________________________________________
 787 void AliTRD::LoadPoints(Int_t )
 788 {
 789   //
 790   // Store x, y, z of all hits in memory.
 791   // Hit originating from TR photons are given a different color
 792   //
 793
 794   //if (!fDrawTR) {
 795   //  AliDetector::LoadPoints(track);
 796   //  return;
 797   //}
 798
 799   if ((fHits == 0) && (fTrackHits == 0)) return;
 800
 801   Int_t nhits;
 802   if (fHitType < 2) {
 803     nhits = fHits->GetEntriesFast();
 804   }
 805   else {
 806     nhits = fTrackHits->GetEntriesFast();
 807   }
 808   if (nhits == 0) return;
 809
 810   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
 811   if (fPoints == 0) fPoints = new TObjArray(tracks);
 812
 813   AliTRDhit *ahit;
 814
 815   Int_t    *ntrkE = new Int_t[tracks];
 816   Int_t    *ntrkT = new Int_t[tracks];
 817   Int_t    *limiE = new Int_t[tracks];
 818   Int_t    *limiT = new Int_t[tracks];
 819   Float_t **coorE = new Float_t*[tracks];
 820   Float_t **coorT = new Float_t*[tracks];
 821   for(Int_t i = 0; i < tracks; i++) {
 822     ntrkE[i] = 0;
 823     ntrkT[i] = 0;
 824     coorE[i] = 0;
 825     coorT[i] = 0;
 826     limiE[i] = 0;
 827     limiT[i] = 0;
 828   }
 829
 830   AliTRDpoints *points = 0;
 831   Float_t      *fp     = 0;
 832   Int_t         trk;
 833   Int_t         chunk  = nhits / 4 + 1;
 834
 835   // Loop over all the hits and store their position
 836   ahit = (AliTRDhit *) FirstHit(-1);
 837   while (ahit) {
 838
 839     // dEdx hits
 840     if (ahit->GetCharge() >= 0) {
 841
 842       trk = ahit->GetTrack();
 843       if (ntrkE[trk] == limiE[trk]) {
 844         // Initialise a new track
 845         fp = new Float_t[3*(limiE[trk]+chunk)];
 846         if (coorE[trk]) {
 847           memcpy(fp,coorE[trk],sizeof(Float_t)*3*limiE[trk]);
 848           delete [] coorE[trk];
 849         }
 850         limiE[trk] += chunk;
 851         coorE[trk]  = fp;
 852       }
 853       else {
 854         fp = coorE[trk];
 855       }
 856       fp[3*ntrkE[trk]  ] = ahit->X();
 857       fp[3*ntrkE[trk]+1] = ahit->Y();
 858       fp[3*ntrkE[trk]+2] = ahit->Z();
 859       ntrkE[trk]++;
 860
 861     }
 862     // TR photon hits
 863     else if ((ahit->GetCharge() < 0) && (fDrawTR)) {
 864
 865       trk = ahit->GetTrack();
 866       if (ntrkT[trk] == limiT[trk]) {
 867         // Initialise a new track
 868         fp = new Float_t[3*(limiT[trk]+chunk)];
 869         if (coorT[trk]) {
 870           memcpy(fp,coorT[trk],sizeof(Float_t)*3*limiT[trk]);
 871           delete [] coorT[trk];
 872         }
 873         limiT[trk] += chunk;
 874         coorT[trk]  = fp;
 875       }
 876       else {
 877         fp = coorT[trk];
 878       }
 879       fp[3*ntrkT[trk]  ] = ahit->X();
 880       fp[3*ntrkT[trk]+1] = ahit->Y();
 881       fp[3*ntrkT[trk]+2] = ahit->Z();
 882       ntrkT[trk]++;
 883
 884     }
 885
 886     ahit = (AliTRDhit *) NextHit();
 887
 888   }
 889
 890   for (trk = 0; trk < tracks; ++trk) {
 891
 892     if (ntrkE[trk] || ntrkT[trk]) {
 893
 894       points = new AliTRDpoints();
 895       points->SetDetector(this);
 896       points->SetParticle(trk);
 897
 898       // Set the dEdx points
 899       if (ntrkE[trk]) {
 900         points->SetMarkerColor(GetMarkerColor());
 901         points->SetMarkerSize(GetMarkerSize());
 902         points->SetPolyMarker(ntrkE[trk],coorE[trk],GetMarkerStyle());
 903         delete [] coorE[trk];
 904         coorE[trk] = 0;
 905       }
 906
 907       // Set the TR photon points
 908       if (ntrkT[trk]) {
 909         points->SetTRpoints(ntrkT[trk],coorT[trk]);
 910         delete [] coorT[trk];
 911         coorT[trk] = 0;
 912       }
 913
 914       fPoints->AddAt(points,trk);
 915
 916     }
 917
 918   }
 919
 920   delete [] coorE;
 921   delete [] coorT;
 922   delete [] ntrkE;
 923   delete [] ntrkT;
 924   delete [] limiE;
 925   delete [] limiT;
 926
 927 }
 928
 929 //_____________________________________________________________________________
 930 void AliTRD::MakeBranch(Option_t* option)
 931 {
 932   //
 933   // Create Tree branches for the TRD digits.
 934   //
 935
 936   Int_t  buffersize = 4000;
 937   Char_t branchname[15];
 938   sprintf(branchname,"%s",GetName());
 939
 940   const char *cD = strstr(option,"D");
 941
 942   AliDetector::MakeBranch(option);
 943
 944   if (fDigits && gAlice->TreeD() && cD) {
 945     MakeBranchInTree(gAlice->TreeD(),branchname,&fDigits,buffersize,0);
 946   }
 947
 948   if (fHitType > 1) {
 949     MakeBranch2(option,0);
 950   }
 951
 952 }
 953
 954 //_____________________________________________________________________________
 955 void AliTRD::ResetDigits()
 956 {
 957   //
 958   // Reset number of digits and the digits array for this detector
 959   //
 960
 961   fNdigits = 0;
 962   if (fDigits) fDigits->Clear();
 963
 964 }
 965
 966 //_____________________________________________________________________________
 967 void AliTRD::SetTreeAddress()
 968 {
 969   //
 970   // Set the branch addresses for the trees.
 971   //
 972
 973   if ( fLoader->TreeH() && (fHits == 0x0))  fHits = new TClonesArray("AliTRDhit",405);
 974   AliDetector::SetTreeAddress();
 975
 976   if (fHitType > 0) {
 977     SetTreeAddress2();
 978   }
 979
 980 }
 981
 982 //_____________________________________________________________________________
 983 void AliTRD::SetGasMix(Int_t imix)
 984 {
 985   //
 986   // Defines the gas mixture (imix=0:  Xe/Isobutane imix=1: Xe/CO2)
 987   //
 988
 989   if ((imix < 0) || (imix > 1)) {
 990     printf("Wrong input value: %d\n",imix);
 991     printf("Use standard setting\n");
 992     fGasMix = 1;
 993     return;
 994   }
 995
 996   fGasMix = imix;
 997
 998 }
 999
1000 //_____________________________________________________________________________
1001 void AliTRD::SetPHOShole()
1002 {
1003   //
1004   // Selects a geometry with a hole in front of the PHOS
1005   //
1006
1007   fGeometry->SetPHOShole();
1008
1009 }
1010
1011 //_____________________________________________________________________________
1012 void AliTRD::SetRICHhole()
1013 {
1014   //
1015   // Selects a geometry with a hole in front of the RICH
1016   //
1017
1018   fGeometry->SetRICHhole();
1019
1020 }
1021
1022 //_____________________________________________________________________________
1023 AliTRD &AliTRD::operator=(const AliTRD &trd)
1024 {
1025   //
1026   // Assignment operator
1027   //
1028
1029   if (this != &trd) ((AliTRD &) trd).Copy(*this);
1030   return *this;
1031
1032 }
1033
1034 //_____________________________________________________________________________
1035 void AliTRD::FinishPrimary()
1036 {
1037   //
1038   // Store the hits in the containers after all primaries are finished
1039   //
1040
1041   if (fTrackHits) {
1042     fTrackHits->FlushHitStack();
1043   }
1044
1045 }
1046
1047 //_____________________________________________________________________________
1048 void AliTRD::RemapTrackHitIDs(Int_t *map)
1049 {
1050   //
1051   // Remap the track IDs
1052   //
1053
1054   if (!fTrackHits) {
1055     return;
1056   }
1057
1058   if (fTrackHits) {
1059     TClonesArray *arr = fTrackHits->GetArray();;
1060     for (Int_t i = 0; i < arr->GetEntriesFast(); i++){
1061       AliTrackHitsParamV2 *info = (AliTrackHitsParamV2 *) (arr->At(i));
1062       info->fTrackID = map[info->fTrackID];
1063     }
1064   }
1065
1066 }
1067
1068 //_____________________________________________________________________________
1069 void AliTRD::ResetHits()
1070 {
1071   //
1072   // Reset the hits
1073   //
1074
1075   AliDetector::ResetHits();
1076   if (fTrackHits) {
1077     fTrackHits->Clear();
1078   }
1079
1080 }
1081
1082 //_____________________________________________________________________________
1083 AliHit* AliTRD::FirstHit(Int_t track)
1084 {
1085   //
1086   // Return the first hit of a track
1087   //
1088
1089   if (fHitType > 1) {
1090     return FirstHit2(track);
1091   }
1092
1093   return AliDetector::FirstHit(track);
1094
1095 }
1096
1097 //_____________________________________________________________________________
1098 AliHit* AliTRD::NextHit()
1099 {
1100   //
1101   // Returns the next hit of a track
1102   //
1103
1104   if (fHitType > 1) {
1105     return NextHit2();
1106   }
1107
1108   return AliDetector::NextHit();
1109
1110 }
1111
1112 //_____________________________________________________________________________
1113 AliHit* AliTRD::FirstHit2(Int_t track)
1114 {
1115   //
1116   // Initializes the hit iterator.
1117   // Returns the address of the first hit of a track.
1118   // If <track> >= 0 the track is read from disk,
1119   // while if <track> < 0 the first hit of the current
1120   // track is returned.
1121   //
1122
1123   if (track >= 0) {
1124     gAlice->ResetHits();
1125     TreeH()->GetEvent(track);
1126   }
1127
1128   if (fTrackHits) {
1129     fTrackHits->First();
1130     return (AliHit*) fTrackHits->GetHit();
1131   }
1132   else {
1133     return 0;
1134   }
1135
1136 }
1137
1138 //_____________________________________________________________________________
1139 AliHit* AliTRD::NextHit2()
1140 {
1141   //
1142   // Returns the next hit of the current track
1143   //
1144
1145   if (fTrackHits) {
1146     fTrackHits->Next();
1147     return (AliHit *) fTrackHits->GetHit();
1148   }
1149   else {
1150     return 0;
1151   }
1152
1153 }
1154
1155 //_____________________________________________________________________________
1156 void AliTRD::MakeBranch2(Option_t *option, const char* )
1157 {
1158   //
1159   // Create a new branch in the current Root tree.
1160   // The branch of fHits is automatically split.
1161   //
1162
1163   if (fHitType < 2) {
1164     return;
1165   }
1166
1167   char branchname[10];
1168   sprintf(branchname,"%s2",GetName());
1169
1170   // Get the pointer to the header
1171   const char *cH = strstr(option,"H");
1172
1173   if (!fTrackHits) {
1174     fTrackHits = new AliTRDtrackHits();
1175    }
1176
1177
1178   if (fTrackHits && TreeH() && cH)
1179    {
1180     TreeH()->Branch(branchname,"AliTRDtrackHits",&fTrackHits,fBufferSize,99);
1181     Info("MakeBranch2","Making Branch %s for trackhits",branchname);
1182    }
1183 }
1184
1185 //_____________________________________________________________________________
1186 void AliTRD::SetTreeAddress2()
1187 {
1188   //
1189   // Set the branch address for the trackHits tree
1190   //
1191   if (GetDebug())  Info("SetTreeAddress2","");
1192
1193   TBranch *branch;
1194   char branchname[20];
1195   sprintf(branchname,"%s2",GetName());
1196
1197   // Branch address for hit tree
1198   TTree *treeH = TreeH();
1199   if ((treeH) && (fHitType > 0)) {
1200     branch = treeH->GetBranch(branchname);
1201     if (branch)
1202      {
1203        branch->SetAddress(&fTrackHits);
1204        if (GetDebug())
1205          Info("SetTreeAddress2","Success.");
1206      }
1207     else
1208      {
1209        if (GetDebug())
1210          Info("SetTreeAddress2","Can NOT get the branch %s",branchname);
1211      }
1212   }
1213
1214 }
1215
1216 //_____________________________________________________________________________
1217 void AliTRD::AddHit2(Int_t track, Int_t det, Float_t *hits, Int_t q
1218                    , Bool_t inDrift)
1219 {
1220   //
1221   // Add a hit to the list
1222   //
1223
1224   Int_t rtrack;
1225
1226   if (fIshunt) {
1227     Int_t primary = gAlice->GetMCApp()->GetPrimary(track);
1228     gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
1229     rtrack = primary;
1230   }
1231   else {
1232     rtrack = track;
1233     gAlice->GetMCApp()->FlagTrack(track);
1234   }
1235
1236   if ((fTrackHits) && (fHitType > 0)) {
1237     fTrackHits->AddHitTRD(det,rtrack,hits[0],hits[1],hits[2],q,inDrift);
1238   }
1239
1240 }
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577