ZDC/ZDCsim/AliZDCv4.cxx

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