Preliminary version of new TPC code
[u/mrichter/AliRoot.git] / TPC / AliTPCv3.cxx
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
21 ClassImp(AliTPCv3)
22  
23 //_____________________________________________________________________________
24 AliTPCv3::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 //_____________________________________________________________________________
35 void 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 //_____________________________________________________________________________
517 void 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 //_____________________________________________________________________________
565 void 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 //_____________________________________________________________________________
580 void 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 //_____________________________________________________________________________
603 void 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 //_____________________________________________________________________________
702 Float_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 }