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