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