Changes in digits IO. Add merging of summable digits
[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 /*
17 $Log$
18 Revision 1.33  2001/11/06 17:19:41  cblume
19 Add detailed geometry and simple simulator
20
21 Revision 1.32  2001/10/08 06:57:33  hristov
22 Branches for  TRD digits are created only during the digitisation
23
24 Revision 1.31  2001/08/30 09:30:30  hristov
25 The split level of branches is set to 99
26
27 Revision 1.30  2001/05/28 17:07:58  hristov
28 Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
29
30 Revision 1.29  2001/05/21 16:45:47  hristov
31 Last minute changes (C.Blume)
32
33 Revision 1.28  2001/05/16 14:57:27  alibrary
34 New files for folders and Stack
35
36 Revision 1.27  2001/05/08 07:05:02  hristov
37 Loop variable declared once (HP, Sun)
38
39 Revision 1.26  2001/05/07 08:03:22  cblume
40 Generate also hits in the amplification region
41
42 Revision 1.25  2001/03/13 09:30:35  cblume
43 Update of digitization. Moved digit branch definition to AliTRD
44
45 Revision 1.24  2001/01/26 19:56:49  hristov
46 Major upgrade of AliRoot code
47
48 Revision 1.23  2000/11/01 14:53:20  cblume
49 Merge with TRD-develop
50
51 Revision 1.17.2.6  2000/10/15 23:29:08  cblume
52 Introduced more detailed geometry for the display
53
54 Revision 1.17.2.5  2000/10/06 16:49:46  cblume
55 Made Getters const
56
57 Revision 1.17.2.4  2000/10/04 16:34:57  cblume
58 Replace include files by forward declarations
59
60 Revision 1.17.2.3  2000/09/22 14:45:17  cblume
61 Included changes for the tracking
62
63 Revision 1.17.2.2  2000/09/18 13:25:13  cblume
64 Included LoadPoints() method to display the TR photons
65
66 Revision 1.22  2000/10/02 21:28:19  fca
67 Removal of useless dependecies via forward declarations
68
69 Revision 1.21  2000/06/09 11:10:07  cblume
70 Compiler warnings and coding conventions, next round
71
72 Revision 1.20  2000/06/08 18:32:57  cblume
73 Make code compliant to coding conventions
74
75 Revision 1.19  2000/06/07 16:25:37  cblume
76 Try to remove compiler warnings on Sun and HP
77
78 Revision 1.18  2000/05/08 16:17:27  cblume
79 Merge TRD-develop
80
81 Revision 1.21  2000/06/09 11:10:07  cblume
82 Compiler warnings and coding conventions, next round
83
84 Revision 1.20  2000/06/08 18:32:57  cblume
85 Make code compliant to coding conventions
86
87 Revision 1.19  2000/06/07 16:25:37  cblume
88 Try to remove compiler warnings on Sun and HP
89
90 Revision 1.18  2000/05/08 16:17:27  cblume
91 Merge TRD-develop
92
93 Revision 1.17.2.1  2000/05/08 14:28:59  cblume
94 Introduced SetPHOShole() and SetRICHhole(). AliTRDrecPoint container is now a TObjArray
95
96 Revision 1.17  2000/02/28 19:10:26  cblume
97 Include the new TRD classes
98
99 Revision 1.16.2.2  2000/02/28 17:53:24  cblume
100 Introduce TRD geometry classes
101
102 Revision 1.16.2.1  2000/02/28 17:04:19  cblume
103 Include functions and data members for AliTRDrecPoint
104
105 Revision 1.16  2000/01/19 17:17:35  fca
106 Introducing a list of lists of hits -- more hits allowed for detector now
107
108 Revision 1.15  1999/11/02 17:04:25  fca
109 Small syntax change for HP compiler
110
111 Revision 1.14  1999/11/02 16:57:02  fca
112 Avoid non ansi warnings on HP compilers
113
114 Revision 1.13  1999/11/02 16:35:56  fca
115 New version of TRD introduced
116
117 Revision 1.12  1999/11/01 20:41:51  fca
118 Added protections against using the wrong version of FRAME
119
120 Revision 1.11  1999/09/29 09:24:34  fca
121 Introduction of the Copyright and cvs Log
122
123 */
124
125 ///////////////////////////////////////////////////////////////////////////////
126 //                                                                           //
127 //  Transition Radiation Detector                                            //
128 //  This class contains the basic functions for the Transition Radiation     //
129 //  Detector.                                                                //
130 //                                                                           //
131 ///////////////////////////////////////////////////////////////////////////////
132
133 #include <stdlib.h>
134 #include <iostream.h>
135
136 #include <TMath.h>
137 #include <TNode.h>
138 #include <TGeometry.h>
139 #include <TTree.h>                                                              
140 #include <TPGON.h> 
141 #include <TFile.h>
142 #include <TROOT.h>
143
144 #include "AliRun.h"
145 #include "AliConst.h"
146 #include "AliDigit.h"
147 #include "AliMagF.h"
148 #include "AliMC.h"                                                              
149  
150 #include "AliTRD.h"
151 #include "AliTRDhit.h"
152 #include "AliTRDpoints.h"
153 #include "AliTRDdigit.h"
154 #include "AliTRDdigitizer.h"
155 #include "AliTRDclusterizer.h"
156 #include "AliTRDgeometryHole.h"
157 #include "AliTRDgeometryFull.h"
158 #include "AliTRDrecPoint.h"
159 #include "AliTRDcluster.h"
160 #include "AliTRDdigitsManager.h"
161 #include "AliTRDdigitsMerger.h"
162
163 ClassImp(AliTRD)
164  
165 //_____________________________________________________________________________
166 AliTRD::AliTRD()
167 {
168   //
169   // Default constructor
170   //
171
172   fIshunt        = 0;
173   fGasMix        = 0;
174   fHits          = 0;
175   fDigits        = 0;
176
177   fRecPoints     = 0;
178   fNRecPoints    = 0;
179
180   fGeometry      = 0;
181
182   fGasDensity    = 0;
183   fFoilDensity   = 0;
184
185   fDrawTR        = 0;
186   fDisplayType   = 0; 
187
188 }
189  
190 //_____________________________________________________________________________
191 AliTRD::AliTRD(const char *name, const char *title)
192        : AliDetector(name,title)
193 {
194   //
195   // Standard constructor for the TRD
196   //
197
198   // Check that FRAME is there otherwise we have no place where to
199   // put TRD
200   AliModule* frame = gAlice->GetModule("FRAME");
201   if (!frame) {
202     Error("Ctor","TRD needs FRAME to be present\n");
203     exit(1);
204   } 
205
206   // Define the TRD geometry according to the FRAME geometry
207   if      (frame->IsVersion() == 0) {
208     // Geometry with hole
209     fGeometry = new AliTRDgeometryHole();
210   }
211   else if (frame->IsVersion() == 1) {
212     // Geometry without hole
213     fGeometry = new AliTRDgeometryFull();
214   }
215   else {
216     Error("Ctor","Could not find valid FRAME version\n");
217     exit(1);
218   }
219
220   // Allocate the hit array
221   fHits           = new TClonesArray("AliTRDhit"     ,405);
222   gAlice->AddHitList(fHits);
223
224   // Allocate the digits array
225   fDigits         = 0;
226
227   // Allocate the rec point array
228   fRecPoints     = new TObjArray(400);
229   fNRecPoints    = 0;
230    
231   fIshunt        = 0;
232   fGasMix        = 1;
233
234   fGasDensity    = 0;
235   fFoilDensity   = 0;
236
237   fDrawTR        = 0;
238   fDisplayType   = 0;
239
240   SetMarkerColor(kWhite);   
241
242 }
243
244 //_____________________________________________________________________________
245 AliTRD::AliTRD(const AliTRD &trd)
246 {
247   //
248   // Copy constructor
249   //
250
251   ((AliTRD &) trd).Copy(*this);
252
253 }
254
255 //_____________________________________________________________________________
256 AliTRD::~AliTRD()
257 {
258   //
259   // TRD destructor
260   //
261
262   fIshunt = 0;
263
264   if (fGeometry) {
265     delete fGeometry;
266   }
267   if (fHits) {
268     delete fHits;
269   }
270   if (fRecPoints) {
271     delete fRecPoints;
272   }
273
274 }
275
276 //_____________________________________________________________________________
277 void AliTRD::AddCluster(Float_t *pos, Int_t *digits, Int_t det, Float_t amp
278                        , Int_t *tracks, Float_t sigmaY2, Int_t iType)
279 {
280   //
281   // Add a cluster for the TRD
282   //
283
284   Int_t   pla = fGeometry->GetPlane(det);
285   Int_t   cha = fGeometry->GetChamber(det);
286   Int_t   sec = fGeometry->GetSector(det);
287
288   Float_t padRow  = pos[0]; 
289   Float_t padCol  = pos[1]; 
290   Float_t row0    = fGeometry->GetRow0(pla,cha,sec);
291   Float_t col0    = fGeometry->GetCol0(pla);
292   Float_t rowSize = fGeometry->GetRowPadSize(pla,cha,sec);
293   Float_t colSize = fGeometry->GetColPadSize(pla);
294
295   AliTRDcluster *c = new AliTRDcluster();
296
297   c->SetDetector(det);
298   c->AddTrackIndex(tracks);
299
300   c->SetQ(amp);
301
302   c->SetLocalTimeBin(((Int_t) pos[2]));
303   c->SetY(- (col0 + padCol * colSize));
304   c->SetZ(   row0 + padRow * rowSize);
305   
306   c->SetSigmaY2((sigmaY2 + 1./12.) * colSize*colSize);   
307   c->SetSigmaZ2(rowSize * rowSize / 12.);
308
309   switch (iType) {
310   case 0:
311     c->Set2pad();
312     break;
313   case 1:
314     c->Set3pad();
315     break;
316   case 2:
317     c->Set4pad();
318     break;
319   case 3:
320     c->Set5pad();
321     break;
322   case 4:
323     c->SetLarge();
324     break;
325   };
326
327   fRecPoints->Add(c);
328
329 }
330
331 //_____________________________________________________________________________
332 void AliTRD::Hits2Digits()
333 {
334   //
335   // Create digits
336   //
337
338   AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
339                                                   ,"TRD digitizer class");
340   digitizer->SetVerbose(GetDebug());
341
342   // Set the parameter
343   digitizer->SetDiffusion();
344   digitizer->SetExB();
345   digitizer->SetEvent(gAlice->GetEvNumber());
346
347   // Initialization
348   digitizer->InitDetector();
349     
350   // Create the digits
351   digitizer->MakeDigits();
352   
353   // Write the digits into the input file
354   if (digitizer->MakeBranch(fDigitsFile)) {
355
356     digitizer->WriteDigits();
357
358     // Save the digitizer class in the AliROOT 
359     digitizer->Write();
360
361   }
362
363 }
364
365 //_____________________________________________________________________________
366 void AliTRD::Hits2SDigits()
367 {
368   //
369   // Create summable digits
370   //
371
372   AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
373                                                   ,"TRD digitizer class");
374   digitizer->SetVerbose(GetDebug());
375
376   // For the summable digits
377   digitizer->SetSDigits(kTRUE);
378
379   // Set the parameter
380   digitizer->SetDiffusion();
381   digitizer->SetExB();
382   digitizer->SetEvent(gAlice->GetEvNumber());
383
384   // Initialization
385   digitizer->InitDetector();
386     
387   // Create the TRD s-digits branch
388   digitizer->MakeDigits();
389   
390   // Write the digits into the input file
391   if (digitizer->MakeBranch(fDigitsFile)) {
392
393     digitizer->WriteDigits();
394
395     // Save the digitizer class in the AliROOT 
396     digitizer->Write();
397
398   }
399
400 }
401
402 //_____________________________________________________________________________
403 void AliTRD::SDigits2Digits()
404 {
405   //
406   // Create final digits from summable digits
407   //
408
409    // Create the TRD digitizer
410   AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
411                                                   ,"TRD digitizer class");  
412   digitizer->SetVerbose(GetDebug());
413
414   // Set the parameter
415   digitizer->SetEvent(gAlice->GetEvNumber());
416
417   // Initialization
418   digitizer->InitDetector();
419
420   // Read the s-digits via digits manager
421   AliTRDdigitsManager *sdigitsManager = new AliTRDdigitsManager();
422   sdigitsManager->SetVerbose(GetDebug());
423   sdigitsManager->SetSDigits(kTRUE);
424   if (fDigitsFile) {
425     sdigitsManager->Open(fDigitsFile);
426   }
427   sdigitsManager->CreateArrays();
428   sdigitsManager->ReadDigits();
429
430   // Add the s-digits to the input list 
431   digitizer->AddSDigitsManager(sdigitsManager);
432
433   // Convert the s-digits to normal digits
434   digitizer->SDigits2Digits();
435
436   // Store the digits
437   if (digitizer->MakeBranch(fDigitsFile)) {
438
439     digitizer->WriteDigits();
440
441   }
442
443 }
444
445 //_____________________________________________________________________________
446 void AliTRD::AddHit(Int_t track, Int_t det, Float_t *hits, Int_t q)
447 {
448   //
449   // Add a hit for the TRD
450   //
451
452   TClonesArray &lhits = *fHits;
453   new(lhits[fNhits++]) AliTRDhit(fIshunt,track,det,hits,q);
454
455 }
456
457 //_____________________________________________________________________________
458 void AliTRD::BuildGeometry()
459 {
460   //
461   // Create the ROOT TNode geometry for the TRD
462   //
463
464   TNode *node, *top;
465   TPGON *pgon;
466
467   Float_t rmin, rmax;
468   Float_t zmax1, zmax2;
469
470   Int_t   iPlan;
471  
472   const Int_t kColorTRD = 46;
473   
474   // Find the top node alice
475   top = gAlice->GetGeometry()->GetNode("alice");
476   
477   if      (fDisplayType == 0) {
478
479     pgon = new TPGON("S_TRD","TRD","void",0,360,AliTRDgeometry::Nsect(),4);
480     rmin = AliTRDgeometry::Rmin();
481     rmax = AliTRDgeometry::Rmax();
482     pgon->DefineSection(0,-AliTRDgeometry::Zmax1(),rmax,rmax);
483     pgon->DefineSection(1,-AliTRDgeometry::Zmax2(),rmin,rmax);
484     pgon->DefineSection(2, AliTRDgeometry::Zmax2(),rmin,rmax);
485     pgon->DefineSection(3, AliTRDgeometry::Zmax1(),rmax,rmax);
486     top->cd();
487     node = new TNode("TRD","TRD","S_TRD",0,0,0,"");
488     node->SetLineColor(kColorTRD);
489     fNodes->Add(node);
490
491   }
492   else if (fDisplayType == 1) {
493
494     Char_t name[7];
495
496     Float_t slope = (AliTRDgeometry::Zmax1() - AliTRDgeometry::Zmax2())
497                   / (AliTRDgeometry::Rmax()  - AliTRDgeometry::Rmin());
498
499     rmin  = AliTRDgeometry::Rmin() + AliTRDgeometry::Ccframe()/2.
500                                    + AliTRDgeometry::DrZpos()
501                                    - AliTRDgeometry::DrThick()/2.;
502     rmax  = rmin + AliTRDgeometry::DrThick();
503
504     Float_t thickness = rmin - AliTRDgeometry::Rmin();
505     zmax2 = AliTRDgeometry::Zmax2() + slope * thickness;
506     zmax1 = zmax2 + slope * AliTRDgeometry::DrThick();
507
508     for (iPlan = 0; iPlan < AliTRDgeometry::Nplan(); iPlan++) {
509
510       sprintf(name,"S_TR1%d",iPlan);
511       pgon  = new TPGON(name,"TRD","void",0,360,AliTRDgeometry::Nsect(),4);
512       pgon->DefineSection(0,-zmax1,rmax,rmax);
513       pgon->DefineSection(1,-zmax2,rmin,rmax);
514       pgon->DefineSection(2, zmax2,rmin,rmax);
515       pgon->DefineSection(3, zmax1,rmax,rmax);
516       top->cd();
517       node = new TNode("TRD","TRD",name,0,0,0,"");
518       node->SetLineColor(kColorTRD);
519       fNodes->Add(node);
520
521       Float_t height = AliTRDgeometry::Cheight() + AliTRDgeometry::Cspace(); 
522       rmin  = rmin  + height;
523       rmax  = rmax  + height;
524       zmax1 = zmax1 + slope * height;
525       zmax2 = zmax2 + slope * height;
526
527     }
528
529     thickness += AliTRDgeometry::DrThick();
530     rmin  = AliTRDgeometry::Rmin() + thickness;
531     rmax  = rmin + AliTRDgeometry::AmThick();
532     zmax2 = AliTRDgeometry::Zmax2() + slope * thickness;
533     zmax1 = zmax2 + slope * AliTRDgeometry::AmThick();
534
535     for (iPlan = 0; iPlan < AliTRDgeometry::Nplan(); iPlan++) {
536
537       sprintf(name,"S_TR2%d",iPlan);
538       pgon  = new TPGON(name,"TRD","void",0,360,AliTRDgeometry::Nsect(),4);
539       pgon->DefineSection(0,-zmax1,rmax,rmax);
540       pgon->DefineSection(1,-zmax2,rmin,rmax);
541       pgon->DefineSection(2, zmax2,rmin,rmax);
542       pgon->DefineSection(3, zmax1,rmax,rmax);
543       top->cd();
544       node = new TNode("TRD","TRD",name,0,0,0,"");
545       node->SetLineColor(kColorTRD);
546       fNodes->Add(node);
547
548       Float_t height = AliTRDgeometry::Cheight() + AliTRDgeometry::Cspace(); 
549       rmin  = rmin  + height;
550       rmax  = rmax  + height;
551       zmax1 = zmax1 + slope * height;
552       zmax2 = zmax2 + slope * height;
553
554     }
555
556   }
557
558 }
559  
560 //_____________________________________________________________________________
561 void AliTRD::Copy(TObject &trd)
562 {
563   //
564   // Copy function
565   //
566
567   ((AliTRD &) trd).fGasMix      = fGasMix;
568   ((AliTRD &) trd).fGeometry    = fGeometry;       
569   ((AliTRD &) trd).fRecPoints   = fRecPoints;
570   ((AliTRD &) trd).fNRecPoints  = fNRecPoints;
571   ((AliTRD &) trd).fGasDensity  = fGasDensity;
572   ((AliTRD &) trd).fFoilDensity = fFoilDensity;
573   ((AliTRD &) trd).fDrawTR      = fDrawTR;
574   ((AliTRD &) trd).fDisplayType = fDisplayType;
575
576   //AliDetector::Copy(trd);
577
578 }
579
580 //_____________________________________________________________________________
581 void AliTRD::CreateGeometry()
582 {
583   //
584   // Creates the volumes for the TRD chambers
585   //
586
587   // Check that FRAME is there otherwise we have no place where to put the TRD
588   AliModule* frame = gAlice->GetModule("FRAME");
589   if (!frame) {
590     printf(" The TRD needs the FRAME to be defined first\n");
591     return;
592   }
593
594   fGeometry->CreateGeometry(fIdtmed->GetArray() - 1299);
595
596 }
597  
598 //_____________________________________________________________________________
599 void AliTRD::CreateMaterials()
600 {
601   //
602   // Create the materials for the TRD
603   // Origin Y.Foka
604   //
605
606   Int_t   isxfld = gAlice->Field()->Integ();
607   Float_t sxmgmx = gAlice->Field()->Max();
608   
609   // For polyethilene (CH2) 
610   Float_t ape[2] = { 12., 1. };
611   Float_t zpe[2] = {  6., 1. };
612   Float_t wpe[2] = {  1., 2. };
613   Float_t dpe    = 0.95;
614
615   // For mylar (C5H4O2) 
616   Float_t amy[3] = { 12., 1., 16. };
617   Float_t zmy[3] = {  6., 1.,  8. };
618   Float_t wmy[3] = {  5., 4.,  2. };
619   Float_t dmy    = 1.39;
620
621   // For CO2 
622   Float_t aco[2] = { 12., 16. };
623   Float_t zco[2] = {  6.,  8. };
624   Float_t wco[2] = {  1.,  2. };
625   Float_t dco    = 0.001977;
626
627   // For water
628   Float_t awa[2] = {  1., 16. };
629   Float_t zwa[2] = {  1.,  8. };
630   Float_t wwa[2] = {  2.,  1. };
631   Float_t dwa    = 1.0;
632
633   // For isobutane (C4H10)
634   Float_t ais[2] = { 12.,  1. };
635   Float_t zis[2] = {  6.,  1. };
636   Float_t wis[2] = {  4., 10. };
637   Float_t dis    = 0.00267;
638
639   // For plexiglas (C5H8O2)
640   Float_t apg[3] = { 12.011 ,  1.0    , 15.9994 };
641   Float_t zpg[3] = {  6.0   ,  1.0    ,  8.0    };
642   Float_t wpg[3] = {  5.0   ,  8.0    ,  2.0    };
643   Float_t dpg    = 1.18; 
644
645   // For Xe/CO2-gas-mixture 
646   // Xe-content of the Xe/CO2-mixture (85% / 15%) 
647   Float_t fxc    = .85;
648   // Xe-content of the Xe/Isobutane-mixture (97% / 3%) 
649   Float_t fxi    = .97;
650   Float_t dxe    = .005858;
651   
652   // General tracking parameter
653   Float_t tmaxfd = -10.;
654   Float_t stemax = -1e10;
655   Float_t deemax = -0.1;
656   Float_t epsil  =  1e-4;
657   Float_t stmin  = -0.001;
658   
659   Float_t absl, radl, d, buf[1];
660   Float_t agm[2], zgm[2], wgm[2];
661   Float_t dgm1, dgm2;
662   Int_t   nbuf;
663   
664   //////////////////////////////////////////////////////////////////////////
665   //     Define Materials 
666   //////////////////////////////////////////////////////////////////////////
667
668   AliMaterial( 1, "Al"   ,  26.98, 13.0, 2.7     ,     8.9 ,    37.2);
669   AliMaterial( 2, "Air"  ,  14.61,  7.3, 0.001205, 30420.0 , 67500.0);
670   AliMaterial( 4, "Xe"   , 131.29, 54.0, dxe     ,  1447.59,     0.0);
671   AliMaterial( 5, "Cu"   ,  63.54, 29.0, 8.96    ,     1.43,    14.8);
672   AliMaterial( 6, "C"    ,  12.01,  6.0, 2.265   ,    18.8 ,    74.4);
673   AliMaterial(12, "G10"  ,  20.00, 10.0, 1.7     ,    19.4 ,   999.0);
674   AliMaterial(15, "Sn"   , 118.71, 50.0, 7.31    ,     1.21,    14.8);
675   AliMaterial(16, "Si"   ,  28.09, 14.0, 2.33    ,     9.36,    37.2);
676   AliMaterial(17, "Epoxy",  17.75,  8.9, 1.8     ,    21.82,   999.0);
677
678   // Mixtures 
679   AliMixture(3, "Polyethilene",   ape, zpe, dpe, -2, wpe);
680   AliMixture(7, "Mylar",          amy, zmy, dmy, -3, wmy);
681   AliMixture(8, "CO2",            aco, zco, dco, -2, wco);
682   AliMixture(9, "Isobutane",      ais, zis, dis, -2, wis);
683   AliMixture(13,"Water",          awa, zwa, dwa, -2, wwa);
684   AliMixture(14,"Plexiglas",      apg, zpg, dpg, -3, wpg);
685
686   // Gas mixtures
687   Char_t namate[21];
688   // Xe/CO2-mixture
689   // Get properties of Xe 
690   gMC->Gfmate((*fIdmate)[4], namate, agm[0], zgm[0], d, radl, absl, buf, nbuf);
691   // Get properties of CO2 
692   gMC->Gfmate((*fIdmate)[8], namate, agm[1], zgm[1], d, radl, absl, buf, nbuf);
693   // Create gas mixture 
694   wgm[0] = fxc;
695   wgm[1] = 1. - fxc;
696   dgm1   = wgm[0] * dxe + wgm[1] * dco;
697   AliMixture(10, "Gas mixture 1", agm, zgm, dgm1,  2, wgm);
698   // Xe/Isobutane-mixture
699   // Get properties of Xe 
700   gMC->Gfmate((*fIdmate)[4], namate, agm[0], zgm[0], d, radl, absl, buf, nbuf);
701   // Get properties of Isobutane
702   gMC->Gfmate((*fIdmate)[9], namate, agm[1], zgm[1], d, radl, absl, buf, nbuf);
703   // Create gas mixture 
704   wgm[0] = fxi;
705   wgm[1] = 1. - fxi;
706   dgm2   = wgm[0] * dxe + wgm[1] * dis;
707   AliMixture(11, "Gas mixture 2", agm, zgm, dgm2,  2, wgm);
708  
709   //////////////////////////////////////////////////////////////////////////
710   //     Tracking Media Parameters 
711   //////////////////////////////////////////////////////////////////////////
712
713   // Al Frame 
714   AliMedium(1, "Al Frame",   1, 0, isxfld, sxmgmx
715                 , tmaxfd, stemax, deemax, epsil, stmin);
716   // Air 
717   AliMedium(2, "Air",        2, 0, isxfld, sxmgmx
718                 , tmaxfd, stemax, deemax, epsil, stmin);
719   // Polyethilene 
720   AliMedium(3, "Radiator",   3, 0, isxfld, sxmgmx
721                 , tmaxfd, stemax, deemax, epsil, stmin);
722   // Xe 
723   AliMedium(4, "Xe",         4, 1, isxfld, sxmgmx
724                 , tmaxfd, stemax, deemax, epsil, stmin);
725   // Cu pads 
726   AliMedium(5, "Padplane",   5, 1, isxfld, sxmgmx
727                 , tmaxfd, stemax, deemax, epsil, stmin);
728   // Fee + cables 
729   AliMedium(6, "Readout",    1, 0, isxfld, sxmgmx
730                 , tmaxfd, stemax, deemax, epsil, stmin);
731   // C frame 
732   AliMedium(7, "C Frame",    6, 0, isxfld, sxmgmx
733                 , tmaxfd, stemax, deemax, epsil, stmin);
734   // Mylar foils 
735   AliMedium(8, "Mylar",      7, 0, isxfld, sxmgmx
736                 , tmaxfd, stemax, deemax, epsil, stmin);
737   if (fGasMix == 1) {
738     // Gas-mixture (Xe/CO2) 
739     AliMedium(9, "Gas-mix",   10, 1, isxfld, sxmgmx
740                   , tmaxfd, stemax, deemax, epsil, stmin);
741   }
742   else {
743     // Gas-mixture (Xe/Isobutane) 
744     AliMedium(9, "Gas-mix",   11, 1, isxfld, sxmgmx
745                   , tmaxfd, stemax, deemax, epsil, stmin);
746   }
747   // Nomex-honeycomb (use carbon for the time being) 
748   AliMedium(10, "Nomex",      6, 0, isxfld, sxmgmx
749                 , tmaxfd, stemax, deemax, epsil, stmin);
750   // Kapton foils (use Mylar for the time being) 
751   AliMedium(11, "Kapton",     7, 0, isxfld, sxmgmx
752                 , tmaxfd, stemax, deemax, epsil, stmin);
753   // Gas-filling of the radiator 
754   AliMedium(12, "CO2",        8, 0, isxfld, sxmgmx
755                 , tmaxfd, stemax, deemax, epsil, stmin);
756   // G10-plates
757   AliMedium(13, "G10-plates",12, 0, isxfld, sxmgmx
758                 , tmaxfd, stemax, deemax, epsil, stmin);
759   // Cooling water
760   AliMedium(14, "Water",     13, 0, isxfld, sxmgmx
761                 , tmaxfd, stemax, deemax, epsil, stmin);
762   // Rohacell (plexiglas) for the radiator
763   AliMedium(15, "Rohacell",  14, 0, isxfld, sxmgmx
764                 , tmaxfd, stemax, deemax, epsil, stmin);
765   // Al layer in MCMs
766   AliMedium(16, "MCM-Al"  ,   1, 0, isxfld, sxmgmx
767                 , tmaxfd, stemax, deemax, epsil, stmin);
768   // Sn layer in MCMs
769   AliMedium(17, "MCM-Sn"  ,  15, 0, isxfld, sxmgmx
770                 , tmaxfd, stemax, deemax, epsil, stmin);
771   // Cu layer in MCMs
772   AliMedium(18, "MCM-Cu"  ,   5, 0, isxfld, sxmgmx
773                 , tmaxfd, stemax, deemax, epsil, stmin);
774   // G10 layer in MCMs
775   AliMedium(19, "MCM-G10" ,  12, 0, isxfld, sxmgmx
776                 , tmaxfd, stemax, deemax, epsil, stmin);
777   // Si in readout chips
778   AliMedium(20, "Chip-Si" ,  16, 0, isxfld, sxmgmx
779                 , tmaxfd, stemax, deemax, epsil, stmin);
780   // Epoxy in readout chips
781   AliMedium(21, "Chip-Ep" ,  17, 0, isxfld, sxmgmx
782                 , tmaxfd, stemax, deemax, epsil, stmin);
783   // PE in connectors
784   AliMedium(22, "Conn-PE" ,   3, 0, isxfld, sxmgmx
785                 , tmaxfd, stemax, deemax, epsil, stmin);
786   // Cu in connectors
787   AliMedium(23, "Chip-Cu" ,   5, 0, isxfld, sxmgmx
788                 , tmaxfd, stemax, deemax, epsil, stmin);
789   // Al of cooling pipes
790   AliMedium(24, "Cooling" ,   1, 0, isxfld, sxmgmx
791                 , tmaxfd, stemax, deemax, epsil, stmin);
792
793   // Save the density values for the TRD absorbtion
794   fFoilDensity = dmy;
795   if (fGasMix == 1)
796     fGasDensity = dgm1;
797   else
798     fGasDensity = dgm2;
799
800 }
801
802 //_____________________________________________________________________________
803 void AliTRD::DrawModule()
804 {
805   //
806   // Draw a shaded view of the Transition Radiation Detector version 0
807   //
808
809   // Set everything unseen
810   gMC->Gsatt("*"   ,"SEEN",-1);
811   
812   // Set ALIC mother transparent
813   gMC->Gsatt("ALIC","SEEN", 0);
814   
815   // Set the volumes visible
816   if (fGeometry->IsVersion() == 0) {
817     gMC->Gsatt("B071","SEEN", 0);
818     gMC->Gsatt("B074","SEEN", 0);
819     gMC->Gsatt("B075","SEEN", 0);
820     gMC->Gsatt("B077","SEEN", 0);
821     gMC->Gsatt("BTR1","SEEN", 0);
822     gMC->Gsatt("BTR2","SEEN", 0);
823     gMC->Gsatt("BTR3","SEEN", 0);
824     gMC->Gsatt("TRD1","SEEN", 0);
825     gMC->Gsatt("TRD2","SEEN", 0);
826     gMC->Gsatt("TRD3","SEEN", 0);
827   }
828   else {
829     gMC->Gsatt("B071","SEEN", 0);
830     gMC->Gsatt("B074","SEEN", 0);
831     gMC->Gsatt("B075","SEEN", 0);
832     gMC->Gsatt("B077","SEEN", 0);
833     gMC->Gsatt("BTR1","SEEN", 0);
834     gMC->Gsatt("BTR2","SEEN", 0);
835     gMC->Gsatt("BTR3","SEEN", 0);
836     gMC->Gsatt("TRD1","SEEN", 0);
837     if (fGeometry->GetPHOShole())
838       gMC->Gsatt("TRD2","SEEN", 0);
839     if (fGeometry->GetRICHhole())
840       gMC->Gsatt("TRD3","SEEN", 0);
841   }
842   gMC->Gsatt("UCII","SEEN", 0);
843   gMC->Gsatt("UCIM","SEEN", 0);
844   gMC->Gsatt("UCIO","SEEN", 0);
845   gMC->Gsatt("UL02","SEEN", 1);
846   gMC->Gsatt("UL05","SEEN", 1);
847   gMC->Gsatt("UL06","SEEN", 1);
848   
849   gMC->Gdopt("hide", "on");
850   gMC->Gdopt("shad", "on");
851   gMC->Gsatt("*", "fill", 7);
852   gMC->SetClipBox(".");
853   gMC->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
854   gMC->DefaultRange();
855   gMC->Gdraw("alic", 40, 30, 0, 12, 9.4, .021, .021);
856   gMC->Gdhead(1111, "Transition Radiation Detector");
857   gMC->Gdman(18, 4, "MAN");
858
859 }
860
861 //_____________________________________________________________________________
862 Int_t AliTRD::DistancetoPrimitive(Int_t , Int_t )
863 {
864   //
865   // Distance between the mouse and the TRD detector on the screen
866   // Dummy routine
867   
868   return 9999;
869
870 }
871  
872 //_____________________________________________________________________________
873 void AliTRD::Init()
874 {
875   //
876   // Initialize the TRD detector after the geometry has been created
877   //
878
879   Int_t i;
880
881   if (fDebug) {
882     printf("\n%s: ",ClassName());
883     for (i = 0; i < 35; i++) printf("*");
884     printf(" TRD_INIT ");
885     for (i = 0; i < 35; i++) printf("*");
886     printf("\n");
887   }
888
889   if      (fGeometry->IsVersion() == 0) {
890     printf("%s: Geometry for spaceframe with holes initialized\n",ClassName());
891   }
892   else if (fGeometry->IsVersion() == 1) {
893     printf("%s: Geometry for spaceframe without holes initialized\n",ClassName());
894     if (fGeometry->GetPHOShole())
895       printf("%s: Leave space in front of PHOS free\n",ClassName());
896     if (fGeometry->GetRICHhole())
897       printf("%s: Leave space in front of RICH free\n",ClassName());
898   }
899   
900   if (fGasMix == 1) {
901     printf("%s: Gas Mixture: 85%% Xe + 15%% CO2\n",ClassName());
902   }
903   else {
904     printf("%s: Gas Mixture: 97%% Xe + 3%% Isobutane\n",ClassName());
905   }
906
907 }
908
909 //_____________________________________________________________________________
910 void AliTRD::LoadPoints(Int_t track)
911 {
912   //
913   // Store x, y, z of all hits in memory.
914   // Hit originating from TR photons are given a different color
915   //
916
917   if (!fDrawTR) {
918     AliDetector::LoadPoints(track);
919     return;
920   }
921
922   if (fHits == 0) return;
923
924   Int_t nhits  = fHits->GetEntriesFast();
925   if (nhits == 0) return;
926
927   Int_t tracks = gAlice->GetNtrack();
928   if (fPoints == 0) fPoints = new TObjArray(tracks);
929
930   AliTRDhit *ahit;
931   
932   Int_t    *ntrkE = new Int_t[tracks];
933   Int_t    *ntrkT = new Int_t[tracks];
934   Int_t    *limiE = new Int_t[tracks];
935   Int_t    *limiT = new Int_t[tracks];
936   Float_t **coorE = new Float_t*[tracks];
937   Float_t **coorT = new Float_t*[tracks];
938   for(Int_t i = 0; i < tracks; i++) {
939     ntrkE[i] = 0;
940     ntrkT[i] = 0;
941     coorE[i] = 0;
942     coorT[i] = 0;
943     limiE[i] = 0;
944     limiT[i] = 0;
945   }
946   
947   AliTRDpoints  *points = 0;
948   Float_t       *fp     = 0;
949   Int_t          trk;
950   Int_t          chunk  = nhits / 4 + 1;
951
952   // Loop over all the hits and store their position
953   for (Int_t hit = 0; hit < nhits; hit++) {
954
955     ahit = (AliTRDhit *) fHits->UncheckedAt(hit);
956
957     // dEdx hits
958     if (ahit->FromDrift() || ahit->FromAmplification()) {
959
960       trk = ahit->GetTrack();
961       if (ntrkE[trk] == limiE[trk]) {
962         // Initialise a new track
963         fp = new Float_t[3*(limiE[trk]+chunk)];
964         if (coorE[trk]) {
965           memcpy(fp,coorE[trk],sizeof(Float_t)*3*limiE[trk]);
966           delete [] coorE[trk];
967         }
968         limiE[trk] += chunk;
969         coorE[trk]  = fp;
970       } 
971       else {
972         fp = coorE[trk];
973       }
974       fp[3*ntrkE[trk]  ] = ahit->X();
975       fp[3*ntrkE[trk]+1] = ahit->Y();
976       fp[3*ntrkE[trk]+2] = ahit->Z();
977       ntrkE[trk]++;
978
979     }
980     // TR photon hits
981     else if (ahit->FromTRphoton()) {
982
983       trk = ahit->GetTrack();
984       if (ntrkT[trk] == limiT[trk]) {
985         // Initialise a new track
986         fp = new Float_t[3*(limiT[trk]+chunk)];
987         if (coorT[trk]) {
988           memcpy(fp,coorT[trk],sizeof(Float_t)*3*limiT[trk]);
989           delete [] coorT[trk];
990         }
991         limiT[trk] += chunk;
992         coorT[trk]  = fp;
993       } 
994       else {
995         fp = coorT[trk];
996       }
997       fp[3*ntrkT[trk]  ] = ahit->X();
998       fp[3*ntrkT[trk]+1] = ahit->Y();
999       fp[3*ntrkT[trk]+2] = ahit->Z();
1000       ntrkT[trk]++;
1001
1002     }
1003
1004   }
1005
1006   for (trk = 0; trk < tracks; ++trk) {
1007
1008     if (ntrkE[trk] || ntrkT[trk]) {
1009
1010       points = new AliTRDpoints();
1011       points->SetDetector(this);
1012       points->SetParticle(trk);
1013
1014       // Set the dEdx points
1015       if (ntrkE[trk]) {
1016         points->SetMarkerColor(GetMarkerColor());
1017         points->SetMarkerSize(GetMarkerSize());
1018         points->SetPolyMarker(ntrkE[trk],coorE[trk],GetMarkerStyle());
1019         delete [] coorE[trk];
1020         coorE[trk] = 0;
1021       }
1022
1023       // Set the TR photon points
1024       if (ntrkT[trk]) {
1025         points->SetTRpoints(ntrkT[trk],coorT[trk]);
1026         delete [] coorT[trk];
1027         coorT[trk] = 0;
1028       }
1029
1030       fPoints->AddAt(points,trk);
1031
1032     }
1033
1034   }
1035
1036   delete [] coorE;
1037   delete [] coorT;
1038   delete [] ntrkE;
1039   delete [] ntrkT;
1040   delete [] limiE;
1041   delete [] limiT;
1042
1043 }
1044
1045 //_____________________________________________________________________________
1046 void AliTRD::MakeBranch(Option_t* option, const char *file)
1047 {
1048   //
1049   // Create Tree branches for the TRD digits.
1050   //
1051
1052   Int_t  buffersize = 4000;
1053   Char_t branchname[15];
1054   sprintf(branchname,"%s",GetName());
1055
1056   const char *cD = strstr(option,"D");
1057
1058   AliDetector::MakeBranch(option,file);
1059
1060   if (fDigits && gAlice->TreeD() && cD) {
1061     MakeBranchInTree(gAlice->TreeD(),branchname,&fDigits,buffersize,file);
1062   }     
1063
1064 }
1065
1066 //_____________________________________________________________________________
1067 void AliTRD::ResetDigits()
1068 {
1069   //
1070   // Reset number of digits and the digits array for this detector
1071   //
1072
1073   fNdigits = 0;
1074   if (fDigits) fDigits->Clear();
1075
1076 }
1077
1078 //_____________________________________________________________________________
1079 void AliTRD::ResetRecPoints()
1080 {
1081   //
1082   // Reset number of reconstructed points and the point array
1083   //
1084
1085   if (fRecPoints) {
1086     fNRecPoints = 0;
1087     Int_t nentr = fRecPoints->GetEntriesFast();
1088     for (Int_t i = 0; i < nentr; i++) delete fRecPoints->RemoveAt(i);
1089   }
1090
1091 }
1092
1093 //_____________________________________________________________________________
1094 void AliTRD::SetTreeAddress()
1095 {
1096   //
1097   // Set the branch addresses for the trees.
1098   //
1099
1100   Char_t branchname[15];
1101
1102   AliDetector::SetTreeAddress();
1103
1104   TBranch *branch;
1105   TTree   *treeR = gAlice->TreeR();
1106
1107   if (treeR) {
1108     sprintf(branchname,"%scluster",GetName());
1109     if (fRecPoints) {
1110       branch = treeR->GetBranch(branchname);
1111       if (branch) {
1112         branch->SetAddress(&fRecPoints);
1113       }
1114     }
1115   }
1116
1117 }
1118
1119 //_____________________________________________________________________________
1120 void AliTRD::SetGasMix(Int_t imix)
1121 {
1122   //
1123   // Defines the gas mixture (imix=0:  Xe/Isobutane imix=1: Xe/CO2)
1124   //
1125   
1126   if ((imix < 0) || (imix > 1)) {
1127     printf("Wrong input value: %d\n",imix);
1128     printf("Use standard setting\n");
1129     fGasMix = 1;
1130     return;
1131   }
1132
1133   fGasMix = imix;
1134
1135 }
1136
1137 //_____________________________________________________________________________
1138 void AliTRD::SetPHOShole()
1139 {
1140   //
1141   // Selects a geometry with a hole in front of the PHOS
1142   //
1143
1144   fGeometry->SetPHOShole();
1145
1146 }
1147
1148 //_____________________________________________________________________________
1149 void AliTRD::SetRICHhole()
1150 {
1151   //
1152   // Selects a geometry with a hole in front of the RICH
1153   //
1154
1155   fGeometry->SetRICHhole();
1156
1157 }
1158
1159 //_____________________________________________________________________________
1160 AliTRD &AliTRD::operator=(const AliTRD &trd)
1161 {
1162   //
1163   // Assignment operator
1164   //
1165
1166   if (this != &trd) ((AliTRD &) trd).Copy(*this);
1167   return *this;
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
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