LHAPDF veraion 5.9.1
[u/mrichter/AliRoot.git] / ZDC / AliZDCv4.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 ///////////////////////////////////////////////////////////////////////
18 //                                                                   //
19 //              AliZDCv4 --- new ZDC geometry                        //
20 //          with both ZDC arms geometry implemented                  //
21 //                                                                   //  
22 ///////////////////////////////////////////////////////////////////////
23
24 // --- Standard libraries
25 #include "stdio.h"
26
27 // --- ROOT system
28 #include <TMath.h>
29 #include <TRandom.h>
30 #include <TSystem.h>
31 #include <TTree.h>
32 #include <TVirtualMC.h>
33 #include <TGeoManager.h>
34 #include <TGeoMatrix.h>
35 #include <TGeoTube.h>
36 #include <TGeoCone.h>
37 #include <TGeoShape.h>
38 #include <TGeoScaledShape.h>
39 #include <TGeoCompositeShape.h>
40 #include <TParticle.h>
41
42 // --- AliRoot classes
43 #include "AliLog.h"
44 #include "AliConst.h"
45 #include "AliMagF.h"
46 #include "AliRun.h"
47 #include "AliZDCv4.h"
48 #include "AliMC.h"
49 #include "AliMCParticle.h"
50  
51 class  AliZDCHit;
52 class  AliPDG;
53 class  AliDetector;
54  
55  
56 ClassImp(AliZDCv4)
57
58 //_____________________________________________________________________________
59 AliZDCv4::AliZDCv4() : 
60   AliZDC(),
61   fMedSensF1(0),
62   fMedSensF2(0),
63   fMedSensZP(0),
64   fMedSensZN(0),
65   fMedSensZEM(0),
66   fMedSensGR(0),
67   fMedSensPI(0),
68   fMedSensTDI(0),
69   fMedSensVColl(0),
70   fMedSensLumi(0),
71   fNalfan(0),
72   fNalfap(0),
73   fNben(0),  
74   fNbep(0),
75   fZEMLength(0),
76   fpLostITC(0), 
77   fpLostD1C(0), 
78   fpcVCollC(0),
79   fpDetectedC(0),
80   fnDetectedC(0),
81   fpLostITA(0), 
82   fpLostD1A(0), 
83   fpLostTDI(0), 
84   fpcVCollA(0),
85   fpDetectedA(0),
86   fnDetectedA(0),
87   fVCollSideCAperture(7./2.),
88   fVCollSideCApertureNeg(7./2.),
89   fVCollSideCCentreY(0.),
90   fTCDDAperturePos(2.0),
91   fTCDDApertureNeg(2.2),
92   fTDIAperturePos(5.5),
93   fTDIApertureNeg(5.5),
94   fLumiLength(15.)
95 {
96   //
97   // Default constructor for Zero Degree Calorimeter
98   //
99   for(Int_t i=0; i<3; i++){
100      fDimZN[i] = fDimZP[i] = 0.;
101      fPosZNC[i] = fPosZNA[i] = fPosZPC[i]= fPosZPA[i] = fPosZEM[i] = 0.;
102      fFibZN[i] = fFibZP[i] = 0.;
103   }
104 }
105  
106 //_____________________________________________________________________________
107 AliZDCv4::AliZDCv4(const char *name, const char *title) : 
108   AliZDC(name,title),
109   fMedSensF1(0),
110   fMedSensF2(0),
111   fMedSensZP(0),
112   fMedSensZN(0),
113   fMedSensZEM(0),
114   fMedSensGR(0),
115   fMedSensPI(0),
116   fMedSensTDI(0),
117   fMedSensVColl(0),
118   fMedSensLumi(0),
119   fNalfan(90),
120   fNalfap(90),
121   fNben(18),  
122   fNbep(28), 
123   fZEMLength(0),
124   fpLostITC(0), 
125   fpLostD1C(0), 
126   fpcVCollC(0),
127   fpDetectedC(0),
128   fnDetectedC(0),
129   fpLostITA(0), 
130   fpLostD1A(0), 
131   fpLostTDI(0), 
132   fpcVCollA(0),
133   fpDetectedA(0),
134   fnDetectedA(0),
135   fVCollSideCAperture(7./2.),
136   fVCollSideCApertureNeg(7./2.),
137   fVCollSideCCentreY(0.),
138   fTCDDAperturePos(2.0),
139   fTCDDApertureNeg(2.2),
140   fTDIAperturePos(5.5),
141   fTDIApertureNeg(5.5),
142   fLumiLength(15.)  
143 {
144   //
145   // Standard constructor for Zero Degree Calorimeter 
146   //
147   //
148   // Check that DIPO, ABSO, DIPO and SHIL is there (otherwise tracking is wrong!!!)
149   
150   AliModule* pipe=gAlice->GetModule("PIPE");
151   AliModule* abso=gAlice->GetModule("ABSO");
152   AliModule* dipo=gAlice->GetModule("DIPO");
153   AliModule* shil=gAlice->GetModule("SHIL");
154   if((!pipe) || (!abso) || (!dipo) || (!shil)) {
155     Error("Constructor","ZDC needs PIPE, ABSO, DIPO and SHIL!!!\n");
156     exit(1);
157   } 
158   //
159   Int_t ip,jp,kp;
160   for(ip=0; ip<4; ip++){
161      for(kp=0; kp<fNalfap; kp++){
162         for(jp=0; jp<fNbep; jp++){
163            fTablep[ip][kp][jp] = 0;
164         } 
165      }
166   }
167   Int_t in,jn,kn;
168   for(in=0; in<4; in++){
169      for(kn=0; kn<fNalfan; kn++){
170         for(jn=0; jn<fNben; jn++){
171            fTablen[in][kn][jn] = 0;
172         } 
173      }
174   }
175   //
176   // Parameters for hadronic calorimeters geometry
177   // Positions updated after post-installation measurements
178   fDimZN[0] = 3.52;
179   fDimZN[1] = 3.52;
180   fDimZN[2] = 50.;  
181   fDimZP[0] = 11.2;
182   fDimZP[1] = 6.;
183   fDimZP[2] = 75.;    
184   fPosZNC[0] = 0.;
185   fPosZNC[1] = 0.;
186   fPosZNC[2] = -11397.3+136; 
187   fPosZPC[0] = 24.35;
188   fPosZPC[1] = 0.;
189   fPosZPC[2] = -11389.3+136; 
190   fPosZNA[0] = 0.;
191   fPosZNA[1] = 0.;
192   fPosZNA[2] = 11395.8-136;  
193   fPosZPA[0] = 24.35;
194   fPosZPA[1] = 0.;
195   fPosZPA[2] = 11387.8-136; 
196   fFibZN[0] = 0.;
197   fFibZN[1] = 0.01825;
198   fFibZN[2] = 50.;
199   fFibZP[0] = 0.;
200   fFibZP[1] = 0.0275;
201   fFibZP[2] = 75.;
202   // Parameters for EM calorimeter geometry
203   fPosZEM[0] = 8.5;
204   fPosZEM[1] = 0.;
205   fPosZEM[2] = 735.;
206   Float_t kDimZEMPb  = 0.15*(TMath::Sqrt(2.));  // z-dimension of the Pb slice
207   Float_t kDimZEMAir = 0.001;                   // scotch
208   Float_t kFibRadZEM = 0.0315;                  // External fiber radius (including cladding)
209   Int_t   kDivZEM[3] = {92, 0, 20};             // Divisions for EM detector
210   Float_t kDimZEM0 = 2*kDivZEM[2]*(kDimZEMPb+kDimZEMAir+kFibRadZEM*(TMath::Sqrt(2.)));
211   fZEMLength = kDimZEM0;
212   
213 }
214  
215 //_____________________________________________________________________________
216 void AliZDCv4::CreateGeometry()
217 {
218   //
219   // Create the geometry for the Zero Degree Calorimeter version 2
220   //* Initialize COMMON block ZDC_CGEOM
221   //*
222
223   CreateBeamLine();
224   CreateZDC();
225 }
226   
227 //_____________________________________________________________________________
228 void AliZDCv4::CreateBeamLine()
229 {
230   //
231   // Create the beam line elements
232   //
233   if(fOnlyZEM) printf("\n  Only ZEM configuration requested: no side-C beam pipe, no side-A hadronic ZDCs\n\n");
234   
235   Double_t zd1, zd2, zCorrDip, zInnTrip, zD1;
236   Double_t conpar[9], tubpar[3], tubspar[5], boxpar[3];
237
238   //-- rotation matrices for the legs
239   Int_t irotpipe1, irotpipe2;
240   TVirtualMC::GetMC()->Matrix(irotpipe1,90.-1.0027,0.,90.,90.,1.0027,180.);      
241   TVirtualMC::GetMC()->Matrix(irotpipe2,90.+1.0027,0.,90.,90.,1.0027,0.);
242
243   Int_t *idtmed = fIdtmed->GetArray();
244   Double_t dx=0., dy=0., dz=0.;
245   Double_t thx=0., thy=0., thz=0.;
246   Double_t phx=0., phy=0., phz=0.;
247   
248   TGeoMedium *medZDCFe = gGeoManager->GetMedium("ZDC_ZIRONT");
249   TGeoMedium *medZDCvoid = gGeoManager->GetMedium("ZDC_ZVOID");
250     
251   ////////////////////////////////////////////////////////////////
252   //                                                            //
253   //                SIDE C - RB26 (dimuon side)                 //
254   //                                                            //
255   ////////////////////////////////////////////////////////////////
256   
257 if(!fOnlyZEM){  
258   // -- Mother of the ZDCs (Vacuum PCON)
259   zd1 = 1921.6;
260   
261   conpar[0] = 0.;
262   conpar[1] = 360.;
263   conpar[2] = 2.;
264   conpar[3] = -13500.;
265   conpar[4] = 0.;
266   conpar[5] = 55.;
267   conpar[6] = -zd1;
268   conpar[7] = 0.;
269   conpar[8] = 55.;
270   TVirtualMC::GetMC()->Gsvolu("ZDCC", "PCON", idtmed[10], conpar, 9);
271   TVirtualMC::GetMC()->Gspos("ZDCC", 1, "ALIC", 0., 0., 0., 0, "ONLY");
272   
273
274   // -- BEAM PIPE from compensator dipole to the beginning of D1) 
275   tubpar[0] = 6.3/2.;
276   tubpar[1] = 6.7/2.;
277   // From beginning of ZDC volumes to beginning of D1
278   tubpar[2] = (5838.3-zd1)/2.;
279   TVirtualMC::GetMC()->Gsvolu("QT01", "TUBE", idtmed[7], tubpar, 3);
280   TVirtualMC::GetMC()->Gspos("QT01", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
281   // Ch.debug
282   //printf("  QT01 TUBE pipe from z = %1.2f to z = %1.2f (D1 begin)\n",-zd1,-2*tubpar[2]-zd1);
283   
284   //-- BEAM PIPE from the end of D1 to the beginning of D2) 
285   
286   //-- FROM MAGNETIC BEGINNING OF D1 TO MAGNETIC END OF D1
287   //--  Cylindrical pipe (r = 3.47) + conical flare  
288   // -> Beginning of D1
289   zd1 += 2.*tubpar[2];
290   
291   tubpar[0] = 6.94/2.;
292   tubpar[1] = 7.34/2.;
293   tubpar[2] = (6909.8-zd1)/2.;
294   TVirtualMC::GetMC()->Gsvolu("QT02", "TUBE", idtmed[7], tubpar, 3);
295   TVirtualMC::GetMC()->Gspos("QT02", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
296   // Ch.debug
297   //printf("    QT02 TUBE pipe from z = %1.2f to z = %1.2f (D1 magnetic end)\n",-zd1,-2*tubpar[2]-zd1);
298
299   zd1 += 2.*tubpar[2];
300   
301   tubpar[0] = 8./2.;
302   tubpar[1] = 8.6/2.;
303   tubpar[2] = (6958.3-zd1)/2.;
304   TVirtualMC::GetMC()->Gsvolu("QT0B", "TUBE", idtmed[7], tubpar, 3);
305   TVirtualMC::GetMC()->Gspos("QT0B", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
306   // Ch.debug
307   //printf("    QT0B TUBE pipe from z = %1.2f to z = %1.2f \n",-zd1,-2*tubpar[2]-zd1);
308  
309   zd1 += 2.*tubpar[2];
310   
311   tubpar[0] = 9./2.;
312   tubpar[1] = 9.6/2.;
313   tubpar[2] = (7022.8-zd1)/2.;
314   TVirtualMC::GetMC()->Gsvolu("QT03", "TUBE", idtmed[7], tubpar, 3);
315   TVirtualMC::GetMC()->Gspos("QT03", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
316   // Ch.debug
317   //printf("    QT03 TUBE pipe from z = %1.2f to z = %1.2f (D1 end)\n",-zd1,-2*tubpar[2]-zd1);
318
319   zd1 += 2.*tubpar[2];
320   
321   conpar[0] = 39.2/2.;
322   conpar[1] = 18./2.;
323   conpar[2] = 18.6/2.;
324   conpar[3] = 9./2.;
325   conpar[4] = 9.6/2.;
326   TVirtualMC::GetMC()->Gsvolu("QC01", "CONE", idtmed[7], conpar, 5);
327   TVirtualMC::GetMC()->Gspos("QC01", 1, "ZDCC", 0., 0., -conpar[0]-zd1, 0, "ONLY");
328   // Ch.debug
329   //printf("    QC01 CONE pipe from z = %1.2f to z= %1.2f (VCTCQ-I)\n",-zd1,-2*conpar[0]-zd1);
330   
331   zd1 += conpar[0] * 2.;
332   
333   // ******************************************************
334   // N.B.-> according to last vacuum layout 
335   // private communication by D. Macina, mail 27/1/2009
336   // updated to new ZDC installation (Janiary 2012) 
337   // ****************************************************** 
338   // 2nd section of    VCTCQ+VAMTF+TCLIA+VAMTF+1st part of VCTCP
339   Float_t totLength1 = 160.8 + 78. + 148. + 78. + 9.3;
340   //
341   tubpar[0] = 18.6/2.;
342   tubpar[1] = 7.6/2.;
343   tubpar[2] = totLength1/2.;
344 //  TVirtualMC::GetMC()->Gsvolu("QE01", "ELTU", idtmed[7], tubpar, 3);  
345   // temporary replace with a scaled tube (AG)
346   TGeoTube *tubeQE01 = new TGeoTube(0.,tubpar[0],tubpar[2]);
347   TGeoScale *scaleQE01 = new TGeoScale(1., tubpar[1]/tubpar[0], 1.);
348   TGeoScaledShape *sshapeQE01 = new TGeoScaledShape(tubeQE01, scaleQE01);
349   new TGeoVolume("QE01", sshapeQE01, gGeoManager->GetMedium(idtmed[7]));
350
351   tubpar[0] = 18.0/2.;
352   tubpar[1] = 7.0/2.;
353   tubpar[2] = totLength1/2.;
354 //  TVirtualMC::GetMC()->Gsvolu("QE02", "ELTU", idtmed[10], tubpar, 3);  
355   // temporary replace with a scaled tube (AG)
356   TGeoTube *tubeQE02 = new TGeoTube(0.,tubpar[0],tubpar[2]);
357   TGeoScale *scaleQE02 = new TGeoScale(1., tubpar[1]/tubpar[0], 1.);
358   TGeoScaledShape *sshapeQE02 = new TGeoScaledShape(tubeQE02, scaleQE02);
359   new TGeoVolume("QE02", sshapeQE02, gGeoManager->GetMedium(idtmed[10]));
360
361   TVirtualMC::GetMC()->Gspos("QE01", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY"); 
362   TVirtualMC::GetMC()->Gspos("QE02", 1, "QE01", 0., 0., 0., 0, "ONLY");  
363   // Ch.debug
364   //printf("    QE01 ELTU from z = %1.2f to z = %1.2f (VCTCQ-II+VAMTF+TCLIA+VAMTF+VCTCP-I)\n",-zd1,-2*tubpar[2]-zd1);
365   
366   // TCLIA collimator jaws (defined ONLY if fVCollAperture<3.5!)
367   if(fVCollSideCAperture<3.5){
368     boxpar[0] = 5.4/2.;
369     boxpar[1] = (3.5-fVCollSideCAperture-fVCollSideCCentreY-0.7)/2.;
370     if(boxpar[1]<0.) boxpar[1]=0.;
371     boxpar[2] = 124.4/2.;
372     printf("  AliZDCv4 -> C side injection collimator jaws: apertures +%1.2f/-%1.2f center %1.2f [cm]\n", 
373         fVCollSideCAperture, fVCollSideCApertureNeg,fVCollSideCCentreY);
374     TVirtualMC::GetMC()->Gsvolu("QCVC" , "BOX ", idtmed[13], boxpar, 3); 
375     TVirtualMC::GetMC()->Gspos("QCVC", 1, "QE02", -boxpar[0],  fVCollSideCAperture+fVCollSideCCentreY+boxpar[1], -totLength1/2.+160.8+78.+148./2., 0, "ONLY");  
376     TVirtualMC::GetMC()->Gspos("QCVC", 2, "QE02", -boxpar[0], -fVCollSideCApertureNeg+fVCollSideCCentreY-boxpar[1], -totLength1/2.+160.8+78.+148./2., 0, "ONLY");  
377   }
378   
379   zd1 += tubpar[2] * 2.;
380   
381   // 2nd part of VCTCP
382   conpar[0] = 31.5/2.;
383   conpar[1] = 21.27/2.;
384   conpar[2] = 21.87/2.;
385   conpar[3] = 18.0/2.;
386   conpar[4] = 18.6/2.;
387   TVirtualMC::GetMC()->Gsvolu("QC02", "CONE", idtmed[7], conpar, 5);
388   TVirtualMC::GetMC()->Gspos("QC02", 1, "ZDCC", 0., 0., -conpar[0]-zd1, 0, "ONLY");
389   // Ch.debug
390   //printf("    QC02 CONE pipe from z = %1.2f to z= %1.2f (VCTCP-II)\n",-zd1,-2*conpar[0]-zd1);
391   
392   zd1 += conpar[0] * 2.;
393
394   // 3rd section of VCTCP+VCDWC+VMLGB   
395   //Float_t totLenght2 = 9.2 + 530.5+40.;
396   Float_t totLenght2 = (8373.3-zd1);
397   tubpar[0] = 21.2/2.;
398   tubpar[1] = 21.9/2.;
399   tubpar[2] = totLenght2/2.;
400   TVirtualMC::GetMC()->Gsvolu("QT04", "TUBE", idtmed[7], tubpar, 3);
401   TVirtualMC::GetMC()->Gspos("QT04", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
402   // Ch.debug
403   //printf("    QT04 TUBE pipe from z = %1.2f to z= %1.2f (VCTCP-III)\n",-zd1,-2*tubpar[2]-zd1);
404   
405   zd1 += tubpar[2] * 2.;
406   
407   // First part of VCTCD
408   // skewed transition cone from ID=212.7 mm to ID=797 mm
409   conpar[0] = 121./2.;
410   conpar[1] = 79.7/2.;
411   conpar[2] = 81.3/2.;
412   conpar[3] = 21.27/2.;
413   conpar[4] = 21.87/2.;
414   TVirtualMC::GetMC()->Gsvolu("QC03", "CONE", idtmed[7], conpar, 5);
415   TVirtualMC::GetMC()->Gspos("QC03", 1, "ZDCC", 0., 0., -conpar[0]-zd1, 0, "ONLY");
416   // Ch.debug
417   //printf("    QC03 CONE pipe from z = %1.2f to z = %1.2f (VCTCD-I)\n",-zd1,-2*conpar[0]-zd1);
418   
419   zd1 += 2.*conpar[0];
420   
421   // VCDGB + 1st part of VCTCH
422   // Modified according to 2012 ZDC installation
423   tubpar[0] = 79.7/2.;
424   tubpar[1] = 81.3/2.;
425   tubpar[2] = (5*475.2+97.-136)/2.;
426   TVirtualMC::GetMC()->Gsvolu("QT05", "TUBE", idtmed[7], tubpar, 3);
427   TVirtualMC::GetMC()->Gspos("QT05", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
428   // Ch.debug
429   //printf("    QT05 TUBE pipe from z = %1.2f to z = %1.2f (VCDGB+VCTCH-I)\n",-zd1,-2*tubpar[2]-zd1);
430   
431   zd1 += 2.*tubpar[2];
432      
433   // 2nd part of VCTCH
434   // Transition from ID=797 mm to ID=196 mm:
435   // in order to simulate the thin window opened in the transition cone
436   // we divide the transition cone in three cones:
437   // (1) 8 mm thick (2) 3 mm thick (3) the third 8 mm thick
438   
439   // (1) 8 mm thick
440   conpar[0] = 9.09/2.; // 15 degree
441   conpar[1] = 74.82868/2.;
442   conpar[2] = 76.42868/2.; // thickness 8 mm 
443   conpar[3] = 79.7/2.;
444   conpar[4] = 81.3/2.; // thickness 8 mm  
445   TVirtualMC::GetMC()->Gsvolu("QC04", "CONE", idtmed[7], conpar, 5);
446   TVirtualMC::GetMC()->Gspos("QC04", 1, "ZDCC", 0., 0., -conpar[0]-zd1, 0, "ONLY");
447   // Ch.debug
448   //printf("    QC04 CONE pipe from z = %1.2f to z = %1.2f (VCTCH-II)\n",-zd1,-2*conpar[0]-zd1);
449
450   zd1 += 2.*conpar[0];  
451
452   // (2) 3 mm thick
453   conpar[0] = 96.2/2.; // 15 degree
454   conpar[1] = 23.19588/2.;
455   conpar[2] = 23.79588/2.; // thickness 3 mm 
456   conpar[3] = 74.82868/2.;
457   conpar[4] = 75.42868/2.; // thickness 3 mm  
458   TVirtualMC::GetMC()->Gsvolu("QC05", "CONE", idtmed[7], conpar, 5);
459   TVirtualMC::GetMC()->Gspos("QC05", 1, "ZDCC", 0., 0., -conpar[0]-zd1, 0, "ONLY");  
460   // Ch.debug
461   //printf("    QC05 CONE pipe from z = %1.2f to z = %1.2f (VCTCH-III)\n",-zd1,-2*conpar[0]-zd1);
462
463   zd1 += 2.*conpar[0];
464   
465   // (3) 8 mm thick
466   conpar[0] = 6.71/2.; // 15 degree
467   conpar[1] = 19.6/2.;
468   conpar[2] = 21.2/2.;// thickness 8 mm 
469   conpar[3] = 23.19588/2.;
470   conpar[4] = 24.79588/2.;// thickness 8 mm 
471   TVirtualMC::GetMC()->Gsvolu("QC06", "CONE", idtmed[7], conpar, 5);
472   TVirtualMC::GetMC()->Gspos("QC06", 1, "ZDCC", 0., 0., -conpar[0]-zd1, 0, "ONLY");
473   // Ch.debug
474   //printf("    QC06 CONE pipe from z = %1.2f to z = %1.2f (VCTCH-III)\n",-zd1,-2*conpar[0]-zd1);
475
476   zd1 += 2.*conpar[0];
477   
478   // VMZAR (5 volumes)  
479   tubpar[0] = 20.2/2.;
480   tubpar[1] = 20.6/2.;
481   tubpar[2] = 2.15/2.;
482   TVirtualMC::GetMC()->Gsvolu("QT06", "TUBE", idtmed[7], tubpar, 3);
483   TVirtualMC::GetMC()->Gspos("QT06", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
484   // Ch.debug
485   //printf("    QT06 TUBE pipe from z = %1.2f to z = %1.2f (VMZAR-I)\n",-zd1,-2*tubpar[2]-zd1);
486
487   zd1 += 2.*tubpar[2];
488   
489   conpar[0] = 6.9/2.;
490   conpar[1] = 23.9/2.;
491   conpar[2] = 24.3/2.;
492   conpar[3] = 20.2/2.;
493   conpar[4] = 20.6/2.;
494   TVirtualMC::GetMC()->Gsvolu("QC07", "CONE", idtmed[7], conpar, 5);
495   TVirtualMC::GetMC()->Gspos("QC07", 1, "ZDCC", 0., 0., -conpar[0]-zd1, 0, "ONLY");
496   // Ch.debug
497   //printf("    QC07 CONE pipe from z = %1.2f to z = %1.2f (VMZAR-II)\n",-zd1,-2*conpar[0]-zd1);
498
499   zd1 += 2.*conpar[0];
500
501   tubpar[0] = 23.9/2.;
502   tubpar[1] = 25.5/2.;
503   tubpar[2] = 17.0/2.;
504   TVirtualMC::GetMC()->Gsvolu("QT07", "TUBE", idtmed[7], tubpar, 3);
505   TVirtualMC::GetMC()->Gspos("QT07", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
506   // Ch.debug
507   //printf("    QT07 TUBE pipe from z = %1.2f to z = %1.2f (VMZAR-III)\n",-zd1,-2*tubpar[2]-zd1);
508  
509   zd1 += 2.*tubpar[2];
510   
511   conpar[0] = 6.9/2.;
512   conpar[1] = 20.2/2.;
513   conpar[2] = 20.6/2.;
514   conpar[3] = 23.9/2.;
515   conpar[4] = 24.3/2.;
516   TVirtualMC::GetMC()->Gsvolu("QC08", "CONE", idtmed[7], conpar, 5);
517   TVirtualMC::GetMC()->Gspos("QC08", 1, "ZDCC", 0., 0., -conpar[0]-zd1, 0, "ONLY");
518   // Ch.debug
519   //printf("    QC08 CONE pipe from z = %1.2f to z = %1.2f (VMZAR-IV)\n",-zd1,-2*conpar[0]-zd1);
520
521   zd1 += 2.*conpar[0];
522   
523   tubpar[0] = 20.2/2.;
524   tubpar[1] = 20.6/2.;
525   tubpar[2] = 2.15/2.;
526   TVirtualMC::GetMC()->Gsvolu("QT08", "TUBE", idtmed[7], tubpar, 3);
527   TVirtualMC::GetMC()->Gspos("QT08", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
528   // Ch.debug
529   //printf("    QT08 TUBE pipe from z = %1.2f to z = %1.2f (VMZAR-V)\n",-zd1,-2*tubpar[2]-zd1);
530
531   zd1 += 2.*tubpar[2];
532   
533   // Flange (ID=196 mm)(last part of VMZAR and first part of VCTYB)
534   tubpar[0] = 19.6/2.;
535   tubpar[1] = 25.3/2.;
536   tubpar[2] = 4.9/2.;
537   TVirtualMC::GetMC()->Gsvolu("QT09", "TUBE", idtmed[7], tubpar, 3);
538   TVirtualMC::GetMC()->Gspos("QT09", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
539   // Ch.debug
540   //printf("    QT09 TUBE pipe from z = %1.2f to z = %1.2f (VMZAR-VI+VCTYB-I)\n",-zd1,-2*tubpar[2]-zd1);
541  
542   zd1 += 2.*tubpar[2];
543   // Ch.debug
544   ////printf("  Beginning of VCTYB volume @ z = %1.2f \n",-zd1);
545   
546   // simulation of the trousers (VCTYB)     
547   tubpar[0] = 19.6/2.;
548   tubpar[1] = 20.0/2.;
549   tubpar[2] = 3.9/2.;
550   TVirtualMC::GetMC()->Gsvolu("QT10", "TUBE", idtmed[7], tubpar, 3);
551   TVirtualMC::GetMC()->Gspos("QT10", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
552   // Ch.debug
553   //printf("    QT10 TUBE pipe from z = %1.2f to z = %1.2f (VCTYB-II)\n",-zd1,-2*tubpar[2]-zd1);
554
555   zd1 += 2.*tubpar[2];
556
557   // transition cone from ID=196. to ID=216.6
558   conpar[0] = 32.55/2.;
559   conpar[1] = 21.66/2.;
560   conpar[2] = 22.06/2.;
561   conpar[3] = 19.6/2.;
562   conpar[4] = 20.0/2.;
563   TVirtualMC::GetMC()->Gsvolu("QC09", "CONE", idtmed[7], conpar, 5);
564   TVirtualMC::GetMC()->Gspos("QC09", 1, "ZDCC", 0., 0., -conpar[0]-zd1, 0, "ONLY");
565   // Ch.debug
566   //printf("    QC09 CONE pipe from z = %1.2f to z= %1.2f\n",-zd1,-2*conpar[0]-zd1);
567
568   zd1 += 2.*conpar[0]; 
569   
570   // tube  
571   tubpar[0] = 21.66/2.;
572   tubpar[1] = 22.06/2.;
573   tubpar[2] = 28.6/2.;
574   TVirtualMC::GetMC()->Gsvolu("QT11", "TUBE", idtmed[7], tubpar, 3);
575   TVirtualMC::GetMC()->Gspos("QT11", 1, "ZDCC", 0., 0., -tubpar[2]-zd1, 0, "ONLY");
576   // Ch.debug
577   //printf("    QT11 TUBE pipe from z = %1.2f to z= %1.2f\n",-zd1,-2*tubpar[2]-zd1);
578
579   zd1 += 2.*tubpar[2];
580   // Ch.debug
581   //printf("    Beginning of C side recombination chamber @ z = %f \n",-zd1);
582
583   // --------------------------------------------------------
584   // RECOMBINATION CHAMBER IMPLEMENTED USING TGeo CLASSES!!!!
585   // author: Chiara (August 2008)
586   // --------------------------------------------------------
587   // TRANSFORMATION MATRICES
588   // Combi transformation: 
589   dx = -3.970000;
590   dy = 0.000000;
591   dz = 0.0;
592   // Rotation: 
593   thx = 84.989100;   phx = 180.000000;
594   thy = 90.000000;   phy = 90.000000;
595   thz = 185.010900;  phz = 0.000000;
596   TGeoRotation *rotMatrix1c = new TGeoRotation("c",thx,phx,thy,phy,thz,phz);
597   // Combi transformation: 
598   dx = -3.970000;
599   dy = 0.000000;
600   dz = 0.0;
601   TGeoCombiTrans *rotMatrix2c = new TGeoCombiTrans("ZDCC_c1", dx,dy,dz,rotMatrix1c);
602   rotMatrix2c->RegisterYourself();
603   // Combi transformation: 
604   dx = 3.970000;
605   dy = 0.000000;
606   dz = 0.0;
607   // Rotation: 
608   thx = 95.010900;   phx = 180.000000;
609   thy = 90.000000;   phy = 90.000000;
610   thz = 180.-5.010900;    phz = 0.000000;
611   TGeoRotation *rotMatrix3c = new TGeoRotation("",thx,phx,thy,phy,thz,phz);
612   TGeoCombiTrans *rotMatrix4c = new TGeoCombiTrans("ZDCC_c2", dx,dy,dz,rotMatrix3c);
613   rotMatrix4c->RegisterYourself();
614
615   // VOLUMES DEFINITION
616   // Volume: ZDCC
617   TGeoVolume *pZDCC = gGeoManager->GetVolume("ZDCC");
618   
619   conpar[0] = (90.1-0.95-0.26-0.0085)/2.;
620   conpar[1] = 0.0/2.;
621   conpar[2] = 21.6/2.;
622   conpar[3] = 0.0/2.;
623   conpar[4] = 5.8/2.;
624   new TGeoCone("QCLext", conpar[0],conpar[1],conpar[2],conpar[3],conpar[4]);
625   
626   conpar[0] = (90.1-0.95-0.26-0.0085)/2.;
627   conpar[1] = 0.0/2.;
628   conpar[2] = 21.2/2.;
629   conpar[3] = 0.0/2.;
630   conpar[4] = 5.4/2.;
631   new TGeoCone("QCLint", conpar[0],conpar[1],conpar[2],conpar[3],conpar[4]);
632
633   // Outer trousers
634   TGeoCompositeShape *pOutTrousersC = new TGeoCompositeShape("outTrousersC", "QCLext:ZDCC_c1+QCLext:ZDCC_c2");
635   
636   // Volume: QCLext
637   TGeoVolume *pQCLext = new TGeoVolume("QCLext",pOutTrousersC, medZDCFe);
638   pQCLext->SetLineColor(kGreen);
639   pQCLext->SetVisLeaves(kTRUE);
640   //
641   TGeoTranslation *tr1c = new TGeoTranslation(0., 0., (Double_t) -conpar[0]-0.95-zd1);
642   //printf("    C side recombination chamber from z = %1.2f to z= %1.2f\n",-zd1,-2*conpar[0]-0.95-zd1);
643   //
644   pZDCC->AddNode(pQCLext, 1, tr1c);
645   // Inner trousers
646   TGeoCompositeShape *pIntTrousersC = new TGeoCompositeShape("intTrousersC", "QCLint:ZDCC_c1+QCLint:ZDCC_c2");
647   // Volume: QCLint
648   TGeoVolume *pQCLint = new TGeoVolume("QCLint",pIntTrousersC, medZDCvoid);
649   pQCLint->SetLineColor(kTeal);
650   pQCLint->SetVisLeaves(kTRUE);
651   pQCLext->AddNode(pQCLint, 1);
652     
653   zd1 += 90.1;
654   Double_t offset = 0.5;
655   zd1 = zd1+offset;
656   
657   //  second section : 2 tubes (ID = 54. OD = 58.)  
658   tubpar[0] = 5.4/2.;
659   tubpar[1] = 5.8/2.;
660   tubpar[2] = 40.0/2.;
661   TVirtualMC::GetMC()->Gsvolu("QT12", "TUBE", idtmed[7], tubpar, 3);
662   TVirtualMC::GetMC()->Gspos("QT12", 1, "ZDCC", -15.8/2., 0., -tubpar[2]-zd1, 0, "ONLY");
663   TVirtualMC::GetMC()->Gspos("QT12", 2, "ZDCC",  15.8/2., 0., -tubpar[2]-zd1, 0, "ONLY");  
664   // Ch.debug
665   //printf("    QT12 TUBE from z = %1.2f to z = %1.2f (separate beam pipes)\n",-zd1,-2*tubpar[2]-zd1);
666   
667   zd1 += 2.*tubpar[2];
668   
669   // transition x2zdc to recombination chamber : skewed cone  
670   conpar[0] = (10.-0.2-offset)/2.;
671   conpar[1] = 6.3/2.;
672   conpar[2] = 7.0/2.;
673   conpar[3] = 5.4/2.;
674   conpar[4] = 5.8/2.;
675   TVirtualMC::GetMC()->Gsvolu("QC10", "CONE", idtmed[7], conpar, 5); 
676   TVirtualMC::GetMC()->Gspos("QC10", 1, "ZDCC", -7.9-0.175, 0., -conpar[0]-0.1-zd1, irotpipe1, "ONLY");
677   TVirtualMC::GetMC()->Gspos("QC10", 2, "ZDCC", 7.9+0.175, 0., -conpar[0]-0.1-zd1, irotpipe2, "ONLY");
678   //printf("    QC10 CONE from z = %1.2f to z = %1.2f (transition X2ZDC)\n",-zd1,-2*conpar[0]-0.2-zd1);
679
680   zd1 += 2.*conpar[0]+0.2;
681   
682   // 2 tubes (ID = 63 mm OD=70 mm)      
683   tubpar[0] = 6.3/2.;
684   tubpar[1] = 7.0/2.;
685   tubpar[2] = 639.8/2.;
686   TVirtualMC::GetMC()->Gsvolu("QT13", "TUBE", idtmed[7], tubpar, 3);
687   TVirtualMC::GetMC()->Gspos("QT13", 1, "ZDCC", -16.5/2., 0., -tubpar[2]-zd1, 0, "ONLY");
688   TVirtualMC::GetMC()->Gspos("QT13", 2, "ZDCC",  16.5/2., 0., -tubpar[2]-zd1, 0, "ONLY");
689   //printf("    QT13 TUBE from z = %1.2f to z = %1.2f (separate beam pipes)\n",-zd1,-2*tubpar[2]-zd1);  
690
691   zd1 += 2.*tubpar[2];
692   printf("      END OF C SIDE BEAM PIPE DEFINITION @ z = %f m from IP2\n\n",-zd1/100.);
693
694            
695   // -- Luminometer (Cu box) in front of ZN - side C
696   if(fLumiLength>0.){
697     boxpar[0] = 8.0/2.;
698     boxpar[1] = 8.0/2.;
699     boxpar[2] = fLumiLength/2.;
700     TVirtualMC::GetMC()->Gsvolu("QLUC", "BOX ", idtmed[9], boxpar, 3);
701     TVirtualMC::GetMC()->Gspos("QLUC", 1, "ZDCC", 0., 0.,  fPosZNC[2]+66.+boxpar[2], 0, "ONLY");
702     printf("    C SIDE LUMINOMETER %1.2f < z < %1.2f\n",  fPosZNC[2]+66., fPosZNC[2]+66.+2*boxpar[2]);
703   }
704 }                
705   // --  END OF BEAM PIPE VOLUME DEFINITION FOR SIDE C (RB26 SIDE) 
706   // ----------------------------------------------------------------
707
708   ////////////////////////////////////////////////////////////////
709   //                                                            //
710   //                SIDE A - RB24                               //
711   //                                                            //
712   ///////////////////////////////////////////////////////////////
713
714   // Rotation Matrices definition
715   Int_t irotpipe3, irotpipe4, irotpipe5;
716   //-- rotation matrices for the tilted cone after the TDI to recenter vacuum chamber      
717   TVirtualMC::GetMC()->Matrix(irotpipe3,90.-1.8934,0.,90.,90.,1.8934,180.);    
718   //-- rotation matrices for the tilted tube before and after the TDI 
719   TVirtualMC::GetMC()->Matrix(irotpipe4,90.-3.8,0.,90.,90.,3.8,180.);       
720   //-- rotation matrix for the tilted cone after the TDI
721   TVirtualMC::GetMC()->Matrix(irotpipe5,90.+9.8,0.,90.,90.,9.8,0.);     
722
723   // -- Mother of the ZDCs (Vacuum PCON)                
724   zd2 = 1910.22;// zd2 initial value
725   
726   conpar[0] = 0.;
727   conpar[1] = 360.;
728   conpar[2] = 2.;
729   conpar[3] = zd2;
730   conpar[4] = 0.;
731   conpar[5] = 55.;
732   conpar[6] = 13500.;
733   conpar[7] = 0.;
734   conpar[8] = 55.;
735   TVirtualMC::GetMC()->Gsvolu("ZDCA", "PCON", idtmed[10], conpar, 9);
736   TVirtualMC::GetMC()->Gspos("ZDCA", 1, "ALIC", 0., 0., 0., 0, "ONLY");
737   
738   // To avoid overlaps 1 micron are left between certain volumes!
739   Double_t dxNoOverlap = 0.0;
740   //zd2 += dxNoOverlap;  
741   
742   // BEAM PIPE from 19.10 m to inner triplet beginning (22.965 m)  
743   tubpar[0] = 6.0/2.;
744   tubpar[1] = 6.4/2.;
745   tubpar[2] = 386.28/2. - dxNoOverlap; 
746   TVirtualMC::GetMC()->Gsvolu("QA01", "TUBE", idtmed[7], tubpar, 3);
747   TVirtualMC::GetMC()->Gspos("QA01", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
748   // Ch.debug
749   //printf("    QA01 TUBE centred in %f from z = %1.2f to z = %1.2f (IT begin)\n",tubpar[2]+zd2,zd2,2*tubpar[2]+zd2);
750   
751   zd2 += 2.*tubpar[2];  
752
753   // -- FIRST SECTION OF THE BEAM PIPE (from beginning of inner triplet to
754   //    beginning of D1)  
755   tubpar[0] = 6.3/2.;
756   tubpar[1] = 6.7/2.;
757   tubpar[2] = 3541.8/2. - dxNoOverlap;
758   TVirtualMC::GetMC()->Gsvolu("QA02", "TUBE", idtmed[7], tubpar, 3);
759   TVirtualMC::GetMC()->Gspos("QA02", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
760   // Ch.debug
761   //printf("    QA02 TUBE from z = %1.2f to z= %1.2f (D1 begin)\n",zd2,2*tubpar[2]+zd2);
762   
763   zd2 += 2.*tubpar[2]; 
764   
765     
766   // -- SECOND SECTION OF THE BEAM PIPE (from the beginning of D1 to the beginning of D2)
767   //
768   //  FROM (MAGNETIC) BEGINNING OF D1 TO THE (MAGNETIC) END OF D1 + 126.5 cm
769   //  CYLINDRICAL PIPE of diameter increasing from 6.75 cm up to 8.0 cm
770   //  from magnetic end :
771   //  1) 80.1 cm still with ID = 6.75 radial beam screen
772   //  2) 2.5 cm conical section from ID = 6.75 to ID = 8.0 cm
773   //  3) 43.9 cm straight section (tube) with ID = 8.0 cm
774
775   tubpar[0] = 6.75/2.;
776   tubpar[1] = 7.15/2.;
777   tubpar[2] = (945.0+80.1)/2.;
778   TVirtualMC::GetMC()->Gsvolu("QA03", "TUBE", idtmed[7], tubpar, 3);
779   TVirtualMC::GetMC()->Gspos("QA03", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
780   // Ch.debug
781   //printf("    QA03 TUBE from z = %1.2f to z = %1.2f (D1 end)\n",zd2,2*tubpar[2]+zd2);
782   
783   zd2 += 2.*tubpar[2];
784
785   // Transition Cone from ID=67.5 mm  to ID=80 mm
786   conpar[0] = 2.5/2.;
787   conpar[1] = 6.75/2.;
788   conpar[2] = 7.15/2.;
789   conpar[3] = 8.0/2.;
790   conpar[4] = 8.4/2.;
791   TVirtualMC::GetMC()->Gsvolu("QA04", "CONE", idtmed[7], conpar, 5);
792   TVirtualMC::GetMC()->Gspos("QA04", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
793   //printf("    QA04 CONE from z = %1.2f to z = %1.2f (transition cone)\n",zd2,2*conpar[0]+zd2);
794
795   zd2 += 2.*conpar[0];
796   
797   tubpar[0] = 8.0/2.;
798   tubpar[1] = 8.4/2.;
799   tubpar[2] = (43.9+20.+28.5+28.5)/2.;
800   TVirtualMC::GetMC()->Gsvolu("QA05", "TUBE", idtmed[7], tubpar, 3);
801   TVirtualMC::GetMC()->Gspos("QA05", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
802   // Ch.debug
803   //printf("    QA05 TUBE from z = %1.2f to z = %1.2f\n",zd2,2*tubpar[2]+zd2);
804   
805   zd2 += 2.*tubpar[2];
806
807   // Second section of VAEHI (transition cone from ID=80mm to ID=98mm)
808   conpar[0] = 4.0/2.;
809   conpar[1] = 8.0/2.;
810   conpar[2] = 8.4/2.;
811   conpar[3] = 9.8/2.;
812   conpar[4] = 10.2/2.;
813   TVirtualMC::GetMC()->Gsvolu("QAV1", "CONE", idtmed[7], conpar, 5);
814   TVirtualMC::GetMC()->Gspos("QAV1", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
815   //printf("    QAV1 CONE from z = %1.2f to z = %1.2f (VAEHI-I)\n",zd2,2*conpar[0]+zd2);
816
817   zd2 += 2.*conpar[0];
818   
819   //Third section of VAEHI (transition cone from ID=98mm to ID=90mm)
820   conpar[0] = 1.0/2.;
821   conpar[1] = 9.8/2.;
822   conpar[2] = 10.2/2.;
823   conpar[3] = 9.0/2.;
824   conpar[4] = 9.4/2.;
825   TVirtualMC::GetMC()->Gsvolu("QAV2", "CONE", idtmed[7], conpar, 5);
826   TVirtualMC::GetMC()->Gspos("QAV2", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
827   //printf("    QAV2 CONE from z = %1.2f to z = %1.2f (VAEHI-II)\n",zd2,2*conpar[0]+zd2);
828
829   zd2 += 2.*conpar[0];
830  
831   // Fourth section of VAEHI (tube ID=90mm)    
832   tubpar[0] = 9.0/2.;
833   tubpar[1] = 9.4/2.;
834   tubpar[2] = 31.0/2.;
835   TVirtualMC::GetMC()->Gsvolu("QAV3", "TUBE", idtmed[7], tubpar, 3);
836   TVirtualMC::GetMC()->Gspos("QAV3", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
837   // Ch.debug
838   //printf("    QAV3 TUBE from z = %1.2f to z = %1.2f (VAEHI-III)\n",zd2,2*tubpar[2]+zd2);
839   
840   zd2 += 2.*tubpar[2]; 
841
842   //---------------------------- TCDD beginning ----------------------------------    
843   // space for the insertion of the collimator TCDD (2 m)
844   // TCDD ZONE - 1st volume
845   conpar[0] = 1.3/2.;
846   conpar[1] = 9.0/2.;
847   conpar[2] = 13.0/2.;
848   conpar[3] = 9.6/2.;
849   conpar[4] = 13.0/2.;
850   TVirtualMC::GetMC()->Gsvolu("Q01T", "CONE", idtmed[7], conpar, 5);
851   TVirtualMC::GetMC()->Gspos("Q01T", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
852   //printf("    Q01T CONE from z = %1.2f to z = %1.2f (TCDD-I)\n",zd2,2*conpar[0]+zd2);
853
854   zd2 += 2.*conpar[0];  
855
856   // TCDD ZONE - 2nd volume    
857   tubpar[0] = 9.6/2.;
858   tubpar[1] = 10.0/2.;
859   tubpar[2] = 1.0/2.;
860   TVirtualMC::GetMC()->Gsvolu("Q02T", "TUBE", idtmed[7], tubpar, 3);
861   TVirtualMC::GetMC()->Gspos("Q02T", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
862   // Ch.debug
863   //printf("    Q02T TUBE from z = %1.2f to z= %1.2f (TCDD-II)\n",zd2,2*tubpar[2]+zd2);
864   
865   zd2 += 2.*tubpar[2]; 
866
867   // TCDD ZONE - third volume
868   conpar[0] = 9.04/2.;
869   conpar[1] = 9.6/2.;
870   conpar[2] = 10.0/2.;
871   conpar[3] = 13.8/2.;
872   conpar[4] = 14.2/2.;
873   TVirtualMC::GetMC()->Gsvolu("Q03T", "CONE", idtmed[7], conpar, 5);
874   TVirtualMC::GetMC()->Gspos("Q03T", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
875   //printf("    Q03T CONE from z = %1.2f to z= %1.2f (TCDD-III)\n",zd2,2*conpar[0]+zd2);
876
877   zd2 += 2.*conpar[0];  
878
879   // TCDD ZONE - 4th volume    
880   tubpar[0] = 13.8/2.;
881   tubpar[1] = 14.2/2.;
882   tubpar[2] = 38.6/2.;
883   TVirtualMC::GetMC()->Gsvolu("Q04T", "TUBE", idtmed[7], tubpar, 3);
884   TVirtualMC::GetMC()->Gspos("Q04T", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
885   // Ch.debug
886   //printf("    Q04T TUBE from z = %1.2f to z= %1.2f (TCDD-IV)\n",zd2,2*tubpar[2]+zd2);
887   
888   zd2 += 2.*tubpar[2]; 
889
890   // TCDD ZONE - 5th volume    
891   tubpar[0] = 21.0/2.;
892   tubpar[1] = 21.4/2.;
893   tubpar[2] = 100.12/2.;
894   TVirtualMC::GetMC()->Gsvolu("Q05T", "TUBE", idtmed[7], tubpar, 3);
895   TVirtualMC::GetMC()->Gspos("Q05T", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
896   // Ch.debug
897   //printf("    Q05T TUBE from z = %1.2f to z= %1.2f (TCDD-V)\n",zd2,2*tubpar[2]+zd2);
898
899   zd2 += 2.*tubpar[2]; 
900  
901   // TCDD ZONE - 6th volume    
902   tubpar[0] = 13.8/2.;
903   tubpar[1] = 14.2/2.;
904   tubpar[2] = 38.6/2.;
905   TVirtualMC::GetMC()->Gsvolu("Q06T", "TUBE", idtmed[7], tubpar, 3);
906   TVirtualMC::GetMC()->Gspos("Q06T", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
907   // Ch.debug
908   //printf("    Q06T TUBE from z = %1.2f to z= %1.2f (TCDD-VI)\n",zd2,2*tubpar[2]+zd2);
909   
910   zd2 += 2.*tubpar[2];
911
912   // TCDD ZONE - 7th volume
913   conpar[0] = 11.34/2.;
914   conpar[1] = 13.8/2.;
915   conpar[2] = 14.2/2.;
916   conpar[3] = 18.0/2.;
917   conpar[4] = 18.4/2.;
918   TVirtualMC::GetMC()->Gsvolu("Q07T", "CONE", idtmed[7], conpar, 5);
919   TVirtualMC::GetMC()->Gspos("Q07T", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
920   //printf("    Q07T CONE from z = %1.2f to z= %1.2f (TCDD-VII)\n",zd2,2*conpar[0]+zd2);
921
922   zd2 += 2.*conpar[0];
923
924   // Upper section : one single phi segment of a tube 
925   //  5 parameters for tubs: inner radius = 0.,
926   //    outer radius = 7. cm, half length = 50 cm
927   //    phi1 = 0., phi2 = 180. 
928   tubspar[0] = 0.0/2.;
929   tubspar[1] = 14.0/2.;
930   tubspar[2] = 100.0/2.;
931   tubspar[3] = 0.;
932   tubspar[4] = 180.;  
933   TVirtualMC::GetMC()->Gsvolu("Q08T", "TUBS", idtmed[7], tubspar, 5);
934   
935   // rectangular beam pipe inside TCDD upper section (Vacuum)  
936   boxpar[0] = 7.0/2.;
937   boxpar[1] = 2.2/2.;
938   boxpar[2] = 100./2.;
939   TVirtualMC::GetMC()->Gsvolu("Q09T", "BOX ", idtmed[10], boxpar, 3);
940   // positioning vacuum box in the upper section of TCDD
941   TVirtualMC::GetMC()->Gspos("Q09T", 1, "Q08T", 0., 1.1,  0., 0, "ONLY");
942   
943   // lower section : one single phi segment of a tube       
944   tubspar[0] = 0.0/2.;
945   tubspar[1] = 14.0/2.;
946   tubspar[2] = 100.0/2.;
947   tubspar[3] = 180.;
948   tubspar[4] = 360.;  
949   TVirtualMC::GetMC()->Gsvolu("Q10T", "TUBS", idtmed[7], tubspar, 5);
950   // rectangular beam pipe inside TCDD lower section (Vacuum)  
951   boxpar[0] = 7.0/2.;
952   boxpar[1] = 2.2/2.;
953   boxpar[2] = 100./2.;
954   TVirtualMC::GetMC()->Gsvolu("Q11T", "BOX ", idtmed[10], boxpar, 3);
955   // positioning vacuum box in the lower section of TCDD
956   TVirtualMC::GetMC()->Gspos("Q11T", 1, "Q10T", 0., -1.1,  0., 0, "ONLY");  
957   
958   // positioning  TCDD elements in ZDCA, (inside TCDD volume)
959   TVirtualMC::GetMC()->Gspos("Q08T", 1, "ZDCA", 0., fTCDDAperturePos, -100.+zd2, 0, "ONLY");  
960   TVirtualMC::GetMC()->Gspos("Q10T", 1, "ZDCA", 0., -fTCDDApertureNeg, -100.+zd2, 0, "ONLY");  
961   printf("  AliZDCv4 -> TCDD apertures +%1.2f/-%1.2f cm\n", 
962         fTCDDAperturePos, fTCDDApertureNeg);
963     
964   // RF screen 
965   boxpar[0] = 0.2/2.;
966   boxpar[1] = 4.0/2.;
967   boxpar[2] = 100./2.;
968   TVirtualMC::GetMC()->Gsvolu("Q12T", "BOX ", idtmed[7], boxpar, 3);  
969   // positioning RF screen at both sides of TCDD
970   TVirtualMC::GetMC()->Gspos("Q12T", 1, "ZDCA", tubspar[1]+boxpar[0], 0., -100.+zd2, 0, "ONLY");  
971   TVirtualMC::GetMC()->Gspos("Q12T", 2, "ZDCA", -tubspar[1]-boxpar[0], 0., -100.+zd2, 0, "ONLY");      
972   //---------------------------- TCDD end ---------------------------------------    
973
974   // The following elliptical tube 180 mm x 70 mm
975   // (obtained positioning the void QA06 in QA07)
976   // represents VAMTF + first part of VCTCP (93 mm)
977   // updated according to 2012 new ZDC installation
978
979   tubpar[0] = 18.4/2.;
980   tubpar[1] = 7.4/2.;
981   tubpar[2] = (78+9.3)/2.;
982 //  TVirtualMC::GetMC()->Gsvolu("QA06", "ELTU", idtmed[7], tubpar, 3);  
983   // temporary replace with a scaled tube (AG)
984   TGeoTube *tubeQA06 = new TGeoTube(0.,tubpar[0],tubpar[2]);
985   TGeoScale *scaleQA06 = new TGeoScale(1., tubpar[1]/tubpar[0], 1.);
986   TGeoScaledShape *sshapeQA06 = new TGeoScaledShape(tubeQA06, scaleQA06);
987   new TGeoVolume("QA06", sshapeQA06, gGeoManager->GetMedium(idtmed[7]));
988   //printf("    QA06 TUBE from z = %1.2f to z = %1.2f (VAMTF+VCTCP-I)\n",zd2,2*tubpar[2]+zd2);
989
990   tubpar[0] = 18.0/2.;
991   tubpar[1] = 7.0/2.;
992   tubpar[2] = (78+9.3)/2.;
993 //  TVirtualMC::GetMC()->Gsvolu("QA07", "ELTU", idtmed[10], tubpar, 3);  
994   // temporary replace with a scaled tube (AG)
995   TGeoTube *tubeQA07 = new TGeoTube(0.,tubpar[0],tubpar[2]);
996   TGeoScale *scaleQA07 = new TGeoScale(1., tubpar[1]/tubpar[0], 1.);
997   TGeoScaledShape *sshapeQA07 = new TGeoScaledShape(tubeQA07, scaleQA07);
998   new TGeoVolume("QA07", sshapeQA07, gGeoManager->GetMedium(idtmed[10]));
999   ////printf("  QA07 TUBE from z = %1.2f to z= %1.2f\n",zd2,2*tubpar[2]+zd2);
1000   TVirtualMC::GetMC()->Gspos("QA06", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY"); 
1001   TVirtualMC::GetMC()->Gspos("QA07", 1, "QA06", 0., 0., 0., 0, "ONLY");  
1002     
1003   zd2 += 2.*tubpar[2];
1004       
1005   // VCTCP second part: transition cone from ID=180 to ID=212.7 
1006   conpar[0] = 31.5/2.;
1007   conpar[1] = 18.0/2.;
1008   conpar[2] = 18.6/2.;
1009   conpar[3] = 21.27/2.;
1010   conpar[4] = 21.87/2.;
1011   TVirtualMC::GetMC()->Gsvolu("QA08", "CONE", idtmed[7], conpar, 5);
1012   TVirtualMC::GetMC()->Gspos("QA08", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
1013   // Ch.debug  
1014   //printf("    QA08 CONE from z = %f to z = %f (VCTCP-II)\n",zd2,2*conpar[0]+zd2);
1015
1016   zd2 += 2.*conpar[0];
1017   
1018   // Tube ID 212.7 mm
1019   // Represents VCTCP third part (92 mm) + VCDWB (765 mm) + VMBGA (400 mm) +
1020   //            VCDWE (300 mm) + VMBGA (400 mm)
1021   // + TCTVB space + VAMTF space (new installation Jan 2012)
1022   tubpar[0] = 21.27/2.;
1023   tubpar[1] = 21.87/2.;
1024   tubpar[2] = (195.7+148.+78.)/2.;
1025   TVirtualMC::GetMC()->Gsvolu("QA09", "TUBE", idtmed[7], tubpar, 3);
1026   TVirtualMC::GetMC()->Gspos("QA09", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1027   //printf("    QA09 TUBE from z = %1.2f to z= %1.2f (VCTCP-III+VCDWB+VMBGA+VCDWE+VMBGA)\n",zd2,2*tubpar[2]+zd2);
1028
1029   zd2 += 2.*tubpar[2];
1030
1031   // skewed transition piece (ID=212.7 mm to 332 mm) (before TDI)   
1032   conpar[0] = (50.0-0.73-1.13)/2.;
1033   conpar[1] = 21.27/2.;
1034   conpar[2] = 21.87/2.;
1035   conpar[3] = 33.2/2.;
1036   conpar[4] = 33.8/2.;
1037   TVirtualMC::GetMC()->Gsvolu("QA10", "CONE", idtmed[7], conpar, 5);
1038   TVirtualMC::GetMC()->Gspos("QA10", 1, "ZDCA", -1.66, 0., conpar[0]+0.73+zd2, irotpipe4, "ONLY");
1039   // Ch.debug  
1040   //printf("    QA10 skewed CONE from z = %1.2f to z= %1.2f\n",zd2,2*conpar[0]+0.73+1.13+zd2);
1041
1042   zd2 += 2.*conpar[0]+0.73+1.13;
1043       
1044   // Vacuum chamber containing TDI  
1045   tubpar[0] = 0.;
1046   tubpar[1] = 54.6/2.;
1047   tubpar[2] = 540.0/2.;
1048   TVirtualMC::GetMC()->Gsvolu("Q13TM", "TUBE", idtmed[10], tubpar, 3);
1049   TVirtualMC::GetMC()->Gspos("Q13TM", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1050   tubpar[0] = 54.0/2.;
1051   tubpar[1] = 54.6/2.;
1052   tubpar[2] = 540.0/2.;
1053   TVirtualMC::GetMC()->Gsvolu("Q13T", "TUBE", idtmed[7], tubpar, 3);
1054   TVirtualMC::GetMC()->Gspos("Q13T", 1, "Q13TM", 0., 0., 0., 0, "ONLY");
1055   // Ch.debug
1056   //printf("    Q13T TUBE from z = %1.2f to z= %1.2f (TDI vacuum chamber)\n",zd2,2*tubpar[2]+zd2);
1057
1058   zd2 += 2.*tubpar[2];
1059   
1060   //---------------- INSERT TDI INSIDE Q13T -----------------------------------    
1061   boxpar[0] = 11.0/2.;
1062   boxpar[1] = 9.0/2.;
1063   boxpar[2] = 540.0/2.;
1064   TVirtualMC::GetMC()->Gsvolu("QTD1", "BOX ", idtmed[7], boxpar, 3);
1065   TVirtualMC::GetMC()->Gspos("QTD1", 1, "Q13TM", -3.8, boxpar[1]+fTDIAperturePos,  0., 0, "ONLY");
1066   boxpar[0] = 11.0/2.;
1067   boxpar[1] = 9.0/2.;
1068   boxpar[2] = 540.0/2.;
1069   TVirtualMC::GetMC()->Gsvolu("QTD2", "BOX ", idtmed[7], boxpar, 3);
1070   TVirtualMC::GetMC()->Gspos("QTD2", 1, "Q13TM", -3.8, -boxpar[1]-fTDIApertureNeg,  0., 0, "ONLY");  
1071   boxpar[0] = 5.1/2.;
1072   boxpar[1] = 0.2/2.;
1073   boxpar[2] = 540.0/2.;
1074   TVirtualMC::GetMC()->Gsvolu("QTD3", "BOX ", idtmed[7], boxpar, 3);
1075   TVirtualMC::GetMC()->Gspos("QTD3", 1, "Q13TM", -3.8+5.5+boxpar[0], fTDIAperturePos,  0., 0, "ONLY");  
1076   TVirtualMC::GetMC()->Gspos("QTD3", 2, "Q13TM", -3.8+5.5+boxpar[0], -fTDIApertureNeg,  0., 0, "ONLY"); 
1077   TVirtualMC::GetMC()->Gspos("QTD3", 3, "Q13TM", -3.8-5.5-boxpar[0], fTDIAperturePos,  0., 0, "ONLY");  
1078   TVirtualMC::GetMC()->Gspos("QTD3", 4, "Q13TM", -3.8-5.5-boxpar[0], -fTDIApertureNeg,  0., 0, "ONLY");  
1079   printf("  AliZDCv4 -> TDI apertures +%1.2f/-%1.2f cm\n", 
1080         fTDIAperturePos, fTDIApertureNeg);
1081   //
1082   tubspar[0] = 12.0/2.;
1083   tubspar[1] = 12.4/2.;
1084   tubspar[2] = 540.0/2.;
1085   tubspar[3] = 90.;
1086   tubspar[4] = 270.;  
1087   TVirtualMC::GetMC()->Gsvolu("QTD4", "TUBS", idtmed[7], tubspar, 5);
1088   TVirtualMC::GetMC()->Gspos("QTD4", 1, "Q13TM", -3.8-10.6, 0.,  0., 0, "ONLY");
1089   tubspar[0] = 12.0/2.;
1090   tubspar[1] = 12.4/2.;
1091   tubspar[2] = 540.0/2.;
1092   tubspar[3] = -90.;
1093   tubspar[4] = 90.;  
1094   TVirtualMC::GetMC()->Gsvolu("QTD5", "TUBS", idtmed[7], tubspar, 5);
1095   TVirtualMC::GetMC()->Gspos("QTD5", 1, "Q13TM", -3.8+10.6, 0.,  0., 0, "ONLY"); 
1096   //---------------- END DEFINING TDI INSIDE Q13T -------------------------------
1097   
1098   // VCTCG skewed transition piece (ID=332 mm to 212.7 mm) (after TDI)
1099   conpar[0] = (50.0-2.92-1.89)/2.;
1100   conpar[1] = 33.2/2.;
1101   conpar[2] = 33.8/2.;
1102   conpar[3] = 21.27/2.;
1103   conpar[4] = 21.87/2.;
1104   TVirtualMC::GetMC()->Gsvolu("QA11", "CONE", idtmed[7], conpar, 5);
1105   TVirtualMC::GetMC()->Gspos("QA11", 1, "ZDCA", 4.32-3.8, 0., conpar[0]+2.92+zd2, irotpipe5, "ONLY");
1106   // Ch.debug  
1107   //printf("    QA11 skewed CONE from z = %f to z =%f (VCTCG)\n",zd2,2*conpar[0]+2.92+1.89+zd2);
1108
1109   zd2 += 2.*conpar[0]+2.92+1.89;
1110   
1111   // The following tube ID 212.7 mm  
1112   // represents VMBGA (400 mm) + VCDWE (300 mm) + VMBGA (400 mm) +
1113   //            BTVTS (600 mm) + VMLGB (400 mm)  
1114   tubpar[0] = 21.27/2.;
1115   tubpar[1] = 21.87/2.;
1116   tubpar[2] = 210.0/2.;
1117   TVirtualMC::GetMC()->Gsvolu("QA12", "TUBE", idtmed[7], tubpar, 3);
1118   TVirtualMC::GetMC()->Gspos("QA12", 1, "ZDCA", 4., 0., tubpar[2]+zd2, 0, "ONLY");
1119   // Ch.debug
1120   //printf("    QA12 TUBE from z = %1.2f to z= %1.2f (VMBGA+VCDWE+VMBGA+BTVTS+VMLGB)\n",zd2,2*tubpar[2]+zd2);
1121
1122   zd2 += 2.*tubpar[2];  
1123   
1124   // First part of VCTCC
1125   // skewed transition cone from ID=212.7 mm to ID=797 mm
1126   conpar[0] = (121.0-0.37-1.35)/2.;
1127   conpar[1] = 21.27/2.;
1128   conpar[2] = 21.87/2.;
1129   conpar[3] = 79.7/2.;
1130   conpar[4] = 81.3/2.;
1131   TVirtualMC::GetMC()->Gsvolu("QA13", "CONE", idtmed[7], conpar, 5);
1132   TVirtualMC::GetMC()->Gspos("QA13", 1, "ZDCA", 4.-2., 0., conpar[0]+0.37+zd2, irotpipe3, "ONLY");
1133   // Ch.debug  
1134   //printf("    QA13 CONE from z = %1.2f to z = %1.2f (VCTCC-I)\n",zd2,2*conpar[0]+0.37+1.35+zd2);
1135
1136   zd2 += 2.*conpar[0]+0.37+1.35;
1137   
1138   // The following tube ID 797 mm  
1139   // represents the second part of VCTCC (4272 mm) + 
1140   //            4 x VCDGA (4 x 4272 mm) + 
1141   //            the first part of VCTCR (850 mm)
1142   // updated according to 2012 ZDC installation
1143   tubpar[0] = 79.7/2.;
1144   tubpar[1] = 81.3/2.;
1145   tubpar[2] = (2221.-136.)/2.;
1146   TVirtualMC::GetMC()->Gsvolu("QA14", "TUBE", idtmed[7], tubpar, 3);
1147   TVirtualMC::GetMC()->Gspos("QA14", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1148   // Ch.debug  
1149   //printf("    QA14 TUBE from z = %1.2f to z = %1.2f (VCTCC-II)\n",zd2,2*tubpar[2]+zd2);
1150
1151   zd2 += 2.*tubpar[2];
1152         
1153   // Second part of VCTCR
1154   // Transition from ID=797 mm to ID=196 mm:
1155   // in order to simulate the thin window opened in the transition cone
1156   // we divide the transition cone in three cones:
1157   // (1) 8 mm thick (2) 3 mm thick (3) the third 8 mm thick
1158   
1159   // (1) 8 mm thick
1160   conpar[0] = 9.09/2.; // 15 degree
1161   conpar[1] = 79.7/2.;
1162   conpar[2] = 81.3/2.; // thickness 8 mm  
1163   conpar[3] = 74.82868/2.;
1164   conpar[4] = 76.42868/2.; // thickness 8 mm 
1165   TVirtualMC::GetMC()->Gsvolu("QA15", "CONE", idtmed[7], conpar, 5);
1166   TVirtualMC::GetMC()->Gspos("QA15", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
1167   //printf("    QA15 CONE from z = %1.2f to z= %1.2f (VCTCR-I)\n",zd2,2*conpar[0]+zd2);
1168
1169   zd2 += 2.*conpar[0];  
1170
1171   // (2) 3 mm thick
1172   conpar[0] = 96.2/2.; // 15 degree
1173   conpar[1] = 74.82868/2.;
1174   conpar[2] = 75.42868/2.; // thickness 3 mm  
1175   conpar[3] = 23.19588/2.;
1176   conpar[4] = 23.79588/2.; // thickness 3 mm 
1177   TVirtualMC::GetMC()->Gsvolu("QA16", "CONE", idtmed[7], conpar, 5);
1178   TVirtualMC::GetMC()->Gspos("QA16", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");  
1179   //printf("    QA16 CONE from z = %1.2f to z= %1.2f\n",zd2,2*conpar[0]+zd2);
1180
1181   zd2 += 2.*conpar[0];
1182   
1183   // (3) 8 mm thick
1184   conpar[0] = 6.71/2.; // 15 degree
1185   conpar[1] = 23.19588/2.;
1186   conpar[2] = 24.79588/2.;// thickness 8 mm 
1187   conpar[3] = 19.6/2.;
1188   conpar[4] = 21.2/2.;// thickness 8 mm 
1189   TVirtualMC::GetMC()->Gsvolu("QA17", "CONE", idtmed[7], conpar, 5);
1190   TVirtualMC::GetMC()->Gspos("QA17", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
1191   //printf("    QA17 CONE from z = %1.2f to z= %1.2f (VCTCR-II)\n",zd2,2*conpar[0]+zd2);
1192
1193   zd2 += 2.*conpar[0];
1194  
1195   // Third part of VCTCR: tube (ID=196 mm)  
1196   tubpar[0] = 19.6/2.;
1197   tubpar[1] = 21.2/2.;
1198   tubpar[2] = 9.55/2.;
1199   TVirtualMC::GetMC()->Gsvolu("QA18", "TUBE", idtmed[7], tubpar, 3);
1200   TVirtualMC::GetMC()->Gspos("QA18", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1201   // Ch.debug  
1202   //printf("    QA18 TUBE from z = %1.2f to z= %1.2f (VCTCR-III)\n",zd2,2*tubpar[2]+zd2);
1203
1204   zd2 += 2.*tubpar[2];  
1205   
1206   // Flange (ID=196 mm) (last part of VCTCR and first part of VMZAR)
1207   tubpar[0] = 19.6/2.;
1208   tubpar[1] = 25.3/2.;
1209   tubpar[2] = 4.9/2.;
1210   TVirtualMC::GetMC()->Gsvolu("QF01", "TUBE", idtmed[7], tubpar, 3);
1211   TVirtualMC::GetMC()->Gspos("QF01", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1212   // Ch.debug  
1213   //printf("    QF01  TUBE from z = %1.2f to z= %1.2f (VMZAR-I)\n",zd2,2*tubpar[2]+zd2);
1214
1215   zd2 += 2.*tubpar[2];
1216   
1217   // VMZAR (5 volumes)  
1218   tubpar[0] = 20.2/2.;
1219   tubpar[1] = 20.6/2.;
1220   tubpar[2] = 2.15/2.;
1221   TVirtualMC::GetMC()->Gsvolu("QA19", "TUBE", idtmed[7], tubpar, 3);
1222   TVirtualMC::GetMC()->Gspos("QA19", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1223   // Ch.debug  
1224   //printf("    QA19  TUBE from z = %1.2f to z = %1.2f (VMZAR-II)\n",zd2,2*tubpar[2]+zd2);
1225
1226   zd2 += 2.*tubpar[2];
1227   
1228   conpar[0] = 6.9/2.;
1229   conpar[1] = 20.2/2.;
1230   conpar[2] = 20.6/2.;
1231   conpar[3] = 23.9/2.;
1232   conpar[4] = 24.3/2.;
1233   TVirtualMC::GetMC()->Gsvolu("QA20", "CONE", idtmed[7], conpar, 5);
1234   TVirtualMC::GetMC()->Gspos("QA20", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
1235   // Ch.debug  
1236   //printf("    QA20 CONE from z = %1.2f to z = %1.2f (VMZAR-III)\n",zd2,2*conpar[0]+zd2);
1237
1238   zd2 += 2.*conpar[0];
1239
1240   tubpar[0] = 23.9/2.;
1241   tubpar[1] = 25.5/2.;
1242   tubpar[2] = 17.0/2.;
1243   TVirtualMC::GetMC()->Gsvolu("QA21", "TUBE", idtmed[7], tubpar, 3);
1244   TVirtualMC::GetMC()->Gspos("QA21", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1245   // Ch.debug  
1246   //printf("    QA21  TUBE from z = %1.2f to z = %1.2f  (VMZAR-IV)\n",zd2,2*tubpar[2]+zd2);
1247
1248   zd2 += 2.*tubpar[2];
1249   
1250   conpar[0] = 6.9/2.;
1251   conpar[1] = 23.9/2.;
1252   conpar[2] = 24.3/2.;
1253   conpar[3] = 20.2/2.;
1254   conpar[4] = 20.6/2.;
1255   TVirtualMC::GetMC()->Gsvolu("QA22", "CONE", idtmed[7], conpar, 5);
1256   TVirtualMC::GetMC()->Gspos("QA22", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
1257   // Ch.debug  
1258   //printf("    QA22 CONE from z = %1.2f to z = %1.2f (VMZAR-V)\n",zd2,2*conpar[0]+zd2);
1259
1260   zd2 += 2.*conpar[0];
1261   
1262   tubpar[0] = 20.2/2.;
1263   tubpar[1] = 20.6/2.;
1264   tubpar[2] = 2.15/2.;
1265   TVirtualMC::GetMC()->Gsvolu("QA23", "TUBE", idtmed[7], tubpar, 3);
1266   TVirtualMC::GetMC()->Gspos("QA23", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1267   // Ch.debug  
1268   //printf("    QA23  TUBE from z = %1.2f to z= %1.2f (VMZAR-VI)\n",zd2,2*tubpar[2]+zd2);
1269
1270   zd2 += 2.*tubpar[2];
1271   
1272   // Flange (ID=196 mm)(last part of VMZAR and first part of VCTYD)
1273   tubpar[0] = 19.6/2.;
1274   tubpar[1] = 25.3/2.;
1275   tubpar[2] = 4.9/2.;
1276   TVirtualMC::GetMC()->Gsvolu("QF02", "TUBE", idtmed[7], tubpar, 3);
1277   TVirtualMC::GetMC()->Gspos("QF02", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1278   // Ch.debug  
1279   //printf("    QF02 TUBE from z = %1.2f to z= %1.2f (VMZAR-VII)\n",zd2,2*tubpar[2]+zd2);
1280
1281   zd2 += 2.*tubpar[2];
1282   
1283   // simulation of the trousers (VCTYB)     
1284   tubpar[0] = 19.6/2.;
1285   tubpar[1] = 20.0/2.;
1286   tubpar[2] = 3.9/2.;
1287   TVirtualMC::GetMC()->Gsvolu("QA24", "TUBE", idtmed[7], tubpar, 3);
1288   TVirtualMC::GetMC()->Gspos("QA24", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1289   // Ch.debug
1290   //printf("    QA24  TUBE from z = %1.2f to z= %1.2f (VCTYB)\n",zd2,2*tubpar[2]+zd2);
1291
1292   zd2 += 2.*tubpar[2];
1293
1294   // transition cone from ID=196. to ID=216.6
1295   conpar[0] = 32.55/2.;
1296   conpar[1] = 19.6/2.;
1297   conpar[2] = 20.0/2.;
1298   conpar[3] = 21.66/2.;
1299   conpar[4] = 22.06/2.;
1300   TVirtualMC::GetMC()->Gsvolu("QA25", "CONE", idtmed[7], conpar, 5);
1301   TVirtualMC::GetMC()->Gspos("QA25", 1, "ZDCA", 0., 0., conpar[0]+zd2, 0, "ONLY");
1302   // Ch.debug  
1303   //printf("    QA25 CONE from z = %1.2f to z= %1.2f (transition cone)\n",zd2,2*conpar[0]+zd2);
1304
1305   zd2 += 2.*conpar[0]; 
1306   
1307   // tube  
1308   tubpar[0] = 21.66/2.;
1309   tubpar[1] = 22.06/2.;
1310   tubpar[2] = 28.6/2.;
1311   TVirtualMC::GetMC()->Gsvolu("QA26", "TUBE", idtmed[7], tubpar, 3);
1312   TVirtualMC::GetMC()->Gspos("QA26", 1, "ZDCA", 0., 0., tubpar[2]+zd2, 0, "ONLY");
1313   // Ch.debug 
1314   //printf("    QA26  TUBE from z = %1.2f to z= %1.2f\n",zd2,2*tubpar[2]+zd2);
1315
1316   zd2 += 2.*tubpar[2];
1317   // Ch.debug
1318   //printf("    Begin of recombination chamber z = %1.2f\n",zd2);
1319
1320   // --------------------------------------------------------
1321   // RECOMBINATION CHAMBER IMPLEMENTED USING TGeo CLASSES!!!!
1322   // author: Chiara (June 2008)
1323   // --------------------------------------------------------
1324   // TRANSFORMATION MATRICES
1325   // Combi transformation: 
1326   dx = -3.970000;
1327   dy = 0.000000;
1328   dz = 0.0;
1329   // Rotation: 
1330   thx = 84.989100;   phx = 0.000000;
1331   thy = 90.000000;   phy = 90.000000;
1332   thz = 5.010900;    phz = 180.000000;
1333   TGeoRotation *rotMatrix1 = new TGeoRotation("",thx,phx,thy,phy,thz,phz);
1334   // Combi transformation: 
1335   dx = -3.970000;
1336   dy = 0.000000;
1337   dz = 0.0;
1338   TGeoCombiTrans *rotMatrix2 = new TGeoCombiTrans("ZDC_c1", dx,dy,dz,rotMatrix1);
1339   rotMatrix2->RegisterYourself();
1340   // Combi transformation: 
1341   dx = 3.970000;
1342   dy = 0.000000;
1343   dz = 0.0;
1344   // Rotation: 
1345   thx = 95.010900;   phx = 0.000000;
1346   thy = 90.000000;   phy = 90.000000;
1347   thz = 5.010900;    phz = 0.000000;
1348   TGeoRotation *rotMatrix3 = new TGeoRotation("",thx,phx,thy,phy,thz,phz);
1349   TGeoCombiTrans *rotMatrix4 = new TGeoCombiTrans("ZDC_c2", dx,dy,dz,rotMatrix3);
1350   rotMatrix4->RegisterYourself();
1351   
1352   
1353   // VOLUMES DEFINITION
1354   // Volume: ZDCA
1355   TGeoVolume *pZDCA = gGeoManager->GetVolume("ZDCA");
1356   
1357   conpar[0] = (90.1-0.95-0.26)/2.;
1358   conpar[1] = 0.0/2.;
1359   conpar[2] = 21.6/2.;
1360   conpar[3] = 0.0/2.;
1361   conpar[4] = 5.8/2.;
1362   new TGeoCone("QALext", conpar[0],conpar[1],conpar[2],conpar[3],conpar[4]);
1363   
1364   conpar[0] = (90.1-0.95-0.26)/2.;
1365   conpar[1] = 0.0/2.;
1366   conpar[2] = 21.2/2.;
1367   conpar[3] = 0.0/2.;
1368   conpar[4] = 5.4/2.;
1369   new TGeoCone("QALint", conpar[0],conpar[1],conpar[2],conpar[3],conpar[4]);
1370
1371   // Outer trousers
1372   TGeoCompositeShape *pOutTrousers = new TGeoCompositeShape("outTrousers", "QALext:ZDC_c1+QALext:ZDC_c2");
1373   
1374   // Volume: QALext
1375   //TGeoMedium *medZDCFe = gGeoManager->GetMedium("ZDC_ZIRON");
1376   TGeoVolume *pQALext = new TGeoVolume("QALext",pOutTrousers, medZDCFe);
1377   pQALext->SetLineColor(kBlue);
1378   pQALext->SetVisLeaves(kTRUE);
1379   //
1380   TGeoTranslation *tr1 = new TGeoTranslation(0., 0., (Double_t) conpar[0]+0.95+zd2);
1381   pZDCA->AddNode(pQALext, 1, tr1);
1382   // Inner trousers
1383   TGeoCompositeShape *pIntTrousers = new TGeoCompositeShape("intTrousers", "QALint:ZDC_c1+QALint:ZDC_c2");
1384   // Volume: QALint
1385   //TGeoMedium *medZDCvoid = gGeoManager->GetMedium("ZDC_ZVOID");
1386   TGeoVolume *pQALint = new TGeoVolume("QALint",pIntTrousers, medZDCvoid);
1387   pQALint->SetLineColor(kAzure);
1388   pQALint->SetVisLeaves(kTRUE);
1389   pQALext->AddNode(pQALint, 1);
1390     
1391   zd2 += 90.1;
1392   // Ch.debug
1393   //printf("    End of recombination chamber z = %1.2f\n",zd2);
1394   
1395   
1396   //  second section : 2 tubes (ID = 54. OD = 58.)  
1397   tubpar[0] = 5.4/2.;
1398   tubpar[1] = 5.8/2.;
1399   tubpar[2] = 40.0/2.;
1400   TVirtualMC::GetMC()->Gsvolu("QA27", "TUBE", idtmed[7], tubpar, 3);
1401   TVirtualMC::GetMC()->Gspos("QA27", 1, "ZDCA", -15.8/2., 0., tubpar[2]+zd2, 0, "ONLY");
1402   TVirtualMC::GetMC()->Gspos("QA27", 2, "ZDCA",  15.8/2., 0., tubpar[2]+zd2, 0, "ONLY");  
1403   // Ch.debug
1404   //printf("    QA27 TUBE from z = %1.2f to z= %1.2f (separate pipes)\n",zd2,2*tubpar[2]+zd2);
1405   
1406   zd2 += 2.*tubpar[2];
1407  
1408   // transition x2zdc to recombination chamber : skewed cone  
1409   conpar[0] = (10.-1.)/2.;
1410   conpar[1] = 5.4/2.;
1411   conpar[2] = 5.8/2.;
1412   conpar[3] = 6.3/2.;
1413   conpar[4] = 7.0/2.;
1414   TVirtualMC::GetMC()->Gsvolu("QA28", "CONE", idtmed[7], conpar, 5); 
1415   TVirtualMC::GetMC()->Gspos("QA28", 1, "ZDCA", -7.9-0.175, 0., conpar[0]+0.5+zd2, irotpipe1, "ONLY");
1416   TVirtualMC::GetMC()->Gspos("QA28", 2, "ZDCA", 7.9+0.175, 0., conpar[0]+0.5+zd2, irotpipe2, "ONLY");
1417   //printf("    QA28 CONE from z = %1.2f to z= %1.2f (transition X2ZDC)\n",zd2,2*conpar[0]+0.2+zd2);
1418
1419   zd2 += 2.*conpar[0]+1.;
1420   
1421   // 2 tubes (ID = 63 mm OD=70 mm)      
1422   tubpar[0] = 6.3/2.;
1423   tubpar[1] = 7.0/2.;
1424   tubpar[2] = (342.5+498.3)/2.;
1425   TVirtualMC::GetMC()->Gsvolu("QA29", "TUBE", idtmed[7], tubpar, 3);
1426   TVirtualMC::GetMC()->Gspos("QA29", 1, "ZDCA", -16.5/2., 0., tubpar[2]+zd2, 0, "ONLY");
1427   TVirtualMC::GetMC()->Gspos("QA29", 2, "ZDCA",  16.5/2., 0., tubpar[2]+zd2, 0, "ONLY");
1428   //printf("    QA29 TUBE from z = %1.2f to z= %1.2f (separate pipes)\n",zd2,2*tubpar[2]+zd2);  
1429
1430   zd2 += 2.*tubpar[2];
1431            
1432   // -- Luminometer (Cu box) in front of ZN - side A
1433   if(fLumiLength>0.){
1434     boxpar[0] = 8.0/2.;
1435     boxpar[1] = 8.0/2.;
1436     boxpar[2] = fLumiLength/2.;
1437     TVirtualMC::GetMC()->Gsvolu("QLUA", "BOX ", idtmed[9], boxpar, 3);
1438     TVirtualMC::GetMC()->Gspos("QLUA", 1, "ZDCA", 0., 0.,  fPosZNA[2]-66.-boxpar[2], 0, "ONLY");
1439     printf("    A SIDE LUMINOMETER %1.2f < z < %1.2f\n\n",  fPosZNA[2]-66., fPosZNA[2]-66.-2*boxpar[2]);
1440   }
1441   printf("      END OF A SIDE BEAM PIPE VOLUME DEFINITION AT z = %f m from IP2\n",zd2/100.);
1442   
1443
1444   // ----------------------------------------------------------------
1445   // --  MAGNET DEFINITION  -> LHC OPTICS 6.5  
1446   // ----------------------------------------------------------------      
1447   // ***************************************************************  
1448   //            SIDE C - RB26  (dimuon side) 
1449   // ***************************************************************   
1450   // --  COMPENSATOR DIPOLE (MBXW)
1451   zCorrDip = 1972.5;   
1452   
1453   // --  GAP (VACUUM WITH MAGNETIC FIELD)
1454   tubpar[0] = 0.;
1455   tubpar[1] = 3.14;
1456   tubpar[2] = 153./2.;
1457   TVirtualMC::GetMC()->Gsvolu("MBXW", "TUBE", idtmed[11], tubpar, 3);
1458
1459   // --  YOKE 
1460   tubpar[0] = 4.5;
1461   tubpar[1] = 55.;
1462   tubpar[2] = 153./2.;
1463   TVirtualMC::GetMC()->Gsvolu("YMBX", "TUBE", idtmed[7], tubpar, 3);
1464
1465   TVirtualMC::GetMC()->Gspos("MBXW", 1, "ZDCC", 0., 0., -tubpar[2]-zCorrDip, 0, "ONLY");
1466   TVirtualMC::GetMC()->Gspos("YMBX", 1, "ZDCC", 0., 0., -tubpar[2]-zCorrDip, 0, "ONLY");
1467   
1468   
1469   // -- INNER TRIPLET 
1470   zInnTrip = 2296.5; 
1471
1472   // -- DEFINE MQXL AND MQX QUADRUPOLE ELEMENT 
1473   // --  MQXL 
1474   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
1475   tubpar[0] = 0.;
1476   tubpar[1] = 3.14;
1477   tubpar[2] = 637./2.;
1478   TVirtualMC::GetMC()->Gsvolu("MQXL", "TUBE", idtmed[11], tubpar, 3);
1479     
1480   // --  YOKE 
1481   tubpar[0] = 3.5;
1482   tubpar[1] = 22.;
1483   tubpar[2] = 637./2.;
1484   TVirtualMC::GetMC()->Gsvolu("YMQL", "TUBE", idtmed[7], tubpar, 3);
1485   
1486   TVirtualMC::GetMC()->Gspos("MQXL", 1, "ZDCC", 0., 0., -tubpar[2]-zInnTrip, 0, "ONLY");
1487   TVirtualMC::GetMC()->Gspos("YMQL", 1, "ZDCC", 0., 0., -tubpar[2]-zInnTrip, 0, "ONLY");
1488   
1489   TVirtualMC::GetMC()->Gspos("MQXL", 2, "ZDCC", 0., 0., -tubpar[2]-zInnTrip-2400., 0, "ONLY");
1490   TVirtualMC::GetMC()->Gspos("YMQL", 2, "ZDCC", 0., 0., -tubpar[2]-zInnTrip-2400., 0, "ONLY");
1491   
1492   // --  MQX 
1493   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
1494   tubpar[0] = 0.;
1495   tubpar[1] = 3.14;
1496   tubpar[2] = 550./2.;
1497   TVirtualMC::GetMC()->Gsvolu("MQX ", "TUBE", idtmed[11], tubpar, 3);
1498   
1499   // --  YOKE 
1500   tubpar[0] = 3.5;
1501   tubpar[1] = 22.;
1502   tubpar[2] = 550./2.;
1503   TVirtualMC::GetMC()->Gsvolu("YMQ ", "TUBE", idtmed[7], tubpar, 3);
1504   
1505   TVirtualMC::GetMC()->Gspos("MQX ", 1, "ZDCC", 0., 0., -tubpar[2]-zInnTrip-908.5,  0, "ONLY");
1506   TVirtualMC::GetMC()->Gspos("YMQ ", 1, "ZDCC", 0., 0., -tubpar[2]-zInnTrip-908.5,  0, "ONLY");
1507   
1508   TVirtualMC::GetMC()->Gspos("MQX ", 2, "ZDCC", 0., 0., -tubpar[2]-zInnTrip-1558.5, 0, "ONLY");
1509   TVirtualMC::GetMC()->Gspos("YMQ ", 2, "ZDCC", 0., 0., -tubpar[2]-zInnTrip-1558.5, 0, "ONLY");
1510   
1511   // -- SEPARATOR DIPOLE D1 
1512   zD1 = 5838.3001;
1513   
1514   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
1515   tubpar[0] = 0.;
1516   tubpar[1] = 3.46;
1517   tubpar[2] = 945./2.;
1518   TVirtualMC::GetMC()->Gsvolu("MD1 ", "TUBE", idtmed[11], tubpar, 3);
1519   
1520   // --  Insert horizontal Cu plates inside D1 
1521   // --   (to simulate the vacuum chamber)
1522   boxpar[0] = TMath::Sqrt(tubpar[1]*tubpar[1]-(2.98+0.2)*(2.98+0.2)) - 0.05;
1523   boxpar[1] = 0.2/2.;
1524   boxpar[2] = 945./2.;
1525   TVirtualMC::GetMC()->Gsvolu("MD1V", "BOX ", idtmed[6], boxpar, 3);
1526   TVirtualMC::GetMC()->Gspos("MD1V", 1, "MD1 ", 0., 2.98+boxpar[1], 0., 0, "ONLY");
1527   TVirtualMC::GetMC()->Gspos("MD1V", 2, "MD1 ", 0., -2.98-boxpar[1], 0., 0, "ONLY");
1528     
1529   // --  YOKE 
1530   tubpar[0] = 3.68;
1531   tubpar[1] = 110./2.;
1532   tubpar[2] = 945./2.;
1533   TVirtualMC::GetMC()->Gsvolu("YD1 ", "TUBE", idtmed[7], tubpar, 3);
1534   
1535   TVirtualMC::GetMC()->Gspos("YD1 ", 1, "ZDCC", 0., 0., -tubpar[2]-zD1, 0, "ONLY");
1536   TVirtualMC::GetMC()->Gspos("MD1 ", 1, "ZDCC", 0., 0., -tubpar[2]-zD1, 0, "ONLY");
1537   // Ch debug
1538   //printf("    MD1 from z = %1.2f to z= %1.2f cm\n",-zD1, -zD1-2*tubpar[2]); 
1539   
1540   // -- DIPOLE D2 
1541 /*  zD2 = 12167.8;
1542   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
1543   tubpar[0] = 0.;
1544   tubpar[1] = 7.5/2.;
1545   tubpar[2] = 945./2.;
1546   TVirtualMC::GetMC()->Gsvolu("MD2 ", "TUBE", idtmed[11], tubpar, 3);
1547   
1548   // --  YOKE 
1549   tubpar[0] = 0.;
1550   tubpar[1] = 55.;
1551   tubpar[2] = 945./2.;
1552   TVirtualMC::GetMC()->Gsvolu("YD2 ", "TUBE", idtmed[7], tubpar, 3);
1553   
1554   TVirtualMC::GetMC()->Gspos("YD2 ", 1, "ZDCC", 0., 0., -tubpar[2]-zD2, 0, "ONLY");
1555   // Ch debug
1556   //printf("    YD2 from z = %1.2f to z= %1.2f cm\n",-zD2, -zD2-2*tubpar[2]); 
1557   
1558   TVirtualMC::GetMC()->Gspos("MD2 ", 1, "YD2 ", -9.4, 0., 0., 0, "ONLY");
1559   TVirtualMC::GetMC()->Gspos("MD2 ", 2, "YD2 ",  9.4, 0., 0., 0, "ONLY");
1560 */  
1561   // ***************************************************************  
1562   //            SIDE A - RB24 
1563   // ***************************************************************
1564   
1565   // COMPENSATOR DIPOLE (MCBWA) (2nd compensator)
1566   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
1567   tubpar[0] = 0.;
1568   tubpar[1] = 3.;  
1569   tubpar[2] = 153./2.;
1570   TVirtualMC::GetMC()->Gsvolu("MCBW", "TUBE", idtmed[11], tubpar, 3);  
1571   TVirtualMC::GetMC()->Gspos("MCBW", 1, "ZDCA", 0., 0., tubpar[2]+zCorrDip, 0, "ONLY");
1572     
1573    // --  YOKE 
1574   tubpar[0] = 4.5;
1575   tubpar[1] = 55.;
1576   tubpar[2] = 153./2.;
1577   TVirtualMC::GetMC()->Gsvolu("YMCB", "TUBE", idtmed[7], tubpar, 3);
1578   TVirtualMC::GetMC()->Gspos("YMCB", 1, "ZDCA", 0., 0., tubpar[2]+zCorrDip, 0, "ONLY");  
1579   
1580    // -- INNER TRIPLET 
1581   // -- DEFINE MQX1 AND MQX2 QUADRUPOLE ELEMENT 
1582   // --  MQX1 
1583   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
1584   tubpar[0] = 0.;
1585   tubpar[1] = 3.14;
1586   tubpar[2] = 637./2.;
1587   TVirtualMC::GetMC()->Gsvolu("MQX1", "TUBE", idtmed[11], tubpar, 3);
1588   TVirtualMC::GetMC()->Gsvolu("MQX4", "TUBE", idtmed[11], tubpar, 3);
1589     
1590   // --  YOKE 
1591   tubpar[0] = 3.5;
1592   tubpar[1] = 22.;
1593   tubpar[2] = 637./2.;
1594   TVirtualMC::GetMC()->Gsvolu("YMQ1", "TUBE", idtmed[7], tubpar, 3);
1595
1596   // -- Q1
1597   TVirtualMC::GetMC()->Gspos("MQX1", 1, "ZDCA", 0., 0., tubpar[2]+zInnTrip, 0, "ONLY");
1598   TVirtualMC::GetMC()->Gspos("YMQ1", 1, "ZDCA", 0., 0., tubpar[2]+zInnTrip, 0, "ONLY");
1599
1600    // -- BEAM SCREEN FOR Q1
1601    tubpar[0] = 4.78/2.;
1602    tubpar[1] = 5.18/2.;
1603    tubpar[2] = 637./2.;
1604    TVirtualMC::GetMC()->Gsvolu("QBS1", "TUBE", idtmed[6], tubpar, 3);
1605    TVirtualMC::GetMC()->Gspos("QBS1", 1, "MQX1", 0., 0., 0., 0, "ONLY");
1606    // INSERT VERTICAL PLATE INSIDE Q1
1607    boxpar[0] = 0.2/2.0;
1608    boxpar[1] = TMath::Sqrt(tubpar[0]*tubpar[0]-(1.9+0.2)*(1.9+0.2));
1609    boxpar[2] =637./2.;
1610    TVirtualMC::GetMC()->Gsvolu("QBS2", "BOX ", idtmed[6], boxpar, 3);
1611    TVirtualMC::GetMC()->Gspos("QBS2", 1, "MQX1", 1.9+boxpar[0], 0., 0., 0, "ONLY");
1612    TVirtualMC::GetMC()->Gspos("QBS2", 2, "MQX1", -1.9-boxpar[0], 0., 0., 0, "ONLY");
1613
1614    // -- Q3   
1615    TVirtualMC::GetMC()->Gspos("MQX4", 1, "ZDCA", 0., 0., tubpar[2]+zInnTrip+2400., 0, "ONLY");
1616    TVirtualMC::GetMC()->Gspos("YMQ1", 2, "ZDCA", 0., 0., tubpar[2]+zInnTrip+2400., 0, "ONLY");
1617
1618    // -- BEAM SCREEN FOR Q3
1619    tubpar[0] = 5.79/2.;
1620    tubpar[1] = 6.14/2.;
1621    tubpar[2] = 637./2.;
1622    TVirtualMC::GetMC()->Gsvolu("QBS3", "TUBE", idtmed[6], tubpar, 3);
1623    TVirtualMC::GetMC()->Gspos("QBS3", 1, "MQX4", 0., 0., 0., 0, "ONLY");
1624    // INSERT VERTICAL PLATE INSIDE Q3
1625    boxpar[0] = 0.2/2.0;
1626    boxpar[1] = TMath::Sqrt(tubpar[0]*tubpar[0]-(2.405+0.2)*(2.405+0.2));
1627    boxpar[2] =637./2.;
1628    TVirtualMC::GetMC()->Gsvolu("QBS4", "BOX ", idtmed[6], boxpar, 3);
1629    TVirtualMC::GetMC()->Gspos("QBS4", 1, "MQX4", 2.405+boxpar[0], 0., 0., 0, "ONLY");
1630    TVirtualMC::GetMC()->Gspos("QBS4", 2, "MQX4", -2.405-boxpar[0], 0., 0., 0, "ONLY");
1631     
1632   
1633   
1634   // --  MQX2
1635   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
1636   tubpar[0] = 0.;
1637   tubpar[1] = 3.14;
1638   tubpar[2] = 550./2.;
1639   TVirtualMC::GetMC()->Gsvolu("MQX2", "TUBE", idtmed[11], tubpar, 3);
1640   TVirtualMC::GetMC()->Gsvolu("MQX3", "TUBE", idtmed[11], tubpar, 3);
1641   
1642   // --  YOKE 
1643   tubpar[0] = 3.5;
1644   tubpar[1] = 22.;
1645   tubpar[2] = 550./2.;
1646   TVirtualMC::GetMC()->Gsvolu("YMQ2", "TUBE", idtmed[7], tubpar, 3);
1647
1648    // -- BEAM SCREEN FOR Q2
1649    tubpar[0] = 5.79/2.;
1650    tubpar[1] = 6.14/2.;
1651    tubpar[2] = 550./2.;
1652    TVirtualMC::GetMC()->Gsvolu("QBS5", "TUBE", idtmed[6], tubpar, 3);
1653    //    VERTICAL PLATE INSIDE Q2
1654    boxpar[0] = 0.2/2.0;
1655    boxpar[1] = TMath::Sqrt(tubpar[0]*tubpar[0]-(2.405+0.2)*(2.405+0.2));
1656    boxpar[2] =550./2.;
1657    TVirtualMC::GetMC()->Gsvolu("QBS6", "BOX ", idtmed[6], boxpar, 3);
1658
1659   // -- Q2A
1660   TVirtualMC::GetMC()->Gspos("MQX2", 1, "ZDCA", 0., 0., tubpar[2]+zInnTrip+908.5,  0, "ONLY");
1661   TVirtualMC::GetMC()->Gspos("QBS5", 1, "MQX2", 0., 0., 0., 0, "ONLY");  
1662   TVirtualMC::GetMC()->Gspos("QBS6", 1, "MQX2", 2.405+boxpar[0], 0., 0., 0, "ONLY");
1663   TVirtualMC::GetMC()->Gspos("QBS6", 2, "MQX2", -2.405-boxpar[0], 0., 0., 0, "ONLY");  
1664   TVirtualMC::GetMC()->Gspos("YMQ2", 1, "ZDCA", 0., 0., tubpar[2]+zInnTrip+908.5,  0, "ONLY");
1665
1666   
1667   // -- Q2B
1668   TVirtualMC::GetMC()->Gspos("MQX3", 1, "ZDCA", 0., 0., tubpar[2]+zInnTrip+1558.5, 0, "ONLY");
1669   TVirtualMC::GetMC()->Gspos("QBS5", 2, "MQX3", 0., 0., 0., 0, "ONLY");  
1670   TVirtualMC::GetMC()->Gspos("QBS6", 3, "MQX3", 2.405+boxpar[0], 0., 0., 0, "ONLY");
1671   TVirtualMC::GetMC()->Gspos("QBS6", 4, "MQX3", -2.405-boxpar[0], 0., 0., 0, "ONLY");
1672   TVirtualMC::GetMC()->Gspos("YMQ2", 2, "ZDCA", 0., 0., tubpar[2]+zInnTrip+1558.5, 0, "ONLY");
1673
1674   // -- SEPARATOR DIPOLE D1 
1675   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
1676   tubpar[0] = 0.;
1677   tubpar[1] = 6.75/2.;//3.375
1678   tubpar[2] = 945./2.;
1679   TVirtualMC::GetMC()->Gsvolu("MD1L", "TUBE", idtmed[11], tubpar, 3);
1680
1681   // --  The beam screen tube is provided by the beam pipe in D1 (QA03 volume)
1682   // --  Insert the beam screen horizontal Cu plates inside D1  
1683   // --   (to simulate the vacuum chamber)
1684   boxpar[0] = TMath::Sqrt(tubpar[1]*tubpar[1]-(2.885+0.2)*(2.885+0.2));
1685   boxpar[1] = 0.2/2.;
1686   boxpar[2] =945./2.;  
1687   TVirtualMC::GetMC()->Gsvolu("QBS7", "BOX ", idtmed[6], boxpar, 3);
1688   TVirtualMC::GetMC()->Gspos("QBS7", 1, "MD1L", 0., 2.885+boxpar[1],0., 0, "ONLY");
1689   TVirtualMC::GetMC()->Gspos("QBS7", 2, "MD1L", 0., -2.885-boxpar[1],0., 0, "ONLY");  
1690     
1691   // --  YOKE 
1692   tubpar[0] = 3.68;
1693   tubpar[1] = 110./2;
1694   tubpar[2] = 945./2.;
1695   TVirtualMC::GetMC()->Gsvolu("YD1L", "TUBE", idtmed[7], tubpar, 3);
1696   
1697   TVirtualMC::GetMC()->Gspos("YD1L", 1, "ZDCA", 0., 0., tubpar[2]+zD1, 0, "ONLY");  
1698   TVirtualMC::GetMC()->Gspos("MD1L", 1, "ZDCA", 0., 0., tubpar[2]+zD1, 0, "ONLY");  
1699   
1700   // -- DIPOLE D2 
1701   // --  GAP (VACUUM WITH MAGNETIC FIELD) 
1702 /*  tubpar[0] = 0.;
1703   tubpar[1] = 7.5/2.; // this has to be checked
1704   tubpar[2] = 945./2.;
1705   TVirtualMC::GetMC()->Gsvolu("MD2L", "TUBE", idtmed[11], tubpar, 3);
1706   
1707   // --  YOKE 
1708   tubpar[0] = 0.;
1709   tubpar[1] = 55.;
1710   tubpar[2] = 945./2.;
1711   TVirtualMC::GetMC()->Gsvolu("YD2L", "TUBE", idtmed[7], tubpar, 3);
1712   
1713   TVirtualMC::GetMC()->Gspos("YD2L", 1, "ZDCA", 0., 0., tubpar[2]+zD2, 0, "ONLY");
1714   
1715   TVirtualMC::GetMC()->Gspos("MD2L", 1, "YD2L", -9.4, 0., 0., 0, "ONLY");
1716   TVirtualMC::GetMC()->Gspos("MD2L", 2, "YD2L",  9.4, 0., 0., 0, "ONLY");
1717 */  
1718   // -- END OF MAGNET DEFINITION     
1719 }
1720   
1721 //_____________________________________________________________________________
1722 void AliZDCv4::CreateZDC()
1723 {
1724  //
1725  // Create the various ZDCs (ZN + ZP)
1726  //
1727   
1728   Float_t dimPb[6], dimVoid[6];
1729   
1730   Int_t *idtmed = fIdtmed->GetArray();
1731
1732   // Parameters for EM calorimeter geometry
1733   // NB -> parameters used ONLY in CreateZDC()
1734   Float_t kDimZEMPb  = 0.15*(TMath::Sqrt(2.));  // z-dimension of the Pb slice
1735   Float_t kFibRadZEM = 0.0315;                  // External fiber radius (including cladding)
1736   Int_t   fDivZEM[3] = {92, 0, 20};             // Divisions for EM detector
1737   Float_t fDimZEM[6] = {fZEMLength, 3.5, 3.5, 45., 0., 0.}; // Dimensions of EM detector
1738   Float_t fFibZEM2 = fDimZEM[2]/TMath::Sin(fDimZEM[3]*kDegrad)-kFibRadZEM;
1739   Float_t fFibZEM[3] = {0., 0.0275, fFibZEM2};  // Fibers for EM calorimeter
1740
1741 if(!fOnlyZEM){
1742   // Parameters for hadronic calorimeters geometry
1743   // NB -> parameters used ONLY in CreateZDC()
1744   Float_t fGrvZN[3] = {0.03, 0.03, 50.};  // Grooves for neutron detector
1745   Float_t fGrvZP[3] = {0.04, 0.04, 75.};  // Grooves for proton detector
1746   Int_t   fDivZN[3] = {11, 11, 0};        // Division for neutron detector
1747   Int_t   fDivZP[3] = {7, 15, 0};         // Division for proton detector
1748   Int_t   fTowZN[2] = {2, 2};             // Tower for neutron detector
1749   Int_t   fTowZP[2] = {4, 1};             // Tower for proton detector
1750
1751
1752   
1753   //-- Create calorimeters geometry
1754   
1755   // -------------------------------------------------------------------------------
1756   //--> Neutron calorimeter (ZN) 
1757   
1758   TVirtualMC::GetMC()->Gsvolu("ZNEU", "BOX ", idtmed[1], fDimZN, 3); // Passive material  
1759   TVirtualMC::GetMC()->Gsvolu("ZNF1", "TUBE", idtmed[3], fFibZN, 3); // Active material
1760   TVirtualMC::GetMC()->Gsvolu("ZNF2", "TUBE", idtmed[4], fFibZN, 3); 
1761   TVirtualMC::GetMC()->Gsvolu("ZNF3", "TUBE", idtmed[4], fFibZN, 3); 
1762   TVirtualMC::GetMC()->Gsvolu("ZNF4", "TUBE", idtmed[3], fFibZN, 3); 
1763   TVirtualMC::GetMC()->Gsvolu("ZNG1", "BOX ", idtmed[12], fGrvZN, 3); // Empty grooves 
1764   TVirtualMC::GetMC()->Gsvolu("ZNG2", "BOX ", idtmed[12], fGrvZN, 3); 
1765   TVirtualMC::GetMC()->Gsvolu("ZNG3", "BOX ", idtmed[12], fGrvZN, 3); 
1766   TVirtualMC::GetMC()->Gsvolu("ZNG4", "BOX ", idtmed[12], fGrvZN, 3); 
1767   
1768   // Divide ZNEU in towers (for hits purposes) 
1769   
1770   TVirtualMC::GetMC()->Gsdvn("ZNTX", "ZNEU", fTowZN[0], 1); // x-tower 
1771   TVirtualMC::GetMC()->Gsdvn("ZN1 ", "ZNTX", fTowZN[1], 2); // y-tower
1772   
1773   //-- Divide ZN1 in minitowers 
1774   //  fDivZN[0]= NUMBER OF FIBERS PER TOWER ALONG X-AXIS, 
1775   //  fDivZN[1]= NUMBER OF FIBERS PER TOWER ALONG Y-AXIS
1776   //  (4 fibres per minitower) 
1777   
1778   TVirtualMC::GetMC()->Gsdvn("ZNSL", "ZN1 ", fDivZN[1], 2); // Slices 
1779   TVirtualMC::GetMC()->Gsdvn("ZNST", "ZNSL", fDivZN[0], 1); // Sticks
1780   
1781   // --- Position the empty grooves in the sticks (4 grooves per stick)
1782   Float_t dx = fDimZN[0] / fDivZN[0] / 4.;
1783   Float_t dy = fDimZN[1] / fDivZN[1] / 4.;
1784   
1785   TVirtualMC::GetMC()->Gspos("ZNG1", 1, "ZNST", 0.-dx, 0.+dy, 0., 0, "ONLY");
1786   TVirtualMC::GetMC()->Gspos("ZNG2", 1, "ZNST", 0.+dx, 0.+dy, 0., 0, "ONLY");
1787   TVirtualMC::GetMC()->Gspos("ZNG3", 1, "ZNST", 0.-dx, 0.-dy, 0., 0, "ONLY");
1788   TVirtualMC::GetMC()->Gspos("ZNG4", 1, "ZNST", 0.+dx, 0.-dy, 0., 0, "ONLY");
1789   
1790   // --- Position the fibers in the grooves 
1791   TVirtualMC::GetMC()->Gspos("ZNF1", 1, "ZNG1", 0., 0., 0., 0, "ONLY");
1792   TVirtualMC::GetMC()->Gspos("ZNF2", 1, "ZNG2", 0., 0., 0., 0, "ONLY");
1793   TVirtualMC::GetMC()->Gspos("ZNF3", 1, "ZNG3", 0., 0., 0., 0, "ONLY");
1794   TVirtualMC::GetMC()->Gspos("ZNF4", 1, "ZNG4", 0., 0., 0., 0, "ONLY");
1795   
1796   // --- Position the neutron calorimeter in ZDC 
1797   // -- Rotation of ZDCs
1798   Int_t irotzdc;
1799   TVirtualMC::GetMC()->Matrix(irotzdc, 90., 180., 90., 90., 180., 0.);
1800   //
1801   TVirtualMC::GetMC()->Gspos("ZNEU", 1, "ZDCC", fPosZNC[0], fPosZNC[1], fPosZNC[2]-fDimZN[2], irotzdc, "ONLY");
1802   //Ch debug
1803   //printf("\n ZN -> %f < z < %f cm\n",fPosZN[2],fPosZN[2]-2*fDimZN[2]);
1804
1805   // --- Position the neutron calorimeter in ZDC2 (left line) 
1806   // -- No Rotation of ZDCs
1807   TVirtualMC::GetMC()->Gspos("ZNEU", 2, "ZDCA", fPosZNA[0], fPosZNA[1], fPosZNA[2]+fDimZN[2], 0, "ONLY");
1808   //Ch debug
1809   //printf("\n ZN left -> %f < z < %f cm\n",fPosZNl[2],fPosZNl[2]+2*fDimZN[2]);
1810
1811
1812   // -------------------------------------------------------------------------------
1813   //--> Proton calorimeter (ZP)  
1814   
1815   TVirtualMC::GetMC()->Gsvolu("ZPRO", "BOX ", idtmed[2], fDimZP, 3); // Passive material
1816   TVirtualMC::GetMC()->Gsvolu("ZPF1", "TUBE", idtmed[3], fFibZP, 3); // Active material
1817   TVirtualMC::GetMC()->Gsvolu("ZPF2", "TUBE", idtmed[4], fFibZP, 3); 
1818   TVirtualMC::GetMC()->Gsvolu("ZPF3", "TUBE", idtmed[4], fFibZP, 3); 
1819   TVirtualMC::GetMC()->Gsvolu("ZPF4", "TUBE", idtmed[3], fFibZP, 3); 
1820   TVirtualMC::GetMC()->Gsvolu("ZPG1", "BOX ", idtmed[12], fGrvZP, 3); // Empty grooves 
1821   TVirtualMC::GetMC()->Gsvolu("ZPG2", "BOX ", idtmed[12], fGrvZP, 3); 
1822   TVirtualMC::GetMC()->Gsvolu("ZPG3", "BOX ", idtmed[12], fGrvZP, 3); 
1823   TVirtualMC::GetMC()->Gsvolu("ZPG4", "BOX ", idtmed[12], fGrvZP, 3); 
1824     
1825   //-- Divide ZPRO in towers(for hits purposes) 
1826   
1827   TVirtualMC::GetMC()->Gsdvn("ZPTX", "ZPRO", fTowZP[0], 1); // x-tower 
1828   TVirtualMC::GetMC()->Gsdvn("ZP1 ", "ZPTX", fTowZP[1], 2); // y-tower
1829   
1830   
1831   //-- Divide ZP1 in minitowers 
1832   //  fDivZP[0]= NUMBER OF FIBERS ALONG X-AXIS PER MINITOWER, 
1833   //  fDivZP[1]= NUMBER OF FIBERS ALONG Y-AXIS PER MINITOWER
1834   //  (4 fiber per minitower) 
1835   
1836   TVirtualMC::GetMC()->Gsdvn("ZPSL", "ZP1 ", fDivZP[1], 2); // Slices 
1837   TVirtualMC::GetMC()->Gsdvn("ZPST", "ZPSL", fDivZP[0], 1); // Sticks
1838   
1839   // --- Position the empty grooves in the sticks (4 grooves per stick)
1840   dx = fDimZP[0] / fTowZP[0] / fDivZP[0] / 2.;
1841   dy = fDimZP[1] / fTowZP[1] / fDivZP[1] / 2.;
1842   
1843   TVirtualMC::GetMC()->Gspos("ZPG1", 1, "ZPST", 0.-dx, 0.+dy, 0., 0, "ONLY");
1844   TVirtualMC::GetMC()->Gspos("ZPG2", 1, "ZPST", 0.+dx, 0.+dy, 0., 0, "ONLY");
1845   TVirtualMC::GetMC()->Gspos("ZPG3", 1, "ZPST", 0.-dx, 0.-dy, 0., 0, "ONLY");
1846   TVirtualMC::GetMC()->Gspos("ZPG4", 1, "ZPST", 0.+dx, 0.-dy, 0., 0, "ONLY");
1847   
1848   // --- Position the fibers in the grooves 
1849   TVirtualMC::GetMC()->Gspos("ZPF1", 1, "ZPG1", 0., 0., 0., 0, "ONLY");
1850   TVirtualMC::GetMC()->Gspos("ZPF2", 1, "ZPG2", 0., 0., 0., 0, "ONLY");
1851   TVirtualMC::GetMC()->Gspos("ZPF3", 1, "ZPG3", 0., 0., 0., 0, "ONLY");
1852   TVirtualMC::GetMC()->Gspos("ZPF4", 1, "ZPG4", 0., 0., 0., 0, "ONLY");
1853   
1854
1855   // --- Position the proton calorimeter in ZDCC
1856   TVirtualMC::GetMC()->Gspos("ZPRO", 1, "ZDCC", fPosZPC[0], fPosZPC[1], fPosZPC[2]-fDimZP[2], irotzdc, "ONLY");
1857   //Ch debug
1858   //printf("\n ZP -> %f < z < %f cm\n",fPosZP[2],fPosZP[2]-2*fDimZP[2]);
1859   
1860   // --- Position the proton calorimeter in ZDCA
1861   // --- No rotation 
1862   TVirtualMC::GetMC()->Gspos("ZPRO", 2, "ZDCA", fPosZPA[0], fPosZPA[1], fPosZPA[2]+fDimZP[2], 0, "ONLY");
1863   //Ch debug
1864   //printf("\n ZP left -> %f < z < %f cm\n",fPosZPl[2],fPosZPl[2]+2*fDimZP[2]);  
1865 }    
1866   
1867   // -------------------------------------------------------------------------------
1868   // -> EM calorimeter (ZEM)  
1869   
1870   TVirtualMC::GetMC()->Gsvolu("ZEM ", "PARA", idtmed[10], fDimZEM, 6);
1871
1872   Int_t irot1, irot2;
1873   TVirtualMC::GetMC()->Matrix(irot1,0.,0.,90.,90.,-90.,0.);                    // Rotation matrix 1  
1874   TVirtualMC::GetMC()->Matrix(irot2,180.,0.,90.,fDimZEM[3]+90.,90.,fDimZEM[3]);// Rotation matrix 2
1875   //printf("irot1 = %d, irot2 = %d \n", irot1, irot2);
1876   
1877   TVirtualMC::GetMC()->Gsvolu("ZEMF", "TUBE", idtmed[3], fFibZEM, 3);   // Active material
1878
1879   TVirtualMC::GetMC()->Gsdvn("ZETR", "ZEM ", fDivZEM[2], 1);            // Tranches 
1880   
1881   dimPb[0] = kDimZEMPb;                                 // Lead slices 
1882   dimPb[1] = fDimZEM[2];
1883   dimPb[2] = fDimZEM[1];
1884   //dimPb[3] = fDimZEM[3]; //controllare
1885   dimPb[3] = 90.-fDimZEM[3]; //originale
1886   dimPb[4] = 0.;
1887   dimPb[5] = 0.;
1888   TVirtualMC::GetMC()->Gsvolu("ZEL0", "PARA", idtmed[5], dimPb, 6);
1889   TVirtualMC::GetMC()->Gsvolu("ZEL1", "PARA", idtmed[5], dimPb, 6);
1890   TVirtualMC::GetMC()->Gsvolu("ZEL2", "PARA", idtmed[5], dimPb, 6);
1891   
1892   // --- Position the lead slices in the tranche 
1893   Float_t zTran = fDimZEM[0]/fDivZEM[2]; 
1894   Float_t zTrPb = -zTran+kDimZEMPb;
1895   TVirtualMC::GetMC()->Gspos("ZEL0", 1, "ZETR", zTrPb, 0., 0., 0, "ONLY");
1896   TVirtualMC::GetMC()->Gspos("ZEL1", 1, "ZETR", kDimZEMPb, 0., 0., 0, "ONLY");
1897   
1898   // --- Vacuum zone (to be filled with fibres)
1899   dimVoid[0] = (zTran-2*kDimZEMPb)/2.;
1900   dimVoid[1] = fDimZEM[2];
1901   dimVoid[2] = fDimZEM[1];
1902   dimVoid[3] = 90.-fDimZEM[3];
1903   dimVoid[4] = 0.;
1904   dimVoid[5] = 0.;
1905   TVirtualMC::GetMC()->Gsvolu("ZEV0", "PARA", idtmed[10], dimVoid,6);
1906   TVirtualMC::GetMC()->Gsvolu("ZEV1", "PARA", idtmed[10], dimVoid,6);
1907   
1908   // --- Divide the vacuum slice into sticks along x axis
1909   TVirtualMC::GetMC()->Gsdvn("ZES0", "ZEV0", fDivZEM[0], 3); 
1910   TVirtualMC::GetMC()->Gsdvn("ZES1", "ZEV1", fDivZEM[0], 3); 
1911   
1912   // --- Positioning the fibers into the sticks
1913   TVirtualMC::GetMC()->Gspos("ZEMF", 1,"ZES0", 0., 0., 0., irot2, "ONLY");
1914   TVirtualMC::GetMC()->Gspos("ZEMF", 1,"ZES1", 0., 0., 0., irot2, "ONLY");
1915   
1916   // --- Positioning the vacuum slice into the tranche
1917   //Float_t displFib = fDimZEM[1]/fDivZEM[0];
1918   TVirtualMC::GetMC()->Gspos("ZEV0", 1,"ZETR", -dimVoid[0], 0., 0., 0, "ONLY");
1919   TVirtualMC::GetMC()->Gspos("ZEV1", 1,"ZETR", -dimVoid[0]+zTran, 0., 0., 0, "ONLY");
1920
1921   // --- Positioning the ZEM into the ZDC - rotation for 90 degrees  
1922   // NB -> ZEM is positioned in ALIC (instead of in ZDC) volume
1923   TVirtualMC::GetMC()->Gspos("ZEM ", 1,"ALIC", -fPosZEM[0], fPosZEM[1], fPosZEM[2]+fDimZEM[0], irot1, "ONLY");
1924   
1925   // Second EM ZDC (same side w.r.t. IP, just on the other side w.r.t. beam pipe)
1926   TVirtualMC::GetMC()->Gspos("ZEM ", 2,"ALIC", fPosZEM[0], fPosZEM[1], fPosZEM[2]+fDimZEM[0], irot1, "ONLY");
1927   
1928   // --- Adding last slice at the end of the EM calorimeter 
1929   Float_t zLastSlice = fPosZEM[2]+kDimZEMPb+2*fDimZEM[0];
1930   TVirtualMC::GetMC()->Gspos("ZEL2", 1,"ALIC", fPosZEM[0], fPosZEM[1], zLastSlice, irot1, "ONLY");
1931   //Ch debug
1932   //printf("\n ZEM lenght = %f cm\n",2*fZEMLength);
1933   //printf("\n ZEM -> %f < z < %f cm\n",fPosZEM[2],fPosZEM[2]+2*fZEMLength+zLastSlice+kDimZEMPb);
1934   
1935 }
1936  
1937 //_____________________________________________________________________________
1938 void AliZDCv4::CreateMaterials()
1939 {
1940   //
1941   // Create Materials for the Zero Degree Calorimeter
1942   //
1943   Float_t dens, ubuf[1], wmat[3], a[3], z[3];
1944
1945   // --- W alloy -> ZN passive material
1946   dens = 17.6;
1947   a[0] = 183.85;
1948   a[1] = 55.85;
1949   a[2] = 58.71;
1950   z[0] = 74.;
1951   z[1] = 26.;
1952   z[2] = 28.;
1953   wmat[0] = .93;
1954   wmat[1] = .03;
1955   wmat[2] = .04;
1956   AliMixture(1, "WALL", a, z, dens, 3, wmat);
1957
1958   // --- Brass (CuZn)  -> ZP passive material
1959   dens = 8.48;
1960   a[0] = 63.546;
1961   a[1] = 65.39;
1962   z[0] = 29.;
1963   z[1] = 30.;
1964   wmat[0] = .63;
1965   wmat[1] = .37;
1966   AliMixture(2, "BRASS", a, z, dens, 2, wmat);
1967   
1968   // --- SiO2 
1969   dens = 2.64;
1970   a[0] = 28.086;
1971   a[1] = 15.9994;
1972   z[0] = 14.;
1973   z[1] = 8.;
1974   wmat[0] = 1.;
1975   wmat[1] = 2.;
1976   AliMixture(3, "SIO2", a, z, dens, -2, wmat);  
1977   
1978   // --- Lead 
1979   ubuf[0] = 1.12;
1980   AliMaterial(5, "LEAD", 207.19, 82., 11.35, .56, 0., ubuf, 1);
1981
1982   // --- Copper (energy loss taken into account)
1983   ubuf[0] = 1.10;
1984   AliMaterial(6, "COPP0", 63.54, 29., 8.96, 1.4, 0., ubuf, 1);
1985
1986   // --- Copper 
1987   ubuf[0] = 1.10;
1988   AliMaterial(9, "COPP1", 63.54, 29., 8.96, 1.4, 0., ubuf, 1);
1989   
1990   // --- Iron (energy loss taken into account)
1991   ubuf[0] = 1.1;
1992   AliMaterial(7, "IRON0", 55.85, 26., 7.87, 1.76, 0., ubuf, 1);
1993   
1994   // --- Iron (no energy loss)
1995   ubuf[0] = 1.1;
1996   AliMaterial(8, "IRON1", 55.85, 26., 7.87, 1.76, 0., ubuf, 1);
1997   
1998   // --- Tatalum 
1999   ubuf[0] = 1.1;
2000   AliMaterial(13, "TANT", 183.84, 74., 19.3, 0.35, 0., ubuf, 1);
2001     
2002   // ---------------------------------------------------------  
2003   Float_t aResGas[3]={1.008,12.0107,15.9994};
2004   Float_t zResGas[3]={1.,6.,8.};
2005   Float_t wResGas[3]={0.28,0.28,0.44};
2006   Float_t dResGas = 3.2E-14;
2007
2008   // --- Vacuum (no magnetic field) 
2009   AliMixture(10, "VOID", aResGas, zResGas, dResGas, 3, wResGas);
2010   
2011   // --- Vacuum (with magnetic field) 
2012   AliMixture(11, "VOIM", aResGas, zResGas, dResGas, 3, wResGas);
2013   
2014   // --- Air (no magnetic field)
2015   Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
2016   Float_t zAir[4]={6.,7.,8.,18.};
2017   Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
2018   Float_t dAir = 1.20479E-3;
2019   //
2020   AliMixture(12, "Air    $", aAir, zAir, dAir, 4, wAir);
2021   
2022   // ---  Definition of tracking media: 
2023   
2024   // --- Tantalum = 1 ; 
2025   // --- Brass = 2 ; 
2026   // --- Fibers (SiO2) = 3 ; 
2027   // --- Fibers (SiO2) = 4 ; 
2028   // --- Lead = 5 ; 
2029   // --- Copper (with high thr.)= 6 ;
2030   // --- Copper (with low thr.)=  9;
2031   // --- Iron (with energy loss) = 7 ; 
2032   // --- Iron (without energy loss) = 8 ; 
2033   // --- Vacuum (no field) = 10 
2034   // --- Vacuum (with field) = 11 
2035   // --- Air (no field) = 12 
2036   
2037   // **************************************************** 
2038   //     Tracking media parameters
2039   //
2040   Float_t epsil  = 0.01;   // Tracking precision, 
2041   Float_t stmin  = 0.01;   // Min. value 4 max. step (cm)
2042   Float_t stemax = 1.;     // Max. step permitted (cm) 
2043   Float_t tmaxfd = 0.;     // Maximum angle due to field (degrees) 
2044   Float_t tmaxfdv = 0.1;   // Maximum angle due to field (degrees) 
2045   Float_t deemax = -1.;    // Maximum fractional energy loss
2046   Float_t nofieldm = 0.;   // Max. field value (no field)
2047   Float_t fieldm = 45.;    // Max. field value (with field)
2048   Int_t isvol = 0;         // ISVOL =0 -> not sensitive volume
2049   Int_t isvolActive = 1;   // ISVOL =1 -> sensitive volume
2050   Int_t inofld = 0;        // IFIELD=0 -> no magnetic field
2051   Int_t ifield =2;         // IFIELD=2 -> magnetic field defined in AliMagFC.h
2052   // *****************************************************
2053   
2054   AliMedium(1, "ZWALL", 1, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2055   AliMedium(2, "ZBRASS",2, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2056   AliMedium(3, "ZSIO2", 3, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2057   AliMedium(4, "ZQUAR", 3, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2058   AliMedium(5, "ZLEAD", 5, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2059   AliMedium(6, "ZCOPP", 6, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2060   AliMedium(7, "ZIRON", 7, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2061   AliMedium(8, "ZIRONN",8, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2062   AliMedium(9, "ZCOPL", 6, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2063   AliMedium(10,"ZVOID",10, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2064   AliMedium(11,"ZVOIM",11, isvol, ifield, fieldm, tmaxfdv, stemax, deemax, epsil, stmin);
2065   AliMedium(12,"ZAIR", 12, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2066   AliMedium(13,"ZTANT",13, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2067   AliMedium(14, "ZIRONT", 7, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2068
2069
2070
2071 //_____________________________________________________________________________
2072 void AliZDCv4::AddAlignableVolumes() const
2073 {
2074  //
2075  // Create entries for alignable volumes associating the symbolic volume
2076  // name with the corresponding volume path. Needs to be syncronized with
2077  // eventual changes in the geometry.
2078  //
2079  if(fOnlyZEM) return;
2080  
2081  TString volpath1 = "ALIC_1/ZDCC_1/ZNEU_1";
2082  TString volpath2 = "ALIC_1/ZDCC_1/ZPRO_1";
2083  TString volpath3 = "ALIC_1/ZDCA_1/ZNEU_2";
2084  TString volpath4 = "ALIC_1/ZDCA_1/ZPRO_2";
2085
2086  TString symname1="ZDC/NeutronZDC_C";
2087  TString symname2="ZDC/ProtonZDC_C";
2088  TString symname3="ZDC/NeutronZDC_A";
2089  TString symname4="ZDC/ProtonZDC_A";
2090
2091  if(!gGeoManager->SetAlignableEntry(symname1.Data(),volpath1.Data()))
2092      AliFatal(Form("Alignable entry %s not created. Volume path %s not valid",   symname1.Data(),volpath1.Data()));
2093
2094  if(!gGeoManager->SetAlignableEntry(symname2.Data(),volpath2.Data()))
2095      AliFatal(Form("Alignable entry %s not created. Volume path %s not valid",   symname2.Data(),volpath2.Data()));
2096
2097  if(!gGeoManager->SetAlignableEntry(symname3.Data(),volpath3.Data()))
2098      AliFatal(Form("Alignable entry %s not created. Volume path %s not valid",   symname1.Data(),volpath1.Data()));
2099
2100  if(!gGeoManager->SetAlignableEntry(symname4.Data(),volpath4.Data()))
2101      AliFatal(Form("Alignable entry %s not created. Volume path %s not valid",   symname2.Data(),volpath2.Data()));
2102
2103 }
2104
2105
2106 //_____________________________________________________________________________
2107 void AliZDCv4::Init()
2108 {
2109  InitTables();
2110   Int_t *idtmed = fIdtmed->GetArray();  
2111   //
2112   fMedSensZN     = idtmed[1];  // Sensitive volume: ZN passive material
2113   fMedSensZP     = idtmed[2];  // Sensitive volume: ZP passive material
2114   fMedSensF1     = idtmed[3];  // Sensitive volume: fibres type 1
2115   fMedSensF2     = idtmed[4];  // Sensitive volume: fibres type 2
2116   fMedSensZEM    = idtmed[5];  // Sensitive volume: ZEM passive material
2117   fMedSensTDI    = idtmed[6];  // Sensitive volume: TDI Cu shield
2118   fMedSensPI     = idtmed[7];  // Sensitive volume: beam pipes
2119   fMedSensLumi   = idtmed[9];  // Sensitive volume: luminometer
2120   fMedSensGR     = idtmed[12]; // Sensitive volume: air into the grooves
2121   fMedSensVColl  = idtmed[13]; // Sensitive volume: collimator jaws
2122 }
2123
2124 //_____________________________________________________________________________
2125 void AliZDCv4::InitTables()
2126 {
2127  //
2128  // Read light tables for Cerenkov light production parameterization 
2129  //
2130
2131   Int_t k, j;
2132   int read=1;
2133
2134   //  --- Reading light tables for ZN 
2135   char *lightfName1 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620362207s");
2136   FILE *fp1 = fopen(lightfName1,"r");
2137   if(fp1 == NULL){
2138      printf("Cannot open file fp1 \n");
2139      return;
2140   }
2141   else{
2142     for(k=0; k<fNalfan; k++){
2143       for(j=0; j<fNben; j++){
2144        read = fscanf(fp1,"%f",&fTablen[0][k][j]);
2145        if(read==0) AliDebug(3, " Error in reading light table 1");
2146       }
2147     }
2148     fclose(fp1);
2149   }
2150   char *lightfName2 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620362208s");
2151   FILE *fp2 = fopen(lightfName2,"r");
2152   if(fp2 == NULL){
2153      printf("Cannot open file fp2 \n");
2154      return;
2155   }  
2156   else{
2157     for(k=0; k<fNalfan; k++){
2158       for(j=0; j<fNben; j++){
2159        read = fscanf(fp2,"%f",&fTablen[1][k][j]);
2160        if(read==0) AliDebug(3, " Error in reading light table 2");
2161       }
2162     }
2163     fclose(fp2);
2164   }
2165   char *lightfName3 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620362209s");
2166   FILE *fp3 = fopen(lightfName3,"r");
2167   if(fp3 == NULL){
2168      printf("Cannot open file fp3 \n");
2169      return;
2170   }
2171   else{
2172     for(k=0; k<fNalfan; k++){
2173       for(j=0; j<fNben; j++){
2174        read = fscanf(fp3,"%f",&fTablen[2][k][j]);
2175        if(read==0) AliDebug(3, " Error in reading light table 3");
2176       }
2177     }
2178     fclose(fp3);
2179   }
2180   char *lightfName4 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620362210s");
2181   FILE *fp4 = fopen(lightfName4,"r");
2182   if(fp4 == NULL){
2183      printf("Cannot open file fp4 \n");
2184      return;
2185   }
2186   else{
2187     for(k=0; k<fNalfan; k++){
2188       for(j=0; j<fNben; j++){
2189        read = fscanf(fp4,"%f",&fTablen[3][k][j]);
2190        if(read==0) AliDebug(3, " Error in reading light table 4");
2191       }
2192     }
2193     fclose(fp4);
2194   }
2195     
2196   //  --- Reading light tables for ZP and ZEM
2197   char *lightfName5 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620552207s");
2198   FILE *fp5 = fopen(lightfName5,"r");
2199   if(fp5 == NULL){
2200      printf("Cannot open file fp5 \n");
2201      return;
2202   }
2203   else{
2204     for(k=0; k<fNalfap; k++){
2205       for(j=0; j<fNbep; j++){
2206        read = fscanf(fp5,"%f",&fTablep[0][k][j]);
2207        if(read==0) AliDebug(3, " Error in reading light table 5");
2208       }
2209     }
2210     fclose(fp5);
2211   }
2212   char *lightfName6 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620552208s");
2213   FILE *fp6 = fopen(lightfName6,"r");
2214   if(fp6 == NULL){
2215      printf("Cannot open file fp6 \n");
2216      return;
2217   }
2218   else{
2219     for(k=0; k<fNalfap; k++){
2220       for(j=0; j<fNbep; j++){
2221        read = fscanf(fp6,"%f",&fTablep[1][k][j]);
2222        if(read==0) AliDebug(3, " Error in reading light table 6");
2223       }
2224     }
2225     fclose(fp6);
2226   }
2227   char *lightfName7 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620552209s");
2228   FILE *fp7 = fopen(lightfName7,"r");
2229   if(fp7 == NULL){
2230      printf("Cannot open file fp7 \n");
2231      return;
2232   }
2233   else{
2234     for(k=0; k<fNalfap; k++){
2235       for(j=0; j<fNbep; j++){
2236        read = fscanf(fp7,"%f",&fTablep[2][k][j]);
2237        if(read==0) AliDebug(3, " Error in reading light table 7");
2238       }
2239     }
2240    fclose(fp7);
2241   }
2242   char *lightfName8 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620552210s");
2243   FILE *fp8 = fopen(lightfName8,"r");
2244   if(fp8 == NULL){
2245      printf("Cannot open file fp8 \n");
2246      return;
2247   }
2248   else{
2249     for(k=0; k<fNalfap; k++){
2250       for(j=0; j<fNbep; j++){
2251        read = fscanf(fp8,"%f",&fTablep[3][k][j]);
2252        if(read==0) AliDebug(3, " Error in reading light table 8");
2253       }
2254     }
2255    fclose(fp8);
2256   }
2257
2258 }
2259 //_____________________________________________________________________________
2260 void AliZDCv4::StepManager()
2261 {
2262   //
2263   // Routine called at every step in the Zero Degree Calorimeters
2264   //
2265   Int_t   j, vol[2]={0,0}, ibeta=0, ialfa=0, ibe=0, nphe=0;
2266   Float_t hits[14], x[3], xdet[3]={999.,999.,999.}, um[3], ud[3];
2267   Float_t destep=0., be=0., out=0.;
2268   Double_t s[3], p[4];
2269   const char *knamed;
2270   //
2271   for(j=0;j<14;j++) hits[j]=-999.;
2272   //
2273   // --- This part is for no shower developement in beam pipe, TDI, VColl
2274   // If particle interacts with beam pipe, TDI, VColl -> return
2275   if(fNoShower==1 && ((TVirtualMC::GetMC()->CurrentMedium() == fMedSensPI) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensTDI) ||  
2276      (TVirtualMC::GetMC()->CurrentMedium() == fMedSensVColl || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensLumi)))){ 
2277     
2278     // If option NoShower is set -> StopTrack
2279
2280     Int_t ipr = 0; 
2281       TVirtualMC::GetMC()->TrackPosition(s[0],s[1],s[2]);
2282       if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensPI){
2283         knamed = TVirtualMC::GetMC()->CurrentVolName();
2284         if(!strncmp(knamed,"YMQ",3)){
2285           if(s[2]<0) fpLostITC += 1;
2286           else fpLostITA += 1;
2287           ipr=1;
2288         }
2289         else if(!strncmp(knamed,"YD1",3)){
2290           if(s[2]<0) fpLostD1C += 1;
2291           else fpLostD1A += 1;
2292           ipr=1;
2293         }
2294       }
2295       else if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensTDI){ 
2296         knamed = TVirtualMC::GetMC()->CurrentVolName();
2297         if(!strncmp(knamed,"MD1",3)){
2298           if(s[2]<0) fpLostD1C += 1;
2299           else  fpLostD1A += 1;
2300           ipr=1;
2301         }
2302         else if(!strncmp(knamed,"QTD",3)) fpLostTDI += 1;
2303       }
2304       else if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensVColl){ 
2305         knamed = TVirtualMC::GetMC()->CurrentVolName();
2306         if(!strncmp(knamed,"QCVC",4)) fpcVCollC++;
2307         else if(!strncmp(knamed,"QCVA",4))  fpcVCollA++;
2308         ipr=1;
2309       }
2310       //
2311       //TVirtualMC::GetMC()->TrackMomentum(p[0], p[1], p[2], p[3]);
2312       //printf("\t Particle: mass = %1.3f, E = %1.3f GeV, pz = %1.2f GeV -> stopped in volume %s\n", 
2313       //     TVirtualMC::GetMC()->TrackMass(), p[3], p[2], TVirtualMC::GetMC()->CurrentVolName());
2314       //
2315       if(ipr<0){
2316         printf("\n\t **********************************\n");
2317         printf("\t ********** Side C **********\n");
2318         printf("\t # of particles in IT = %d\n",fpLostITC);
2319         printf("\t # of particles in D1 = %d\n",fpLostD1C);
2320         printf("\t # of particles in VColl = %d\n",fpcVCollC);
2321         printf("\t ********** Side A **********\n");
2322         printf("\t # of particles in IT = %d\n",fpLostITA);
2323         printf("\t # of particles in D1 = %d\n",fpLostD1A);
2324         printf("\t # of particles in TDI = %d\n",fpLostTDI);
2325         printf("\t # of particles in VColl = %d\n",fpcVCollA);
2326         printf("\t **********************************\n");
2327       }
2328       TVirtualMC::GetMC()->StopTrack();
2329       return;
2330   }
2331   
2332   if((TVirtualMC::GetMC()->CurrentMedium() == fMedSensZN) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensZP) ||
2333      (TVirtualMC::GetMC()->CurrentMedium() == fMedSensGR) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensF1) ||
2334      (TVirtualMC::GetMC()->CurrentMedium() == fMedSensF2) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensZEM)){
2335
2336     
2337   //Particle coordinates 
2338     TVirtualMC::GetMC()->TrackPosition(s[0],s[1],s[2]);
2339     for(j=0; j<=2; j++) x[j] = s[j];
2340     hits[0] = x[0];
2341     hits[1] = x[1];
2342     hits[2] = x[2];
2343
2344   // Determine in which ZDC the particle is
2345     knamed = TVirtualMC::GetMC()->CurrentVolName();
2346     if(!strncmp(knamed,"ZN",2)){
2347           if(x[2]<0.) vol[0]=1; // ZNC (dimuon side)
2348           else if(x[2]>0.) vol[0]=4; //ZNA
2349     }
2350     else if(!strncmp(knamed,"ZP",2)){ 
2351           if(x[2]<0.) vol[0]=2; //ZPC (dimuon side)
2352           else if(x[2]>0.) vol[0]=5; //ZPA  
2353     }
2354     else if(!strncmp(knamed,"ZE",2)) vol[0]=3; //ZEM
2355   
2356   // Determine in which quadrant the particle is
2357     if(vol[0]==1){      //Quadrant in ZNC
2358       // Calculating particle coordinates inside ZNC
2359       xdet[0] = x[0]-fPosZNC[0];
2360       xdet[1] = x[1]-fPosZNC[1];
2361       // Calculating quadrant in ZN
2362       if(xdet[0]<=0.){
2363         if(xdet[1]<=0.) vol[1]=1;
2364         else vol[1]=3;
2365       }
2366       else if(xdet[0]>0.){
2367         if(xdet[1]<=0.) vol[1]=2;
2368         else vol[1]=4;
2369       }
2370     }
2371     
2372     else if(vol[0]==2){ //Quadrant in ZPC
2373       // Calculating particle coordinates inside ZPC
2374       xdet[0] = x[0]-fPosZPC[0];
2375       xdet[1] = x[1]-fPosZPC[1];
2376       if(xdet[0]>=fDimZP[0])  xdet[0]=fDimZP[0]-0.01;
2377       if(xdet[0]<=-fDimZP[0]) xdet[0]=-fDimZP[0]+0.01;
2378       // Calculating tower in ZP
2379       Float_t xqZP = xdet[0]/(fDimZP[0]/2.);
2380       for(int i=1; i<=4; i++){
2381          if(xqZP>=(i-3) && xqZP<(i-2)){
2382            vol[1] = i;
2383            break;
2384          }
2385       }
2386     }
2387     //
2388     // Quadrant in ZEM: vol[1] = 1 -> particle in 1st ZEM (placed at x = 8.5 cm)
2389     //                  vol[1] = 2 -> particle in 2nd ZEM (placed at x = -8.5 cm)
2390     else if(vol[0] == 3){       
2391       if(x[0]>0.){
2392         vol[1] = 1;
2393         // Particle x-coordinate inside ZEM1
2394         xdet[0] = x[0]-fPosZEM[0];
2395       }
2396       else{
2397         vol[1] = 2;
2398         // Particle x-coordinate inside ZEM2
2399         xdet[0] = x[0]+fPosZEM[0];
2400       }
2401       xdet[1] = x[1]-fPosZEM[1];
2402     }
2403     //
2404     else if(vol[0]==4){ //Quadrant in ZNA
2405       // Calculating particle coordinates inside ZNA
2406       xdet[0] = x[0]-fPosZNA[0];
2407       xdet[1] = x[1]-fPosZNA[1];
2408       // Calculating quadrant in ZNA
2409       if(xdet[0]>=0.){
2410         if(xdet[1]<=0.) vol[1]=1;
2411         else vol[1]=3;
2412       }
2413       else if(xdet[0]<0.){
2414         if(xdet[1]<=0.) vol[1]=2;
2415         else vol[1]=4;
2416       }
2417     }    
2418     //
2419     else if(vol[0]==5){ //Quadrant in ZPA
2420       // Calculating particle coordinates inside ZPA
2421       xdet[0] = x[0]-fPosZPA[0];
2422       xdet[1] = x[1]-fPosZPA[1];
2423       if(xdet[0]>=fDimZP[0])  xdet[0]=fDimZP[0]-0.01;
2424       if(xdet[0]<=-fDimZP[0]) xdet[0]=-fDimZP[0]+0.01;
2425       // Calculating tower in ZP
2426       Float_t xqZP = -xdet[0]/(fDimZP[0]/2.);
2427       for(int i=1; i<=4; i++){
2428          if(xqZP>=(i-3) && xqZP<(i-2)){
2429            vol[1] = i;
2430            break;
2431          }
2432       }
2433     }    
2434     if((vol[1]!=1) && (vol[1]!=2) && (vol[1]!=3) && (vol[1]!=4))
2435       AliError(Form(" WRONG tower for det %d: tow %d with xdet=(%f, %f)\n",
2436                 vol[0], vol[1], xdet[0], xdet[1]));
2437     // Ch. debug
2438     //printf("\t *** det %d vol %d xdet(%f, %f)\n",vol[0], vol[1], xdet[0], xdet[1]);
2439     
2440     
2441     // Store impact point and kinetic energy of the ENTERING particle
2442     
2443     if(TVirtualMC::GetMC()->IsTrackEntering()){
2444       //Particle energy
2445       TVirtualMC::GetMC()->TrackMomentum(p[0],p[1],p[2],p[3]);
2446       hits[3] = p[3];
2447       
2448       // Impact point on ZDC
2449       // X takes into account the LHC x-axis sign
2450       // which is opposite to positive x on detector front face
2451       // for side A detectors (ZNA and ZPA)  
2452       if(vol[0]==4 || vol[0]==5){
2453         hits[4] = -xdet[0];
2454       }
2455       else{
2456         hits[4] = xdet[0];
2457       }
2458       hits[5] = xdet[1];
2459       hits[6] = 0;
2460       hits[7] = 0;
2461       hits[8] = 0;
2462       hits[9] = 0;
2463       //
2464       Int_t curTrackN = gAlice->GetMCApp()->GetCurrentTrackNumber();
2465       TParticle *part = gAlice->GetMCApp()->Particle(curTrackN);
2466       hits[10] = part->GetPdgCode();
2467       hits[11] = 0;
2468       hits[12] = 1.0e09*TVirtualMC::GetMC()->TrackTime(); // in ns!
2469       hits[13] = part->Eta();
2470       //
2471       if(fFindMother){
2472          Int_t imo = part->GetFirstMother();
2473          //printf(" tracks: pc %d -> mother %d \n", curTrackN,imo); 
2474          
2475          int trmo = imo;
2476          TParticle *pmot = 0x0;
2477          Bool_t isChild = kFALSE;
2478          if(imo>-1){
2479            pmot = gAlice->GetMCApp()->Particle(imo);
2480            trmo = pmot->GetFirstMother();
2481            isChild = kTRUE;
2482            while(trmo!=-1){
2483               pmot = gAlice->GetMCApp()->Particle(trmo);
2484               //printf("  **** pc %d -> mother %d \n", trch,trmo); 
2485               trmo = pmot->GetFirstMother();
2486            }
2487          }
2488       
2489          if(isChild && pmot){
2490              hits[6]  = 1;
2491              hits[11] = pmot->GetPdgCode();
2492              hits[13] = pmot->Eta();
2493          }
2494       }
2495       
2496
2497       AddHit(curTrackN, vol, hits);
2498
2499       if(fNoShower==1){
2500         if(vol[0]==1){
2501           fnDetectedC += 1;
2502           //if(fnDetectedC==1) printf(" ### Particle in ZNC\n\n");
2503         }
2504         else if(vol[0]==2){
2505           fpDetectedC += 1;
2506           //if(fpDetectedC==1) printf(" ### Particle in ZPC\n\n");
2507         }
2508         //else if(vol[0]==3) printf("   ### Particle in ZEM\n\n");        
2509         else if(vol[0]==4){
2510           fnDetectedA += 1;
2511           //if(fnDetectedA==1) printf(" ### Particle in ZNA\n\n");        
2512         }
2513         else if(vol[0]==5){
2514           fpDetectedA += 1;
2515           //if(fpDetectedA==1) printf(" ### Particle in ZPA\n\n");       
2516         }
2517         //
2518         //printf("\t Pc: x %1.2f y %1.2f z %1.2f  E %1.2f GeV pz = %1.2f GeV in volume %s\n", 
2519         //   x[0],x[1],x[3],p[3],p[2],TVirtualMC::GetMC()->CurrentVolName());
2520         //
2521         TVirtualMC::GetMC()->StopTrack();
2522         return;
2523       }
2524     }
2525            
2526     // Particle energy loss
2527     if(TVirtualMC::GetMC()->Edep() != 0){
2528       hits[9] = TVirtualMC::GetMC()->Edep();
2529       hits[7] = 0.;
2530       hits[8] = 0.;
2531       AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2532     }
2533   }
2534  
2535
2536   // *** Light production in fibres 
2537   if((TVirtualMC::GetMC()->CurrentMedium() == fMedSensF1) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensF2)){
2538
2539      //Select charged particles
2540      if((destep=TVirtualMC::GetMC()->Edep())){
2541
2542        // Particle velocity
2543        Float_t beta = 0.;
2544        TVirtualMC::GetMC()->TrackMomentum(p[0],p[1],p[2],p[3]);
2545        Float_t ptot=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
2546        if(p[3] > 0.00001) beta =  ptot/p[3];
2547        else return;
2548        if(beta<0.67)return;
2549        else if((beta>=0.67) && (beta<=0.75)) ibeta = 0;
2550        else if((beta>0.75)  && (beta<=0.85)) ibeta = 1;
2551        else if((beta>0.85)  && (beta<=0.95)) ibeta = 2;
2552        else if(beta>0.95) ibeta = 3;
2553  
2554        // Angle between particle trajectory and fibre axis
2555        // 1 -> Momentum directions
2556        um[0] = p[0]/ptot;
2557        um[1] = p[1]/ptot;
2558        um[2] = p[2]/ptot;
2559        TVirtualMC::GetMC()->Gmtod(um,ud,2);
2560        // 2 -> Angle < limit angle
2561        Double_t alfar = TMath::ACos(ud[2]);
2562        Double_t alfa = alfar*kRaddeg;
2563        if(alfa>=110.) return;
2564        //
2565        ialfa = Int_t(1.+alfa/2.);
2566  
2567        // Distance between particle trajectory and fibre axis
2568        TVirtualMC::GetMC()->TrackPosition(s[0],s[1],s[2]);
2569        for(j=0; j<=2; j++){
2570           x[j] = s[j];
2571        }
2572        TVirtualMC::GetMC()->Gmtod(x,xdet,1);
2573        if(TMath::Abs(ud[0])>0.00001){
2574          Float_t dcoeff = ud[1]/ud[0];
2575          be = TMath::Abs((xdet[1]-dcoeff*xdet[0])/TMath::Sqrt(dcoeff*dcoeff+1.));
2576        }
2577        else{
2578          be = TMath::Abs(ud[0]);
2579        }
2580  
2581        ibe = Int_t(be*1000.+1);
2582   
2583        //Looking into the light tables 
2584        Float_t charge = 0.;
2585        Int_t curTrackN = gAlice->GetMCApp()->GetCurrentTrackNumber();
2586        TParticle *part = gAlice->GetMCApp()->Particle(curTrackN);
2587        Int_t pdgCode = part->GetPdgCode();
2588        if(pdgCode<10000) charge = TVirtualMC::GetMC()->TrackCharge();
2589        else{
2590           float z = (pdgCode/10000-100000);
2591           charge = TMath::Abs(z);
2592           //printf(" PDG %d   charge %f\n",pdgCode,charge);
2593        } 
2594        
2595        if(vol[0]==1 || vol[0]==4) {     // (1)  ZN fibres
2596          if(ibe>fNben) ibe=fNben;
2597          out =  charge*charge*fTablen[ibeta][ialfa][ibe];
2598          nphe = gRandom->Poisson(out);
2599          // Ch. debug
2600          //if(ibeta==3) printf("\t %f \t %f \t %f\n",alfa, be, out);
2601          //printf("\t ibeta = %d, ialfa = %d, ibe = %d -> nphe = %d\n\n",ibeta,ialfa,ibe,nphe);
2602          if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensF1){
2603            hits[7] = nphe;      //fLightPMQ
2604            hits[8] = 0;
2605            hits[9] = 0;
2606            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2607          }
2608          else{
2609            hits[7] = 0;
2610            hits[8] = nphe;      //fLightPMC
2611            hits[9] = 0;
2612            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2613          }
2614        } 
2615        else if(vol[0]==2 || vol[0]==5) {// (2) ZP fibres
2616          if(ibe>fNbep) ibe=fNbep;
2617          out =  charge*charge*fTablep[ibeta][ialfa][ibe];
2618          nphe = gRandom->Poisson(out);
2619          if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensF1){
2620            hits[7] = nphe;      //fLightPMQ
2621            hits[8] = 0;
2622            hits[9] = 0;
2623            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2624          }
2625          else{
2626            hits[7] = 0;
2627            hits[8] = nphe;      //fLightPMC
2628            hits[9] = 0;
2629            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2630          }
2631        } 
2632        else if(vol[0]==3) {     // (3) ZEM fibres
2633          if(ibe>fNbep) ibe=fNbep;
2634          out =  charge*charge*fTablep[ibeta][ialfa][ibe];
2635          TVirtualMC::GetMC()->TrackPosition(s[0],s[1],s[2]);
2636          Float_t xalic[3];
2637          for(j=0; j<3; j++){
2638             xalic[j] = s[j];
2639          }
2640          // z-coordinate from ZEM front face 
2641          // NB-> fPosZEM[2]+fZEMLength = -1000.+2*10.3 = 979.69 cm
2642          Float_t z = -xalic[2]+fPosZEM[2]+2*fZEMLength-xalic[1];
2643          //z = xalic[2]-fPosZEM[2]-fZEMLength-xalic[1]*(TMath::Tan(45.*kDegrad));
2644          //printf("     fPosZEM[2]+2*fZEMLength = %f", fPosZEM[2]+2*fZEMLength);
2645          //
2646          // Parametrization for light guide uniformity
2647          // NEW!!! Light guide tilted @ 51 degrees
2648          Float_t guiPar[4]={0.31,-0.0006305,0.01337,0.8895};
2649          Float_t guiEff = guiPar[0]*(guiPar[1]*z*z+guiPar[2]*z+guiPar[3]);
2650          out = out*guiEff;
2651          nphe = gRandom->Poisson(out);
2652          //printf("     out*guiEff = %f nphe = %d", out, nphe);
2653          if(vol[1] == 1){
2654            hits[7] = 0;         
2655            hits[8] = nphe;      //fLightPMC (ZEM1)
2656            hits[9] = 0;
2657            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2658          }
2659          else{
2660            hits[7] = nphe;      //fLightPMQ (ZEM2)
2661            hits[8] = 0;         
2662            hits[9] = 0;
2663            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2664          }
2665        }
2666      }
2667    }
2668 }