]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCv3.cxx
Added the entry for the IRST code checking utility
[u/mrichter/AliRoot.git] / TPC / AliTPCv3.cxx
CommitLineData
4b0fdcad 1///////////////////////////////////////////////////////////////////////////////
2// //
3// Time Projection Chamber version 3 -- detailed TPC and slow simulation //
4// //
5//Begin_Html
6/*
7<img src="picts/AliTPCv3Class.gif">
8*/
9//End_Html
10// //
11// //
12///////////////////////////////////////////////////////////////////////////////
13
14#include <stdlib.h>
15#include <TMath.h>
16
17#include "AliTPCv3.h"
18#include "AliRun.h"
19#include "AliConst.h"
20
21ClassImp(AliTPCv3)
22
23//_____________________________________________________________________________
24AliTPCv3::AliTPCv3(const char *name, const char *title) :
25 AliTPC(name, title)
26{
27 //
28 // Standard constructor for Time Projection Chamber version 2
29 //
30
31 SetBufferSize(128000);
32}
33
34//_____________________________________________________________________________
35void AliTPCv3::CreateGeometry()
36{
37 //
38 // Create the geometry of Time Projection Chamber version 3
39 //
40 //Begin_Html
41 /*
42 <img src="picts/AliTPCv3.gif">
43 */
44 //End_Html
45 //Begin_Html
46 /*
47 <img src="picts/AliTPCv3Tree.gif">
48 */
49 //End_Html
50
51 Int_t *idtmed = fIdtmed->GetArray()-399;
52
53 Float_t tana, rlsl, wlsl, rssl, rlsu, wssl, wlsu,
54 rssu, wssu, alpha, x, y, sec_thick;
55
56 Float_t x1, z0, z1, x2, theta1, theta2, theta3, dm[21];
57 Int_t il, iu;
58 Float_t z_side;
59 Int_t idrotm[100];
60
61 Float_t x0l, x0u;
62 Int_t idr;
63 //Float_t thl, thu;
64 Float_t opl, opu, phi1, phi2, phi3;
65
66 // ----------------------------------------------------
67 // FIELD CAGE WITH ENDCAPS - CARBON FIBER
68 // THIS IS ALSO A TPC MOTHER VOLUME
69 // ----------------------------------------------------
70 dm[0] = 76.;
71 dm[1] = 278.;
72 dm[2] = 275.;
73
74 gMC->Gsvolu("TPC ", "TUBE", idtmed[407], dm, 3);
75 // -------------------------------------------------------
76 // drift gas Ne/CO2 (90/10 volume) - sensitive
77 // field cage thickness = 0.52% X0
78 // ----------------------------------------------------
79 dm[0] = 76.+0.09776;
80 dm[1] = 257.;
81 dm[2] = 250.;
82
83 gMC->Gsvolu("TGAS", "TUBE", idtmed[403], dm, 3);
84 // ------------------------------------------------------
85 // "side" gas volume (the same as drift gas)
86 // here the readout chambers are positioned
87 // ------------------------------------------------------
88 dm[2] = 0.5*(275.-250.);
89 z_side = dm[2];
90
91 gMC->Gsvolu("TPSG", "TUBE", idtmed[401], dm, 3);
92 // ------------------------------------------------------
93 // HV midplane - 20 microns of mylar
94 // -----------------------------------------------------
95 dm[2] = .001;
96
97 gMC->Gsvolu("TPHV", "TUBE", idtmed[405], dm, 3);
98
99 // ====================================================
100 // lower and upper readout chambers
101 // ====================================================
102 // sectros opening angles in degrees
103 // ---------------------------------------------------
104 opl = 30.;
105 opu = 15.;
106 //thl = TMath::Tan(opl * .5 * kDegrad);
107 //thu = TMath::Tan(opu * .5 * kDegrad);
108 // ---------------------------------------------------
109 // S and L-sectors radii
110 // ---------------------------------------------------
111 rssl = 88.;
112 rssu = 136.;
113 rlsl = 142.;
114 rlsu = 250.;
115 // --------------------------------------------------
116 // Sectors widths
117 // --------------------------------------------------
118 wssl = 46.5;
119 wssu = 72.2;
120 wlsl = 37.;
121 wlsu = 65.4;
122 // ---------------------------------------------------
123 // Sector thickness 25% of X0 (Al)
124 // ---------------------------------------------------
125 sec_thick = 2.225;
126 // ---------------------------------------------------
127 // S-sectors readout chambers (lower sectors)
128 // ---------------------------------------------------
129 dm[0] = wssl * .5;
130 dm[1] = wssu * .5;
131 dm[2] = sec_thick * .5;
132 dm[3] = (rssu - rssl) * .5;
133
134 x0l = rssl + dm[3];
135
136 gMC->Gsvolu("TRCS", "TRD1", idtmed[399], dm, 4);
137 // ---------------------------------------------------
138 // L-sectors readout chambers (upper sectors)
139 // ---------------------------------------------------
140 dm[0] = wlsl * .5;
141 dm[1] = wlsu * .5;
142 dm[2] = sec_thick * .5;
143 dm[3] = (rlsu - rlsl) * .5;
144
145 x0u = rlsl + dm[3];
146
147 gMC->Gsvolu("TRCL", "TRD1", idtmed[399], dm, 4);
148 // ----------------------------------------------------
149 // positioning of the S-sector readout chambers
150 // rotation matices 1-12
151 // ----------------------------------------------------
152 z1 = -z_side + sec_thick * .5;
153
154 for (il = 1; il <= 12; ++il) {
155 phi1 = (il - 1) * opl + 270.;
156 if (phi1 > 360.) {
157 phi1 += -360.;
158 }
159 theta1 = 90.;
160 phi2 = 90.;
161 theta2 = 180.;
162 phi3 = (il - 1) * opl;
163 theta3 = 90.;
164
165 idr = il;
166
167 alpha = (il - 1) * opl * kDegrad;
168 x = x0l * TMath::Cos(alpha);
169 y = x0l * TMath::Sin(alpha);
170
171 AliMatrix(idrotm[idr], theta1, phi1, theta2, phi2, theta3, phi3);
172 gMC->Gspos("TRCS", il, "TPSG", x, y, z1, idrotm[idr], "ONLY");
173 }
174 // ----------------------------------------------------
175 // positioning of the L-sector readout chambers
176 // rotation matices 13-36
177 // ----------------------------------------------------
178 for (iu = 1; iu <= 24; ++iu) {
179 phi1 = (iu - 1) * opu + 270.;
180 if (phi1 > 360.) {
181 phi1 += -360.;
182 }
183 theta1 = 90.;
184 phi2 = 90.;
185 theta2 = 180.;
186 phi3 = (iu - 1) * opu;
187 theta3 = 90.;
188
189 idr = iu + 12;
190 AliMatrix(idrotm[idr], theta1, phi1, theta2, phi2, theta3, phi3);
191
192 alpha = (iu - 1) * opu * kDegrad;
193 x = x0u * TMath::Cos(alpha);
194 y = x0u * TMath::Sin(alpha);
195
196 gMC->Gspos("TRCL", iu, "TPSG", x, y, z1, idrotm[idr], "ONLY");
197 }
198 // --------------------------------------------------------
199 // Spoke wheel structures
200 // --------------------------------------------------------
201 gMC->Gsvolu("TSWS", "TUBE", idtmed[399], dm, 0);
202
203 z0 = -z_side + 2.;
204
205 dm[0] = 82.;
206 dm[1] = 86.;
207 dm[2] = 1.;
208
209 gMC->Gsposp("TSWS", 1, "TPSG", 0, 0, z0, 0, "ONLY", dm, 3);
210
211 dm[0] = 253.;
212 dm[1] = 257.;
213
214 gMC->Gsposp("TSWS", 2, "TPSG", 0, 0, z0, 0, "ONLY", dm, 3);
215
216 dm[0] = 140.9;
217 dm[1] = 141.9;
218
219 gMC->Gsposp("TSWS", 3, "TPSG", 0, 0, z0, 0, "ONLY", dm, 3);
220
221 // -------------------------------------------------------
222 // this volumes are to avoid overlaping
223 // -------------------------------------------------------
224 z0 = 253.;
225
226 dm[0] = 76.;
227 dm[1] = 76.+0.09776;
228
229 gMC->Gsposp("TSWS", 4, "TPC ", 0, 0, z0, 0, "ONLY", dm, 3);
230 gMC->Gsposp("TSWS", 5, "TPC ", 0, 0, -z0,0, "ONLY", dm, 3);
231
232 z0 += 21.;
233
234 gMC->Gsposp("TSWS", 6, "TPC ", 0, 0, z0, 0, "ONLY", dm, 3);
235 gMC->Gsposp("TSWS", 7, "TPC ", 0, 0, -z0, 0, "ONLY", dm, 3);
236
237 dm[0] = 257.;
238 dm[1] = 257.+0.09776;
239 dm[2] = 11.5;
240
241 z0 = 263.5;
242
243 gMC->Gsposp("TSWS", 8, "TPC ", 0, 0, z0, 0, "ONLY", dm, 3);
244 gMC->Gsposp("TSWS", 9, "TPC ", 0, 0, -z0, 0, "ONLY", dm, 3);
245 // ==========================================================
246 // wheels
247 // ==========================================================
248 // ----------------------------------------------------------
249 // Large wheel -> positioned in the TPC
250 // ----------------------------------------------------------
251 dm[0] = 257.+0.09776;
252 dm[1] = 278.;
253 dm[2] = 11.5;
254 gMC->Gsvolu("TPW1", "TUBE", idtmed[399], dm, 3);
255
256 dm[0] = 259.;
257 dm[1] = 278.;
258 dm[2] = 9.5;
259
260 gMC->Gsvolu("TPW2", "TUBE", idtmed[498], dm, 3);
261
262 gMC->Gspos("TPW2", 1, "TPW1", 0, 0, 0, 0, "ONLY");
263
264 gMC->Gspos("TPW1", 1, "TPC ", 0, 0, z0, 0, "ONLY");
265 gMC->Gspos("TPW1", 2, "TPC ", 0, 0, -z0, 0, "ONLY");
266 // -----------------------------------------------------------
267 // Small wheel -> positioned in the TPSG
268 // -----------------------------------------------------------
269 dm[0] = 76.+0.09776;
270 dm[1] = 82.;
271 dm[2] = 11.5;
272
273 gMC->Gsvolu("TPW3", "TUBE", idtmed[399], dm, 3);
274
275 dm[0] = 76.+0.09776;
276 dm[1] = 80.;
277 dm[2] = 9.5;
278
279 gMC->Gsvolu("TPW4", "TUBE", idtmed[401], dm, 3);
280
281 gMC->Gspos("TPW4", 1, "TPW3", 0, 0, 0, 0, "ONLY");
282
283 z0 = 1.;
284
285 gMC->Gspos("TPW3", 1, "TPSG", 0, 0, z0, 0, "ONLY");
286 // ---------------------------------------------------------
287 // spokes, inner and outer, also the inner ring
288 // ---------------------------------------------------------
289 dm[0] = 0.5*(135.9-82.1);
290 dm[1] = 3.;
291 dm[2] = 2.;
292
293 x1 = dm[0] + 82.;
294
295 gMC->Gsvolu("TSPI", "BOX ", idtmed[399], dm, 3);
296
297 dm[1] = 2.;
298 dm[2] = 1.;
299
300 gMC->Gsvolu("TSP1", "BOX ", idtmed[498], dm, 3);
301
302 gMC->Gspos("TSP1", 1, "TSPI", 0, 0, 0, 0, "ONLY");
303
304 dm[0] = 0.5*(256.9-142.1);
305 dm[1] = 3.;
306 dm[2] = 2.;
307
308 x2 = dm[0] + 142.;
309
310 gMC->Gsvolu("TSPO", "BOX ", idtmed[399], dm, 3);
311
312 dm[1] = 2.;
313 dm[2] = 1.;
314
315 gMC->Gsvolu("TSP2", "BOX ", idtmed[498], dm, 3);
316
317 gMC->Gspos("TSP2", 1, "TSPO", 0, 0, 0, 0, "ONLY");
318 // --------------------------------------------------------
319 dm[0] = 136.;
320 dm[1] = 142.;
321 dm[2] = 2.;
322
323 gMC->Gsvolu("TSWH", "TUBE", idtmed[399], dm, 3);
324
325 dm[0] = 137.;
326 dm[1] = 141.;
327 dm[2] = 1.;
328
329 gMC->Gsvolu("TSW1", "TUBE", idtmed[498], dm, 3);
330
331 gMC->Gspos("TSW1", 1, "TSWH", 0, 0, 0, 0, "ONLY");
332
333 z0 = z_side - .16168 - 2.;
334 // --------------------------------------------------------
335 gMC->Gspos("TSWH", 1, "TPSG", 0, 0, z0, 0, "ONLY");
336 // -------------------------------------------------------
337 // posiioning of the inner spokes
338 // -------------------------------------------------------
339 for (il = 1; il <= 6; ++il) {
340 phi1 = opl * .5 + (il - 1) * 2. * opl;
341 theta1 = 90.;
342 phi2 = opl * .5 + 90. + (il - 1) * 2. * opl;
343 if (phi2 > 360.) {
344 phi2 += -360.;
345 }
346 theta2 = 90.;
347 phi3 = 0.;
348 theta3 = 0.;
349
350 alpha = phi1 * kDegrad;
351 x = x1 * TMath::Cos(alpha);
352 y = x1 * TMath::Sin(alpha);
353
354 idr = il + 36;
355
356 AliMatrix(idrotm[idr], theta1, phi1, theta2, phi2, theta3, phi3);
357 gMC->Gspos("TSPI", il, "TPSG", x, y, z0, idrotm[idr], "ONLY");
358
359 }
360
361 for (iu = 1; iu <= 12; ++iu) {
362 phi1 = opu * .5 + (iu - 1) * 2. * opu;
363 theta1 = 90.;
364 phi2 = opu * .5 + 90. + (iu - 1) * 2. * opu;
365 if (phi2 > 360.) {
366 phi2 += -360.;
367 }
368 theta2 = 90.;
369 phi3 = 0.;
370 theta3 = 0.;
371
372 alpha = phi1 * kDegrad;
373 x = x2 * TMath::Cos(alpha);
374 y = x2 * TMath::Sin(alpha);
375
376 idr = iu + 42;
377
378 AliMatrix(idrotm[idr], theta1, phi1, theta2, phi2, theta3, phi3);
379 gMC->Gspos("TSPO", iu, "TPSG", x, y, z0, idrotm[idr], "ONLY");
380 }
381 // --------------------------------------------------------
382 // endcap cover (C, 0.86% X0)
383 // --------------------------------------------------------
384 dm[0] = 76.+0.09776;
385 dm[1] = 257.;
386 dm[2] = 0.16168*0.5;
387
388 gMC->Gsvolu("TCOV", "TUBE", idtmed[407], dm, 3);
389
390 z0 = z_side - dm[2];
391
392 gMC->Gspos("TCOV", 1, "TPSG", 0, 0, z0, 0, "ONLY");
393 // --------------------------------------------------------
394 // put the readout chambers into the TPC
395 // --------------------------------------------------------
396 theta1 = 90.;
397 phi1 = 0.;
398 theta2 = 90.;
399 phi2 = 270.;
400 theta3 = 180.;
401 phi3 = 0.;
402
403 AliMatrix(idrotm[55], theta1, phi1, theta2, phi2, theta3, phi3);
404
405 z0 = z_side + 250.;
406
407 gMC->Gspos("TPSG", 1, "TPC ", 0, 0, z0, 0, "ONLY");
408 gMC->Gspos("TPSG", 2, "TPC ", 0, 0, -z0, idrotm[55], "ONLY");
409 // ---------------------------------------------------------
410 // outer gas insulation (CO2)
411 // ---------------------------------------------------------
412 dm[0] = 257.+0.09776;
413 dm[1] = 278.-0.25004;
414 dm[2] = 275.-23.;
415
416 gMC->Gsvolu("TPOI", "TUBE", idtmed[406], dm, 3);
417
418 gMC->Gspos("TPHV", 1, "TGAS", 0, 0, 0, 0, "ONLY");
419 gMC->Gspos("TGAS", 1, "TPC ", 0, 0, 0, 0, "ONLY");
420 gMC->Gspos("TPOI", 1, "TPC ", 0, 0, 0, 0, "ONLY");
421
422 gMC->Gspos("TPC ", 1, "ALIC", 0, 0, 0, 0, "ONLY");
423 // ======================================================
424 // all volumes below are positioned in ALIC
425 // ======================================================
426 // ------------------------------------------------------
427 // the last parts of the smaller wheel (TSWS)
428 // ------------------------------------------------------
429 dm[0] = 74.;
430 dm[1] = 76.;
431 dm[2] = 1.;
432
433 z0 = 253.;
434
435 gMC->Gsposp("TSWS", 10, "TPC ", 0, 0, z0, 0, "ONLY", dm, 3);
436 gMC->Gsposp("TSWS", 11, "TPC ", 0, 0, -z0, 0, "ONLY", dm, 3);
437
438 dm[0] = 70.;
439
440 z0 += 21.;
441
442 gMC->Gsposp("TSWS", 12, "TPC ", 0, 0, z0, 0, "ONLY", dm, 3);
443 gMC->Gsposp("TSWS", 13, "TPC ", 0, 0, -z0, 0, "ONLY", dm, 3);
444 // ----------------------------------------------------
445 // Inner vessel (PCON)
446 // This volume is to be positioned directly in ALIC
447 // ----------------------------------------------------
448 dm[0] = 0.;
449 dm[1] = 360.;
450 dm[2] = 4.;
451
452 dm[3] = -250.;
453 dm[4] = 75.;
454 dm[5] = 76.;
455
456 dm[6] = -64.5;
457 dm[7] = 50.;
458 dm[8] = 76.;
459
460 dm[9] = 64.5;
461 dm[10] = 50.;
462 dm[11] = 76.;
463
464 dm[12] = 250.;
465 dm[13] = 75.;
466 dm[14] = 76.;
467
468 gMC->Gsvolu("TPIV", "PCON", idtmed[407], dm, 15);
469 // --------------------------------------------------------
470 // fill the inner vessel with CO2, (HV kDegrader)
471 // cone parts have different thickness
472 // than the central barrel, according to the TP
473 // --------------------------------------------------------
474 tana = 75./185.5;
475
476 dm[0] = 0.;
477 dm[1] = 360.;
478 dm[2] = 6.;
479
480 dm[3] = -(250.-0.2162);
481 dm[4] = (185.5-0.2126)*tana+0.2126;
482 dm[5] = 76-0.001;
483
484 dm[6] = -64.5;
485 dm[7] = 50.+0.2162;
486 dm[8] = 76-0.001;
487
488 dm[9] = -64.5;
489 dm[10] = 50+0.05076;
490 dm[11] = 76-0.001;
491
492 dm[12] = 64.5;
493 dm[13] = 50+0.05076;
494 dm[14] = 76-0.001;
495
496 dm[15] = 64.5;
497 dm[16] = 50.+0.2162;
498 dm[17] = 76-0.001;
499
500 dm[18] = (250.-0.2162);
501 dm[19] = (185.5-0.2126)*tana+0.2126;
502 dm[20] = 76-0.001;
503
504 gMC->Gsvolu("TPVD", "PCON", idtmed[406], dm, 21);
505
506 gMC->Gspos("TPVD", 1, "TPIV", 0, 0, 0, 0, "ONLY");
507
508 gMC->Gspos("TPIV", 1, "ALIC", 0, 0, 0, 0, "ONLY");
509 // ---------------------------------------------------
510 // volumes ordering
511 // ---------------------------------------------------
512 gMC->Gsord("TPSG", 6);
513}
514
515
516//_____________________________________________________________________________
517void AliTPCv3::DrawDetector()
518{
519 //
520 // Draw a shaded view of the Time Projection Chamber version 1
521 //
522
523
524 // Set everything unseen
525 gMC->Gsatt("*", "seen", -1);
526 //
527 // Set ALIC mother transparent
528 gMC->Gsatt("ALIC","SEEN",0);
529 //
530 // Set the volumes visible
531 gMC->Gsatt("TPC","SEEN",0);
532 gMC->Gsatt("TGAS","SEEN",0);
533 gMC->Gsatt("TPSG","SEEN",0);
534 gMC->Gsatt("TPHV","SEEN",1);
535 gMC->Gsatt("TRCS","SEEN",1);
536 gMC->Gsatt("TRCL","SEEN",1);
537 gMC->Gsatt("TSWS","SEEN",1);
538 gMC->Gsatt("TPW1","SEEN",1);
539 gMC->Gsatt("TPW2","SEEN",1);
540 gMC->Gsatt("TPW3","SEEN",1);
541 gMC->Gsatt("TPW4","SEEN",1);
542 gMC->Gsatt("TSPI","SEEN",1);
543 gMC->Gsatt("TSP1","SEEN",0);
544 gMC->Gsatt("TSPO","SEEN",1);
545 gMC->Gsatt("TSP2","SEEN",0);
546 gMC->Gsatt("TSWH","SEEN",1);
547 gMC->Gsatt("TSW1","SEEN",1);
548 gMC->Gsatt("TPOI","SEEN",1);
549 gMC->Gsatt("TPIV","SEEN",1);
550 gMC->Gsatt("TPVD","SEEN",1);
551 //
552 gMC->Gdopt("hide", "on");
553 gMC->Gdopt("shad", "on");
554 gMC->Gsatt("*", "fill", 7);
555 gMC->SetClipBox(".");
556 gMC->SetClipBox("*", 0, 1000, -1000, 1000, -1000, 1000);
557 gMC->DefaultRange();
558 gMC->Gdraw("alic", 40, 30, 0, 12, 9.5, .025, .025);
559 gMC->Gdhead(1111, "Time Projection Chamber");
560 gMC->Gdman(18, 4, "MAN");
561 gMC->Gdopt("hide","off");
562}
563
564//_____________________________________________________________________________
565void AliTPCv3::CreateMaterials()
566{
567 //
568 // Define materials for version 2 of the Time Projection Chamber
569 //
570
571
572 //
573 // Increase maximum number of steps
574 gMC->SetMaxNStep(30000);
575 //
576 AliTPC::CreateMaterials();
577}
578
579//_____________________________________________________________________________
580void AliTPCv3::Init()
581{
582 //
583 // Initialises version 3 of the TPC after that it has been built
584 //
585 Int_t *idtmed = fIdtmed->GetArray()-399;
586
587 AliTPC::Init();
588
589 fIdSens1=gMC->VolId("TGAS"); // drift gas as a sensitive volume
590
591 gMC->SetMaxNStep(30000); // max. number of steps increased
592
593 gMC->Gstpar(idtmed[403],"LOSS",5);
594
595 printf("*** TPC version 3 initialized ***\n");
596 printf("Maximum number of steps = %d\n",gMC->GetMaxNStep());
597
598 //
599
600}
601
602//_____________________________________________________________________________
603void AliTPCv3::StepManager()
604{
605 //
606 // Called for every step in the Time Projection Chamber
607 //
608
609 //
610 // parameters used for the energy loss calculations
611 //
612 const Float_t prim = 14.35; // number of primary collisions per 1 cm
613 const Float_t poti = 20.77e-9; // first ionization potential for Ne/CO2
614 const Float_t w_ion = 35.97e-9; // energy for the ion-electron pair creation
615
616
617 const Float_t big = 1.e10;
618
619 Int_t id,copy;
620 TLorentzVector pos;
621 Float_t hits[4];
622 Int_t vol[2];
623 TClonesArray &lhits = *fHits;
624
625 vol[1]=0;
626 vol[0]=0;
627
628 //
629
630 gMC->SetMaxStep(big);
631
632 if(!gMC->IsTrackAlive()) return; // particle has disappeared
633
634 Float_t charge = gMC->TrackCharge();
635
636 if(TMath::Abs(charge)<=0.) return; // take only charged particles
637
638
639 id=gMC->CurrentVolID(copy);
640
641 // Check the sensitive volume
642
643 if (id != fIdSens1) return;
644
645 //
646 // charged particle is in the sensitive volume
647 //
648
649 if(gMC->TrackStep() > 0) {
650
651
652 Int_t nel = (Int_t)(((gMC->Edep())-poti)/w_ion) + 1;
653 nel=TMath::Min(nel,300); // 300 electrons corresponds to 10 keV
654
655 gMC->TrackPosition(pos);
656 hits[0]=pos[0];
657 hits[1]=pos[1];
658 hits[2]=pos[2];
659
660 //
661 // check the selected side of the TPC
662 //
663
664 if(fSide && fSide*hits[2]<=0.) return;
665
666 hits[3]=(Float_t)nel;
667
668 // Add this hit
669
670 new(lhits[fNhits++]) AliTPChit(fIshunt,gAlice->CurrentTrack(),vol,hits);
671
672 }
673
674 // Stemax calculation for the next step
675
676 Float_t pp;
677 TLorentzVector mom;
678 gMC->TrackMomentum(mom);
679 Float_t ptot=mom.Rho();
680 Float_t beta_gamma = ptot/gMC->TrackMass();
681
682 if(gMC->TrackPid() <= 3 && ptot > 0.002)
683 {
684 pp = prim*1.58; // electrons above 20 MeV/c are on the plateau!
685 }
686 else
687 {
688 pp=prim*BetheBloch(beta_gamma);
689 if(TMath::Abs(charge) > 1.) pp *= (charge*charge);
690 }
691
692 Float_t random[1];
693 gMC->Rndm(random,1); // good, old GRNDM from Geant3
694
695 Double_t rnd = (Double_t)random[0];
696
697 gMC->SetMaxStep(-TMath::Log(rnd)/pp);
698
699}
700
701//_____________________________________________________________________________
702Float_t AliTPCv3::BetheBloch(Float_t bg)
703{
704 //
705 // Bethe-Bloch energy loss formula
706 //
707 const Double_t p1=0.76176e-1;
708 const Double_t p2=10.632;
709 const Double_t p3=0.13279e-4;
710 const Double_t p4=1.8631;
711 const Double_t p5=1.9479;
712
713 Double_t dbg = (Double_t) bg;
714
715 Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
716
717 Double_t aa = TMath::Power(beta,p4);
718 Double_t bb = TMath::Power(1./dbg,p5);
719
720 bb=TMath::Log(p3+bb);
721
722 return ((Float_t)((p2-aa-bb)*p1/aa));
723}