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