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