CASTOR/AliCASTOR.cxx

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