Bug fixes, warnings only in case of debug option, more comments (P.Skowronski)
[u/mrichter/AliRoot.git] / TRD / AliTRD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  Transition Radiation Detector                                            //
21 //  This class contains the basic functions for the Transition Radiation     //
22 //  Detector.                                                                //
23 //                                                                           //
24 ///////////////////////////////////////////////////////////////////////////////
25
26 #include <stdlib.h>
27 #include <Riostream.h>
28
29 #include <TFile.h>
30 #include <TGeometry.h>
31 #include <TLorentzVector.h>
32 #include <TMath.h>
33 #include <TNode.h>
34 #include <TPGON.h> 
35 #include <TParticle.h>
36 #include <TROOT.h>
37 #include <TTree.h>
38 #include <TVirtualMC.h>
39  
40 #include "AliConst.h"
41 #include "AliDigit.h"
42 #include "AliLoader.h"
43 #include "AliMagF.h"
44 #include "AliRun.h"
45 #include "AliTRD.h"
46 #include "AliTRDcluster.h"
47 #include "AliTRDclusterizer.h"
48 #include "AliTRDdigit.h"
49 #include "AliTRDdigitizer.h"
50 #include "AliTRDdigitsManager.h"
51 #include "AliTRDgeometryFull.h"
52 #include "AliTRDgeometryHole.h"
53 #include "AliTRDhit.h"
54 #include "AliTRDpoints.h"
55 #include "AliTRDrecPoint.h"
56 #include "AliTRDtrackHits.h"  
57 #include "AliTrackReference.h"
58
59 ClassImp(AliTRD)
60  
61 //_____________________________________________________________________________
62 AliTRD::AliTRD()
63 {
64   //
65   // Default constructor
66   //
67
68   fIshunt        = 0;
69   fGasMix        = 0;
70   fHits          = 0;
71   fDigits        = 0;
72
73   fRecPoints     = 0;
74   fNRecPoints    = 0;
75
76   fGeometry      = 0;
77
78   fGasDensity    = 0;
79   fFoilDensity   = 0;
80
81   fDrawTR        = 0;
82   fDisplayType   = 0;
83  
84   fTrackHits     = 0; 
85   fHitType       = 0; 
86
87 }
88  
89 //_____________________________________________________________________________
90 AliTRD::AliTRD(const char *name, const char *title)
91        : AliDetector(name,title)
92 {
93   //
94   // Standard constructor for the TRD
95   //
96
97   // Check that FRAME is there otherwise we have no place where to
98   // put TRD
99   AliModule* frame = gAlice->GetModule("FRAME");
100   if (!frame) {
101     Error("Ctor","TRD needs FRAME to be present\n");
102     exit(1);
103   } 
104
105   // Define the TRD geometry according to the FRAME geometry
106   if      (frame->IsVersion() == 0) {
107     // Geometry with hole
108     fGeometry = new AliTRDgeometryHole();
109   }
110   else if (frame->IsVersion() == 1) {
111     // Geometry without hole
112     fGeometry = new AliTRDgeometryFull();
113   }
114   else {
115     Error("Ctor","Could not find valid FRAME version\n");
116     exit(1);
117   }
118
119   // Allocate the hit array
120   fHits           = new TClonesArray("AliTRDhit"     ,405);
121   gAlice->AddHitList(fHits);
122
123   // Allocate the digits array
124   fDigits         = 0;
125
126   // Allocate the rec point array
127   fRecPoints     = new TObjArray(400);
128   fNRecPoints    = 0;
129    
130   fIshunt        = 0;
131   fGasMix        = 1;
132
133   fGasDensity    = 0;
134   fFoilDensity   = 0;
135
136   fDrawTR        = 0;
137   fDisplayType   = 0;
138
139   fTrackHits     = 0;
140   fHitType       = 2;
141
142   SetMarkerColor(kWhite);   
143
144 }
145
146 //_____________________________________________________________________________
147 AliTRD::AliTRD(const AliTRD &trd)
148 {
149   //
150   // Copy constructor
151   //
152
153   ((AliTRD &) trd).Copy(*this);
154
155 }
156
157 //_____________________________________________________________________________
158 AliTRD::~AliTRD()
159 {
160   //
161   // TRD destructor
162   //
163
164   fIshunt = 0;
165
166   if (fGeometry) {
167     delete fGeometry;
168     fGeometry  = 0;
169   }
170   if (fHits) {
171     delete fHits;
172     fHits      = 0;
173   }
174   if (fRecPoints) {
175     delete fRecPoints;
176     fRecPoints = 0;
177   }
178   if (fTrackHits) {
179     delete fTrackHits;
180     fTrackHits = 0;
181   }
182
183 }
184
185 //_____________________________________________________________________________
186 void AliTRD::AddCluster(Float_t *pos, Int_t det, Float_t amp
187                       , Int_t *tracks, Float_t *sig, Int_t iType)
188 {
189   //
190   // Add a cluster for the TRD
191   //
192
193   AliTRDcluster *c = new AliTRDcluster();
194
195   c->SetDetector(det);
196   c->AddTrackIndex(tracks);
197   c->SetQ(amp);
198   c->SetY(pos[0]);
199   c->SetZ(pos[1]);
200   c->SetSigmaY2(sig[0]);   
201   c->SetSigmaZ2(sig[1]);
202   c->SetLocalTimeBin(((Int_t) pos[2]));
203
204   switch (iType) {
205   case 0:
206     c->Set2pad();
207     break;
208   case 1:
209     c->Set3pad();
210     break;
211   case 2:
212     c->Set4pad();
213     break;
214   case 3:
215     c->Set5pad();
216     break;
217   case 4:
218     c->SetLarge();
219     break;
220   };
221
222   fRecPoints->Add(c);
223
224 }
225
226 //_____________________________________________________________________________
227 void AliTRD::Hits2Digits()
228 {
229   //
230   // Create digits
231   //
232   if (!fLoader->TreeH()) fLoader->LoadHits("read");
233   SetTreeAddress();
234   
235   AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
236                                                   ,"TRD digitizer class");
237   digitizer->SetDebug(GetDebug());
238   digitizer->SetEvent(gAlice->GetEvNumber());
239   
240   // Initialization
241   digitizer->InitDetector();
242     
243   // Create the digits
244   digitizer->MakeDigits();
245   
246   // Write the digits into the input file
247
248   if (!fLoader->TreeD()) fLoader->MakeTree("D");
249   if (digitizer->MakeBranch(fLoader->TreeD())) {
250
251     digitizer->WriteDigits();
252
253     // Save the digitizer class in the AliROOT 
254     digitizer->Write();
255
256   }
257
258 }
259
260 //_____________________________________________________________________________
261 void AliTRD::Hits2SDigits()
262 {
263   //
264   // Create summable digits
265   //
266   if (!fLoader->TreeH()) fLoader->LoadHits("read");
267   SetTreeAddress();
268
269   AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
270                                                   ,"TRD digitizer class");
271   digitizer->SetDebug(GetDebug());
272
273   // For the summable digits
274   digitizer->SetSDigits(kTRUE);
275   digitizer->SetEvent(gAlice->GetEvNumber());
276
277   // Initialization
278   digitizer->InitDetector();
279     
280   // Create the TRD s-digits branch
281   digitizer->MakeDigits();
282   
283   // Write the digits into the input file
284   if (!fLoader->TreeS()) fLoader->MakeTree("S");
285   
286   if (digitizer->MakeBranch(fLoader->TreeS())) {
287
288     digitizer->WriteDigits();
289
290     // Save the digitizer class in the AliROOT 
291     digitizer->Write();
292
293   }
294
295 }
296
297 //_____________________________________________________________________________
298 void AliTRD::SDigits2Digits()
299 {
300   //
301   // Create final digits from summable digits
302   //
303
304    // Create the TRD digitizer
305   AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
306                                                   ,"TRD digitizer class");  
307   digitizer->SetDebug(GetDebug());
308
309   // Set the parameter
310   digitizer->SetEvent(gAlice->GetEvNumber());
311
312   // Initialization
313   digitizer->InitDetector();
314
315   // Read the s-digits via digits manager
316   AliTRDdigitsManager *sdigitsManager = new AliTRDdigitsManager();
317  
318   sdigitsManager->SetDebug(GetDebug());
319   sdigitsManager->SetSDigits(kTRUE);
320   sdigitsManager->CreateArrays();
321   
322   if (!fLoader->TreeS()) 
323     if (fLoader->LoadSDigits("read"))
324      {
325        Error("SDigits2Digits","Error while reading SDigits for event %d",gAlice->GetEvNumber());
326        return;
327      }
328   if (!fLoader->TreeS()) return;
329   
330   sdigitsManager->ReadDigits(fLoader->TreeS());
331
332   // Add the s-digits to the input list 
333   digitizer->AddSDigitsManager(sdigitsManager);
334
335   // Convert the s-digits to normal digits
336   digitizer->SDigits2Digits();
337
338   // Store the digits
339   if (!fLoader->TreeD()) fLoader->MakeTree("D");
340   if (digitizer->MakeBranch(fLoader->TreeD())){
341     digitizer->WriteDigits();
342   }
343
344 }
345
346 //_____________________________________________________________________________
347 void AliTRD::AddHit(Int_t track, Int_t det, Float_t *hits, Int_t q
348                   , Bool_t inDrift)
349 {
350   //
351   // Add a hit for the TRD
352   // 
353   // The data structure is set according to fHitType:
354   //   bit0: standard TClonesArray
355   //   bit1: compressed trackHits structure
356   //
357
358   if (fHitType & 1) {
359     TClonesArray &lhits = *fHits;
360     new(lhits[fNhits++]) AliTRDhit(fIshunt,track,det,hits,q);
361   }
362
363   if (fHitType > 1) {
364     AddHit2(track,det,hits,q,inDrift);
365   }
366
367 }
368
369 //_____________________________________________________________________________
370 void AliTRD::BuildGeometry()
371 {
372   //
373   // Create the ROOT TNode geometry for the TRD
374   //
375
376   TNode *node, *top;
377   TPGON *pgon;
378
379   Float_t rmin, rmax;
380   Float_t zmax1, zmax2;
381
382   Int_t   iPlan;
383  
384   const Int_t kColorTRD = 46;
385   
386   // Find the top node alice
387   top = gAlice->GetGeometry()->GetNode("alice");
388   
389   if      (fDisplayType == 0) {
390
391     pgon = new TPGON("S_TRD","TRD","void",0,360,AliTRDgeometry::Nsect(),4);
392     rmin = AliTRDgeometry::Rmin();
393     rmax = AliTRDgeometry::Rmax();
394     pgon->DefineSection(0,-AliTRDgeometry::Zmax1(),rmax,rmax);
395     pgon->DefineSection(1,-AliTRDgeometry::Zmax2(),rmin,rmax);
396     pgon->DefineSection(2, AliTRDgeometry::Zmax2(),rmin,rmax);
397     pgon->DefineSection(3, AliTRDgeometry::Zmax1(),rmax,rmax);
398     top->cd();
399     node = new TNode("TRD","TRD","S_TRD",0,0,0,"");
400     node->SetLineColor(kColorTRD);
401     fNodes->Add(node);
402
403   }
404   else if (fDisplayType == 1) {
405
406     Char_t name[7];
407
408     Float_t slope = (AliTRDgeometry::Zmax1() - AliTRDgeometry::Zmax2())
409                   / (AliTRDgeometry::Rmax()  - AliTRDgeometry::Rmin());
410
411     rmin  = AliTRDgeometry::Rmin() + AliTRDgeometry::CraHght();
412     rmax  = rmin                   + AliTRDgeometry::CdrHght();
413
414     Float_t thickness = rmin - AliTRDgeometry::Rmin();
415     zmax2 = AliTRDgeometry::Zmax2() + slope * thickness;
416     zmax1 = zmax2 + slope * AliTRDgeometry::DrThick();
417
418     for (iPlan = 0; iPlan < AliTRDgeometry::Nplan(); iPlan++) {
419
420       sprintf(name,"S_TR1%d",iPlan);
421       pgon  = new TPGON(name,"TRD","void",0,360,AliTRDgeometry::Nsect(),4);
422       pgon->DefineSection(0,-zmax1,rmax,rmax);
423       pgon->DefineSection(1,-zmax2,rmin,rmax);
424       pgon->DefineSection(2, zmax2,rmin,rmax);
425       pgon->DefineSection(3, zmax1,rmax,rmax);
426       top->cd();
427       node = new TNode("TRD","TRD",name,0,0,0,"");
428       node->SetLineColor(kColorTRD);
429       fNodes->Add(node);
430
431       Float_t height = AliTRDgeometry::Cheight() + AliTRDgeometry::Cspace(); 
432       rmin  = rmin  + height;
433       rmax  = rmax  + height;
434       zmax1 = zmax1 + slope * height;
435       zmax2 = zmax2 + slope * height;
436
437     }
438
439     thickness += AliTRDgeometry::DrThick();
440     rmin  = AliTRDgeometry::Rmin() + thickness;
441     rmax  = rmin + AliTRDgeometry::AmThick();
442     zmax2 = AliTRDgeometry::Zmax2() + slope * thickness;
443     zmax1 = zmax2 + slope * AliTRDgeometry::AmThick();
444
445     for (iPlan = 0; iPlan < AliTRDgeometry::Nplan(); iPlan++) {
446
447       sprintf(name,"S_TR2%d",iPlan);
448       pgon  = new TPGON(name,"TRD","void",0,360,AliTRDgeometry::Nsect(),4);
449       pgon->DefineSection(0,-zmax1,rmax,rmax);
450       pgon->DefineSection(1,-zmax2,rmin,rmax);
451       pgon->DefineSection(2, zmax2,rmin,rmax);
452       pgon->DefineSection(3, zmax1,rmax,rmax);
453       top->cd();
454       node = new TNode("TRD","TRD",name,0,0,0,"");
455       node->SetLineColor(kColorTRD);
456       fNodes->Add(node);
457
458       Float_t height = AliTRDgeometry::Cheight() + AliTRDgeometry::Cspace(); 
459       rmin  = rmin  + height;
460       rmax  = rmax  + height;
461       zmax1 = zmax1 + slope * height;
462       zmax2 = zmax2 + slope * height;
463
464     }
465
466   }
467
468 }
469  
470 //_____________________________________________________________________________
471 void AliTRD::Copy(TObject &trd)
472 {
473   //
474   // Copy function
475   //
476
477   ((AliTRD &) trd).fGasMix      = fGasMix;
478   ((AliTRD &) trd).fGeometry    = fGeometry;       
479   ((AliTRD &) trd).fRecPoints   = fRecPoints;
480   ((AliTRD &) trd).fNRecPoints  = fNRecPoints;
481   ((AliTRD &) trd).fGasDensity  = fGasDensity;
482   ((AliTRD &) trd).fFoilDensity = fFoilDensity;
483   ((AliTRD &) trd).fDrawTR      = fDrawTR;
484   ((AliTRD &) trd).fDisplayType = fDisplayType;
485   ((AliTRD &) trd).fHitType     = fHitType;
486
487   //AliDetector::Copy(trd);
488
489 }
490
491 //_____________________________________________________________________________
492 void AliTRD::CreateGeometry()
493 {
494   //
495   // Creates the volumes for the TRD chambers
496   //
497
498   // Check that FRAME is there otherwise we have no place where to put the TRD
499   AliModule* frame = gAlice->GetModule("FRAME");
500   if (!frame) {
501     printf(" The TRD needs the FRAME to be defined first\n");
502     return;
503   }
504
505   fGeometry->CreateGeometry(fIdtmed->GetArray() - 1299);
506
507 }
508  
509 //_____________________________________________________________________________
510 void AliTRD::CreateMaterials()
511 {
512   //
513   // Create the materials for the TRD
514   // Origin Y.Foka
515   //
516
517   Int_t   isxfld = gAlice->Field()->Integ();
518   Float_t sxmgmx = gAlice->Field()->Max();
519   
520   // For polyethilene (CH2) 
521   Float_t ape[2] = { 12., 1. };
522   Float_t zpe[2] = {  6., 1. };
523   Float_t wpe[2] = {  1., 2. };
524   Float_t dpe    = 0.95;
525
526   // For mylar (C5H4O2) 
527   Float_t amy[3] = { 12., 1., 16. };
528   Float_t zmy[3] = {  6., 1.,  8. };
529   Float_t wmy[3] = {  5., 4.,  2. };
530   Float_t dmy    = 1.39;
531
532   // For CO2 
533   Float_t aco[2] = { 12., 16. };
534   Float_t zco[2] = {  6.,  8. };
535   Float_t wco[2] = {  1.,  2. };
536   Float_t dco    = 0.001977;
537
538   // For water
539   Float_t awa[2] = {  1., 16. };
540   Float_t zwa[2] = {  1.,  8. };
541   Float_t wwa[2] = {  2.,  1. };
542   Float_t dwa    = 1.0;
543
544   // For isobutane (C4H10)
545   Float_t ais[2] = { 12.,  1. };
546   Float_t zis[2] = {  6.,  1. };
547   Float_t wis[2] = {  4., 10. };
548   Float_t dis    = 0.00267;
549
550   // For plexiglas (C5H8O2)
551   Float_t apg[3] = { 12.011 ,  1.0    , 15.9994 };
552   Float_t zpg[3] = {  6.0   ,  1.0    ,  8.0    };
553   Float_t wpg[3] = {  5.0   ,  8.0    ,  2.0    };
554   Float_t dpg    = 1.18; 
555
556   // For Xe/CO2-gas-mixture 
557   // Xe-content of the Xe/CO2-mixture (85% / 15%) 
558   Float_t fxc    = .85;
559   // Xe-content of the Xe/Isobutane-mixture (97% / 3%) 
560   Float_t fxi    = .97;
561   Float_t dxe    = .005858;
562   
563   // General tracking parameter
564   Float_t tmaxfd = -10.;
565   Float_t stemax = -1e10;
566   Float_t deemax = -0.1;
567   Float_t epsil  =  1e-4;
568   Float_t stmin  = -0.001;
569   
570   Float_t absl, radl, d, buf[1];
571   Float_t agm[2], zgm[2], wgm[2];
572   Float_t dgm1, dgm2;
573   Int_t   nbuf;
574   
575   //////////////////////////////////////////////////////////////////////////
576   //     Define Materials 
577   //////////////////////////////////////////////////////////////////////////
578
579   AliMaterial( 1, "Al"   ,  26.98, 13.0, 2.7     ,     8.9 ,    37.2);
580   AliMaterial( 2, "Air"  ,  14.61,  7.3, 0.001205, 30420.0 , 67500.0);
581   AliMaterial( 4, "Xe"   , 131.29, 54.0, dxe     ,  1447.59,     0.0);
582   AliMaterial( 5, "Cu"   ,  63.54, 29.0, 8.96    ,     1.43,    14.8);
583   AliMaterial( 6, "C"    ,  12.01,  6.0, 2.265   ,    18.8 ,    74.4);
584   AliMaterial(12, "G10"  ,  20.00, 10.0, 1.7     ,    19.4 ,   999.0);
585   AliMaterial(15, "Sn"   , 118.71, 50.0, 7.31    ,     1.21,    14.8);
586   AliMaterial(16, "Si"   ,  28.09, 14.0, 2.33    ,     9.36,    37.2);
587   AliMaterial(17, "Epoxy",  17.75,  8.9, 1.8     ,    21.82,   999.0);
588
589   // Mixtures 
590   AliMixture(3, "Polyethilene",   ape, zpe, dpe, -2, wpe);
591   AliMixture(7, "Mylar",          amy, zmy, dmy, -3, wmy);
592   AliMixture(8, "CO2",            aco, zco, dco, -2, wco);
593   AliMixture(9, "Isobutane",      ais, zis, dis, -2, wis);
594   AliMixture(13,"Water",          awa, zwa, dwa, -2, wwa);
595   AliMixture(14,"Plexiglas",      apg, zpg, dpg, -3, wpg);
596
597   // Gas mixtures
598   Char_t namate[21];
599   // Xe/CO2-mixture
600   // Get properties of Xe 
601   gMC->Gfmate((*fIdmate)[4], namate, agm[0], zgm[0], d, radl, absl, buf, nbuf);
602   // Get properties of CO2 
603   gMC->Gfmate((*fIdmate)[8], namate, agm[1], zgm[1], d, radl, absl, buf, nbuf);
604   // Create gas mixture 
605   wgm[0] = fxc;
606   wgm[1] = 1. - fxc;
607   dgm1   = wgm[0] * dxe + wgm[1] * dco;
608   AliMixture(10, "Gas mixture 1", agm, zgm, dgm1,  2, wgm);
609   // Xe/Isobutane-mixture
610   // Get properties of Xe 
611   gMC->Gfmate((*fIdmate)[4], namate, agm[0], zgm[0], d, radl, absl, buf, nbuf);
612   // Get properties of Isobutane
613   gMC->Gfmate((*fIdmate)[9], namate, agm[1], zgm[1], d, radl, absl, buf, nbuf);
614   // Create gas mixture 
615   wgm[0] = fxi;
616   wgm[1] = 1. - fxi;
617   dgm2   = wgm[0] * dxe + wgm[1] * dis;
618   AliMixture(11, "Gas mixture 2", agm, zgm, dgm2,  2, wgm);
619  
620   //////////////////////////////////////////////////////////////////////////
621   //     Tracking Media Parameters 
622   //////////////////////////////////////////////////////////////////////////
623
624   // Al Frame 
625   AliMedium(1, "Al Frame",   1, 0, isxfld, sxmgmx
626                 , tmaxfd, stemax, deemax, epsil, stmin);
627   // Air 
628   AliMedium(2, "Air",        2, 0, isxfld, sxmgmx
629                 , tmaxfd, stemax, deemax, epsil, stmin);
630   // Polyethilene 
631   AliMedium(3, "Radiator",   3, 0, isxfld, sxmgmx
632                 , tmaxfd, stemax, deemax, epsil, stmin);
633   // Xe 
634   AliMedium(4, "Xe",         4, 1, isxfld, sxmgmx
635                 , tmaxfd, stemax, deemax, epsil, stmin);
636   // Cu pads 
637   AliMedium(5, "Padplane",   5, 1, isxfld, sxmgmx
638                 , tmaxfd, stemax, deemax, epsil, stmin);
639   // Fee + cables 
640   AliMedium(6, "Readout",    1, 0, isxfld, sxmgmx
641                 , tmaxfd, stemax, deemax, epsil, stmin);
642   // C frame 
643   AliMedium(7, "C Frame",    6, 0, isxfld, sxmgmx
644                 , tmaxfd, stemax, deemax, epsil, stmin);
645   // Mylar foils 
646   AliMedium(8, "Mylar",      7, 0, isxfld, sxmgmx
647                 , tmaxfd, stemax, deemax, epsil, stmin);
648   if (fGasMix == 1) {
649     // Gas-mixture (Xe/CO2) 
650     AliMedium(9, "Gas-mix",   10, 1, isxfld, sxmgmx
651                   , tmaxfd, stemax, deemax, epsil, stmin);
652   }
653   else {
654     // Gas-mixture (Xe/Isobutane) 
655     AliMedium(9, "Gas-mix",   11, 1, isxfld, sxmgmx
656                   , tmaxfd, stemax, deemax, epsil, stmin);
657   }
658   // Nomex-honeycomb (use carbon for the time being) 
659   AliMedium(10, "Nomex",      6, 0, isxfld, sxmgmx
660                 , tmaxfd, stemax, deemax, epsil, stmin);
661   // Kapton foils (use Mylar for the time being) 
662   AliMedium(11, "Kapton",     7, 0, isxfld, sxmgmx
663                 , tmaxfd, stemax, deemax, epsil, stmin);
664   // Gas-filling of the radiator 
665   AliMedium(12, "CO2",        8, 0, isxfld, sxmgmx
666                 , tmaxfd, stemax, deemax, epsil, stmin);
667   // G10-plates
668   AliMedium(13, "G10-plates",12, 0, isxfld, sxmgmx
669                 , tmaxfd, stemax, deemax, epsil, stmin);
670   // Cooling water
671   AliMedium(14, "Water",     13, 0, isxfld, sxmgmx
672                 , tmaxfd, stemax, deemax, epsil, stmin);
673   // Rohacell (plexiglas) for the radiator
674   AliMedium(15, "Rohacell",  14, 0, isxfld, sxmgmx
675                 , tmaxfd, stemax, deemax, epsil, stmin);
676   // Al layer in MCMs
677   AliMedium(16, "MCM-Al"  ,   1, 0, isxfld, sxmgmx
678                 , tmaxfd, stemax, deemax, epsil, stmin);
679   // Sn layer in MCMs
680   AliMedium(17, "MCM-Sn"  ,  15, 0, isxfld, sxmgmx
681                 , tmaxfd, stemax, deemax, epsil, stmin);
682   // Cu layer in MCMs
683   AliMedium(18, "MCM-Cu"  ,   5, 0, isxfld, sxmgmx
684                 , tmaxfd, stemax, deemax, epsil, stmin);
685   // G10 layer in MCMs
686   AliMedium(19, "MCM-G10" ,  12, 0, isxfld, sxmgmx
687                 , tmaxfd, stemax, deemax, epsil, stmin);
688   // Si in readout chips
689   AliMedium(20, "Chip-Si" ,  16, 0, isxfld, sxmgmx
690                 , tmaxfd, stemax, deemax, epsil, stmin);
691   // Epoxy in readout chips
692   AliMedium(21, "Chip-Ep" ,  17, 0, isxfld, sxmgmx
693                 , tmaxfd, stemax, deemax, epsil, stmin);
694   // PE in connectors
695   AliMedium(22, "Conn-PE" ,   3, 0, isxfld, sxmgmx
696                 , tmaxfd, stemax, deemax, epsil, stmin);
697   // Cu in connectors
698   AliMedium(23, "Chip-Cu" ,   5, 0, isxfld, sxmgmx
699                 , tmaxfd, stemax, deemax, epsil, stmin);
700   // Al of cooling pipes
701   AliMedium(24, "Cooling" ,   1, 0, isxfld, sxmgmx
702                 , tmaxfd, stemax, deemax, epsil, stmin);
703
704   // Save the density values for the TRD absorbtion
705   fFoilDensity = dmy;
706   if (fGasMix == 1)
707     fGasDensity = dgm1;
708   else
709     fGasDensity = dgm2;
710
711 }
712
713 //_____________________________________________________________________________
714 void AliTRD::DrawModule() const
715 {
716   //
717   // Draw a shaded view of the Transition Radiation Detector version 0
718   //
719
720   // Set everything unseen
721   gMC->Gsatt("*"   ,"SEEN",-1);
722   
723   // Set ALIC mother transparent
724   gMC->Gsatt("ALIC","SEEN", 0);
725   
726   // Set the volumes visible
727   if (fGeometry->IsVersion() == 0) {
728     gMC->Gsatt("B071","SEEN", 0);
729     gMC->Gsatt("B074","SEEN", 0);
730     gMC->Gsatt("B075","SEEN", 0);
731     gMC->Gsatt("B077","SEEN", 0);
732     gMC->Gsatt("BTR1","SEEN", 0);
733     gMC->Gsatt("BTR2","SEEN", 0);
734     gMC->Gsatt("BTR3","SEEN", 0);
735     gMC->Gsatt("UTR1","SEEN", 0);
736     gMC->Gsatt("UTR2","SEEN", 0);
737     gMC->Gsatt("UTR3","SEEN", 0);
738   }
739   else {
740     gMC->Gsatt("B071","SEEN", 0);
741     gMC->Gsatt("B074","SEEN", 0);
742     gMC->Gsatt("B075","SEEN", 0);
743     gMC->Gsatt("B077","SEEN", 0);
744     gMC->Gsatt("BTR1","SEEN", 0);
745     gMC->Gsatt("BTR2","SEEN", 0);
746     gMC->Gsatt("BTR3","SEEN", 0);
747     gMC->Gsatt("UTR1","SEEN", 0);
748     if (fGeometry->GetPHOShole())
749       gMC->Gsatt("UTR2","SEEN", 0);
750     if (fGeometry->GetRICHhole())
751       gMC->Gsatt("UTR3","SEEN", 0);
752   }
753 //   gMC->Gsatt("UCII","SEEN", 0);
754 //   gMC->Gsatt("UCIM","SEEN", 0);
755 //   gMC->Gsatt("UCIO","SEEN", 0);
756 //   gMC->Gsatt("UL02","SEEN", 1);
757 //   gMC->Gsatt("UL05","SEEN", 1);
758 //   gMC->Gsatt("UL06","SEEN", 1);
759   
760   gMC->Gdopt("hide", "on");
761   gMC->Gdopt("shad", "on");
762   gMC->Gsatt("*", "fill", 7);
763   gMC->SetClipBox(".");
764   gMC->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
765   gMC->DefaultRange();
766   gMC->Gdraw("alic", 40, 30, 0, 12, 9.4, .021, .021);
767   gMC->Gdhead(1111, "Transition Radiation Detector");
768   gMC->Gdman(18, 4, "MAN");
769
770 }
771
772 //_____________________________________________________________________________
773 Int_t AliTRD::DistancetoPrimitive(Int_t , Int_t ) const
774 {
775   //
776   // Distance between the mouse and the TRD detector on the screen
777   // Dummy routine
778   
779   return 9999;
780
781 }
782  
783 //_____________________________________________________________________________
784 void AliTRD::Init()
785 {
786   //
787   // Initialize the TRD detector after the geometry has been created
788   //
789
790   Int_t i;
791
792   if (fDebug) {
793     printf("\n%s: ",ClassName());
794     for (i = 0; i < 35; i++) printf("*");
795     printf(" TRD_INIT ");
796     for (i = 0; i < 35; i++) printf("*");
797     printf("\n");
798   }
799
800   if      (fGeometry->IsVersion() == 0) {
801     printf("%s: Geometry for spaceframe with holes initialized\n",ClassName());
802   }
803   else if (fGeometry->IsVersion() == 1) {
804     printf("%s: Geometry for spaceframe without holes initialized\n",ClassName());
805     if (fGeometry->GetPHOShole())
806       printf("%s: Leave space in front of PHOS free\n",ClassName());
807     if (fGeometry->GetRICHhole())
808       printf("%s: Leave space in front of RICH free\n",ClassName());
809   }
810   
811   if (fGasMix == 1) {
812     printf("%s: Gas Mixture: 85%% Xe + 15%% CO2\n",ClassName());
813   }
814   else {
815     printf("%s: Gas Mixture: 97%% Xe + 3%% Isobutane\n",ClassName());
816   }
817
818 }
819
820 //_____________________________________________________________________________
821 void AliTRD::LoadPoints(Int_t track)
822 {
823   //
824   // Store x, y, z of all hits in memory.
825   // Hit originating from TR photons are given a different color
826   //
827
828   //if (!fDrawTR) {
829   //  AliDetector::LoadPoints(track);
830   //  return;
831   //}
832
833   if ((fHits == 0) && (fTrackHits == 0)) return;
834
835   Int_t nhits;
836   if (fHitType < 2) {
837     nhits = fHits->GetEntriesFast();
838   }
839   else {
840     nhits = fTrackHits->GetEntriesFast();
841   } 
842   if (nhits == 0) return;
843
844   Int_t tracks = gAlice->GetNtrack();
845   if (fPoints == 0) fPoints = new TObjArray(tracks);
846
847   AliTRDhit *ahit;
848   
849   Int_t    *ntrkE = new Int_t[tracks];
850   Int_t    *ntrkT = new Int_t[tracks];
851   Int_t    *limiE = new Int_t[tracks];
852   Int_t    *limiT = new Int_t[tracks];
853   Float_t **coorE = new Float_t*[tracks];
854   Float_t **coorT = new Float_t*[tracks];
855   for(Int_t i = 0; i < tracks; i++) {
856     ntrkE[i] = 0;
857     ntrkT[i] = 0;
858     coorE[i] = 0;
859     coorT[i] = 0;
860     limiE[i] = 0;
861     limiT[i] = 0;
862   }
863   
864   AliTRDpoints *points = 0;
865   Float_t      *fp     = 0;
866   Int_t         trk;
867   Int_t         chunk  = nhits / 4 + 1;
868
869   // Loop over all the hits and store their position
870   ahit = (AliTRDhit *) FirstHit(-1);
871   while (ahit) {
872
873     // dEdx hits
874     if (ahit->GetCharge() >= 0) {
875
876       trk = ahit->GetTrack();
877       if (ntrkE[trk] == limiE[trk]) {
878         // Initialise a new track
879         fp = new Float_t[3*(limiE[trk]+chunk)];
880         if (coorE[trk]) {
881           memcpy(fp,coorE[trk],sizeof(Float_t)*3*limiE[trk]);
882           delete [] coorE[trk];
883         }
884         limiE[trk] += chunk;
885         coorE[trk]  = fp;
886       } 
887       else {
888         fp = coorE[trk];
889       }
890       fp[3*ntrkE[trk]  ] = ahit->X();
891       fp[3*ntrkE[trk]+1] = ahit->Y();
892       fp[3*ntrkE[trk]+2] = ahit->Z();
893       ntrkE[trk]++;
894
895     }
896     // TR photon hits
897     else if ((ahit->GetCharge() < 0) && (fDrawTR)) {
898
899       trk = ahit->GetTrack();
900       if (ntrkT[trk] == limiT[trk]) {
901         // Initialise a new track
902         fp = new Float_t[3*(limiT[trk]+chunk)];
903         if (coorT[trk]) {
904           memcpy(fp,coorT[trk],sizeof(Float_t)*3*limiT[trk]);
905           delete [] coorT[trk];
906         }
907         limiT[trk] += chunk;
908         coorT[trk]  = fp;
909       } 
910       else {
911         fp = coorT[trk];
912       }
913       fp[3*ntrkT[trk]  ] = ahit->X();
914       fp[3*ntrkT[trk]+1] = ahit->Y();
915       fp[3*ntrkT[trk]+2] = ahit->Z();
916       ntrkT[trk]++;
917
918     }
919
920     ahit = (AliTRDhit *) NextHit();
921
922   }
923
924   for (trk = 0; trk < tracks; ++trk) {
925
926     if (ntrkE[trk] || ntrkT[trk]) {
927
928       points = new AliTRDpoints();
929       points->SetDetector(this);
930       points->SetParticle(trk);
931
932       // Set the dEdx points
933       if (ntrkE[trk]) {
934         points->SetMarkerColor(GetMarkerColor());
935         points->SetMarkerSize(GetMarkerSize());
936         points->SetPolyMarker(ntrkE[trk],coorE[trk],GetMarkerStyle());
937         delete [] coorE[trk];
938         coorE[trk] = 0;
939       }
940
941       // Set the TR photon points
942       if (ntrkT[trk]) {
943         points->SetTRpoints(ntrkT[trk],coorT[trk]);
944         delete [] coorT[trk];
945         coorT[trk] = 0;
946       }
947
948       fPoints->AddAt(points,trk);
949
950     }
951
952   }
953
954   delete [] coorE;
955   delete [] coorT;
956   delete [] ntrkE;
957   delete [] ntrkT;
958   delete [] limiE;
959   delete [] limiT;
960
961 }
962
963 //_____________________________________________________________________________
964 void AliTRD::MakeBranch(Option_t* option)
965 {
966   //
967   // Create Tree branches for the TRD digits.
968   //
969
970   Int_t  buffersize = 4000;
971   Char_t branchname[15];
972   sprintf(branchname,"%s",GetName());
973
974   const char *cD = strstr(option,"D");
975
976   AliDetector::MakeBranch(option);
977
978   if (fDigits && gAlice->TreeD() && cD) {
979     MakeBranchInTree(gAlice->TreeD(),branchname,&fDigits,buffersize,0);
980   }     
981
982   if (fHitType > 1) {
983     MakeBranch2(option,0); 
984   }
985
986 }
987
988 //_____________________________________________________________________________
989 void AliTRD::ResetDigits()
990 {
991   //
992   // Reset number of digits and the digits array for this detector
993   //
994
995   fNdigits = 0;
996   if (fDigits) fDigits->Clear();
997
998 }
999
1000 //_____________________________________________________________________________
1001 void AliTRD::ResetRecPoints()
1002 {
1003   //
1004   // Reset number of reconstructed points and the point array
1005   //
1006
1007   if (fRecPoints) {
1008     fNRecPoints = 0;
1009     Int_t nentr = fRecPoints->GetEntriesFast();
1010     for (Int_t i = 0; i < nentr; i++) delete fRecPoints->RemoveAt(i);
1011   }
1012
1013 }
1014
1015 //_____________________________________________________________________________
1016 void AliTRD::SetTreeAddress()
1017 {
1018   //
1019   // Set the branch addresses for the trees.
1020   //
1021
1022   Char_t branchname[15];
1023
1024   if ( fLoader->TreeH() && (fHits == 0x0))  fHits = new TClonesArray("AliTRDhit",405);
1025   AliDetector::SetTreeAddress();
1026
1027   TBranch *branch;
1028   TTree   *treeR = fLoader->TreeR();
1029
1030   if (treeR) {
1031     sprintf(branchname,"%scluster",GetName());
1032     if (fRecPoints == 0x0) fRecPoints  = new TObjArray(400);
1033     if (fRecPoints) {
1034       branch = treeR->GetBranch(branchname);
1035       if (branch) {
1036         branch->SetAddress(&fRecPoints);
1037       }
1038     }
1039   }
1040
1041   if (fHitType > 0) {
1042     SetTreeAddress2();    
1043   }
1044
1045 }
1046
1047 //_____________________________________________________________________________
1048 void AliTRD::SetGasMix(Int_t imix)
1049 {
1050   //
1051   // Defines the gas mixture (imix=0:  Xe/Isobutane imix=1: Xe/CO2)
1052   //
1053   
1054   if ((imix < 0) || (imix > 1)) {
1055     printf("Wrong input value: %d\n",imix);
1056     printf("Use standard setting\n");
1057     fGasMix = 1;
1058     return;
1059   }
1060
1061   fGasMix = imix;
1062
1063 }
1064
1065 //_____________________________________________________________________________
1066 void AliTRD::SetPHOShole()
1067 {
1068   //
1069   // Selects a geometry with a hole in front of the PHOS
1070   //
1071
1072   fGeometry->SetPHOShole();
1073
1074 }
1075
1076 //_____________________________________________________________________________
1077 void AliTRD::SetRICHhole()
1078 {
1079   //
1080   // Selects a geometry with a hole in front of the RICH
1081   //
1082
1083   fGeometry->SetRICHhole();
1084
1085 }
1086
1087 //_____________________________________________________________________________
1088 AliTRD &AliTRD::operator=(const AliTRD &trd)
1089 {
1090   //
1091   // Assignment operator
1092   //
1093
1094   if (this != &trd) ((AliTRD &) trd).Copy(*this);
1095   return *this;
1096
1097
1098
1099 //_____________________________________________________________________________
1100 void AliTRD::FinishPrimary()
1101 {
1102   //
1103   // Store the hits in the containers after all primaries are finished
1104   //
1105
1106   if (fTrackHits) { 
1107     fTrackHits->FlushHitStack();
1108   }
1109
1110 }
1111
1112 //_____________________________________________________________________________
1113 void AliTRD::RemapTrackHitIDs(Int_t *map)
1114 {
1115   //
1116   // Remap the track IDs
1117   //
1118
1119   if (!fTrackHits) {
1120     return;
1121   }
1122
1123   if (fTrackHits) {
1124     TClonesArray *arr = fTrackHits->GetArray();;
1125     for (Int_t i = 0; i < arr->GetEntriesFast(); i++){
1126       AliTrackHitsParamV2 *info = (AliTrackHitsParamV2 *) (arr->At(i));
1127       info->fTrackID = map[info->fTrackID];
1128     }
1129   }
1130
1131 }
1132
1133 //_____________________________________________________________________________
1134 void AliTRD::ResetHits()
1135 {
1136   //
1137   // Reset the hits
1138   //
1139
1140   AliDetector::ResetHits();
1141   if (fTrackHits) {
1142     fTrackHits->Clear();
1143   }
1144
1145 }
1146
1147 //_____________________________________________________________________________
1148 AliHit* AliTRD::FirstHit(Int_t track)
1149 {
1150   //
1151   // Return the first hit of a track
1152   //
1153
1154   if (fHitType > 1) {
1155     return FirstHit2(track);
1156   }
1157
1158   return AliDetector::FirstHit(track);
1159
1160 }
1161
1162 //_____________________________________________________________________________
1163 AliHit* AliTRD::NextHit()
1164 {
1165   //
1166   // Returns the next hit of a track
1167   //
1168
1169   if (fHitType > 1) {
1170     return NextHit2();
1171   }
1172
1173   return AliDetector::NextHit();
1174
1175 }
1176
1177 //_____________________________________________________________________________
1178 AliHit* AliTRD::FirstHit2(Int_t track) 
1179 {
1180   //
1181   // Initializes the hit iterator.
1182   // Returns the address of the first hit of a track.
1183   // If <track> >= 0 the track is read from disk,
1184   // while if <track> < 0 the first hit of the current
1185   // track is returned.
1186   //
1187
1188   if (track >= 0) {
1189     gAlice->ResetHits();
1190     TreeH()->GetEvent(track);
1191   }
1192   
1193   if (fTrackHits) {
1194     fTrackHits->First();
1195     return (AliHit*) fTrackHits->GetHit();
1196   }
1197   else {
1198     return 0;
1199   }
1200
1201 }
1202
1203 //_____________________________________________________________________________
1204 AliHit* AliTRD::NextHit2()
1205 {
1206   //
1207   // Returns the next hit of the current track
1208   //
1209
1210   if (fTrackHits) {
1211     fTrackHits->Next();
1212     return (AliHit *) fTrackHits->GetHit();
1213   }
1214   else {
1215     return 0;
1216   }
1217
1218 }
1219
1220 //_____________________________________________________________________________
1221 void AliTRD::MakeBranch2(Option_t *option, const char *file)
1222 {
1223   //
1224   // Create a new branch in the current Root tree.
1225   // The branch of fHits is automatically split.
1226   //
1227
1228   if (fHitType < 2) {
1229     return;
1230   }
1231
1232   char branchname[10];
1233   sprintf(branchname,"%s2",GetName());
1234
1235   // Get the pointer to the header
1236   const char *cH = strstr(option,"H");
1237  
1238   if (!fTrackHits) {
1239     fTrackHits = new AliTRDtrackHits();
1240    }
1241
1242
1243   if (fTrackHits && TreeH() && cH) 
1244    {
1245     TreeH()->Branch(branchname,"AliTRDtrackHits",&fTrackHits,fBufferSize,99);
1246     Info("MakeBranch2","Making Branch %s for trackhits",branchname);
1247    }
1248 }
1249
1250 //_____________________________________________________________________________
1251 void AliTRD::SetTreeAddress2()
1252 {
1253   //
1254   // Set the branch address for the trackHits tree
1255   //
1256   if (GetDebug())  Info("SetTreeAddress2","");
1257
1258   TBranch *branch;
1259   char branchname[20];
1260   sprintf(branchname,"%s2",GetName());
1261   
1262   // Branch address for hit tree
1263   TTree *treeH = TreeH();
1264   if ((treeH) && (fHitType > 0)) {
1265     branch = treeH->GetBranch(branchname);
1266     if (branch) 
1267      {
1268        branch->SetAddress(&fTrackHits);
1269        if (GetDebug())
1270          Info("SetTreeAddress2","Success.");
1271      }
1272     else
1273      {
1274        if (GetDebug())
1275          Info("SetTreeAddress2","Can NOT get the branch %s",branchname);
1276      }
1277   }
1278
1279 }
1280
1281 //_____________________________________________________________________________
1282 void AliTRD::AddHit2(Int_t track, Int_t det, Float_t *hits, Int_t q
1283                    , Bool_t inDrift)
1284 {
1285   //
1286   // Add a hit to the list
1287   //
1288
1289   Int_t rtrack;
1290
1291   if (fIshunt) {
1292     Int_t primary = gAlice->GetPrimary(track);
1293     gAlice->Particle(primary)->SetBit(kKeepBit);
1294     rtrack = primary;
1295   } 
1296   else {
1297     rtrack = track;
1298     gAlice->FlagTrack(track);
1299   }
1300
1301   if ((fTrackHits) && (fHitType > 0)) {
1302     fTrackHits->AddHitTRD(det,rtrack,hits[0],hits[1],hits[2],q,inDrift);
1303   }
1304
1305 }
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644