868564393fa614c0fe0974ee863c99d3141edddf
[u/mrichter/AliRoot.git] / TRD / AliTRDv2.cxx
1 ///////////////////////////////////////////////////////////////////////////////
2 //                                                                           //
3 //  Transition Radiation Detector version 2 -- detailed simulation           //
4 //                                                                           //
5 //Begin_Html
6 /*
7 <img src="gif/AliTRDv2Class.gif">
8 */
9 //End_Html
10 //                                                                           //
11 //                                                                           //
12 ///////////////////////////////////////////////////////////////////////////////
13
14 #include <TMath.h>
15 #include <TVector.h>
16
17 #include "AliTRDv2.h"
18 #include "AliRun.h"
19 #include "AliMC.h"
20 #include "AliConst.h"
21
22 ClassImp(AliTRDv2)
23
24 //_____________________________________________________________________________
25 AliTRDv2::AliTRDv2(const char *name, const char *title) 
26          :AliTRD(name, title) 
27 {
28   //
29   // Standard constructor for Transition Radiation Detector version 2
30   //
31   for (Int_t icham = 0; icham < ncham; ++icham) {
32     fIdSensI[icham] = 0;
33     fIdSensN[icham] = 0;
34     fIdSensO[icham] = 0;
35   }
36   fDeltaE = NULL;
37
38   SetBufferSize(128000);
39 }
40
41 AliTRDv2::~AliTRDv2()
42 {
43    if (fDeltaE)  delete fDeltaE;
44 }
45  
46 //_____________________________________________________________________________
47 void AliTRDv2::CreateGeometry()
48 {
49   //
50   // Create geometry for the Transition Radiation Detector version 2
51   // This version covers the full azimuth. 
52   // --- Author :  Christoph Blume (GSI) 20/5/99 
53   //
54   // --- Volume names : 
55   //       TRD         --> Mother TRD volume                              (Al) 
56   //       UTRS        --> Sectors of the sub-detector                    (Al)
57   //       UTRI        --> Inner part of the detector frame               (Air) 
58   //     The chambers 
59   //       UCI1-6      --> The frame of the inner chambers                (C) 
60   //       UCN1-6      --> The frame of the neighbouring chambers         (C) 
61   //       UCO1-6      --> The frame of the outer chambers                (C) 
62   //       UII1-6      --> The inner part of the inner chambers           (Air) 
63   //       UIN1-6      --> The inner part of the neighbouring chambers    (Air) 
64   //       UIO1-6      --> The inner part of the outer chambers           (Air) 
65   //     The layers inside a chamber 
66   //       UT0I(N,O)   --> Radiator seal                                  (G10)
67   //       UT1I(N,O)   --> Radiator                                       (CO2)
68   //       UT2I(N,O)   --> Polyethylene of radiator                       (PE)
69   //       UT3I(N,O)   --> Entrance window                                (Mylar)
70   //       UXI(N,O)1-6 --> Gas volume (sensitive)                         (Xe/Isobutane)
71   //       UT5I(N,O)   --> Pad plane                                      (Cu)
72   //       UT6I(N,O)   --> Support structure                              (G10)
73   //       UT7I(N,O)   --> FEE + signal lines                             (Cu)
74   //       UT8I(N,O)   --> Polyethylene of cooling device                 (PE)
75   //       UT9I(N,O)   --> Cooling water                                  (Water)
76   //
77   //Begin_Html
78   /*
79     <img src="gif/AliTRDv2.gif">
80   */
81   //End_Html
82   //Begin_Html
83   /*
84     <img src="gif/AliTRDv2Tree.gif">
85   */
86   //End_Html
87   
88   Float_t xpos, ypos, zpos;
89   Int_t   idmat[2];
90
91   const Int_t nparmo = 10;
92   const Int_t nparfr =  4;
93   const Int_t nparch =  3;
94   const Int_t nparic =  4;
95   const Int_t nparnc =  4;
96   const Int_t nparoc = 11;
97
98   Float_t par_mo[nparmo];
99   Float_t par_fr[nparfr];
100   Float_t par_ch[nparch];
101   Float_t par_ic[nparic];
102   Float_t par_nc[nparnc];
103   Float_t par_oc[nparoc];
104
105   Int_t *idtmed = gAlice->Idtmed();
106   
107   AliMC* pMC = AliMC::GetMC();
108
109   //////////////////////////////////////////////////////////////////////////
110   //     Definition of Volumes 
111   //////////////////////////////////////////////////////////////////////////
112   
113   // Definition of the mother volume for the TRD (Al) 
114   par_mo[0] =   0.;
115   par_mo[1] = 360.;
116   par_mo[2] = nsect;
117   par_mo[3] = 2.;
118   par_mo[4] = -zmax1;
119   par_mo[5] = rmin;
120   par_mo[6] = rmax;
121   par_mo[7] =  zmax1;
122   par_mo[8] = rmin;
123   par_mo[9] = rmax;
124   pMC->Gsvolu("TRD ", "PGON", idtmed[1301-1], par_mo, nparmo);
125   pMC->Gsdvn("UTRS", "TRD ", nsect, 2);
126
127   // The minimal width of a sector in rphi-direction
128   Float_t widmi = rmin * TMath::Sin(kPI/nsect);
129   // The maximal width of a sector in rphi-direction
130   Float_t widma = rmax * TMath::Sin(kPI/nsect);
131   // The total thickness of the spaceframe (Al + Air)
132   Float_t frame = widmi - (widpl1 / 2);
133
134   // Definition of the inner part of the detector frame (Air) 
135   par_fr[0] = widmi - alframe / 2.;
136   par_fr[1] = widma - alframe / 2.;
137   par_fr[2] = zmax1;
138   par_fr[3] = (rmax - rmin) / 2;
139   pMC->Gsvolu("UTRI", "TRD1", idtmed[1302-1], par_fr, nparfr); 
140
141   // Some parameter for the chambers 
142   Float_t lendifc = (zmax1 - zmax2) / nmodul;
143   Float_t heightc = (rmax  - rmin ) / nmodul;
144   Float_t widdifc = (widma - widmi) / nmodul;
145
146   // Definition of the chambers 
147   Char_t ctagc[5], ctagi[5];
148   for (Int_t icham = 1; icham <= ncham; ++icham) {
149
150     // Carbon frame of the inner chambers (C) 
151     par_ch[0] = widmi + (icham-1) * widdifc - frame;
152     par_ch[1] = zleni   / 2.;
153     par_ch[2] = heightc / 2.;
154     sprintf(ctagc,"UCI%1d",icham);
155     pMC->Gsvolu(ctagc, "BOX ", idtmed[1307-1], par_ch, nparch);
156     // Inner part of the inner chambers (Air) 
157     par_ch[0] -= ccframe;
158     par_ch[1] -= ccframe;
159     sprintf(ctagc,"UII%1d",icham);
160     pMC->Gsvolu(ctagc, "BOX ", idtmed[1302-1], par_ch, nparch);
161
162     // Carbon frame of the neighbouring chambers (C) 
163     par_ch[0] = widmi + (icham-1) * widdifc - frame;
164     par_ch[1] = zlenn   / 2.;
165     par_ch[2] = heightc / 2.;
166     sprintf(ctagc,"UCN%1d",icham);
167     pMC->Gsvolu(ctagc, "BOX ", idtmed[1307-1], par_ch, nparch);
168     // Inner part of the neighbouring chambers (Air) 
169     par_ch[0] -= ccframe;
170     par_ch[1] -= ccframe;
171     sprintf(ctagc,"UIN%1d",icham);
172     pMC->Gsvolu(ctagc, "BOX ", idtmed[1302-1], par_ch, nparch);
173
174     // Carbon frame of the outer chambers (C) 
175     par_ch[0] = widmi + (icham-1) * widdifc - frame;
176     par_ch[1] = (icham - 6) * lendifc / 2. + zleno   / 2.;
177     par_ch[2] = heightc / 2.;
178     sprintf(ctagc,"UCO%1d",icham);
179     pMC->Gsvolu(ctagc, "BOX ", idtmed[1307-1], par_ch, nparch);
180     // Inner part of the outer chambers (Air) 
181     par_ch[0] -= ccframe;
182     par_ch[1] -= ccframe;
183     sprintf(ctagc,"UIO%1d",icham);
184     pMC->Gsvolu(ctagc, "BOX ", idtmed[1302-1], par_ch, nparch);
185
186   }
187
188   // Definition of the layers in each inner chamber 
189   par_ic[0] = -1.;
190   par_ic[1] = -1.;
191   // G10 layer (radiator layer)
192   par_ic[2] = sethick / 2;
193   pMC->Gsvolu("UT0I", "BOX ", idtmed[1313-1], par_ic, nparic);
194   // CO2 layer (radiator)
195   par_ic[2] = rathick / 2;
196   pMC->Gsvolu("UT1I", "BOX ", idtmed[1312-1], par_ic, nparic);
197   // PE layer (radiator)
198   par_ic[2] = pethick / 2;
199   pMC->Gsvolu("UT2I", "BOX ", idtmed[1303-1], par_ic, nparic);
200   // Mylar layer (entrance window + HV cathode) 
201   par_ic[2] = mythick / 2;
202   pMC->Gsvolu("UT3I", "BOX ", idtmed[1308-1], par_ic, nparic);
203   // Xe/Isobutane layer (gasvolume) 
204   par_ic[2] = xethick / 2.;
205   for (Int_t icham = 1; icham <= 6; ++icham) {
206     sprintf(ctagc,"UXI%1d",icham);
207     pMC->Gsvolu(ctagc, "BOX ", idtmed[1309-1], par_ic, nparic);
208   }
209   // Cu layer (pad plane)
210   par_ic[2] = cuthick / 2;
211   pMC->Gsvolu("UT5I", "BOX ", idtmed[1305-1], par_ic, nparic);
212   // G10 layer (support structure)
213   par_ic[2] = suthick / 2;
214   pMC->Gsvolu("UT6I", "BOX ", idtmed[1313-1], par_ic, nparic);
215   // Cu layer (FEE + signal lines)
216   par_ic[2] = fethick / 2;
217   pMC->Gsvolu("UT7I", "BOX ", idtmed[1305-1], par_ic, nparic);
218   // PE layer (cooling devices)
219   par_ic[2] = cothick / 2;
220   pMC->Gsvolu("UT8I", "BOX ", idtmed[1303-1], par_ic, nparic);
221   // Water layer (cooling)
222   par_ic[2] = wathick / 2;
223   pMC->Gsvolu("UT9I", "BOX ", idtmed[1314-1], par_ic, nparic);
224
225   // Definition of the layers in each neighbouring chamber 
226   par_nc[0] = -1.;
227   par_nc[1] = -1.;
228   // G10 layer (radiator layer)
229   par_nc[2] = sethick / 2;
230   pMC->Gsvolu("UT0N", "BOX ", idtmed[1313-1], par_nc, nparnc);
231   // CO2 layer (radiator)
232   par_nc[2] = rathick / 2;
233   pMC->Gsvolu("UT1N", "BOX ", idtmed[1312-1], par_nc, nparnc);
234   // PE layer (radiator)
235   par_nc[2] = pethick / 2;
236   pMC->Gsvolu("UT2N", "BOX ", idtmed[1303-1], par_nc, nparnc);
237   // Mylar layer (entrance window + HV cathode) 
238   par_nc[2] = mythick / 2;
239   pMC->Gsvolu("UT3N", "BOX ", idtmed[1308-1], par_nc, nparnc);
240   // Xe/Isobutane layer (gasvolume) 
241   par_nc[2] = xethick / 2.;
242   for (Int_t icham = 1; icham <= 6; ++icham) {
243     sprintf(ctagc,"UXN%1d",icham);
244     pMC->Gsvolu(ctagc, "BOX ", idtmed[1309-1], par_nc, nparnc);
245   }
246   // Cu layer (pad plane)
247   par_nc[2] = cuthick / 2;
248   pMC->Gsvolu("UT5N", "BOX ", idtmed[1305-1], par_nc, nparnc);
249   // G10 layer (support structure)
250   par_nc[2] = suthick / 2;
251   pMC->Gsvolu("UT6N", "BOX ", idtmed[1313-1], par_nc, nparnc);
252   // Cu layer (FEE + signal lines)
253   par_nc[2] = fethick / 2;
254   pMC->Gsvolu("UT7N", "BOX ", idtmed[1305-1], par_nc, nparnc);
255   // PE layer (cooling devices)
256   par_nc[2] = cothick / 2;
257   pMC->Gsvolu("UT8N", "BOX ", idtmed[1303-1], par_nc, nparnc);
258   // Water layer (cooling)
259   par_nc[2] = wathick / 2;
260   pMC->Gsvolu("UT9N", "BOX ", idtmed[1314-1], par_nc, nparnc);
261
262   // Definition of the layers in each outer chamber 
263   par_oc[0] = -1.;
264   par_oc[1] = -1.;
265   // G10 layer (radiator layer)
266   par_oc[2] = sethick / 2;
267   pMC->Gsvolu("UT0O", "BOX ", idtmed[1313-1], par_oc, nparoc);
268   // CO2 layer (radiator)
269   par_oc[2] = rathick / 2;
270   pMC->Gsvolu("UT1O", "BOX ", idtmed[1312-1], par_oc, nparoc);
271   // PE layer (radiator)
272   par_oc[2] = pethick / 2;
273   pMC->Gsvolu("UT2O", "BOX ", idtmed[1303-1], par_oc, nparoc);
274   // Mylar layer (entrance window + HV cathode) 
275   par_oc[2] = mythick / 2;
276   pMC->Gsvolu("UT3O", "BOX ", idtmed[1308-1], par_oc, nparoc);
277   // Xe/Isobutane layer (gasvolume) 
278   par_oc[2] = xethick / 2.;
279   for (Int_t icham = 1; icham <= 6; ++icham) {
280     sprintf(ctagc,"UXO%1d",icham);
281     pMC->Gsvolu(ctagc, "BOX ", idtmed[1309-1], par_oc, nparoc);
282   }
283   // Cu layer (pad plane)
284   par_oc[2] = cuthick / 2;
285   pMC->Gsvolu("UT5O", "BOX ", idtmed[1305-1], par_oc, nparoc);
286   // G10 layer (support structure)
287   par_oc[2] = suthick / 2;
288   pMC->Gsvolu("UT6O", "BOX ", idtmed[1313-1], par_oc, nparoc);
289   // Cu layer (FEE + signal lines)
290   par_oc[2] = fethick / 2;
291   pMC->Gsvolu("UT7O", "BOX ", idtmed[1305-1], par_oc, nparoc);
292   // PE layer (cooling devices)
293   par_oc[2] = cothick / 2;
294   pMC->Gsvolu("UT8O", "BOX ", idtmed[1303-1], par_oc, nparoc);
295   // Water layer (cooling)
296   par_oc[2] = wathick / 2;
297   pMC->Gsvolu("UT9O", "BOX ", idtmed[1314-1], par_oc, nparoc);
298
299   //////////////////////////////////////////////////////////////////////////
300   //     Positioning of Volumes   
301   //////////////////////////////////////////////////////////////////////////
302
303   // The rotation matrices 
304   AliMatrix(idmat[0],  90.,  90., 180.,   0.,  90.,   0.);
305   AliMatrix(idmat[1],  90.,  90.,   0.,   0.,  90.,   0.);
306
307   // Position of the layers in a chamber 
308   pMC->Gspos("UT2I", 1, "UT1I", 0., 0., pezpos, 0, "ONLY");
309   pMC->Gspos("UT2N", 1, "UT1N", 0., 0., pezpos, 0, "ONLY");
310   pMC->Gspos("UT2O", 1, "UT1O", 0., 0., pezpos, 0, "ONLY");
311   for (Int_t icham = 1; icham <= ncham; ++icham) {
312     // The inner chambers 
313     sprintf(ctagi,"UII%1d",icham);
314     sprintf(ctagc,"UXI%1d",icham);
315     pMC->Gspos("UT9I", icham, ctagi, 0., 0., wazpos, 0, "ONLY");
316     pMC->Gspos("UT8I", icham, ctagi, 0., 0., cozpos, 0, "ONLY");
317     pMC->Gspos("UT7I", icham, ctagi, 0., 0., fezpos, 0, "ONLY");
318     pMC->Gspos("UT6I", icham, ctagi, 0., 0., suzpos, 0, "ONLY");
319     pMC->Gspos("UT5I", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
320     pMC->Gspos(ctagc ,     1, ctagi, 0., 0., xezpos, 0, "ONLY");
321     pMC->Gspos("UT3I", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
322     pMC->Gspos("UT1I", icham, ctagi, 0., 0., razpos, 0, "ONLY");
323     pMC->Gspos("UT0I", icham, ctagi, 0., 0., sezpos, 0, "ONLY");
324     // The neighbouring chambers 
325     sprintf(ctagi,"UIN%1d",icham);
326     sprintf(ctagc,"UXN%1d",icham);
327     pMC->Gspos("UT9N", icham, ctagi, 0., 0., wazpos, 0, "ONLY");
328     pMC->Gspos("UT8N", icham, ctagi, 0., 0., cozpos, 0, "ONLY");
329     pMC->Gspos("UT7N", icham, ctagi, 0., 0., fezpos, 0, "ONLY");
330     pMC->Gspos("UT6N", icham, ctagi, 0., 0., suzpos, 0, "ONLY");
331     pMC->Gspos("UT5N", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
332     pMC->Gspos(ctagc ,     1, ctagi, 0., 0., xezpos, 0, "ONLY");
333     pMC->Gspos("UT3N", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
334     pMC->Gspos("UT1N", icham, ctagi, 0., 0., razpos, 0, "ONLY");
335     pMC->Gspos("UT0N", icham, ctagi, 0., 0., sezpos, 0, "ONLY");
336     // The outer chambers 
337     sprintf(ctagi,"UIO%1d",icham);
338     sprintf(ctagc,"UXO%1d",icham);
339     pMC->Gspos("UT9O", icham, ctagi, 0., 0., wazpos, 0, "ONLY");
340     pMC->Gspos("UT8O", icham, ctagi, 0., 0., cozpos, 0, "ONLY");
341     pMC->Gspos("UT7O", icham, ctagi, 0., 0., fezpos, 0, "ONLY");
342     pMC->Gspos("UT6O", icham, ctagi, 0., 0., suzpos, 0, "ONLY");
343     pMC->Gspos("UT5O", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
344     pMC->Gspos(ctagc ,     1, ctagi, 0., 0., xezpos, 0, "ONLY");
345     pMC->Gspos("UT3O", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
346     pMC->Gspos("UT1O", icham, ctagi, 0., 0., razpos, 0, "ONLY");
347     pMC->Gspos("UT0O", icham, ctagi, 0., 0., sezpos, 0, "ONLY");
348   }
349
350   // Position of the inner part of the chambers in the carbon-frames 
351   for (Int_t icham = 1; icham <= ncham; ++icham) {
352     xpos = 0.;
353     ypos = 0.;
354     zpos = 0.;
355     // The inner chambers 
356     sprintf(ctagi,"UII%1d",icham);
357     sprintf(ctagc,"UCI%1d",icham);
358     pMC->Gspos(ctagi, 1, ctagc, xpos, ypos, zpos, 0, "ONLY");
359     // The neighbouring chambers 
360     sprintf(ctagi,"UIN%1d",icham);
361     sprintf(ctagc,"UCN%1d",icham);
362     pMC->Gspos(ctagi, 1, ctagc, xpos, ypos, zpos, 0, "ONLY");
363     // The outer chambers 
364     sprintf(ctagi,"UIO%1d",icham);
365     sprintf(ctagc,"UCO%1d",icham);
366     pMC->Gspos(ctagi, 1, ctagc, xpos, ypos, zpos, 0, "ONLY");
367   }
368
369   // Position of the chambers in the full TRD-setup 
370   for (Int_t icham = 1; icham <= ncham; ++icham) {
371     // The inner chambers
372     xpos = 0.;
373     ypos = 0.;
374     zpos = (icham-0.5) * heightc - (rmax - rmin) / 2;
375     sprintf(ctagc,"UCI%1d",icham);
376     pMC->Gspos(ctagc, 1, "UTRI", xpos, ypos, zpos, 0, "ONLY");
377     // The neighbouring chambers
378     xpos = 0.;
379     ypos = (zleni + zlenn) / 2.;
380     zpos = (icham-0.5) * heightc - (rmax - rmin) / 2;
381     sprintf(ctagc,"UCN%1d",icham);
382     pMC->Gspos(ctagc, 1, "UTRI", xpos, ypos, zpos, 0, "ONLY");
383     ypos = -ypos;
384     sprintf(ctagc,"UCN%1d",icham);
385     pMC->Gspos(ctagc, 2, "UTRI", xpos, ypos, zpos, 0, "ONLY");
386     // The outer chambers 
387     xpos = 0.;
388     ypos = (zleni / 2. + zlenn + zmax2 + (icham-1) * lendifc) / 2.;
389     zpos = (icham-0.5) * heightc - (rmax-rmin)/2;
390     sprintf(ctagc,"UCO%1d",icham);
391     pMC->Gspos(ctagc, 1, "UTRI", xpos, ypos, zpos, 0, "ONLY");
392     ypos = -ypos;
393     sprintf(ctagc,"UCO%1d",icham);
394     pMC->Gspos(ctagc, 2, "UTRI", xpos, ypos, zpos, 0, "ONLY");
395   }
396
397   // Position of the inner part of the detector frame
398   xpos = (rmax + rmin) / 2;
399   ypos = 0.;
400   zpos = 0.;
401   pMC->Gspos("UTRI", 1, "UTRS", xpos, ypos, zpos, idmat[0], "ONLY");
402
403   // Position of the TRD mother volume in the ALICE experiment 
404   xpos = 0.;
405   ypos = 0.;
406   zpos = 0.;
407   pMC->Gspos("TRD ", 1, "ALIC", xpos, ypos, zpos,        0, "ONLY");
408
409 }
410
411 //_____________________________________________________________________________
412 void AliTRDv2::DrawModule()
413 {
414   //
415   // Draw a shaded view of the Transition Radiation Detector version 1
416   //
417
418   AliMC* pMC = AliMC::GetMC();
419   
420   // Set everything unseen
421   pMC->Gsatt("*", "seen", -1);
422    
423   // Set ALIC mother transparent
424   pMC->Gsatt("ALIC","SEEN",0);
425   
426   // Set the volumes visible
427   pMC->Gsatt("TRD ","SEEN",0);
428   pMC->Gsatt("UTRS","SEEN",0);
429   pMC->Gsatt("UTRI","SEEN",0);
430   Char_t ctag[5];
431   for (Int_t icham = 0; icham < ncham; ++icham) {
432     sprintf(ctag,"UCI%1d",icham+1);
433     pMC->Gsatt(ctag,"SEEN",0);
434     sprintf(ctag,"UCN%1d",icham+1);
435     pMC->Gsatt(ctag,"SEEN",0);
436     sprintf(ctag,"UCO%1d",icham+1);
437     pMC->Gsatt(ctag,"SEEN",0);
438     sprintf(ctag,"UII%1d",icham+1);
439     pMC->Gsatt(ctag,"SEEN",0);
440     sprintf(ctag,"UIN%1d",icham+1);
441     pMC->Gsatt(ctag,"SEEN",0);
442     sprintf(ctag,"UIO%1d",icham+1);
443     pMC->Gsatt(ctag,"SEEN",0);
444     sprintf(ctag,"UXI%1d",icham+1);
445     pMC->Gsatt(ctag,"SEEN",1);
446     sprintf(ctag,"UXN%1d",icham+1);
447     pMC->Gsatt(ctag,"SEEN",1);
448     sprintf(ctag,"UXO%1d",icham+1);
449     pMC->Gsatt(ctag,"SEEN",1);
450   }
451   pMC->Gsatt("UT1I","SEEN",1);
452   pMC->Gsatt("UT1N","SEEN",1);
453   pMC->Gsatt("UT1O","SEEN",1);
454
455   pMC->Gdopt("hide", "on");
456   pMC->Gdopt("shad", "on");
457   pMC->Gsatt("*", "fill", 7);
458   pMC->SetClipBox(".");
459   pMC->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
460   pMC->DefaultRange();
461   pMC->Gdraw("alic", 40, 30, 0, 12, 9.4, .021, .021);
462   pMC->Gdhead(1111, "Transition Radiation Detector Version 2");
463   pMC->Gdman(18, 4, "MAN");
464   pMC->Gdopt("hide", "off");
465 }
466
467 //_____________________________________________________________________________
468 void AliTRDv2::CreateMaterials()
469 {
470   //
471   // Create materials for the Transition Radiation Detector version 2
472   //
473   AliTRD::CreateMaterials();
474 }
475
476 //_____________________________________________________________________________
477 void AliTRDv2::Init() 
478 {
479   //
480   // Initialise Transition Radiation Detector after geometry has been built
481   //
482
483   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
484   const Float_t kPoti = 12.1;
485   // Maximum energy (50 keV);
486   const Float_t kEend = 50000.0;
487
488   AliTRD::Init();
489
490   AliMC* pMC = AliMC::GetMC();
491
492   // Get the sensitive volumes
493   Char_t ctag[5];
494   for (Int_t icham = 0; icham < ncham; ++icham) {
495     sprintf(ctag,"UXI%1d",icham+1);
496     fIdSensI[icham] = pMC->VolId(ctag);
497     sprintf(ctag,"UXN%1d",icham+1);
498     fIdSensN[icham] = pMC->VolId(ctag);
499     sprintf(ctag,"UXO%1d",icham+1);
500     fIdSensO[icham] = pMC->VolId(ctag);
501   }
502
503   Float_t Poti = TMath::Log(kPoti);
504   Float_t Eend = TMath::Log(kEend);
505
506   // Ermilova distribution for the delta-ray spectrum
507   fDeltaE  = new TF1("deltae",Ermilova,Poti,Eend,0);
508
509 }
510
511 //_____________________________________________________________________________
512 void AliTRDv2::StepManager()
513 {
514   //
515   // Called at every step in the Transition Radiation Detector version 2
516   //
517
518   Int_t          idSens, icSens, id;
519   Int_t          iPla, iCha, iSec;
520   Int_t          iOut;
521   Int_t          vol[3]; 
522   Int_t          iPid;
523
524   const Double_t kBig = 1.0E+12;
525
526   Float_t        hits[4];
527   Float_t        mom[4];
528   Float_t        random[1];
529   Float_t        charge;
530   Float_t        aMass;
531
532   Double_t       pTot;
533   Double_t       qTot;
534   Double_t       eDelta;
535   Double_t       betaGamma, pp;
536
537   TClonesArray  &lhits = *fHits;
538
539   AliMC* pMC = AliMC::GetMC();
540
541   // Ionization energy
542   const Float_t kWion    = 22.04;
543   // Maximum energy for e+ e- g for the step-size calculation
544   const Float_t kPTotMax = 0.002;
545   // Plateau value of the energy-loss for electron in xenon
546   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
547   //const Double_t kPlateau = 1.70;
548   // the averaged value (26/3/99)
549   const Float_t kPlateau = 1.55;
550   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
551   const Float_t kPrim    = 48.0;
552   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
553   const Float_t kPoti    = 12.1;
554
555   // Set the maximum step size to a very large number for all 
556   // neutral particles and those outside the driftvolume
557   pMC->SetMaxStep(kBig); 
558
559   // Use only charged tracks 
560   if (( pMC->TrackCharge()   ) &&
561       (!pMC->TrackStop()     ) && 
562       (!pMC->TrackDisappear())) {
563
564     // Find the sensitive volume
565     idSens = pMC->CurrentVol(0,icSens);
566     iPla   = 0;
567     iOut   = 0;
568     for (Int_t icham = 0; icham < ncham; ++icham) {
569       if (idSens == fIdSensI[icham]) {
570         iOut = 0;
571         iPla = icham + 1;
572       }
573       if (idSens == fIdSensN[icham]) {
574         iOut = 1;
575         iPla = icham + 1;
576       }
577       if (idSens == fIdSensO[icham]) {
578         iOut = 2;
579         iPla = icham + 1;
580       }
581     }
582
583     // Inside a sensitive volume?
584     if (iPla) {
585
586       // Calculate the energy of the delta-electrons
587       eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
588       eDelta = TMath::Max(eDelta,0.0);
589
590       // The number of secondary electrons created
591       qTot = (Double_t) ((Int_t) (eDelta / kWion) + 1);
592
593       // The sector number
594       id = pMC->CurrentVolOff(4,0,iSec);
595
596       // The chamber number
597       //   1: outer left
598       //   2: neighbouring left
599       //   3: inner
600       //   4: neighbouring right
601       //   5: outer right
602       id = pMC->CurrentVolOff(2,0,iCha);
603       if (iCha == 1) 
604         iCha = 3 + iOut;
605       else
606         iCha = 3 - iOut;
607
608       vol[0]  = iSec;
609       vol[1]  = iCha;
610       vol[2]  = iPla;
611
612       // Check on selected volumes
613       Int_t addthishit = 1;
614       if (fSensSelect) {
615         if ((fSensPlane)   && (vol[2] != fSensPlane  )) addthishit = 0;
616         if ((fSensChamber) && (vol[1] != fSensChamber)) addthishit = 0;
617         if ((fSensSector)  && (vol[0] != fSensSector )) addthishit = 0;
618       }
619
620       if (addthishit) {
621
622         // Add this hit
623         pMC->TrackPosition(hits);
624         hits[3] = qTot;
625         new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
626
627         // The energy loss according to Bethe Bloch
628         pMC->TrackMomentum(mom);
629         pTot = mom[3];
630         iPid = pMC->TrackPid();
631         if ( (iPid >  3) ||
632             ((iPid <= 3) && (pTot < kPTotMax))) {
633           aMass     = pMC->TrackMass();
634           betaGamma = pTot / aMass;
635           pp        = kPrim * BetheBloch(betaGamma);
636           // Take charge > 1 into account
637           charge = pMC->TrackCharge();
638           if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
639         }
640         // Electrons above 20 Mev/c are at the plateau
641         else {
642           pp = kPrim * kPlateau;
643         }
644       
645         // Calculate the maximum step size for the next tracking step
646         if (pp > 0) {
647           do 
648             pMC->Rndm(random,1);
649           while ((random[0] == 1.) || (random[0] == 0.));
650           pMC->SetMaxStep( - TMath::Log(random[0]) / pp);
651         }
652
653       }
654       else {
655         // set step size to maximal value
656         pMC->SetMaxStep(kBig); 
657       }
658
659     }
660
661   }
662
663 }
664
665 //_____________________________________________________________________________
666 Double_t AliTRDv2::BetheBloch(Double_t bg) 
667 {
668   //
669   // Parametrization of the Bethe-Bloch-curve
670   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
671   //
672
673   // The parameters have been adjusted to Xe-data found in:
674   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
675   //const Double_t kP1 = 0.76176E-1;
676   //const Double_t kP2 = 10.632;
677   //const Double_t kP3 = 3.17983E-6;
678   //const Double_t kP4 = 1.8631;
679   //const Double_t kP5 = 1.9479;
680
681   // This parameters have been adjusted to averaged values from GEANT
682   const Double_t kP1 = 7.17960e-02;
683   const Double_t kP2 = 8.54196;
684   const Double_t kP3 = 1.38065e-06;
685   const Double_t kP4 = 5.30972;
686   const Double_t kP5 = 2.83798;
687
688   if (bg > 0) {
689     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
690     Double_t aa = TMath::Power(yy,kP4);
691     Double_t bb = TMath::Power((1./bg),kP5);
692              bb = TMath::Log(kP3 + bb);
693     return ((kP2 - aa - bb)*kP1 / aa);
694   }
695   else
696     return 0;
697
698 }
699
700 //_____________________________________________________________________________
701 Double_t Ermilova(Double_t *x, Double_t *par)
702 {
703   //
704   // Calculates the delta-ray energy distribution according to Ermilova
705   // Logarithmic scale !
706   //
707
708   Double_t energy;
709   Double_t dpos;
710   Double_t dnde;
711
712   Int_t    pos1, pos2;
713
714   const Int_t nV = 31;
715
716   Float_t vxe[nV] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120  
717                     , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
718                     , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
719                     , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
720                     , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
721                     , 9.4727, 9.9035,10.3735,10.5966,10.8198
722                     ,11.5129 };
723
724   Float_t vye[nV] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
725                     , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
726                     ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
727                     ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
728                     ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
729                     ,  0.04 ,  0.023,  0.015,  0.011,  0.01
730                     ,  0.004 };
731
732   energy = x[0];
733
734   // Find the position 
735   pos1 = pos2 = 0;
736   dpos = 0;
737   do {
738     dpos = energy - vxe[pos2++];
739   } 
740   while (dpos > 0);
741   pos2--; 
742   if (pos2 > nV) pos2 = nV;
743   pos1 = pos2 - 1;
744
745   // Differentiate between the sampling points
746   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
747
748   return dnde;
749
750 }