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