ZDC/ZDCsim/AliZDCv4.cxx

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