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