CASTOR/AliCASTOR.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 $Log$
  18 Revision 1.9  2000/12/04 08:48:15  alibrary
  19 Fixing problems in the HEAD
  20
  21 Revision 1.8  2000/10/02 21:28:05  fca
  22 Removal of useless dependecies via forward declarations
  23
  24 Revision 1.7  2000/01/19 17:16:41  fca
  25 Introducing a list of lists of hits -- more hits allowed for detector now
  26
  27 Revision 1.6  1999/09/29 09:24:07  fca
  28 Introduction of the Copyright and cvs Log
  29
  30 */
  31
  32 ///////////////////////////////////////////////////////////////////////////////
  33 //                                                                           //
  34 //  CASTOR                                                                   //
  35 //  This class contains the description of the CASTOR detector               //
  36 //                                                                           //
  37 //Begin_Html
  38 /*
  39 <img src="picts/AliCASTORClass.gif">
  40 </pre>
  41 <br clear=left>
  42 <font size=+2 color=red>
  43 <p>The responsible person for this module is
  44 <a href="mailto:aris.angelis@cern.ch">Aris Angelis</a>.
  45 </font>
  46 <pre>
  47 */
  48 //End_Html
  49 //                                                                           //
  50 //                                                                           //
  51 ///////////////////////////////////////////////////////////////////////////////
  52
  53
  54 #include "AliCASTOR.h"
  55 #include <TNode.h>
  56 #include <TPGON.h>
  57 #include "TGeometry.h"
  58 #include "AliMagF.h"
  59 #include "AliRun.h"
  60 #include "AliMC.h"
  61 #include "AliConst.h"
  62 #include "TMath.h"
  63 static const Double_t kPI=TMath::Pi();
  64
  65 ClassImp(AliCASTOR)
  66
  67 //_____________________________________________________________________________
  68 AliCASTOR::AliCASTOR()
  69 {
  70   //
  71   // Default constructor for CASTOR
  72   //
  73   fIshunt   = 0;
  74 }
  75
  76 //_____________________________________________________________________________
  77 AliCASTOR::AliCASTOR(const char *name, const char *title)
  78        : AliDetector(name,title)
  79 {
  80   //
  81   // Standard constructor for CASTOR
  82   //
  83
  84   //
  85   // Create a tree of castor hits
  86   fHits   = new TClonesArray("AliCASTORhit",  405);
  87   gAlice->AddHitList(fHits);
  88
  89   fIshunt     =  0;
  90
  91   SetMarkerColor(7);
  92   SetMarkerStyle(2);
  93   SetMarkerSize(0.4);
  94 }
  95
  96 //_____________________________________________________________________________
  97 void AliCASTOR::AddHit(Int_t track, Int_t *vol, Float_t *hits)
  98 {
  99   //
 100   // Add a CASTOR hit
 101   //
 102   TClonesArray &lhits = *fHits;
 103   new(lhits[fNhits++]) AliCASTORhit(fIshunt,track,vol,hits);
 104 }
 105
 106 //_____________________________________________________________________________
 107 void AliCASTOR::BuildGeometry()
 108 {
 109   //
 110   // Build CASTOR ROOT TNode geometry for event display
 111   TNode *Node, *Top;
 112   TPGON *pgon;
 113   const int kColorCASTOR  = 4;
 114   //
 115   Top=gAlice->GetGeometry()->GetNode("alice");
 116
 117   // CASTOR
 118   pgon = new TPGON("S_CASTOR","S_CASTOR","void",22.5,360,8,2);
 119   pgon->DefineSection(0,-69.05885,2.598121,12.86874);
 120   pgon->DefineSection(1,69.05885,2.787778,13.88912);
 121   new TRotMatrix("rotcas","rotcas",90,180,90,90,180,0);
 122
 123   Top->cd();
 124   Node = new TNode("CASTOR","CASTOR","S_CASTOR",0,0,-1809.59,"rotcas");
 125   Node->SetLineColor(kColorCASTOR);
 126   fNodes->Add(Node);
 127 }
 128
 129 //_____________________________________________________________________________
 130 Int_t AliCASTOR::DistancetoPrimitive(Int_t , Int_t )
 131 {
 132    return 9999;
 133 }
 134
 135
 136 ClassImp(AliCASTORv1)
 137
 138 //_____________________________________________________________________________
 139 AliCASTORv1::AliCASTORv1() : AliCASTOR()
 140 {
 141   //
 142   // Default constructor for CASTOR version 1
 143   //
 144   fOdFiber = 0;
 145   fOdCladding = 0;
 146   fOdAbsorber = 0;
 147   fOctants = 0;
 148   fLayersEM = 0;
 149   fLayersHad = 0;
 150   fPhiOct = 0;
 151   fRadCore = 0;
 152   fRadFactor = 0;
 153 }
 154
 155 //_____________________________________________________________________________
 156 AliCASTORv1::AliCASTORv1(const char *name, const char *title)
 157        : AliCASTOR(name,title)
 158 {
 159   //
 160   // Standard constructor for CASTOR version 1
 161   //
 162   fOdFiber = 0;
 163   fOdCladding = 0;
 164   fOdAbsorber = 0;
 165   fOctants = 0;
 166   fLayersEM = 0;
 167   fLayersHad = 0;
 168   fPhiOct = 0;
 169   fRadCore = 0;
 170   fRadFactor = 0;
 171 }
 172
 173 //_____________________________________________________________________________
 174 void AliCASTORv1::CreateGeometry()
 175 {
 176   //
 177   // Creation of the geometry of the CASTOR detector
 178   //
 179   //Begin_Html
 180   /*
 181     <img src="picts/AliCASTORTree.gif">
 182   */
 183   //End_Html
 184   //Begin_Html
 185   /*
 186     <img src="picts/AliCASTOR.gif">
 187   */
 188   //End_Html
 189   //
 190   //   28 March 1997   23 February 1998              Aris L. S. Angelis   *
 191   // >--------------------------------------------------------------------<*
 192
 193
 194   Float_t dhad[11], dcal[3], beta, doct[11], alfa1, alfa2, fact1, fact2,fact3;
 195   Float_t dclha[3], dcoha[3], dclem[3], dbxha[3], dcoem[3], dcalt[5], dcalv[5], dbxem[3];
 196   Float_t rzhig;
 197   Float_t s1, s2, s3, rxyin, rzlow, rxyut, facemd, facein, dlayha, dlayem, doctem, doctha, faceut, zendha, phicov;
 198   Float_t doctnt;
 199   Float_t zemhad;
 200   Int_t idrotm[100];
 201   Float_t thecen, xp, xxmdhi, zp, yp, rinbeg;
 202   Float_t rutbeg, xxinhi, rinend, rutend, xxmdlo;
 203   Float_t dztotl, xxinlo, xxuthi;
 204   Float_t xxutlo, dem[11], ang;
 205   Int_t nfx;
 206   Float_t rxy;
 207   // Angle (deg) of inclination of quartz fibres w.r.t. to beam (Cerenkov angle).
 208   const Float_t kBetaD = 45;
 209   //Rapidity range covered by the calorimeter.
 210   const Float_t kEtaLow  = 5.6;
 211   const Float_t kEtaHigh = 7.2;
 212   // Z position (cm) of beginning of calorimeter EM section (the tip.
 213   const Float_t kZbegem = 1740;
 214   // Number of azimuthal calorimeter sectors: octants.
 215   fOctants = 8;
 216   // Number of e-m and hadronic layers (each layer comprises a slice
 217   // of absorber material followed by a slice of active quartz fibres).
 218   //     DATA NLAYEM,NLAYHA /9,69/  ! 0.64 + 9.73 lambda_i
 219   fLayersEM  = 8;
 220   fLayersHad = 72;  // 0.57 + 10.15 lambda_i
 221   // Number of planes of quartz fibres within each active slice for
 222   // e-m and hadronic sections.
 223   const Int_t kFibersEM  = 2;
 224   const Int_t kFibersHad = 4;
 225   // Thickness (cm) of absorber material for e-m and hadronic layers.
 226   const Float_t kAbsorberEM  = 0.5;
 227   const Float_t kAbsorberHad = 1;
 228   // Diameter (cm) of fibre core and of fibre with cladding.
 229   const Float_t kDiamCore     = 0.043;
 230   const Float_t kDiamCladding = 0.045;
 231
 232   Int_t i;
 233   static Int_t debugFlag = fDebug-1;
 234
 235   Int_t *idtmed = fIdtmed->GetArray()-1499;
 236
 237
 238   // >--------------------------------------------------------------------<*
 239   // **> Note: ALICE frame XYZ, proper ref. frame of a trapezoid X'Y'Z'.
 240   // --- Common which contains debug flags for the various detectors ---
 241   // --- Also control flags (JPAWF,JOUTF) for each detector added ---
 242
 243   // **> Common containing some of the Castor FCAL geometry data.
 244
 245   //**> Angle (deg) of inclination of quartz fibres w.r.t. to beam
 246   //**> (Cerenkovangle).
 247   // **> Rapidity range covered by the calorimeter.
 248   // **> Z position (cm) of beginning of calorimeter EM section (the tip.
 249   // **> Number of planes of quartz fibres within each active slice for
 250   // **> e-m and hadronic sections.
 251   // **> Thickness (cm) of absorber material for e-m and hadronic layers.
 252   // **> Diameter (cm) of fibre core and of fibre with cladding.
 253   // **> E-M and hadronic sections of an octant and complete octant module
 254   // **> (general trapezoids).
 255   // **> Imaginary box to hold the complete calorimeter.
 256   // **> Imaginary rectangular boxes containing the trapezoids of the
 257   // **> EM and Hadronic sections of an Octant.
 258   // **> Cylindrical volumes for clad fibres and fibre cores in the
 259   // **> EM and Had sections.
 260   //**> Narrow stainless steel conical beam tube traversing the calorimeter.
 261   // **> Print calorimeter parameters.
 262   // **> Number of azimuthal calorimeter sectors: octants.
 263   //      DATA NOCTS / 16 /
 264   // **> Number of e-m and hadronic layers (each layer comprises a slice
 265   // **> of absorber material followed by a slice of active quartz fibres).
 266   //      DATA NLAYEM,NLAYHA /9,69/  ! 0.64 + 9.73 lambda_i
 267   // 0.57 + 10.15 lambda_i
 268   if (debugFlag > 0) {
 269     printf("----------------------------------\n");
 270     printf(" EtaLo = %f, EtaHigh = %f, ZbegEM =%f\n",kEtaLow, kEtaHigh,kZbegem);
 271     printf(" Nocts =%d, NlayEM=%d, NlayHad = %d\n",fOctants,fLayersEM,fLayersHad);
 272     printf("----------------------------------\n");
 273   }
 274   // **> Radius of sensitive fibre core.
 275   fRadCore = kDiamCore/2;
 276   // **> Radius normalised to radius of 0.5 mm used in the calculation of
 277   // **> the Cherenkov tables.
 278   fRadFactor = fRadCore / .05;
 279   // **> Total number of sensitive QF plane layers.
 280   //nqemly = fLayersEM*kFibersEM;
 281   //nqhaly = fLayersHad*kFibersHad;
 282   beta   = kBetaD*kDegrad; // **> Conversions to radians.
 283   // **> Thickness of e-m and hadronic layers:
 284   // **> Thickness = Thickness_of_Absorber + Thickness_of_N_Fibre_Planes
 285   // **> For N pair: Thickness_of_N_Fibre_Planes = N/2 * [2+TMath::Sqrt(3)]*R_fibre
 286   // **> taking into account staggering of fibres in adjacent planes.
 287   //**> For simplicity staggering not yet introduced, use TMath::Sqrt(4) temporarily.
 288   dlayem = kAbsorberEM +(0.5*kFibersEM )*(2+TMath::Sqrt(4.))*kDiamCladding/2;
 289   dlayha = kAbsorberHad+(0.5*kFibersHad)*(2+TMath::Sqrt(4.))*kDiamCladding/2;
 290   if (debugFlag > 0) {
 291     printf(" Layer Thickness. EM = %f, Had = %f\n",dlayem,dlayha);
 292   }
 293   // **> Thickness of complete octant, along the line perpendicular
 294   // **> to the layers.
 295   // **> Thickness = NlayerEM*DlayerEM + NlayerHad*DlayerHad (DeltaZ').
 296   doctem = fLayersEM*dlayem;
 297   doctha = fLayersHad*dlayha;
 298   doctnt = doctem + doctha;
 299   if (debugFlag > 0) {
 300     printf(" Octant Thickness. EM = %f, Had = %f, Total = %f\n",doctem,doctha,doctnt);
 301   }
 302   // **> Construct one octant module: general trapezoid, rotated such
 303   // **> that the fibre planes are perpenicular to the Z axis of the
 304   // **> proper reference frame (X'Y'Z' frame).
 305   // **> Calculation of the length of the faces at +/- DeltaZ'/2 of an
 306   // **> octant, projected onto the Y'Z' plane (see notes dated 4/4/97).
 307   alfa1 = TMath::ATan(exp(-kEtaLow)) * 2.;
 308   alfa2 = TMath::ATan(exp(-kEtaHigh)) * 2.;
 309   fact1 = (TMath::Tan(alfa1) - TMath::Tan(alfa2)) * TMath::Cos(alfa1) / TMath::Sin(beta - alfa1);
 310   if (debugFlag > 0) {
 311     printf(" Beta =%f,Fact1 =%f\n",kBetaD, fact1);
 312     printf(" EtaLow=%f, EtaHigh=%f, Alfa1=%f, Alfa2=%f\n",kEtaLow,kEtaHigh,alfa1*kRaddeg,alfa2*kRaddeg);
 313   }
 314   // **> Face at entrance to E-M section (-DeltaZ'/2).
 315   facein = fact1 * kZbegem;
 316   // **> Face at interface from E-M to Hadronic section.
 317   facemd = (doctem / TMath::Sin(beta) + kZbegem) * fact1;
 318   // **> Face at exit of Hadronic section (+DeltaZ'/2).
 319   faceut = (doctnt / TMath::Sin(beta) + kZbegem) * fact1;
 320   if (debugFlag > 0) {
 321     printf(" Octant Face Length. Front: %f, Back: %f, EM-Had: %f\n",facein,faceut,facemd);
 322   }
 323   // **> Angular coverage of octant (360./8) projected onto plane
 324   // **> tilted at angle Beta (see notes dated 28/3/97).
 325   //**> PhiTilted = 2*atan[TMath::Tan(phi/2)TMath::Cos(beta)] = 32.65 deg for beta=45,phi=22.5.
 326   fPhiOct = k2PI / fOctants;
 327   phicov = TMath::ATan(TMath::Tan(fPhiOct / 2.) * TMath::Cos(beta)) * 2.;
 328   if (debugFlag > 0) {
 329     printf(" FPhiOct =%f, PhiCov =%f\n",fPhiOct * kRaddeg,phicov * kRaddeg);
 330   }
 331   // **> Dimensions along X' of front and back faces of calorimeter
 332   // **> (see notes dated 8/4/97).
 333   fact2  = TMath::Tan(alfa2) / TMath::Sin(beta);
 334   fact3  = TMath::Cos(alfa2) / TMath::Sin(beta - alfa2);
 335   zendha = doctnt * fact3 + kZbegem;
 336   zemhad = doctem * fact3 + kZbegem;
 337   if (debugFlag > 0) {
 338     printf(" ZbegEM =%f, ZendHA =%f, ZEMHad =%f\n",kZbegem,zendha, zemhad);
 339     printf(" Fact2 =%f, Fact3 =%f\n",fact2,fact3);
 340   }
 341   // **> DeltaX' at -DeltaY'/2, -DeltaZ'/2.
 342   xxinlo = fact2 * 2*kZbegem * TMath::Tan(phicov / 2.);
 343   // **> DeltaX' at +DeltaY'/2, -DeltaZ'/2.
 344   xxinhi = (fact2 + fact1) * 2*kZbegem * TMath::Tan(phicov / 2.);
 345   // **> DeltaX' at -DeltaY'/2, +DeltaZ'/2.
 346   xxutlo = zendha * 2. * fact2 * TMath::Tan(phicov / 2.);
 347   // **> DeltaX' at +DeltaY'/2, +DeltaZ'/2.
 348   xxuthi = zendha * 2. * (fact2 + fact1) * TMath::Tan(phicov / 2.);
 349   // **> DeltaX' at -DeltaY'/2, at EM/Had interface.
 350   xxmdlo = zemhad * 2. * fact2 * TMath::Tan(phicov / 2.);
 351   // **> DeltaX' at +DeltaY'/2, at EM/Had interface.
 352   xxmdhi = zemhad * 2. * (fact2 + fact1) * TMath::Tan(phicov / 2.);
 353   if (debugFlag > 0) {
 354     printf(" XXinLo=%f, XXinHi=%f, XXutLo=%f, XXutHi=%f, XXmdLo=%f, XXmdHi=%f\n",
 355            xxinlo,xxinhi,xxutlo,xxuthi,xxmdlo,xxmdhi);
 356   }
 357   //**> Calculate the polar angle in the X'Y'Z' frame of the line joining the
 358   //**> centres of the front and back faces of the octant (see notes dated 9/4/97).
 359   s1  = (1. - fact2 * TMath::Cos(beta)) * kZbegem;
 360   s2  = (fact2 + fact1 / 2.) * kZbegem;
 361   s3  = TMath::Sqrt(s1 * s1 + s2 * s2 - s1 * s2 * TMath::Cos(kPI - beta));
 362   ang = TMath::ASin(sin(kPI - beta) * s2 / s3);
 363   thecen = kPI/2 - beta + ang;
 364   if (debugFlag > 0) {
 365     printf(" S1=%f, S2=%f, S3=%f, Ang=%f, TheCen=%f\n",s1,s2,s3,ang*kRaddeg,thecen*kRaddeg);
 366   }
 367   // **> Construct the octant volume.
 368   doct[0] = 180*0.125;
 369   doct[1] = 360.;
 370   doct[2] = 8.;
 371   doct[3] = 2.;
 372   doct[4] = -(zendha - kZbegem + faceut * TMath::Cos(beta)) / 2.;
 373   doct[5] = TMath::Tan(alfa2) * kZbegem;
 374   doct[6] = TMath::Tan(alfa1) * kZbegem;
 375   doct[7] = (zendha - kZbegem + faceut * TMath::Cos(beta)) / 2.;
 376   doct[8] = zendha * TMath::Tan(alfa2);
 377   doct[9] = (faceut + zendha * fact2) * TMath::Sin(beta);
 378
 379   if (debugFlag > 0) {
 380     printf("\n Doct(1-10) = ");
 381     for (i = 1; i <= 10; ++i) {
 382       printf("%f, ",doct[i - 1]);
 383     }
 384     printf("   \n");
 385   }
 386   gMC->Gsvolu("OCTA", "PGON", idtmed[fOdAbsorber - 1], doct, 10);
 387   gMC->Gsdvn("OCT ", "OCTA", 8, 2);
 388   // absorber material.
 389   // **> Construct the E-M section volume.
 390   dem[0]  = doctem / 2.;      // DeltaZ'/2
 391   dem[1]  = thecen *kRaddeg;  // Theta[(Centre(-DeltaZ')--Centre(+DeltaZ'
 392   dem[2]  = 90.;              // Phi[(Centre(-DeltaZ')--Centre(+DeltaZ')]
 393   dem[3]  = facein / 2.;      // DeltaY'/2 at -DeltaZ'/2.
 394   dem[4]  = xxinlo / 2.;      // DeltaX'/2 at -DeltaY'/2 at -DeltaZ'/2.
 395   dem[5]  = xxinhi / 2.;      // DeltaX'/2 at +DeltaY'/2 at -DeltaZ'/2.
 396   dem[6]  = 0.;               // Angle w.r.t. Y axis of line joining cent
 397                                 // at +/- DeltaY at -DeltaZ. // Angle w.r.t. Y axis of line joining cent
 398   dem[7]  = facemd / 2.;      // DeltaY'/2 at +DeltaZ'.
 399   dem[8]  = xxmdlo / 2.;      // DeltaX'/2 at -DeltaY'/2 at +DeltaZ'/2.
 400   dem[9]  = xxmdhi / 2.;      // DeltaX'/2 at +DeltaY'/2 at +DeltaZ'/2.
 401   dem[10] = 0.;               // Angle w.r.t. Y axis of line joining cent
 402                                 // at +/- DeltaY at +DeltaZ.
 403
 404   if (debugFlag > 0) {
 405     printf("\n De-m(1-11) =");
 406     for (i = 1; i <= 11; ++i) {
 407       printf("%f, ",dem[i - 1]);
 408     }
 409     printf("   \n");
 410   }
 411   gMC->Gsvolu("EM  ", "TRAP", idtmed[fOdAbsorber - 1], dem, 11);
 412   // absorber material.
 413   // **> Construct the Hadronic section volume.
 414   // Fill with s
 415   dhad[0]  = doctha / 2.;      // DeltaZ'/2
 416   dhad[1]  = thecen *kRaddeg;  // Theta[(Centre(-DeltaZ')--Centre(+DeltaZ'
 417   dhad[2]  = 90.;              // Phi[(Centre(-DeltaZ')--Centre(+DeltaZ')]
 418   dhad[3]  = facemd / 2.;      // DeltaY'/2 at -DeltaZ'/2.
 419   dhad[4]  = xxmdlo / 2.;      // DeltaX'/2 at -DeltaY'/2 at -DeltaZ'/2.
 420   dhad[5]  = xxmdhi / 2.;      // DeltaX'/2 at +DeltaY'/2 at -DeltaZ'/2.
 421   dhad[6]  = 0.;               // Angle w.r.t. Y axis of line joining cent
 422   // at +/- DeltaY at -DeltaZ.
 423   dhad[7]  = faceut / 2.;      // DeltaY'/2 at +DeltaZ'.
 424   dhad[8]  = xxutlo / 2.;      // DeltaX'/2 at -DeltaY'/2 at +DeltaZ'/2.
 425   dhad[9]  = xxuthi / 2.;      // DeltaX'/2 at +DeltaY'/2 at +DeltaZ'/2.
 426   dhad[10] = 0.;               // Angle w.r.t. Y axis of line joining cent
 427   // at +/- DeltaY at +DeltaZ.
 428
 429   if (debugFlag > 0) {
 430     printf("\n Dhad(1-11) = ");
 431     for (i = 1; i <= 11; ++i) {
 432       printf("%f, ",dhad[i - 1]);
 433     }
 434     printf("   \n");
 435   }
 436   gMC->Gsvolu("HAD ", "TRAP", idtmed[fOdAbsorber - 1], dhad, 11); // absorber material.
 437   // **> Rotation matrix to rotate fibres verticaly to fit into holes.
 438   // Fill with
 439   AliMatrix(idrotm[0], 90., 0., 180., 0., 90., 90.);
 440   // **> Internal structure of the EM section starts here.  <---
 441   // **> Construct one sampling module
 442   gMC->Gsdvn("SLEM", "EM  ", fLayersEM, 3);
 443   gMC->Gsatt("SLEM", "SEEN", 0);
 444   // **> Construct the (imaginary) rectangular box embedding the fibres
 445   // **> Fill with air, make it invisible on the drawings.
 446   dbxem[0] = xxmdhi / 2.;
 447   dbxem[2] = kFibersEM*kDiamCladding/2;
 448   dbxem[1] = facemd / 2. + dbxem[2] * TMath::Tan(thecen);
 449   if (debugFlag > 0) {
 450     printf(" DbxEM(1-3) =");
 451     for (i = 1; i <= 3; ++i) {
 452       printf("%f, ",dbxem[i - 1]);
 453     }
 454     printf("   \n");
 455   }
 456   gMC->Gsvolu("BXEM", "BOX ", idtmed[1501], dbxem, 3);
 457   gMC->Gsatt("BXEM", "SEEN", 0);
 458   // **> Divide along Z to obtain one layer
 459   gMC->Gsdvn("RWEM", "BXEM", 2, 3);
 460   gMC->Gsatt("RWEM", "SEEN", 0);
 461   // **> Divide along X' to accomodate the maximum number of individual
 462   //**> fibres packed along X', make the divisions invisible on the drawings.
 463   nfx = Int_t(xxmdhi / .045);
 464   if (debugFlag > 0) {
 465     printf(" NfxEM = %d\n",nfx);
 466   }
 467   gMC->Gsdvn("FXEM", "RWEM", nfx, 1);
 468   gMC->Gsatt("FXEM", "SEEN", 0);
 469   // **> Construct the fiber cladding
 470   dclem[0] = 0.;
 471   dclem[1] = kDiamCladding/2;
 472   dclem[2] = dbxem[1];
 473   if (debugFlag > 0) {
 474     printf(" DclEM(1-3) = \n");
 475     for (i = 1; i <= 3; ++i) {
 476       printf("%f, ",dclem[i - 1]);
 477     }
 478     printf("   \n");
 479   }
 480   gMC->Gsvolu("CLEM", "TUBE", idtmed[fOdCladding - 1], dclem,3);
 481   gMC->Gsatt("CLEM", "SEEN", 0);
 482   //**> Construct the cylindrical volume for a fibre core in the EM section.
 483   //**> Fill with selected fibre material, make it invisible on the drawings.
 484   dcoem[0] = 0.;
 485   dcoem[1] = kDiamCore/2;
 486   dcoem[2] = dbxem[1];
 487   if (debugFlag > 0) {
 488     printf(" DcoEM(1-3) = ");
 489     for (i = 1; i <= 3; ++i) {
 490       printf("%f, ",dcoem[i - 1]);
 491     }
 492     printf("   \n");
 493   }
 494   gMC->Gsvolu("COEM", "TUBE", idtmed[fOdFiber - 1], dcoem,3);
 495   gMC->Gsatt("COEM", "SEEN", 0);
 496   // **> Position the volumes
 497   // **> Put the air section inside one sampling module
 498   // **> Use MANY to obtain clipping of protruding edges.
 499   xp = 0.;
 500   zp = dlayem / 2. - 0.5*kFibersEM*kDiamCladding;
 501   yp = zp * TMath::Tan(thecen);
 502   gMC->Gspos("BXEM", 1, "SLEM", xp, yp, zp, 0, "MANY");
 503   // **> Place the core fibre in the clad
 504   xp = 0.;
 505   yp = 0.;
 506   zp = 0.;
 507   gMC->Gspos("COEM", 1, "CLEM", xp, yp, zp, 0, "MANY");
 508   // **> Put the fiber in its air box
 509   gMC->Gspos("CLEM", 1, "FXEM", xp, yp, zp, idrotm[0], "MANY");
 510   // **> Internal structure of the Hadronic section starts here.  <---
 511   gMC->Gsdvn("SLHA", "HAD ", fLayersHad, 3);
 512   gMC->Gsatt("SLHA", "SEEN", 0);
 513   // **> Construct the air section where the fibers are
 514   dhad[0] = 0.5*kFibersEM*kDiamCladding;
 515   gMC->Gsvolu("AIHA", "TRAP", idtmed[1501], dhad, 11);
 516   // **> Divide along z in the appropriate number of layers
 517   gMC->Gsdvn("SAHA", "AIHA", 4, 3);
 518   //**> Construct the (imaginary) rectangular box embedding one lauer of fibres
 519   // **> Fill with air, make it invisible on the drawings.
 520   dbxha[0] = xxuthi / 2.;
 521   dbxha[2] = 0.5*kFibersHad*kDiamCladding;
 522   dbxha[1] = faceut / 2. + dbxha[2] * TMath::Tan(thecen);
 523   if (debugFlag > 0) {
 524     printf(" DbxHa(1-3) = ");
 525     for (i = 1; i <= 3; ++i) {
 526       printf("%f, ",dbxem[i - 1]);
 527     }
 528     printf("   \n");
 529   }
 530   gMC->Gsvolu("BXHA", "BOX ", idtmed[1501], dbxha, 3);
 531   gMC->Gsatt("BXHA", "SEEN", 0);
 532   // **> Divide along Z to obtain one layer
 533   gMC->Gsdvn("RWHA", "BXHA", 4, 3);
 534   gMC->Gsatt("RWHA", "SEEN", 0);
 535   // **> Divide along X' to accomodate the maximum number of individual
 536   //**> fibres packed along X', make the divisions invisible on the drawings.
 537   nfx = Int_t(xxuthi / .045);
 538   if (debugFlag > 0) {
 539     printf(" NfxHad = %d\n",nfx);
 540   }
 541   gMC->Gsdvn("FXHA", "RWHA", nfx, 1);
 542   gMC->Gsatt("FXHA", "SEEN", 0);
 543   // **> Construct one fiber cladding
 544   dclha[0] = 0.;
 545   dclha[1] = 0.5*kDiamCladding;
 546   dclha[2] = dbxha[1];
 547   if (debugFlag > 0) {
 548     printf(" DclHa(1-3) = ");
 549     for (i = 1; i <= 3; ++i) {
 550       printf("%f, ",dclha[i - 1]);
 551     }
 552     printf("   \n");
 553   }
 554   gMC->Gsvolu("CLHA", "TUBE", idtmed[fOdCladding - 1], dclha,3);
 555   gMC->Gsatt("CLHA", "SEEN", 0);
 556   //**> Construct the cylindrical volume for a fibre core in the Had section.
 557   //**> Fill with selected fibre material, make it invisible on the drawings.
 558   dcoha[0] = 0.;
 559   dcoha[1] = 0.5*kDiamCore;
 560   dcoha[2] = dbxha[1];
 561   if (debugFlag > 0) {
 562     printf(" DcoHa(1-3) = ");
 563     for (i = 1; i <= 3; ++i) {
 564       printf("%f, ",dcoha[i - 1]);
 565     }
 566     printf("   \n");
 567   }
 568   gMC->Gsvolu("COHA", "TUBE", idtmed[fOdFiber - 1], dcoha,3);
 569   gMC->Gsatt("COHA", "SEEN", 0);
 570   // **> Position the volumes
 571   // **> Put the air section inside one sampling module
 572   // **> Use MANY to obtain clipping of protruding edges.
 573   xp = 0.;
 574   zp = dlayha / 2. - 0.5*kFibersHad*kDiamCladding;
 575   yp = zp * TMath::Tan(thecen);
 576   gMC->Gspos("BXHA", 1, "SLHA", xp, yp, zp, 0, "MANY");
 577   // **> Place the core fibre in the clad
 578   xp = 0.;
 579   yp = 0.;
 580   zp = 0.;
 581   gMC->Gspos("COHA", 1, "CLHA", xp, yp, zp, 0, "MANY");
 582   // **> Place the fibre in its air box
 583   gMC->Gspos("CLHA", 1, "FXHA", xp, yp, zp, idrotm[0], "MANY");
 584   // **> Rotation matrices for consecutive calorimeter octants
 585   // **> filling the imaginary box.
 586   AliMatrix(idrotm[1], 90., -90., 45., 0., 45., 180.);
 587   // **> Place the EM and Hadronic sections inside the Octant.
 588   rzlow = (doct[5] + doct[6]) * .5;
 589   rzhig = (doct[8] + doct[9]) * .5;
 590   zp = doct[7] - (faceut * TMath::Cos(beta) + doctha * fact3) * .5;
 591   yp = 0.;
 592   xp = rzlow + (rzhig - rzlow) * .5 * (zp - doct[4]) / doct[7];
 593   gMC->Gspos("HAD ", 1, "OCT ", xp, yp, zp, idrotm[1], "ONLY");
 594   yp = 0.;
 595   zp = doct[7] - faceut * TMath::Cos(beta) * .5 - doctha * fact3 - doctem * fact3 * .5;
 596   xp = rzlow + (rzhig - rzlow) * .5 * (zp - doct[4]) / doct[7];
 597   gMC->Gspos("EM  ", 1, "OCT ", xp, yp, zp, idrotm[1], "ONLY");
 598   // **> An imaginary box to hold the complete calorimeter.
 599   dcal[0] = (faceut + zendha * fact2) * TMath::Sin(beta);
 600   dcal[1] = dcal[0];
 601   dcal[2] = (zendha - kZbegem + faceut * TMath::Cos(beta)) / 2.;
 602   if (debugFlag > 0) {
 603     printf(" Dcal(1-3) = ");
 604     for (i = 1; i <= 3; ++i) {
 605       printf("%f, ",dcal[i - 1]);
 606     }
 607     printf("   \n");
 608   }
 609   gMC->Gsvolu("CAL ", "BOX ", idtmed[1501], dcal, 3);
 610   // Fill with air
 611   rinbeg = TMath::Tan(alfa2) * kZbegem;
 612   rutbeg = TMath::Tan(alfa1) * kZbegem;
 613   dztotl = dcal[2] * 2.;
 614   rinend = (dztotl + kZbegem) * TMath::Tan(alfa2);
 615   rutend = (dztotl + kZbegem) * TMath::Tan(alfa1);
 616   if (debugFlag > 0) {
 617     printf(" RinBeg=%f, RoutBeg=%f\n",rinbeg,rutbeg);
 618     printf(" RinEnd=%f, RoutEnd=%f\n",rinend,rutend);
 619     printf(" DeltaZtotal = %f\n",dztotl);
 620   }
 621   // **> Build the calorimeter inside the imaginary box.
 622   rxyin = (fact2 + fact1 / 2.) * kZbegem; // Radius to centre of octant in X'Y'
 623   // plane at calorimeter entrance.
 624   rxyut = zendha * (fact2 + fact1 / 2.);  // Radius to centre of octant in X'Y'
 625   // plane at calorimeter exit.
 626   rxy   = (rxyin + rxyut) / 2.;           // Radius to geometrical centre of octant in
 627   rxy  *= TMath::Sin(beta);               // projected to the XY plane.
 628   if (debugFlag > 0) {
 629     printf(" \n");
 630   }
 631   gMC->Gspos("OCTA", 1, "CAL ", 0., 0., 0., 0, "ONLY");
 632   //**> Construct the narrow stainless steel conical beam tube traversing the
 633   // **> calorimeter and its vacuum filling:  WallThickness = 0.1 cm,
 634   // **> Router = touching the inner side of the calorimeter,
 635   // **> DeltaZ = all through the calorimeter box.
 636   dcalt[0] = dcal[2];
 637   dcalt[2] = TMath::Tan(alfa2) * kZbegem;
 638   dcalt[1] = dcalt[2] - .1 / TMath::Cos(alfa2);
 639   dcalt[4] = (dcalt[0] * 2. + kZbegem) * TMath::Tan(alfa2);
 640   dcalt[3] = dcalt[4] - .1 / TMath::Cos(alfa2);
 641   dcalv[0] = dcalt[0];
 642   dcalv[2] = dcalt[1];
 643   dcalv[1] = 0.;
 644   dcalv[4] = dcalt[3];
 645   dcalv[3] = 0.;
 646   gMC->Gsvolu("CALT", "CONE", idtmed[1506], dcalt, 5);
 647   // Fe (steel a
 648   gMC->Gsvolu("CALV", "CONE", idtmed[1500], dcalv, 5);
 649   // Vacuum.
 650   gMC->Gsatt("CALV", "SEEN", 0);
 651   // **> Position at centre of calorimeter box.
 652   zp = 0.;
 653   gMC->Gspos("CALT", 1, "CAL ", 0., 0., zp, 0, "ONLY");
 654   gMC->Gspos("CALV", 1, "CAL ", 0., 0., zp, 0, "ONLY");
 655   if (debugFlag > 0) {
 656     printf(" Dcalt,Zp,-/+ = ");
 657     for (i = 1; i <= 5; ++i) {
 658       printf("%f, ",dcalt[i - 1]);
 659     }
 660     printf("%f, %f, %f\n",zp, zp - dcalt[0], zp + dcalt[0]);
 661     printf(" Dcalt,Zp,-/+ = ");
 662     for (i = 1; i <= 5; ++i) {
 663       printf("%f, ",dcalt[i - 1]);
 664     }
 665     printf("%f, %f, %f\n",zp, zp - dcalt[0], zp + dcalt[0]);
 666   }
 667   // **> Rotate the imaginary box carrying the calorimeter and place it
 668   // **> in the ALICE volume on the -Z side.
 669   xp = 0.;
 670   yp = 0.;
 671   zp = dcal[2] + kZbegem;
 672   AliMatrix(idrotm[2], 90., 180., 90., 90., 180., 0.);
 673   // -X theta and phi w.r.t. to box XYZ.
 674   //  Y theta and phi w.r.t. to box XYZ.
 675   // -Z theta and phi w.r.t. to box XYZ.
 676   gMC->Gspos("CAL ", 1, "ALIC", xp, yp, -zp, idrotm[2], "ONLY");
 677   if (debugFlag > 0) {
 678     printf(" Dcal,Zp,-/+ = ");
 679     for (i = 1; i <= 3; ++i) {
 680       printf("%f, ",dcal[i - 1]);
 681     }
 682     printf("%f, %f, %f\n",zp, zp - dcal[2], zp + dcal[2]);
 683   }
 684 }
 685
 686 //_____________________________________________________________________________
 687 void AliCASTORv1::DrawModule()
 688 {
 689   //
 690   // Draw a shaded view of CASTOR version 1
 691   //
 692
 693
 694   gMC->Gsatt("*", "seen", -1);
 695   gMC->Gsatt("alic", "seen", 0);
 696   //
 697   // Set visibility of elements
 698   gMC->Gsatt("OCTA","seen",0);
 699   gMC->Gsatt("EM  ","seen",0);
 700   gMC->Gsatt("HAD ","seen",0);
 701   gMC->Gsatt("CAL ","seen",0);
 702   gMC->Gsatt("CALT","seen",1);
 703   gMC->Gsatt("OCT ","seen",0);
 704   gMC->Gsatt("SLEM","seen",1);
 705   gMC->Gsatt("SLHA","seen",1);
 706   gMC->Gsatt("SAHA","seen",1);
 707   //
 708   gMC->Gdopt("hide", "on");
 709   gMC->Gdopt("shad", "on");
 710   gMC->Gsatt("*", "fill", 7);
 711   gMC->SetClipBox(".");
 712   gMC->SetClipBox("*", 0, 20, -20, 20, -1900, -1700);
 713   gMC->DefaultRange();
 714   gMC->Gdraw("alic", 40, 30, 0, -191.5, -78, .19, .19);
 715   gMC->Gdhead(1111, "CASTOR Version 1");
 716   gMC->Gdman(15,-2, "MAN");
 717   gMC->Gdopt("hide", "off");
 718 }
 719
 720 //_____________________________________________________________________________
 721 void AliCASTORv1::CreateMaterials()
 722 {
 723   //
 724   // Create materials for CASTOR version 1
 725   //
 726   //   30 March 1997   27 November 1997              Aris L. S. Angelis   *
 727   // >--------------------------------------------------------------------<*
 728   Int_t   ISXFLD = gAlice->Field()->Integ();
 729   Float_t SXMGMX = gAlice->Field()->Max();
 730
 731   Int_t *idtmed = fIdtmed->GetArray()-1499;
 732
 733   Float_t cute, ubuf[1], cutg, epsil, awmix[3], dwmix, stmin;
 734   Int_t isvol;
 735   Float_t wwmix[3], zwmix[3], aq[2], dq, zq[2], wq[2];
 736   Float_t tmaxfd, stemax, deemax;
 737   Int_t kod;
 738
 739
 740   // **> Quartz and Wmixture.
 741   // **> UBUF is the value of r0, used for calculation of the radii of
 742   // **> the nuclei and the Woods-Saxon potential.
 743   ubuf[0] = .68;
 744   AliMaterial(1, "Vacuum$", 1e-16, 1e-16, 1e-16, 1e16, 1e16, ubuf, 1);
 745   ubuf[0] = .68;
 746   AliMaterial(2, "Air   $", 14.61, 7.3, .001205, 30420., 67500., ubuf, 1);
 747   //**> Quartz (SiO2) and fluorinated (?) quartz for cladding (insensitive).
 748   dq    = 2.64;
 749   aq[0] = 28.086;
 750   aq[1] = 15.9994;
 751   zq[0] = 14.;
 752   zq[1] = 8.;
 753   wq[0] = 1.;
 754   wq[1] = 2.;
 755   AliMixture(3, "Quartz$", aq, zq, dq, -2, wq);
 756   // After a call with ratios by number (negative number of elements),
 757   // the ratio array is changed to the ratio by weight, so all successive
 758   // calls with the same array must specify the number of elements as
 759   // positive
 760   AliMixture(4, "FQuartz$", aq, zq, dq, 2, wq);
 761   // **> W mixture (90% W + 7.5% Ni + 2.5% Cu).
 762   awmix[0] = 183.85;
 763   zwmix[0] = 74.;
 764   wwmix[0] = .9;
 765   awmix[1] = 58.69;
 766   zwmix[1] = 28.;
 767   wwmix[1] = .075;
 768   awmix[2] = 63.55;
 769   zwmix[2] = 29.;
 770   wwmix[2] = .025;
 771   dwmix    = 17.2;
 772   // **> (Pure W and W mixture are given the same material number
 773   // **> so that they can be used interchangeably).
 774   ubuf[0] = 1.1;
 775   AliMixture(5, "W Mix $", awmix, zwmix, dwmix, 3, wwmix);
 776   // **> Lead.
 777   ubuf[0] = 1.12;
 778   AliMaterial(6, "Pb208 $", 207.19, 82., 11.35, .56, 18.5, ubuf, 1);
 779   // **> Iron.
 780   ubuf[0] = .99;
 781   AliMaterial(7, "Fe56  $", 55.85, 26., 7.87, 1.76, 16.7, ubuf, 1);
 782   // **> Copper.
 783   ubuf[0] = 1.01;
 784   AliMaterial(8, "Cu63  $", 63.54, 29., 8.96, 1.43, 15., ubuf, 1);
 785   // **> Debug Printout.
 786   //      CALL GPRINT('MATE',0)
 787   // **> (Negative values for automatic calculation in case of AUTO=0).
 788   isvol  = 0;    // Sensitive volume flag.
 789   tmaxfd = .1;   // Max allowed angular deviation in 1 step due to field
 790   stemax = -.5;  // Maximum permitted step size (cm).
 791   deemax = -.2;  // Maximum permitted fractional energy loss.
 792   epsil  = .01;  // Boundary crossing precision (cm).
 793   stmin  = -.1;  // Minimum permitted step size inside absorber (cm).
 794   AliMedium(1, "Vacuum$", 1, isvol, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 795   AliMedium(2, "Air   $", 2, isvol, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 796
 797   // **> Options for Cherenkov fibres and cladding.
 798   isvol = 1;    // Declare fibre core as sensitive.
 799   AliMedium(3, "Quartz$", 3, isvol, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 800   isvol = 0;    // Declare fibre cladding as not sensitive.
 801   AliMedium(4, "FQuartz$", 4, isvol, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 802
 803   // **> Options for absorber material (not sensitive).
 804   isvol  = 0;   // Sensitive volume flag.
 805   stemax = .5;  // Maximum permitted step size (cm).
 806   deemax = .5;  // Maximum permitted fractional energy loss.
 807   stmin  = .1;  // Minimum permitted step size inside absorber (cm).
 808   AliMedium(5, "W Mix $",  5, isvol, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 809   AliMedium(6, "Pb208 $",  6, isvol, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 810   AliMedium(7, "Fe56  $ ", 7, isvol, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 811   AliMedium(8, "Cu63  $ ", 8, isvol, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
 812
 813   // **> Select material for the Cherenkov fibres.
 814   fOdFiber    = 1503;
 815   //      CALL GPTMED(IDTMED(KODFBR))
 816   // **> Select material for the fibre cladding.
 817   // Quartz.
 818   fOdCladding = 1504;
 819   //      CALL GPTMED(IDTMED(KODCLD))
 820   // **> Select absorber material.
 821   // FQuartz.
 822   fOdAbsorber = 1505;  // W184/Mix
 823   //      KODABS=1506   ! Pb208.
 824   //      KODABS=1507   ! Fe56.
 825   //      KODABS=1508   ! Cu63.
 826   //      CALL GPTMED(IDTMED(KODABS))
 827   // **> Set by default all interactions and decays explicitly ON
 828   // **> and redefine the kinetic energy cutoffs:
 829   //      CUTE=0.0031       ! Allow beta >= 0.99 only.
 830   cute = 7e-4;  // Allow beta >= 0.67 only.
 831   cutg = cute * 1.33;
 832
 833   // **> Inside the absorber material,
 834   for (kod = 1505; kod <= 1508; ++kod) {
 835     Int_t absorber = idtmed[kod - 1];
 836     gMC->Gstpar(absorber, "CUTELE", cute);  // Allow beta >= 0.xx
 837     gMC->Gstpar(absorber, "CUTGAM", cutg);  // = 1.33 cutele.
 838     gMC->Gstpar(absorber, "CUTNEU", .01);   // Default.
 839     gMC->Gstpar(absorber, "CUTHAD", .01);   // Default.
 840     gMC->Gstpar(absorber, "CUTMUO", .01);   // Default.
 841     gMC->Gstpar(absorber, "BCUTE", cutg);   // = cutgam.
 842     gMC->Gstpar(absorber, "BCUTM", cutg);   // = cutgam.
 843     gMC->Gstpar(absorber, "DCUTE", cute);   // = cutele.
 844     gMC->Gstpar(absorber, "DCUTM", cute);   // = cutele.
 845     gMC->Gstpar(absorber, "PPCUTM", cutg);  // = 1.33 cutele.
 846     gMC->Gstpar(absorber, "DCAY", 1.);
 847     gMC->Gstpar(absorber, "MULS", 1.);
 848     gMC->Gstpar(absorber, "PFIS", 1.);
 849     gMC->Gstpar(absorber, "MUNU", 1.);
 850     gMC->Gstpar(absorber, "LOSS", 1.);
 851     gMC->Gstpar(absorber, "PHOT", 1.);
 852     gMC->Gstpar(absorber, "COMP", 1.);
 853     gMC->Gstpar(absorber, "PAIR", 1.);
 854     gMC->Gstpar(absorber, "BREM", 1.);
 855     gMC->Gstpar(absorber, "RAYL", 1.);
 856     gMC->Gstpar(absorber, "DRAY", 1.);
 857     gMC->Gstpar(absorber, "ANNI", 1.);
 858     gMC->Gstpar(absorber, "HADR", 1.);
 859     gMC->Gstpar(absorber, "LABS", 1.);
 860   }
 861   // **> Inside the cladding,
 862   Int_t cladding = idtmed[fOdCladding - 1];
 863   gMC->Gstpar(cladding, "CUTELE", cute);  // Allow beta >= 0.xx
 864   gMC->Gstpar(cladding, "CUTGAM", cutg);  // = 1.33 cutele.
 865   gMC->Gstpar(cladding, "CUTNEU", .01);   // Default.
 866   gMC->Gstpar(cladding, "CUTHAD", .01);   // Default.
 867   gMC->Gstpar(cladding, "CUTMUO", .01);   // Default.
 868   gMC->Gstpar(cladding, "BCUTE", cutg);   // = cutgam.
 869   gMC->Gstpar(cladding, "BCUTM", cutg);   // = cutgam.
 870   gMC->Gstpar(cladding, "DCUTE", cute);   // = cutele.
 871   gMC->Gstpar(cladding, "DCUTM", cute);   // = cutele.
 872   gMC->Gstpar(cladding, "PPCUTM", cutg);  // = 1.33 cutele.
 873   gMC->Gstpar(cladding, "DCAY", 1.);
 874   gMC->Gstpar(cladding, "MULS", 1.);
 875   gMC->Gstpar(cladding, "PFIS", 1.);
 876   gMC->Gstpar(cladding, "MUNU", 1.);
 877   gMC->Gstpar(cladding, "LOSS", 1.);
 878   gMC->Gstpar(cladding, "PHOT", 1.);
 879   gMC->Gstpar(cladding, "COMP", 1.);
 880   gMC->Gstpar(cladding, "PAIR", 1.);
 881   gMC->Gstpar(cladding, "BREM", 1.);
 882   gMC->Gstpar(cladding, "RAYL", 1.);
 883   gMC->Gstpar(cladding, "DRAY", 1.);
 884   gMC->Gstpar(cladding, "ANNI", 1.);
 885   gMC->Gstpar(cladding, "HADR", 1.);
 886   gMC->Gstpar(cladding, "LABS", 1.);
 887
 888   // **> and Inside the Cherenkov fibres,
 889   Int_t fiber = idtmed[fOdFiber - 1];
 890   gMC->Gstpar(fiber, "CUTELE", cute);  // Allow beta >= 0.xx
 891   gMC->Gstpar(fiber, "CUTGAM", cutg);  // = 1.33 cutele.
 892   gMC->Gstpar(fiber, "CUTNEU", .01);   // Default.
 893   gMC->Gstpar(fiber, "CUTHAD", .01);   // Default.
 894   gMC->Gstpar(fiber, "CUTMUO", .01);   // Default.
 895   gMC->Gstpar(fiber, "BCUTE", cutg);   // = cutgam.
 896   gMC->Gstpar(fiber, "BCUTM", cutg);   // = cutgam.
 897   gMC->Gstpar(fiber, "DCUTE", cute);   // = cutele.
 898   gMC->Gstpar(fiber, "DCUTM", cute);   // = cutele.
 899   gMC->Gstpar(fiber, "PPCUTM", cutg);  // = 1.33 cutele.
 900   gMC->Gstpar(fiber, "DCAY", 1.);
 901   gMC->Gstpar(fiber, "MULS", 1.);
 902   gMC->Gstpar(fiber, "PFIS", 1.);
 903   gMC->Gstpar(fiber, "MUNU", 1.);
 904   gMC->Gstpar(fiber, "LOSS", 1.);
 905   gMC->Gstpar(fiber, "PHOT", 1.);
 906   gMC->Gstpar(fiber, "COMP", 1.);
 907   gMC->Gstpar(fiber, "PAIR", 1.);
 908   gMC->Gstpar(fiber, "BREM", 1.);
 909   gMC->Gstpar(fiber, "RAYL", 1.);
 910   gMC->Gstpar(fiber, "DRAY", 1.);
 911   gMC->Gstpar(fiber, "ANNI", 1.);
 912   gMC->Gstpar(fiber, "HADR", 1.);
 913   gMC->Gstpar(fiber, "LABS", 1.);
 914 }
 915
 916 //_____________________________________________________________________________
 917 void AliCASTORv1::StepManager()
 918 {
 919   //
 920   // Called at every step in CASTOR
 921   //
 922 }
 923
 924 //_____________________________________________________________________________
 925 void AliCASTORv1::Init()
 926 {
 927   //
 928   // Initialise CASTOR detector after it has been built
 929   //
 930   Int_t i;
 931   //
 932   if(fDebug) {
 933     printf("\n%s: ",ClassName());
 934     for(i=0;i<35;i++) printf("*");
 935     printf(" CASTOR_INIT ");
 936     for(i=0;i<35;i++) printf("*");
 937     printf("\n%s: ",ClassName());
 938     //
 939     // Here the ABSO initialisation code (if any!)
 940     for(i=0;i<80;i++) printf("*");
 941     printf("\n");
 942   }
 943 }
 944
 945 ClassImp(AliCASTORhit)
 946
 947 //_____________________________________________________________________________
 948 AliCASTORhit::AliCASTORhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
 949 AliHit(shunt, track)
 950 {
 951   //
 952   // Store a CASTOR hit
 953   //
 954   fVolume  = vol[0];
 955   fX=hits[0];
 956   fY=hits[1];
 957   fZ=hits[2];
 958 }
 959
 960