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