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