Move Gstpar calls to Init()
[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 <TClonesArray.h>
30 #include <TFile.h>
31 #include <TGeometry.h>
32 #include <TLorentzVector.h>
33 #include <TMath.h>
34 #include <TNode.h>
35 #include <TPGON.h> 
36 #include <TParticle.h>
37 #include <TROOT.h>
38 #include <TTree.h>
39 #include <TVirtualMC.h>
40  
41 #include "AliConst.h"
42 #include "AliDigit.h"
43 #include "AliLoader.h"
44 #include "AliLog.h"
45 #include "AliMC.h"
46 #include "AliMagF.h"
47 #include "AliRun.h"
48 #include "AliTrackReference.h"
49 #include "AliRawReader.h"
50
51 #include "AliTRD.h"
52 #include "AliTRDdigit.h"
53 #include "AliTRDdigitizer.h"
54 #include "AliTRDdigitsManager.h"
55 #include "AliTRDgeometry.h"
56 #include "AliTRDhit.h"
57 #include "AliTRDpoints.h"
58 #include "AliTRDrawData.h"
59 #include "AliTRDSimParam.h"
60 #include "AliTRDRecParam.h"
61 #include "AliTRDCommonParam.h"
62 #include "AliTRDcalibDB.h"
63
64 ClassImp(AliTRD)
65  
66 //_____________________________________________________________________________
67 AliTRD::AliTRD()
68   :AliDetector()
69   ,fGeometry(0)
70   ,fGasDensity(0)
71   ,fFoilDensity(0)
72   ,fDrawTR(0)
73   ,fDisplayType(0)
74 {
75   //
76   // Default constructor
77   //
78  
79 }
80  
81 //_____________________________________________________________________________
82 AliTRD::AliTRD(const char *name, const char *title)
83   :AliDetector(name,title)
84   ,fGeometry(0)
85   ,fGasDensity(0)
86   ,fFoilDensity(0)
87   ,fDrawTR(0)
88   ,fDisplayType(0)
89 {
90   //
91   // Standard constructor for the TRD
92   //
93
94   // Check that FRAME is there otherwise we have no place where to put TRD
95   AliModule *frame = gAlice->GetModule("FRAME");
96   if (!frame) {
97     AliError("TRD needs FRAME to be present\n");
98     exit(1);
99   } 
100
101   // Define the TRD geometry
102   if ((frame->IsVersion() == 0) ||
103       (frame->IsVersion() == 1)) {
104     fGeometry = new AliTRDgeometry();
105   }
106   else {
107     AliError("Could not find valid FRAME version\n");
108     exit(1);
109   }
110
111   // Save the geometry
112   TDirectory *saveDir = gDirectory;
113   gAlice->GetRunLoader()->CdGAFile();
114   fGeometry->Write("TRDgeometry");
115   saveDir->cd();
116
117   // Allocate the hit array
118   fHits = new TClonesArray("AliTRDhit",405);
119   gAlice->GetMCApp()->AddHitList(fHits);
120
121   //PH SetMarkerColor(kWhite);   
122
123 }
124
125 //_____________________________________________________________________________
126 AliTRD::~AliTRD()
127 {
128   //
129   // TRD destructor
130   //
131
132   if (fGeometry) {
133     delete fGeometry;
134     fGeometry = 0;
135   }
136
137   if (fHits) {
138     delete fHits;
139     fHits     = 0;
140   }
141
142 }
143
144 //_____________________________________________________________________________
145 void AliTRD::Hits2Digits()
146 {
147   //
148   // Create digits
149   //
150
151   AliTRDdigitizer digitizer("TRDdigitizer","TRD digitizer class");
152   AliLog::SetClassDebugLevel("TRDdigitizer",AliDebugLevel());
153   
154   // Initialization
155   digitizer.InitDetector();
156     
157   if (!fLoader->TreeH()) {
158     fLoader->LoadHits("read");
159   }
160   fLoader->LoadDigits("recreate");
161
162   AliRunLoader *runLoader = fLoader->GetRunLoader(); 
163
164   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
165     runLoader->GetEvent(iEvent);
166     digitizer.Open(runLoader,iEvent);
167     digitizer.MakeDigits();
168     digitizer.WriteDigits();
169   }
170
171   fLoader->UnloadHits();
172   fLoader->UnloadDigits();
173
174 }
175
176 //_____________________________________________________________________________
177 void AliTRD::Hits2SDigits()
178 {
179   //
180   // Create summable digits
181   //
182
183   AliTRDdigitizer digitizer("TRDdigitizer","TRD digitizer class");
184   // For the summable digits
185   digitizer.SetSDigits(kTRUE);
186   AliLog::SetClassDebugLevel("TRDdigitizer",AliDebugLevel());
187
188   // Initialization
189   digitizer.InitDetector();
190     
191   if (!fLoader->TreeH()) {
192     fLoader->LoadHits("read");
193   }
194   fLoader->LoadSDigits("recreate");
195
196   AliRunLoader *runLoader = fLoader->GetRunLoader(); 
197
198   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
199     runLoader->GetEvent(iEvent);
200     digitizer.Open(runLoader,iEvent);
201     digitizer.MakeDigits();
202     digitizer.WriteDigits();
203   }
204
205   fLoader->UnloadHits();
206   fLoader->UnloadSDigits();
207   
208 }
209
210 //_____________________________________________________________________________
211 AliDigitizer *AliTRD::CreateDigitizer(AliRunDigitizer *manager) const
212 {
213   //
214   // Creates a new digitizer object
215   //
216
217   return new AliTRDdigitizer(manager);
218
219 }
220
221 //_____________________________________________________________________________
222 void AliTRD::SDigits2Digits()
223 {
224   //
225   // Create final digits from summable digits
226   //
227
228   // Create the TRD digitizer
229   AliTRDdigitizer digitizer("TRDdigitizer","TRD digitizer class");  
230   AliLog::SetClassDebugLevel("TRDdigitizer",AliDebugLevel());
231
232   // Set the parameter
233   digitizer.SetEvent(gAlice->GetEvNumber());
234
235   // Initialization
236   digitizer.InitDetector();
237
238   // Read the s-digits via digits manager
239   AliTRDdigitsManager sdigitsManager;
240  
241   AliLog::SetClassDebugLevel("TRDdigitisManager",AliDebugLevel());
242   sdigitsManager.SetSDigits(kTRUE);
243   sdigitsManager.CreateArrays();
244   
245   if (!fLoader->TreeS()) { 
246     if (fLoader->LoadSDigits("read")) {
247       return;
248     }
249   }
250   if (!fLoader->TreeS()) {
251     AliError(Form("Error while reading SDigits for event %d",gAlice->GetEvNumber()));
252     return;
253   }
254   
255   sdigitsManager.ReadDigits(fLoader->TreeS());
256
257   // Add the s-digits to the input list 
258   digitizer.AddSDigitsManager(&sdigitsManager);
259
260   // Convert the s-digits to normal digits
261   digitizer.SDigits2Digits();
262
263   // Store the digits
264   if (!fLoader->TreeD()) {
265     fLoader->MakeTree("D");
266   }
267   if (digitizer.MakeBranch(fLoader->TreeD())){
268     digitizer.WriteDigits();
269   }
270
271 }
272
273 //_____________________________________________________________________________
274 void AliTRD::Digits2Raw() 
275 {
276   //
277   // Convert digits of the current event to raw data
278   //
279
280   fLoader->LoadDigits();
281   TTree *digits = fLoader->TreeD();
282   if (!digits) {
283     AliError("No digits tree");
284     return;
285   }
286
287   AliTRDrawData rawWriter;
288   if (!rawWriter.Digits2Raw(digits)) {
289     AliError("The raw writer could not load the digits tree");
290   }
291
292   fLoader->UnloadDigits();
293
294 }
295
296 //_____________________________________________________________________________
297 void AliTRD::AddHit(Int_t track, Int_t det, Float_t *hits, Int_t q
298                   , Float_t time, Bool_t inDrift)
299 {
300   //
301   // Add a hit for the TRD
302   // 
303
304   TClonesArray &lhits = *fHits;
305   AliTRDhit *hit = new(lhits[fNhits++]) AliTRDhit(fIshunt
306                                                  ,track
307                                                  ,det
308                                                  ,hits
309                                                  ,q
310                                                  ,time);
311
312   if (inDrift) {
313     hit->SetDrift();
314   }
315   else {
316     hit->SetAmplification();
317   }
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;
333   TNode *top;
334   TPGON *pgon;
335
336   // The dimensions of the TRD super module
337   const Float_t kRmin  = 291.0;
338   const Float_t kRmax  = 369.0;
339   const Float_t kZmax1 = 378.35;
340   const Float_t kZmax2 = 302.0;
341
342   Float_t rmin;
343   Float_t rmax;
344   Float_t zmax1;
345   Float_t zmax2;
346
347   Int_t   iPlan;
348  
349   const Int_t kColorTRD = 46;
350   
351   // Find the top node alice
352   top = gAlice->GetGeometry()->GetNode("alice");
353   
354   if      (fDisplayType == 0) {
355
356     pgon = new TPGON("S_TRD","TRD","void",0,360,AliTRDgeometry::Nsect(),4);
357     rmin = kRmin;
358     rmax = kRmax;
359     pgon->DefineSection(0,-kZmax1,rmax,rmax);
360     pgon->DefineSection(1,-kZmax2,rmin,rmax);
361     pgon->DefineSection(2, kZmax2,rmin,rmax);
362     pgon->DefineSection(3, kZmax1,rmax,rmax);
363     top->cd();
364     node = new TNode("TRD","TRD","S_TRD",0,0,0,"");
365     node->SetLineColor(kColorTRD);
366     fNodes->Add(node);
367
368   }
369   else if (fDisplayType == 1) {
370
371     Char_t name[7];
372
373     Float_t slope = (kZmax1 - kZmax2) / (kRmax  - kRmin);
374
375     rmin  = kRmin + AliTRDgeometry::CraHght();
376     rmax  = rmin  + AliTRDgeometry::CdrHght();
377
378     Float_t thickness = rmin - kRmin;
379     zmax2 = kZmax2 + slope * thickness;
380     zmax1 = zmax2 + slope * AliTRDgeometry::DrThick();
381
382     for (iPlan = 0; iPlan < AliTRDgeometry::Nplan(); iPlan++) {
383
384       sprintf(name,"S_TR1%d",iPlan);
385       pgon  = new TPGON(name,"TRD","void",0,360,AliTRDgeometry::Nsect(),4);
386       pgon->DefineSection(0,-zmax1,rmax,rmax);
387       pgon->DefineSection(1,-zmax2,rmin,rmax);
388       pgon->DefineSection(2, zmax2,rmin,rmax);
389       pgon->DefineSection(3, zmax1,rmax,rmax);
390       top->cd();
391       node = new TNode("TRD","TRD",name,0,0,0,"");
392       node->SetLineColor(kColorTRD);
393       fNodes->Add(node);
394
395       Float_t height = AliTRDgeometry::Cheight() + AliTRDgeometry::Cspace(); 
396       rmin  = rmin  + height;
397       rmax  = rmax  + height;
398       zmax1 = zmax1 + slope * height;
399       zmax2 = zmax2 + slope * height;
400
401     }
402
403     thickness += AliTRDgeometry::DrThick();
404     rmin       = kRmin  + thickness;
405     rmax       = rmin   + AliTRDgeometry::AmThick();
406     zmax2      = kZmax2 + slope * thickness;
407     zmax1      = zmax2  + slope * AliTRDgeometry::AmThick();
408
409     for (iPlan = 0; iPlan < AliTRDgeometry::Nplan(); iPlan++) {
410
411       sprintf(name,"S_TR2%d",iPlan);
412       pgon  = new TPGON(name,"TRD","void",0,360,AliTRDgeometry::Nsect(),4);
413       pgon->DefineSection(0,-zmax1,rmax,rmax);
414       pgon->DefineSection(1,-zmax2,rmin,rmax);
415       pgon->DefineSection(2, zmax2,rmin,rmax);
416       pgon->DefineSection(3, zmax1,rmax,rmax);
417       top->cd();
418       node = new TNode("TRD","TRD",name,0,0,0,"");
419       node->SetLineColor(kColorTRD);
420       fNodes->Add(node);
421
422       Float_t height = AliTRDgeometry::Cheight() + AliTRDgeometry::Cspace(); 
423       rmin  = rmin  + height;
424       rmax  = rmax  + height;
425       zmax1 = zmax1 + slope * height;
426       zmax2 = zmax2 + slope * height;
427
428     }
429
430   }
431
432 }
433  
434 //_____________________________________________________________________________
435 void AliTRD::CreateGeometry()
436 {
437   //
438   // Creates the volumes for the TRD chambers
439   //
440
441   // Check that FRAME is there otherwise we have no place where to put the TRD
442   AliModule *frame = gAlice->GetModule("FRAME");
443   if (!frame) {
444     AliFatal("The TRD needs the FRAME to be defined first");
445   }
446
447   fGeometry->CreateGeometry(fIdtmed->GetArray() - 1299);
448
449 }
450  
451 //_____________________________________________________________________________
452 void AliTRD::CreateMaterials()
453 {
454   //
455   // Create the materials for the TRD
456   //
457
458   Int_t   isxfld = gAlice->Field()->Integ();
459   Float_t sxmgmx = gAlice->Field()->Max();
460   
461   // For polyethilene (CH2) 
462   Float_t ape[2]    = { 12.011 ,  1.0079 };
463   Float_t zpe[2]    = {  6.0   ,  1.0    };
464   Float_t wpe[2]    = {  1.0   ,  2.0    };
465   Float_t dpe       = 0.95;
466
467   // For CO2 
468   Float_t aco[2]    = { 12.011 , 15.9994 };
469   Float_t zco[2]    = {  6.0   ,  8.0    };
470   Float_t wco[2]    = {  1.0   ,  2.0    };
471   Float_t dco       = 0.00186;
472
473   // For water
474   Float_t awa[2]    = {  1.0079, 15.9994 };
475   Float_t zwa[2]    = {  1.0   ,  8.0    };
476   Float_t wwa[2]    = {  2.0   ,  1.0    };
477   Float_t dwa       = 1.0;
478
479   // For isobutane (C4H10)
480   Float_t ais[2]    = { 12.011 ,  1.0079 };
481   Float_t zis[2]    = {  6.0   ,  1.0    };
482   Float_t wis[2]    = {  4.0   , 10.0    };
483   Float_t dis       = 0.00267;
484
485   // For plexiglas (C5H8O2)
486   Float_t apg[3]    = { 12.011 ,  1.0079, 15.9994 };
487   Float_t zpg[3]    = {  6.0   ,  1.0   ,  8.0    };
488   Float_t wpg[3]    = {  5.0   ,  8.0   ,  2.0    };
489   Float_t dpg       = 1.18; 
490   
491   // For epoxy (C18H19O3)
492   Float_t aEpoxy[3] = { 15.9994,  1.0079, 12.011  }; 
493   Float_t zEpoxy[3] = {  8.0   ,  1.0   ,  6.0    }; 
494   Float_t wEpoxy[3] = {  3.0   , 19.0   , 18.0    }; 
495   Float_t dEpoxy    = 1.8 ; 
496
497   // For Araldite, low density epoxy (C18H19O3)
498   Float_t aAral[3]  = { 15.9994,  1.0079, 12.011  }; 
499   Float_t zAral[3]  = {  8.0   ,  1.0   ,  6.0    }; 
500   Float_t wAral[3]  = {  3.0   , 19.0   , 18.0    }; 
501   Float_t dAral     = 1.05; 
502
503   // For air  
504   Float_t aAir[4]   = { 12.011   , 14.0     , 15.9994  , 36.0      };
505   Float_t zAir[4]   = {  6.0     ,  7.0     ,  8.0     , 18.0      };
506   Float_t wAir[4]   = {  0.000124,  0.755267,  0.231781,  0.012827 };
507   Float_t dAir      = 1.20479e-03;
508
509   // For G10
510   Float_t aG10[4]   = {  1.0079  , 12.011   , 15.9994  , 28.086    };
511   Float_t zG10[4]   = {  1.0     ,  6.0     ,  8.0     , 14.0      };
512   Float_t wG10[4]   = {  0.15201 ,  0.10641 ,  0.49444 ,  0.24714  };
513   Float_t dG10      = 1.7;
514
515   // For Xe/CO2-gas-mixture 
516   Float_t aXeCO2[3] = { 131.29   ,  12.0107 ,  15.9994  };
517   Float_t zXeCO2[3] = {  54.0    ,   6.0    ,   8.0     };
518   Float_t wXeCO2[3] = {   8.5    ,   1.5    ,   3.0     }; 
519   // Xe-content of the Xe/CO2-mixture (85% / 15%) 
520   Float_t fxc       = 0.85;
521   Float_t dxe       = 0.00549;
522   Float_t dgm       = fxc * dxe + (1.0 - fxc) * dco;
523   
524   // General tracking parameter
525   Float_t tmaxfd    = -10.0;
526   Float_t stemax    = -1.0e10;
527   Float_t deemax    = -0.1;
528   Float_t epsil     =  1.0e-4;
529   Float_t stmin     = -0.001;
530   
531   //////////////////////////////////////////////////////////////////////////
532   //     Define Materials 
533   //////////////////////////////////////////////////////////////////////////
534
535   AliMaterial( 1, "Al"   ,  26.98, 13.0, 2.7     ,     8.9 ,    37.2);
536   AliMaterial( 4, "Xe"   , 131.29, 54.0, dxe     ,  1546.16,     0.0);
537   AliMaterial( 5, "Cu"   ,  63.54, 29.0, 8.96    ,     1.43,    14.8);
538   AliMaterial( 6, "C"    ,  12.01,  6.0, 2.265   ,    18.8 ,    74.4);
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(18, "Fe"   ,  55.85, 26.0, 7.87    ,     1.76,    14.8);
542
543   // Mixtures 
544   AliMixture(2, "Air"         , aAir,   zAir,   dAir,    4, wAir  );
545   AliMixture(3, "Polyethilene", ape,    zpe,    dpe,    -2, wpe   );
546   AliMixture(8, "CO2",          aco,    zco,    dco,    -2, wco   );
547   AliMixture(9, "Isobutane",    ais,    zis,    dis,    -2, wis   );
548   AliMixture(10,"Gas mixture",  aXeCO2, zXeCO2, dgm,    -3, wXeCO2);
549   AliMixture(12,"G10",          aG10,   zG10,   dG10,    4, wG10  );
550   AliMixture(13,"Water",        awa,    zwa,    dwa,    -2, wwa   );
551   AliMixture(14,"Plexiglas",    apg,    zpg,    dpg,    -3, wpg   );
552   AliMixture(17,"Epoxy",        aEpoxy, zEpoxy, dEpoxy, -3, wEpoxy);
553   AliMixture(19,"Araldite",     aAral,  zAral,  dAral,  -3, wAral );
554
555   //////////////////////////////////////////////////////////////////////////
556   //     Tracking Media Parameters 
557   //////////////////////////////////////////////////////////////////////////
558
559   // Al Frame 
560   AliMedium( 1,"Al Frame"   , 1,0,isxfld,sxmgmx
561               ,tmaxfd,stemax,deemax,epsil,stmin);
562   // Air 
563   AliMedium( 2,"Air"        , 2,0,isxfld,sxmgmx
564               ,tmaxfd,stemax,deemax,epsil,stmin);
565   // Wires
566   AliMedium( 3,"Wires"      , 5,0,isxfld,sxmgmx
567               ,tmaxfd,stemax,deemax,epsil,stmin);
568   // All other ROB materials (caps, etc.)
569   AliMedium( 4,"ROB Other"  , 5,0,isxfld,sxmgmx
570               ,tmaxfd,stemax,deemax,epsil,stmin);
571   // Cu pads 
572   AliMedium( 5,"Padplane"   , 5,1,isxfld,sxmgmx
573               ,tmaxfd,stemax,deemax,epsil,stmin);
574   // Fee + cables 
575   AliMedium( 6,"Readout"    , 5,0,isxfld,sxmgmx
576               ,tmaxfd,stemax,deemax,epsil,stmin);
577   // C frame 
578   AliMedium( 7,"C Frame"    , 6,0,isxfld,sxmgmx
579               ,tmaxfd,stemax,deemax,epsil,stmin);
580   // INOX of cooling bus bars
581   AliMedium( 8,"Cooling bus",18,0,isxfld,sxmgmx
582               ,tmaxfd,stemax,deemax,epsil,stmin);
583   // Gas-mixture (Xe/CO2) 
584   AliMedium( 9,"Gas-mix"    ,10,1,isxfld,sxmgmx
585               ,tmaxfd,stemax,deemax,epsil,stmin);
586   // Nomex-honeycomb
587   AliMedium(10,"Nomex"      ,12,0,isxfld,sxmgmx
588               ,tmaxfd,stemax,deemax,epsil,stmin);
589   // Araldite glue
590   AliMedium(11,"Glue"       ,19,0,isxfld,sxmgmx
591               ,tmaxfd,stemax,deemax,epsil,stmin);
592   // G10-plates
593   AliMedium(13,"G10-plates" ,12,0,isxfld,sxmgmx
594               ,tmaxfd,stemax,deemax,epsil,stmin);
595   // Cooling water
596   AliMedium(14,"Water"      ,13,0,isxfld,sxmgmx
597               ,tmaxfd,stemax,deemax,epsil,stmin);
598   // Rohacell (plexiglas) for the radiator
599   AliMedium(15,"Rohacell"   ,14,0,isxfld,sxmgmx
600               ,tmaxfd,stemax,deemax,epsil,stmin);
601   // Al layer in MCMs
602   AliMedium(16,"MCM-Al"     , 1,0,isxfld,sxmgmx
603               ,tmaxfd,stemax,deemax,epsil,stmin);
604   // Sn layer in MCMs
605   AliMedium(17,"MCM-Sn"     ,15,0,isxfld,sxmgmx
606               ,tmaxfd,stemax,deemax,epsil,stmin);
607   // Cu layer in MCMs
608   AliMedium(18,"MCM-Cu"     , 5,0,isxfld,sxmgmx
609               ,tmaxfd,stemax,deemax,epsil,stmin);
610   // G10 layer in MCMs
611   AliMedium(19,"MCM-G10"    ,12,0,isxfld,sxmgmx
612               ,tmaxfd,stemax,deemax,epsil,stmin);
613   // Si in readout chips
614   AliMedium(20,"Chip-Si"    ,16,0,isxfld,sxmgmx
615               ,tmaxfd,stemax,deemax,epsil,stmin);
616   // Epoxy in readout chips
617   AliMedium(21,"Chip-Ep"    ,17,0,isxfld,sxmgmx
618               ,tmaxfd,stemax,deemax,epsil,stmin);
619   // PE in connectors
620   AliMedium(22,"Conn-PE"    , 3,0,isxfld,sxmgmx
621               ,tmaxfd,stemax,deemax,epsil,stmin);
622   // Cu in connectors
623   AliMedium(23,"Chip-Cu"    , 5,0,isxfld,sxmgmx
624               ,tmaxfd,stemax,deemax,epsil,stmin);
625   // Al of cooling pipes
626   AliMedium(24,"Cooling"    , 1,0,isxfld,sxmgmx
627               ,tmaxfd,stemax,deemax,epsil,stmin);
628   // Cu in services
629   AliMedium(25,"Serv-Cu"    , 5,0,isxfld,sxmgmx
630               ,tmaxfd,stemax,deemax,epsil,stmin);
631
632   // Save the density values for the TRD absorbtion
633   Float_t dmy  = 1.39;
634   fFoilDensity = dmy;
635   fGasDensity  = dgm;
636
637 }
638
639 //_____________________________________________________________________________
640 void AliTRD::DrawModule() const
641 {
642   //
643   // Draw a shaded view of the Transition Radiation Detector version 0
644   //
645
646   // Set everything unseen
647   gMC->Gsatt("*"   ,"SEEN",-1);
648   
649   // Set ALIC mother transparent
650   gMC->Gsatt("ALIC","SEEN", 0);
651   
652   // Set the volumes visible
653   if (fGeometry->IsVersion() == 0) {
654     gMC->Gsatt("B071","SEEN", 0);
655     gMC->Gsatt("B074","SEEN", 0);
656     gMC->Gsatt("B075","SEEN", 0);
657     gMC->Gsatt("B077","SEEN", 0);
658     gMC->Gsatt("BTR1","SEEN", 0);
659     gMC->Gsatt("BTR2","SEEN", 0);
660     gMC->Gsatt("BTR3","SEEN", 0);
661     gMC->Gsatt("UTR1","SEEN", 0);
662     gMC->Gsatt("UTR2","SEEN", 0);
663     gMC->Gsatt("UTR3","SEEN", 0);
664   }
665   else {
666     gMC->Gsatt("B071","SEEN", 0);
667     gMC->Gsatt("B074","SEEN", 0);
668     gMC->Gsatt("B075","SEEN", 0);
669     gMC->Gsatt("B077","SEEN", 0);
670     gMC->Gsatt("BTR1","SEEN", 0);
671     gMC->Gsatt("BTR2","SEEN", 0);
672     gMC->Gsatt("BTR3","SEEN", 0);
673     gMC->Gsatt("UTR1","SEEN", 0);
674   }
675   
676   gMC->Gdopt("hide", "on");
677   gMC->Gdopt("shad", "on");
678   gMC->Gsatt("*", "fill", 7);
679   gMC->SetClipBox(".");
680   gMC->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
681   gMC->DefaultRange();
682   gMC->Gdraw("alic", 40, 30, 0, 12, 9.4, .021, .021);
683   gMC->Gdhead(1111, "Transition Radiation Detector");
684   gMC->Gdman(18, 4, "MAN");
685
686 }
687
688 //_____________________________________________________________________________
689 Int_t AliTRD::DistancetoPrimitive(Int_t , Int_t )
690 {
691   //
692   // Distance between the mouse and the TRD detector on the screen
693   // Dummy routine
694   //
695   
696   return 9999;
697
698 }
699  
700 //_____________________________________________________________________________
701 void AliTRD::Init()
702 {
703   //
704   // Initialize the TRD detector after the geometry has been created
705   //
706
707   AliDebug(1,"++++++++++++++++++++++++++++++++++++++++++++++");
708
709   if (fGeometry->IsVersion() != 1) {
710     AliError("Not a valid geometry");
711   }
712
713   // Special tracking options for charged particles for XeCO2
714   gMC->Gstpar((* fIdtmed)[9],"DRAY"    , 1.0);
715   gMC->Gstpar((* fIdtmed)[9],"STRA"    , 1.0); 
716   gMC->Gstpar((* fIdtmed)[9],"LOSS"    ,13.0);      // Specific energy loss
717   gMC->Gstpar((* fIdtmed)[9],"PRIMIO_E",23.53);     // 1st ionisation potential
718   gMC->Gstpar((* fIdtmed)[9],"PRIMIO_N",19.344431); // Number of primaries
719
720 }
721
722 //_____________________________________________________________________________
723 void AliTRD::LoadPoints(Int_t )
724 {
725   //
726   // Store x, y, z of all hits in memory.
727   // Hit originating from TR photons are given a different color
728   //
729
730   if (fHits == 0) {
731     return;
732   }
733
734   Int_t nhits  = fHits->GetEntriesFast();
735   if (nhits == 0) {
736     return;
737   }
738
739   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
740   if (fPoints == 0) {
741     fPoints = new TObjArray(tracks);
742   }
743
744   AliTRDhit *ahit;
745   
746   Int_t    *ntrkE = new Int_t[tracks];
747   Int_t    *ntrkT = new Int_t[tracks];
748   Int_t    *limiE = new Int_t[tracks];
749   Int_t    *limiT = new Int_t[tracks];
750   Float_t **coorE = new Float_t*[tracks];
751   Float_t **coorT = new Float_t*[tracks];
752   for(Int_t i = 0; i < tracks; i++) {
753     ntrkE[i] = 0;
754     ntrkT[i] = 0;
755     coorE[i] = 0;
756     coorT[i] = 0;
757     limiE[i] = 0;
758     limiT[i] = 0;
759   }
760   
761   AliTRDpoints *points = 0;
762   Float_t      *fp     = 0;
763   Int_t         trk;
764   Int_t         chunk  = nhits / 4 + 1;
765
766   // Loop over all the hits and store their position
767   ahit = (AliTRDhit *) FirstHit(-1);
768   while (ahit) {
769
770     // dEdx hits
771     if (ahit->GetCharge() >= 0) {
772
773       trk = ahit->GetTrack();
774       if (ntrkE[trk] == limiE[trk]) {
775         // Initialise a new track
776         fp = new Float_t[3*(limiE[trk]+chunk)];
777         if (coorE[trk]) {
778           memcpy(fp,coorE[trk],sizeof(Float_t)*3*limiE[trk]);
779           delete [] coorE[trk];
780         }
781         limiE[trk] += chunk;
782         coorE[trk]  = fp;
783       } 
784       else {
785         fp = coorE[trk];
786       }
787       fp[3*ntrkE[trk]  ] = ahit->X();
788       fp[3*ntrkE[trk]+1] = ahit->Y();
789       fp[3*ntrkE[trk]+2] = ahit->Z();
790       ntrkE[trk]++;
791
792     }
793     // TR photon hits
794     else if ((ahit->GetCharge() < 0) && 
795              (fDrawTR)) {
796
797       trk = ahit->GetTrack();
798       if (ntrkT[trk] == limiT[trk]) {
799         // Initialise a new track
800         fp = new Float_t[3*(limiT[trk]+chunk)];
801         if (coorT[trk]) {
802           memcpy(fp,coorT[trk],sizeof(Float_t)*3*limiT[trk]);
803           delete [] coorT[trk];
804         }
805         limiT[trk] += chunk;
806         coorT[trk]  = fp;
807       } 
808       else {
809         fp = coorT[trk];
810       }
811       fp[3*ntrkT[trk]  ] = ahit->X();
812       fp[3*ntrkT[trk]+1] = ahit->Y();
813       fp[3*ntrkT[trk]+2] = ahit->Z();
814       ntrkT[trk]++;
815
816     }
817
818     ahit = (AliTRDhit *) NextHit();
819
820   }
821
822   for (trk = 0; trk < tracks; ++trk) {
823
824     if (ntrkE[trk] || ntrkT[trk]) {
825
826       points = new AliTRDpoints();
827       points->SetDetector(this);
828       points->SetParticle(trk);
829
830       // Set the dEdx points
831       if (ntrkE[trk]) {
832         points->SetMarkerColor(kWhite); //PH This is the default color in TRD
833         points->SetMarkerSize(1); //PH Default size=1
834         points->SetPolyMarker(ntrkE[trk],coorE[trk],1); //PH Default style=1
835         delete [] coorE[trk];
836         coorE[trk] = 0;
837       }
838
839       // Set the TR photon points
840       if (ntrkT[trk]) {
841         points->SetTRpoints(ntrkT[trk],coorT[trk]);
842         delete [] coorT[trk];
843         coorT[trk] = 0;
844       }
845
846       fPoints->AddAt(points,trk);
847
848     }
849
850   }
851
852   delete [] coorE;
853   delete [] coorT;
854   delete [] ntrkE;
855   delete [] ntrkT;
856   delete [] limiE;
857   delete [] limiT;
858
859 }
860
861 //_____________________________________________________________________________
862 void AliTRD::MakeBranch(Option_t *option)
863 {
864   //
865   // Create Tree branches for the TRD digits.
866   //
867
868   Int_t  buffersize = 4000;
869   Char_t branchname[15];
870   sprintf(branchname,"%s",GetName());
871
872   const Char_t *cD = strstr(option,"D");
873
874   AliDetector::MakeBranch(option);
875
876   if (fDigits         && 
877       gAlice->TreeD() && 
878       cD) {
879     MakeBranchInTree(gAlice->TreeD(),branchname,&fDigits,buffersize,0);
880   }
881
882 }
883
884 //_____________________________________________________________________________
885 void AliTRD::ResetDigits()
886 {
887   //
888   // Reset number of digits and the digits array for this detector
889   //
890
891   fNdigits = 0;
892
893   if (fDigits) {
894     fDigits->Clear();
895   }
896
897 }
898
899 //_____________________________________________________________________________
900 void AliTRD::SetTreeAddress()
901 {
902   //
903   // Set the branch addresses for the trees.
904   //
905
906   if (fLoader->TreeH() && 
907       (fHits == 0x0)) {
908     fHits = new TClonesArray("AliTRDhit",405);
909   }
910   AliDetector::SetTreeAddress();
911
912 }
913
914 //_____________________________________________________________________________
915 Bool_t AliTRD::Raw2SDigits(AliRawReader *rawReader)
916 {
917   //
918   // Converts RAW data to SDigits
919   //
920
921   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
922   if (!loader) {
923     AliError("Can not get TRD loader from Run Loader");
924     return kFALSE;
925   }
926     
927   TTree *tree = 0;
928   tree = loader->TreeS();
929   if (!tree) {
930     loader->MakeTree("S");
931     tree = loader->TreeS();
932   }
933
934   AliTRDrawData       *rawdata        = new AliTRDrawData();
935   AliTRDdigitsManager *sdigitsManager = rawdata->Raw2Digits(rawReader);
936   if (sdigitsManager) {
937     sdigitsManager->SetSDigits(kTRUE);
938     sdigitsManager->MakeBranch(tree);
939     sdigitsManager->WriteDigits();
940     return kTRUE;
941   } 
942   else {
943     return kFALSE;
944   }
945
946 }
947
948 //_____________________________________________________________________________
949 AliTRD &AliTRD::operator=(const AliTRD &trd)
950 {
951   //
952   // Assignment operator
953   //
954
955   if (this != &trd) {
956     ((AliTRD &) trd).Copy(*this);
957   }
958
959   return *this;
960
961 }