New version from M.Kowalski
[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 #include "AliTPCD.h"
21 #include"AliTPCParam.h"
22
23 ClassImp(AliTPCv3)
24  
25 //_____________________________________________________________________________
26 AliTPCv3::AliTPCv3(const char *name, const char *title) :
27   AliTPC(name, title) 
28 {
29   //
30   // Standard constructor for Time Projection Chamber version 2
31   //
32
33   SetBufferSize(128000);
34 }
35  
36 //_____________________________________________________________________________
37 void AliTPCv3::CreateGeometry()
38 {
39   //
40   // Creation of the TPC coarse geometry (version 0)
41   // Origin Marek Kowalski Crakow
42   //
43   //Begin_Html
44   /*
45     <img src="picts/AliTPCv0.gif">
46   */
47   //End_Html
48   //Begin_Html
49   /*
50     <img src="picts/AliTPCv0Tree.gif">
51   */
52   //End_Html
53
54   AliTPCParam * fTPCParam = &(fDigParam->GetParam());
55
56   Int_t *idtmed = fIdtmed->GetArray();
57
58   Float_t dm[21];
59   Int_t idrotm[100];
60
61   Int_t nRotMat = 0;
62
63
64   // ---------------------------------------------------- 
65   //          FIELD CAGE WITH ENDCAPS - G10
66   //          THIS IS ALSO A TPC MOTHER VOLUME 
67   // ---------------------------------------------------- 
68
69   dm[0] = 76.;
70   dm[1] = 278.;
71   dm[2] = 275.;
72
73   gMC->Gsvolu("TPC ", "TUBE", idtmed[8], dm, 3); 
74
75   //-----------------------------------------------------
76   //  Endcap cover c-fibre 0.86% X0
77   //-----------------------------------------------------
78
79   dm[0] = 78.;
80   dm[1] = 258.;
81   dm[2] = 0.95;
82
83   gMC->Gsvolu("TPEC","TUBE",idtmed[10],dm,3);
84
85   //-----------------------------------------------------
86   // Drift gas , leave 2 cm at the outer radius
87   // and inner raddius
88   //-----------------------------------------------------
89
90   dm[0] = 78.;
91   dm[1] = 258.;
92   dm[2] = 250.;
93
94   gMC->Gsvolu("TGAS", "TUBE", idtmed[4], dm, 3);
95
96
97   //------------------------------------------------------
98   //  membrane holder - carbon fiber
99   //------------------------------------------------------
100
101
102   gMC->Gsvolu("TPMH","TUBE",idtmed[6],dm,0);
103
104   dm[0] = 252.;
105   dm[1] = 258.;
106   dm[2] = 0.2;
107
108   gMC->Gsposp("TPMH",1,"TGAS",0.,0.,0.,0,"ONLY",dm,3);
109  
110   dm[0] = 78.;
111   dm[2] = 82.;
112   dm[2] = 0.1;
113
114   gMC->Gsposp("TPMH",2,"TGAS",0.,0.,0.,0,"ONLY",dm,3);
115
116   //----------------------------------------------------------
117   //  HV membrane - 25 microns of mylar
118   //----------------------------------------------------------
119
120   dm[0] = 82.;
121   dm[1] = 252.;
122   dm[2] = 0.00125;
123
124   gMC->Gsvolu("TPHV","TUBE",idtmed[5],dm,3);
125
126   gMC->Gspos("TPHV",1,"TGAS",0.,0.,0.,0,"ONLY");
127
128   gMC->Gspos("TGAS",1,"TPC ",0.,0.,0.,0,"ONLY");
129
130   //----------------------------------------------------------
131   // "side" gas volume, the same as the drift gas
132   // the readout chambers are placed there.  
133   //----------------------------------------------------------
134
135   dm[0] = 78.;
136   dm[1] = 258.;
137   dm[2] = 0.5*(275. - 250.);
138    
139   gMC->Gsvolu("TPSG", "TUBE", idtmed[2], dm, 3);
140
141   Float_t z_side = dm[2]; // 1/2 of the side gas thickness
142
143   //-----------------------------------------------------------
144   //   Readout chambers , 25% of X0, I use Al as the material
145   //-----------------------------------------------------------
146
147   Float_t InnerOpenAngle = fTPCParam->GetInnerAngle();
148   Float_t OuterOpenAngle = fTPCParam->GetOuterAngle();
149
150   Float_t InnerAngleShift = fTPCParam->GetInnerAngleShift();
151   Float_t OuterAngleShift = fTPCParam->GetOuterAngleShift();
152
153
154   Int_t nInnerSector = fTPCParam->GetNInnerSector()/2;
155   Int_t nOuterSector = fTPCParam->GetNOuterSector()/2;
156
157
158   Float_t InSecLowEdge = fTPCParam->GetInSecLowEdge();
159   Float_t InSecUpEdge =  fTPCParam->GetInSecUpEdge();
160
161   Float_t OuSecLowEdge = fTPCParam->GetOuSecLowEdge();
162   Float_t OuSecUpEdge = fTPCParam->GetOuSecUpEdge();
163
164   Float_t SecThick = 2.225; // Al
165
166   Float_t edge = fTPCParam->GetEdge();
167
168   //  S (Inner) sectors
169
170   dm[0] = InSecLowEdge*TMath::Tan(0.5*InnerOpenAngle)-edge;
171   dm[1] = InSecUpEdge*TMath::Tan(0.5*InnerOpenAngle)-edge;
172   dm[2] = SecThick;
173   dm[3] = 0.5*(InSecUpEdge-InSecLowEdge);
174
175   Float_t xCenterS = InSecLowEdge+dm[3];
176
177   gMC->Gsvolu("TRCS", "TRD1", idtmed[0], dm, 4); 
178
179   //  L (Outer) sectors
180
181   dm[0] = OuSecLowEdge*TMath::Tan(0.5*OuterOpenAngle)-edge;
182   dm[1] = OuSecUpEdge*TMath::Tan(0.5*OuterOpenAngle)-edge;
183   dm[2] = SecThick;
184   dm[3] = 0.5*(OuSecUpEdge-OuSecLowEdge);
185
186   Float_t xCenterL = OuSecLowEdge+dm[3];  
187
188   gMC->Gsvolu("TRCL", "TRD1", idtmed[0], dm, 4);
189
190   Float_t z1 = -z_side + SecThick*0.5;
191
192   //------------------------------------------------------------------
193   // Positioning of the S-sector readout chambers
194   //------------------------------------------------------------------
195
196   Int_t ns;
197   Float_t theta1,theta2,theta3;
198   Float_t phi1,phi2,phi3;
199   Float_t alpha;
200   Float_t x,y;
201
202   for(ns=0;ns<nInnerSector;ns++){
203     
204     phi1 = ns * InnerOpenAngle + 270.*kDegrad + InnerAngleShift;
205     phi1 *= kRaddeg; // in degrees
206
207     phi1 = (Float_t)TMath::Nint(phi1);
208
209     if (phi1 > 360.) phi1 -= 360.;
210       
211     theta1 = 90.;
212     phi2   = 90.;
213     theta2 = 180.;
214     phi3   = ns * InnerOpenAngle + InnerAngleShift;
215     phi3 *= kRaddeg; // in degrees
216
217     phi3 = (Float_t)TMath::Nint(phi3);
218       
219     if(phi3 > 360.) phi3 -= 360.;
220
221     theta3 = 90.;
222
223     alpha = phi3*kDegrad;
224
225     x = xCenterS * TMath::Cos(alpha);
226     y = xCenterS * TMath::Sin(alpha); 
227  
228     AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3);  
229      
230     gMC->Gspos("TRCS", ns+1, "TPSG", x, y, z1, idrotm[nRotMat], "ONLY");
231
232     nRotMat++;     
233
234   }
235     
236   //-------------------------------------------------------------------
237   //  Positioning of the L-sectors readout chambers
238   //-------------------------------------------------------------------
239     
240   for(ns=0;ns<nOuterSector;ns++){
241     phi1 = ns * OuterOpenAngle + 270.*kDegrad + OuterAngleShift;
242     phi1 *= kRaddeg; // in degrees
243
244     phi1 = (Float_t)TMath::Nint(phi1);
245     
246
247     if (phi1 > 360.) phi1 -= 360.;
248       
249     theta1 = 90.;
250     phi2   = 90.;
251     theta2 = 180.;
252     phi3   = ns * OuterOpenAngle+OuterAngleShift;
253     phi3 *= kRaddeg; // in degrees
254
255     phi3 = (Float_t)TMath::Nint(phi3);
256
257       
258     if(phi3 > 360.) phi3 -= 360.;
259
260     theta3 = 90.;
261
262     alpha = phi3*kDegrad;
263
264     x = xCenterL * TMath::Cos(alpha);
265     y = xCenterL * TMath::Sin(alpha); 
266  
267     AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3);  
268      
269
270     gMC->Gspos("TRCL", ns+1, "TPSG", x, y, z1, idrotm[nRotMat], "ONLY"); 
271
272     nRotMat++;   
273
274   }
275
276   Float_t z0 = z_side - 0.95;
277
278   gMC->Gspos("TPEC",1,"TPSG",0.,0.,z0,0,"ONLY");
279
280   // ========================================================== 
281   //                  wheels 
282   // ========================================================== 
283
284   //
285   //  auxilary structures
286   //
287
288
289   gMC->Gsvolu("TPWI","TUBE",idtmed[24],dm,0); // "air" 
290
291   // ---------------------------------------------------------- 
292   //       Large wheel -> positioned in the TPC 
293   // ---------------------------------------------------------- 
294   
295
296   z0 = 263.5; // TPC length - 1/2 spoke wheel width
297
298   dm[0] = 258.;
299   dm[1] = 278.;
300   dm[2] = 11.5;
301   
302   gMC->Gsvolu("TPWL", "TUBE", idtmed[0], dm, 3); 
303
304   dm[0] = dm[0]+2.;
305   dm[1] = 278.;
306   dm[2] = dm[2]-2.;
307
308   gMC->Gsposp("TPWI",1,"TPWL",0.,0.,0.,0,"ONLY",dm,3);
309
310   gMC->Gspos("TPWL", 1, "TPC ", 0, 0, z0, 0, "ONLY");
311   gMC->Gspos("TPWL", 2, "TPC ", 0, 0, -z0, 0, "ONLY");
312
313   //
314   //  Outer vessel + CO2 HV degrader
315   //
316
317   dm[0] = 260.;
318   dm[1] = 278.;
319   dm[2] = 252.;
320
321   gMC->Gsvolu("TPCO","TUBE",idtmed[12],dm,3);
322
323   dm[0] = 275.;
324   dm[1] = 278.;
325   
326   gMC->Gsvolu("TPOV","TUBE",idtmed[10],dm,3);
327
328   gMC->Gspos("TPOV",1,"TPCO",0.,0.,0.,0,"ONLY");
329
330
331   // G10 plugs
332
333   dm[0] = 258.;
334   dm[1] = 260.;
335   dm[2] = 1.;
336
337   gMC->Gsvolu("TPG1","TUBE",idtmed[8],dm,3);
338   gMC->Gspos("TPG1",1,"TPCO",0.,0.,251.,0,"ONLY");
339   gMC->Gspos("TPG1",2,"TPCO",0.,0.,-251.,0,"ONLY");  
340
341   gMC->Gspos("TPCO",1,"TPC ",0.,0.,0.,0,"ONLY");
342
343
344   //----------------------------------------------------------
345   //  Small wheel -> positioned in "side gas
346   //----------------------------------------------------------
347
348   dm[0] = 78.;
349   dm[1] = 82.;
350   dm[2] = 11.5;
351
352   gMC->Gsvolu("TPWS", "TUBE", idtmed[0], dm, 3);
353
354   dm[0] = 78.;
355   dm[1] = dm[1]-2;
356   dm[2] = dm[2]-2.;
357
358   gMC->Gsvolu("TPW1", "TUBE", idtmed[2], dm, 3);
359   
360   gMC->Gspos("TPW1", 1, "TPWS", 0., 0., 0., 0, "ONLY");
361
362   z0 = 1.; // spoke wheel is shifted w.r.t. center of the "side gas"
363
364   gMC->Gspos("TPWS", 1, "TPSG", 0, 0, z0, 0, "ONLY");
365
366
367   // to avoid overlaps
368
369   dm[0] = 76.;
370   dm[1] = 78.;
371   dm[2] = 11.5;
372
373   gMC->Gsvolu("TPS1","TUBE",idtmed[0],dm,3);
374
375   dm[2] = 9.5;
376
377   gMC->Gsvolu("TPS2","TUBE",idtmed[24],dm,3);
378
379   gMC->Gspos("TPS2",1,"TPS1",0.,0.,0.,0,"ONLY");
380
381   z0= 263.5;
382   
383   gMC->Gspos("TPS1",1,"TPC ",0.,0.,z0,0,"ONLY");
384   gMC->Gspos("TPS1",2,"TPC ",0.,0.,-z0,0,"ONLY");
385
386   // G10 plug
387
388   dm[0] = 76.;
389   dm[2] = 78.;
390   dm[3] = 1.;
391
392   gMC->Gsvolu("TPG2","TUBE",idtmed[8],dm,3);
393
394   z0 = 251.;
395
396   gMC->Gspos("TPG2",1,"TPC ",0.,0.,z0,0,"ONLY");
397   gMC->Gspos("TPG2",2,"TPC ",0.,0.,-z0,0,"ONLY");
398
399
400   //---------------------------------------------------------
401   //  central wheel  6 (radial direction) x 4 (along z) cm2
402   //---------------------------------------------------------
403
404   dm[0] = 140.;
405   dm[1] = 146.;
406   dm[2] = 2.;
407
408   gMC->Gsvolu("TPWC","TUBE",idtmed[0],dm,3);
409
410   dm[0] = dm[0] + 2.;
411   dm[1] = dm[1] - 2.;
412   dm[2] = dm[2] - 1.;
413
414   gMC->Gsposp("TPWI",2,"TPWC",0.,0.,0.,0,"ONLY",dm,3);
415
416   z0 = z_side - 1.9 - 2.;
417
418   gMC->Gspos("TPWC",1,"TPSG",0.,0.,z0,0,"ONLY");
419
420   //
421
422   gMC->Gsvolu("TPSE","BOX ",idtmed[24],dm,0); // "empty" part of the spoke 
423
424  
425   //---------------------------------------------------------
426   //  inner spokes (nSectorInner)
427   //---------------------------------------------------------
428
429   dm[0] = 0.5*(139.9-82.1);
430   dm[1] = 3.;
431   dm[2] = 2.;
432
433   Float_t x1 = dm[0]+82.;
434
435   gMC->Gsvolu("TPSI","BOX",idtmed[0],dm,3);
436
437   dm[1] = dm[1]-1.;
438   dm[2] = dm[2]-1.;
439
440   gMC->Gsposp("TPSE",1,"TPSI",0.,0.,0.,0,"ONLY",dm,3);
441
442   for(ns=0;ns<nInnerSector;ns++){
443
444     phi1 = 0.5*InnerOpenAngle + ns*InnerOpenAngle + InnerAngleShift;
445     theta1=90.;
446     phi1 *=kRaddeg;
447
448     phi1 = (Float_t)TMath::Nint(phi1);
449
450     phi2 = phi1+90.;
451     if(phi2>360.) phi2 -= 360.;
452     theta2=90.;
453     phi3=0.;
454     theta3=0.;
455
456     alpha = phi1 * kDegrad;
457     x     = x1 * TMath::Cos(alpha);
458     y     = x1 * TMath::Sin(alpha);    
459
460    AliMatrix(idrotm[nRotMat],theta1,phi1,theta2,phi2,theta3,phi3);
461
462    gMC->Gspos("TPSI",ns+1,"TPSG",x,y,z0,idrotm[nRotMat],"ONLY");  
463
464    nRotMat++;
465
466   }
467
468   //-------------------------------------------------------------
469   // outer spokes (nSectorOuter)
470   //-------------------------------------------------------------
471
472   dm[0] = 0.5*(257.9-146.1);
473   dm[1] = 3.;
474   dm[2] = 2.;
475
476   x1 = dm[0] + 146.;
477
478   gMC->Gsvolu("TPSO","BOX ",idtmed[0],dm,3);
479
480   dm[1] = dm[1] - 1.;
481   dm[2] = dm[2] - 1.;
482
483   gMC->Gsposp("TPSE",2,"TPSO",0.,0.,0.,0,"ONLY",dm,3);
484
485   for(ns=0;ns<nOuterSector;ns++){
486
487     phi1 = 0.5*OuterOpenAngle + ns*OuterOpenAngle + OuterAngleShift;
488     theta1=90.;
489     phi1 *=kRaddeg;
490
491     phi1 = (Float_t)TMath::Nint(phi1);
492
493     phi2 = phi1+90.;
494     if(phi2>360.) phi2 -= 360.;
495     theta2=90.;
496     phi3=0.;
497     theta3=0.;
498
499     alpha = phi1 * kDegrad;
500     x     = x1 * TMath::Cos(alpha);
501     y     = x1 * TMath::Sin(alpha);    
502
503    AliMatrix(idrotm[nRotMat],theta1,phi1,theta2,phi2,theta3,phi3);
504
505    gMC->Gspos("TPSO",ns+1,"TPSG",x,y,z0,idrotm[nRotMat],"ONLY");  
506
507    nRotMat++;
508
509   }  
510   
511
512   
513   // -------------------------------------------------------- 
514   //         put the readout chambers into the TPC 
515   // -------------------------------------------------------- 
516
517   theta1 = 90.;
518   phi1   = 0.;
519   theta2 = 90.;
520   phi2   = 270.;
521   theta3 = 180.;
522   phi3   = 0.;
523   
524   AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3);
525   
526   z0 = z_side + 250.;
527   
528   gMC->Gspos("TPSG", 1, "TPC ", 0, 0, z0, 0, "ONLY");
529   gMC->Gspos("TPSG", 2, "TPC ", 0, 0, -z0, idrotm[nRotMat], "ONLY");
530   
531   gMC->Gspos("TPC ", 1, "ALIC", 0, 0, 0, 0, "ONLY");
532
533   //----------------------------------------------------
534   //  Inner vessel and HV degrader
535   //----------------------------------------------------
536
537   dm[0] = 0.;
538   dm[1] = 360.;
539   dm[2] = 4.;
540   
541   dm[3] = -250.;
542   dm[4] = 74.4;
543   dm[5] = 76.;
544
545   dm[6] = -64.5;
546   dm[7] = 50.;
547   dm[8] = 76.;
548
549   dm[9] = -64.5;
550   dm[10] = 50.;
551   dm[11] = 76.;
552
553   dm[12] = 250.;
554   dm[13] = 74.4;
555   dm[14] = 76.;
556
557   gMC->Gsvolu("TPVD", "PCON", idtmed[12], dm, 15); // CO2
558
559   // cone parts
560
561   dm[0] = 0.;
562   dm[1] = 360.;
563   dm[2] = 2.;
564
565   dm[3] = 64.5;
566   dm[4] = 50.;
567   dm[5] = 51.6;
568  
569   dm[6] = 250.;
570   dm[7] = 74.4;
571   dm[8] = 76.;
572
573
574   gMC->Gsvolu("TIVC","PCON",idtmed[11],dm,9); // C-fibre
575
576   gMC->Gspos("TIVC",1,"TPVD",0.,0.,0.,0,"ONLY");
577   gMC->Gspos("TIVC",2,"TPVD",0.,0.,0.,idrotm[nRotMat],"ONLY");
578
579   // barrel part
580
581   dm[0] = 50.;
582   dm[1] = 50.5;
583   dm[2] = 32.25;
584
585   gMC->Gsvolu("TIVB","TUBE",idtmed[9],dm,3);
586
587   gMC->Gspos("TIVB",1,"TPVD",0.,0.,0.,0,"ONLY");
588
589   gMC->Gspos("TPVD",1,"ALIC",0.,0.,0.,0,"ONLY");
590
591   
592
593   
594
595   // --------------------------------------------------- 
596   //               volumes ordering 
597   // --------------------------------------------------- 
598   gMC->Gsord("TPSG", 6);
599  
600 } // end of function
601
602  
603  
604 //_____________________________________________________________________________
605 void AliTPCv3::DrawDetector()
606 {
607   //
608   // Draw a shaded view of the Time Projection Chamber version 1
609   //
610
611
612   // Set everything unseen
613   gMC->Gsatt("*", "seen", -1);
614   // 
615   // Set ALIC mother transparent
616   gMC->Gsatt("ALIC","SEEN",0);
617   //
618   // Set the volumes visible
619   gMC->Gsatt("TPC","SEEN",0);
620   gMC->Gsatt("TGAS","SEEN",0);
621   gMC->Gsatt("TPSG","SEEN",0);
622   gMC->Gsatt("TPHV","SEEN",1);
623   gMC->Gsatt("TPMH","SEEN",1);
624   gMC->Gsatt("TPEC","SEEN",0);
625   gMC->Gsatt("TRCS","SEEN",1);
626   gMC->Gsatt("TRCL","SEEN",1);
627   gMC->Gsatt("TPWL","SEEN",1);
628   gMC->Gsatt("TPWI","SEEN",1);
629   gMC->Gsatt("TPWS","SEEN",1);
630   gMC->Gsatt("TPW1","SEEN",1);
631   gMC->Gsatt("TPS1","SEEN",1);
632   gMC->Gsatt("TPS2","SEEN",1);
633   gMC->Gsatt("TPG1","SEEN",1);
634   gMC->Gsatt("TPG2","SEEN",1);
635   gMC->Gsatt("TPWC","SEEN",1);
636   gMC->Gsatt("TPSI","SEEN",1); 
637   gMC->Gsatt("TPSO","SEEN",1);
638   gMC->Gsatt("TPCO","SEEN",1);
639   gMC->Gsatt("TPOV","SEEN",1);
640   gMC->Gsatt("TPVD","SEEN",1);
641   //
642   gMC->Gdopt("hide", "on");
643   gMC->Gdopt("shad", "on");
644   gMC->Gsatt("*", "fill", 7);
645   gMC->SetClipBox(".");
646   gMC->SetClipBox("*", 0, 1000, -1000, 1000, -1000, 1000);
647   gMC->DefaultRange();
648   gMC->Gdraw("alic", 40, 30, 0, 12, 9.5, .025, .025);
649   gMC->Gdhead(1111, "Time Projection Chamber");
650   gMC->Gdman(18, 4, "MAN");
651   gMC->Gdopt("hide","off");
652 }
653
654 //_____________________________________________________________________________
655 void AliTPCv3::CreateMaterials()
656 {
657   //
658   // Define materials for version 2 of the Time Projection Chamber
659   //
660
661
662   //
663   // Increase maximum number of steps
664   gMC->SetMaxNStep(30000);
665   //
666   AliTPC::CreateMaterials();
667 }
668
669 //_____________________________________________________________________________
670 void AliTPCv3::Init()
671 {
672   //
673   // Initialises version 3 of the TPC after that it has been built
674   //
675   Int_t *idtmed = fIdtmed->GetArray()-399;
676
677   AliTPC::Init();
678
679   fIdSens1=gMC->VolId("TGAS"); // drift gas as a sensitive volume
680
681   gMC->SetMaxNStep(30000); // max. number of steps increased
682
683   gMC->Gstpar(idtmed[403],"LOSS",5);
684
685   printf("*** TPC version 3 initialized ***\n");
686   printf("Maximum number of steps = %d\n",gMC->GetMaxNStep());
687
688   //
689   
690 }
691
692 //_____________________________________________________________________________
693 void AliTPCv3::StepManager()
694 {
695   //
696   // Called for every step in the Time Projection Chamber
697   //
698
699   //
700   // parameters used for the energy loss calculations
701   //
702   const Float_t prim = 14.35; // number of primary collisions per 1 cm
703   const Float_t poti = 20.77e-9; // first ionization potential for Ne/CO2
704   const Float_t w_ion = 35.97e-9; // energy for the ion-electron pair creation 
705  
706  
707   const Float_t big = 1.e10;
708
709   Int_t id,copy;
710   TLorentzVector pos;
711   Float_t hits[4];
712   Int_t vol[2];  
713   TClonesArray &lhits = *fHits;
714   
715   vol[1]=0;
716   vol[0]=0;
717
718   //
719
720   gMC->SetMaxStep(big);
721   
722   if(!gMC->IsTrackAlive()) return; // particle has disappeared
723   
724   Float_t charge = gMC->TrackCharge();
725   
726   if(TMath::Abs(charge)<=0.) return; // take only charged particles
727   
728   
729   id=gMC->CurrentVolID(copy);
730   
731   // Check the sensitive volume
732   
733   if (id != fIdSens1) return;
734   
735   //
736   //  charged particle is in the sensitive volume
737   //
738   
739   if(gMC->TrackStep() > 0) {
740
741     
742     Int_t nel = (Int_t)(((gMC->Edep())-poti)/w_ion) + 1;
743     nel=TMath::Min(nel,300); // 300 electrons corresponds to 10 keV
744     
745     gMC->TrackPosition(pos);
746     hits[0]=pos[0];
747     hits[1]=pos[1];
748     hits[2]=pos[2];
749
750     //
751     // check the selected side of the TPC
752     //
753  
754     if(fSide && fSide*hits[2]<=0.) return;
755
756     hits[3]=(Float_t)nel;
757     
758     // Add this hit
759    
760     new(lhits[fNhits++]) AliTPChit(fIshunt,gAlice->CurrentTrack(),vol,hits);
761     
762   } 
763   
764   // Stemax calculation for the next step
765   
766   Float_t pp;
767   TLorentzVector mom;
768   gMC->TrackMomentum(mom);
769   Float_t ptot=mom.Rho();
770   Float_t beta_gamma = ptot/gMC->TrackMass();
771   
772   if(gMC->IdFromPDG(gMC->TrackPid()) <= 3 && ptot > 0.002)
773     { 
774       pp = prim*1.58; // electrons above 20 MeV/c are on the plateau!
775     }
776   else
777     {
778       pp=prim*BetheBloch(beta_gamma);    
779       if(TMath::Abs(charge) > 1.) pp *= (charge*charge);
780     }
781   
782   Float_t random[1];
783   gMC->Rndm(random,1); // good, old GRNDM from Geant3
784   
785   Double_t rnd = (Double_t)random[0];
786   
787   gMC->SetMaxStep(-TMath::Log(rnd)/pp);
788   
789 }
790
791 //_____________________________________________________________________________
792 Float_t AliTPCv3::BetheBloch(Float_t bg)
793 {
794   //
795   // Bethe-Bloch energy loss formula
796   //
797   const Double_t p1=0.76176e-1;
798   const Double_t p2=10.632;
799   const Double_t p3=0.13279e-4;
800   const Double_t p4=1.8631;
801   const Double_t p5=1.9479;
802
803   Double_t dbg = (Double_t) bg;
804
805   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
806
807   Double_t aa = TMath::Power(beta,p4);
808   Double_t bb = TMath::Power(1./dbg,p5);
809
810   bb=TMath::Log(p3+bb);
811   
812   return ((Float_t)((p2-aa-bb)*p1/aa));
813 }