FairRoot/PandaRoot
create4StationsGem.C
Go to the documentation of this file.
1 #include "iomanip.h"
2 
3 int create4StationsGem()
4 {
5  // Gem disk geometry parameters
6  //-----------------------------
7  const Int_t kNofDisks = 4;
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  const Double_t kDiskInnerRadius[kNofDisks] = { 2.5, 2.5, 2.5, 2.5};
13  const Double_t kDiskOuterRadius[kNofDisks] = { 45.0, 45.0, 56.0, 74.0};
14  const Double_t kDiskZPosition [kNofDisks] = { 84.0,117.0,153.0,189.0};
15  const Int_t kDiskNFoils [kNofDisks] = { 2 , 2 , 2 , 4 };
16  const Double_t kMiddleROBarHfTh[kNofDisks] = { 4.3, 4.3, 4.3, 4.3};
17 
18  const Double_t kHalfStationThickness = 3.00;
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=31 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  const Int_t kSensorStripType [2] = { 0 , 2 };
69 
70  const Double_t kSensorStripAngle[2][2] = { 0. , 0. , 0. , 0. };
71  const Double_t kSensorStripPitch[2][2] = { 0.02, 0.02, 0.02, 0.02};
72 
73  Double_t firstLayerOffset = 0.;
74 
75  for ( Int_t ilayer = 0 ; ilayer < kNofLayers ; ilayer++ ) {
76  cout << kLayerName[ilayer].Data() << " -> " << kLayerThick[ilayer] << endl;
77  firstLayerOffset += kLayerThick[ilayer];
78  }
79 
80  cout << "total thickness is " << firstLayerOffset << endl;
81  firstLayerOffset = firstLayerOffset/2.;
82  firstLayerOffset = -0.001*(TMath::Floor(1000.*firstLayerOffset));
83  cout << "first layer offset is " << firstLayerOffset << endl;
84 
85  //--------------------------------------------------------------------
86  TString vmcWorkdir = getenv("VMCWORKDIR");
87 
88  TString outfile= "../../geometry/gem_4Stations.root";
89  TFile* fi = new TFile(outfile,"RECREATE");
90 
91  cout << "created output file" << endl;
92  FairGeoLoader* geoLoad = new FairGeoLoader("TGeo","FairGeoLoader");
93  FairGeoInterface *geoFace = geoLoad->getGeoInterface();
94  cout << "geoface setmediafile" << endl;
95  geoFace->setMediaFile("../../geometry/media_pnd.geo");
96  cout << "geoface readmedia" << endl;
97  geoFace->readMedia();
98  cout << "geoface print" << endl;
99  geoFace->print();
100  cout << "geoface done" << endl;
101 
102  FairGeoMedia *Media = geoFace->getMedia();
103  FairGeoBuilder *geobuild=geoLoad->getGeoBuilder();
104 
105  FairGeoMedium *CbmMediumAir = Media->getMedium("air");
106  FairGeoMedium *CbmMediumPWO = Media->getMedium("PWO");
107  FairGeoMedium *CbmMediumCarbon = Media->getMedium("carbon");
108  FairGeoMedium *CbmMediumAluminium = Media->getMedium("aluminium");
109  FairGeoMedium *CbmMediumCopper = Media->getMedium("copper");
110  FairGeoMedium *CbmMediumKapton = Media->getMedium("kapton");
111  FairGeoMedium *CbmMediumArCO2 = Media->getMedium("GEMmixture");
112 
113  Int_t nmed=geobuild->createMedium(CbmMediumAir);
114  nmed=geobuild->createMedium(CbmMediumPWO);
115  nmed=geobuild->createMedium(CbmMediumCarbon);
116  nmed=geobuild->createMedium(CbmMediumAluminium);
117  nmed=geobuild->createMedium(CbmMediumCopper);
118  nmed=geobuild->createMedium(CbmMediumKapton);
119  nmed=geobuild->createMedium(CbmMediumArCO2);
120 
121  TGeoManager* gGeoMan = (TGeoManager*)gROOT->FindObject("FAIRGeom");
122 
123  TGeoVolume *top = new TGeoVolumeAssembly("Gem");
124 
125  gGeoMan->SetTopVolume(top);
126 
127  cout << "-------------------------------------------------------------------" << endl;
128  TList* mediaList = (TList*)gGeoMan->GetListOfMedia();
129  //TList* mediaList = (TList*)Media->getListOfMedia();
130  cout << "media known: " << mediaList->GetEntries() << "." << endl;
131  for ( Int_t itemp = 0 ; itemp < mediaList->GetEntries() ; itemp++ ) {
132  TGeoMedium* medium = (TGeoMedium*)mediaList->At(itemp);
133  cout << "medium " << itemp << " name is " << medium->GetName() << endl;
134  }
135  cout << "-------------------------------------------------------------------" << endl;
136 
137  TGeoRotation* dummyrot = new TGeoRotation();
138 
139  TGeoVolumeAssembly* SubunitVol = new TGeoVolumeAssembly("Gem_Disks");
140 
141  TGeoShape *DiskShape[kNofDisks];
142  TGeoVolume *DiskVol [kNofDisks];
143  TGeoTranslation *DiskTrans[kNofDisks];
144  TGeoCombiTrans *DiskCombi[kNofDisks];
145 
146  TGeoShape *DiskLayersShapeA[kNofDisks][kNofLayers][4];
147  TGeoShape *DiskLayersShapeB[kNofDisks][kNofLayers][4];
148  TGeoSubtraction *DiskLayersSubtr [kNofDisks][kNofLayers][4];
149  TGeoShape *DiskLayersShapeC[kNofDisks][kNofLayers][4]; // final, C = A-B
150  TGeoVolume *DiskLayersVol [kNofDisks][kNofLayers][4];
151  TGeoTranslation *DiskLayersTrans[kNofDisks][kNofLayers][4];
152  TGeoCombiTrans *DiskLayersCombi[kNofDisks][kNofLayers][4];
153 
154  TString outParFileName = Form("%s/macro/params/gem_4Stations.digi.par",vmcWorkdir.Data());
155 
156  cout << "parameter file = \"" << outParFileName.Data() << "\"" << endl;
157 
158  ofstream pout(outParFileName.Data());
159  pout.setf(ios::fixed);
160 
161  pout << "#################################################################" << endl;
162  pout << "# Digitization parameters for GEM " << endl;
163  pout << "# with 4 Stations at 90, 120, 150, 180 cm from the target " << endl;
164  pout << "# Format: " << endl;
165  pout << "#Description of parameters: " << endl;
166  pout << "#[PndGemDetectors] " << endl;
167  pout << "#parameters:d station_number, rotation_angle, number_of_sensors, \\" << endl;
168  pout << "# sensor_number, sensor_type, pos_x, pos_y, pos_z, rotAngle, inn_rad, out_rad, thick, str_ang_0, str_ang_1, pitch_0, pitch_1, \\" << endl;
169  pout << "# sensor_number, ...." << endl;
170  pout << "# station_number, ..." << endl;
171  pout << "#################################################################" << endl;
172  pout << "[PndGemDetectors]" << setprecision(4) << endl;
173 
174  pout << "parameters:Double_t \\" << endl;
175 
176  for ( Int_t istat = 0 ; istat < kNofDisks ; istat++ ) {
177 
178  pout << " " << istat+1 << ", "
179  << setw(9) << kDiskZPosition[istat]
180  << ", 0.0, " << 2 << ", \\" << endl;
181 
182  DiskShape[istat] = new TGeoTube (Form("disk%dshape",istat+1),kDiskInnerRadius[istat],kDiskOuterRadius[istat],kHalfStationThickness);
183  DiskVol [istat] = new TGeoVolume(Form("Gem_Disk%d_Volume",istat+1),DiskShape[istat],gGeoMan->GetMedium("GEMmixture"));
184  DiskTrans[istat] = new TGeoTranslation(0.,0.,kDiskZPosition[istat]);
185  cout << "station " << kDiskInnerRadius[istat] << " " << kDiskOuterRadius[istat] << " at " << kDiskZPosition[istat] << endl;
186  DiskCombi[istat] = new TGeoCombiTrans(*DiskTrans[istat],*dummyrot);
187  DiskCombi[istat]->SetName(Form("Gem_Disk%d_Volume",istat+1));
188  DiskCombi[istat]->RegisterYourself();
189 
190  Double_t layerPosition = firstLayerOffset;//-kLayerThick[0]/2.;
191 
192  Int_t sensorNumber = 0;
193 
194  for ( Int_t ilay = 0 ; ilay < kNofLayers ; ilay++ ) {
195  layerPosition += kLayerThick[ilay]/2.;
196  if ( kLayerName[ilay].Contains("space") && kLayerName[ilay].Length() == 5 ) {
197  layerPosition += kLayerThick[ilay]/2.;
198  continue;
199  }
200  // cout << " HAHA, got layer " << kLayerName[ilay].Data() << endl;
201 
202  Double_t segPhiSpan = 360./(Double_t(kDiskNFoils[istat]));
203  Double_t segBegin = 90.;
204  for ( Int_t iseg = 0 ; iseg < kDiskNFoils[istat] ; iseg++ ) {
205  DiskLayersShapeA[istat][ilay][iseg] = new TGeoTubeSeg(Form("disk%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
206  kDiskInnerRadius[istat],kDiskOuterRadius[istat],
207  kLayerThick[ilay]/2.,
208  segBegin,segBegin+segPhiSpan);
209  DiskLayersShapeB[istat][ilay][iseg] = new TGeoBBox (Form("robo%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
210  kMiddleROBarHfTh[istat]/2.,
211  kDiskOuterRadius[istat],
212  kLayerThick[ilay]);
213  DiskLayersSubtr [istat][ilay][iseg] = new TGeoSubtraction(DiskLayersShapeA[istat][ilay][iseg],
214  DiskLayersShapeB[istat][ilay][iseg]);
215  DiskLayersShapeC[istat][ilay][iseg] = new TGeoCompositeShape(Form("comp%dseg%d%sshape",istat+1,iseg+1,kLayerName[ilay].Data()),
216  DiskLayersSubtr[istat][ilay][iseg]);
217  segBegin += segPhiSpan;
218 
219  TString layerMaterial = kLayerName[ilay].Data();
220  layerMaterial.Remove(0,layerMaterial.Last('_')+1);
221  // cout << "THE MATERIAL IS \"" << layerMaterial.Data() << "\"" << endl;
222  DiskLayersVol [istat][ilay][iseg] = new TGeoVolume(Form("Gem_Disk%d_Seg%d_%s",istat+1,iseg+1,kLayerName[ilay].Data()),
223  DiskLayersShapeC[istat][ilay][iseg],
224  gGeoMan->GetMedium(layerMaterial.Data()));
225 
226  // cout << "layer material = " << layerMaterial.Data() << endl;
227  if ( layerMaterial.Contains("copper" ) )
228  DiskLayersVol[istat][ilay][iseg]->SetLineColor(2);
229  if ( layerMaterial.Contains("kapton" ) )
230  DiskLayersVol[istat][ilay][iseg]->SetLineColor(3);
231  if ( layerMaterial.Contains("aluminium" ) )
232  DiskLayersVol[istat][ilay][iseg]->SetLineColor(4);
233  DiskLayersTrans[istat][ilay][iseg] = new TGeoTranslation(0.,0.,layerPosition);
234  DiskLayersCombi[istat][ilay][iseg] = new TGeoCombiTrans(*DiskLayersTrans[istat][ilay][iseg],*dummyrot);
235  DiskLayersCombi[istat][ilay][iseg]->SetName(Form("Gem_Disk%d_Seg%d_%s",istat+1,iseg+1,kLayerName[ilay].Data()));
236  DiskLayersCombi[istat][ilay][iseg]->RegisterYourself();
237  DiskVol[istat]->AddNode(DiskLayersVol [istat][ilay][iseg],0,DiskLayersCombi[istat][ilay][iseg]);
238 
239  }
240 
241  cout << "volume " << kLayerName[ilay] << " from "
242  << setprecision(10) << kDiskZPosition[istat]+layerPosition-kLayerThick[ilay]/2. << " to "
243  << setprecision(10) << kDiskZPosition[istat]+layerPosition+kLayerThick[ilay]/2. << endl;
244 
245  if ( kLayerName[ilay].Contains("Gem") && kLayerName[ilay].Contains("Sensor") ) {
246  Double_t newRadius = kDiskInnerRadius[istat];
247  Double_t nofStrips = 0;
248 
249  cout << "rad = " << kDiskInnerRadius[istat] << " pitch = " << kSensorStripPitch[sensorNumber][0] << " for sensor " << sensorNumber << endl;
250  if ( kSensorStripType[sensorNumber] != 2 ) {
251  nofStrips = TMath::Ceil(2.*TMath::Pi()*kDiskInnerRadius[istat]/kSensorStripPitch[sensorNumber][0]);
252  newRadius = nofStrips*kSensorStripPitch[sensorNumber][0]/2./TMath::Pi();
253  }
254  cout << "!!!! " << istat << " " << ilay << " > there shall be " << nofStrips << " strips here so the radius should be " << newRadius << endl;
255  pout << " " << sensorNumber+1 << ", " << kSensorStripType[sensorNumber] << ", "
256  << setw(9) << 0. << ", "
257  << setw(9) << 0. << ", "
258  << setw(9) << kDiskZPosition[istat]+layerPosition << ", "
259  << setw(9) << 0. << ", "
260  << setw(9) << newRadius << ", "
261  << setw(9) << kDiskOuterRadius[istat] << ", "
262  << setw(9) << kLayerThick[ilay] << ", "
263  << setw(9) << kSensorStripAngle[sensorNumber][0] << ", "
264  // << setw(9) << kSensorStripAngle[sensorNumber][1] << ", "
265  << setw(9) << kMiddleROBarHfTh[istat] << ", "
266  << setw(9) << kSensorStripPitch[sensorNumber][0] << ", "
267  << setw(9) << kSensorStripPitch[sensorNumber][1] << ((istat==kNofDisks-1&&sensorNumber==1)?"":", \\")
268  << endl;
269  sensorNumber++;
270  }
271 
272  layerPosition += kLayerThick[ilay]/2.;
273  }
274  SubunitVol->AddNode(DiskVol[istat],0,DiskCombi[istat]);
275  }
276  pout << "TrackFinderOnHits_ParThetaA: Double_t 59.4" << endl
277  << "TrackFinderOnHits_ParThetaB: Double_t -0.02" << endl
278  << "TrackFinderOnHits_ParTheta0: Double_t 56.1372" << endl
279  << "TrackFinderOnHits_ParTheta1: Double_t -0.000564362" << endl
280  << "TrackFinderOnHits_ParTheta2: Double_t -0.181828" << endl
281  << "TrackFinderOnHits_ParTheta3: Double_t 0.284289" << endl
282  << "TrackFinderOnHits_ParRadPhi0:Double_t 0.9944432" << endl
283  << "TrackFinderOnHits_ParRadPhi2:Double_t -0.000590706" << endl
284  << "TrackFinderOnHits_ParMat0: Double_t \\" << endl
285  << " -2.31333e-6, 0.00067035, 0.10173" << endl
286  << "TrackFinderOnHits_ParMat1: Double_t \\" << endl
287  << " -7.46844e-10, -6.6696e-7, 0.000736672" << endl
288  << "##########################################################################################" << flush;
289 
290  top->AddNode(SubunitVol,0,new TGeoCombiTrans());
291 
292  gGeoMan->CloseGeometry();
293  top->Write();
294  fi->Close();
295  // gGeoManager->Export(outfile);
296  top->Draw("ogl");
297 
298  pout.close();
299 
300  return 0;
301 }
dummyrot
TGeoRotation * dummyrot
Definition: create_gem_rootgeo.C:61
create4StationsGem
int create4StationsGem()
Definition: create4StationsGem.C:3
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
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 Wed Apr 3 2019 10:02:45 for FairRoot/PandaRoot by   doxygen 1.8.5