Remove warnings on SunOS
[u/mrichter/AliRoot.git] / RICH / AliRICHv0.cxx
1 /////////////////////////////////////////////////////////
2 //  Manager and hits classes for set:RICH version 0    //
3 /////////////////////////////////////////////////////////
4
5 #include <TTUBE.h>
6 #include <TNode.h> 
7 #include <TRandom.h> 
8
9 #include "AliRICHv0.h"
10 #include "AliRun.h"
11 #include "AliMC.h"
12 #include "iostream.h"
13 #include "AliCallf77.h"
14 #include "AliConst.h" 
15 #include "TGeant3.h"
16
17 ClassImp(AliRICHv0)
18     
19 //___________________________________________
20 AliRICHv0::AliRICHv0() : AliRICH()
21 {
22     fChambers = 0;
23 }
24
25 //___________________________________________
26 AliRICHv0::AliRICHv0(const char *name, const char *title)
27     : AliRICH(name,title)
28 {
29     
30     fChambers = new TObjArray(7);
31     for (Int_t i=0; i<7; i++) {
32     
33         (*fChambers)[i] = new AliRICHchamber();  
34         
35     }      
36 }
37
38
39 //___________________________________________
40 void AliRICHv0::CreateGeometry()
41 {
42     //
43     // Create the geometry for RICH version 1
44     //
45     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it) 
46     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it) 
47     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it) 
48     //
49     //Begin_Html
50     /*
51       <img src="picts/AliRICHv1.gif">
52     */
53     //End_Html
54     //Begin_Html
55     /*
56       <img src="picts/AliRICHv1Tree.gif">
57     */
58     //End_Html
59     
60     
61     Int_t *idtmed = fIdtmed->GetArray()-999;
62     
63     Int_t i;
64     Float_t zs;
65     Int_t idrotm[1099];
66     Float_t par[3];
67     
68     // --- Define the RICH detector 
69     //     External aluminium box 
70     par[0] = 71.1;
71     par[1] = 11.5;
72     par[2] = 73.15;
73     gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
74     
75     //     Sensitive part of the whole RICH 
76     par[0] = 64.8;
77     par[1] = 11.5;
78     par[2] = 66.55;
79     gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
80     
81     //     Honeycomb 
82     par[0] = 63.1;
83     par[1] = .188;
84     par[2] = 66.55;
85     gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
86     
87     //     Aluminium sheet 
88     par[0] = 63.1;
89     par[1] = .025;
90     par[2] = 66.55;
91     gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
92     
93     //     Quartz 
94     par[0] = 63.1;
95     par[1] = .25;
96     par[2] = 65.5;
97     gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
98     
99     //     Spacers (cylinders) 
100     par[0] = 0.;
101     par[1] = .5;
102     par[2] = .5;
103     gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
104     
105     //     Opaque quartz 
106     par[0] = 61.95;
107     par[1] = .2;
108     par[2] = 66.5;
109     gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
110   
111     //     Frame of opaque quartz 
112     par[0] = 20.65;
113     par[1] = .5;
114     par[2] = 66.5;
115     gMC->Gsvolu("OQUF", "BOX ", idtmed[1007], par, 3);
116     
117     //     Little bar of opaque quartz 
118     par[0] = 63.1;
119     par[1] = .25;
120     par[2] = .275;
121     gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
122     
123     //     Freon 
124     par[0] = 20.15;
125     par[1] = .5;
126     par[2] = 65.5;
127     gMC->Gsvolu("FREO", "BOX ", idtmed[1003], par, 3);
128     
129     //     Methane 
130     par[0] = 64.8;
131     par[1] = 5.;
132     par[2] = 64.8;
133     gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
134     
135     //     Methane gap 
136     par[0] = 64.8;
137     par[1] = .2;
138     par[2] = 64.8;
139     gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
140     
141     //     CsI photocathode 
142     par[0] = 64.8;
143     par[1] = .25;
144     par[2] = 64.8;
145     gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
146     
147     //     Anode grid 
148     par[0] = 0.;
149     par[1] = .0025;
150     par[2] = 20.;
151     gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
152     
153     // --- Places the detectors defined with GSVOLU 
154     //     Place material inside RICH 
155     gMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
156     
157     gMC->Gspos("ALUM", 1, "SRIC", 0., -6.075, 0., 0, "ONLY");
158     gMC->Gspos("HONE", 1, "SRIC", 0., -5.862, 0., 0, "ONLY");
159     gMC->Gspos("ALUM", 2, "SRIC", 0., -5.649, 0., 0, "ONLY");
160     gMC->Gspos("OQUA", 1, "SRIC", 0., -5.424, 0., 0, "ONLY");
161     
162     AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
163     
164     for (i = 1; i <= 9; ++i) {
165         zs = (5 - i) * 14.4;
166         gMC->Gspos("SPAC", i, "FREO", 6.7, 0., zs, idrotm[1019], "ONLY");
167     }
168     for (i = 10; i <= 18; ++i) {
169         zs = (14 - i) * 14.4;
170         gMC->Gspos("SPAC", i, "FREO", -6.7, 0., zs, idrotm[1019], "ONLY");
171     }
172     
173     gMC->Gspos("FREO", 1, "OQUF", 0., 0., 0., 0, "ONLY");
174     gMC->Gspos("OQUF", 1, "SRIC", 41.3, -4.724, 0., 0, "ONLY");
175     gMC->Gspos("OQUF", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
176     gMC->Gspos("OQUF", 3, "SRIC", -41.3, -4.724, 0., 0, "ONLY");
177     gMC->Gspos("BARR", 1, "QUAR", 0., 0., -21.65, 0, "ONLY");
178     gMC->Gspos("BARR", 2, "QUAR", 0., 0., 21.65, 0, "ONLY");
179     gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
180     gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
181     gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
182     gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");
183     
184     //     Place RICH inside ALICE apparatus 
185   
186     AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
187     AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
188     AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
189     AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
190     AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
191     AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
192     AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
193     
194     gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26,     idrotm[1000], "ONLY");
195     gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0.,        idrotm[1001], "ONLY");
196     gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0.,          idrotm[1002], "ONLY");
197     gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0.,       idrotm[1003], "ONLY");
198     gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3,  idrotm[1004], "ONLY");
199     gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3,     idrotm[1005], "ONLY");
200     gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
201     
202 }
203
204
205 //___________________________________________
206 void AliRICHv0::CreateMaterials()
207 {
208     //
209     // *** DEFINITION OF AVAILABLE RICH MATERIALS *** 
210     // ORIGIN    : NICK VAN EIJNDHOVEN 
211     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it) 
212     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it) 
213     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it) 
214     //
215     Int_t   ISXFLD = gAlice->Field()->Integ();
216     Float_t SXMGMX = gAlice->Field()->Max();
217     
218     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,
219                            6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
220     Float_t rindex_quarz[14] = { 1.528309,1.533333,
221                                  1.538243,1.544223,1.550568,1.55777,
222                                  1.565463,1.574765,1.584831,1.597027,
223                                1.611858,1.6277,1.6472,1.6724 };
224     Float_t rindex_quarzo[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
225     Float_t rindex_methane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
226     Float_t rindex_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
227     Float_t absco_freon[14] = { 179.0987,179.0987,
228                                 179.0987,179.0987,179.0987,35.7,12.54,5.92,4.92,3.86,1.42,.336,.134,0. };
229     Float_t absco_quarz[14] = { 20.126,16.27,13.49,11.728,9.224,8.38,7.44,7.17,
230                                 6.324,4.483,1.6,.323,.073,0. };
231     Float_t absco_quarzo[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
232                                  1e-5,1e-5,1e-5,1e-5,1e-5 };
233     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,
234                               1e-4,1e-4,1e-4,1e-4 };
235     Float_t absco_methane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
236                                   1e6,1e6,1e6 };
237     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,
238                               1e-4,1e-4,1e-4,1e-4 };
239     Float_t effic_all[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
240     Float_t effic_csi[14] = { 4.74e-4,.00438,.009,.0182,.0282,.0653,.1141,.163,
241                               .2101,.2554,.293,.376,.3861,.418 };
242     Float_t effic_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
243     
244     Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon, 
245     zmet[2], zqua[2];
246     Int_t nlmatfre;
247     Float_t densquao;
248     Int_t nlmatmet, nlmatqua;
249     Float_t wmatquao[2], rindex_freon[14];
250     Int_t i;
251     Float_t aquao[2], epsil, stmin, zquao[2];
252     Int_t nlmatquao;
253     Float_t radlal, densal, tmaxfd, deemax, stemax;
254     Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
255     
256     Int_t *idtmed = fIdtmed->GetArray()-999;
257     
258     TGeant3 *geant3 = (TGeant3*) gMC;
259     
260     // --- Photon energy (GeV) 
261     // --- Refraction indexes 
262     for (i = 0; i < 14; ++i) {
263         rindex_freon[i] = ppckov[i] * .01095 * 1e9 + 1.2177;
264     }
265     // need to be changed 
266     
267     // --- Absorbtion lenghts (in cm) 
268     //      DATA ABSCO_QUARZ / 
269     //     &    5 * 1000000., 
270     //     &    29.85,    7.34,     4.134,    1.273,    0.722, 
271     //     &    0.365, 0.365, 0.365, 0.  / 
272     // need to be changed 
273     
274     // --- Detection efficiencies (quantum efficiency for CsI) 
275     // --- Define parameters for honeycomb. 
276     //     Used carbon of equivalent rad. lenght 
277     
278     ahon    = 12.01;
279     zhon    = 6.;
280     denshon = 2.265;
281     radlhon = 18.8;
282     
283     // --- Parameters to include in GSMIXT, relative to Quarz (SiO2) 
284     
285     aqua[0]    = 28.09;
286     aqua[1]    = 16.;
287     zqua[0]    = 14.;
288     zqua[1]    = 8.;
289     densqua    = 2.64;
290     nlmatqua   = -2;
291     wmatqua[0] = 1.;
292     wmatqua[1] = 2.;
293     
294     // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2) 
295     
296     aquao[0]    = 28.09;
297     aquao[1]    = 16.;
298     zquao[0]    = 14.;
299     zquao[1]    = 8.;
300     densquao    = 2.64;
301     nlmatquao   = -2;
302     wmatquao[0] = 1.;
303     wmatquao[1] = 2.;
304     
305     // --- Parameters to include in GSMIXT, relative to Freon (C6F14) 
306     
307     afre[0]    = 12.;
308     afre[1]    = 19.;
309     zfre[0]    = 6.;
310     zfre[1]    = 9.;
311     densfre    = 1.7;
312     nlmatfre   = -2;
313     wmatfre[0] = 6.;
314     wmatfre[1] = 14.;
315     
316     // --- Parameters to include in GSMIXT, relative to methane (CH4) 
317     
318     amet[0]    = 12.01;
319     amet[1]    = 1.;
320     zmet[0]    = 6.;
321     zmet[1]    = 1.;
322     densmet    = 7.17e-4;
323     nlmatmet   = -2;
324     wmatmet[0] = 1.;
325     wmatmet[1] = 4.;
326     
327     // --- Parameters to include in GSMIXT, relative to anode grid (Cu) 
328   
329     agri    = 63.54;
330     zgri    = 29.;
331     densgri = 8.96;
332     radlgri = 1.43;
333     
334     // --- Parameters to include in GSMATE related to aluminium sheet 
335     
336     aal    = 26.98;
337     zal    = 13.;
338     densal = 2.7;
339     radlal = 8.9;
340     
341     AliMaterial(1, "Air     $", 14.61, 7.3, .001205, 30420., 67500);
342     AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
343     AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
344     AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
345     AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
346     AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
347     AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
348     AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
349     AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
350     AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
351     
352     tmaxfd = -10.;
353     stemax = -.1;
354     deemax = -.2;
355     epsil  = .001;
356     stmin  = -.001;
357     
358     AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
359     AliMedium(2, "HONEYCOMB$", 6, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
360     AliMedium(3, "QUARZO$", 20, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
361     AliMedium(4, "FREON$", 30, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
362     AliMedium(5, "METANO$", 40, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
363     AliMedium(6, "CSI$", 16, 1, ISXFLD, SXMGMX,tmaxfd, stemax, deemax, epsil, stmin);
364     AliMedium(7, "GRIGLIA$", 11, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
365     AliMedium(8, "QUARZOO$", 21, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
366     AliMedium(9, "GAP$", 41, 1, ISXFLD, SXMGMX,tmaxfd, .1, -deemax, epsil, -stmin);
367     AliMedium(10, "ALUMINUM$", 50, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
368     
369     
370     //     Switch on delta-ray production in the methane and freon gaps 
371     
372     gMC->Gstpar(idtmed[1002], "LOSS", 1.);
373     gMC->Gstpar(idtmed[1003], "LOSS", 1.);
374     gMC->Gstpar(idtmed[1004], "LOSS", 1.);
375     gMC->Gstpar(idtmed[1008], "LOSS", 1.);
376     gMC->Gstpar(idtmed[1005], "LOSS", 1.);
377     gMC->Gstpar(idtmed[1002], "HADR", 1.);
378     gMC->Gstpar(idtmed[1003], "HADR", 1.);
379     gMC->Gstpar(idtmed[1004], "HADR", 1.);
380     gMC->Gstpar(idtmed[1008], "HADR", 1.);
381     gMC->Gstpar(idtmed[1005], "HADR", 1.);
382     gMC->Gstpar(idtmed[1002], "DCAY", 1.);
383     gMC->Gstpar(idtmed[1003], "DCAY", 1.);
384     gMC->Gstpar(idtmed[1004], "DCAY", 1.);
385     gMC->Gstpar(idtmed[1008], "DCAY", 1.);
386     gMC->Gstpar(idtmed[1005], "DCAY", 1.);
387     geant3->Gsckov(idtmed[1000], 14, ppckov, absco_methane, effic_all, rindex_methane);
388     geant3->Gsckov(idtmed[1001], 14, ppckov, absco_methane, effic_all, rindex_methane);
389     geant3->Gsckov(idtmed[1002], 14, ppckov, absco_quarz, effic_all,rindex_quarz);
390     geant3->Gsckov(idtmed[1003], 14, ppckov, absco_freon, effic_all,rindex_freon);
391     geant3->Gsckov(idtmed[1004], 14, ppckov, absco_methane, effic_all, rindex_methane);
392     geant3->Gsckov(idtmed[1005], 14, ppckov, absco_csi, effic_csi, rindex_methane);
393     geant3->Gsckov(idtmed[1006], 14, ppckov, absco_gri, effic_gri, rindex_gri);
394     geant3->Gsckov(idtmed[1007], 14, ppckov, absco_quarzo, effic_all, rindex_quarzo);
395     geant3->Gsckov(idtmed[1008], 14, ppckov, absco_methane, effic_all, rindex_methane);
396     geant3->Gsckov(idtmed[1009], 14, ppckov, absco_gri, effic_gri, rindex_gri);
397 }
398
399 //___________________________________________
400
401 void AliRICHv0::Init()
402 {
403     printf("\n\n\n Start Init for version 0 - CPC chamber type \n\n\n");
404     
405     // 
406     // Initialize Tracking Chambers
407     //
408     for (Int_t i=0; i<7; i++) {
409         ( (AliRICHchamber*) (*fChambers)[i])->Init();
410     }
411     
412     //
413     // Set the chamber (sensitive region) GEANT identifier
414     
415     ((AliRICHchamber*)(*fChambers)[0])->SetGid(1);
416     ((AliRICHchamber*)(*fChambers)[1])->SetGid(2);
417     ((AliRICHchamber*)(*fChambers)[2])->SetGid(3);
418     ((AliRICHchamber*)(*fChambers)[3])->SetGid(4);
419     ((AliRICHchamber*)(*fChambers)[4])->SetGid(5);
420     ((AliRICHchamber*)(*fChambers)[5])->SetGid(6);
421     ((AliRICHchamber*)(*fChambers)[6])->SetGid(7);
422     
423     printf("\n\n\n Finished Init for version 0 - CPC chamber type\n\n\n");
424 }
425
426 //___________________________________________
427 void AliRICHv0::StepManager()
428 {
429     Int_t          copy, id;
430     static Int_t   idvol;
431     static Int_t   vol[2];
432     Int_t          ipart;
433     static Float_t hits[9];
434     TLorentzVector Position;
435     TLorentzVector Momentum;
436     Float_t        pos[3];
437     Float_t        mom[4];
438     Float_t        Localpos[3];
439     Float_t        Localmom[4];
440     Float_t        Localtheta,Localphi;
441     Float_t        theta,phi;
442     Float_t        destep, step;
443     static Float_t eloss, xhit, yhit, tlength;
444     const  Float_t big=1.e10;
445     
446     TClonesArray &lhits = *fHits;
447     TClonesArray &lcerenkovs = *fCerenkovs;
448     
449     // Only gas gap inside chamber
450     // Tag chambers and record hits when track enters 
451     
452     idvol=-1;
453     id=gMC->CurrentVolID(copy);
454     Float_t cherenkov_loss=0.00001;
455     
456     //        Treat photons produced in Freon and Quartz 
457     if (gMC->TrackPid() == 50000050 ) {
458         if (gMC->IsTrackEntering()){
459             if (gMC->VolId("FREO")==gMC->CurrentVolID(copy) || gMC->VolId("QUAR")==gMC->CurrentVolID(copy)){
460                 //printf("GOT ONE! Type:%d    \n",gMC->TrackPid());
461             }
462         }
463     }
464     
465   
466     if (gMC->TrackPid() == 50000050 ) {
467         if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
468           {
469             if (gMC->Edep() > 0.){
470               gMC->TrackPosition(Position);
471               gMC->TrackMomentum(Momentum);
472               pos[0]=Position(0);
473               pos[1]=Position(1);
474               pos[2]=Position(2);
475               mom[0]=Momentum(0);
476               mom[1]=Momentum(1);
477               mom[2]=Momentum(2);
478               mom[3]=Momentum(3);
479               Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
480               Double_t rt = TMath::Sqrt(tc);
481               theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
482               phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
483               gMC->Gmtod(pos,Localpos,1);                                                                    
484               gMC->Gmtod(mom,Localmom,2);
485
486               gMC->CurrentVolOffID(2,copy);
487               vol[0]=copy;
488               idvol=vol[0]-1;
489
490               ((AliRICHchamber*) (*fChambers)[idvol])
491                 ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
492               if(idvol<7) {
493                 hits[0] = 50000050;               // particle type
494                 hits[1] = pos[0];                 // X-position for hit
495                 hits[2] = pos[1];                 // Y-position for hit
496                 hits[3] = pos[2];                 // Z-position for hit
497                 hits[4] = theta;                      // theta angle of incidence
498                 hits[5] = phi;                      // phi angle of incidence 
499                 hits[8] = (Float_t) fNclusters;   // first padhit
500                 hits[9] = -1;                     // last pad hit
501                 
502                 MakePadHits(Localpos[0],Localpos[2],cherenkov_loss,idvol,cerenkov);
503                 if (fNclusters > (Int_t)hits[8]) {
504                   hits[8]= hits[8]+1;
505                   hits[9]= (Float_t) fNclusters;
506                 }
507                   
508                 AddHit(gAlice->CurrentTrack(),vol,hits);
509                 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,gAlice->CurrentTrack(),vol,hits);
510                 
511               }
512             }
513           }
514         
515         
516 // 
517 // treat charged particles
518     } else if (gMC->TrackCharge())
519     {
520         
521 //If MIP
522         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
523 // Get current particle id (ipart), track position (pos)  and momentum (mom) 
524
525           gMC->CurrentVolOffID(3,copy);
526           vol[0]=copy;
527           idvol=vol[0]-1;
528                 
529             gMC->TrackPosition(Position);
530             gMC->TrackMomentum(Momentum);
531             pos[0]=Position(0);
532             pos[1]=Position(1);
533             pos[2]=Position(2);
534             mom[0]=Momentum(0);
535             mom[1]=Momentum(1);
536             mom[2]=Momentum(2);
537             mom[3]=Momentum(3);
538             gMC->Gmtod(pos,Localpos,1);                                                                    
539             gMC->Gmtod(mom,Localmom,2);
540             
541             ipart  = gMC->TrackPid();
542             //
543             // momentum loss and steplength in last step
544             destep = gMC->Edep();
545             step   = gMC->TrackStep();
546   
547             //
548             // record hits when track enters ...
549             if( gMC->IsTrackEntering()) {
550                 gMC->SetMaxStep(fMaxStepGas);
551                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
552                 Double_t rt = TMath::Sqrt(tc);
553                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
554                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
555                 
556                 Localtheta   = Float_t(TMath::ATan2(rt,Double_t(Localmom[2])))*kRaddeg;                       
557                 Localphi     = Float_t(TMath::ATan2(Double_t(Localmom[1]),Double_t(Localmom[0])))*kRaddeg;    
558                 
559                 hits[0] = Float_t(ipart);         // particle type
560                 hits[1] = pos[0];                 // X-position for hit
561                 hits[2] = pos[1];                 // Y-position for hit
562                 hits[3] = pos[2];                 // Z-position for hit
563                 hits[4] = theta;                  // theta angle of incidence
564                 hits[5] = phi;                    // phi angle of incidence 
565                 hits[8] = (Float_t) fNclusters;   // first padhit
566                 hits[9] = -1;                     // last pad hit
567                 // phi angle of incidence
568                 tlength = 0;
569                 eloss   = 0;
570                 
571                 Chamber(idvol).LocaltoGlobal(Localpos,hits+1);
572                 
573                 //To make chamber coordinates x-y had to pass LocalPos[0], LocalPos[2]
574                 xhit    = Localpos[0];
575                 yhit    = Localpos[2];
576                 // Only if not trigger chamber
577                 if(idvol<7) {
578                     //
579                     //  Initialize hit position (cursor) in the segmentation model 
580                     ((AliRICHchamber*) (*fChambers)[idvol])
581                         ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
582                 }
583             }
584             
585             // 
586             // Calculate the charge induced on a pad (disintegration) in case 
587             //
588             // Mip left chamber ...
589             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
590                 gMC->SetMaxStep(big);
591                 eloss   += destep;
592                 tlength += step;
593                 
594                 
595                 
596                 // Only if not trigger chamber
597                 if(idvol<7) {
598                     if (eloss > 0) MakePadHits(xhit,yhit,eloss,idvol,mip);
599                 }
600                 
601                 hits[6]=tlength;
602                 hits[7]=eloss;
603                 if (fNclusters > (Int_t)hits[8]) {
604                     hits[8]= hits[8]+1;
605                     hits[9]= (Float_t) fNclusters;
606                 }
607
608                 new(lhits[fNhits++]) 
609                     AliRICHhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
610                 eloss = 0; 
611                 //
612                 // Check additional signal generation conditions 
613                 // defined by the segmentation
614                 // model (boundary crossing conditions) 
615             } else if 
616                 (((AliRICHchamber*) (*fChambers)[idvol])
617                  ->SigGenCond(Localpos[0], Localpos[2], Localpos[1]))
618             {
619                 ((AliRICHchamber*) (*fChambers)[idvol])
620                     ->SigGenInit(Localpos[0], Localpos[2], Localpos[1]);
621                 if (eloss > 0) MakePadHits(xhit,yhit,eloss,idvol,mip);
622                 xhit     = Localpos[0];
623                 yhit     = Localpos[2]; 
624                 eloss    = destep;
625                 tlength += step ;
626                 //
627                 // nothing special  happened, add up energy loss
628             } else {        
629                 eloss   += destep;
630                 tlength += step ;
631             }
632         }
633     }
634 }
635
636   
637 //___________________________________________
638 void AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, Response_t res)
639 {
640 //
641 //  Calls the charge disintegration method of the current chamber and adds
642 //  the simulated cluster to the root treee 
643 //
644     Int_t clhits[7];
645     Float_t newclust[6][500];
646     Int_t nnew;
647     
648 //
649 //  Integrated pulse height on chamber
650     
651     clhits[0]=fNhits+1;
652     
653     ((AliRICHchamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
654     Int_t ic=0;
655     
656 //
657 //  Add new clusters
658     for (Int_t i=0; i<nnew; i++) {
659         if (Int_t(newclust[3][i]) > 0) {
660             ic++;
661 // Cathode plane
662             clhits[1] = Int_t(newclust[5][i]);
663 //  Cluster Charge
664             clhits[2] = Int_t(newclust[0][i]);
665 //  Pad: ix
666             clhits[3] = Int_t(newclust[1][i]);
667 //  Pad: iy 
668             clhits[4] = Int_t(newclust[2][i]);
669 //  Pad: charge
670             clhits[5] = Int_t(newclust[3][i]);
671 //  Pad: chamber sector
672             clhits[6] = Int_t(newclust[4][i]);
673             
674             AddCluster(clhits);
675         }
676     }
677 }
678
679 ClassImp(AliRICHchamber)        
680     
681     AliRICHchamber::AliRICHchamber() 
682 {
683     fSegmentation = new TObjArray(2);
684     fResponse= new TObjArray(2);
685     fnsec=1;
686     frMin=0.1;
687     frMax=140;
688 }
689
690 //  
691 //  Get reference to response model
692 AliRICHresponse* AliRICHchamber::GetResponseModel(Response_t res)
693 {
694     if (res==mip) {
695         return (AliRICHresponse*) (*fResponse)[0];
696     } else if (res==cerenkov) {
697         return (AliRICHresponse*) (*fResponse)[1];
698     }
699     return (AliRICHresponse*) (*fResponse)[0];
700 }
701
702 // Configure response model
703 void   AliRICHchamber::ResponseModel(Response_t res, AliRICHresponse* thisResponse)
704 {
705     
706     if (res==mip) {
707         (*fResponse)[0]=thisResponse;
708     } else if (res==cerenkov) {
709         (*fResponse)[1]=thisResponse;
710     }
711 }
712
713 void AliRICHchamber::Init()
714 {
715     
716     ((AliRICHsegmentation *) (*fSegmentation)[0])->Init(this);
717     if (fnsec==2) {
718         ((AliRICHsegmentation *) (*fSegmentation)[1])->Init(this);
719     }
720 }
721
722 void AliRICHchamber::LocaltoGlobal(Float_t pos[3],Float_t Localpos[3])
723 {
724
725     Double_t *fMatrix;
726     fMatrix =  fChamberMatrix->GetMatrix();
727     Localpos[0]=pos[0]*fMatrix[0]+pos[1]*fMatrix[3]+pos[2]*fMatrix[6];
728     Localpos[1]=pos[0]*fMatrix[1]+pos[1]*fMatrix[4]+pos[2]*fMatrix[7];
729     Localpos[2]=pos[0]*fMatrix[2]+pos[1]*fMatrix[5]+pos[2]*fMatrix[8];
730     Localpos[0]+=fChamberTrans[0];
731     Localpos[1]+=fChamberTrans[1];
732     Localpos[2]+=fChamberTrans[2];
733 }
734
735
736 void AliRICHchamber::DisIntegration(Float_t eloss, Float_t xhit, Float_t yhit,
737                                     Int_t& nnew,Float_t newclust[6][500],Response_t res) 
738 {
739 //    
740 //  Generates pad hits (simulated cluster) 
741 //  using the segmentation and the response model
742     
743     Float_t dx, dy;
744     Float_t local[3];
745     Float_t source[3];
746     Int_t Nfp=0;
747     //    Float_t dummy=0;
748     
749     
750     //
751     // Width of the integration area
752     //
753     dx=((AliRICHresponse*) (*fResponse)[0])->Nsigma()*((AliRICHresponse*) (*fResponse)[0])->ChwX();
754     dy=((AliRICHresponse*) (*fResponse)[0])->Nsigma()*((AliRICHresponse*) (*fResponse)[0])->ChwY();
755     //
756     // Get pulse height from energy loss
757     Float_t qtot;
758     
759     if (res==mip) {
760         qtot = ((AliRICHresponse*) (*fResponse)[0])->IntPH(eloss);
761         
762         local[0]=xhit;
763         //Z position of the wires relative to the RICH mother volume
764         local[1]=6.026;
765         local[2]=yhit;
766         //Generate feedback photons
767         Nfp  = ((AliRICHresponse*) (*fResponse)[0])->FeedBackPhotons(source,qtot);
768         //printf("\nFeedbacks (Mip)     :%d",Nfp);
769     } else if (res==cerenkov) {
770         qtot = ((AliRICHresponse*) (*fResponse)[1])->IntPH();
771         local[0]=xhit;
772         local[1]=6.026;
773         local[2]=yhit;
774         Nfp  = ((AliRICHresponse*) (*fResponse)[1])->FeedBackPhotons(source,qtot);
775         //printf("\nFeedbacks (Cerenkov):%d",Nfp);
776     }
777
778     //
779     // Loop Over Pads
780     
781     Float_t qcheck=0, qp;
782     nnew=0;
783     for (Int_t i=1; i<=fnsec; i++) {
784         qcheck=0;
785         AliRICHsegmentation * segmentation=(AliRICHsegmentation *) (*fSegmentation)[i-1];
786         for (segmentation->FirstPad(xhit, yhit, dx, dy); 
787              segmentation->MorePads(); 
788              segmentation->NextPad()) 
789         {
790             if (res==mip) {
791                 qp= ((AliRICHresponse*) (*fResponse)[0])->IntXY(segmentation);
792             }
793             if (res==cerenkov) {
794                 qp= ((AliRICHresponse*) (*fResponse)[0])->IntXY(segmentation);
795             }
796             
797             qp= TMath::Abs(qp);
798
799             if (qp > 1.e-4) {
800                 qcheck+=qp;
801                 //
802                 // --- store signal information
803                 newclust[0][nnew]=qtot;
804                 newclust[1][nnew]=segmentation->Ix();
805                 newclust[2][nnew]=segmentation->Iy();
806                 newclust[3][nnew]=qp * qtot;
807                 newclust[4][nnew]=segmentation->ISector();
808                 newclust[5][nnew]=(Float_t) i;
809                 nnew++;
810                 
811                 
812             }
813         } // Pad loop
814     } // Cathode plane loop
815 }
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834