FairRoot/PandaRoot
create3StationsGem_Tube.C
Go to the documentation of this file.
1 #include "iomanip.h"
2 
3 int create3StationsGem_Tube()
4 {
5  // Gem disk geometry parameters
6  //-----------------------------
7  const Int_t kNofDisks = 3;
8 
9 // const Double_t kDiskInnerRadius[kNofDisks] = { 5.0, 5.0, 5.0};
10 // const Double_t kDiskOuterRadius[kNofDisks] = { 42.0, 66.0, 90.0};
11 // const Double_t kDiskZPosition [kNofDisks] = { 120.0, 150.0, 180.0};
12 
13  const Double_t kHalfStationThickness = 3.00;
14  const Double_t kDiskInnerRadius[kNofDisks] = { 4.52, 4.52, 4.52};//Adjusted for XY strips
15  const Double_t kDiskOuterRadius[kNofDisks] = { 37.92, 46.4, 64.4};//Adjusted for XY strips
16  const Double_t kDiskZPosition [kNofDisks] = {119.6,143.1,175.6};//CAD V1312
17  const Int_t kDiskNFoils [kNofDisks] = { 2 , 2 , 4 };
18  const Double_t kMiddleROBarHfTh[kNofDisks] = { 4.3, 4.3, 4.3};
19 
20  const Int_t kNofLayers = 39;
21  const TString kLayerName [kNofLayers] = {"WindowF_kapton","WindowF_aluminium",
22  "space",
23  "CathodeF_kapton","CathodeF_aluminium",
24  "Gem1_Sensor_GEMmixture",
25  "Gem1F_copper","Gem1_kapton","Gem1B_copper",
26  "space",
27  "Gem2F_copper","Gem2_kapton","Gem2B_copper",
28  "space",
29  "Gem3F_copper","Gem3_kapton","Gem3B_copper",
30  "space",
31  "PadF_copper","Pad_kapton","Pad_copper",
32  "space",
33  "Gem4F_copper","Gem4_kapton","Gem4B_copper",
34  "space",
35  "Gem5F_copper","Gem5_kapton","Gem5B_copper",
36  "space",
37  "Gem6F_copper","Gem6_kapton","Gem6B_copper",
38  "Gem6_Sensor_GEMmixture",
39  "CathodeB_aluminium","CathodeB_kapton",
40  "space",
41  "WindowB_aluminium","WindowB_kapton"};
42  const Double_t kLayerThick[kNofLayers] = { 0.0007, 0.0001, // 2 window
43  1.0010, // +1= 3 space
44  0.0007, 0.0001, // +2= 5 cathode
45  1.0020, // +1= 6 active space
46  0.0002, 0.0050, 0.0002, // +3= 9 gemfoil
47  0.2010, // +1=10 space
48  0.0002, 0.0050, 0.0002, // +3=13 gemfoil
49  0.2010, // +1=14 space
50  0.0002, 0.0050, 0.0002, // +3=17 gemfoil
51  0.2010, // +1=18 space
52  0.0012, 0.0050, 0.0012, // +3=21 padplane
53  0.2010, // +1=22 space
54  0.0002, 0.0050, 0.0002, // +3=25 gemfoil
55  0.2010, // +1=26 space
56  0.0002, 0.0050, 0.0002, // +3=29 gemfoil
57  0.2010, // +1=30 space
58  0.0002, 0.0050, 0.0002, // +3=33 gemfoil
59  1.0020, // +1=34 active space
60  0.0001, 0.0007, // +2=36 cathode
61  1.0010, // +1=37 space
62  0.0001, 0.0007}; // +2=39 window
63 
64 
65 // type 0: r,phi strips
66 // type 1: tilted strips
67 // type 2: x,y strips
68 // type 3: r_tree,phi strips
69  const Int_t kSensorStripType [2] = { 3 , 2 };
70 
71  const Double_t kSensorStripAngle[2][2] = { 0. , 0. , 0. , 0. };
72  const Double_t kSensorStripPitch[2][2] = { 0.04, 0.04, 0.04, 0.04};
73 
74  Double_t firstLayerOffset = 0.;
75 
76  for ( Int_t ilayer = 0 ; ilayer < kNofLayers ; ilayer++ ) {
77  cout << kLayerName[ilayer].Data() << " -> " << kLayerThick[ilayer] << endl;
78  firstLayerOffset += kLayerThick[ilayer];
79  }
80 
81  cout << "total thickness is " << firstLayerOffset << endl;
82  firstLayerOffset = firstLayerOffset/2.;
83  firstLayerOffset = -0.001*(TMath::Floor(1000.*firstLayerOffset));
84  cout << "first layer offset is " << firstLayerOffset << endl;
85 
86  //--------------------------------------------------------------------
87  TString vmcWorkdir = getenv("VMCWORKDIR");
88 
89  TString outfile= "../../geometry/gem_3Stations_Tube.root";
90  TFile* fi = new TFile(outfile,"RECREATE");
91 
92  cout << "created output file" << endl;
93  FairGeoLoader* geoLoad = new FairGeoLoader("TGeo","FairGeoLoader");
94  FairGeoInterface *geoFace = geoLoad->getGeoInterface();
95  cout << "geoface setmediafile" << endl;
96  geoFace->setMediaFile("../../geometry/media_pnd.geo");
97  cout << "geoface readmedia" << endl;
98  geoFace->readMedia();
99  cout << "geoface print" << endl;
100  geoFace->print();
101  cout << "geoface done" << endl;
102 
103  FairGeoMedia *Media = geoFace->getMedia();
104  FairGeoBuilder *geobuild=geoLoad->getGeoBuilder();
105 
106  FairGeoMedium *CbmMediumAir = Media->getMedium("air");
107  FairGeoMedium *CbmMediumPWO = Media->getMedium("PWO");
108  FairGeoMedium *CbmMediumCarbon = Media->getMedium("carbon");
109  FairGeoMedium *CbmMediumAluminium = Media->getMedium("aluminium");
110  FairGeoMedium *CbmMediumCopper = Media->getMedium("copper");
111  FairGeoMedium *CbmMediumKapton = Media->getMedium("kapton");
112  FairGeoMedium *CbmMediumArCO2 = Media->getMedium("GEMmixture");
113 
114  Int_t nmed=geobuild->createMedium(CbmMediumAir);
115  nmed=geobuild->createMedium(CbmMediumPWO);
116  nmed=geobuild->createMedium(CbmMediumCarbon);
117  nmed=geobuild->createMedium(CbmMediumAluminium);
118  nmed=geobuild->createMedium(CbmMediumCopper);
119  nmed=geobuild->createMedium(CbmMediumKapton);
120  nmed=geobuild->createMedium(CbmMediumArCO2);
121 
122  TGeoManager* gGeoMan = (TGeoManager*)gROOT->FindObject("FAIRGeom");
123 
124  TGeoVolume *top = new TGeoVolumeAssembly("Gem");
125 
126  gGeoMan->SetTopVolume(top);
127 
128  cout << "-------------------------------------------------------------------" << endl;
129  TList* mediaList = (TList*)gGeoMan->GetListOfMedia();
130  //TList* mediaList = (TList*)Media->getListOfMedia();
131  cout << "media known: " << mediaList->GetEntries() << "." << endl;
132  for ( Int_t itemp = 0 ; itemp < mediaList->GetEntries() ; itemp++ ) {
133  TGeoMedium* medium = (TGeoMedium*)mediaList->At(itemp);
134  cout << "medium " << itemp << " name is " << medium->GetName() << endl;
135  }
136  cout << "-------------------------------------------------------------------" << endl;
137 
138  TGeoRotation* dummyrot = new TGeoRotation();
139 
140  TGeoVolumeAssembly* SubunitVol = new TGeoVolumeAssembly("Gem_Disks");
141 
142  TGeoShape *DiskShape[kNofDisks];
143  TGeoVolume *DiskVol [kNofDisks];
144  TGeoTranslation *DiskTrans[kNofDisks];
145  TGeoCombiTrans *DiskCombi[kNofDisks];
146 
147  TGeoShape *DiskLayersShapeA[kNofDisks][kNofLayers][4];
148  TGeoShape *DiskLayersShapeB[kNofDisks][kNofLayers][4];
149  TGeoSubtraction *DiskLayersSubtr [kNofDisks][kNofLayers][4];
150  TGeoShape *DiskLayersShapeC[kNofDisks][kNofLayers][4]; // final, C = A-B
151  TGeoVolume *DiskLayersVol [kNofDisks][kNofLayers][4];
152  TGeoTranslation *DiskLayersTrans[kNofDisks][kNofLayers][4];
153  TGeoCombiTrans *DiskLayersCombi[kNofDisks][kNofLayers][4];
154 
155  TString outParFileName = Form("%s/macro/params/gem_3Stations.digi.par",vmcWorkdir.Data());
156 
157  cout << "parameter file = \"" << outParFileName.Data() << "\"" << endl;
158 
159  ofstream pout(outParFileName.Data());
160  pout.setf(ios::fixed);
161 
162  pout << "#################################################################" << endl;
163  pout << "# Digitization parameters for GEM " << endl;
164  pout << "# with 3 Stations at "
165  << kDiskZPosition[0] << ", "
166  << kDiskZPosition[1] << ", "
167  << kDiskZPosition[2] << " cm from the target " << endl;
168  pout << "# Format: " << endl;
169  pout << "#Description of parameters: " << endl;
170  pout << "#[PndGemDetectors] " << endl;
171  pout << "#parameters:d station_number, ZPos, rotation_angle, number_of_sensors, \\" << endl;
172  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;
173  pout << "# sensor_number, ...." << endl;
174  pout << "# station_number, ..." << endl;
175  pout << "#################################################################" << endl;
176  pout << "[PndGemDetectors]" << setprecision(4) << endl;
177 
178  pout << "parameters:Double_t \\" << endl;
179 
180  for ( Int_t istat = 0 ; istat < kNofDisks ; istat++ ) {
181 
182  pout << " " << istat+1 << ", "
183  << setw(9) << kDiskZPosition[istat]
184  << ", 0.0, " << 2 << ", \\" << endl;
185 
186  DiskShape[istat] = new TGeoTube (Form("disk%dshape",istat+1),kDiskInnerRadius[istat],kDiskOuterRadius[istat],kHalfStationThickness);
187  DiskVol [istat] = new TGeoVolume(Form("Gem_Disk%d_Volume",istat+1),DiskShape[istat],gGeoMan->GetMedium("GEMmixture"));
188  DiskTrans[istat] = new TGeoTranslation(0.,0.,kDiskZPosition[istat]);
189  cout << "station " << kDiskInnerRadius[istat] << " " << kDiskOuterRadius[istat] << " at " << kDiskZPosition[istat] << endl;
190  DiskCombi[istat] = new TGeoCombiTrans(*DiskTrans[istat],*dummyrot);
191  DiskCombi[istat]->SetName(Form("Gem_Disk%d_Volume",istat+1));
192  DiskCombi[istat]->RegisterYourself();
193 
194  Double_t layerPosition = firstLayerOffset;//-kLayerThick[0]/2.;
195 
196  Int_t sensorNumber = 0;
197 
198  for ( Int_t ilay = 0 ; ilay < kNofLayers ; ilay++ ) {
199  layerPosition += kLayerThick[ilay]/2.;
200  if ( kLayerName[ilay].Contains("space") && kLayerName[ilay].Length() == 5 ) {
201  layerPosition += kLayerThick[ilay]/2.;
202  continue;
203  }
204  // cout << " HAHA, got layer " << kLayerName[ilay].Data() << endl;
205 
206  Double_t segPhiSpan = 360./(Double_t(kDiskNFoils[istat]));
207  Double_t segBegin = 90.;
208  for ( Int_t iseg = 0 ; iseg < 1 ; iseg++ ) {
209  DiskLayersShapeA[istat][ilay][iseg] = new TGeoTube(Form("disk%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
210  kDiskInnerRadius[istat],kDiskOuterRadius[istat],
211  kLayerThick[ilay]/2.);
212  DiskLayersShapeB[istat][ilay][iseg] = new TGeoBBox (Form("robo%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
213  kMiddleROBarHfTh[istat]/2.,
214  kDiskOuterRadius[istat],
215  kLayerThick[ilay]);
216  DiskLayersSubtr [istat][ilay][iseg] = new TGeoSubtraction(DiskLayersShapeA[istat][ilay][iseg],
217  DiskLayersShapeB[istat][ilay][iseg]);
218  DiskLayersShapeC[istat][ilay][iseg] = new TGeoCompositeShape(Form("comp%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
219  DiskLayersSubtr[istat][ilay][iseg]);
220  segBegin += segPhiSpan;
221 
222  TString layerMaterial = kLayerName[ilay].Data();
223  layerMaterial.Remove(0,layerMaterial.Last('_')+1);
224  // cout << "THE MATERIAL IS \"" << layerMaterial.Data() << "\"" << endl;
225  DiskLayersVol [istat][ilay][iseg] = new TGeoVolume(Form("Gem_Disk%d_Seg%d_%s",istat+1,iseg+1,kLayerName[ilay].Data()),
226  DiskLayersShapeC[istat][ilay][iseg],
227  gGeoMan->GetMedium(layerMaterial.Data()));
228 
229  // cout << "layer material = " << layerMaterial.Data() << endl;
230  if ( layerMaterial.Contains("copper" ) )
231  DiskLayersVol[istat][ilay][iseg]->SetLineColor(2);
232  if ( layerMaterial.Contains("kapton" ) )
233  DiskLayersVol[istat][ilay][iseg]->SetLineColor(3);
234  if ( layerMaterial.Contains("aluminium" ) )
235  DiskLayersVol[istat][ilay][iseg]->SetLineColor(4);
236  DiskLayersTrans[istat][ilay][iseg] = new TGeoTranslation(0.,0.,layerPosition);
237  DiskLayersCombi[istat][ilay][iseg] = new TGeoCombiTrans(*DiskLayersTrans[istat][ilay][iseg],*dummyrot);
238  DiskLayersCombi[istat][ilay][iseg]->SetName(Form("Gem_Disk%d_Seg%d_%s",istat+1,iseg+1,kLayerName[ilay].Data()));
239  DiskLayersCombi[istat][ilay][iseg]->RegisterYourself();
240  DiskVol[istat]->AddNode(DiskLayersVol [istat][ilay][iseg],0,DiskLayersCombi[istat][ilay][iseg]);
241 
242  }
243 
244  cout << "volume " << kLayerName[ilay] << " from "
245  << setprecision(10) << kDiskZPosition[istat]+layerPosition-kLayerThick[ilay]/2. << " to "
246  << setprecision(10) << kDiskZPosition[istat]+layerPosition+kLayerThick[ilay]/2. << endl;
247 
248  if ( kLayerName[ilay].Contains("Gem") && kLayerName[ilay].Contains("Sensor") ) {
249  Double_t newRadius = kDiskInnerRadius[istat];
250  Double_t nofStrips = 0;
251 
252  cout << "rad = " << kDiskInnerRadius[istat] << " pitch = " << kSensorStripPitch[sensorNumber][0] << " for sensor " << sensorNumber << endl;
253  if ( kSensorStripType[sensorNumber] != 2 ) {
254  nofStrips = TMath::Ceil(2.*TMath::Pi()*kDiskInnerRadius[istat]/kSensorStripPitch[sensorNumber][0]);
255  newRadius = nofStrips*kSensorStripPitch[sensorNumber][0]/2./TMath::Pi();
256  }
257  cout << "!!!! " << istat << " " << ilay << " > there shall be " << nofStrips << " strips here so the radius should be " << newRadius << endl;
258  pout << " " << sensorNumber+1 << ", " << kSensorStripType[sensorNumber] << ", "
259  << setw(9) << 0. << ", "
260  << setw(9) << 0. << ", "
261  << setw(9) << kDiskZPosition[istat]+layerPosition << ", "
262  << setw(9) << 0. << ", "
263  << setw(9) << newRadius << ", "
264  << setw(9) << kDiskOuterRadius[istat] << ", "
265  << setw(9) << kLayerThick[ilay] << ", "
266  << setw(9) << kSensorStripAngle[sensorNumber][0] << ", "
267  // << setw(9) << kSensorStripAngle[sensorNumber][1] << ", "
268  << setw(9) << kMiddleROBarHfTh[istat] << ", "
269  << setw(9) << kSensorStripPitch[sensorNumber][0] << ", "
270  << setw(9) << kSensorStripPitch[sensorNumber][1] << ((istat==kNofDisks-1&&sensorNumber==1)?"":", \\")
271  << endl;
272  sensorNumber++;
273  }
274 
275  layerPosition += kLayerThick[ilay]/2.;
276  }
277  SubunitVol->AddNode(DiskVol[istat],0,DiskCombi[istat]);
278  }
279  pout << "TrackFinderOnHits_ParThetaA: Double_t 59.4" << endl
280  << "TrackFinderOnHits_ParThetaB: Double_t -0.02" << endl
281  << "TrackFinderOnHits_ParTheta0: Double_t 56.1372" << endl
282  << "TrackFinderOnHits_ParTheta1: Double_t -0.000564362" << endl
283  << "TrackFinderOnHits_ParTheta2: Double_t -0.181828" << endl
284  << "TrackFinderOnHits_ParTheta3: Double_t 0.284289" << endl
285  << "TrackFinderOnHits_ParRadPhi0:Double_t 0.9944432" << endl
286  << "TrackFinderOnHits_ParRadPhi2:Double_t -0.000590706" << endl
287  << "TrackFinderOnHits_ParMat0: Double_t \\" << endl
288  << " -2.31333e-6, 0.00067035, 0.10173" << endl
289  << "TrackFinderOnHits_ParMat1: Double_t \\" << endl
290  << " -7.46844e-10, -6.6696e-7, 0.000736672" << endl
291  << "##########################################################################################" << flush;
292 
293  top->AddNode(SubunitVol,0,new TGeoCombiTrans());
294 
295  gGeoMan->CloseGeometry();
296  top->Write();
297  fi->Close();
298  // gGeoManager->Export(outfile);
299  top->Draw("ogl");
300 
301  pout.close();
302 
303  return 0;
304 }
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
create3StationsGem_Tube
int create3StationsGem_Tube()
Definition: create3StationsGem_Tube.C:3
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
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