ZDC/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.2),
  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.2),
 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] = 540.0/2.;
1064   TVirtualMC::GetMC()->Gsvolu("QTD1", "BOX ", idtmed[7], boxpar, 3);
1065   TVirtualMC::GetMC()->Gspos("QTD1", 1, "Q13TM", -3.8, boxpar[1]+fTDIAperturePos,  0., 0, "ONLY");
1066   boxpar[0] = 11.0/2.;
1067   boxpar[1] = 9.0/2.;
1068   boxpar[2] = 540.0/2.;
1069   TVirtualMC::GetMC()->Gsvolu("QTD2", "BOX ", idtmed[7], boxpar, 3);
1070   TVirtualMC::GetMC()->Gspos("QTD2", 1, "Q13TM", -3.8, -boxpar[1]-fTDIApertureNeg,  0., 0, "ONLY");
1071   boxpar[0] = 5.1/2.;
1072   boxpar[1] = 0.2/2.;
1073   boxpar[2] = 540.0/2.;
1074   TVirtualMC::GetMC()->Gsvolu("QTD3", "BOX ", idtmed[7], boxpar, 3);
1075   TVirtualMC::GetMC()->Gspos("QTD3", 1, "Q13TM", -3.8+5.5+boxpar[0], fTDIAperturePos,  0., 0, "ONLY");
1076   TVirtualMC::GetMC()->Gspos("QTD3", 2, "Q13TM", -3.8+5.5+boxpar[0], -fTDIApertureNeg,  0., 0, "ONLY");
1077   TVirtualMC::GetMC()->Gspos("QTD3", 3, "Q13TM", -3.8-5.5-boxpar[0], fTDIAperturePos,  0., 0, "ONLY");
1078   TVirtualMC::GetMC()->Gspos("QTD3", 4, "Q13TM", -3.8-5.5-boxpar[0], -fTDIApertureNeg,  0., 0, "ONLY");
1079   printf("  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] = 540.0/2.;
1085   tubspar[3] = 90.;
1086   tubspar[4] = 270.;
1087   TVirtualMC::GetMC()->Gsvolu("QTD4", "TUBS", idtmed[7], tubspar, 5);
1088   TVirtualMC::GetMC()->Gspos("QTD4", 1, "Q13TM", -3.8-10.6, 0.,  0., 0, "ONLY");
1089   tubspar[0] = 12.0/2.;
1090   tubspar[1] = 12.4/2.;
1091   tubspar[2] = 540.0/2.;
1092   tubspar[3] = -90.;
1093   tubspar[4] = 90.;
1094   TVirtualMC::GetMC()->Gsvolu("QTD5", "TUBS", idtmed[7], tubspar, 5);
1095   TVirtualMC::GetMC()->Gspos("QTD5", 1, "Q13TM", -3.8+10.6, 0.,  0., 0, "ONLY");
1096   //---------------- END DEFINING TDI INSIDE Q13T -------------------------------
1097
1098   // VCTCG skewed transition piece (ID=332 mm to 212.7 mm) (after TDI)
1099   conpar[0] = (50.0-2.92-1.89)/2.;
1100   conpar[1] = 33.2/2.;
1101   conpar[2] = 33.8/2.;
1102   conpar[3] = 21.27/2.;
1103   conpar[4] = 21.87/2.;
1104   TVirtualMC::GetMC()->Gsvolu("QA11", "CONE", idtmed[7], conpar, 5);
1105   TVirtualMC::GetMC()->Gspos("QA11", 1, "ZDCA", 4.32-3.8, 0., conpar[0]+2.92+zd2, irotpipe5, "ONLY");
1106   // Ch.debug
1107   //printf("    QA11 skewed CONE from z = %f to z =%f (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 ZN left -> %f < z < %f cm\n",fPosZNl[2],fPosZNl[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 ZP left -> %f < z < %f cm\n",fPosZPl[2],fPosZPl[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, ubuf[1], wmat[3], a[3], z[3];
1944
1945   // --- W alloy -> ZN passive material
1946   dens = 17.6;
1947   a[0] = 183.85;
1948   a[1] = 55.85;
1949   a[2] = 58.71;
1950   z[0] = 74.;
1951   z[1] = 26.;
1952   z[2] = 28.;
1953   wmat[0] = .93;
1954   wmat[1] = .03;
1955   wmat[2] = .04;
1956   AliMixture(1, "WALL", a, z, dens, 3, wmat);
1957
1958   // --- Brass (CuZn)  -> ZP passive material
1959   dens = 8.48;
1960   a[0] = 63.546;
1961   a[1] = 65.39;
1962   z[0] = 29.;
1963   z[1] = 30.;
1964   wmat[0] = .63;
1965   wmat[1] = .37;
1966   AliMixture(2, "BRASS", a, z, dens, 2, wmat);
1967
1968   // --- SiO2
1969   dens = 2.64;
1970   a[0] = 28.086;
1971   a[1] = 15.9994;
1972   z[0] = 14.;
1973   z[1] = 8.;
1974   wmat[0] = 1.;
1975   wmat[1] = 2.;
1976   AliMixture(3, "SIO2", a, z, dens, -2, wmat);
1977
1978   // --- Lead
1979   ubuf[0] = 1.12;
1980   AliMaterial(5, "LEAD", 207.19, 82., 11.35, .56, 0., ubuf, 1);
1981
1982   // --- Copper (energy loss taken into account)
1983   ubuf[0] = 1.10;
1984   AliMaterial(6, "COPP0", 63.54, 29., 8.96, 1.4, 0., ubuf, 1);
1985
1986   // --- Copper
1987   ubuf[0] = 1.10;
1988   AliMaterial(9, "COPP1", 63.54, 29., 8.96, 1.4, 0., ubuf, 1);
1989
1990   // --- Iron (energy loss taken into account)
1991   ubuf[0] = 1.1;
1992   AliMaterial(7, "IRON0", 55.85, 26., 7.87, 1.76, 0., ubuf, 1);
1993
1994   // --- Iron (no energy loss)
1995   ubuf[0] = 1.1;
1996   AliMaterial(8, "IRON1", 55.85, 26., 7.87, 1.76, 0., ubuf, 1);
1997
1998   // --- Tatalum
1999   ubuf[0] = 1.1;
2000   AliMaterial(13, "TANT", 183.84, 74., 19.3, 0.35, 0., ubuf, 1);
2001
2002   // ---------------------------------------------------------
2003   Float_t aResGas[3]={1.008,12.0107,15.9994};
2004   Float_t zResGas[3]={1.,6.,8.};
2005   Float_t wResGas[3]={0.28,0.28,0.44};
2006   Float_t dResGas = 3.2E-14;
2007
2008   // --- Vacuum (no magnetic field)
2009   AliMixture(10, "VOID", aResGas, zResGas, dResGas, 3, wResGas);
2010
2011   // --- Vacuum (with magnetic field)
2012   AliMixture(11, "VOIM", aResGas, zResGas, dResGas, 3, wResGas);
2013
2014   // --- Air (no magnetic field)
2015   Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
2016   Float_t zAir[4]={6.,7.,8.,18.};
2017   Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
2018   Float_t dAir = 1.20479E-3;
2019   //
2020   AliMixture(12, "Air    $", aAir, zAir, dAir, 4, wAir);
2021
2022   // ---  Definition of tracking media:
2023
2024   // --- Tantalum = 1 ;
2025   // --- Brass = 2 ;
2026   // --- Fibers (SiO2) = 3 ;
2027   // --- Fibers (SiO2) = 4 ;
2028   // --- Lead = 5 ;
2029   // --- Copper (with high thr.)= 6 ;
2030   // --- Copper (with low thr.)=  9;
2031   // --- Iron (with energy loss) = 7 ;
2032   // --- Iron (without energy loss) = 8 ;
2033   // --- Vacuum (no field) = 10
2034   // --- Vacuum (with field) = 11
2035   // --- Air (no field) = 12
2036
2037   // ****************************************************
2038   //     Tracking media parameters
2039   //
2040   Float_t epsil  = 0.01;   // Tracking precision,
2041   Float_t stmin  = 0.01;   // Min. value 4 max. step (cm)
2042   Float_t stemax = 1.;     // Max. step permitted (cm)
2043   Float_t tmaxfd = 0.;     // Maximum angle due to field (degrees)
2044   Float_t tmaxfdv = 0.1;   // Maximum angle due to field (degrees)
2045   Float_t deemax = -1.;    // Maximum fractional energy loss
2046   Float_t nofieldm = 0.;   // Max. field value (no field)
2047   Float_t fieldm = 45.;    // Max. field value (with field)
2048   Int_t isvol = 0;         // ISVOL =0 -> not sensitive volume
2049   Int_t isvolActive = 1;   // ISVOL =1 -> sensitive volume
2050   Int_t inofld = 0;        // IFIELD=0 -> no magnetic field
2051   Int_t ifield =2;         // IFIELD=2 -> magnetic field defined in AliMagFC.h
2052   // *****************************************************
2053
2054   AliMedium(1, "ZWALL", 1, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2055   AliMedium(2, "ZBRASS",2, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2056   AliMedium(3, "ZSIO2", 3, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2057   AliMedium(4, "ZQUAR", 3, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2058   AliMedium(5, "ZLEAD", 5, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2059   AliMedium(6, "ZCOPP", 6, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2060   AliMedium(7, "ZIRON", 7, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2061   AliMedium(8, "ZIRONN",8, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2062   AliMedium(9, "ZCOPL", 6, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2063   AliMedium(10,"ZVOID",10, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2064   AliMedium(11,"ZVOIM",11, isvol, ifield, fieldm, tmaxfdv, stemax, deemax, epsil, stmin);
2065   AliMedium(12,"ZAIR", 12, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2066   AliMedium(13,"ZTANT",13, isvolActive, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2067   AliMedium(14, "ZIRONT", 7, isvol, inofld, nofieldm, tmaxfd, stemax, deemax, epsil, stmin);
2068
2069 }
2070
2071 //_____________________________________________________________________________
2072 void AliZDCv4::AddAlignableVolumes() const
2073 {
2074  //
2075  // Create entries for alignable volumes associating the symbolic volume
2076  // name with the corresponding volume path. Needs to be syncronized with
2077  // eventual changes in the geometry.
2078  //
2079  if(fOnlyZEM) return;
2080
2081  TString volpath1 = "ALIC_1/ZDCC_1/ZNEU_1";
2082  TString volpath2 = "ALIC_1/ZDCC_1/ZPRO_1";
2083  TString volpath3 = "ALIC_1/ZDCA_1/ZNEU_2";
2084  TString volpath4 = "ALIC_1/ZDCA_1/ZPRO_2";
2085
2086  TString symname1="ZDC/NeutronZDC_C";
2087  TString symname2="ZDC/ProtonZDC_C";
2088  TString symname3="ZDC/NeutronZDC_A";
2089  TString symname4="ZDC/ProtonZDC_A";
2090
2091  if(!gGeoManager->SetAlignableEntry(symname1.Data(),volpath1.Data()))
2092      AliFatal(Form("Alignable entry %s not created. Volume path %s not valid",   symname1.Data(),volpath1.Data()));
2093
2094  if(!gGeoManager->SetAlignableEntry(symname2.Data(),volpath2.Data()))
2095      AliFatal(Form("Alignable entry %s not created. Volume path %s not valid",   symname2.Data(),volpath2.Data()));
2096
2097  if(!gGeoManager->SetAlignableEntry(symname3.Data(),volpath3.Data()))
2098      AliFatal(Form("Alignable entry %s not created. Volume path %s not valid",   symname1.Data(),volpath1.Data()));
2099
2100  if(!gGeoManager->SetAlignableEntry(symname4.Data(),volpath4.Data()))
2101      AliFatal(Form("Alignable entry %s not created. Volume path %s not valid",   symname2.Data(),volpath2.Data()));
2102
2103 }
2104
2105
2106 //_____________________________________________________________________________
2107 void AliZDCv4::Init()
2108 {
2109  InitTables();
2110   Int_t *idtmed = fIdtmed->GetArray();
2111   //
2112   fMedSensZN     = idtmed[1];  // Sensitive volume: ZN passive material
2113   fMedSensZP     = idtmed[2];  // Sensitive volume: ZP passive material
2114   fMedSensF1     = idtmed[3];  // Sensitive volume: fibres type 1
2115   fMedSensF2     = idtmed[4];  // Sensitive volume: fibres type 2
2116   fMedSensZEM    = idtmed[5];  // Sensitive volume: ZEM passive material
2117   fMedSensTDI    = idtmed[6];  // Sensitive volume: TDI Cu shield
2118   fMedSensPI     = idtmed[7];  // Sensitive volume: beam pipes
2119   fMedSensLumi   = idtmed[9];  // Sensitive volume: luminometer
2120   fMedSensGR     = idtmed[12]; // Sensitive volume: air into the grooves
2121   fMedSensVColl  = idtmed[13]; // Sensitive volume: collimator jaws
2122 }
2123
2124 //_____________________________________________________________________________
2125 void AliZDCv4::InitTables()
2126 {
2127  //
2128  // Read light tables for Cerenkov light production parameterization
2129  //
2130
2131   Int_t k, j;
2132   int read=1;
2133
2134   //  --- Reading light tables for ZN
2135   char *lightfName1 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620362207s");
2136   FILE *fp1 = fopen(lightfName1,"r");
2137   if(fp1 == NULL){
2138      printf("Cannot open file fp1 \n");
2139      return;
2140   }
2141   else{
2142     for(k=0; k<fNalfan; k++){
2143       for(j=0; j<fNben; j++){
2144        read = fscanf(fp1,"%f",&fTablen[0][k][j]);
2145        if(read==0) AliDebug(3, " Error in reading light table 1");
2146       }
2147     }
2148     fclose(fp1);
2149   }
2150   char *lightfName2 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620362208s");
2151   FILE *fp2 = fopen(lightfName2,"r");
2152   if(fp2 == NULL){
2153      printf("Cannot open file fp2 \n");
2154      return;
2155   }
2156   else{
2157     for(k=0; k<fNalfan; k++){
2158       for(j=0; j<fNben; j++){
2159        read = fscanf(fp2,"%f",&fTablen[1][k][j]);
2160        if(read==0) AliDebug(3, " Error in reading light table 2");
2161       }
2162     }
2163     fclose(fp2);
2164   }
2165   char *lightfName3 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620362209s");
2166   FILE *fp3 = fopen(lightfName3,"r");
2167   if(fp3 == NULL){
2168      printf("Cannot open file fp3 \n");
2169      return;
2170   }
2171   else{
2172     for(k=0; k<fNalfan; k++){
2173       for(j=0; j<fNben; j++){
2174        read = fscanf(fp3,"%f",&fTablen[2][k][j]);
2175        if(read==0) AliDebug(3, " Error in reading light table 3");
2176       }
2177     }
2178     fclose(fp3);
2179   }
2180   char *lightfName4 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620362210s");
2181   FILE *fp4 = fopen(lightfName4,"r");
2182   if(fp4 == NULL){
2183      printf("Cannot open file fp4 \n");
2184      return;
2185   }
2186   else{
2187     for(k=0; k<fNalfan; k++){
2188       for(j=0; j<fNben; j++){
2189        read = fscanf(fp4,"%f",&fTablen[3][k][j]);
2190        if(read==0) AliDebug(3, " Error in reading light table 4");
2191       }
2192     }
2193     fclose(fp4);
2194   }
2195
2196   //  --- Reading light tables for ZP and ZEM
2197   char *lightfName5 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620552207s");
2198   FILE *fp5 = fopen(lightfName5,"r");
2199   if(fp5 == NULL){
2200      printf("Cannot open file fp5 \n");
2201      return;
2202   }
2203   else{
2204     for(k=0; k<fNalfap; k++){
2205       for(j=0; j<fNbep; j++){
2206        read = fscanf(fp5,"%f",&fTablep[0][k][j]);
2207        if(read==0) AliDebug(3, " Error in reading light table 5");
2208       }
2209     }
2210     fclose(fp5);
2211   }
2212   char *lightfName6 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620552208s");
2213   FILE *fp6 = fopen(lightfName6,"r");
2214   if(fp6 == NULL){
2215      printf("Cannot open file fp6 \n");
2216      return;
2217   }
2218   else{
2219     for(k=0; k<fNalfap; k++){
2220       for(j=0; j<fNbep; j++){
2221        read = fscanf(fp6,"%f",&fTablep[1][k][j]);
2222        if(read==0) AliDebug(3, " Error in reading light table 6");
2223       }
2224     }
2225     fclose(fp6);
2226   }
2227   char *lightfName7 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620552209s");
2228   FILE *fp7 = fopen(lightfName7,"r");
2229   if(fp7 == NULL){
2230      printf("Cannot open file fp7 \n");
2231      return;
2232   }
2233   else{
2234     for(k=0; k<fNalfap; k++){
2235       for(j=0; j<fNbep; j++){
2236        read = fscanf(fp7,"%f",&fTablep[2][k][j]);
2237        if(read==0) AliDebug(3, " Error in reading light table 7");
2238       }
2239     }
2240    fclose(fp7);
2241   }
2242   char *lightfName8 = gSystem->ExpandPathName("$ALICE_ROOT/ZDC/light22620552210s");
2243   FILE *fp8 = fopen(lightfName8,"r");
2244   if(fp8 == NULL){
2245      printf("Cannot open file fp8 \n");
2246      return;
2247   }
2248   else{
2249     for(k=0; k<fNalfap; k++){
2250       for(j=0; j<fNbep; j++){
2251        read = fscanf(fp8,"%f",&fTablep[3][k][j]);
2252        if(read==0) AliDebug(3, " Error in reading light table 8");
2253       }
2254     }
2255    fclose(fp8);
2256   }
2257
2258 }
2259 //_____________________________________________________________________________
2260 void AliZDCv4::StepManager()
2261 {
2262   //
2263   // Routine called at every step in the Zero Degree Calorimeters
2264   //
2265   Int_t   j, vol[2]={0,0}, ibeta=0, ialfa=0, ibe=0, nphe=0;
2266   Float_t hits[14], x[3], xdet[3]={999.,999.,999.}, um[3], ud[3];
2267   Float_t destep=0., be=0., out=0.;
2268   Double_t s[3], p[4];
2269   const char *knamed;
2270   //
2271   for(j=0;j<14;j++) hits[j]=-999.;
2272   //
2273   // --- This part is for no shower developement in beam pipe, TDI, VColl
2274   // If particle interacts with beam pipe, TDI, VColl -> return
2275   if(fNoShower==1 && ((TVirtualMC::GetMC()->CurrentMedium() == fMedSensPI) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensTDI) ||
2276      (TVirtualMC::GetMC()->CurrentMedium() == fMedSensVColl || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensLumi)))){
2277
2278     // If option NoShower is set -> StopTrack
2279
2280     Int_t ipr = 0;
2281       TVirtualMC::GetMC()->TrackPosition(s[0],s[1],s[2]);
2282       if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensPI){
2283         knamed = TVirtualMC::GetMC()->CurrentVolName();
2284         if(!strncmp(knamed,"YMQ",3)){
2285           if(s[2]<0) fpLostITC += 1;
2286           else fpLostITA += 1;
2287           ipr=1;
2288         }
2289         else if(!strncmp(knamed,"YD1",3)){
2290           if(s[2]<0) fpLostD1C += 1;
2291           else fpLostD1A += 1;
2292           ipr=1;
2293         }
2294       }
2295       else if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensTDI){
2296         knamed = TVirtualMC::GetMC()->CurrentVolName();
2297         if(!strncmp(knamed,"MD1",3)){
2298           if(s[2]<0) fpLostD1C += 1;
2299           else  fpLostD1A += 1;
2300           ipr=1;
2301         }
2302         else if(!strncmp(knamed,"QTD",3)) fpLostTDI += 1;
2303       }
2304       else if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensVColl){
2305         knamed = TVirtualMC::GetMC()->CurrentVolName();
2306         if(!strncmp(knamed,"QCVC",4)) fpcVCollC++;
2307         else if(!strncmp(knamed,"QCVA",4))  fpcVCollA++;
2308         ipr=1;
2309       }
2310       //
2311       //TVirtualMC::GetMC()->TrackMomentum(p[0], p[1], p[2], p[3]);
2312       //printf("\t Particle: mass = %1.3f, E = %1.3f GeV, pz = %1.2f GeV -> stopped in volume %s\n",
2313       //     TVirtualMC::GetMC()->TrackMass(), p[3], p[2], TVirtualMC::GetMC()->CurrentVolName());
2314       //
2315       if(ipr<0){
2316         printf("\n\t **********************************\n");
2317         printf("\t ********** Side C **********\n");
2318         printf("\t # of particles in IT = %d\n",fpLostITC);
2319         printf("\t # of particles in D1 = %d\n",fpLostD1C);
2320         printf("\t # of particles in VColl = %d\n",fpcVCollC);
2321         printf("\t ********** Side A **********\n");
2322         printf("\t # of particles in IT = %d\n",fpLostITA);
2323         printf("\t # of particles in D1 = %d\n",fpLostD1A);
2324         printf("\t # of particles in TDI = %d\n",fpLostTDI);
2325         printf("\t # of particles in VColl = %d\n",fpcVCollA);
2326         printf("\t **********************************\n");
2327       }
2328       TVirtualMC::GetMC()->StopTrack();
2329       return;
2330   }
2331
2332   if((TVirtualMC::GetMC()->CurrentMedium() == fMedSensZN) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensZP) ||
2333      (TVirtualMC::GetMC()->CurrentMedium() == fMedSensGR) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensF1) ||
2334      (TVirtualMC::GetMC()->CurrentMedium() == fMedSensF2) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensZEM)){
2335
2336
2337   //Particle coordinates
2338     TVirtualMC::GetMC()->TrackPosition(s[0],s[1],s[2]);
2339     for(j=0; j<=2; j++) x[j] = s[j];
2340     hits[0] = x[0];
2341     hits[1] = x[1];
2342     hits[2] = x[2];
2343
2344   // Determine in which ZDC the particle is
2345     knamed = TVirtualMC::GetMC()->CurrentVolName();
2346     if(!strncmp(knamed,"ZN",2)){
2347           if(x[2]<0.) vol[0]=1; // ZNC (dimuon side)
2348           else if(x[2]>0.) vol[0]=4; //ZNA
2349     }
2350     else if(!strncmp(knamed,"ZP",2)){
2351           if(x[2]<0.) vol[0]=2; //ZPC (dimuon side)
2352           else if(x[2]>0.) vol[0]=5; //ZPA
2353     }
2354     else if(!strncmp(knamed,"ZE",2)) vol[0]=3; //ZEM
2355
2356   // Determine in which quadrant the particle is
2357     if(vol[0]==1){      //Quadrant in ZNC
2358       // Calculating particle coordinates inside ZNC
2359       xdet[0] = x[0]-fPosZNC[0];
2360       xdet[1] = x[1]-fPosZNC[1];
2361       // Calculating quadrant in ZN
2362       if(xdet[0]<=0.){
2363         if(xdet[1]<=0.) vol[1]=1;
2364         else vol[1]=3;
2365       }
2366       else if(xdet[0]>0.){
2367         if(xdet[1]<=0.) vol[1]=2;
2368         else vol[1]=4;
2369       }
2370     }
2371
2372     else if(vol[0]==2){ //Quadrant in ZPC
2373       // Calculating particle coordinates inside ZPC
2374       xdet[0] = x[0]-fPosZPC[0];
2375       xdet[1] = x[1]-fPosZPC[1];
2376       if(xdet[0]>=fDimZP[0])  xdet[0]=fDimZP[0]-0.01;
2377       if(xdet[0]<=-fDimZP[0]) xdet[0]=-fDimZP[0]+0.01;
2378       // Calculating tower in ZP
2379       Float_t xqZP = xdet[0]/(fDimZP[0]/2.);
2380       for(int i=1; i<=4; i++){
2381          if(xqZP>=(i-3) && xqZP<(i-2)){
2382            vol[1] = i;
2383            break;
2384          }
2385       }
2386     }
2387     //
2388     // Quadrant in ZEM: vol[1] = 1 -> particle in 1st ZEM (placed at x = 8.5 cm)
2389     //                  vol[1] = 2 -> particle in 2nd ZEM (placed at x = -8.5 cm)
2390     else if(vol[0] == 3){
2391       if(x[0]>0.){
2392         vol[1] = 1;
2393         // Particle x-coordinate inside ZEM1
2394         xdet[0] = x[0]-fPosZEM[0];
2395       }
2396       else{
2397         vol[1] = 2;
2398         // Particle x-coordinate inside ZEM2
2399         xdet[0] = x[0]+fPosZEM[0];
2400       }
2401       xdet[1] = x[1]-fPosZEM[1];
2402     }
2403     //
2404     else if(vol[0]==4){ //Quadrant in ZNA
2405       // Calculating particle coordinates inside ZNA
2406       xdet[0] = x[0]-fPosZNA[0];
2407       xdet[1] = x[1]-fPosZNA[1];
2408       // Calculating quadrant in ZNA
2409       if(xdet[0]>=0.){
2410         if(xdet[1]<=0.) vol[1]=1;
2411         else vol[1]=3;
2412       }
2413       else if(xdet[0]<0.){
2414         if(xdet[1]<=0.) vol[1]=2;
2415         else vol[1]=4;
2416       }
2417     }
2418     //
2419     else if(vol[0]==5){ //Quadrant in ZPA
2420       // Calculating particle coordinates inside ZPA
2421       xdet[0] = x[0]-fPosZPA[0];
2422       xdet[1] = x[1]-fPosZPA[1];
2423       if(xdet[0]>=fDimZP[0])  xdet[0]=fDimZP[0]-0.01;
2424       if(xdet[0]<=-fDimZP[0]) xdet[0]=-fDimZP[0]+0.01;
2425       // Calculating tower in ZP
2426       Float_t xqZP = -xdet[0]/(fDimZP[0]/2.);
2427       for(int i=1; i<=4; i++){
2428          if(xqZP>=(i-3) && xqZP<(i-2)){
2429            vol[1] = i;
2430            break;
2431          }
2432       }
2433     }
2434     if((vol[1]!=1) && (vol[1]!=2) && (vol[1]!=3) && (vol[1]!=4))
2435       AliError(Form(" WRONG tower for det %d: tow %d with xdet=(%f, %f)\n",
2436                 vol[0], vol[1], xdet[0], xdet[1]));
2437     // Ch. debug
2438     //printf("\t *** det %d vol %d xdet(%f, %f)\n",vol[0], vol[1], xdet[0], xdet[1]);
2439
2440
2441     // Store impact point and kinetic energy of the ENTERING particle
2442
2443     if(TVirtualMC::GetMC()->IsTrackEntering()){
2444       //Particle energy
2445       TVirtualMC::GetMC()->TrackMomentum(p[0],p[1],p[2],p[3]);
2446       hits[3] = p[3];
2447
2448       // Impact point on ZDC
2449       // X takes into account the LHC x-axis sign
2450       // which is opposite to positive x on detector front face
2451       // for side A detectors (ZNA and ZPA)
2452       if(vol[0]==4 || vol[0]==5){
2453         hits[4] = -xdet[0];
2454       }
2455       else{
2456         hits[4] = xdet[0];
2457       }
2458       hits[5] = xdet[1];
2459       hits[6] = 0;
2460       hits[7] = 0;
2461       hits[8] = 0;
2462       hits[9] = 0;
2463       //
2464       Int_t curTrackN = gAlice->GetMCApp()->GetCurrentTrackNumber();
2465       TParticle *part = gAlice->GetMCApp()->Particle(curTrackN);
2466       hits[10] = part->GetPdgCode();
2467       hits[11] = 0;
2468       hits[12] = 1.0e09*TVirtualMC::GetMC()->TrackTime(); // in ns!
2469       hits[13] = part->Eta();
2470       //
2471       if(fFindMother){
2472          Int_t imo = part->GetFirstMother();
2473          //printf(" tracks: pc %d -> mother %d \n", curTrackN,imo);
2474
2475          int trmo = imo;
2476          TParticle *pmot = 0x0;
2477          Bool_t isChild = kFALSE;
2478          if(imo>-1){
2479            pmot = gAlice->GetMCApp()->Particle(imo);
2480            trmo = pmot->GetFirstMother();
2481            isChild = kTRUE;
2482            while(trmo!=-1){
2483               pmot = gAlice->GetMCApp()->Particle(trmo);
2484               //printf("  **** pc %d -> mother %d \n", trch,trmo);
2485               trmo = pmot->GetFirstMother();
2486            }
2487          }
2488
2489          if(isChild && pmot){
2490              hits[6]  = 1;
2491              hits[11] = pmot->GetPdgCode();
2492              hits[13] = pmot->Eta();
2493          }
2494       }
2495
2496
2497       AddHit(curTrackN, vol, hits);
2498
2499       if(fNoShower==1){
2500         if(vol[0]==1){
2501           fnDetectedC += 1;
2502           //if(fnDetectedC==1) printf(" ### Particle in ZNC\n\n");
2503         }
2504         else if(vol[0]==2){
2505           fpDetectedC += 1;
2506           //if(fpDetectedC==1) printf(" ### Particle in ZPC\n\n");
2507         }
2508         //else if(vol[0]==3) printf("   ### Particle in ZEM\n\n");
2509         else if(vol[0]==4){
2510           fnDetectedA += 1;
2511           //if(fnDetectedA==1) printf(" ### Particle in ZNA\n\n");
2512         }
2513         else if(vol[0]==5){
2514           fpDetectedA += 1;
2515           //if(fpDetectedA==1) printf(" ### Particle in ZPA\n\n");
2516         }
2517         //
2518         //printf("\t Pc: x %1.2f y %1.2f z %1.2f  E %1.2f GeV pz = %1.2f GeV in volume %s\n",
2519         //   x[0],x[1],x[3],p[3],p[2],TVirtualMC::GetMC()->CurrentVolName());
2520         //
2521         TVirtualMC::GetMC()->StopTrack();
2522         return;
2523       }
2524     }
2525
2526     // Particle energy loss
2527     if(TVirtualMC::GetMC()->Edep() != 0){
2528       hits[9] = TVirtualMC::GetMC()->Edep();
2529       hits[7] = 0.;
2530       hits[8] = 0.;
2531       AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2532     }
2533   }
2534
2535
2536   // *** Light production in fibres
2537   if((TVirtualMC::GetMC()->CurrentMedium() == fMedSensF1) || (TVirtualMC::GetMC()->CurrentMedium() == fMedSensF2)){
2538
2539      //Select charged particles
2540      if((destep=TVirtualMC::GetMC()->Edep())){
2541
2542        // Particle velocity
2543        Float_t beta = 0.;
2544        TVirtualMC::GetMC()->TrackMomentum(p[0],p[1],p[2],p[3]);
2545        Float_t ptot=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
2546        if(p[3] > 0.00001) beta =  ptot/p[3];
2547        else return;
2548        if(beta<0.67)return;
2549        else if((beta>=0.67) && (beta<=0.75)) ibeta = 0;
2550        else if((beta>0.75)  && (beta<=0.85)) ibeta = 1;
2551        else if((beta>0.85)  && (beta<=0.95)) ibeta = 2;
2552        else if(beta>0.95) ibeta = 3;
2553
2554        // Angle between particle trajectory and fibre axis
2555        // 1 -> Momentum directions
2556        um[0] = p[0]/ptot;
2557        um[1] = p[1]/ptot;
2558        um[2] = p[2]/ptot;
2559        TVirtualMC::GetMC()->Gmtod(um,ud,2);
2560        // 2 -> Angle < limit angle
2561        Double_t alfar = TMath::ACos(ud[2]);
2562        Double_t alfa = alfar*kRaddeg;
2563        if(alfa>=110.) return;
2564        //
2565        ialfa = Int_t(1.+alfa/2.);
2566
2567        // Distance between particle trajectory and fibre axis
2568        TVirtualMC::GetMC()->TrackPosition(s[0],s[1],s[2]);
2569        for(j=0; j<=2; j++){
2570           x[j] = s[j];
2571        }
2572        TVirtualMC::GetMC()->Gmtod(x,xdet,1);
2573        if(TMath::Abs(ud[0])>0.00001){
2574          Float_t dcoeff = ud[1]/ud[0];
2575          be = TMath::Abs((xdet[1]-dcoeff*xdet[0])/TMath::Sqrt(dcoeff*dcoeff+1.));
2576        }
2577        else{
2578          be = TMath::Abs(ud[0]);
2579        }
2580
2581        ibe = Int_t(be*1000.+1);
2582
2583        //Looking into the light tables
2584        Float_t charge = 0.;
2585        Int_t curTrackN = gAlice->GetMCApp()->GetCurrentTrackNumber();
2586        TParticle *part = gAlice->GetMCApp()->Particle(curTrackN);
2587        Int_t pdgCode = part->GetPdgCode();
2588        if(pdgCode<10000) charge = TVirtualMC::GetMC()->TrackCharge();
2589        else{
2590           float z = (pdgCode/10000-100000);
2591           charge = TMath::Abs(z);
2592           //printf(" PDG %d   charge %f\n",pdgCode,charge);
2593        }
2594
2595        if(vol[0]==1 || vol[0]==4) {     // (1)  ZN fibres
2596          if(ibe>fNben) ibe=fNben;
2597          out =  charge*charge*fTablen[ibeta][ialfa][ibe];
2598          nphe = gRandom->Poisson(out);
2599          // Ch. debug
2600          //if(ibeta==3) printf("\t %f \t %f \t %f\n",alfa, be, out);
2601          //printf("\t ibeta = %d, ialfa = %d, ibe = %d -> nphe = %d\n\n",ibeta,ialfa,ibe,nphe);
2602          if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensF1){
2603            hits[7] = nphe;      //fLightPMQ
2604            hits[8] = 0;
2605            hits[9] = 0;
2606            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2607          }
2608          else{
2609            hits[7] = 0;
2610            hits[8] = nphe;      //fLightPMC
2611            hits[9] = 0;
2612            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2613          }
2614        }
2615        else if(vol[0]==2 || vol[0]==5) {// (2) ZP fibres
2616          if(ibe>fNbep) ibe=fNbep;
2617          out =  charge*charge*fTablep[ibeta][ialfa][ibe];
2618          nphe = gRandom->Poisson(out);
2619          if(TVirtualMC::GetMC()->CurrentMedium() == fMedSensF1){
2620            hits[7] = nphe;      //fLightPMQ
2621            hits[8] = 0;
2622            hits[9] = 0;
2623            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2624          }
2625          else{
2626            hits[7] = 0;
2627            hits[8] = nphe;      //fLightPMC
2628            hits[9] = 0;
2629            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2630          }
2631        }
2632        else if(vol[0]==3) {     // (3) ZEM fibres
2633          if(ibe>fNbep) ibe=fNbep;
2634          out =  charge*charge*fTablep[ibeta][ialfa][ibe];
2635          TVirtualMC::GetMC()->TrackPosition(s[0],s[1],s[2]);
2636          Float_t xalic[3];
2637          for(j=0; j<3; j++){
2638             xalic[j] = s[j];
2639          }
2640          // z-coordinate from ZEM front face
2641          // NB-> fPosZEM[2]+fZEMLength = -1000.+2*10.3 = 979.69 cm
2642          Float_t z = -xalic[2]+fPosZEM[2]+2*fZEMLength-xalic[1];
2643          //z = xalic[2]-fPosZEM[2]-fZEMLength-xalic[1]*(TMath::Tan(45.*kDegrad));
2644          //printf("     fPosZEM[2]+2*fZEMLength = %f", fPosZEM[2]+2*fZEMLength);
2645          //
2646          // Parametrization for light guide uniformity
2647          // NEW!!! Light guide tilted @ 51 degrees
2648          Float_t guiPar[4]={0.31,-0.0006305,0.01337,0.8895};
2649          Float_t guiEff = guiPar[0]*(guiPar[1]*z*z+guiPar[2]*z+guiPar[3]);
2650          out = out*guiEff;
2651          nphe = gRandom->Poisson(out);
2652          //printf("     out*guiEff = %f nphe = %d", out, nphe);
2653          if(vol[1] == 1){
2654            hits[7] = 0;
2655            hits[8] = nphe;      //fLightPMC (ZEM1)
2656            hits[9] = 0;
2657            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2658          }
2659          else{
2660            hits[7] = nphe;      //fLightPMQ (ZEM2)
2661            hits[8] = 0;
2662            hits[9] = 0;
2663            AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(), vol, hits);
2664          }
2665        }
2666      }
2667    }
2668 }