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