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