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