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