4f67230787836cf3f1797434a5469761f359f8b1
[u/mrichter/AliRoot.git] / RICH / AliRICHv0.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.8  2000/05/26 17:30:08  jbarbosa
19   Cerenkov angle now stored within cerenkov data structure.
20
21   Revision 1.7  2000/05/18 10:31:36  jbarbosa
22   Fixed positioning of spacers inside freon.
23   Fixed positioning of proximity gap
24   inside methane.
25   Fixed cut on neutral particles in the StepManager.
26
27   Revision 1.6  2000/04/28 11:51:58  morsch
28    Dimensions of arrays hits and Ckov_data corrected.
29
30   Revision 1.5  2000/04/19 13:28:46  morsch
31   Major changes in geometry (parametrised), materials (updated) and
32   step manager (diagnostics) (JB, AM)
33
34 */
35
36
37
38 ////////////////////////////////////////////////////////
39 //  Manager and hits classes for set:RICH version 0    //
40 /////////////////////////////////////////////////////////
41
42 #include <TTUBE.h>
43 #include <TNode.h> 
44 #include <TRandom.h> 
45
46 #include "AliRICHv0.h"
47 #include "AliRun.h"
48 #include "AliMC.h"
49 #include "iostream.h"
50 #include "AliCallf77.h"
51 #include "AliConst.h" 
52 #include "AliPDG.h" 
53 #include "TGeant3.h"
54
55 ClassImp(AliRICHv0)
56     
57 //___________________________________________
58 AliRICHv0::AliRICHv0() : AliRICH()
59 {
60     //fChambers = 0;
61 }
62
63 //___________________________________________
64 AliRICHv0::AliRICHv0(const char *name, const char *title)
65     : AliRICH(name,title)
66 {
67     fCkov_number=0;
68     fFreon_prod=0;
69
70     fChambers = new TObjArray(7);
71     for (Int_t i=0; i<7; i++) {
72     
73         (*fChambers)[i] = new AliRICHChamber();  
74         
75     }
76 }
77
78
79 //___________________________________________
80 void AliRICHv0::CreateGeometry()
81 {
82     //
83     // Create the geometry for RICH version 1
84     //
85     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it) 
86     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it) 
87     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it) 
88     //
89     //Begin_Html
90     /*
91       <img src="picts/AliRICHv1.gif">
92     */
93     //End_Html
94     //Begin_Html
95     /*
96       <img src="picts/AliRICHv1Tree.gif">
97     */
98     //End_Html
99
100   AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH"); 
101   AliRICHSegmentation*  segmentation;
102   AliRICHGeometry*  geometry;
103   AliRICHChamber*       iChamber;
104
105   iChamber = &(RICH->Chamber(0));
106   segmentation=iChamber->GetSegmentationModel(0);
107   geometry=iChamber->GetGeometryModel();
108
109   Float_t distance;
110   distance = geometry->GetFreonThickness()/2 + geometry->GetQuartzThickness() + geometry->GetGapThickness();
111   geometry->SetRadiatorToPads(distance);
112     
113     
114     Int_t *idtmed = fIdtmed->GetArray()-999;
115     
116     Int_t i;
117     Float_t zs;
118     Int_t idrotm[1099];
119     Float_t par[3];
120     
121     // --- Define the RICH detector 
122     //     External aluminium box 
123     par[0] = 71.1;
124     par[1] = 11.5;                 //Original Settings
125     par[2] = 73.15;
126     /*par[0] = 73.15;
127     par[1] = 11.5;
128     par[2] = 71.1;*/
129     gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
130     
131     //     Sensitive part of the whole RICH 
132     par[0] = 64.8;
133     par[1] = 11.5;                 //Original Settings
134     par[2] = 66.55;
135     /*par[0] = 66.55;
136     par[1] = 11.5;
137     par[2] = 64.8;*/
138     gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
139     
140     //     Honeycomb 
141     par[0] = 63.1;
142     par[1] = .188;                 //Original Settings
143     par[2] = 66.55;
144     /*par[0] = 66.55;
145     par[1] = .188;
146     par[2] = 63.1;*/
147     gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
148     
149     //     Aluminium sheet 
150     par[0] = 63.1;
151     par[1] = .025;                 //Original Settings
152     par[2] = 66.55;
153     /*par[0] = 66.5;
154     par[1] = .025;
155     par[2] = 63.1;*/
156     gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
157     
158     //     Quartz 
159     par[0] = geometry->GetQuartzWidth()/2;
160     par[1] = geometry->GetQuartzThickness()/2;
161     par[2] = geometry->GetQuartzLength()/2;
162     /*par[0] = 63.1;
163     par[1] = .25;                  //Original Settings
164     par[2] = 65.5;*/
165     /*par[0] = geometry->GetQuartzWidth()/2;
166     par[1] = geometry->GetQuartzThickness()/2;
167     par[2] = geometry->GetQuartzLength()/2;*/
168     //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f %f %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[0],par[1],par[2]);
169     gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
170     
171     //     Spacers (cylinders) 
172     par[0] = 0.;
173     par[1] = .5;
174     par[2] = geometry->GetFreonThickness()/2;
175     gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
176     
177     //     Opaque quartz 
178     par[0] = 61.95;
179     par[1] = .2;                   //Original Settings
180     par[2] = 66.5;
181     /*par[0] = 66.5;
182     par[1] = .2;
183     par[2] = 61.95;*/
184     gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
185   
186     //     Frame of opaque quartz
187     par[0] = geometry->GetOuterFreonWidth()/2;
188     par[1] = geometry->GetFreonThickness()/2;
189     par[2] = geometry->GetOuterFreonLength()/2 + 1; 
190     /*par[0] = 20.65;
191     par[1] = .5;                   //Original Settings
192     par[2] = 66.5;*/
193     /*par[0] = 66.5;
194     par[1] = .5;
195     par[2] = 20.65;*/
196     gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
197
198     par[0] = geometry->GetInnerFreonWidth()/2;
199     par[1] = geometry->GetFreonThickness()/2;
200     par[2] = geometry->GetInnerFreonLength()/2 + 1; 
201     gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
202     
203     //     Little bar of opaque quartz 
204     par[0] = .275;
205     par[1] = geometry->GetQuartzThickness()/2;
206     par[2] = geometry->GetInnerFreonLength()/2 - 2.4; 
207     /*par[0] = .275;
208     par[1] = .25;                   //Original Settings
209     par[2] = 63.1;*/
210     /*par[0] = 63.1;
211     par[1] = .25;
212     par[2] = .275;*/
213     gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
214     
215     //     Freon 
216     par[0] = geometry->GetOuterFreonWidth()/2;
217     par[1] = geometry->GetFreonThickness()/2;
218     par[2] = geometry->GetOuterFreonLength()/2; 
219     /*par[0] = 20.15;
220     par[1] = .5;                   //Original Settings
221     par[2] = 65.5;*/
222     /*par[0] = 65.5;
223     par[1] = .5;
224     par[2] = 20.15;*/
225     gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
226
227     par[0] = geometry->GetInnerFreonWidth()/2;
228     par[1] = geometry->GetFreonThickness()/2;
229     par[2] = geometry->GetInnerFreonLength()/2; 
230     gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
231     
232     //     Methane 
233     par[0] = 64.8;
234     par[1] = geometry->GetGapThickness()/2;
235     //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
236     par[2] = 64.8;
237     gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
238     
239     //     Methane gap 
240     par[0] = 64.8;
241     par[1] = geometry->GetProximityGapThickness()/2;
242     //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
243     par[2] = 64.8;
244     gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
245     
246     //     CsI photocathode 
247     par[0] = 64.8;
248     par[1] = .25;
249     par[2] = 64.8;
250     gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
251     
252     //     Anode grid 
253     par[0] = 0.;
254     par[1] = .001;
255     par[2] = 20.;
256     gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
257     
258     // --- Places the detectors defined with GSVOLU 
259     //     Place material inside RICH 
260     gMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
261     
262     gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 -.05 - .376 -.025, 0., 0, "ONLY");
263     gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2  - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 -.05 - .188, 0., 0, "ONLY");
264     gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .025, 0., 0, "ONLY");
265     gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
266     
267     AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
268     
269     Int_t nspacers = (Int_t)(TMath::Abs(geometry->GetInnerFreonLength()/14.4));
270     //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n Spacers:%d\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n",nspacers); 
271
272     printf("Nspacers: %d", nspacers);
273     
274     //for (i = 1; i <= 9; ++i) {
275       //zs = (5 - i) * 14.4;                       //Original settings 
276     for (i = 0; i < nspacers; i++) {
277         zs = (TMath::Abs(nspacers/2) - i) * 14.4;
278         gMC->Gspos("SPAC", i, "FRE1", 6.7, 0., zs, idrotm[1019], "ONLY");  //Original settings 
279         //gMC->Gspos("SPAC", i, "FRE1", zs, 0., 6.7, idrotm[1019], "ONLY"); 
280     }
281     //for (i = 10; i <= 18; ++i) {
282       //zs = (14 - i) * 14.4;                       //Original settings 
283     for (i = nspacers; i < nspacers*2; ++i) {
284         zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4;
285         gMC->Gspos("SPAC", i, "FRE1", -6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings  
286         //gMC->Gspos("SPAC", i, "FRE1", zs, 0., -6.7, idrotm[1019], "ONLY");  
287     }
288
289     //for (i = 1; i <= 9; ++i) {
290       //zs = (5 - i) * 14.4;                       //Original settings 
291       for (i = 0; i < nspacers; i++) {
292         zs = (TMath::Abs(nspacers/2) - i) * 14.4;
293         gMC->Gspos("SPAC", i, "FRE2", 6.7, 0., zs, idrotm[1019], "ONLY");  //Original settings 
294         //gMC->Gspos("SPAC", i, "FRE2", zs, 0., 6.7, idrotm[1019], "ONLY");
295     }
296     //for (i = 10; i <= 18; ++i) {
297       //zs = (5 - i) * 14.4;                       //Original settings 
298       for (i = nspacers; i < nspacers*2; ++i) {
299         zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4;
300         gMC->Gspos("SPAC", i, "FRE2", -6.7, 0., zs, idrotm[1019], "ONLY");  //Original settings 
301         //gMC->Gspos("SPAC", i, "FRE2", zs, 0., -6.7, idrotm[1019], "ONLY");
302     }
303     
304     /*gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
305     gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
306     gMC->Gspos("OQF1", 1, "SRIC", 31.3, -4.724, 41.3, 0, "ONLY");
307     gMC->Gspos("OQF2", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
308     gMC->Gspos("OQF1", 3, "SRIC", -31.3, -4.724, -41.3, 0, "ONLY");
309     gMC->Gspos("BARR", 1, "QUAR", -21.65, 0., 0., 0, "ONLY");           //Original settings 
310     gMC->Gspos("BARR", 2, "QUAR", 21.65, 0., 0., 0, "ONLY");            //Original settings 
311     gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
312     gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
313     gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
314     gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");*/
315
316
317     gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
318     gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
319     gMC->Gspos("OQF1", 1, "SRIC", geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2, 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3)
320     gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY");          //Original settings 
321     gMC->Gspos("OQF1", 3, "SRIC", - (geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2), 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY");       //Original settings (-31.3)
322     gMC->Gspos("BARR", 1, "QUAR", -21.65, 0., 0., 0, "ONLY");           //Original settings 
323     gMC->Gspos("BARR", 2, "QUAR", 21.65, 0., 0., 0, "ONLY");            //Original settings 
324     gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
325     gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY");
326     gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
327     gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY");
328
329     printf("Position of the gap: %f to %f\n", 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - .2, 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 + .2);
330     
331     //     Place RICH inside ALICE apparatus 
332   
333     AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
334     AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
335     AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
336     AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
337     AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
338     AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
339     AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
340     
341     gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26,     idrotm[1000], "ONLY");
342     gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0.,        idrotm[1001], "ONLY");
343     gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0.,          idrotm[1002], "ONLY");
344     gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0.,       idrotm[1003], "ONLY");
345     gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3,  idrotm[1004], "ONLY");
346     gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3,     idrotm[1005], "ONLY");
347     gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
348     
349 }
350
351
352 //___________________________________________
353 void AliRICHv0::CreateMaterials()
354 {
355     //
356     // *** DEFINITION OF AVAILABLE RICH MATERIALS *** 
357     // ORIGIN    : NICK VAN EIJNDHOVEN 
358     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it) 
359     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it) 
360     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it) 
361     //
362     Int_t   ISXFLD = gAlice->Field()->Integ();
363     Float_t SXMGMX = gAlice->Field()->Max();
364     Int_t i;
365
366     /************************************Antonnelo's Values (14-vectors)*****************************************/
367     /*
368     Float_t ppckov[14] = { 5.63e-9,5.77e-9,5.9e-9,6.05e-9,6.2e-9,6.36e-9,6.52e-9,
369                            6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
370     Float_t rindex_quarz[14] = { 1.528309,1.533333,
371                                  1.538243,1.544223,1.550568,1.55777,
372                                  1.565463,1.574765,1.584831,1.597027,
373                                1.611858,1.6277,1.6472,1.6724 };
374     Float_t rindex_quarzo[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
375     Float_t rindex_methane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
376     Float_t rindex_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
377     Float_t absco_freon[14] = { 179.0987,179.0987,
378                                 179.0987,179.0987,179.0987,142.92,56.65,13.95,10.43,7.07,2.03,.5773,.33496,0. };
379     //Float_t absco_freon[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
380         //                       1e-5,1e-5,1e-5,1e-5,1e-5 };
381     Float_t absco_quarz[14] = { 64.035,39.98,35.665,31.262,27.527,22.815,21.04,17.52,
382                                 14.177,9.282,4.0925,1.149,.3627,.10857 };
383     Float_t absco_quarzo[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
384                                  1e-5,1e-5,1e-5,1e-5,1e-5 };
385     Float_t absco_csi[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
386                               1e-4,1e-4,1e-4,1e-4 };
387     Float_t absco_methane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
388                                   1e6,1e6,1e6 };
389     Float_t absco_gri[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
390                               1e-4,1e-4,1e-4,1e-4 };
391     Float_t effic_all[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
392     Float_t effic_csi[14] = { 6e-4,.005,.0075,.01125,.045,.117,.135,.16575,
393                               .17425,.1785,.1836,.1904,.1938,.221 };
394     Float_t effic_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
395     */
396    
397     
398     /**********************************End of Antonnelo's Values**********************************/
399     
400     /**********************************Values from rich_media.f (31-vectors)**********************************/
401     
402
403     //Photons energy intervals
404     Float_t ppckov[26];
405     for (i=0;i<26;i++) 
406     {
407         ppckov[i] = (Float_t(i)*0.1+5.5)*1e-9;
408         //printf ("Energy intervals: %e\n",ppckov[i]);
409     }
410     
411     
412     //Refraction index for quarz
413     Float_t rindex_quarz[26];
414     Float_t  e1= 10.666;
415     Float_t  e2= 18.125;
416     Float_t  f1= 46.411;
417     Float_t  f2= 228.71;
418     for (i=0;i<26;i++)
419     {
420         Float_t ene=ppckov[i]*1e9;
421         Float_t a=f1/(e1*e1 - ene*ene);
422         Float_t b=f2/(e2*e2 - ene*ene);
423         rindex_quarz[i] = TMath::Sqrt(1. + a + b );
424         //printf ("Rindex_quarz: %e\n",rindex_quarz[i]);
425     } 
426     
427     //Refraction index for opaque quarz, methane and grid
428     Float_t rindex_quarzo[26];
429     Float_t rindex_methane[26];
430     Float_t rindex_gri[26];
431     for (i=0;i<26;i++)
432     {
433         rindex_quarzo[i]=1;
434         rindex_methane[i]=1.000444;
435         rindex_gri[i]=1;
436         //printf ("Rindex_quarzo , etc: %e, %e, %e\n",rindex_quarzo[i], rindex_methane[i], rindex_gri[i]=1);
437     } 
438     
439     //Absorption index for freon
440     Float_t absco_freon[26] = {179.0987, 179.0987, 179.0987, 179.0987, 179.0987,  179.0987, 179.0987, 179.0987, 
441                                179.0987, 142.9206, 56.64957, 25.58622, 13.95293, 12.03905, 10.42953, 8.804196, 
442                                7.069031, 4.461292, 2.028366, 1.293013, .577267,   .40746,  .334964, 0., 0., 0.};
443     
444     //Absorption index for quarz
445     /*Float_t Qzt [21] = {.0,.0,.005,.04,.35,.647,.769,.808,.829,.844,.853,.858,.869,.887,.903,.902,.902,
446                         .906,.907,.907,.907};
447     Float_t Wavl2[] = {150.,155.,160.0,165.0,170.0,175.0,180.0,185.0,190.0,195.0,200.0,205.0,210.0,
448                        215.0,220.0,225.0,230.0,235.0,240.0,245.0,250.0};                                 
449     Float_t absco_quarz[31];         
450     for (Int_t i=0;i<31;i++)
451     {
452         Float_t Xlam = 1237.79 / (ppckov[i]*1e9);
453         if (Xlam <= 160) absco_quarz[i] = 0;
454         if (Xlam > 250) absco_quarz[i] = 1;
455         else 
456         {
457             for (Int_t j=0;j<21;j++)
458             {
459                 //printf ("Passed\n");
460                 if (Xlam > Wavl2[j] && Xlam < Wavl2[j+1])
461                 {
462                     Float_t Dabs = (Qzt[j+1] - Qzt[j])/(Wavl2[j+1] - Wavl2[j]);
463                     Float_t Abso = Qzt[j] + Dabs*(Xlam - Wavl2[j]);
464                     absco_quarz[i] = -5.0/(TMath::Log(Abso));
465                 } 
466             }
467         }
468         printf ("Absco_quarz: %e Absco_freon: %e for energy: %e\n",absco_quarz[i],absco_freon[i],ppckov[i]);
469     }*/
470
471     /*Float_t absco_quarz[31] = {49.64211, 48.41296, 47.46989, 46.50492, 45.13682, 44.47883, 43.1929 , 41.30922, 40.5943 ,
472                                39.82956, 38.98623, 38.6247 , 38.43448, 37.41084, 36.22575, 33.74852, 30.73901, 24.25086, 
473                                17.94531, 11.88753, 5.99128,  3.83503,  2.36661,  1.53155, 1.30582, 1.08574, .8779708, 
474                                .675275, 0., 0., 0.};
475     
476     for (Int_t i=0;i<31;i++)
477     {
478         absco_quarz[i] = absco_quarz[i]/10;
479     }*/
480
481     Float_t absco_quarz [26] = {105.8, 65.52, 48.58, 42.85, 35.79, 31.262, 28.598, 27.527, 25.007, 22.815, 21.004,
482                                 19.266, 17.525, 15.878, 14.177, 11.719, 9.282, 6.62, 4.0925, 2.601, 1.149, .667, .3627,
483                                 .192, .1497, .10857};
484     
485     //Absorption index for methane
486     Float_t absco_methane[26];
487     for (i=0;i<26;i++) 
488     {
489         absco_methane[i]=AbsoCH4(ppckov[i]*1e9); 
490         //printf("Absco_methane: %e for energy: %e\n", absco_methane[i],ppckov[i]*1e9);
491     }
492     
493     //Absorption index for opaque quarz, csi and grid, efficiency for all and grid
494     Float_t absco_quarzo[26];
495     Float_t absco_csi[26];
496     Float_t absco_gri[26];
497     Float_t effic_all[26];
498     Float_t effic_gri[26];
499     for (i=0;i<26;i++)
500     { 
501         absco_quarzo[i]=1e-5; 
502         absco_csi[i]=1e-4; 
503         absco_gri[i]=1e-4; 
504         effic_all[i]=1; 
505         effic_gri[i]=1;
506         //printf ("All must be 1: %e,  %e,  %e,  %e,  %e\n",absco_quarzo[i],absco_csi[i],absco_gri[i],effic_all[i],effic_gri[i]);
507     } 
508     
509     //Efficiency for csi 
510     
511     Float_t effic_csi[26] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983, 0.010125,
512                              0.0242999997, 0.0405000001, 0.0688500032, 0.105299994, 0.121500008, 0.141749993, 0.157949999,
513                              0.162, 0.166050002, 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992, 0.18592,
514                              0.187579989, 0.189239994, 0.190899998, 0.207499996, 0.215799987};
515         
516     
517
518     //FRESNEL LOSS CORRECTION FOR PERPENDICULAR INCIDENCE AND
519     //UNPOLARIZED PHOTONS
520
521     for (i=0;i<26;i++)
522     {
523         effic_csi[i] = effic_csi[i]/(1.-Fresnel(ppckov[i]*1e9,1.,0)); 
524         //printf ("Fresnel result: %e for energy: %e\n",Fresnel(ppckov[i]*1e9,1.,0),ppckov[i]*1e9);
525     }
526         
527     /*******************************************End of rich_media.f***************************************/
528
529   
530
531     
532     
533     
534     Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon, 
535     zmet[2], zqua[2];
536     Int_t nlmatfre;
537     Float_t densquao;
538     Int_t nlmatmet, nlmatqua;
539     Float_t wmatquao[2], rindex_freon[26];
540     Float_t aquao[2], epsil, stmin, zquao[2];
541     Int_t nlmatquao;
542     Float_t radlal, densal, tmaxfd, deemax, stemax;
543     Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
544     
545     Int_t *idtmed = fIdtmed->GetArray()-999;
546     
547     TGeant3 *geant3 = (TGeant3*) gMC;
548     
549     // --- Photon energy (GeV) 
550     // --- Refraction indexes 
551     for (i = 0; i < 26; ++i) {
552         rindex_freon[i] = ppckov[i] * .0172 * 1e9 + 1.177;
553         //printf ("Rindex_freon: %e \n Effic_csi: %e for energy: %e\n",rindex_freon[i], effic_csi[i], ppckov[i]);
554     }
555             
556     // --- Detection efficiencies (quantum efficiency for CsI) 
557     // --- Define parameters for honeycomb. 
558     //     Used carbon of equivalent rad. lenght 
559     
560     ahon    = 12.01;
561     zhon    = 6.;
562     denshon = 2.265;
563     radlhon = 18.8;
564     
565     // --- Parameters to include in GSMIXT, relative to Quarz (SiO2) 
566     
567     aqua[0]    = 28.09;
568     aqua[1]    = 16.;
569     zqua[0]    = 14.;
570     zqua[1]    = 8.;
571     densqua    = 2.64;
572     nlmatqua   = -2;
573     wmatqua[0] = 1.;
574     wmatqua[1] = 2.;
575     
576     // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2) 
577     
578     aquao[0]    = 28.09;
579     aquao[1]    = 16.;
580     zquao[0]    = 14.;
581     zquao[1]    = 8.;
582     densquao    = 2.64;
583     nlmatquao   = -2;
584     wmatquao[0] = 1.;
585     wmatquao[1] = 2.;
586     
587     // --- Parameters to include in GSMIXT, relative to Freon (C6F14) 
588     
589     afre[0]    = 12.;
590     afre[1]    = 19.;
591     zfre[0]    = 6.;
592     zfre[1]    = 9.;
593     densfre    = 1.7;
594     nlmatfre   = -2;
595     wmatfre[0] = 6.;
596     wmatfre[1] = 14.;
597     
598     // --- Parameters to include in GSMIXT, relative to methane (CH4) 
599     
600     amet[0]    = 12.01;
601     amet[1]    = 1.;
602     zmet[0]    = 6.;
603     zmet[1]    = 1.;
604     densmet    = 7.17e-4;
605     nlmatmet   = -2;
606     wmatmet[0] = 1.;
607     wmatmet[1] = 4.;
608     
609     // --- Parameters to include in GSMIXT, relative to anode grid (Cu) 
610   
611     agri    = 63.54;
612     zgri    = 29.;
613     densgri = 8.96;
614     radlgri = 1.43;
615     
616     // --- Parameters to include in GSMATE related to aluminium sheet 
617     
618     aal    = 26.98;
619     zal    = 13.;
620     densal = 2.7;
621     radlal = 8.9;
622     
623     AliMaterial(1, "Air     $", 14.61, 7.3, .001205, 30420., 67500);
624     AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
625     AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
626     AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
627     AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
628     AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
629     AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
630     AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
631     AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
632     AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
633     
634     tmaxfd = -10.;
635     stemax = -.1;
636     deemax = -.2;
637     epsil  = .001;
638     stmin  = -.001;
639     
640     AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
641     AliMedium(2, "HONEYCOMB$", 6, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
642     AliMedium(3, "QUARZO$", 20, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
643     AliMedium(4, "FREON$", 30, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
644     AliMedium(5, "METANO$", 40, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
645     AliMedium(6, "CSI$", 16, 1, ISXFLD, SXMGMX,tmaxfd, stemax, deemax, epsil, stmin);
646     AliMedium(7, "GRIGLIA$", 11, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
647     AliMedium(8, "QUARZOO$", 21, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
648     AliMedium(9, "GAP$", 41, 1, ISXFLD, SXMGMX,tmaxfd, .1, -deemax, epsil, -stmin);
649     AliMedium(10, "ALUMINUM$", 50, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
650     
651
652     geant3->Gsckov(idtmed[1000], 26, ppckov, absco_methane, effic_all, rindex_methane);
653     geant3->Gsckov(idtmed[1001], 26, ppckov, absco_methane, effic_all, rindex_methane);
654     geant3->Gsckov(idtmed[1002], 26, ppckov, absco_quarz, effic_all,rindex_quarz);
655     geant3->Gsckov(idtmed[1003], 26, ppckov, absco_freon, effic_all,rindex_freon);
656     geant3->Gsckov(idtmed[1004], 26, ppckov, absco_methane, effic_all, rindex_methane);
657     geant3->Gsckov(idtmed[1005], 26, ppckov, absco_csi, effic_csi, rindex_methane);
658     geant3->Gsckov(idtmed[1006], 26, ppckov, absco_gri, effic_gri, rindex_gri);
659     geant3->Gsckov(idtmed[1007], 26, ppckov, absco_quarzo, effic_all, rindex_quarzo);
660     geant3->Gsckov(idtmed[1008], 26, ppckov, absco_methane, effic_all, rindex_methane);
661     geant3->Gsckov(idtmed[1009], 26, ppckov, absco_gri, effic_gri, rindex_gri);
662 }
663
664 //___________________________________________
665
666 Float_t AliRICHv0::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
667 {
668
669     //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
670     
671     Float_t en[36] = {5.0,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6.0,6.1,6.2,
672                       6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,7.7,
673                       7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
674      
675
676     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
677                         2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
678                         2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
679                         1.72,1.53};
680       
681     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
682                         0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
683                         0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
684                         1.714,1.498};
685     Float_t xe=ene;
686     Int_t  j=Int_t(xe*10)-49;
687     Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
688     Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
689
690     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
691     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
692
693     Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
694     Float_t tanin=sinin/pdoti;
695
696     Float_t c1=cn*cn-ck*ck-sinin*sinin;
697     Float_t c2=4*cn*cn*ck*ck;
698     Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
699     Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
700     
701     Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
702     Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
703     
704
705     //CORRECTION FACTOR FOR SURFACE ROUGHNESS
706     //B.J. STAGG  APPLIED OPTICS, 30(1991),4113
707
708     Float_t sigraf=18.;
709     Float_t lamb=1240/ene;
710     Float_t fresn;
711  
712     Float_t  rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
713
714     if(pola)
715     {
716         Float_t pdotr=0.8;                                 //DEGREE OF POLARIZATION : 1->P , -1->S
717         fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
718     }
719     else
720         fresn=0.5*(rp+rs);
721       
722     fresn = fresn*rO;
723     return(fresn);
724 }
725
726 //__________________________________________
727
728 Float_t AliRICHv0::AbsoCH4(Float_t x)
729 {
730
731     //LOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
732     Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.};              //MB X 10^22
733     //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
734     Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
735     const Float_t losch=2.686763E19;                                      // LOSCHMIDT NUMBER IN CM-3
736     const Float_t igas1=100, igas2=0, oxy=10., wat=5., pre=750.,tem=283.;                                      
737     Float_t pn=pre/760.;
738     Float_t tn=tem/273.16;
739     
740         
741 // ------- METHANE CROSS SECTION -----------------
742 // ASTROPH. J. 214, L47 (1978)
743         
744     Float_t sm=0;
745     if (x<7.75) 
746         sm=.06e-22;
747     
748     if(x>=7.75 && x<=8.1)
749     {
750         Float_t c0=-1.655279e-1;
751         Float_t c1=6.307392e-2;
752         Float_t c2=-8.011441e-3;
753         Float_t c3=3.392126e-4;
754         sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
755     }
756     
757     if (x> 8.1)
758     {
759         Int_t j=0;
760         while (x<=em[j] && x>=em[j+1])
761         {
762             j++;
763             Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
764             sm=(sch4[j]+a*(x-em[j]))*1e-22;
765         }
766     }
767     
768     Float_t dm=(igas1/100.)*(1.-((oxy+wat)/1.e6))*losch*pn/tn;
769     Float_t abslm=1./sm/dm;
770     
771 //    ------- ISOBUTHANE CROSS SECTION --------------
772 //     i-C4H10 (ai) abs. length from curves in
773 //     Lu-McDonald paper for BARI RICH workshop .
774 //     -----------------------------------------------------------
775     
776     Float_t ai;
777     Float_t absli;
778     if (igas2 != 0) 
779     {
780         if (x<7.25)
781             ai=100000000.;
782         
783         if(x>=7.25 && x<7.375)
784             ai=24.3;
785         
786         if(x>=7.375)
787             ai=.0000000001;
788         
789         Float_t si = 1./(ai*losch*273.16/293.);                    // ISOB. CRO.SEC.IN CM2
790         Float_t di=(igas2/100.)*(1.-((oxy+wat)/1.e6))*losch*pn/tn;
791         absli =1./si/di;
792     }
793     else
794         absli=1.e18;
795 //    ---------------------------------------------------------
796 //
797 //       transmission of O2
798 //
799 //       y= path in cm, x=energy in eV
800 //       so= cross section for UV absorption in cm2
801 //       do= O2 molecular density in cm-3
802 //    ---------------------------------------------------------
803     
804     Float_t abslo;
805     Float_t so=0;
806     if(x>=6.0)
807     {
808         if(x>=6.0 && x<6.5)
809         {
810             so=3.392709e-13 * TMath::Exp(2.864104 *x);
811             so=so*1e-18;
812         }
813         
814         if(x>=6.5 && x<7.0) 
815         {
816             so=2.910039e-34 * TMath::Exp(10.3337*x);
817             so=so*1e-18;
818         }
819             
820
821         if (x>=7.0) 
822         {
823             Float_t a0=-73770.76;
824             Float_t a1=46190.69;
825             Float_t a2=-11475.44;
826             Float_t a3=1412.611;
827             Float_t a4=-86.07027;
828             Float_t a5=2.074234;
829             so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
830             so=so*1e-18;
831         }
832         
833         Float_t dox=(oxy/1e6)*losch*pn/tn;
834         abslo=1./so/dox;
835     }
836     else
837         abslo=1.e18;
838 //     ---------------------------------------------------------
839 //
840 //       transmission of H2O
841 //
842 //       y= path in cm, x=energy in eV
843 //       sw= cross section for UV absorption in cm2
844 //       dw= H2O molecular density in cm-3
845 //     ---------------------------------------------------------
846     
847     Float_t abslw;
848     
849     Float_t b0=29231.65;
850     Float_t b1=-15807.74;
851     Float_t b2=3192.926;
852     Float_t b3=-285.4809;
853     Float_t b4=9.533944;
854     
855     if(x>6.75)
856     {    
857         Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
858         sw=sw*1e-18;
859         Float_t dw=(wat/1e6)*losch*pn/tn;
860         abslw=1./sw/dw;
861     }
862     else
863         abslw=1.e18;
864             
865 //    ---------------------------------------------------------
866     
867     Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
868     return (alength);
869 }
870
871
872
873
874 //___________________________________________
875
876 void AliRICHv0::Init()
877 {
878     printf("\n\n\n Start Init for version 0 - CPC chamber type \n\n\n");
879     
880     // 
881     // Initialize Tracking Chambers
882     //
883     for (Int_t i=1; i<7; i++) {
884         //printf ("i:%d",i);
885         ( (AliRICHChamber*) (*fChambers)[i])->Init();  
886     }  
887     
888     //
889     // Set the chamber (sensitive region) GEANT identifier
890     
891     ((AliRICHChamber*)(*fChambers)[0])->SetGid(1);  
892     ((AliRICHChamber*)(*fChambers)[1])->SetGid(2);  
893     ((AliRICHChamber*)(*fChambers)[2])->SetGid(3);  
894     ((AliRICHChamber*)(*fChambers)[3])->SetGid(4);  
895     ((AliRICHChamber*)(*fChambers)[4])->SetGid(5);  
896     ((AliRICHChamber*)(*fChambers)[5])->SetGid(6);  
897     ((AliRICHChamber*)(*fChambers)[6])->SetGid(7); 
898
899     Float_t pos1[3]={0,471.8999,165.2599};
900     Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90));
901
902     Float_t pos2[3]={171,470,0};
903     Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],new TRotMatrix("rot994","rot994",90,-20,90,70,0,0));
904
905     Float_t pos3[3]={0,500,0};
906     Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],new TRotMatrix("rot995","rot995",90,0,90,90,0,0));
907     
908     Float_t pos4[3]={-171,470,0};
909     Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2], new TRotMatrix("rot996","rot996",90,20,90,110,0,0));  
910
911     Float_t pos5[3]={161.3999,443.3999,-165.3};
912     Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70));
913
914     Float_t pos6[3]={0., 471.9, -165.3,};
915     Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90));
916
917     Float_t pos7[3]={-161.399,443.3999,-165.3};
918     Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110));
919     
920     printf("\n\n\n Finished Init for version 0 - CPC chamber type\n\n\n");
921 }
922
923 //___________________________________________
924 void AliRICHv0::StepManager()
925 {
926     Int_t          copy, id;
927     static Int_t   idvol;
928     static Int_t   vol[2];
929     Int_t          ipart;
930     static Float_t hits[18];
931     static Float_t Ckov_data[19];
932     TLorentzVector Position;
933     TLorentzVector Momentum;
934     Float_t        pos[3];
935     Float_t        mom[4];
936     Float_t        Localpos[3];
937     Float_t        Localmom[4];
938     Float_t        Localtheta,Localphi;
939     Float_t        theta,phi;
940     Float_t        destep, step;
941     Float_t        ranf[2];
942     Int_t          NPads;
943     Float_t        coscerenkov;
944     static Float_t eloss, xhit, yhit, tlength;
945     const  Float_t big=1.e10;
946        
947     TClonesArray &lhits = *fHits;
948     TGeant3 *geant3 = (TGeant3*) gMC;
949     TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
950
951  //if (current->Energy()>1)
952    //{
953         
954     // Only gas gap inside chamber
955     // Tag chambers and record hits when track enters 
956     
957     idvol=-1;
958     id=gMC->CurrentVolID(copy);
959     Float_t cherenkov_loss=0;
960     //gAlice->KeepTrack(gAlice->CurrentTrack());
961     
962     gMC->TrackPosition(Position);
963     pos[0]=Position(0);
964     pos[1]=Position(1);
965     pos[2]=Position(2);
966     Ckov_data[1] = pos[0];                 // X-position for hit
967     Ckov_data[2] = pos[1];                 // Y-position for hit
968     Ckov_data[3] = pos[2];                 // Z-position for hit
969     //Ckov_data[11] = gAlice->CurrentTrack();
970
971     AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH"); 
972     
973     /********************Store production parameters for Cerenkov photons************************/ 
974 //is it a Cerenkov photon? 
975     if (gMC->TrackPid() == 50000050) {          
976
977       //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
978         //{                    
979           Float_t Ckov_energy = current->Energy();
980           //energy interval for tracking
981           if  (Ckov_energy > 5.6e-09 && Ckov_energy < 7.8e-09 )       
982             //if (Ckov_energy > 0)
983             {
984               if (gMC->IsTrackEntering()){                                     //is track entering?
985                 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
986                   {                                                          //is it in freo?
987                     if (geant3->Gctrak()->nstep<1){                          //is it the first step?
988                       Int_t mother = current->GetFirstMother(); 
989                       
990                       //printf("Second Mother:%d\n",current->GetSecondMother());
991                       
992                       Ckov_data[10] = mother;
993                       Ckov_data[11] = gAlice->CurrentTrack();
994                       Ckov_data[12] = 1;             //Media where photon was produced 1->Freon, 2->Quarz
995                       fCkov_number++;
996                       fFreon_prod=1;
997                       //printf("Index: %d\n",fCkov_number);
998                     }    //first step question
999                   }        //freo question
1000                 
1001                 if (geant3->Gctrak()->nstep<1){                                  //is it first step?
1002                   if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy))             //is it in quarz?
1003                     {
1004                       Ckov_data[12] = 2;
1005                     }    //quarz question
1006                 }        //first step question
1007                 
1008                 //printf("Before %d\n",fFreon_prod);
1009               }   //track entering question
1010               
1011               if (Ckov_data[12] == 1)                                        //was it produced in Freon?
1012                 //if (fFreon_prod == 1)
1013                 {
1014                   if (gMC->IsTrackEntering()){                                     //is track entering?
1015                     //printf("Got in");
1016                     if (gMC->VolId("META")==gMC->CurrentVolID(copy))                //is it in gap?      
1017                       {
1018                         //printf("Got in\n");
1019                         gMC->TrackMomentum(Momentum);
1020                         mom[0]=Momentum(0);
1021                         mom[1]=Momentum(1);
1022                         mom[2]=Momentum(2);
1023                         mom[3]=Momentum(3);
1024                         // Z-position for hit
1025                         
1026                         
1027                         /**************** Photons lost in second grid have to be calculated by hand************/ 
1028                         
1029                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
1030                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
1031                         gMC->Rndm(ranf, 1);
1032                         //printf("grid calculation:%f\n",t);
1033                         if (ranf[0] > t) {
1034                           geant3->StopTrack();
1035                           Ckov_data[13] = 5;
1036                           AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1037                           //printf("Lost one in grid\n");
1038                         }
1039                         /**********************************************************************************/
1040                       }    //gap
1041                     
1042                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))             //is it in csi?      
1043                       {
1044                         gMC->TrackMomentum(Momentum);
1045                         mom[0]=Momentum(0);
1046                         mom[1]=Momentum(1);
1047                         mom[2]=Momentum(2);
1048                         mom[3]=Momentum(3);
1049                         
1050                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/ 
1051                         /***********************Cerenkov phtons (always polarised)*************************/
1052                         
1053                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
1054                         Float_t t = Fresnel(Ckov_energy*1e9,cophi,1);
1055                         gMC->Rndm(ranf, 1);
1056                         if (ranf[0] < t) {
1057                           geant3->StopTrack();
1058                           Ckov_data[13] = 6;
1059                           AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1060                           //printf("Lost by Fresnel\n");
1061                         }
1062                         /**********************************************************************************/
1063                       }
1064                   } //track entering?
1065                   
1066                   
1067                   /********************Evaluation of losses************************/
1068                   /******************still in the old fashion**********************/
1069                   
1070                   Int_t i1 = geant3->Gctrak()->nmec;            //number of physics mechanisms acting on the particle
1071                   for (Int_t i = 0; i < i1; ++i) {
1072                     //        Reflection loss 
1073                     if (geant3->Gctrak()->lmec[i] == 106) {        //was it reflected
1074                       Ckov_data[13]=10;
1075                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
1076                         Ckov_data[13]=1;
1077                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
1078                         Ckov_data[13]=2;
1079                       geant3->StopTrack();
1080                       AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1081                     } //reflection question
1082                     
1083                     
1084                     //        Absorption loss 
1085                     else if (geant3->Gctrak()->lmec[i] == 101) {              //was it absorbed?
1086                       Ckov_data[13]=20;
1087                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy)) 
1088                         Ckov_data[13]=11;
1089                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR")) 
1090                         Ckov_data[13]=12;
1091                       if (gMC->CurrentVolID(copy) == gMC->VolId("META")) 
1092                         Ckov_data[13]=13;
1093                       if (gMC->CurrentVolID(copy) == gMC->VolId("GAP ")) 
1094                         Ckov_data[13]=13;
1095                       
1096                       if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC")) 
1097                         Ckov_data[13]=15;
1098                       
1099                       //        CsI inefficiency 
1100                       if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
1101                         Ckov_data[13]=16;
1102                       }
1103                       geant3->StopTrack();
1104                       AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1105                       //printf("Added cerenkov %d\n",fCkov_number);
1106                     } //absorption question 
1107                     
1108                     
1109                     //        Photon goes out of tracking scope 
1110                     else if (geant3->Gctrak()->lmec[i] == 30) {                 //is it below energy treshold?
1111                       Ckov_data[13]=21;
1112                       geant3->StopTrack();
1113                       AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1114                     }   // energy treshold question         
1115                   }  //number of mechanisms cycle
1116                   /**********************End of evaluation************************/
1117                 } //freon production question
1118             } //energy interval question
1119         //}//inside the proximity gap question
1120     } //cerenkov photon question
1121       
1122     /**************************************End of Production Parameters Storing*********************/ 
1123     
1124     
1125     /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/ 
1126     
1127     if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
1128       //printf("Cerenkov\n");
1129         if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
1130         {
1131             
1132           if (gMC->Edep() > 0.){
1133                 gMC->TrackPosition(Position);
1134                 gMC->TrackMomentum(Momentum);
1135                 pos[0]=Position(0);
1136                 pos[1]=Position(1);
1137                 pos[2]=Position(2);
1138                 mom[0]=Momentum(0);
1139                 mom[1]=Momentum(1);
1140                 mom[2]=Momentum(2);
1141                 mom[3]=Momentum(3);
1142                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
1143                 Double_t rt = TMath::Sqrt(tc);
1144                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
1145                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
1146                 gMC->Gmtod(pos,Localpos,1);                                                                    
1147                 gMC->Gmtod(mom,Localmom,2);
1148                 
1149                 gMC->CurrentVolOffID(2,copy);
1150                 vol[0]=copy;
1151                 idvol=vol[0]-1;
1152
1153                 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
1154                         //->Sector(Localpos[0], Localpos[2]);
1155                 //printf("Sector:%d\n",sector);
1156
1157                 /*if (gMC->TrackPid() == 50000051){
1158                   fFeedbacks++;
1159                   printf("Feedbacks:%d\n",fFeedbacks);
1160                 }*/     
1161                 
1162                 ((AliRICHChamber*) (*fChambers)[idvol])
1163                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
1164                 if(idvol<7) {   
1165                     Ckov_data[0] = gMC->TrackPid();        // particle type
1166                     Ckov_data[1] = pos[0];                 // X-position for hit
1167                     Ckov_data[2] = pos[1];                 // Y-position for hit
1168                     Ckov_data[3] = pos[2];                 // Z-position for hit
1169                     Ckov_data[4] = theta;                      // theta angle of incidence
1170                     Ckov_data[5] = phi;                      // phi angle of incidence 
1171                     Ckov_data[8] = (Float_t) fNPadHits;      // first padhit
1172                     Ckov_data[9] = -1;                       // last pad hit
1173                     Ckov_data[13] = 4;                       // photon was detected
1174                     Ckov_data[14] = mom[0];
1175                     Ckov_data[15] = mom[1];
1176                     Ckov_data[16] = mom[2];
1177                     
1178                     destep = gMC->Edep();
1179                     gMC->SetMaxStep(big);
1180                     cherenkov_loss  += destep;
1181                     Ckov_data[7]=cherenkov_loss;
1182                     
1183                     NPads = MakePadHits(Localpos[0],Localpos[2],cherenkov_loss,idvol,cerenkov);
1184                     if (fNPadHits > (Int_t)Ckov_data[8]) {
1185                         Ckov_data[8]= Ckov_data[8]+1;
1186                         Ckov_data[9]= (Float_t) fNPadHits;
1187                     }
1188
1189                     Ckov_data[17] = NPads;
1190                     //printf("Npads:%d",NPads);
1191                     
1192                     //TClonesArray *Hits = RICH->Hits();
1193                     AliRICHHit *mipHit =  (AliRICHHit*) (fHits->UncheckedAt(0));
1194                     if (mipHit)
1195                       {
1196                         mom[0] = current->Px();
1197                         mom[1] = current->Py();
1198                         mom[2] = current->Pz();
1199                         Float_t Mip_px = mipHit->fMomX;
1200                         Float_t Mip_py = mipHit->fMomY;
1201                         Float_t Mip_pz = mipHit->fMomZ;
1202                         
1203                         Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
1204                         Float_t rt = TMath::Sqrt(r);
1205                         Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz;  
1206                         Float_t Mip_rt = TMath::Sqrt(Mip_r);
1207                         if ((rt*Mip_rt) > 0)
1208                           {
1209                             coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt);
1210                           }
1211                         else
1212                           {
1213                             coscerenkov = 0;
1214                           }
1215                         Float_t cherenkov = TMath::ACos(coscerenkov);
1216                         Ckov_data[18]=cherenkov;
1217                       }
1218                     //if (sector != -1)
1219                     //{
1220                     AddHit(gAlice->CurrentTrack(),vol,Ckov_data);
1221                     AddCerenkov(gAlice->CurrentTrack(),vol,Ckov_data);
1222                     //}
1223                 }
1224             }
1225         }
1226     }
1227     
1228     /***********************************************End of photon hits*********************************************/
1229     
1230
1231     /**********************************************Charged particles treatment*************************************/
1232
1233     else if (gMC->TrackCharge())
1234     //else if (1 == 1)
1235       {
1236 //If MIP
1237         /*if (gMC->IsTrackEntering())
1238           {                
1239             hits[13]=20;//is track entering?
1240           }*/
1241         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1242           {
1243             fFreon_prod=1;
1244           }
1245
1246         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
1247 // Get current particle id (ipart), track position (pos)  and momentum (mom)
1248             
1249             gMC->CurrentVolOffID(3,copy);
1250             vol[0]=copy;
1251             idvol=vol[0]-1;
1252
1253             //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
1254                         //->Sector(Localpos[0], Localpos[2]);
1255             //printf("Sector:%d\n",sector);
1256             
1257             gMC->TrackPosition(Position);
1258             gMC->TrackMomentum(Momentum);
1259             pos[0]=Position(0);
1260             pos[1]=Position(1);
1261             pos[2]=Position(2);
1262             mom[0]=Momentum(0);
1263             mom[1]=Momentum(1);
1264             mom[2]=Momentum(2);
1265             mom[3]=Momentum(3);
1266             gMC->Gmtod(pos,Localpos,1);                                                                    
1267             gMC->Gmtod(mom,Localmom,2);
1268             
1269             ipart  = gMC->TrackPid();
1270             //
1271             // momentum loss and steplength in last step
1272             destep = gMC->Edep();
1273             step   = gMC->TrackStep();
1274   
1275             //
1276             // record hits when track enters ...
1277             if( gMC->IsTrackEntering()) {
1278 //              gMC->SetMaxStep(fMaxStepGas);
1279                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
1280                 Double_t rt = TMath::Sqrt(tc);
1281                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
1282                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
1283                 
1284
1285                 Double_t Localtc = Localmom[0]*Localmom[0]+Localmom[2]*Localmom[2];
1286                 Double_t Localrt = TMath::Sqrt(Localtc);
1287                 Localtheta   = Float_t(TMath::ATan2(Localrt,Double_t(Localmom[1])))*kRaddeg;                       
1288                 Localphi     = Float_t(TMath::ATan2(Double_t(Localmom[2]),Double_t(Localmom[0])))*kRaddeg;    
1289                 
1290                 hits[0] = Float_t(ipart);         // particle type
1291                 hits[1] = Localpos[0];                 // X-position for hit
1292                 hits[2] = Localpos[1];                 // Y-position for hit
1293                 hits[3] = Localpos[2];                 // Z-position for hit
1294                 hits[4] = Localtheta;                  // theta angle of incidence
1295                 hits[5] = Localphi;                    // phi angle of incidence 
1296                 hits[8] = (Float_t) fNPadHits;    // first padhit
1297                 hits[9] = -1;                     // last pad hit
1298                 hits[13] = fFreon_prod;           // did id hit the freon?
1299                 hits[14] = mom[0];
1300                 hits[15] = mom[1];
1301                 hits[16] = mom[2];
1302
1303                 tlength = 0;
1304                 eloss   = 0;
1305                 fFreon_prod = 0;
1306         
1307                 Chamber(idvol).LocaltoGlobal(Localpos,hits+1);
1308            
1309                 
1310                 //To make chamber coordinates x-y had to pass LocalPos[0], LocalPos[2]
1311                 xhit    = Localpos[0];
1312                 yhit    = Localpos[2];
1313                 // Only if not trigger chamber
1314                 if(idvol<7) {
1315                     //
1316                     //  Initialize hit position (cursor) in the segmentation model 
1317                     ((AliRICHChamber*) (*fChambers)[idvol])
1318                         ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
1319                 }
1320             }
1321             
1322             // 
1323             // Calculate the charge induced on a pad (disintegration) in case 
1324             //
1325             // Mip left chamber ...
1326             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
1327                 gMC->SetMaxStep(big);
1328                 eloss   += destep;
1329                 tlength += step;
1330                 
1331                                 
1332                 // Only if not trigger chamber
1333                 if(idvol<7) {
1334                   if (eloss > 0) 
1335                     {
1336                       if(gMC->TrackPid() == kNeutron)
1337                         printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
1338                       NPads = MakePadHits(xhit,yhit,eloss,idvol,mip);
1339                       hits[17] = NPads;
1340                       //printf("Npads:%d",NPads);
1341                     }
1342                 }
1343                 
1344                 hits[6]=tlength;
1345                 hits[7]=eloss;
1346                 if (fNPadHits > (Int_t)hits[8]) {
1347                     hits[8]= hits[8]+1;
1348                     hits[9]= (Float_t) fNPadHits;
1349                 }
1350                 
1351                 //if(sector !=-1)
1352                 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
1353                 eloss = 0; 
1354                 //
1355                 // Check additional signal generation conditions 
1356                 // defined by the segmentation
1357                 // model (boundary crossing conditions) 
1358             } else if 
1359                 (((AliRICHChamber*) (*fChambers)[idvol])
1360                  ->SigGenCond(Localpos[0], Localpos[2], Localpos[1]))
1361             {
1362                 ((AliRICHChamber*) (*fChambers)[idvol])
1363                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
1364                 if (eloss > 0) 
1365                   {
1366                     if(gMC->TrackPid() == kNeutron)
1367                       printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
1368                     NPads = MakePadHits(xhit,yhit,eloss,idvol,mip);
1369                     hits[17] = NPads;
1370                     //printf("Npads:%d",NPads);
1371                   }
1372                 xhit     = Localpos[0];
1373                 yhit     = Localpos[2]; 
1374                 eloss    = destep;
1375                 tlength += step ;
1376                 //
1377                 // nothing special  happened, add up energy loss
1378             } else {        
1379                 eloss   += destep;
1380                 tlength += step ;
1381             }
1382         }
1383       }
1384     /*************************************************End of MIP treatment**************************************/
1385    //}
1386 }
1387
1388   
1389 //___________________________________________
1390 Int_t AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, Response_t res)
1391 {
1392 //
1393 //  Calls the charge disintegration method of the current chamber and adds
1394 //  the simulated cluster to the root treee 
1395 //
1396     Int_t clhits[7];
1397     Float_t newclust[6][500];
1398     Int_t nnew;
1399     
1400 //
1401 //  Integrated pulse height on chamber
1402     
1403     clhits[0]=fNhits+1;
1404     
1405     ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
1406     Int_t ic=0;
1407     
1408 //
1409 //  Add new clusters
1410     for (Int_t i=0; i<nnew; i++) {
1411         if (Int_t(newclust[3][i]) > 0) {
1412             ic++;
1413 // Cathode plane
1414             clhits[1] = Int_t(newclust[5][i]);
1415 //  Cluster Charge
1416             clhits[2] = Int_t(newclust[0][i]);
1417 //  Pad: ix
1418             clhits[3] = Int_t(newclust[1][i]);
1419 //  Pad: iy 
1420             clhits[4] = Int_t(newclust[2][i]);
1421 //  Pad: charge
1422             clhits[5] = Int_t(newclust[3][i]);
1423 //  Pad: chamber sector
1424             clhits[6] = Int_t(newclust[4][i]);
1425             
1426             AddPadHit(clhits);
1427         }
1428     }
1429 return nnew;
1430 }