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