72bf6b32c333706da7890c050de6417ab4a044b8
[u/mrichter/AliRoot.git] / ZDC / AliZDCv1.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.7  2000/01/19 17:17:40  fca
19
20 Revision 1.6  1999/09/29 09:24:35  fca
21 Introduction of the Copyright and cvs Log
22
23 */
24
25 ///////////////////////////////////////////////////////////////////////////////
26 //                                                                           //
27 //  Zero Degree Calorimeter                                                  //
28 //  This class contains the basic functions for the Time Of Flight           //
29 //  detector. Functions specific to one particular geometry are              //
30 //  contained in the derived classes                                         //
31 //                                                                           //
32 ///////////////////////////////////////////////////////////////////////////////
33
34 #include <TBRIK.h>
35 #include <TNode.h>
36 #include <TMath.h>
37
38 #include "stdio.h"
39 #include "AliZDCv1.h"
40 #include "AliRun.h"
41 #include "AliMC.h"
42 #include "AliCallf77.h"
43 #include "AliConst.h"
44 #include "AliPDG.h"
45  
46  
47 ClassImp(AliZDCv1)
48  
49
50 ///////////////////////////////////////////////////////////////////////////////
51 //                                                                           //
52 //  Zero Degree Calorimeter version 1                                        //
53 //                                                                           //
54 ///////////////////////////////////////////////////////////////////////////////
55
56 //_____________________________________________________________________________
57 AliZDCv1::AliZDCv1() : AliZDC()
58 {
59   //
60   // Default constructor for Zero Degree Calorimeter
61   //
62   fMedSensF1 = 0;
63   fMedSensF2 = 0;
64   fMedSensZN = 0;
65   fMedSensZP = 0;
66   fMedSensGR = 0;
67 }
68  
69 //_____________________________________________________________________________
70 AliZDCv1::AliZDCv1(const char *name, const char *title)
71   : AliZDC(name,title)
72 {
73   //
74   // Standard constructor for Zero Degree Calorimeter 
75   //
76   fMedSensF1 = 0;
77   fMedSensF2 = 0;
78   fMedSensZN = 0;
79   fMedSensZP = 0;
80   fMedSensGR = 0;
81 }
82  
83 //_____________________________________________________________________________
84 void AliZDCv1::CreateGeometry()
85 {
86   //
87   // Create the geometry for the Zero Degree Calorimeter version 1
88   //* Initialize COMMON block ZDC_CGEOM
89   //*
90
91   CreateBeamLine();
92   CreateZDC();
93 }
94   
95 //_____________________________________________________________________________
96 void AliZDCv1::CreateBeamLine()
97 {
98   
99   Float_t angle;
100   Float_t zq, conpar[9], elpar[3], tubpar[3];
101   Int_t im1, im2;
102   Float_t zd1, zd2;
103   
104   
105   Int_t *idtmed = fIdtmed->GetArray();
106   
107   // -- Mother of the ZDC 
108   
109   conpar[0] = 0.;
110   conpar[1] = 360.;
111   conpar[2] = 2.;
112   conpar[3] = 1921.6;
113   conpar[4] = 0.;
114   conpar[5] = 55.;
115   conpar[6] = 13060.;
116   conpar[7] = 0.;
117   conpar[8] = 55.;
118   gMC->Gsvolu("ZDC ", "PCON", idtmed[10], conpar, 9);
119   gMC->Gspos("ZDC ", 1, "ALIC", 0., 0., 0., 0, "ONLY");
120
121   // -- FIRST SECTION OF THE BEAM PIPE (from compensator dipole to 
122   //    beginning of D1) 
123   
124   zd1 = 1921.6;
125   
126   tubpar[0] = 6.3/2.;
127   tubpar[1] = 6.7/2.;
128   tubpar[2] = 3916.7/2.;
129   gMC->Gsvolu("P001", "TUBE", idtmed[5], tubpar, 3);
130   gMC->Gspos("P001", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
131   
132   //-- SECOND SECTION OF THE BEAM PIPE (FROM THE END OF D1 TO THE BEGINNING OF
133   //    D2) 
134   
135   //-- FROM MAGNETIC BEGINNING OG D1 TO MAGNETIC END OF D1 + 23.5 cm
136   //-- Elliptic pipe
137   
138   zd1 = 6310.8-472.5;
139   
140   elpar[0] = 6.84/2.;
141   elpar[1] = 5.86/2.;
142   elpar[2] = 945./2.;
143   gMC->Gsvolu("E001", "ELTU", idtmed[5], elpar, 3);
144   gMC->Gspos("E001", 1, "ZDC ", 0., 0., elpar[2] + zd1, 0, "ONLY");
145   
146   elpar[0] = 6.44/2.;
147   elpar[1] = 5.46/2.;
148   elpar[2] = 945./2.;
149   gMC->Gsvolu("E002", "ELTU", idtmed[10], elpar, 3);
150   gMC->Gspos("E002", 1, "E001", 0., 0., 0., 0, "ONLY");
151
152   zd1 += 2.*elpar[2];
153   
154   elpar[0] = 6.84/2.;
155   elpar[1] = 5.86/2.;
156   elpar[2] = 13.5/2.;
157   gMC->Gsvolu("E003", "ELTU", idtmed[5], elpar, 3);
158   gMC->Gspos("E002", 1, "ZDC ", 0., 0., elpar[2] + zd1, 0, "ONLY");
159   
160   elpar[0] = 6.44/2.;
161   elpar[1] = 5.46/2.;
162   elpar[2] = 13.5/2.;
163   gMC->Gsvolu("E004", "ELTU", idtmed[10], elpar, 3);
164   gMC->Gspos("E004", 1, "E003", 0., 0., 0., 0, "ONLY");
165
166   zd1 += 2.*elpar[2];
167   
168   conpar[0] = 25./2.;
169   conpar[1] = 6.44/2.;
170   conpar[2] = 6.84/2.;
171   conpar[3] = 10./2.;
172   conpar[4] = 10.4/2.;
173   gMC->Gsvolu("C001", "CONE", idtmed[5], conpar, 5);
174   gMC->Gspos("C001", 1, "ZDC ", 0., 0., conpar[0] + zd1, 0, "ONLY");
175
176   zd1 += 2.*conpar[0];
177   
178   tubpar[0] = 10./2.;
179   tubpar[1] = 10.4/2.;
180   tubpar[2] = 50./2.;
181   gMC->Gsvolu("P002", "TUBE", idtmed[5], tubpar, 3);
182   gMC->Gspos("P002", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
183   
184   zd1 += tubpar[2] * 2.;
185   
186   tubpar[0] = 10./2.;
187   tubpar[1] = 10.4/2.;
188   tubpar[2] = 10./2.;
189   gMC->Gsvolu("P003", "TUBE", idtmed[5], tubpar, 3);
190   gMC->Gspos("P003", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
191   
192   zd1 += tubpar[2] * 2.;
193   
194   tubpar[0] = 10./2.;
195   tubpar[1] = 10.4/2.;
196   tubpar[2] = 3.16/2.;
197   gMC->Gsvolu("P004", "TUBE", idtmed[5], tubpar, 3);
198   gMC->Gspos("P004", 1, "ZDC ", 0., 0., tubpar[0] + zd1, 0, "ONLY");
199   
200   zd1 += tubpar[2] * 2.;
201   
202   tubpar[0] = 10.0/2.;
203   tubpar[1] = 10.4/2;
204   tubpar[2] = 190./2.;
205   gMC->Gsvolu("P005", "TUBE", idtmed[5], tubpar, 3);
206   gMC->Gspos("P005", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
207   
208   zd1 += tubpar[2] * 2.;
209   
210   conpar[0] = 30./2.;
211   conpar[1] = 10./2.;
212   conpar[2] = 10.4/2.;
213   conpar[3] = 20.6/2.;
214   conpar[4] = 21./2.;
215   gMC->Gsvolu("P006", "CONE", idtmed[5], conpar, 5);
216   gMC->Gspos("P006", 1, "ZDC ", 0., 0., conpar[0] + zd1, 0, "ONLY");
217   
218   zd1 += conpar[0] * 2.;
219   
220   tubpar[0] = 20.6/2.;
221   tubpar[1] = 21./2.;
222   tubpar[2] = 450./2.;
223   gMC->Gsvolu("P007", "TUBE", idtmed[5], tubpar, 3);
224   gMC->Gspos("P007", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
225   
226   zd1 += tubpar[2] * 2.;
227   
228   conpar[0] = 13.6/2.;
229   conpar[1] = 20.6/2.;
230   conpar[2] = 21./2.;
231   conpar[3] = 25.4/2.;
232   conpar[4] = 25.8/2.;
233   gMC->Gsvolu("P008", "CONE", idtmed[5], conpar, 5);
234   gMC->Gspos("P008", 1, "ZDC ", 0., 0., conpar[0] + zd1, 0, "ONLY");
235   
236   zd1 += conpar[0] * 2.;
237   
238   tubpar[0] = 25.4/2.;
239   tubpar[1] = 25.8/2.;
240   tubpar[2] = 205.8/2.;
241   gMC->Gsvolu("P009", "TUBE", idtmed[5], tubpar, 3);
242   gMC->Gspos("P009", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
243   
244   zd1 += tubpar[2] * 2.;
245   
246   tubpar[0] = 50./2.;
247   tubpar[1] = 50.4/2.;
248   tubpar[2] = 505.4/2.;
249   gMC->Gsvolu("P010", "TUBE", idtmed[5], tubpar, 3);
250   gMC->Gspos("P010", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
251   
252   zd1 += tubpar[2] * 2.;
253   
254   tubpar[0] = 50./2.;
255   tubpar[1] = 50.4/2.;
256   tubpar[2] = 700./2.;
257   gMC->Gsvolu("P011", "TUBE", idtmed[5], tubpar, 3);
258   gMC->Gspos("P011", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
259   
260   zd1 += tubpar[2] * 2.;
261   
262   tubpar[0] = 50./2.;
263   tubpar[1] = 50.4/2.;
264   tubpar[2] = 778.5/2.;
265   gMC->Gsvolu("P012", "TUBE", idtmed[5], tubpar, 3);
266   gMC->Gspos("P012", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
267   
268   zd1 += tubpar[2] * 2.;
269   
270   conpar[0] = 14.18/2.;
271   conpar[1] = 50./2.;
272   conpar[2] = 50.4/2.;
273   conpar[3] = 55./2.;
274   conpar[4] = 55.4/2.;
275   gMC->Gsvolu("P013", "CONE", idtmed[5], conpar, 5);
276   gMC->Gspos("P013", 1, "ZDC ", 0., 0., conpar[0] + zd1, 0, "ONLY");
277   
278   zd1 += conpar[0] * 2.;
279   
280   tubpar[0] = 55./2.;
281   tubpar[1] = 55.4/2.;
282   tubpar[2] = 730./2.;
283   gMC->Gsvolu("P014", "TUBE", idtmed[5], tubpar, 3);
284   gMC->Gspos("P014", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
285   
286   zd1 += tubpar[2] * 2.;
287   
288   conpar[0] = 36.86/2.;
289   conpar[1] = 55./2.;
290   conpar[2] = 55.4/2.;
291   conpar[3] = 68./2.;
292   conpar[4] = 68.4/2.;
293   gMC->Gsvolu("P015", "CONE", idtmed[5], conpar, 5);
294   gMC->Gspos("P015", 1, "ZDC ", 0., 0., conpar[0] + zd1, 0, "ONLY");
295   
296   zd1 += conpar[0] * 2.;
297   
298   tubpar[0] = 68./2.;
299   tubpar[1] = 68.4/2.;
300   tubpar[2] = 927.3/2.;
301   gMC->Gsvolu("P016", "TUBE", idtmed[5], tubpar, 3);
302   gMC->Gspos("P016", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
303   
304   zd1 += tubpar[2] * 2.;
305   
306   tubpar[0] = 0./2.;
307   tubpar[1] = 68.4/2.;
308   tubpar[2] = 0.2/2.;
309   gMC->Gsvolu("P017", "TUBE", idtmed[5], tubpar, 3);
310   gMC->Gspos("P017", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
311   
312   zd1 += tubpar[2] * 2.;
313   
314   tubpar[0] = 0./2.;
315   tubpar[1] = 5./2.;
316   tubpar[2] = 0.2/2.;
317   gMC->Gsvolu("Q017", "TUBE", idtmed[10], tubpar, 3);
318   
319   //-- Position Q017 inside P017
320   gMC->Gspos("Q017", 1, "P017", -7.7, 0., 0., 0, "ONLY");
321   
322   tubpar[0] = 0./2.;
323   tubpar[1] = 7./2.;
324   tubpar[2] = 0.2/2.;
325   gMC->Gsvolu("R017", "TUBE", idtmed[10], tubpar, 3);
326   
327   //-- Position R017 inside P017
328   gMC->Gspos("R017", 1, "P017", 7.7, 0., 0., 0, "ONLY");
329   
330   //-- BEAM PIPE BETWEEN END OF CONICAL PIPE AND BEGINNING OF D2
331   
332   tubpar[0] = 5./2.;
333   tubpar[1] = 5.4/2.;
334   tubpar[2] = 678./2.;
335   gMC->Gsvolu("P018", "TUBE", idtmed[5], tubpar, 3);
336
337   tubpar[0] = 7./2.;
338   tubpar[1] = 7.4/2.;
339   tubpar[2] = 678./2.;
340   gMC->Gsvolu("P019", "TUBE", idtmed[5], tubpar, 3);
341   
342   // -- ROTATE PIPES 
343
344   AliMatrix(im1, 90.-0.071, 0., 90., 90., .071, 180.);
345   angle = .071*kDegrad;
346   gMC->Gspos("P018", 1, "ZDC ", TMath::Sin(angle) * 645. / 2. - 9.7 + 
347                TMath::Sin(angle) * 945. / 2., 0., tubpar[2] + zd1, im1, "ONLY");
348   AliMatrix(im2, 90.+0.071, 0., 90., 90., .071, 0.);
349   gMC->Gspos("P019", 1, "ZDC ", 9.7 - TMath::Sin(angle) * 645. / 2., 0., 
350                tubpar[2] + zd1, im2, "ONLY");
351   
352   // --  END OF BEAM PIPE VOLUME DEFINITION. MAGNET DEFINITION FOLLOWS 
353   //     (LHC OPTICS 6) 
354   
355   // -- COMPENSATOR DIPOLE (MBXW) 
356   //     GAP (VACUUM WITH MAGNETIC FIELD) 
357   
358   tubpar[0] = 0.;
359   tubpar[1] = 4.5;
360   tubpar[2] = 340./2.;
361   gMC->Gsvolu("MBXW", "TUBE", idtmed[11], tubpar, 3);
362   gMC->Gspos("MBXW", 1, "ZDC ", 0., 0., tubpar[2] + 805., 0, "ONLY");
363   
364   // --  YOKE (IRON WITHOUT MAGNETIC FIELD) 
365   
366   tubpar[0] = 4.5;
367   tubpar[1] = 55.;
368   tubpar[2] = 340./2.;
369   gMC->Gsvolu("YMBX", "TUBE", idtmed[5], tubpar, 3);
370   gMC->Gspos("YMBX", 1, "ZDC ", 0., 0., tubpar[2] + 805., 0, "ONLY");
371   
372   // -- COMPENSATOR DIPOLE (MCBWA) 
373   //     GAP (VACUUM WITH MAGNETIC FIELD) 
374   
375   tubpar[0] = 0.;
376   tubpar[1] = 4.5;
377   tubpar[2] = 170./2.;
378   gMC->Gsvolu("MCBW", "TUBE", idtmed[11], tubpar, 3);
379   gMC->Gspos("MCBW", 1, "ZDC ", 0., 0., tubpar[2] + 1921.6, 0, "ONLY");
380   
381   // --  YOKE (IRON WITHOUT MAGNETIC FIELD) 
382   
383   tubpar[0] = 4.5;
384   tubpar[1] = 55.;
385   tubpar[2] = 170./2.;
386   gMC->Gsvolu("YMCB", "TUBE", idtmed[5], tubpar, 3);
387   gMC->Gspos("YMCB", 1, "ZDC ", 0., 0., tubpar[2] + 1921.6, 0, "ONLY");
388   
389   // -- INNER TRIPLET 
390   
391   zq = 2300.;
392   
393   // -- DEFINE MQXL AND MQX QUADRUPOLE ELEMENT 
394   
395   //     MQXL 
396   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
397   
398   tubpar[0] = 0.;
399   tubpar[1] = 3.5;
400   tubpar[2] = 630./2.;
401   gMC->Gsvolu("MQXL", "TUBE", idtmed[11], tubpar, 3);
402   
403   // --  YOKE 
404   
405   tubpar[0] = 3.5;
406   tubpar[1] = 22.;
407   tubpar[2] = 630./2.;
408   gMC->Gsvolu("YMQL", "TUBE", idtmed[5], tubpar, 3);
409   
410   gMC->Gspos("MQXL", 1, "ZDC ", 0., 0., tubpar[2] + zq, 0, "ONLY");
411   gMC->Gspos("YMQL", 1, "ZDC ", 0., 0., tubpar[2] + zq, 0, "ONLY");
412   
413   gMC->Gspos("MQXL", 2, "ZDC ", 0., 0., tubpar[2] + zq + 2430., 0, "ONLY");
414   gMC->Gspos("YMQL", 2, "ZDC ", 0., 0., tubpar[2] + zq + 2430., 0, "ONLY");
415   
416   // --  MQX 
417   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
418   
419   tubpar[0] = 0.;
420   tubpar[1] = 3.5;
421   tubpar[2] = 550./2.;
422   gMC->Gsvolu("MQX ", "TUBE", idtmed[11], tubpar, 3);
423   
424   // --  YOKE 
425   
426   tubpar[0] = 3.5;
427   tubpar[1] = 22.;
428   tubpar[2] = 550./2.;
429   gMC->Gsvolu("YMQ ", "TUBE", idtmed[5], tubpar, 3);
430   
431   gMC->Gspos("MQX ", 1, "ZDC ", 0., 0., tubpar[2] + zq + 880.,  0, "ONLY");
432   gMC->Gspos("YMQ ", 1, "ZDC ", 0., 0., tubpar[2] + zq + 880.,  0, "ONLY");
433   
434   gMC->Gspos("MQX ", 2, "ZDC ", 0., 0., tubpar[2] + zq + 1530., 0, "ONLY");
435   gMC->Gspos("YMQ ", 2, "ZDC ", 0., 0., tubpar[2] + zq + 1530., 0, "ONLY");
436   
437   // -- SEPARATOR DIPOLE D1 
438   
439   zd1 = 5838.3;
440   
441   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
442   
443   tubpar[0] = 0.;
444   tubpar[1] = 7.5/2.;
445   tubpar[2] = 945./2.;
446   gMC->Gsvolu("D1  ", "TUBE", idtmed[11], tubpar, 3);
447   
448   // --  YOKE 
449   
450   tubpar[0] = 0.;
451   tubpar[1] = 110./2;
452   tubpar[2] = 945./2.;
453   gMC->Gsvolu("YD1 ", "TUBE", idtmed[5], tubpar, 3);
454   
455   gMC->Gspos("YD1 ", 1, "ZDC ", 0., 0., tubpar[2] + zd1, 0, "ONLY");
456   gMC->Gspos("D1  ", 1, "YD1 ", 0., 0., 0., 0, "ONLY");
457   
458   // -- DIPOLE D2 
459   
460   zd2 = 12147.6;
461   
462   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
463   
464   tubpar[0] = 0.;
465   tubpar[1] = 7.5/2.;
466   tubpar[2] = 945./2.;
467   gMC->Gsvolu("D2  ", "TUBE", idtmed[11], tubpar, 3);
468   
469   // --  YOKE 
470   
471   tubpar[0] = 0.;
472   tubpar[1] = 55.;
473   tubpar[2] = 945./2.;
474   gMC->Gsvolu("YD2 ", "TUBE", idtmed[5], tubpar, 3);
475   
476   gMC->Gspos("YD2 ", 1, "ZDC ", 0., 0., tubpar[2] + zd2, 0, "ONLY");
477   
478   gMC->Gspos("D2  ", 1, "YD2 ", -9.7, 0., 0., 0, "ONLY");
479   gMC->Gspos("D2  ", 2, "YD2 ",  9.7, 0., 0., 0, "ONLY");
480   
481   // -- END OF MAGNET DEFINITION 
482 }
483   
484 //_____________________________________________________________________________
485 void AliZDCv1::CreateZDC()
486 {
487   
488   Int_t *idtmed = fIdtmed->GetArray();
489   
490   //-- Create calorimeters geometry
491   
492   //--> Neutron calorimeter (ZN) 
493   
494   gMC->Gsvolu("ZNEU", "BOX ", idtmed[1], fDimZN, 3); // Passive material  
495   gMC->Gsvolu("ZNF1", "TUBE", idtmed[3], fFibZN, 3); // Active material
496   gMC->Gsvolu("ZNF2", "TUBE", idtmed[4], fFibZN, 3); 
497   gMC->Gsvolu("ZNF3", "TUBE", idtmed[4], fFibZN, 3); 
498   gMC->Gsvolu("ZNF4", "TUBE", idtmed[3], fFibZN, 3); 
499   gMC->Gsvolu("ZNG1", "BOX ", idtmed[12], fGrvZN, 3); // Empty grooves 
500   gMC->Gsvolu("ZNG2", "BOX ", idtmed[12], fGrvZN, 3); 
501   gMC->Gsvolu("ZNG3", "BOX ", idtmed[12], fGrvZN, 3); 
502   gMC->Gsvolu("ZNG4", "BOX ", idtmed[12], fGrvZN, 3); 
503   
504   // Divide ZNEU in towers (for hits purposes) 
505   
506   gMC->Gsdvn("ZNTX", "ZNEU", fTowZN[0], 1); // x-tower 
507   gMC->Gsdvn("ZN1 ", "ZNTX", fTowZN[1], 2); // y-tower
508   
509   //-- Divide ZN1 in minitowers 
510   //  fDivZN[0]= NUMBER OF FIBERS PER TOWER ALONG X-AXIS, 
511   //  fDivZN[1]= NUMBER OF FIBERS PER TOWER ALONG Y-AXIS
512   //  (4 fibres per minitower) 
513   
514   gMC->Gsdvn("ZNSL", "ZN1 ", fDivZN[1], 2); // Slices 
515   gMC->Gsdvn("ZNST", "ZNSL", fDivZN[0], 1); // Sticks
516   
517   // --- Position the empty grooves in the sticks (4 grooves per stick)
518   Float_t dx = fDimZN[0] / fDivZN[0] / 4.;
519   Float_t dy = fDimZN[1] / fDivZN[1] / 4.;
520   
521   gMC->Gspos("ZNG1", 1, "ZNST", 0.-dx, 0.+dy, 0., 0, "ONLY");
522   gMC->Gspos("ZNG2", 1, "ZNST", 0.+dx, 0.+dy, 0., 0, "ONLY");
523   gMC->Gspos("ZNG3", 1, "ZNST", 0.-dx, 0.-dy, 0., 0, "ONLY");
524   gMC->Gspos("ZNG4", 1, "ZNST", 0.+dx, 0.-dy, 0., 0, "ONLY");
525   
526   // --- Position the fibers in the grooves 
527   gMC->Gspos("ZNF1", 1, "ZNG1", 0., 0., 0., 0, "ONLY");
528   gMC->Gspos("ZNF2", 1, "ZNG2", 0., 0., 0., 0, "ONLY");
529   gMC->Gspos("ZNF3", 1, "ZNG3", 0., 0., 0., 0, "ONLY");
530   gMC->Gspos("ZNF4", 1, "ZNG4", 0., 0., 0., 0, "ONLY");
531   
532   // --- Position the neutron calorimeter in ZDC 
533   gMC->Gspos("ZNEU", 1, "ZDC ", fPosZN[0], fPosZN[1], fPosZN[2] + fDimZN[2], 0, "ONLY");
534   
535   //--> Proton calorimeter 
536   
537   gMC->Gsvolu("ZPRO", "BOX ", idtmed[2], fDimZP, 3); // Passive material
538   gMC->Gsvolu("ZPF1", "TUBE", idtmed[3], fFibZP, 3); // Active material
539   gMC->Gsvolu("ZPF2", "TUBE", idtmed[4], fFibZP, 3); 
540   gMC->Gsvolu("ZPF3", "TUBE", idtmed[4], fFibZP, 3); 
541   gMC->Gsvolu("ZPF4", "TUBE", idtmed[3], fFibZP, 3); 
542   gMC->Gsvolu("ZPG1", "BOX ", idtmed[12], fGrvZP, 3); // Empty grooves 
543   gMC->Gsvolu("ZPG2", "BOX ", idtmed[12], fGrvZP, 3); 
544   gMC->Gsvolu("ZPG3", "BOX ", idtmed[12], fGrvZP, 3); 
545   gMC->Gsvolu("ZPG4", "BOX ", idtmed[12], fGrvZP, 3); 
546     
547   //-- Divide ZPRO in towers(for hits purposes) 
548   
549   gMC->Gsdvn("ZPTX", "ZPRO", fTowZP[0], 1); // x-tower 
550   gMC->Gsdvn("ZP1 ", "ZPTX", fTowZP[1], 2); // y-tower
551   
552   
553   //-- Divide ZP1 in minitowers 
554   //  fDivZP[0]= NUMBER OF FIBERS ALONG X-AXIS PER MINITOWER, 
555   //  fDivZP[1]= NUMBER OF FIBERS ALONG Y-AXIS PER MINITOWER
556   //  (4 fiber per minitower) 
557   
558   gMC->Gsdvn("ZPSL", "ZP1 ", fDivZP[1], 2); // Slices 
559   gMC->Gsdvn("ZPST", "ZPSL", fDivZP[0], 1); // Sticks
560   
561   // --- Position the empty grooves in the sticks (4 grooves per stick)
562   dx = fDimZP[0] / fTowZP[0] / fDivZP[0] / 2.;
563   dy = fDimZP[1] / fTowZP[1] / fDivZP[1] / 2.;
564   
565   gMC->Gspos("ZPG1", 1, "ZPST", 0.-dx, 0.+dy, 0., 0, "ONLY");
566   gMC->Gspos("ZPG2", 1, "ZPST", 0.+dx, 0.+dy, 0., 0, "ONLY");
567   gMC->Gspos("ZPG3", 1, "ZPST", 0.-dx, 0.-dy, 0., 0, "ONLY");
568   gMC->Gspos("ZPG4", 1, "ZPST", 0.+dx, 0.-dy, 0., 0, "ONLY");
569   
570   // --- Position the fibers in the grooves 
571   gMC->Gspos("ZPF1", 1, "ZPG1", 0., 0., 0., 0, "ONLY");
572   gMC->Gspos("ZPF2", 1, "ZPG2", 0., 0., 0., 0, "ONLY");
573   gMC->Gspos("ZPF3", 1, "ZPG3", 0., 0., 0., 0, "ONLY");
574   gMC->Gspos("ZPF4", 1, "ZPG4", 0., 0., 0., 0, "ONLY");
575   
576
577   // --- Position the proton calorimeter in ZDC 
578   gMC->Gspos("ZPRO", 1, "ZDC ", fPosZP[0], fPosZP[1], fPosZP[2] + fDimZP[2], 0, "ONLY");
579   
580 }
581  
582 //_____________________________________________________________________________
583 void AliZDCv1::DrawModule()
584 {
585   //
586   // Draw a shaded view of the Zero Degree Calorimeter version 1
587   //
588
589   // Set everything unseen
590   gMC->Gsatt("*", "seen", -1);
591   // 
592   // Set ALIC mother transparent
593   gMC->Gsatt("ALIC","SEEN",0);
594   //
595   // Set the volumes visible
596   gMC->Gsatt("ZDC ","SEEN",0);
597   gMC->Gsatt("P001","SEEN",1);
598   gMC->Gsatt("E001","SEEN",1);
599   gMC->Gsatt("E002","SEEN",1);
600   gMC->Gsatt("E003","SEEN",1);
601   gMC->Gsatt("E004","SEEN",1);
602   gMC->Gsatt("C001","SEEN",1);
603   gMC->Gsatt("P002","SEEN",1);
604   gMC->Gsatt("P003","SEEN",1);
605   gMC->Gsatt("P004","SEEN",1);
606   gMC->Gsatt("P005","SEEN",1);
607   gMC->Gsatt("P006","SEEN",1);
608   gMC->Gsatt("P007","SEEN",1);
609   gMC->Gsatt("P008","SEEN",1);
610   gMC->Gsatt("P009","SEEN",1);
611   gMC->Gsatt("P010","SEEN",1);
612   gMC->Gsatt("P011","SEEN",1);
613   gMC->Gsatt("P012","SEEN",1);
614   gMC->Gsatt("P013","SEEN",1);
615   gMC->Gsatt("P014","SEEN",1);
616   gMC->Gsatt("P015","SEEN",1);
617   gMC->Gsatt("P016","SEEN",1);
618   gMC->Gsatt("P017","SEEN",1);
619   gMC->Gsatt("Q017","SEEN",1);
620   gMC->Gsatt("R017","SEEN",1);
621   gMC->Gsatt("P018","SEEN",1);
622   gMC->Gsatt("P019","SEEN",1);
623   gMC->Gsatt("MBXW","SEEN",1);
624   gMC->Gsatt("YMBX","SEEN",1);
625   gMC->Gsatt("MCBW","SEEN",1);
626   gMC->Gsatt("YMCB","SEEN",1);
627   gMC->Gsatt("MQXL","SEEN",1);
628   gMC->Gsatt("YMQL","SEEN",1);
629   gMC->Gsatt("MQX ","SEEN",1);
630   gMC->Gsatt("YMQ ","SEEN",1);
631   gMC->Gsatt("D1  ","SEEN",1);
632   gMC->Gsatt("YD1 ","SEEN",1);
633   gMC->Gsatt("D2  ","SEEN",1);
634   gMC->Gsatt("YD2 ","SEEN",1);
635   gMC->Gsatt("ZNEU","SEEN",0);
636   gMC->Gsatt("ZNF1","SEEN",0);
637   gMC->Gsatt("ZNF2","SEEN",0);
638   gMC->Gsatt("ZNF3","SEEN",0);
639   gMC->Gsatt("ZNF4","SEEN",0);
640   gMC->Gsatt("ZNG1","SEEN",0);
641   gMC->Gsatt("ZNG2","SEEN",0);
642   gMC->Gsatt("ZNG3","SEEN",0);
643   gMC->Gsatt("ZNG4","SEEN",0);
644   gMC->Gsatt("ZNTX","SEEN",0);
645   gMC->Gsatt("ZN1 ","COLO",2); 
646   gMC->Gsatt("ZN1 ","SEEN",1);
647   gMC->Gsatt("ZNSL","SEEN",0);
648   gMC->Gsatt("ZNST","SEEN",0);
649   gMC->Gsatt("ZPRO","SEEN",0);
650   gMC->Gsatt("ZPF1","SEEN",0);
651   gMC->Gsatt("ZPF2","SEEN",0);
652   gMC->Gsatt("ZPF3","SEEN",0);
653   gMC->Gsatt("ZPF4","SEEN",0);
654   gMC->Gsatt("ZPG1","SEEN",0);
655   gMC->Gsatt("ZPG2","SEEN",0);
656   gMC->Gsatt("ZPG3","SEEN",0);
657   gMC->Gsatt("ZPG4","SEEN",0);
658   gMC->Gsatt("ZPTX","SEEN",0);
659   gMC->Gsatt("ZP1 ","COLO",2); 
660   gMC->Gsatt("ZP1 ","SEEN",1);
661   gMC->Gsatt("ZPSL","SEEN",0);
662   gMC->Gsatt("ZPST","SEEN",0);
663   
664   //
665   gMC->Gdopt("hide", "on");
666   gMC->Gdopt("shad", "on");
667   gMC->Gsatt("*", "fill", 7);
668   gMC->SetClipBox(".");
669   gMC->SetClipBox("*", 0, 100, -100, 100, 12000, 16000);
670   gMC->DefaultRange();
671   gMC->Gdraw("alic", 40, 30, 0, 488, 220, .07, .07);
672   gMC->Gdhead(1111, "Zero Degree Calorimeter Version 1");
673   gMC->Gdman(18, 4, "MAN");
674 }
675
676 //_____________________________________________________________________________
677 void AliZDCv1::CreateMaterials()
678 {
679   //
680   // Create Materials for the Zero Degree Calorimeter
681   //
682   // Origin    : E. Scomparin 
683   
684   Int_t *idtmed = fIdtmed->GetArray();
685   
686   Float_t dens, ubuf[1], wmat[2];
687   Int_t isvolActive;
688   Float_t a[2];
689   Int_t i;
690   Float_t z[2], epsil=0.001, stmin=0.01;
691   Int_t isvol;
692   Float_t fieldm = gAlice->Field()->Max();
693   Int_t inofld;
694   Float_t deemax=-1;
695   Float_t tmaxfd=gAlice->Field()->Max();
696   Int_t isxfld = gAlice->Field()->Integ();
697   Float_t stemax;
698   
699   // --- Store in UBUF r0 for nuclear radius calculation R=r0*A**1/3 
700
701   // --- Tantalum -> ZN passive material
702   ubuf[0] = 1.1;
703   AliMaterial(1, "TANT", 180.95, 73., 16.65, .4, 11.9, ubuf, 1);
704     
705   // --- Tungsten 
706 //  ubuf[0] = 1.11;
707 //  AliMaterial(1, "TUNG", 183.85, 74., 19.3, .35, 10.3, ubuf, 1);
708   
709   // --- Brass (CuZn)  -> ZP passive material
710   dens = 8.48;
711   a[0] = 63.546;
712   a[1] = 65.39;
713   z[0] = 29.;
714   z[1] = 30.;
715   wmat[0] = .63;
716   wmat[1] = .37;
717   AliMixture(2, "BRASS               ", a, z, dens, 2, wmat);
718   
719   // --- SiO2 
720   dens = 2.64;
721   a[0] = 28.086;
722   a[1] = 15.9994;
723   z[0] = 14.;
724   z[1] = 8.;
725   wmat[0] = 1.;
726   wmat[1] = 2.;
727   AliMixture(3, "SIO2                ", a, z, dens, -2, wmat);
728
729   // --- Copper 
730 //  ubuf[0] = 1.1;
731 //  AliMaterial(7, "COPP", 63.54, 29., 8.96, 1.4, 0., ubuf, 1);
732   
733   
734   // --- Lead 
735 //  ubuf[0] = 1.12;
736 //  AliMaterial(6, "LEAD", 207.19, 82., 11.35, .56, 18.5, ubuf, 1);
737   
738   // --- Iron 
739   ubuf[0] = 1.1;
740   AliMaterial(5, "IRON", 55.85, 26., 7.87, 1.76, 0., ubuf, 1);
741   
742   // --- Vacuum (no magnetic field) 
743   AliMaterial(10, "VOID", 1e-16, 1e-16, 1e-16, 1e16, 1e16, ubuf,0);
744   
745   // --- Vacuum (with magnetic field) 
746   AliMaterial(11, "VOIM", 1e-16, 1e-16, 1e-16, 1e16, 1e16, ubuf,0);
747   
748   // --- Air (no magnetic field)
749   AliMaterial(12, "Air    $", 14.61, 7.3, .001205, 30420., 67500., ubuf, 0);
750   
751   // ---  Definition of tracking media: 
752   
753   // --- Tantalum = 1 ; 
754   // --- Brass = 2 ; 
755   // --- Fibers (SiO2) = 3 ; 
756   // --- Fibers (SiO2) = 4 ; 
757   // --- Iron = 5 ; 
758   // --- Lead = 6 ; 
759   // --- Vacuum (no field) = 10 
760   // --- Vacuum (with field) = 11 
761   // --- Air (no field) = 12 
762   
763   
764   // --- Tracking media parameters 
765   epsil  = .01;
766   stemax = 1.;
767   isvol  = 0;
768   isvolActive = 1;
769   inofld = 0;
770   fieldm = 0.;
771   
772   AliMedium(1, "ZTANT", 1, isvolActive, inofld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
773 //  AliMedium(1, "ZW", 1, isvolActive, inofld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
774   AliMedium(2, "ZBRASS", 2, isvolActive, inofld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
775   AliMedium(3, "ZSIO2", 3, isvolActive, inofld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
776   AliMedium(4, "ZQUAR", 3, isvolActive, inofld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
777 //  AliMedium(7, "ZCOPP", 7, isvolActive, inofld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
778 //  AliMedium(6, "ZLEAD", 6, isvolActive, inofld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
779   AliMedium(5, "ZIRON", 5, isvol, inofld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
780   AliMedium(10, "ZVOID", 10, isvol, inofld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
781   AliMedium(12, "ZAIR", 12, 0, inofld, fieldm, tmaxfd, stemax,deemax, epsil, stmin);
782   
783   fieldm = 45.;
784   AliMedium(11, "ZVOIM", 11, isvol, isxfld, fieldm, tmaxfd, stemax, deemax, epsil, stmin);
785   
786   // Thresholds for showering in the ZDCs 
787   
788   i = 1;
789   gMC->Gstpar(idtmed[i], "CUTGAM", .001);
790   gMC->Gstpar(idtmed[i], "CUTELE", .001);
791   gMC->Gstpar(idtmed[i], "CUTNEU", .01);
792   gMC->Gstpar(idtmed[i], "CUTHAD", .01);
793   i = 2;
794   gMC->Gstpar(idtmed[i], "CUTGAM", .001);
795   gMC->Gstpar(idtmed[i], "CUTELE", .001);
796   gMC->Gstpar(idtmed[i], "CUTNEU", .01);
797   gMC->Gstpar(idtmed[i], "CUTHAD", .01);
798   
799   // Avoid too detailed showering along the beam line 
800   
801   i = 5;
802   gMC->Gstpar(idtmed[i], "CUTGAM", .1);
803   gMC->Gstpar(idtmed[i], "CUTELE", .1);
804   gMC->Gstpar(idtmed[i], "CUTNEU", 1.);
805   gMC->Gstpar(idtmed[i], "CUTHAD", 1.);
806   
807   // Avoid interaction in fibers (only energy loss allowed) 
808   i = 3;
809   gMC->Gstpar(idtmed[i], "DCAY", 0.);
810   gMC->Gstpar(idtmed[i], "MULS", 0.);
811   gMC->Gstpar(idtmed[i], "PFIS", 0.);
812   gMC->Gstpar(idtmed[i], "MUNU", 0.);
813   gMC->Gstpar(idtmed[i], "LOSS", 1.);
814   gMC->Gstpar(idtmed[i], "PHOT", 0.);
815   gMC->Gstpar(idtmed[i], "COMP", 0.);
816   gMC->Gstpar(idtmed[i], "PAIR", 0.);
817   gMC->Gstpar(idtmed[i], "BREM", 0.);
818   gMC->Gstpar(idtmed[i], "DRAY", 0.);
819   gMC->Gstpar(idtmed[i], "ANNI", 0.);
820   gMC->Gstpar(idtmed[i], "HADR", 0.);
821   i = 4;
822   gMC->Gstpar(idtmed[i], "DCAY", 0.);
823   gMC->Gstpar(idtmed[i], "MULS", 0.);
824   gMC->Gstpar(idtmed[i], "PFIS", 0.);
825   gMC->Gstpar(idtmed[i], "MUNU", 0.);
826   gMC->Gstpar(idtmed[i], "LOSS", 1.);
827   gMC->Gstpar(idtmed[i], "PHOT", 0.);
828   gMC->Gstpar(idtmed[i], "COMP", 0.);
829   gMC->Gstpar(idtmed[i], "PAIR", 0.);
830   gMC->Gstpar(idtmed[i], "BREM", 0.);
831   gMC->Gstpar(idtmed[i], "DRAY", 0.);
832   gMC->Gstpar(idtmed[i], "ANNI", 0.);
833   gMC->Gstpar(idtmed[i], "HADR", 0.);
834   //
835   fMedSensF1 = idtmed[3];  // Sensitive volume: fibres type 1
836   fMedSensF2 = idtmed[4];  // Sensitive volume: fibres type 2
837   fMedSensZN = idtmed[1];  // Sensitive volume: ZN passive material
838   fMedSensZP = idtmed[2];  // Sensitive volume: ZP passive material
839   fMedSensGR = idtmed[12]; // Sensitive volume: air into the grooves
840
841
842 //_____________________________________________________________________________
843 void AliZDCv1::Init()
844 {
845  InitTables();
846
847 }
848
849 //_____________________________________________________________________________
850 void AliZDCv1::InitTables()
851 {
852   //Initialize parameters for light tables and read them
853   fNalfan = 90;
854   fNalfap = 90;
855   fNben = 18;
856   fNbep = 28;
857   
858   FILE *fp1, *fp2, *fp3, *fp4, *fp5, *fp6, *fp7, *fp8;
859
860   if((fp1 = fopen("light22620362207s","r")) == NULL){
861      printf("Cannot open file fp1 \n");
862      return;
863   }
864   if((fp2 = fopen("light22620362208s","r")) == NULL){
865      printf("Cannot open file fp2 \n");
866      return;
867   }  
868   if((fp3 = fopen("light22620362209s","r")) == NULL){
869      printf("Cannot open file fp3 \n");
870      return;
871   }
872   if((fp4 = fopen("light22620362210s","r")) == NULL){
873      printf("Cannot open file fp4 \n");
874      return;
875   }
876 //  printf(" --- Reading light tables for ZN \n");
877   for(int k=0; k<fNalfan; k++){
878      for(int j=0; j<fNben; j++){
879        fscanf(fp1,"%f",&fTablen[0][k][j]);
880        fscanf(fp2,"%f",&fTablen[1][k][j]);
881        fscanf(fp3,"%f",&fTablen[2][k][j]);
882        fscanf(fp4,"%f",&fTablen[3][k][j]);
883      } 
884   }
885   fclose(fp1);
886   fclose(fp2);
887   fclose(fp3);
888   fclose(fp4);
889   
890   if((fp5 = fopen("light22620552207s","r")) == NULL){
891      printf("Cannot open file fp5 \n");
892      return;
893   }
894   if((fp6 = fopen("light22620552208s","r")) == NULL){
895      printf("Cannot open file fp6 \n");
896      return;
897   }
898   if((fp7 = fopen("light22620552209s","r")) == NULL){
899      printf("Cannot open file fp7 \n");
900      return;
901   }
902   if((fp8 = fopen("light22620552210s","r")) == NULL){
903      printf("Cannot open file fp8 \n");
904      return;
905   }
906 //  printf(" --- Reading light tables for ZP \n");
907   for(int k=0; k<fNalfap; k++){
908      for(int j=0; j<fNbep; j++){
909        fscanf(fp5,"%f",&fTablep[0][k][j]);
910        fscanf(fp6,"%f",&fTablep[1][k][j]);
911        fscanf(fp7,"%f",&fTablep[2][k][j]);
912        fscanf(fp8,"%f",&fTablep[3][k][j]);
913      } 
914   }
915   fclose(fp5);
916   fclose(fp6);
917   fclose(fp7);
918   fclose(fp8);
919 }
920
921 //_____________________________________________________________________________
922 void AliZDCv1::StepManager()
923 {
924   //
925   // Routine called at every step in the Zero Degree Calorimeters
926   //
927
928   Int_t vol[2], ibeta, ialfa, ibe;
929   Float_t x[3], xdet[3], destep, hits[9], m, ekin, um[3], ud[3], be, radius, out;
930   TLorentzVector s, p;
931   const char *knamed;
932   
933   if((gMC->GetMedium() == fMedSensZN) || (gMC->GetMedium() == fMedSensZP) ||
934      (gMC->GetMedium() == fMedSensGR) || (gMC->GetMedium() == fMedSensF1) ||
935      (gMC->GetMedium() == fMedSensF2)){
936         
937   //Particle coordinates 
938     gMC->TrackPosition(s);
939     for(Int_t j=0; j<=2; j++){
940        x[j] = s[j];
941     }
942     hits[0] = x[0];
943     hits[1] = x[1];
944     hits[2] = x[2];
945
946   // Determine in which ZDC the particle is
947     knamed = gMC->CurrentVolName();
948     if(!strncmp(knamed,"ZN",2))vol[0]=1;
949     if(!strncmp(knamed,"ZP",2))vol[0]=2;
950   
951   // Determine in which quadrant the particle is
952     
953     //Quadrant in ZN
954     gMC->Gmtod(x,xdet,1);
955     if(vol[0]==1){
956       if((xdet[0]<0.) && (xdet[1]>0.)) vol[1]=1;
957       if((xdet[0]>0.) && (xdet[1]>0.)) vol[1]=2;
958       if((xdet[0]<0.) && (xdet[1]<0.)) vol[1]=3;
959       if((xdet[0]>0.) && (xdet[1]<0.)) vol[1]=4;
960     }
961     
962     //Quadrant in ZP
963     if(vol[0]==2){
964       Float_t xqZP = xdet[0]/(fDimZP[0]/2);
965       for(int i=1; i<=4; i++){
966          if(xqZP>(i-3) && xqZP<(i-2)){
967            vol[1] = i;
968            break;
969         }
970      }
971     }
972 //    printf("  -> Det. %d Quad. %d \n", vol[0], vol[1]);
973
974   // Store impact point and kinetic energy of the ENTERING particle
975     
976 //    Int_t Curtrack = gAlice->CurrentTrack();
977 //    Int_t Prim = gAlice->GetPrimary(Curtrack);
978 //    printf ("Primary: %d, Current Track: %d \n", Prim, Curtrack); 
979     
980 //    if(Curtrack==Prim){
981       if(gMC->IsTrackEntering()){
982         //Particle energy
983         gMC->TrackMomentum(p);
984 //       printf("p[0] = %f, p[1] = %f, p[2] = %f, p[3] = %f \n", 
985 //                 p[0], p[1], p[2], p[3]);
986         hits[3] = p[3];
987
988         // Impact point on ZN  
989         hits[4] = xdet[0];
990         hits[5] = xdet[1];
991         hits[7] = 0;
992         hits[8] = 0;
993         hits[9] = 0;
994
995 //        printf(" hits[2] = %f \n",hits[2]);
996         AddHit(gAlice->CurrentTrack(), vol, hits);
997       }
998 //    }
999              
1000       // Charged particles -> Energy loss
1001       if((destep=gMC->Edep())){
1002          if(gMC->IsTrackStop()){
1003            gMC->TrackMomentum(p);
1004            m = gMC->TrackMass();
1005            ekin = p[3]-m;
1006            if(ekin<0.) printf("ATTENTION!!!!!!!!!!!!!!! ->      ekin = %f <0 (?)",ekin);
1007            hits[9] = ekin;
1008            hits[7] = 0.;
1009            hits[8] = 0.;
1010            AddHit(gAlice->CurrentTrack(), vol, hits);
1011            }
1012          else{
1013            hits[9] = destep;
1014            hits[7] = 0.;
1015            hits[8] = 0.;
1016            AddHit(gAlice->CurrentTrack(), vol, hits);
1017            }
1018 //       printf("       -> Charged particle -> Dep. E = %f eV \n",hits[8]);
1019          }
1020 //       printf(" \n");
1021   }
1022
1023
1024   // *** Light production in fibres 
1025   if((gMC->GetMedium() == fMedSensF1) || (gMC->GetMedium() == fMedSensF2)){
1026 //    printf("%%%%%%%%%%%%%%%% Particle in fibre %%%%%%%%%%%%%%%%\n");
1027
1028      //Select charged particles
1029      if((destep=gMC->Edep())){
1030 //       printf("               -> CHARGED particle!!! \n");
1031
1032        // Particle velocity
1033        gMC->TrackMomentum(p);
1034        Float_t ptot=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
1035        Float_t beta =  ptot/p[3];
1036 //       printf("p[0] = %f, p[1] = %f, p[2] = %f, p[3] = %f, ptot = %f \n",
1037 //              p[0], p[1], p[2], p[3], ptot);
1038 //       Int_t pcID = gMC->TrackPid();
1039 //       printf("       Pc %d in quadrant %d -> beta = %f \n", pcID, vol[1], beta);
1040        if(beta<0.67) return;
1041        if((beta>=0.67) && (beta<=0.75)) ibeta = 0;
1042        if((beta>0.75)  && (beta<=0.85)) ibeta = 1;
1043        if((beta>0.85)  && (beta<=0.95)) ibeta = 2;
1044        if((beta>0.95)  && (beta<=1.00)) ibeta = 3;
1045  
1046        // Angle between particle trajectory and fibre axis
1047        // 1 -> Momentum directions
1048        um[0] = p[0]/ptot;
1049        um[1] = p[1]/ptot;
1050        um[2] = p[2]/ptot;
1051        gMC->Gmtod(um,ud,2);
1052        // 2 -> Angle < limit angle
1053        Double_t alfar = TMath::ACos(ud[2]);
1054        Double_t alfa = alfar*kRaddeg;
1055        if(alfa>110.) return;
1056        ialfa = Int_t(1.+alfa/2.);
1057  
1058        // Distance between particle trajectory and fibre axis
1059        gMC->TrackPosition(s);
1060        for(Int_t j=0; j<=2; j++){
1061           x[j] = s[j];
1062        }
1063        gMC->Gmtod(x,xdet,1);
1064        if(TMath::Abs(ud[0])>0.00001){
1065          Float_t dcoeff = ud[1]/ud[0];
1066          be = TMath::Abs((xdet[1]-dcoeff*xdet[0])/TMath::Sqrt(dcoeff*dcoeff+1.));
1067        }
1068        else{
1069          be = TMath::Abs(ud[0]);
1070        }
1071  
1072        if((vol[0]==1)) radius = fFibZN[1];
1073        if((vol[0]==2)) radius = fFibZP[1];
1074        ibe = Int_t(be*1000.+1);
1075  
1076        //Looking into the light tables 
1077        Float_t charge = gMC->TrackCharge();
1078        
1079        // (1)  ZN
1080        if((vol[0]==1)) {
1081          if(ibe>fNben) ibe=fNben;
1082          out =  charge*charge*fTablen[ibeta][ialfa][ibe];
1083 //       printf("       -> fTablen [%d][%d][%d] = %f \n", 
1084 //              ibeta, ialfa, ibe, fTablen[ibeta][ialfa][ibe]);
1085          if(gMC->GetMedium() == fMedSensF1){
1086            hits[7] = out;       //fLightPMQ
1087            hits[8] = 0;
1088            hits[9] = 0;
1089            AddHit(gAlice->CurrentTrack(), vol, hits);
1090          }
1091          else{
1092            hits[7] = 0;
1093            hits[8] = out;       //fLightPMC
1094            hits[9] = 0;
1095            AddHit(gAlice->CurrentTrack(), vol, hits);
1096          }
1097        } 
1098        
1099        // (2) ZP
1100        if((vol[0]==2)) {
1101          if(ibe>fNbep) ibe=fNbep;
1102          out =  charge*charge*fTablep[ibeta][ialfa][ibe];
1103 //       printf("       -> fTablep [%d][%d][%d] = %f\n", 
1104 //              ibeta, ialfa, ibe, fTablen[ibeta][ialfa][ibe]);
1105          if(gMC->GetMedium() == fMedSensF1){
1106            hits[7] = out;       //fLightPMQ
1107            hits[8] = 0;
1108            hits[9] = 0;
1109            AddHit(gAlice->CurrentTrack(), vol, hits);
1110          }
1111          else{
1112            hits[7] = 0;
1113            hits[8] = out;       //fLightPMC
1114            hits[9] = 0;
1115            AddHit(gAlice->CurrentTrack(), vol, hits);
1116          }
1117        } 
1118      }
1119 //    printf("\n");
1120    }
1121 }