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