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