Only use PDG codes and not GEANT ones
[u/mrichter/AliRoot.git] / TPC / AliTPCv2.cxx
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 /*
17 $Log$
18 Revision 1.17  1999/10/08 06:27:23  fca
19 Corrected bug in the HV degrader geometry, thanks to G.Tabary
20
21 Revision 1.16  1999/10/04 13:39:54  fca
22 Correct array index problem
23
24 Revision 1.15  1999/09/29 09:24:34  fca
25 Introduction of the Copyright and cvs Log
26
27 */
28
29 ///////////////////////////////////////////////////////////////////////////////
30 //                                                                           //
31 //  Time Projection Chamber version 2 -- detailed TPC and slow simulation    //
32 //                                                                           //
33 //Begin_Html
34 /*
35 <img src="picts/AliTPCv2Class.gif">
36 */
37 //End_Html
38 //                                                                           //
39 //                                                                           //
40 ///////////////////////////////////////////////////////////////////////////////
41 #include <stdlib.h>
42
43 #include <TMath.h>
44
45 #include "AliTPCv2.h"
46 #include "AliTPCD.h"
47 #include "AliRun.h"
48 #include "AliConst.h"
49 #include "AliPDG.h"
50
51 ClassImp(AliTPCv2)
52  
53 //_____________________________________________________________________________
54 AliTPCv2::AliTPCv2(const char *name, const char *title) :
55   AliTPC(name, title) 
56 {
57   //
58   // Standard constructor for Time Projection Chamber version 2
59   //
60   fIdSens1=0;
61   fIdSens2=0;
62   SetBufferSize(128000);
63 }
64  
65 //_____________________________________________________________________________
66 void AliTPCv2::CreateGeometry()
67 {
68   //
69   // Create the geometry of Time Projection Chamber version 2
70   //
71   //Begin_Html
72   /*
73     <img src="picts/AliTPCv2.gif">
74   */
75   //End_Html
76   //Begin_Html
77   /*
78     <img src="picts/AliTPCv2Tree.gif">
79   */
80   //End_Html
81
82   AliTPCParam * fTPCParam = &(fDigParam->GetParam());
83
84   Int_t *idtmed = fIdtmed->GetArray();
85
86   Float_t dm[21];
87   Int_t idrotm[120];
88
89   Int_t nRotMat = 0;
90
91   Int_t i,ifl1,ifl2;
92
93   Int_t nInnerSector = fTPCParam->GetNInnerSector()/2;
94   Int_t nOuterSector = fTPCParam->GetNOuterSector()/2;
95   
96   // --------------------------------------------------- 
97   //        sector specification check 
98   // --------------------------------------------------- 
99   if (fSecAL >= 0) {
100     ifl1 = 0;
101     
102     for (i = 0; i < 6; ++i) {
103       if (fSecLows[i] >= 0 && fSecLows[i] < 2*nInnerSector) {
104         ifl1 = 1;
105         printf("*** SECTOR %d selected\n",fSecLows[i]);
106       }
107     }
108
109   } else {
110     printf("*** ALL LOWER SECTORS SELECTED ***\n");
111     ifl1 = 1;
112   }
113
114   if (fSecAU >= 0) {
115     ifl2 = 0;
116     
117     for (i = 0; i < 12; ++i) {
118       if (fSecUps[i] > 2*nInnerSector-1 && 
119           fSecUps[i] < 2*(nInnerSector+nOuterSector)) {
120         ifl2 = 1;
121         printf("*** SECTOR %d selected\n",fSecUps[i]);
122       }
123     }
124     
125   } else {
126     printf("*** ALL UPPER SECTORS SELECTED ***\n");
127     ifl1 = 1;
128   }
129   
130   if (ifl1 == 0 && ifl2 == 0) {
131     printf("*** ERROR: AT LEAST ONE SECTOR MUST BE SPECIFIED ***\n");
132     printf("!!! PROGRAM STOPPED !!!\n");
133     exit(1);
134   }
135   
136   if ((fSecAL < 0 || fSecAU < 0) && fSens >= 0) {
137     printf("** ERROR: STRIPS CANNOT BE SPECIFIED FOR ALL SECTORS **\n");
138     printf("!!! PROGRAM STOPPED !!!\n");
139     exit(1);
140   }
141   
142   // ---------------------------------------------------- 
143   //          FIELD CAGE WITH ENDCAPS - G10
144   //          THIS IS ALSO A TPC MOTHER VOLUME 
145   // ---------------------------------------------------- 
146
147   dm[0] = 76.;
148   dm[1] = 278.;
149   dm[2] = 275.;
150
151   gMC->Gsvolu("TPC ", "TUBE", idtmed[8], dm, 3); 
152
153   //-----------------------------------------------------
154   //  Endcap cover c-fibre 0.86% X0
155   //-----------------------------------------------------
156
157   dm[0] = 78.;
158   dm[1] = 258.;
159   dm[2] = 0.95;
160
161   gMC->Gsvolu("TPEC","TUBE",idtmed[10],dm,3);
162
163   //-----------------------------------------------------
164   // Drift gas , leave 2 cm at the outer radius
165   // and inner raddius
166   //-----------------------------------------------------
167
168   dm[0] = 78.;
169   dm[1] = 258.;
170   dm[2] = 250.;
171
172   gMC->Gsvolu("TGAS", "TUBE", idtmed[3], dm, 3);
173
174   //------------------------------------------------------
175   //  membrane holder - carbon fiber
176   //------------------------------------------------------
177
178
179   gMC->Gsvolu("TPMH","TUBE",idtmed[6],dm,0);
180
181   dm[0] = 252.;
182   dm[1] = 258.;
183   dm[2] = 0.2;
184
185   gMC->Gsposp("TPMH",1,"TGAS",0.,0.,0.,0,"ONLY",dm,3);
186  
187   dm[0] = 78.;
188   dm[1] = 82.;
189   dm[2] = 0.1;
190
191   gMC->Gsposp("TPMH",2,"TGAS",0.,0.,0.,0,"ONLY",dm,3);
192
193   //----------------------------------------------------------
194   //  HV membrane - 25 microns of mylar
195   //----------------------------------------------------------
196
197   dm[0] = 82.;
198   dm[1] = 252.;
199   dm[2] = 0.00125;
200
201   gMC->Gsvolu("TPHV","TUBE",idtmed[5],dm,3);
202
203   gMC->Gspos("TPHV",1,"TGAS",0.,0.,0.,0,"ONLY");
204
205   gMC->Gspos("TGAS",1,"TPC ",0.,0.,0.,0,"ONLY");
206
207   //----------------------------------------------------------
208   // "side" gas volume, the same as the drift gas
209   // the readout chambers are placed there.  
210   //----------------------------------------------------------
211
212   dm[0] = 78.;
213   dm[1] = 258.;
214   dm[2] = 0.5*(275. - 250.);
215    
216   gMC->Gsvolu("TPSG", "TUBE", idtmed[2], dm, 3);
217
218   Float_t z_side = dm[2]; // 1/2 of the side gas thickness
219
220   //-----------------------------------------------------------
221   //   Readout chambers , 25% of X0, I use Al as the material
222   //-----------------------------------------------------------
223
224   Float_t InnerOpenAngle = fTPCParam->GetInnerAngle();
225   Float_t OuterOpenAngle = fTPCParam->GetOuterAngle();
226
227   Float_t InnerAngleShift = fTPCParam->GetInnerAngleShift();
228   Float_t OuterAngleShift = fTPCParam->GetOuterAngleShift();
229
230   Float_t InSecLowEdge = fTPCParam->GetInSecLowEdge();
231   Float_t InSecUpEdge =  fTPCParam->GetInSecUpEdge();
232
233   Float_t OuSecLowEdge = fTPCParam->GetOuSecLowEdge();
234   Float_t OuSecUpEdge = fTPCParam->GetOuSecUpEdge();
235
236
237   Float_t SecThick = 2.225; // Al
238
239   Float_t edge = fTPCParam->GetEdge();
240
241   //  S (Inner) sectors
242
243   dm[0] = InSecLowEdge*TMath::Tan(0.5*InnerOpenAngle)-edge;
244   dm[1] = InSecUpEdge*TMath::Tan(0.5*InnerOpenAngle)-edge;
245   dm[2] = 0.5*SecThick;
246   dm[3] = 0.5*(InSecUpEdge-InSecLowEdge);
247
248   Float_t xCenterS = InSecLowEdge+dm[3];
249
250   gMC->Gsvolu("TRCS", "TRD1", idtmed[0], dm, 4); 
251
252   //  L (Outer) sectors
253
254   dm[0] = OuSecLowEdge*TMath::Tan(0.5*OuterOpenAngle)-edge;
255   dm[1] = OuSecUpEdge*TMath::Tan(0.5*OuterOpenAngle)-edge;
256   dm[2] = 0.5*SecThick;
257   dm[3] = 0.5*(OuSecUpEdge-OuSecLowEdge);
258
259   Float_t xCenterL = OuSecLowEdge+dm[3];  
260
261   gMC->Gsvolu("TRCL", "TRD1", idtmed[0], dm, 4);
262
263   Float_t z1 = -z_side + SecThick*0.5;
264
265   //------------------------------------------------------------------
266   // S sectors - "gas sectors" (TRD1)
267   //------------------------------------------------------------------
268
269   dm[0] = InSecLowEdge*TMath::Tan(0.5*InnerOpenAngle)-0.01;
270   dm[1] = InSecUpEdge*TMath::Tan(0.5*InnerOpenAngle)-0.01;
271   dm[2] = 0.5*(250. - 0.001);
272   dm[3] = 0.5*(InSecUpEdge-InSecLowEdge);  
273
274   gMC->Gsvolu("TSGA", "TRD1", idtmed[4], dm, 4); // sensitive
275
276   // ------------------------------------------------------------- 
277   //  Only for the debugging purpose and resolution calculation 
278   //  Sensitive strips at the pad-row center 
279   // ------------------------------------------------------------- 
280
281   Int_t ns;
282
283   if(fSens>=0){
284
285     Float_t r1,r2,zz;
286
287     Float_t StripThick = 0.01; // 100 microns
288     Float_t dead = fTPCParam->GetDeadZone();
289
290     gMC->Gsvolu("TSST", "TRD1", idtmed[4], dm, 0);
291
292     dm[2] = 0.5*(250. - 0.002);
293     dm[3] = 0.5 * StripThick;
294
295
296     for (ns = 0; ns < fTPCParam->GetNRowLow(); ns++) {
297
298       r1 = fTPCParam->GetPadRowRadiiLow(ns);
299       r2 = r1 + StripThick;     
300       dm[0] = r1 * TMath::Tan(0.5*InnerOpenAngle) - dead;
301       dm[1] = r2 * TMath::Tan(0.5*InnerOpenAngle) - dead;
302
303       zz = -InSecLowEdge -0.5*(InSecUpEdge-InSecLowEdge);
304       zz += r1;
305       zz += dm[3];
306
307       gMC->Gsposp("TSST", ns+1, "TSGA", 0., 0., zz, 0, "ONLY", dm, 4);
308
309     
310     }   
311
312     gMC->Gsord("TSGA", 3); 
313   
314   } // if strips selected
315
316   
317   //-----------------------------------------------------------------
318   //  L sectors - "gas sectors" (PGON to avoid overlaps)
319   //-----------------------------------------------------------------
320
321   dm[0] = 360.*kDegrad - 0.5*OuterOpenAngle;
322   dm[0] *= kRaddeg;
323   dm[0] = (Float_t)TMath::Nint(dm[0]);
324
325   dm[1] = OuterOpenAngle*kRaddeg;
326   dm[1] = (Float_t)TMath::Nint(dm[1]);
327
328   dm[2] = 1.;
329   dm[3] = 4.;
330
331   dm[4] = 0.002;
332   dm[5] = OuSecLowEdge;
333   dm[6] = 252.*TMath::Cos(0.5*OuterOpenAngle)-0.002;
334
335   dm[7] = dm[4]+0.2;
336   dm[8] = dm[5];
337   dm[9] = dm[6];
338
339   dm[10] = dm[7];
340   dm[11] = OuSecLowEdge;
341   dm[12] = OuSecUpEdge;
342
343   dm[13] = 250.;
344   dm[14] = dm[11];
345   dm[15] = dm[12];
346
347   gMC->Gsvolu("TLGA","PGON",idtmed[4],dm,16);  
348
349   if (fSens >= 0) {
350     
351     Float_t rmax = dm[6];
352     Float_t r1,r2;
353     Float_t dead = fTPCParam->GetDeadZone();
354
355     Float_t StripThick = 0.01; // 100 microns
356
357     gMC->Gsvolu("TLST", "PGON", idtmed[4], dm, 0);
358
359     dm[0] = 360.*kDegrad - 0.5*OuterOpenAngle;
360     dm[0] *= kRaddeg;
361     dm[0] = (Float_t)TMath::Nint(dm[0]);
362
363     dm[1] = OuterOpenAngle*kRaddeg;
364     dm[1] = (Float_t)TMath::Nint(dm[1]);
365
366     dm[2] = 1.;
367     dm[3] = 2.;
368
369     dm[7] = 250.;
370
371     Float_t xx = dead/TMath::Tan(0.5*OuterOpenAngle);
372
373     for(ns=0;ns<fTPCParam->GetNRowUp();ns++){
374
375       r1 = fTPCParam->GetPadRowRadiiUp(ns)-xx;
376       r2 = r1 + StripThick;
377
378       dm[5] = r1;
379       dm[6] = r2;
380
381       dm[8] = r1;
382       dm[9] = r2;
383
384       if(r2+xx < rmax){
385         dm[4] = 0.002;
386       }
387       else{
388         dm[4] = 0.202;
389       }
390
391       gMC->Gsposp("TLST",ns+1,"TLGA",xx,0.,0.,0,"ONLY",dm,10);
392
393     }    
394
395         gMC->Gsord("TLGA", 4);
396
397   } // if strips selected
398
399   //------------------------------------------------------------------
400   // Positioning of the S-sector readout chambers
401   //------------------------------------------------------------------
402
403   Float_t zs = 0.5*(250.+0.002);
404
405   Float_t theta1,theta2,theta3;
406   Float_t phi1,phi2,phi3;
407   Float_t alpha;
408   Float_t x,y;
409
410   for(ns=0;ns<nInnerSector;ns++){
411     
412     phi1 = ns * InnerOpenAngle + 270.*kDegrad + InnerAngleShift;
413     phi1 *= kRaddeg; // in degrees
414
415     phi1 = (Float_t)TMath::Nint(phi1);
416
417     if (phi1 > 360.) phi1 -= 360.;
418       
419     theta1 = 90.;
420     phi2   = 90.;
421     theta2 = 180.;
422     phi3   = ns * InnerOpenAngle + InnerAngleShift;
423     phi3 *= kRaddeg; // in degrees
424
425     phi3 = (Float_t)TMath::Nint(phi3);
426       
427     if(phi3 > 360.) phi3 -= 360.;
428
429     theta3 = 90.;
430
431     alpha = phi3*kDegrad;
432
433     x = xCenterS * TMath::Cos(alpha);
434     y = xCenterS * TMath::Sin(alpha); 
435  
436     AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3);  
437      
438     gMC->Gspos("TRCS", ns+1, "TPSG", x, y, z1, idrotm[nRotMat], "ONLY");
439
440     if(fSecAL < 0){
441
442       //---------------------------------------------------------------
443       //  position all sectors
444       //---------------------------------------------------------------
445
446       gMC->Gspos("TSGA",ns+1,"TGAS",x,y,zs,idrotm[nRotMat], "ONLY");
447       gMC->Gspos("TSGA",ns+1+nInnerSector,"TGAS",x,y,-zs,idrotm[nRotMat], "ONLY");
448     }
449
450     else{
451
452       //---------------------------------------------------------------
453       //  position selected sectors
454       //---------------------------------------------------------------
455
456       for(Int_t sel=0;sel<6;sel++){
457
458         if(fSecLows[sel] == ns){
459         gMC->Gspos("TSGA", ns+1, "TGAS", x, y, zs, idrotm[nRotMat], "ONLY");
460         }
461         else if(fSecLows[sel] == ns+nInnerSector){
462         gMC->
463           Gspos("TSGA",ns+1+nInnerSector,"TGAS", x, y,-zs,idrotm[nRotMat],"ONLY");
464         }
465       }
466     }
467
468     nRotMat++;     
469
470   }
471     
472   //-------------------------------------------------------------------
473   //  Positioning of the L-sectors readout chambers
474   //-------------------------------------------------------------------
475     
476   for(ns=0;ns<nOuterSector;ns++){
477     phi1 = ns * OuterOpenAngle + 270.*kDegrad + OuterAngleShift;
478     phi1 *= kRaddeg; // in degrees
479
480     phi1 = (Float_t)TMath::Nint(phi1);
481     
482
483     if (phi1 > 360.) phi1 -= 360.;
484       
485     theta1 = 90.;
486     phi2   = 90.;
487     theta2 = 180.;
488     phi3   = ns * OuterOpenAngle+OuterAngleShift;
489     phi3 *= kRaddeg; // in degrees
490
491     phi3 = (Float_t)TMath::Nint(phi3);
492
493       
494     if(phi3 > 360.) phi3 -= 360.;
495
496     theta3 = 90.;
497
498     alpha = phi3*kDegrad;
499
500     x = xCenterL * TMath::Cos(alpha);
501     y = xCenterL * TMath::Sin(alpha); 
502  
503     AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3);  
504      
505
506     gMC->Gspos("TRCL", ns+1, "TPSG", x, y, z1, idrotm[nRotMat], "ONLY"); 
507
508     nRotMat++;   
509
510   }
511
512   //-------------------------------------------------------------------
513   // Positioning of the L-sectors (gas sectors)
514   //-------------------------------------------------------------------
515
516   for(ns=0;ns<nOuterSector;ns++){
517
518      phi1 = ns*OuterOpenAngle + OuterAngleShift;
519      phi1 *= kRaddeg;
520     
521      phi1 = (Float_t)TMath::Nint(phi1);
522      if(phi1>360.) phi1 -= 360.;
523
524      theta1 = 90.;
525
526      phi2 = 90. + phi1;
527      if(phi2>360.) phi2 -= 360.;
528
529      theta2 = 90.; 
530
531      phi3 = 0.;
532      theta3 = 0.;
533
534      if(fSecAU < 0) {
535
536        //--------------------------------------------------------------
537        //  position all sectors
538        //--------------------------------------------------------------
539
540        AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3); 
541
542        gMC->Gspos("TLGA",ns+1,"TGAS" ,0.,0.,0.,idrotm[nRotMat],"ONLY");
543
544        nRotMat++;
545    
546        // reflection !!
547
548        phi3 = 0.;
549        theta3 = 180.;
550      
551        AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3);
552   
553        gMC->Gspos("TLGA",ns+1+nOuterSector,"TGAS" ,0.,0.,0.,idrotm[nRotMat],"ONLY");
554           
555        nRotMat++;
556      }
557
558      else{
559
560        //---------------------------------------------------------------
561        //  position selected sectors
562        //---------------------------------------------------------------
563
564        for(Int_t sel=0;sel<12;sel++){
565
566          if(fSecUps[sel] == ns+2*nInnerSector){
567            
568           AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3); 
569           gMC->Gspos("TLGA",ns+1,"TGAS" ,0.,0.,0.,idrotm[nRotMat],"ONLY");
570           nRotMat++; 
571
572          }
573          else if(fSecUps[sel] == ns+2*nInnerSector+nOuterSector){
574
575            // reflection
576
577            phi3 = 0.;
578            theta3 = 180.;
579
580            AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3);
581            gMC->
582            Gspos("TLGA",ns+1+nOuterSector,"TGAS" ,0.,0.,0.,idrotm[nRotMat],"ONLY");
583            nRotMat++; 
584
585          }
586
587        }
588
589      }
590
591   }
592   
593   Float_t z0 = z_side - 0.95;
594
595   gMC->Gspos("TPEC",1,"TPSG",0.,0.,z0,0,"ONLY");
596
597   // ========================================================== 
598   //                  wheels 
599   // ========================================================== 
600
601   //
602   //  auxilary structures
603   //
604
605
606   gMC->Gsvolu("TPWI","TUBE",idtmed[24],dm,0); // "air" 
607
608   // ---------------------------------------------------------- 
609   //       Large wheel -> positioned in the TPC 
610   // ---------------------------------------------------------- 
611   
612
613   z0 = 263.5; // TPC length - 1/2 spoke wheel width
614
615   dm[0] = 258.;
616   dm[1] = 278.;
617   dm[2] = 11.5;
618   
619   gMC->Gsvolu("TPWL", "TUBE", idtmed[0], dm, 3); 
620
621   dm[0] = dm[0]+2.;
622   dm[1] = 278.;
623   dm[2] = dm[2]-2.;
624
625   gMC->Gsposp("TPWI",1,"TPWL",0.,0.,0.,0,"ONLY",dm,3);
626
627   gMC->Gspos("TPWL", 1, "TPC ", 0, 0, z0, 0, "ONLY");
628   gMC->Gspos("TPWL", 2, "TPC ", 0, 0, -z0, 0, "ONLY");
629
630   //
631   //  Outer vessel + CO2 HV degrader
632   //
633
634   dm[0] = 260.;
635   dm[1] = 278.;
636   dm[2] = 252.;
637
638   gMC->Gsvolu("TPCO","TUBE",idtmed[12],dm,3);
639
640   dm[0] = 275.;
641   dm[1] = 278.;
642   
643   gMC->Gsvolu("TPOV","TUBE",idtmed[10],dm,3);
644
645   gMC->Gspos("TPOV",1,"TPCO",0.,0.,0.,0,"ONLY");
646
647
648   // G10 plugs
649
650   dm[0] = 258.;
651   dm[1] = 260.;
652   dm[2] = 1.;
653
654   gMC->Gsvolu("TPG1","TUBE",idtmed[8],dm,3);
655   gMC->Gspos("TPG1",1,"TPCO",0.,0.,251.,0,"ONLY");
656   gMC->Gspos("TPG1",2,"TPCO",0.,0.,-251.,0,"ONLY");  
657
658   gMC->Gspos("TPCO",1,"TPC ",0.,0.,0.,0,"ONLY");
659
660
661   //----------------------------------------------------------
662   //  Small wheel -> positioned in "side gas
663   //----------------------------------------------------------
664
665   dm[0] = 78.;
666   dm[1] = 82.;
667   dm[2] = 11.5;
668
669   gMC->Gsvolu("TPWS", "TUBE", idtmed[0], dm, 3);
670
671   dm[0] = 78.;
672   dm[1] = dm[1]-2;
673   dm[2] = dm[2]-2.;
674
675   gMC->Gsvolu("TPW1", "TUBE", idtmed[2], dm, 3);
676   
677   gMC->Gspos("TPW1", 1, "TPWS", 0., 0., 0., 0, "ONLY");
678
679   z0 = 1.; // spoke wheel is shifted w.r.t. center of the "side gas"
680
681   gMC->Gspos("TPWS", 1, "TPSG", 0, 0, z0, 0, "ONLY");
682
683
684   // to avoid overlaps
685
686   dm[0] = 76.;
687   dm[1] = 78.;
688   dm[2] = 11.5;
689
690   gMC->Gsvolu("TPS1","TUBE",idtmed[0],dm,3);
691
692   dm[2] = 9.5;
693
694   gMC->Gsvolu("TPS2","TUBE",idtmed[24],dm,3);
695
696   gMC->Gspos("TPS2",1,"TPS1",0.,0.,0.,0,"ONLY");
697
698   z0= 263.5;
699   
700   gMC->Gspos("TPS1",1,"TPC ",0.,0.,z0,0,"ONLY");
701   gMC->Gspos("TPS1",2,"TPC ",0.,0.,-z0,0,"ONLY");
702
703   // G10 plug
704
705   dm[0] = 76.;
706   dm[1] = 78.;
707   dm[2] = 1.;
708
709   gMC->Gsvolu("TPG2","TUBE",idtmed[8],dm,3);
710
711   z0 = 251.;
712
713   gMC->Gspos("TPG2",1,"TPC ",0.,0.,z0,0,"ONLY");
714   gMC->Gspos("TPG2",2,"TPC ",0.,0.,-z0,0,"ONLY");
715
716
717   //---------------------------------------------------------
718   //  central wheel  6 (radial direction) x 4 (along z) cm2
719   //---------------------------------------------------------
720
721   dm[0] = 140.;
722   dm[1] = 146.;
723   dm[2] = 2.;
724
725   gMC->Gsvolu("TPWC","TUBE",idtmed[0],dm,3);
726
727   dm[0] = dm[0] + 2.;
728   dm[1] = dm[1] - 2.;
729   dm[2] = dm[2] - 1.;
730
731   gMC->Gsposp("TPWI",2,"TPWC",0.,0.,0.,0,"ONLY",dm,3);
732
733   z0 = z_side - 1.9 - 2.;
734
735   gMC->Gspos("TPWC",1,"TPSG",0.,0.,z0,0,"ONLY");
736
737   //
738
739   gMC->Gsvolu("TPSE","BOX ",idtmed[24],dm,0); // "empty" part of the spoke 
740
741  
742   //---------------------------------------------------------
743   //  inner spokes (nSectorInner)
744   //---------------------------------------------------------
745
746   dm[0] = 0.5*(139.9-82.1);
747   dm[1] = 3.;
748   dm[2] = 2.;
749
750   Float_t x1 = dm[0]+82.;
751
752   gMC->Gsvolu("TPSI","BOX",idtmed[0],dm,3);
753
754   dm[1] = dm[1]-1.;
755   dm[2] = dm[2]-1.;
756
757   gMC->Gsposp("TPSE",1,"TPSI",0.,0.,0.,0,"ONLY",dm,3);
758
759   for(ns=0;ns<nInnerSector;ns++){
760
761     phi1 = 0.5*InnerOpenAngle + ns*InnerOpenAngle + InnerAngleShift;
762     theta1=90.;
763     phi1 *=kRaddeg;
764
765     phi1 = (Float_t)TMath::Nint(phi1);
766     if(phi1>360.) phi1 -= 360.;    
767
768     phi2 = phi1+90.;
769     if(phi2>360.) phi2 -= 360.;
770     theta2=90.;
771     phi3=0.;
772     theta3=0.;
773
774     alpha = phi1 * kDegrad;
775     x     = x1 * TMath::Cos(alpha);
776     y     = x1 * TMath::Sin(alpha);    
777
778    AliMatrix(idrotm[nRotMat],theta1,phi1,theta2,phi2,theta3,phi3);
779
780    gMC->Gspos("TPSI",ns+1,"TPSG",x,y,z0,idrotm[nRotMat],"ONLY");  
781
782    nRotMat++;
783
784   }
785
786   //-------------------------------------------------------------
787   // outer spokes (nSectorOuter)
788   //-------------------------------------------------------------
789
790   dm[0] = 0.5*(257.9-146.1);
791   dm[1] = 3.;
792   dm[2] = 2.;
793
794   x1 = dm[0] + 146.;
795
796   gMC->Gsvolu("TPSO","BOX ",idtmed[0],dm,3);
797
798   dm[1] = dm[1] - 1.;
799   dm[2] = dm[2] - 1.;
800
801   gMC->Gsposp("TPSE",2,"TPSO",0.,0.,0.,0,"ONLY",dm,3);
802
803   for(ns=0;ns<nOuterSector;ns++){
804
805     phi1 = 0.5*OuterOpenAngle + ns*OuterOpenAngle + OuterAngleShift;
806     theta1=90.;
807     phi1 *=kRaddeg;
808
809     phi1 = (Float_t)TMath::Nint(phi1);
810     if(phi1>360.) phi1 -= 360.;
811
812     phi2 = phi1+90.;
813     if(phi2>360.) phi2 -= 360.;
814     theta2=90.;
815     phi3=0.;
816     theta3=0.;
817
818     alpha = phi1 * kDegrad;
819     x     = x1 * TMath::Cos(alpha);
820     y     = x1 * TMath::Sin(alpha);    
821
822    AliMatrix(idrotm[nRotMat],theta1,phi1,theta2,phi2,theta3,phi3);
823
824    gMC->Gspos("TPSO",ns+1,"TPSG",x,y,z0,idrotm[nRotMat],"ONLY");  
825
826    nRotMat++;
827
828   }  
829   
830
831   
832   // -------------------------------------------------------- 
833   //         put the readout chambers into the TPC 
834   // -------------------------------------------------------- 
835
836   theta1 = 90.;
837   phi1   = 0.;
838   theta2 = 90.;
839   phi2   = 270.;
840   theta3 = 180.;
841   phi3   = 0.;
842   
843   AliMatrix(idrotm[nRotMat], theta1, phi1, theta2, phi2, theta3, phi3);
844   
845   z0 = z_side + 250.;
846   
847   gMC->Gspos("TPSG", 1, "TPC ", 0, 0, z0, 0, "ONLY");
848   gMC->Gspos("TPSG", 2, "TPC ", 0, 0, -z0, idrotm[nRotMat], "ONLY");
849   
850   gMC->Gspos("TPC ", 1, "ALIC", 0, 0, 0, 0, "ONLY");
851
852   //----------------------------------------------------
853   //  Inner vessel and HV degrader
854   //----------------------------------------------------
855
856   dm[0] = 0.;
857   dm[1] = 360.;
858   dm[2] = 4.;
859   
860   dm[3] = -250.;
861   dm[4] = 74.4;
862   dm[5] = 76.;
863
864   dm[6] = -64.5;
865   dm[7] = 50.;
866   dm[8] = 76.;
867
868   dm[9] = 64.5;
869   dm[10] = 50.;
870   dm[11] = 76.;
871
872   dm[12] = 250.;
873   dm[13] = 74.4;
874   dm[14] = 76.;
875
876   gMC->Gsvolu("TPVD", "PCON", idtmed[12], dm, 15); // CO2
877
878   // cone parts
879
880   dm[0] = 0.;
881   dm[1] = 360.;
882   dm[2] = 2.;
883
884   dm[3] = 64.5;
885   dm[4] = 50.;
886   dm[5] = 51.6;
887  
888   dm[6] = 250.;
889   dm[7] = 74.4;
890   dm[8] = 76.;
891
892
893   gMC->Gsvolu("TIVC","PCON",idtmed[11],dm,9); // C-fibre
894
895   gMC->Gspos("TIVC",1,"TPVD",0.,0.,0.,0,"ONLY");
896   gMC->Gspos("TIVC",2,"TPVD",0.,0.,0.,idrotm[nRotMat],"ONLY");
897
898   // barrel part
899
900   dm[0] = 50.;
901   dm[1] = 50.5;
902   dm[2] = 32.25;
903
904   gMC->Gsvolu("TIVB","TUBE",idtmed[9],dm,3);
905
906   gMC->Gspos("TIVB",1,"TPVD",0.,0.,0.,0,"ONLY");
907
908   gMC->Gspos("TPVD",1,"ALIC",0.,0.,0.,0,"ONLY");
909
910
911   // --------------------------------------------------- 
912   //               volumes ordering 
913   // ---------------------------------------------------
914
915   gMC->Gsord("TGAS", 6);
916   gMC->Gsord("TPSG", 6); 
917
918   
919
920 } // end of function
921  
922 //_____________________________________________________________________________
923 void AliTPCv2::DrawDetector()
924 {
925   //
926   // Draw a shaded view of the Time Projection Chamber version 1
927   //
928
929   // Set everything unseen
930   gMC->Gsatt("*", "seen", -1);
931   // 
932   // Set ALIC mother transparent
933   gMC->Gsatt("ALIC","SEEN",0);
934   //
935   // Set the volumes visible
936   gMC->Gsatt("TPC","SEEN",0);
937   gMC->Gsatt("TGAS","SEEN",0);
938   gMC->Gsatt("TPSG","SEEN",0);
939   gMC->Gsatt("TPHV","SEEN",1);
940   gMC->Gsatt("TPMH","SEEN",1);
941   gMC->Gsatt("TPEC","SEEN",0);
942   gMC->Gsatt("TRCS","SEEN",1);
943   gMC->Gsatt("TRCL","SEEN",1);
944   gMC->Gsatt("TPWL","SEEN",1);
945   gMC->Gsatt("TPWI","SEEN",1);
946   gMC->Gsatt("TPWS","SEEN",1);
947   gMC->Gsatt("TPW1","SEEN",1);
948   gMC->Gsatt("TPS1","SEEN",1);
949   gMC->Gsatt("TPS2","SEEN",1);
950   gMC->Gsatt("TPG1","SEEN",1);
951   gMC->Gsatt("TPG2","SEEN",1);
952   gMC->Gsatt("TPWC","SEEN",1);
953   gMC->Gsatt("TPSI","SEEN",1); 
954   gMC->Gsatt("TPSO","SEEN",1);
955   gMC->Gsatt("TPCO","SEEN",1);
956   gMC->Gsatt("TPOV","SEEN",1);
957   gMC->Gsatt("TPVD","SEEN",1);
958   //
959   gMC->Gdopt("hide", "on");
960   gMC->Gdopt("shad", "on");
961   gMC->Gsatt("*", "fill", 7);
962   gMC->SetClipBox(".");
963   gMC->SetClipBox("*", 0, 1000, -1000, 1000, -1000, 1000);
964   gMC->DefaultRange();
965   gMC->Gdraw("alic", 40, 30, 0, 12, 9.5, .025, .025);
966   gMC->Gdhead(1111, "Time Projection Chamber");
967   gMC->Gdman(18, 4, "MAN");
968   gMC->Gdopt("hide","off");
969 }
970
971 //_____________________________________________________________________________
972 void AliTPCv2::CreateMaterials()
973 {
974   //
975   // Define materials for version 2 of the Time Projection Chamber
976   //
977
978   //
979   // Increase maximum number of steps
980   gMC->SetMaxNStep(30000);
981   //
982   AliTPC::CreateMaterials();
983 }
984
985 //_____________________________________________________________________________
986 void AliTPCv2::Init()
987 {
988   //
989   // Initialises version 2 of the TPC after that it has been built
990   //
991   Int_t *idtmed = fIdtmed->GetArray()-399;
992   AliTPC::Init();
993   fIdSens1=gMC->VolId("TLGA"); // L-sector
994   fIdSens2=gMC->VolId("TSGA"); // S-sector 
995   fIdSens3=gMC->VolId("TSST"); // strip - S-sector (not always used)
996   fIdSens4=gMC->VolId("TLST"); // strip - S-sector (not always used)
997
998   gMC->SetMaxNStep(30000); // max. number of steps increased
999
1000   gMC->Gstpar(idtmed[403],"LOSS",5);
1001
1002   printf("*** TPC version 2 initialized ***\n");
1003   printf("Maximum number of steps = %d\n",gMC->GetMaxNStep());
1004
1005   //
1006   
1007 }
1008
1009 //_____________________________________________________________________________
1010 void AliTPCv2::StepManager()
1011 {
1012   //
1013   // Called for every step in the Time Projection Chamber
1014   //
1015
1016   //
1017   // parameters used for the energy loss calculations
1018   //
1019   const Float_t prim = 14.35; // number of primary collisions per 1 cm
1020   const Float_t poti = 20.77e-9; // first ionization potential for Ne/CO2
1021   const Float_t w_ion = 35.97e-9; // energy for the ion-electron pair creation 
1022
1023   //  const Float_t prim = 17.65;
1024   //  const Float_t poti = 19.02e-9;
1025   // const Float_t w_ion = 33.06e-9;
1026  
1027  
1028   const Float_t big = 1.e10;
1029
1030   Int_t id,copy;
1031   Float_t hits[4];
1032   Int_t vol[2];  
1033   TClonesArray &lhits = *fHits;
1034   TLorentzVector pos;
1035
1036   AliTPCParam *fTPCParam = &(fDigParam->GetParam());
1037   
1038   vol[1]=0;
1039
1040   //
1041
1042   gMC->SetMaxStep(big);
1043   
1044   if(!gMC->IsTrackAlive()) return; // particle has disappeared
1045   
1046   Float_t charge = gMC->TrackCharge();
1047   
1048   if(TMath::Abs(charge)<=0.) return; // take only charged particles
1049   
1050   
1051   id=gMC->CurrentVolID(copy);
1052   
1053   // Check the sensitive volume
1054   
1055   if(id == fIdSens1) 
1056     {
1057       vol[0] = copy + fTPCParam->GetNInnerSector()-1; // L-sector number
1058     }
1059   else if(id == fIdSens2) 
1060     {
1061       vol[0] = copy-1; // S-sector number 
1062     }
1063   else if(id == fIdSens3 && gMC->IsTrackEntering())
1064     {
1065       vol[1] = copy-1;  // row number  
1066       id = gMC->CurrentVolOffID(1,copy);
1067       vol[0] = copy-1; // sector number (S-sector)
1068       
1069       gMC->TrackPosition(pos);
1070       hits[0]=pos[0];
1071       hits[1]=pos[1];
1072       hits[2]=pos[2];
1073       hits[3]=0.; // this hit has no energy loss
1074       new(lhits[fNhits++]) AliTPChit(fIshunt,gAlice->CurrentTrack(),vol,hits);
1075     }
1076   else if(id == fIdSens4 && gMC->IsTrackEntering())
1077     {
1078       vol[1] = copy-1; // row number 
1079       id = gMC->CurrentVolOffID(1,copy);
1080       vol[0] = copy+fTPCParam->GetNInnerSector()-1; // sector number (L-sector)
1081       
1082       gMC->TrackPosition(pos);
1083       hits[0]=pos[0];
1084       hits[1]=pos[1];
1085       hits[2]=pos[2];
1086       hits[3]=0.; // this hit has no energy loss
1087       new(lhits[fNhits++]) AliTPChit(fIshunt,gAlice->CurrentTrack(),vol,hits);
1088     }
1089   else return;
1090   
1091   //
1092   //  charged particle is in the sensitive volume
1093   //
1094   
1095   if(gMC->TrackStep() > 0) {
1096     
1097     Int_t nel = (Int_t)(((gMC->Edep())-poti)/w_ion) + 1;
1098     nel=TMath::Min(nel,300); // 300 electrons corresponds to 10 keV
1099     
1100     gMC->TrackPosition(pos);
1101     hits[0]=pos[0];
1102     hits[1]=pos[1];
1103     hits[2]=pos[2];
1104     hits[3]=(Float_t)nel;
1105     
1106     // Add this hit
1107     
1108     new(lhits[fNhits++]) AliTPChit(fIshunt,gAlice->CurrentTrack(),vol,hits);
1109     
1110   } 
1111   
1112   // Stemax calculation for the next step
1113   
1114   Float_t pp;
1115   TLorentzVector mom;
1116   gMC->TrackMomentum(mom);
1117   Float_t ptot=mom.Rho();
1118   Float_t beta_gamma = ptot/gMC->TrackMass();
1119   
1120   Int_t pid=gMC->TrackPid();
1121   if((pid==kElectron || pid==kPositron || pid==kGamma) && ptot > 0.002)
1122     { 
1123       pp = prim*1.58; // electrons above 20 MeV/c are on the plateau!
1124     }
1125   else
1126     {
1127       pp=prim*BetheBloch(beta_gamma);    
1128       if(TMath::Abs(charge) > 1.) pp *= (charge*charge);
1129     }
1130   
1131   Float_t random[1];
1132   gMC->Rndm(random,1); // good, old GRNDM from Geant3
1133   
1134   Double_t rnd = (Double_t)random[0];
1135   
1136   gMC->SetMaxStep(-TMath::Log(rnd)/pp);
1137   
1138 }
1139
1140 //_____________________________________________________________________________
1141 Float_t AliTPCv2::BetheBloch(Float_t bg)
1142 {
1143   //
1144   // Bethe-Bloch energy loss formula
1145   //
1146   const Double_t p1=0.76176e-1;
1147   const Double_t p2=10.632;
1148   const Double_t p3=0.13279e-4;
1149   const Double_t p4=1.8631;
1150   const Double_t p5=1.9479;
1151
1152   Double_t dbg = (Double_t) bg;
1153
1154   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
1155
1156   Double_t aa = TMath::Power(beta,p4);
1157   Double_t bb = TMath::Power(1./dbg,p5);
1158
1159   bb=TMath::Log(p3+bb);
1160   
1161   return ((Float_t)((p2-aa-bb)*p1/aa));
1162 }