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