FairRoot/PandaRoot
create3StationsGem.C
Go to the documentation of this file.
1 #include "iomanip.h"
2 
3 int create3StationsGem()
4 {
5  // Gem disk geometry parameters (R.Karabowicz)
6  //----------------------------------------------------------------------------------------------------------------------------------------------
7  // Done Some Changes By Nazila Divani-Veis ( since May 2015 ... )
8  //----------------------------------------------------------------------------------------------------------------------------------------------
9  // Try to get close to the CAD Geometry ( After could not use CadConverter ) //Some parameters from DL.
10  //----------------------------------------------------------------------------------------------------------------------------------------------
11  // Units are in [cm]
12  //---------------------------------------------------------------------------------------------------------------------------------------------
13 
14  const Int_t kNofDisks = 3; //number of GEM stations in form of disk
15 
16  const Double_t kDiskVolInnerRadius[kNofDisks] = { 4.50, 4.50, 4.50 }; // InnerRadius for 3 Stations in form of disk
17  const Double_t kDiskVolOuterRadius[kNofDisks] = { 45.00, 56.00, 74.0 }; // OuterRadius for 3 Stations in form of disk
18  const Double_t kDiskZPosition [kNofDisks] = { 119.40, 155.40, 188.50 }; // ZPosition for 3 Stations in form of disk
19 
20 
21  const Int_t kDiskNFoils [kNofDisks] = { 2 , 2 , 4 }; // For modifying misalignment
22  const Double_t kMiddleROBarHfTh[kNofDisks] = { 2.35, 2.35, 2.35 }; // half thickness of space in the middle of foils [cm]
23  const Double_t kHalfStationThickness = 6.00;
24 
25  const Double_t carbonRingInnerRadius[kNofDisks] = { 40.0, 51.0, 69.0 };
26  const Double_t carbonRingOuterRadius[kNofDisks] = { 40.2, 51.2, 69.2 };
27  const Double_t carbonRingHalfThickness = 1.5;
28 
29  const Double_t copperRingInnerRadius[kNofDisks] = { 39.9, 50.9, 68.9 };
30  const Double_t copperRingOuterRadius[kNofDisks] = { 40.1, 51.1, 69.1 };
31  const Double_t copperRingHalfThickness = 3.75;
32 
33  const Double_t AlRingInnerRadius[kNofDisks] = { 40.1, 51.1, 69.1 };
34  const Double_t AlRingOuterRadius[kNofDisks] = { 44.7, 55.7, 73.7 };
35  const Double_t AlRingHalfThickness = 3.75;
36 
37  const Double_t SegmentHalfThickness = 0.25;
38  const Int_t NofSegments = 24;
39  const Double_t SegmentDeltaAngle = 360./(Double_t(NofSegments));
40  const Double_t FirstSegmentAngle = 7.5;
41  const TString newsegment [NofSegments] = { "seg1" ,"seg2" ,"seg3" ,"seg4" ,"seg5" ,"seg6" ,"seg7" ,"seg8" ,
42  "seg9" ,"seg10","seg11","seg12","seg13","seg14","seg15","seg16",
43  "seg17","seg18","seg19","seg20","seg21","seg22","seg23","seg24" };
44 
45  const Double_t moduleRingInnerRadius[kNofDisks] = { 40.15, 51.15, 69.15 };
46  const Double_t moduleRingOuterRadius[kNofDisks] = { 44.65, 55.65, 73.65 };
47  const Double_t moduleRingHalfThickness = 3.75;
48 
49  const Double_t moduleSegmentHalfThickness = 1.30;
50  const Int_t NofmoduleSegments = 24;
51  const Double_t moduleSegmentDeltaAngle = 360./(Double_t(NofmoduleSegments));
52  const Double_t FirstmoduleSegmentAngle = 0.0 ;
53  const TString newmodulesegment [NofmoduleSegments] = { "moduleseg1" ,"moduleseg2" ,"moduleseg3" ,"moduleseg4" ,"moduleseg5" ,"moduleseg6" ,"moduleseg7" ,"moduleseg8" ,
54  "moduleseg9" ,"moduleseg10","moduleseg11","moduleseg12","moduleseg13","moduleseg14","moduleseg15","moduleseg16",
55  "moduleseg17","moduleseg18","moduleseg19","moduleseg20","moduleseg21","moduleseg22","moduleseg23","moduleseg24" };
56 
57  const Double_t AlumiRingInnerRadius[kNofDisks] = { 44.7, 55.7, 73.7 };
58  const Double_t AlumiRingOuterRadius[kNofDisks] = { 45.0, 56.0, 74.0 };
59  const Double_t AlumiRingHalfThickness = 3.75;
60 
61  const Double_t coverRingInnerRadius[kNofDisks] = { 39.9, 50.9, 68.9 };
62  const Double_t coverRingOuterRadius[kNofDisks] = { 45.0, 56.0, 74.0 };
63  const Double_t coverRingHalfThickness = 0.2;
64 
65  const Double_t rcopperbarx = 5.00;
66  const Double_t rcopperbary = 1.40;
67  const Double_t rcopperbarHalfThickness = 17.0;
68 
69  const Double_t lcopperbarx = 5.00;
70  const Double_t lcopperbary = 1.00;
71  const Double_t lcopperbarHalfThickness = 16.0;
72 
73  //------------------------------------------------------------ main layers in shape of Disk----------------------------------------------------------------------
74  const Int_t kNofLayers = 47;
75  const TString kLayerName [kNofLayers] = { "window_foil_out_aluminium","WindowF_kapton",
76  "space", // these "spaces" belong to holding structure layers
77  "WindowF_aluminium","CathodeF_kapton","CathodeF_aluminium",
78  "space",
79  "Gem1_Sensor_GEMmixture", // sensitive layer
80  "space",
81  "Gem1F_copper","Gem1_kapton","Gem1B_copper",
82  "space",
83  "Gem2F_copper","Gem2_kapton","Gem2B_copper",
84  "space",
85  "Gem3F_copper","Gem3_kapton","Gem3B_copper",
86  "space",
87  "space",
88  "PadF_copper","Pad_kapton","PadB_copper",
89  "space",
90  "space",
91  "Gem4F_copper","Gem4_kapton","Gem4B_copper",
92  "space",
93  "Gem5F_copper","Gem5_kapton","Gem5B_copper",
94  "space",
95  "Gem6F_copper","Gem6_kapton","Gem6B_copper",
96  "space",
97  "Gem6_Sensor_GEMmixture", // sensitive layer
98  "space",
99  "CathodeB_aluminium","CathodeB_kapton","WindowB_aluminium",
100  "space",
101  "WindowB_kapton","window_foil_in_aluminium" };
102  const Double_t kLayerThick[kNofLayers] = { 0.0001, 0.0007, // 2 = 2 window
103  1.00, // windowF_ring_carbon=holding structure // +1 = 3 space
104  0.0001, 0.0007, 0.0001, // +3 = 6 window + cathode
105  0.80, // cathodeF_ring_GlassFiber=holding structure // +1 = 7 space
106  1.0020, // +1 = 8 SENSOR
107  0.0050, // gem_ring1_GlassFiber=holding structure // +1 = 9 space
108  0.0005, 0.0050, 0.0005, // +3 = 12 GEM Foil
109  0.0050, // gem_ring2_GlassFiber=holding structure // +1 = 13 space
110  0.0005, 0.0050, 0.0005, // +3 = 16 GEM Foil
111  0.0050, // gem_ring3_GlassFiber=holding structure // +1 = 17 space
112  0.0005, 0.0050, 0.0005, // +3 = 20 GEM Foil
113  0.0050, // gem_ring4_GlassFiber=holding structure // +1 = 21 space
114  0.10, // padplaneF_support_GlassFiber=holding structure // +1 = 22 space
115  0.001, 0.001, 0.001, // +3 = 25 PAD plane
116  0.10, // padplaneB_support_GlassFiber=holding structure // +1 = 26 space
117  0.0050, // gem_ring5_GlassFiber=holding structure // +1 = 27 space
118  0.0005, 0.0050, 0.0005, // +3 = 30 GEM Foil
119  0.0050, // gem_ring6_GlassFiber=holding structure // +1 = 31 space
120  0.0005, 0.0050, 0.0005, // +3 = 34 GEM Foil
121  0.0050, // gem_ring7_GlassFiber=holding structure // +1 = 35 space
122  0.0005, 0.0050, 0.0005, // +3 = 38 GEM Foil
123  0.0050, // gem_ring8_GlassFiber=holding structure // +1 = 39 space
124  1.0020, // +1 = 40 SENSOR
125  0.80, // cathodeB_ring_GlassFiber=holding structure // +1 = 41 space
126  0.0001, 0.0007, 0.0001, // +3 = 44 cathode + window
127  1.00, // windowB_ring_carbon=holding structure // +1 = 45 space
128  0.0007, 0.0001 }; // +2 = 47 window
129  const Double_t kDiskOuterRadius[kNofLayers][kNofDisks] = { 38.95, 49.95, 67.95, 39.40, 50.40, 68.40,
130  39.90, 50.90, 68.90,
131  38.45, 49.45, 67.45, 39.40, 50.40, 68.40, 38.45, 49.45, 67.45,
132  39.40, 50.40, 68.40,
133  38.45, 49.45, 67.45,
134  39.05, 50.40, 68.40,
135  38.10, 49.45, 67.45, 39.05, 50.40, 68.40, 38.10, 49.45, 67.45,
136  39.05, 50.40, 68.40,
137  38.10, 49.45, 67.45, 39.05, 50.40, 68.40, 38.10, 49.45, 67.45,
138  39.05, 50.40, 68.40,
139  38.10, 49.45, 67.45, 39.05, 50.40, 68.40, 38.10, 49.45, 67.45,
140  39.05, 50.40, 68.40,
141  45.00, 56.00, 74.00,
142  44.90, 55.90, 73.90, 44.90, 55.90, 73.90, 44.90, 55.90, 73.90,
143  45.00, 56.00, 74.00,
144  39.05, 50.40, 68.40,
145  38.10, 49.45, 67.45, 39.05, 50.40, 68.40, 38.10, 49.45, 67.45,
146  39.05, 50.40, 68.40,
147  38.10, 49.45, 67.45, 39.05, 50.40, 68.40, 38.10, 49.45, 67.45,
148  39.05, 50.40, 68.40,
149  38.10, 49.45, 67.45, 39.05, 50.40, 68.40, 38.10, 49.45, 67.45,
150  39.05, 50.40, 68.40,
151  38.45, 49.45, 67.45,
152  39.40, 50.40, 64.40,
153  38.45, 49.45, 67.45, 39.40, 50.40, 68.40, 38.45, 49.45, 67.45,
154  39.40, 50.40, 68.40,
155  39.40, 50.40, 68.40, 38.95, 49.45, 67.95 };
156  const Double_t kDiskInnerRadius[kNofLayers][kNofDisks] = { 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
157  4.50, 4.50, 4.50,
158  4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
159  4.50, 4.50, 4.50,
160  4.50, 4.50, 4.50,
161  4.50, 4.50, 4.50,
162  4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
163  4.50, 4.50, 4.50,
164  4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
165  4.50, 4.50, 4.50,
166  4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
167  4.50, 4.50, 4.50,
168  4.50, 4.50, 4.50,
169  4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
170  4.50, 4.50, 4.50,
171  4.50, 4.50, 4.50,
172  4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
173  4.50, 4.50, 4.50,
174  4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
175  4.50, 4.50, 4.50,
176  4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
177  4.50, 4.50, 4.50,
178  4.50, 4.50, 4.50,
179  4.50, 4.50, 4.50,
180  4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50, 4.50,
181  4.50, 4.50, 4.50,
182  4.50, 4.50, 4.50, 4.50, 4.50, 4.50 };
183 
184  const Double_t HoleTZ = 0.0; // (top and down) holes Translation parameters
185  const Double_t HoleTX = 0.0;
186  const Double_t HoleTY = 33.50;
187 
188 //----------------------- holding structure layers ------------------------------------------------------------------------------------------------
189  const Int_t NofHLayers = 14;
190  const TString HLayersName [NofHLayers] = { "windowF_ring_carbon",
191  "cathodeF_ring_GlassFiber",
192  "gem_ring1_GlassFiber",
193  "gem_ring2_GlassFiber",
194  "gem_ring3_GlassFiber",
195  "gem_ring4_GlassFiber",
196  "padplaneF_support_GlassFiber",
197  "padplaneB_support_GlassFiber",
198  "gem_ring5_GlassFiber",
199  "gem_ring6_GlassFiber",
200  "gem_ring7_GlassFiber",
201  "gem_ring8_GlassFiber",
202  "cathodeB_ring_GlassFiber",
203  "windowB_ring_carbon" };
204  const Double_t HLayersThick[NofHLayers] = { 1.00,
205  0.80,
206  0.0050,
207  0.0050,
208  0.0050,
209  0.0050,
210  0.10,
211  0.10,
212  0.0050,
213  0.0050,
214  0.0050,
215  0.0050,
216  0.80,
217  1.00 };
218  const Double_t HZPosition[NofHLayers][kNofDisks] = { -1.9426, -1.9426, -1.9426,
219  -1.5417, -1.5417, -1.5417,
220  -0.1372, -0.1372, -0.1372,
221  -0.1262, -0.1262, -0.1262,
222  -0.1152, -0.1152, -0.1152,
223  -0.1042, -0.1042, -0.1042,
224  -0.0517, -0.0517, -0.0517,
225  0.0517, 0.0517, 0.0517,
226  0.1042, 0.1042, 0.1042,
227  0.1152, 0.1152, 0.1152,
228  0.1262, 0.1262, 0.1262,
229  0.1372, 0.1372, 0.1372,
230  1.5417, 1.5417, 1.5417,
231  1.9426, 1.9426, 1.9426 };
232  const Double_t HOuterRadius[NofHLayers][kNofDisks] = { 39.90, 50.90, 68.90,
233  39.40, 50.40, 68.40,
234  39.05, 50.05, 68.05,
235  39.05, 50.05, 68.05,
236  39.05, 50.05, 68.05,
237  39.05, 50.05, 68.05,
238  45.00, 56.00, 74.00,
239  45.00, 56.00, 74.00,
240  39.05, 50.05, 68.05,
241  39.05, 50.05, 68.05,
242  39.05, 50.05, 68.05,
243  39.05, 50.05, 68.05,
244  39.40, 50.40, 68.40,
245  39.90, 50.90, 68.90 };
246  const Double_t HInnerRadius[NofHLayers][kNofDisks] = { 39.65, 50.65, 68.65,
247  38.50, 49.50, 67.50,
248  38.05, 49.05, 67.05,
249  38.05, 49.05, 67.05,
250  38.05, 49.05, 67.05,
251  38.05, 49.05, 67.05,
252  38.50, 49.50, 67.50,
253  38.50, 49.50, 67.50,
254  38.05, 49.05, 67.05,
255  38.05, 49.05, 67.05,
256  38.05, 49.05, 67.05,
257  38.05, 49.05, 67.05,
258  38.50, 49.50, 67.50,
259  39.65, 50.65, 68.65 };
260 
261  const Double_t HXBoxWidth = 2.30; // Using to define holding structure layers holes
262  const Double_t HXPlateWidth = 1.90;
263  const Double_t HYPlateWidth = 10.0;
264 
265  const Double_t HTZ = 0.0; // Translation parameters for vertical holes
266  const Double_t HTX = 0.0;
267  const Double_t HTY[NofHLayers][kNofDisks] = { 25.0, 45.0, 65.0,
268  25.0, 45.0, 65.0,
269  25.0, 45.0, 65.0,
270  25.0, 45.0, 65.0,
271  25.0, 45.0, 65.0,
272  25.0, 45.0, 65.0,
273  25.0, 45.0, 65.0,
274  25.0, 45.0, 65.0,
275  25.0, 45.0, 65.0,
276  25.0, 45.0, 65.0,
277  25.0, 45.0, 65.0,
278  25.0, 45.0, 60.0,
279  25.0, 45.0, 65.0,
280  25.0, 45.0, 65.0 };
281 //-------------------------------------------------------------------------------------------------------------------------------------------------
282  const Int_t kSensorStripType [2] = { 3 , 2 };
283 
284  const Double_t kSensorStripAngle[2][2] = { 0. , 0. , 0. , 0. };
285  const Double_t kSensorStripPitch[2][2] = { 0.04, 0.04, 0.04, 0.04};
286 
287  Double_t firstLayerOffset = 0.;
288 
289  for ( Int_t ilayer = 0 ; ilayer < kNofLayers ; ilayer++ ) {
290  cout << kLayerName[ilayer].Data() << " -> " << kLayerThick[ilayer] << endl;
291  firstLayerOffset += kLayerThick[ilayer];
292  }
293 
294  cout << "total thickness is " << firstLayerOffset << endl;
295  firstLayerOffset = firstLayerOffset/2.;
296  firstLayerOffset = -0.001*(TMath::Floor(1000.*firstLayerOffset));
297  cout << "first layer offset is " << firstLayerOffset << endl;
298  //-----------------------------------------------------------------------------------------------
299 
300  TString vmcWorkdir = getenv("VMCWORKDIR");
301 
302  TString outfile= "../../geometry/gem_3Stations_realistic_v1.root";
303  // TString outfile= "../../geometry/gem_3Stations.root";
304 
305  TFile* fi = new TFile(outfile,"RECREATE");
306 
307  cout << "created output file" << endl;
308  FairGeoLoader* geoLoad = new FairGeoLoader("TGeo","FairGeoLoader");
309  FairGeoInterface *geoFace = geoLoad->getGeoInterface();
310  cout << "geoface setmediafile" << endl;
311  geoFace->setMediaFile("../../geometry/media_pnd.geo");
312  cout << "geoface readmedia" << endl;
313  geoFace->readMedia();
314  cout << "geoface print" << endl;
315  geoFace->print();
316  cout << "geoface done" << endl;
317 
318  FairGeoMedia *Media = geoFace->getMedia();
319  FairGeoBuilder *geobuild=geoLoad->getGeoBuilder();
320 
321  FairGeoMedium *CbmMediumAir = Media->getMedium("air");
322  FairGeoMedium *CbmMediumPWO = Media->getMedium("PWO");
323  FairGeoMedium *CbmMediumCarbon = Media->getMedium("carbon");
324  FairGeoMedium *CbmMediumAluminium = Media->getMedium("aluminium");
325  FairGeoMedium *CbmMediumCopper = Media->getMedium("copper");
326  FairGeoMedium *CbmMediumKapton = Media->getMedium("kapton");
327  FairGeoMedium *CbmMediumArCO2 = Media->getMedium("GEMmixture");
328  FairGeoMedium *CbmMediumUranium = Media->getMedium("uranium");
329  FairGeoMedium *CbmMediumGlassFiber= Media->getMedium("GlassFiber");
330 
331  Int_t nmed=geobuild->createMedium(CbmMediumAir);
332  nmed=geobuild->createMedium(CbmMediumPWO);
333  nmed=geobuild->createMedium(CbmMediumCarbon);
334  nmed=geobuild->createMedium(CbmMediumAluminium);
335  nmed=geobuild->createMedium(CbmMediumCopper);
336  nmed=geobuild->createMedium(CbmMediumKapton);
337  nmed=geobuild->createMedium(CbmMediumArCO2);
338  nmed=geobuild->createMedium(CbmMediumUranium);
339  nmed=geobuild->createMedium(CbmMediumGlassFiber);
340 
341 
342  TGeoManager* gGeoMan = (TGeoManager*)gROOT->FindObject("FAIRGeom");
343 
344  TGeoVolume *top = new TGeoVolumeAssembly("Gem");
345 
346  gGeoMan->SetTopVolume(top);
347 
348  cout << "-------------------------------------------------------------------" << endl;
349  TList* mediaList = (TList*)gGeoMan->GetListOfMedia();
350  //TList* mediaList = (TList*)Media->getListOfMedia();
351  cout << "media known: " << mediaList->GetEntries() << "." << endl;
352  for ( Int_t itemp = 0 ; itemp < mediaList->GetEntries() ; itemp++ ) {
353  TGeoMedium* medium = (TGeoMedium*)mediaList->At(itemp);
354  cout << "medium " << itemp << " name is " << medium->GetName() << endl;
355  }
356  cout << "-------------------------------------------------------------------" << endl;
357 
358  TGeoRotation* dummyrot = new TGeoRotation();
359 
360  TGeoVolumeAssembly* SubunitVol = new TGeoVolumeAssembly("Gem_Disks");
361 
362  //----------------------------------------------------------------------------------------
363  TGeoShape *DiskShape[kNofDisks];
364  TGeoVolume *DiskVol [kNofDisks];
365  TGeoTranslation *DiskTrans[kNofDisks];
366  TGeoCombiTrans *DiskCombi[kNofDisks];
367  //--------------------------------------------------------------------------------------------
368  TGeoShape *carbonRingShape[kNofDisks];
369  TGeoVolume *carbonRingVol [kNofDisks];
370  TGeoTranslation *carbonRingTrans[kNofDisks];
371  TGeoCombiTrans *carbonRingCombi[kNofDisks];
372  //------------------------------------------------------------------------------------------
373  TGeoShape *copperRingShape[kNofDisks];
374  TGeoVolume *copperRingVol [kNofDisks];
375  TGeoTranslation *copperRingTrans[kNofDisks];
376  TGeoCombiTrans *copperRingCombi[kNofDisks];
377  //-----------------------------------------------------------------------------------------
378  TGeoShape *coverRingShape[kNofDisks];
379  TGeoVolume *coverRingVol [kNofDisks];
380  TGeoTranslation *coverRingTrans[kNofDisks];
381  TGeoCombiTrans *coverRingCombi[kNofDisks];
382  //-------------------------------------------------------------------------------------------
383  TGeoShape *tAcopperbarShape;
384  TGeoVolume *tAcopperbarVol ;
385  TGeoTranslation *tAcopperbarTrans;
386  TGeoRotation *tAcopperbarRot;
387  TGeoCombiTrans *tAcopperbarCombi;
388  //-----------------------------------------------------------------------------------------
389  TGeoShape *tBcopperbarShape;
390  TGeoVolume *tBcopperbarVol ;
391  TGeoTranslation *tBcopperbarTrans;
392  TGeoRotation *tBcopperbarRot;
393  TGeoCombiTrans *tBcopperbarCombi;
394  //-----------------------------------------------------------------------------------------
395  TGeoShape *dAcopperbarShape;
396  TGeoVolume *dAcopperbarVol ;
397  TGeoTranslation *dAcopperbarTrans;
398  TGeoRotation *dAcopperbarRot;
399  TGeoCombiTrans *dAcopperbarCombi;
400  //-----------------------------------------------------------------------------------------
401  TGeoShape *dBcopperbarShape;
402  TGeoVolume *dBcopperbarVol ;
403  TGeoTranslation *dBcopperbarTrans;
404  TGeoRotation *dBcopperbarRot;
405  TGeoCombiTrans *dBcopperbarCombi;
406  //-----------------------------------------------------------------------------------------
407  TGeoShape *AlRingShape[kNofDisks][NofSegments];
408  TGeoVolume *AlRingVol [kNofDisks][NofSegments];
409  TGeoTranslation *AlRingTrans[kNofDisks][NofSegments];
410  TGeoCombiTrans *AlRingCombi[kNofDisks][NofSegments];
411  //-----------------------------------------------------------------------------------------
412  TGeoShape *moduleRingShape[kNofDisks][NofmoduleSegments];
413  TGeoVolume *moduleRingVol [kNofDisks][NofmoduleSegments];
414  TGeoTranslation *moduleRingTrans[kNofDisks][NofmoduleSegments];
415  TGeoCombiTrans *moduleRingCombi[kNofDisks][NofmoduleSegments];
416  //-----------------------------------------------------------------------------------------
417  TGeoShape *AlumiRingShape[kNofDisks];
418  TGeoVolume *AlumiRingVol [kNofDisks];
419  TGeoTranslation *AlumiRingTrans[kNofDisks];
420  TGeoCombiTrans *AlumiRingCombi[kNofDisks];
421  //----------------------------------------------------------------------------------------------
422  TGeoShape *DiskLayersShapeA [kNofLayers][kNofDisks][4]; // 4 is number of seg
423  TGeoShape *DiskLayersShapeB [kNofLayers][kNofDisks][4];
424  TGeoSubtraction *DiskLayersSubtr [kNofLayers][kNofDisks][4];
425  TGeoShape *DiskLayersShapeC [kNofLayers][kNofDisks][4]; // final, C = A-B
426  TGeoShape *DiskLayersShapeHole [kNofLayers][kNofDisks][4];
427  TGeoCompositeShape *DiskLayersShapecompos[kNofLayers][kNofDisks][4];
428  TGeoVolume *DiskLayersVol [kNofLayers][kNofDisks][4];
429  TGeoTranslation *DiskLayersTranshA [kNofLayers][kNofDisks][4];
430  TGeoTranslation *DiskLayersTranshB [kNofLayers][kNofDisks][4];
431  TGeoTranslation *DiskLayersTrans [kNofLayers][kNofDisks][4];
432  TGeoCombiTrans *DiskLayersCombi [kNofLayers][kNofDisks][4];
433  //------------------------------------------------------------------------------------------------
434  TGeoShape *HLayersShapeTube [NofHLayers][kNofDisks];
435  TGeoShape *HLayersShapeBox [NofHLayers][kNofDisks];
436  TGeoShape *HLayersShapePlate [NofHLayers][kNofDisks];
437  TGeoShape *HLayersShapeHT [NofHLayers][kNofDisks];
438  TGeoShape *HLayersShapeHTM [NofHLayers][kNofDisks];
439  TGeoShape *HLayersShapeHTD [NofHLayers][kNofDisks];
440  TGeoCompositeShape *HLayersShapecompos[NofHLayers][kNofDisks];
441  TGeoVolume *HLayersVolcomp [NofHLayers][kNofDisks];
442  TGeoTranslation *HLayersTranstA [NofHLayers][kNofDisks];
443  TGeoTranslation *HLayersTranstB [NofHLayers][kNofDisks];
444  TGeoTranslation *HLayersTranstC [NofHLayers][kNofDisks];
445  TGeoTranslation *HLayersTranstD [NofHLayers][kNofDisks];
446  TGeoTranslation *HLayersTranstE [NofHLayers][kNofDisks];
447  TGeoTranslation *HLayersTranstF [NofHLayers][kNofDisks];
448  TGeoRotation *HLayersRotrA [NofHLayers][kNofDisks];
449  TGeoShape *HLayersShape [NofHLayers][kNofDisks];
450  TGeoVolume *HLayersVol [NofHLayers][kNofDisks];
451  TGeoTranslation *HLayersTrans [NofHLayers][kNofDisks];
452  TGeoCombiTrans *HLayersCombi [NofHLayers][kNofDisks];
453  //------------------------------------------------------------------------------------------------
454  TGeoShape *RiddleShapeTubeA ;
455  TGeoShape *RiddleShapeTubeB ;
456  TGeoShape *RiddleShapeCone ;
457  TGeoShape *RiddleShapeTubeC ;
458  TGeoShape *RiddleShapeTubeD ;
459  TGeoShape *RiddleShapeTubeE ;
460  TGeoShape *RiddleShapeTubeF ;
461  TGeoSubtraction *RiddleSubtr ;
462  TGeoCompositeShape *RiddleShapecompos;
463  TGeoVolume *RiddleVolcomp ;
464  TGeoTranslation *RiddleTrans ;
465  TGeoTranslation *RiddleTransTubeA ;
466  TGeoTranslation *RiddleTransCone ;
467  TGeoTranslation *RiddleTransTubeC ;
468  TGeoTranslation *RiddleTransTubeD ;
469  TGeoTranslation *RiddleTransTubeE ;
470  TGeoCombiTrans *RiddleCombi ;
471  TGeoTranslation *RiddleTransTubeF[4][100] ;
472  TGeoCombiTrans *RiddleCombiTranshole[4][100] ;
473  TGeoRotation *RiddleRothole[4][100] ;
474  //----------------------------------------------------------------------------------------------
475  TString outParFileName = Form("%s/macro/params/gem_3Stations_realistic_v1.digi.par",vmcWorkdir.Data());
476  //TString outParFileName = Form("%s/macro/params/gem_3Stations.digi.par",vmcWorkdir.Data());
477 
478  cout << "parameter file = \"" << outParFileName.Data() << "\"" << endl;
479 
480  ofstream pout(outParFileName.Data());
481  pout.setf(ios::fixed);
482 
483  pout << "#################################################################" << endl;
484  pout << "# Digitization parameters for GEM " << endl;
485  pout << "# with 3 Stations " << endl;
486  pout << "# Format: " << endl;
487  pout << "# Description of parameters: " << endl;
488  pout << "# [PndGemDetectors] " << endl;
489  pout << "# parameters:d station_number, ZPos, rotation_angle, number_of_sensors, \\" << endl;
490  pout << "# sensor_number, sensor_type, pos_x, pos_y, pos_z, rotAngle, inn_rad, out_rad, thick, str_ang_0, str_ang_1, barWidth, pitch_0, pitch_1, \\" << endl;
491  pout << "# sensor_number, ...." << endl;
492  pout << "# station_number, ..." << endl;
493  pout << "#################################################################" << endl;
494  pout << "[PndGemDetectors]" << setprecision(4) << endl;
495 
496  pout << "parameters:Double_t \\" << endl;
497 
498  for ( Int_t istat = 0 ; istat < kNofDisks ; istat++ ) {
499 
500  pout << " " << istat+1 << ", "
501  << setw(9) << kDiskZPosition[istat]
502  << ", 0.0, " << 2 << ", \\" << endl;
503  //-----------------------------------GEM Disk------------------------------------------------------------------------------------------------------
504  DiskShape[istat] = new TGeoTube (Form("disk%dshape",istat+1),kDiskVolInnerRadius[istat],kDiskVolOuterRadius[istat],kHalfStationThickness);
505  DiskVol [istat] = new TGeoVolume(Form("Gem_Disk%d_Volume",istat+1),DiskShape[istat],gGeoMan->GetMedium("GEMmixture"));
506  DiskTrans[istat] = new TGeoTranslation(0.,0.,kDiskZPosition[istat]);
507  cout << "station " << kDiskVolInnerRadius[istat] << " " << kDiskVolOuterRadius[istat] << " at " << kDiskZPosition[istat] << endl;
508  DiskCombi[istat] = new TGeoCombiTrans(*DiskTrans[istat],*dummyrot);
509  DiskCombi[istat]->SetName(Form("Gem_Disk%d_Volume",istat+1));
510  DiskCombi[istat]->RegisterYourself();
511  DiskVol[istat]->SetLineColor(kYellow);
512  //------------------------------------------------------------------------------------------------------------------------------------------
513  //------------------------- Gas container Ring Bottom -------------------------------------------------------------------------------------
514  carbonRingShape[istat] = new TGeoTube (Form("carbonRing%dshape",istat+1),carbonRingInnerRadius[istat],carbonRingOuterRadius[istat],carbonRingHalfThickness);
515  carbonRingVol [istat] = new TGeoVolume(Form("Gem_carbonRing%d_Volume",istat+1),carbonRingShape[istat],gGeoMan->GetMedium("carbon"));
516  carbonRingTrans[istat] = new TGeoTranslation(0.,0.,-1.5);
517  cout << "carbonRing " << carbonRingInnerRadius[istat] << " " << carbonRingOuterRadius[istat] << endl;
518  carbonRingCombi[istat] = new TGeoCombiTrans(*carbonRingTrans[istat],*dummyrot);
519  carbonRingCombi[istat]->SetName(Form("Gem_carbonRing%d_Volume",istat+1));
520  carbonRingCombi[istat]->RegisterYourself();
521  carbonRingVol[istat]->SetLineColor(kMagenta);
522  //----------------------------------------------------------------------------------------------------------------------------------------
523  DiskVol[istat]->AddNode(carbonRingVol [istat],0,carbonRingCombi[istat]);
524  //------------------------------------------------------------------------------------------------------------------------------------------
525  //------------------------- Gas container Ring top -----------------------------------------------------------------------------------------
526  copperRingShape[istat] = new TGeoTube (Form("copperRing%dshape",istat+1),copperRingInnerRadius[istat],copperRingOuterRadius[istat],copperRingHalfThickness);
527  copperRingVol [istat] = new TGeoVolume(Form("Gem_copperRing%d_Volume",istat+1),copperRingShape[istat],gGeoMan->GetMedium("copper"));
528  copperRingTrans[istat] = new TGeoTranslation(0.,0.,3.75);
529  cout << "copperRing " << copperRingInnerRadius[istat] << " " << copperRingOuterRadius[istat] << endl;
530  copperRingCombi[istat] = new TGeoCombiTrans(*copperRingTrans[istat],*dummyrot);
531  copperRingCombi[istat]->SetName(Form("Gem_copperRing%d_Volume",istat+1));
532  copperRingCombi[istat]->RegisterYourself();
533  copperRingVol[istat]->SetLineColor(kYellow-5);
534  //----------------------------------------------------------------------------------------------------------------------------------------
535  DiskVol[istat]->AddNode(copperRingVol [istat],0,copperRingCombi[istat]);
536  //----------------------------------------------------------------------------------------------------------------------------------------
537  //----------------------segments for electronics --------------------------------------------------------------------------------------------------------
538  Double_t segmentAngularSize = SegmentHalfThickness/AlRingInnerRadius[istat]*360.;
539  for ( Int_t isegm = 0 ; isegm < NofSegments ; isegm++ ) {
540  cout << " Segment " << isegm << " with name " << newsegment[isegm] << " will be put at " << FirstSegmentAngle+isegm*SegmentDeltaAngle << " (in fact from " << FirstSegmentAngle+isegm*SegmentDeltaAngle-segmentAngularSize/2. << " to " << FirstSegmentAngle+isegm*SegmentDeltaAngle+segmentAngularSize/2. << ")" << endl;
541  AlRingShape[istat][isegm] = new TGeoTubeSeg (Form("AlRing%dshape",istat+1,isegm+1),AlRingInnerRadius[istat],AlRingOuterRadius[istat],AlRingHalfThickness,FirstSegmentAngle+isegm*SegmentDeltaAngle-segmentAngularSize/2.,FirstSegmentAngle+isegm*SegmentDeltaAngle+segmentAngularSize/2. );
542  AlRingVol[istat][isegm] = new TGeoVolume (Form("Gem_AlRing%d_Volume",istat+1,isegm+1),AlRingShape[istat][isegm],gGeoMan->GetMedium("aluminium"));
543  AlRingTrans[istat][isegm] = new TGeoTranslation (0.,0.,3.75);
544  cout << "AlRing " << AlRingInnerRadius[istat] << " " << AlRingOuterRadius[istat] << endl;
545  AlRingCombi[istat][isegm] = new TGeoCombiTrans (*AlRingTrans[istat][isegm],*dummyrot);
546  AlRingCombi[istat][isegm]->SetName(Form("Gem_AlRing%d_Volume",istat+1,isegm+1));
547  AlRingCombi[istat][isegm]->RegisterYourself();
548  AlRingVol[istat][isegm]->SetLineColor(kCyan-9);
549  //-----------------------------------------------------------------------------------------------------------------------------------------
550  DiskVol[istat]->AddNode(AlRingVol[istat][isegm],0,AlRingCombi[istat][isegm]);
551  //------------------------------------------------------------------------------------------------------------------------------------------
552  }
553  //----------------------------------------------------------------------------------------------------------------------------------------
554  //---------------------- electronic module --------------------------------------------------------------------------------------------------------
555  Double_t modulesegmentAngularSize = moduleSegmentHalfThickness/moduleRingInnerRadius[istat]*360.;
556  for ( Int_t imodulesegm = 0 ; imodulesegm < NofmoduleSegments ; imodulesegm++ ) {
557  cout << " moduleSegment " << imodulesegm << " with name " << newmodulesegment[imodulesegm] << " will be put at " << FirstmoduleSegmentAngle+imodulesegm*moduleSegmentDeltaAngle << " (in fact from " << FirstmoduleSegmentAngle+imodulesegm*moduleSegmentDeltaAngle-modulesegmentAngularSize/2. << " to " << FirstmoduleSegmentAngle+imodulesegm*moduleSegmentDeltaAngle+modulesegmentAngularSize/2. << ")" << endl;
558  moduleRingShape[istat][imodulesegm] = new TGeoTubeSeg (Form("moduleRing%dshape",istat+1,imodulesegm+1),moduleRingInnerRadius[istat],moduleRingOuterRadius[istat],moduleRingHalfThickness,FirstmoduleSegmentAngle+imodulesegm*moduleSegmentDeltaAngle-modulesegmentAngularSize/2.,FirstmoduleSegmentAngle+imodulesegm*moduleSegmentDeltaAngle+modulesegmentAngularSize/2. );
559  moduleRingVol[istat][imodulesegm] = new TGeoVolume (Form("Gem_moduleRing%d_Volume",istat+1,imodulesegm+1),moduleRingShape[istat][imodulesegm],gGeoMan->GetMedium("copper"));
560  moduleRingTrans[istat][imodulesegm] = new TGeoTranslation (0.,0.,3.65);
561  cout << "moduleRing " << moduleRingInnerRadius[istat] << " " << moduleRingOuterRadius[istat] << endl;
562  moduleRingCombi[istat][imodulesegm] = new TGeoCombiTrans (*moduleRingTrans[istat][imodulesegm],*dummyrot);
563  moduleRingCombi[istat][imodulesegm]->SetName(Form("Gem_moduleRing%d_Volume",istat+1,imodulesegm+1));
564  moduleRingCombi[istat][imodulesegm]->RegisterYourself();
565  moduleRingVol[istat][imodulesegm]->SetLineColor(kGreen);
566  //-----------------------------------------------------------------------------------------------------------------------------------------
567  DiskVol[istat]->AddNode(moduleRingVol[istat][imodulesegm],0,moduleRingCombi[istat][imodulesegm]);
568  //------------------------------------------------------------------------------------------------------------------------------------------
569  }
570  //---------------------------------Cooling Ring-------------------------------------------------------------------------------------------------------------
571  AlumiRingShape[istat] = new TGeoTube (Form("AlumiRing%dshape",istat+1),AlumiRingInnerRadius[istat],AlumiRingOuterRadius[istat],AlumiRingHalfThickness);
572  AlumiRingVol [istat] = new TGeoVolume(Form("Gem_AlumiRing%d_Volume",istat+1),AlumiRingShape[istat],gGeoMan->GetMedium("aluminium"));
573  AlumiRingTrans[istat] = new TGeoTranslation(0.,0.,3.75);
574  cout << "AlumiRing " << AlumiRingInnerRadius[istat] << " " << AlumiRingOuterRadius[istat] << endl;
575  AlumiRingCombi[istat] = new TGeoCombiTrans(*AlumiRingTrans[istat],*dummyrot);
576  AlumiRingCombi[istat]->SetName(Form("Gem_AlumiRing%d_Volume",istat+1));
577  AlumiRingCombi[istat]->RegisterYourself();
578  AlumiRingVol[istat]->SetLineColor(kCyan-9);
579  //----------------------------------------------------------------------------------------------------------------------------------------
580  DiskVol[istat]->AddNode(AlumiRingVol[istat],0,AlumiRingCombi[istat]);
581  //----------------------------------------------------------------------------------------------------------------------------------------
582  //------------------------- GEM tracker cover electronic module -----------------------------------------------------------------------------------------
583  coverRingShape[istat] = new TGeoTube (Form("coverRing%dshape",istat+1),coverRingInnerRadius[istat],coverRingOuterRadius[istat],coverRingHalfThickness);
584  coverRingVol [istat] = new TGeoVolume(Form("Gem_coverRing%d_Volume",istat+1),coverRingShape[istat],gGeoMan->GetMedium("GlassFiber"));
585  coverRingTrans[istat] = new TGeoTranslation(0.,0.,7.75);
586  cout << "coverRing " << coverRingInnerRadius[istat] << " " << coverRingOuterRadius[istat] << endl;
587  coverRingCombi[istat] = new TGeoCombiTrans(*coverRingTrans[istat],*dummyrot);
588  coverRingCombi[istat]->SetName(Form("Gem_coverRing%d_Volume",istat+1));
589  coverRingCombi[istat]->RegisterYourself();
590  coverRingVol[istat]->SetLineColor(kGreen+3);
591  //----------------------------------------------------------------------------------------------------------------------------------------
592  DiskVol[istat]->AddNode(coverRingVol [istat],0,coverRingCombi[istat]);
593  //----------------------------------------------------------------------------------------------------------------------------------------
595  for ( Int_t jlay = 0 ; jlay < NofHLayers ; jlay++ ) {
596  // for ( Int_t jlay = 7 ; jlay < 8 ; jlay++ ) {
597  cout << "doing Hlayers " << jlay << endl;
598 
599  HLayersShapeTube[jlay][istat] = new TGeoTube (Form("T%dT%s",istat+1,HLayersName[jlay].Data()),HInnerRadius[jlay][istat],HOuterRadius[jlay][istat],HLayersThick[jlay]);
600  HLayersShapeHT[jlay][istat] = new TGeoTube (Form("H%dH%s",istat+1,HLayersName[jlay].Data()),0.0,4.60,HLayersThick[jlay]);
601  HLayersShapeHTM[jlay][istat] = new TGeoTube (Form("HTM%dHTM%s",istat+1,HLayersName[jlay].Data()),0.0,4.40,HLayersThick[jlay]+0.0002);
602  HLayersShapeHTD[jlay][istat] = new TGeoTube (Form("HTD%dHTD%s",istat+1,HLayersName[jlay].Data()),0.0,1.90,HLayersThick[jlay]+0.0001);
603  cout << "Tube name is " << HLayersShapeTube[jlay][istat]->GetName() << endl;
604  // cout << "TubeHTM name is " << HLayersShapeHTM[jlay][istat]->GetName() << endl;
605  HLayersShapeBox[jlay][istat] = new TGeoBBox(Form("B%dB%s",istat+1,HLayersName[jlay].Data()),HXBoxWidth[jlay][istat],HOuterRadius[jlay][istat],HLayersThick[jlay]);
606  HLayersShapePlate[jlay][istat] = new TGeoBBox(Form("P%dP%s",istat+1,HLayersName[jlay].Data()),HXPlateWidth[jlay][istat],HYPlateWidth[jlay][istat],HLayersThick[jlay]+0.0001);
607 
608  HLayersTranstA[jlay][istat] = new TGeoTranslation("tA",HTX,HTY[jlay][istat],HTZ);
609  HLayersTranstA[jlay][istat] ->RegisterYourself();
610  HLayersTranstB[jlay][istat] = new TGeoTranslation("tB",HTX,-HTY[jlay][istat],HTZ);
611  HLayersTranstB[jlay][istat] ->RegisterYourself();
612  HLayersTranstC[jlay][istat] = new TGeoTranslation("tC",HTX,HTY[jlay][istat]+9.0,HTZ);
613  HLayersTranstC[jlay][istat] ->RegisterYourself();
614  HLayersTranstD[jlay][istat] = new TGeoTranslation("tD",HTX,HTY[jlay][istat]-9.0,HTZ);
615  HLayersTranstD[jlay][istat] ->RegisterYourself();
616  HLayersTranstE[jlay][istat] = new TGeoTranslation("tE",HTX,-HTY[jlay][istat]+9.0,HTZ);
617  HLayersTranstE[jlay][istat] ->RegisterYourself();
618  HLayersTranstF[jlay][istat] = new TGeoTranslation("tF",HTX,-HTY[jlay][istat]-9.0,HTZ);
619  HLayersTranstF[jlay][istat] ->RegisterYourself();
620 
621  HLayersShapecompos[jlay][istat] = new TGeoCompositeShape(Form("compos%dcompos%s",istat+1,HLayersName[jlay].Data()),
622  Form("T%dT%s+B%dB%s+H%dH%s-HTM%dHTM%s-(P%dP%s:tA)-(HTD%dHTD%s:tC)-(HTD%dHTD%s:tD)-(P%dP%s:tB)-(HTD%dHTD%s:tE)-(HTD%dHTD%s:tF)",
623  istat+1,HLayersName[jlay].Data(),istat+1,HLayersName[jlay].Data(),istat+1,HLayersName[jlay].Data(),istat+1,HLayersName[jlay].Data(),
624  istat+1,HLayersName[jlay].Data(),istat+1,HLayersName[jlay].Data(),istat+1,HLayersName[jlay].Data(),istat+1,HLayersName[jlay].Data(),
625  istat+1,HLayersName[jlay].Data(),istat+1,HLayersName[jlay].Data()));
626 
627  cout << "composite name is " << HLayersShapecompos[jlay][istat] -> GetName() << endl;
628 
629  TString HlayersMaterial = HLayersName[jlay].Data();
630  HlayersMaterial.Remove(0,HlayersMaterial.Last('_')+1);
631  cout << "THE HMATERIAL IS \"" << HlayersMaterial.Data() << "\"" << endl;
632  HLayersVolcomp[jlay][istat] = new TGeoVolume(Form("GEMHLayersCOMP%dGEMHLayersCOMP%s",istat+1,HLayersName[jlay].Data()),HLayersShapecompos[jlay][istat],gGeoMan->GetMedium(HlayersMaterial.Data()));
633  cout << "COMP name is " << HLayersVolcomp[jlay][istat] -> GetName() << endl;
634  cout << "Hlayersmaterial = " << HlayersMaterial.Data() << endl;
635 
636  if ( HlayersMaterial.Contains("carbon" ) )
637  HLayersVolcomp[jlay][istat]->SetLineColor(kPink);
638  if ( HlayersMaterial.Contains("GlassFiber" ) )
639  HLayersVolcomp[jlay][istat]->SetLineColor(kGreen+3);
640 
641  cout << "STATION " << istat << " LAYER " << jlay << " POSITION " << HZPosition[jlay][istat] << endl;
642 
643  HLayersTrans[jlay][istat] = new TGeoTranslation(0.0,0.0,HZPosition[jlay][istat]);
644  HLayersCombi[jlay][istat] = new TGeoCombiTrans(*HLayersTrans[jlay][istat],*dummyrot);
645  HLayersCombi[jlay][istat]->SetName(Form("GEMHLayersCOMP%dGEMHLayersCOMP%s",istat+1,HLayersName[jlay].Data()));
646  HLayersCombi[jlay][istat]->RegisterYourself();
647  DiskVol[istat]->AddNode( HLayersVolcomp[jlay][istat],0,HLayersCombi[jlay][istat] );
648  }
650 
651  Double_t layerPosition = firstLayerOffset;//-kLayerThick[0]/2.;
652 
653  Int_t sensorNumber = 0;
654 
655  for ( Int_t ilay = 0 ; ilay < kNofLayers ; ilay++ ) {
656  cout << "doing layer " << ilay << endl;
657  layerPosition += kLayerThick[ilay]/2.;
658  if ( kLayerName[ilay].Contains("space") && kLayerThick[ilay] > 0.7 ) {
659  cout << "***** THE THICK SPACE LAYER IS AT " << layerPosition << endl;
660  }
661  if ( kLayerName[ilay].Contains("space") && kLayerName[ilay].Length() == 5 ) {
662  layerPosition += kLayerThick[ilay]/2.;
663  continue;
664  }
665 
666  cout << " HAHA, got layer " << kLayerName[ilay].Data() << endl;
667 
668  Double_t segPhiSpan = 360./(Double_t(kDiskNFoils[istat]));
669  Double_t segBegin = 90.;
670  cout << "will do loop over segments" << endl;
671  for ( Int_t iseg = 0 ; iseg < 1 ; iseg++ ) {
672  cout << "segment " << iseg << endl;
673  DiskLayersShapeA[ilay][istat][iseg] = new TGeoTube(Form("disk%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
674  kDiskInnerRadius[ilay][istat],kDiskOuterRadius[ilay][istat],
675  kLayerThick[ilay]/2.);
676  DiskLayersShapeB[ilay][istat][iseg] = new TGeoBBox (Form("robo%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
677  kMiddleROBarHfTh[istat]/2.,
678  kDiskOuterRadius[ilay][istat],
679  kLayerThick[ilay]);
680  DiskLayersSubtr[ilay][istat][iseg] = new TGeoSubtraction(DiskLayersShapeA[ilay][istat][iseg],
681  DiskLayersShapeB[ilay][istat][iseg]);
682 
683 
684  DiskLayersShapeC[ilay][istat][iseg] = new TGeoCompositeShape(Form("comp%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
685  DiskLayersSubtr[ilay][istat][iseg]);
686  segBegin += segPhiSpan;
687  cout << " segBegin " << segBegin << endl;
688 
689  DiskLayersShapeHole[ilay][istat][iseg] = new TGeoTube(Form("Hole%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),0.0,3.80,kLayerThick[ilay]+0.001);
690 
691  DiskLayersTranshA[ilay][istat][iseg] = new TGeoTranslation("hA",HoleTX,HoleTY,HoleTZ);
692  DiskLayersTranshA[ilay][istat][iseg] ->RegisterYourself();
693  DiskLayersTranshB[ilay][istat][iseg] = new TGeoTranslation("hB",HoleTX,-HoleTY,HoleTZ);
694  DiskLayersTranshB[ilay][istat][iseg] ->RegisterYourself();
695 
696  DiskLayersShapecompos[ilay][istat][iseg] = new TGeoCompositeShape(Form("compos%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
697  Form("comp%dseg%d%sshape-(Hole%dseg%d%sshape:hA)-(Hole%dseg%d%sshape:hB)",
698  istat+1,iseg+1,kLayerName[ilay].Data(),istat+1,iseg+1,kLayerName[ilay].Data(),istat+1,iseg+1,kLayerName[ilay].Data()));
699 
700 
701  TString layerMaterial = kLayerName[ilay].Data();
702  layerMaterial.Remove(0,layerMaterial.Last('_')+1);
703  cout << "THE MATERIAL IS \"" << layerMaterial.Data() << "\"" << endl;
704  DiskLayersVol[ilay][istat][iseg] = new TGeoVolume(Form("Gem_Disk%d_Seg%d_%s",istat+1,iseg+1,kLayerName[ilay].Data()),
705  DiskLayersShapecompos[ilay][istat][iseg],
706  gGeoMan->GetMedium(layerMaterial.Data()));
707 
708 
709  cout << "layer material = " << layerMaterial.Data() << endl;
710  if ( layerMaterial.Contains("air" ) )
711  DiskLayersVol[ilay][istat][iseg]->SetLineColor(kGray+1);
712  if ( layerMaterial.Contains("copper" ) )
713  DiskLayersVol[ilay][istat][iseg]->SetLineColor(kOrange+1);
714  if ( layerMaterial.Contains("kapton" ) )
715  DiskLayersVol[ilay][istat][iseg]->SetLineColor(kOrange+2);
716  if ( layerMaterial.Contains("aluminium" ) )
717  DiskLayersVol[ilay][istat][iseg]->SetLineColor(kCyan-9);
718  if ( layerMaterial.Contains("GEMmixture" ) )
719  DiskLayersVol[ilay][istat][iseg]->SetLineColor(kYellow);
720  if ( layerMaterial.Contains("carbon" ) )
721  DiskLayersVol[jlay][istat][iseg]->SetLineColor(kPink);
722  if ( layerMaterial.Contains("GlassFiber" ) )
723  DiskLayersVol[jlay][istat][iseg]->SetLineColor(kGreen+3);
724 
725  DiskLayersTrans[ilay][istat][iseg] = new TGeoTranslation(0.,0.,layerPosition);
726  DiskLayersCombi[ilay][istat][iseg] = new TGeoCombiTrans(*DiskLayersTrans[ilay][istat][iseg],*dummyrot);
727  DiskLayersCombi[ilay][istat][iseg]->SetName(Form("Gem_Disk%d_Seg%d_%s",istat+1,iseg+1,kLayerName[ilay].Data()));
728  DiskLayersCombi[ilay][istat][iseg]->RegisterYourself();
729  DiskVol[istat]->AddNode(DiskLayersVol[ilay][istat][iseg],0,DiskLayersCombi[ilay][istat][iseg]);
730 
731  }
732  //-------------------------------------------------------------------------------------------------------------------------------
733  cout << "volume " << kLayerName[ilay] << " from "
734  << setprecision(10) << kDiskZPosition[istat]+layerPosition-kLayerThick[ilay]/2. << " to "
735  << setprecision(10) << kDiskZPosition[istat]+layerPosition+kLayerThick[ilay]/2. << endl;
736 
737  if ( kLayerName[ilay].Contains("Gem") && kLayerName[ilay].Contains("Sensor") ) {
738  Double_t newRadius = kDiskInnerRadius[ilay][istat];
739  Double_t nofStrips = 0;
740 
741  cout << "rad = " << kDiskInnerRadius[ilay][istat] << " pitch = " << kSensorStripPitch[sensorNumber][0] << " for sensor " << sensorNumber << endl;
742  if ( kSensorStripType[sensorNumber] != 2 ) {
743  nofStrips = TMath::Ceil(2.*TMath::Pi()*kDiskInnerRadius[ilay][istat]/kSensorStripPitch[sensorNumber][0]);
744  newRadius = nofStrips*kSensorStripPitch[sensorNumber][0]/2./TMath::Pi();
745  }
746  cout << "!!!! " << istat << " " << ilay << " > there shall be " << nofStrips << " strips here so the radius should be " << newRadius << endl;
747  pout << " " << sensorNumber+1 << ", " << kSensorStripType[sensorNumber] << ", "
748  << setw(9) << 0. << ", "
749  << setw(9) << 0. << ", "
750  << setw(9) << kDiskZPosition[istat]+layerPosition << ", "
751  << setw(9) << 0. << ", "
752  << setw(9) << newRadius << ", "
753  << setw(9) << kDiskOuterRadius[ilay][istat] << ", "
754  << setw(9) << kLayerThick[ilay] << ", "
755  << setw(9) << kSensorStripAngle[sensorNumber][0] << ", "
756  // << setw(9) << kSensorStripAngle[sensorNumber][1] << ", "
757  << setw(9) << kMiddleROBarHfTh[istat] << ", "
758  << setw(9) << kSensorStripPitch[sensorNumber][0] << ", "
759  << setw(9) << kSensorStripPitch[sensorNumber][1] << ((istat==kNofDisks-1&&sensorNumber==1)?"":", \\")
760  << endl;
761  sensorNumber++;
762  }
763 
764  layerPosition += kLayerThick[ilay]/2.;
765  }
766 
767  SubunitVol->AddNode(DiskVol[istat],0,DiskCombi[istat]);
768 
769  }
770 
771  pout << "TrackFinderOnHits_ParThetaA: Double_t 59.4" << endl
772  << "TrackFinderOnHits_ParThetaB: Double_t -0.02" << endl
773  << "TrackFinderOnHits_ParTheta0: Double_t 56.1372" << endl
774  << "TrackFinderOnHits_ParTheta1: Double_t -0.000564362" << endl
775  << "TrackFinderOnHits_ParTheta2: Double_t -0.181828" << endl
776  << "TrackFinderOnHits_ParTheta3: Double_t 0.284289" << endl
777  << "TrackFinderOnHits_ParRadPhi0:Double_t 0.9944432" << endl
778  << "TrackFinderOnHits_ParRadPhi2:Double_t -0.000590706" << endl
779  << "TrackFinderOnHits_ParMat0: Double_t \\" << endl
780  // << " -2.31333e-6, 0.00067035, 0.10173" << endl
781  << " -2.32143e-6, 0.00067270, 0.10209" << endl
782  << "TrackFinderOnHits_ParMat1: Double_t \\" << endl
783  // << " -7.46844e-10, -6.6696e-7, 0.000736672" << endl
784  << " -7.49458e-10, -6.6929e-7, 0.000739251" << endl
785  << "##########################################################################################" << flush;
787  // ----------------------- Riddle shell -------------------------------------------------------------------------------------------------------------------------------
788  RiddleShapeTubeA = new TGeoTube ("TubeA" , 44.0, 47.0, 0.5 );
789  RiddleShapeTubeB = new TGeoTube ("TubeB" , 47.0, 47.1, 10.0 );
790  RiddleShapeCone = new TGeoCone ("Cone" , 20.0, 47.0, 47.1, 71.0, 71.1 );
791  RiddleShapeTubeC = new TGeoTube ("TubeC" , 71.0, 74.5, 0.5 );
792  RiddleShapeTubeD = new TGeoTube ("TubeD" , 74.5, 75.0, 10.1 );
793  RiddleShapeTubeE = new TGeoTube ("TubeE" , 71.5, 75.0, 0.5 );
794  RiddleShapeTubeF = new TGeoTube ("TubeF" , 0.0, 8.0, 2.0 );
795 
796  RiddleTransTubeA = new TGeoTranslation("trA",0.0,0.0,9.50);
797  RiddleTransTubeA ->RegisterYourself();
798  RiddleTransCone = new TGeoTranslation("trB",0.0,0.0,39.50);
799  RiddleTransCone ->RegisterYourself();
800  RiddleTransTubeC = new TGeoTranslation("trC",0.0,0.0,60.0);
801  RiddleTransTubeC ->RegisterYourself();
802  RiddleTransTubeD = new TGeoTranslation("trD",0.0,0.0,69.50);
803  RiddleTransTubeD ->RegisterYourself();
804  RiddleTransTubeE = new TGeoTranslation("trE",0.0,0.0,79.50);
805  RiddleTransTubeE ->RegisterYourself();
806 
807  //RiddleShapecompos = new TGeoCompositeShape("riddle", "TubeA+(TubeB:trA)+(Cone:trB)+(TubeC:trC)+(TubeD:trD)+(TubeE:trE)" );
808 
809  // ---------putting holes on the Riddle--------------------------------------
810 
811  const Int_t NofHoles = 16.0;
812  const Double_t rotDeltaAngle = 360.0/(Double_t(NofHoles));
813  const Int_t NofHolesRows = 4;
814 
815  Double_t holePosR[4] = {46,52,64,75}; // holes X position
816  Double_t holePosD[4] = {90,-60,-60,90}; // holes angles
817  Double_t holePosZ[4] = {10,30,50,70}; // holes Z position
818 
819  TString RiddleshapeFormula = "TubeA+(TubeB:trA)+(Cone:trB)+(TubeC:trC)+(TubeD:trD)+(TubeE:trE)" ;
820 
821 
822  for ( Int_t irow = 0 ; irow < NofHolesRows ; irow++ ) {
823  for ( Int_t ihole= 0 ; ihole < NofHoles ; ihole++ ) {
824 
825  Double_t holePosA = ihole*rotDeltaAngle;
826  Double_t holePosX = holePosR[irow]*TMath::Cos(holePosA*TMath::DegToRad());
827  Double_t holePosY = holePosR[irow]*TMath::Sin(holePosA*TMath::DegToRad());
828 
829  RiddleTransTubeF[irow][ihole] = new TGeoTranslation(Form("trF%d_%d",irow+1, ihole+1) , holePosX, holePosY,holePosZ[irow]);
830  RiddleTransTubeF[irow][ihole] ->RegisterYourself();
831 
832  RiddleRothole[irow][ihole] = new TGeoRotation(Form("ro%d_%d" , irow+1, ihole+1) ,holePosA+90, holePosD[irow],20);
833  RiddleRothole[irow][ihole] ->RegisterYourself();
834 
835  RiddleCombiTranshole[irow][ihole] = new TGeoCombiTrans( *RiddleTransTubeF[irow][ihole] , *RiddleRothole[irow][ihole] );
836 
837  RiddleCombiTranshole[irow][ihole] ->SetName(Form("TR%d_%d",irow+1, ihole+1));
838  RiddleCombiTranshole[irow][ihole] ->RegisterYourself();
839 
840  RiddleshapeFormula += Form("-(TubeF:TR%d_%d)" , irow+1,ihole+1 );
841 
842  }
843  }
844  cout<< " RiddleshapeFormula = \"" <<RiddleshapeFormula.Data()<< "\""<< endl;
845 
846  RiddleShapecompos = new TGeoCompositeShape("riddle", RiddleshapeFormula );
847 
849 
850  RiddleVolcomp = new TGeoVolume("GEMriddleCOMP", RiddleShapecompos,gGeoMan->GetMedium("carbon"));
851 
852 
853  RiddleTrans = new TGeoTranslation(0.0,0.0,117.0);
854  RiddleCombi = new TGeoCombiTrans(*RiddleTrans,*dummyrot);
855  RiddleCombi->SetName("GEMriddleCOMP");
856  RiddleCombi->RegisterYourself();
857  RiddleVolcomp->SetLineColor(kGray+1);
858  SubunitVol->AddNode(RiddleVolcomp,0,RiddleCombi);
859  //---------------------------------------------------------------------------------------------------------------------------------------------------------
860 //------------------ cables top 1 --------------------------------------------------------------------------------------------------------------------------
861  tAcopperbarShape = new TGeoBBox ("tAcopperbarshape",lcopperbarx,lcopperbary,lcopperbarHalfThickness);
862  tAcopperbarVol = new TGeoVolume("GEM_tAcopperbarVolume",tAcopperbarShape,gGeoMan->GetMedium("copper"));
863  tAcopperbarTrans = new TGeoTranslation(0.,46.1,141.);
864  tAcopperbarRot = new TGeoRotation("tAcopperbarrot",0.0,-25.0,0.0);
865  tAcopperbarRot ->RegisterYourself();
866  tAcopperbarCombi = new TGeoCombiTrans(*tAcopperbarTrans,*tAcopperbarRot);
867  tAcopperbarCombi->SetName("GEM_tAcopperbar_Volume");
868  tAcopperbarCombi->RegisterYourself();
869  tAcopperbarVol->SetLineColor(kOrange+1);
870  SubunitVol->AddNode(tAcopperbarVol,0,tAcopperbarCombi);
871  //------------------------ cables top 2 -------------------------------------------------------------------------------------------------------------------------
872  tBcopperbarShape = new TGeoBBox ("tBcopperbarshape",rcopperbarx,rcopperbary,rcopperbarHalfThickness);
873  tBcopperbarVol = new TGeoVolume("GEM_tBcopperbarVolume",tBcopperbarShape,gGeoMan->GetMedium("copper"));
874  tBcopperbarTrans = new TGeoTranslation(0.,62.5,175.);
875  tBcopperbarRot = new TGeoRotation("tBcopperbarrot",0.0,-35.0,0.0);
876  tBcopperbarRot ->RegisterYourself();
877  tBcopperbarCombi = new TGeoCombiTrans(*tBcopperbarTrans,*tBcopperbarRot);
878  tBcopperbarCombi->SetName("GEM_tBcopperbar_Volume");
879  tBcopperbarCombi->RegisterYourself();
880  tBcopperbarVol->SetLineColor(kOrange+1);
881  SubunitVol->AddNode(tBcopperbarVol,0,tBcopperbarCombi);
882  //------------------ cables down 1 --------------------------------------------------------------------------------------------------------------------------
883  dAcopperbarShape = new TGeoBBox ("dAcopperbarshape",lcopperbarx,lcopperbary,lcopperbarHalfThickness);
884  dAcopperbarVol = new TGeoVolume("GEM_dAcopperbarVolume",dAcopperbarShape,gGeoMan->GetMedium("copper"));
885  dAcopperbarTrans = new TGeoTranslation(0.,-46.1,141.);
886  dAcopperbarRot = new TGeoRotation("dAcopperbarrot",0.0,25.0,0.0);
887  dAcopperbarRot ->RegisterYourself();
888  dAcopperbarCombi = new TGeoCombiTrans(*dAcopperbarTrans,*dAcopperbarRot);
889  dAcopperbarCombi->SetName("GEM_dAcopperbar_Volume");
890  dAcopperbarCombi->RegisterYourself();
891  dAcopperbarVol->SetLineColor(kOrange+1);
892  SubunitVol->AddNode(dAcopperbarVol,0,dAcopperbarCombi);
893  //------------------------ cables down 2 -------------------------------------------------------------------------------------------------------------------------
894  dBcopperbarShape = new TGeoBBox ("dBcopperbarshape",rcopperbarx,rcopperbary,rcopperbarHalfThickness);
895  dBcopperbarVol = new TGeoVolume("GEM_dBcopperbarVolume",dBcopperbarShape,gGeoMan->GetMedium("copper"));
896  dBcopperbarTrans = new TGeoTranslation(0.,-62.5,175.);
897  dBcopperbarRot = new TGeoRotation("dBcopperbarrot",0.0,35.0,0.0);
898  dBcopperbarRot ->RegisterYourself();
899  dBcopperbarCombi = new TGeoCombiTrans(*dBcopperbarTrans,*dBcopperbarRot);
900  dBcopperbarCombi->SetName("GEM_dBcopperbar_Volume");
901  dBcopperbarCombi->RegisterYourself();
902  dBcopperbarVol->SetLineColor(kOrange+1);
903  SubunitVol->AddNode(dBcopperbarVol,0,dBcopperbarCombi);
905 
906  top->AddNode(SubunitVol,0,new TGeoCombiTrans());
907 
908  // top->CheckOverlaps(0.0001, "");
909  // gGeoManager->CheckOverlaps(0.0001,""); // [cm]
910  // gGeoManager->CheckGeometryFull();
911 
912  // TObjArray *listOfOverlaps = gGeoManager->GetListOfOverlaps();
913  // cout << "************************************************" << endl;
914  // cout<<listOfOverlaps->GetEntries()<<endl;
915  // listOfOverlaps->Print();
916  // cout << "************************************************" << endl;
917 
918  // gGeoManager->CheckOverlaps();
919  // gGeoManager->PrintOverlaps();
920 
921  gGeoMan->CloseGeometry();
922  top->Write();
923  fi->Close();
924  //gGeoManager->Export(outfile);
925 
926  //top->Raytrace();
927  top->Draw("ogl");
928  //top->Draw();
929 
930  pout.close();
931 
932  return 0;
933 }
934 
dummyrot
TGeoRotation * dummyrot
Definition: create_gem_rootgeo.C:61
geoLoad
FairGeoLoader * geoLoad
Definition: createRootGeoFile.C:16
Media
FairGeoMedia * Media
Definition: createRootFscGeometryFile.C:68
gGeoMan
TGeoManager * gGeoMan
Definition: createRootFscGeometryFile.C:85
CbmMediumCarbon
FairGeoMedium * CbmMediumCarbon
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:149
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
top
TGeoVolume * top
Definition: createRootFscGeometryFile.C:89
CbmMediumAluminium
FairGeoMedium * CbmMediumAluminium
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:150
SubunitVol
TGeoVolume * SubunitVol
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:184
geobuild
FairGeoBuilder * geobuild
Definition: createRootFscGeometryFile.C:69
fi
TFile * fi
Definition: createRootGeoFile.C:238
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
medium
TString medium
Definition: createRootGeoFile.C:83
create3StationsGem
int create3StationsGem()
Definition: create3StationsGem.C:3
CbmMediumAir
FairGeoMedium * CbmMediumAir
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:147
CbmMediumPWO
FairGeoMedium * CbmMediumPWO
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:148
geoFace
FairGeoInterface * geoFace
Definition: createRootGeoFile.C:17
SetLineColor
vDisk SetLineColor(colYellow)
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
outfile
TString outfile
Definition: createRootFscGeometryFile.C:57
Generated on Sun Aug 18 2019 03:01:08 for FairRoot/PandaRoot by   doxygen 1.8.8