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