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