Put headers before libraries in the make.
[u/mrichter/AliRoot.git] / TRD / AliTRDv2.cxx
CommitLineData
fe4da5cc 1///////////////////////////////////////////////////////////////////////////////
2// //
3// Transition Radiation Detector version 2 -- detailed simulation //
4// //
5//Begin_Html
6/*
7<img src="gif/AliTRDv2Class.gif">
8*/
9//End_Html
10// //
11// //
12///////////////////////////////////////////////////////////////////////////////
13
14#include <TMath.h>
15#include <TRandom.h>
16#include <TVector.h>
17#include <TGeometry.h>
18#include <TNode.h>
19#include <TPGON.h>
20
21#include "GParticle.h"
22#include "AliTRDv2.h"
23#include "AliRun.h"
24#include "AliConst.h"
25#include "AliMC.h"
26
27ClassImp(AliTRDv2)
28
29//_____________________________________________________________________________
30AliTRDv2::AliTRDv2(const char *name, const char *title)
31 :AliTRD(name, title)
32{
33 //
34 // Standard constructor for Transition Radiation Detector version 2
35 //
36 fIdSenO1 = fIdSenO2 = fIdSenO3 = fIdSenO4 = fIdSenO5 = fIdSenO6 = 0;
37 fIdSenI1 = fIdSenI2 = fIdSenI3 = fIdSenI4 = fIdSenI5 = fIdSenI6 = 0;
38 SetBufferSize(128000);
39}
40
41//_____________________________________________________________________________
42void AliTRDv2::CreateGeometry()
43{
44 //
45 // Create geometry for the Transition Radiation Detector version 2
46 // This version covers the full azimuth.
47 //
48 //Begin_Html
49 /*
50 <img src="gif/AliTRDv2.gif">
51 */
52 //End_Html
53 //Begin_Html
54 /*
55 <img src="gif/AliTRDv2Tree.gif">
56 */
57 //End_Html
58
59
60 // --- Name Conventions :
61
62 // The mother volume and the support structure
63 // TRD --> Mother TRD volume (Air)
64 // UTRS --> The sectors of the detector (Air)
65 // UTSP --> The aluminum part of the support structure (Al)
66 // UTII(O) --> The inner parts of the support structure (Air)
67
68 // The chambers
69 // UCI1-6 --> The frame of the inner chambers (C)
70 // UCO1-6 --> The frame of the outer chambers (C)
71 // UII1-6 --> The inner part of the inner chambers (Air)
72 // UIO1-6 --> The inner part of the outer chambers (Air)
73
74 // The layers inside a chamber
75 // UT1I(O) --> Radiator layer (CO2)
76 // UT2I(O) --> Polyethylene layer (PE)
77 // UT3I(O) --> Mylar layer (Mylar)
78 // UXI1-6 --> Xe/C02 layer in the inner chambers (Xe/C02)
79 // UXO1-6 --> Xe/C02 layer in the outer chambers (Xe/C02)
80 // UT5I(O) --> Cu layer (pads/sensitive) (Cu)
81 // UT6I(O) --> Kapton layer (Kapton)
82 // UT7I(O) --> NOMEX layer (C)
83 // UT8I(O) --> Readout layer (Al)
84
85
86 // --- Contains geometry information
87
88 // --- Number of sectors in the full detector
89 // --- Number of modules in each sector
90 // --- z-Coordinates of the TRD-frame
91 // --- r-Coordinates of the TRD-frame
92 // --- Thickness of the aluminium of the support frame
93 // --- Thickness of the interior of the support frame
94 // --- Thickness of the carbon chamber frame
95 // --- Thickness and z-position of the PE-layer in the radiator
96 // --- Thickness and z-position of the radiator
97 // --- Thickness and z-position of the mylar-layer
98 // --- Thickness and z-position of the Xe/C02-layer
99 // --- Thickness and z-position of the Cu-layer (Pads)
100 // --- Thickness and z-position of the kapton-layer
101 // --- Thickness and z-position of the NOMEX-layer
102 // Simple C-layer for the time being
103 // --- Thickness and z-position of the readout-layer
104 // --- Parameter for the arrays
105 // --- Number of different chambers
106 // --- Number of rotation matrices
107
108 AliMC* pMC = AliMC::GetMC();
109
110 Float_t xpos, ypos, zpos;
111 Int_t icham;
112 Int_t idmat[2];
113 Float_t widma, widmi;
114 Float_t lendifc, widdifc, heightc;
115 Float_t par_ic[3], par_oc[3], par_mo[10], par_sp[4], par_ch[3];
116
117 Int_t *idtmed = gAlice->Idtmed();
118
119 //************************************************************************
120
121 // Definition of Volumes
122
123 //************************************************************************
124
125 const Int_t ncham = 6; //Number of different chambers
126
127 widmi = rmin*TMath::Sin(kPI/nsect);
128 widma = rmax*TMath::Sin(kPI/nsect);
129 // --- Some parameter for the chambers
130 lendifc = (zmax1-zmax2)/nmodul;
131 heightc = (rmax-rmin)/nmodul;
132 widdifc = (widma - widmi)/nmodul;
133 // --- Definition of the Mother volume for the TRD (Air)
134 par_mo[0] = 0.;
135 par_mo[1] = 360.;
136 par_mo[2] = nsect;
137 par_mo[3] = 2.;
138 par_mo[4] = -zmax1;
139 par_mo[5] = rmin;
140 par_mo[6] = rmax;
141 par_mo[7] = zmax1;
142 par_mo[8] = rmin;
143 par_mo[9] = rmax;
144 pMC->Gsvolu("TRD ", "PGON", idtmed[1301], par_mo, 10);
145 // --- Divide the mother volume into sectors
146 pMC->Gsdvn("UTRS", "TRD ", nsect, 2);
147 // --- Definition of the aluminum part of the support structure (Al)
148 par_sp[0] = widmi - inframe/2;
149 par_sp[1] = widma - inframe/2;
150 par_sp[2] = zmax1/2;
151 par_sp[3] = (rmax-rmin)/2;
152 pMC->Gsvolu("UTSP", "TRD1", idtmed[1300], par_sp, 4);
153 // --- Definition of the inner part of the support structure (Air)
154 par_sp[0] = widmi - inframe/2 - alfram1/2;
155 par_sp[1] = widma - inframe/2 - alfram1/2;
156 par_sp[2] = zmax1/4 -alfram2/2;
157 par_sp[3] = (rmax-rmin)/2;
158 pMC->Gsvolu("UTII", "TRD1", idtmed[1301], par_sp, 4);
159 pMC->Gsvolu("UTIO", "TRD1", idtmed[1301], par_sp, 4);
160 // --- Definition of the chambers
161 char ctagc[5], ctagi[5], ctagx[5];
162 for (icham = 1; icham <= ncham; ++icham) {
163 // --- Carbon frame of the inner chambers (C)
164 par_ch[0] = widmi + (icham-1) * widdifc - (inframe+alfram1)/2;
165 par_ch[1] = zmax1/4 -alfram2/2;
166 par_ch[2] = heightc/2.;
167 sprintf(ctagc,"UCI%1d",icham);
168 pMC->Gsvolu(ctagc, "BOX ", idtmed[1306], par_ch, 3);
169 // --- Inner part of the inner chambers (Air)
170 par_ch[0] -= ccframe;
171 par_ch[1] -= ccframe;
172 sprintf(ctagi,"UII%1d",icham);
173 pMC->Gsvolu(ctagi, "BOX ", idtmed[1301], par_ch, 3);
174 // --- Carbon frame of the outer chambers (C)
175 par_ch[0] = widmi + (icham - 1) * widdifc - 2.;
176 par_ch[1] = (icham - 6) * lendifc / 2. + zmax1/4 -alfram2/2;
177 par_ch[2] = heightc / 2.;
178 sprintf(ctagc,"UCO%1d",icham);
179 pMC->Gsvolu(ctagc, "BOX ", idtmed[1306], par_ch, 3);
180 // --- Inner part of the outer chambers (Air)
181 par_ch[0] -= ccframe;
182 par_ch[1] -= ccframe;
183 sprintf(ctagi,"UIO%1d",icham);
184 pMC->Gsvolu(ctagi, "BOX ", idtmed[1301], par_ch, 3);
185 }
186 // --- Definition of the layers in each inner chamber
187 par_ic[0] = -1.;
188 par_ic[1] = -1.;
189 // --- Radiator layer
190 par_ic[2] = rathick/2;;
191 pMC->Gsvolu("UT1I", "BOX ", idtmed[1311], par_ic, 3);
192 // --- Polyethylene layer
193 par_ic[2] = pethick/2;
194 pMC->Gsvolu("UT2I", "BOX ", idtmed[1302], par_ic, 3);
195 // --- Mylar layer
196 par_ic[2] = mythick/2;
197 pMC->Gsvolu("UT3I", "BOX ", idtmed[1307], par_ic, 3);
198 // --- Xe/CO2 layer
199 par_ic[2] = 1.8;
200 for (icham = 1; icham <= 6; ++icham) {
201 sprintf(ctagx,"UXI%1d",icham);
202 pMC->Gsvolu(ctagx, "BOX ", idtmed[1308], par_ic, 3);
203 }
204 // --- Cu layer
205 par_ic[2] = cuthick/2;
206 pMC->Gsvolu("UT5I", "BOX ", idtmed[1304], par_ic, 3);
207 // --- Kapton layer
208 par_ic[2] = kathick/2;
209 pMC->Gsvolu("UT6I", "BOX ", idtmed[1310], par_ic, 3);
210 // --- NOMEX layer
211 par_ic[2] = nothick/2;
212 pMC->Gsvolu("UT7I", "BOX ", idtmed[1309], par_ic, 3);
213 // --- Read out layer
214 par_ic[2] = rothick/2;
215 pMC->Gsvolu("UT8I", "BOX ", idtmed[1305], par_ic, 3);
216 // --- Definition of the layers in each outer chamber
217 par_oc[0] = -1.;
218 par_oc[1] = -1.;
219 // --- Radiator layer
220 par_oc[2] = rathick/2;;
221 pMC->Gsvolu("UT1O", "BOX ", idtmed[1311], par_oc, 3);
222 // --- Polyethylene layer
223 par_oc[2] = pethick/2;
224 pMC->Gsvolu("UT2O", "BOX ", idtmed[1302], par_oc, 3);
225 // --- Mylar layer
226 par_oc[2] = mythick/2;
227 pMC->Gsvolu("UT3O", "BOX ", idtmed[1307], par_oc, 3);
228 // --- Xe/CO2 layer
229 par_oc[2] = xethick/2;
230 for (icham = 1; icham <= ncham; ++icham) {
231 sprintf(ctagx,"UXO%1d",icham);
232 pMC->Gsvolu(ctagx, "BOX ", idtmed[1308], par_oc, 3);
233 }
234 // --- Cu layer
235 par_oc[2] = cuthick/2;
236 pMC->Gsvolu("UT5O", "BOX ", idtmed[1304], par_oc, 3);
237 // --- Kapton layer
238 par_oc[2] = kathick/2;
239 pMC->Gsvolu("UT6O", "BOX ", idtmed[1310], par_oc, 3);
240 // --- NOMEX layer
241 par_oc[2] = nothick/2;
242 pMC->Gsvolu("UT7O", "BOX ", idtmed[1309], par_oc, 3);
243 // --- Read out layer
244 par_oc[2] = rothick/2;
245 pMC->Gsvolu("UT8O", "BOX ", idtmed[1305], par_oc, 3);
246 //*************************************************************************
247
248 // Positioning of Volumes
249
250 //************************************************************************
251 // --- The rotation matrices
252 AliMatrix(idmat[0], 90., 90., 180., 0., 90., 0.);
253 AliMatrix(idmat[1], 90., 90., 0., 0., 90., 0.);
254 // --- Position of the layers in a chamber
255 for (icham = 1; icham <= ncham; ++icham) {
256 // --- The inner chambers
257 sprintf(ctagi,"UII%1d",icham);
258 sprintf(ctagx,"UXI%1d",icham);
259 pMC->Gspos("UT8I", icham, ctagi, 0., 0., rozpos, 0, "ONLY");
260 pMC->Gspos("UT7I", icham, ctagi, 0., 0., nozpos, 0, "ONLY");
261 pMC->Gspos("UT6I", icham, ctagi, 0., 0., kazpos, 0, "ONLY");
262 pMC->Gspos("UT5I", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
263 pMC->Gspos(ctagx, 1, ctagi, 0., 0., xezpos, 0, "ONLY");
264 pMC->Gspos("UT3I", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
265 pMC->Gspos("UT1I", icham, ctagi, 0., 0., razpos, 0, "ONLY");
266 // --- The outer chambers
267 sprintf(ctagi,"UIO%d",icham);
268 sprintf(ctagx,"UXO%d",icham);
269 pMC->Gspos("UT8O", icham, ctagi, 0., 0., rozpos, 0, "ONLY");
270 pMC->Gspos("UT7O", icham, ctagi, 0., 0., nozpos, 0, "ONLY");
271 pMC->Gspos("UT6O", icham, ctagi, 0., 0., kazpos, 0, "ONLY");
272 pMC->Gspos("UT5O", icham, ctagi, 0., 0., cuzpos, 0, "ONLY");
273 pMC->Gspos(ctagx, 1, ctagi, 0., 0., xezpos, 0, "ONLY");
274 pMC->Gspos("UT3O", icham, ctagi, 0., 0., myzpos, 0, "ONLY");
275 pMC->Gspos("UT1O", icham, ctagi, 0., 0., razpos, 0, "ONLY");
276 }
277 pMC->Gspos("UT2I", 1, "UT1I", 0., 0., 0., 0, "ONLY");
278 pMC->Gspos("UT2O", 1, "UT1O", 0., 0., 0., 0, "ONLY");
279 // --- Position of the inner part of the chambers in the carbon-frames
280 for (icham = 1; icham <= ncham; ++icham) {
281 // --- The inner chambers
282 sprintf(ctagi,"UII%1d",icham);
283 sprintf(ctagc,"UCI%1d",icham);
284 pMC->Gspos(ctagi, 1, ctagc, 0., 0., 0., 0, "ONLY");
285 // --- The outer chambers
286 sprintf(ctagi,"UIO%1d",icham);
287 sprintf(ctagc,"UCO%1d",icham);
288 pMC->Gspos(ctagi, 1, ctagc, 0., 0., 0., 0, "ONLY");
289 }
290 // --- Position of the chambers in the full TRD-setup
291 for (icham = 1; icham <= ncham; ++icham) {
292 // --- The inner chambers
293 xpos = 0.;
294 ypos = 0.;
295 zpos = (icham - .5) * heightc - (rmax-rmin)/2;
296 sprintf(ctagc,"UCI%1d",icham);
297 pMC->Gspos(ctagc, 1, "UTII", xpos, ypos, zpos, 0, "ONLY");
298 // --- The outer chambers
299 xpos = 0.;
300 ypos = 0. - (icham - 6) * lendifc / 2.;
301 zpos = (icham - .5) * heightc - (rmax-rmin)/2;
302 sprintf(ctagc,"UCO%1d",icham);
303 pMC->Gspos(ctagc, 1, "UTIO", xpos, ypos, zpos, 0, "ONLY");
304 }
305 // --- Position of the inner parts of the support frame
306 xpos = 0.;
307 ypos = zmax1/4;
308 zpos = 0.;
309 pMC->Gspos("UTIO", 1, "UTSP", xpos,-ypos, zpos, 0, "ONLY");
310 pMC->Gspos("UTII", 1, "UTSP", xpos, ypos, zpos, 0, "ONLY");
311 // --- Position of the support frame in the TRD-sectors
312 xpos = (rmax+rmin)/2;
313 ypos = 0.;
314 zpos = zmax1/2;
315 pMC->Gspos("UTSP", 1, "UTRS", xpos, ypos, zpos, idmat[0], "ONLY");
316 pMC->Gspos("UTSP", 2, "UTRS", xpos, ypos,-zpos, idmat[1], "ONLY");
317 // --- Position of TRD mother volume in ALICE experiment
318 pMC->Gspos("TRD ", 1, "ALIC", 0., 0., 0., 0, "ONLY");
319}
320
321//_____________________________________________________________________________
05e51f55 322void AliTRDv2::DrawModule()
fe4da5cc 323{
324 //
325 // Draw a shaded view of the Transition Radiation Detector version 1
326 //
327
328 AliMC* pMC = AliMC::GetMC();
329
330 // Set everything unseen
331 pMC->Gsatt("*", "seen", -1);
332 //
333 // Set ALIC mother transparent
334 pMC->Gsatt("ALIC","SEEN",0);
335 //
336 // Set the volumes visible
337 pMC->Gsatt("TRD","SEEN",0);
338 pMC->Gsatt("UTRS","SEEN",0);
339 pMC->Gsatt("UTSP","SEEN",0);
340 pMC->Gsatt("UTII","SEEN",0);
341 pMC->Gsatt("UTIO","SEEN",0);
342 pMC->Gsatt("UCI1","SEEN",0);
343 pMC->Gsatt("UII1","SEEN",0);
344 pMC->Gsatt("UCO1","SEEN",0);
345 pMC->Gsatt("UIO1","SEEN",0);
346 pMC->Gsatt("UCI2","SEEN",0);
347 pMC->Gsatt("UII2","SEEN",0);
348 pMC->Gsatt("UCO2","SEEN",0);
349 pMC->Gsatt("UIO2","SEEN",0);
350 pMC->Gsatt("UCI3","SEEN",0);
351 pMC->Gsatt("UII3","SEEN",0);
352 pMC->Gsatt("UCO3","SEEN",0);
353 pMC->Gsatt("UIO3","SEEN",0);
354 pMC->Gsatt("UCI4","SEEN",0);
355 pMC->Gsatt("UII4","SEEN",0);
356 pMC->Gsatt("UCO4","SEEN",0);
357 pMC->Gsatt("UIO4","SEEN",0);
358 pMC->Gsatt("UCI5","SEEN",0);
359 pMC->Gsatt("UII5","SEEN",0);
360 pMC->Gsatt("UCO5","SEEN",0);
361 pMC->Gsatt("UIO5","SEEN",0);
362 pMC->Gsatt("UCI6","SEEN",0);
363 pMC->Gsatt("UII6","SEEN",0);
364 pMC->Gsatt("UCO6","SEEN",0);
365 pMC->Gsatt("UIO6","SEEN",0);
366 pMC->Gsatt("UT1I","SEEN",1);
367 pMC->Gsatt("UXI1","SEEN",1);
368 pMC->Gsatt("UXI2","SEEN",1);
369 pMC->Gsatt("UXI3","SEEN",1);
370 pMC->Gsatt("UXI4","SEEN",1);
371 pMC->Gsatt("UXI5","SEEN",1);
372 pMC->Gsatt("UXI6","SEEN",1);
373 pMC->Gsatt("UT1O","SEEN",1);
374 pMC->Gsatt("UXO1","SEEN",1);
375 pMC->Gsatt("UXO2","SEEN",1);
376 pMC->Gsatt("UXO3","SEEN",1);
377 pMC->Gsatt("UXO4","SEEN",1);
378 pMC->Gsatt("UXO5","SEEN",1);
379 pMC->Gsatt("UXO6","SEEN",1);
380 //
381 pMC->Gdopt("hide", "on");
382 pMC->Gdopt("shad", "on");
383 pMC->Gsatt("*", "fill", 7);
384 pMC->SetClipBox(".");
385 pMC->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
386 pMC->DefaultRange();
387 pMC->Gdraw("alic", 40, 30, 0, 12, 9.4, .021, .021);
388 pMC->Gdhead(1111, "Transition Radiation Detector Version 2");
389 pMC->Gdman(18, 4, "MAN");
390 pMC->Gdopt("hide", "off");
391}
392
393//_____________________________________________________________________________
394void AliTRDv2::CreateMaterials()
395{
396 //
397 // Create materials for the Transition Radiation Detector version 2
398 //
399 printf("TRD: Slow simulation with fine geometry\n");
400 AliTRD::CreateMaterials();
401}
402
403//_____________________________________________________________________________
404void AliTRDv2::Init()
405{
406 //
407 // Initialise Transition Radiation Detector after geometry has been built
408 //
409 AliTRD::Init();
410 AliMC* pMC = AliMC::GetMC();
411 fIdSenI1 = pMC->VolId("UXI1");
412 fIdSenI2 = pMC->VolId("UXI2");
413 fIdSenI3 = pMC->VolId("UXI3");
414 fIdSenI4 = pMC->VolId("UXI4");
415 fIdSenI5 = pMC->VolId("UXI5");
416 fIdSenI6 = pMC->VolId("UXI6");
417
418 fIdSenO1 = pMC->VolId("UXO1");
419 fIdSenO2 = pMC->VolId("UXO2");
420 fIdSenO3 = pMC->VolId("UXO3");
421 fIdSenO4 = pMC->VolId("UXO4");
422 fIdSenO5 = pMC->VolId("UXO5");
423 fIdSenO6 = pMC->VolId("UXO6");
424}
425
426//_____________________________________________________________________________
427void AliTRDv2::StepManager()
428{
429 //
430 // Called at every step in the Transition Radiation Detector version 2
431 //
432 Int_t idSens, icSens;
433 Int_t iPla, iCha, iSec;
434 Int_t iOut;
435 Int_t vol[3];
436 Int_t iPid;
437
438 const Float_t kBig = 1.0E+12;
439
440 Float_t random[1];
441 Float_t charge;
442 Float_t betaGamma, pp;
443 Float_t aMass;
444 Float_t eLos, qTot;
445 Float_t hits[4];
446 Float_t mom[4];
447 Float_t pTot;
448
449 TClonesArray &lhits = *fHits;
450
451 AliMC* pMC = AliMC::GetMC();
452
453 // Ionization energy
454 // taken from: Ionization Measurements in High Energy Physics, Springer
455 const Double_t kWion = 23.0E-9;
456 // Maximum energy for e+ e- g for the step-size calculation
457 const Float_t kPTotMax = 0.002;
458 // Plateau value of the energy-loss for electron in xenon
459 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
460 const Float_t kPlateau = 1.70;
461 // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
462 // taken from: Ionization Measurements in High Energy Physics, Springer
463 const Float_t kPrim = 43.68;
464
465 // Set the maximum step size to a very large number for all
466 // neutral particles and those outside the driftvolume
467 pMC->SetMaxStep(kBig);
468
469 // Use only charged tracks
470 if (( pMC->TrackCharge() ) &&
471 (!pMC->TrackStop() ) &&
472 (!pMC->TrackDisappear())) {
473
474 // Find the sensitive volume
475 idSens = pMC->CurrentVol(0,icSens);
476 iPla = 0;
477 iOut = 0;
478 if (idSens == fIdSenI1) iPla = 1;
479 else if (idSens == fIdSenI2) iPla = 2;
480 else if (idSens == fIdSenI3) iPla = 3;
481 else if (idSens == fIdSenI4) iPla = 4;
482 else if (idSens == fIdSenI5) iPla = 5;
483 else if (idSens == fIdSenI6) iPla = 6;
484 else if (idSens == fIdSenO1) {
485 iPla = 1;
486 iOut = 1; }
487 else if (idSens == fIdSenO2) {
488 iPla = 2;
489 iOut = 1; }
490 else if (idSens == fIdSenO3) {
491 iPla = 3;
492 iOut = 1; }
493 else if (idSens == fIdSenO4) {
494 iPla = 4;
495 iOut = 1; }
496 else if (idSens == fIdSenO5) {
497 iPla = 5;
498 iOut = 1; }
499 else if (idSens == fIdSenO6) {
500 iPla = 6;
501 iOut = 1; }
502
503 // Inside a sensitive volume?
504 if (iPla) {
505
506 // Calculate the energy loss
507 // ((1/E^2.0) distribution for the fluctuations)
508 pMC->Rndm(random,1);
509 eLos = Eloss(random[0]);
510
511 // The amount of charge created
512 qTot = (Float_t) ((Int_t) (eLos / kWion) + 1);
513
514 // The sector number
515 pMC->CurrentVolOff(5,0,iSec);
516
517 // The chamber number
518 pMC->CurrentVolOff(4,0,iCha);
519 if (iOut)
520 iCha = (iCha - 1) * 3 + 1;
521 else
522 iCha = iCha + 1;
523
524 vol[0] = iSec;
525 vol[1] = iCha;
526 vol[2] = iPla;
527
528 pMC->TrackPosition(hits);
529 hits[3] = qTot;
530
531 // Add the hit
532 new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
533
534 pMC->TrackMomentum(mom);
535 pTot = mom[3];
536
537 // New step size for electrons only if momentum is small enough
538 iPid = pMC->TrackPid();
539 if ( (iPid > 3) ||
540 ((iPid <= 3) && (pTot < kPTotMax))) {
541 aMass = pMC->TrackMass();
542 betaGamma = pTot / aMass;
543 pp = kPrim * BetheBloch(betaGamma);
544 // Take charge > 1 into account
545 charge = pMC->TrackCharge();
546 if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
547 }
548 // Electrons above 20 Mev/c are at the plateau
549 else {
550 pp = kPrim * kPlateau;
551 }
552
553 // Calculate the maximum step size for the next tracking step
554 do
555 pMC->Rndm(random,1);
556 while ((random[0] == 1.) || (random[0] == 0.));
557 pMC->SetMaxStep( - TMath::Log(random[0]) / pp);
558
559 }
560 }
561}