FairRoot/PandaRoot
create_HPGe_rootgeo.C
Go to the documentation of this file.
1 void HPGe_geo(){
2  // Forward tof geometry parameters
3  //-----------------------------
4  //-- macro created by Alicia S. based on
5  //-- Panda TPR and on previous geant3 geometry for forward tof(Vladimir V.)
6  const Double_t kBCentX = 2.5; // half length(cm) //from EMC TDR
7  const Double_t kBCentY = 70.0; //half length (cm) //?
8  const Double_t kBCentZ = 0.75; //half length (cm) //?
9 
10  const Double_t kBBeamX = 2.5; // half length(cm) //from EMC TDR
11  const Double_t kBBeamY = 32.5; //alf length (cm) //?
12  const Double_t kBBeamZ = 0.75; //half length (cm) //?
13 
14  const Double_t kBVertX = 5.0; //half length (cm) //from EMC TDR
15  const Double_t kBVertY = 70.0; //half length (cm) //?
16  const Double_t kBVertZ = 0.75; //half length (cm) //from EMC TDR
17 
18  //--------------------------------------------------------------------
19  int fcrl; // each cluster consists of several crystals
20  int fclr;
21 
22  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
23 
24  // Load this libraries
25  gSystem->Load("libGeoBase");
26  gSystem->Load("libParBase");
27  gSystem->Load("libBase");
28  gSystem->Load("libPndData");
29  gSystem->Load("libPassive");
30 
31  //TString outfile= "../../geometry/HPGeCluster.root";
32  TString outfile= "./HPGeCluster.root";
33  TFile* fi = new TFile(outfile,"RECREATE");
34 
35  FairGeoLoader* geoLoad = new FairGeoLoader("TGeo","FairGeoLoader");
36  FairGeoInterface *geoFace = geoLoad->getGeoInterface();
37  geoFace->setMediaFile("../../geometry/media_pnd.geo");
38  geoFace->readMedia();
39  geoFace->print();
40 
41  FairGeoMedia *Media = geoFace->getMedia();
42  FairGeoBuilder *geobuild=geoLoad->getGeoBuilder();
43 
44  FairGeoMedium *medGe = Media->getMedium("germanium");
45  Int_t nmedGe=geobuild->createMedium(medGe);
46 
47  //Scintillator plate
48  FairGeoMedium *medsci = Media->getMedium("polypropylene");
49  Int_t nmedsci=geobuild->createMedium(medsci);
50 
51  // aluminum capsule
52  FairGeoMedium *medcap = Media->getMedium("HYPaluminium");
53  Int_t nmedcap=geobuild->createMedium(medcap);
54 
55 
56 
57  TGeoManager* gGeoMan = (TGeoManager*)gROOT->FindObject("FAIRGeom");
58 
59  TGeoVolume *top = new TGeoVolumeAssembly("hpge");
60 
61  gGeoMan->SetTopVolume(top);
62  Double_t deg = TMath::Pi()/180.;
63 
64 
65  TGeoVolume *sphere = new TGeoVolumeAssembly("sphere");
66 
67  TGeoTube* Crystal_Tube = new TGeoTube("tube",0.,7.0/2.,7.8/2.);
68 
69  TGeoPgon* Crystal_Polyhedra = new TGeoPgon("Pgon",0., 360., 6,2);
70 
71  Crystal_Polyhedra->DefineSection(0,-3.9,0.,3.50*cos(30.*deg));
72  Crystal_Polyhedra->DefineSection(1,3.9,0.,(3.50));
73 
74  // scintillator plate
75  TGeoPgon* lay_sci = new TGeoPgon("PgLay",0., 360., 6,2);
76  //first thickness 0.5 cm,
77  //second 1.cm
78  lay_sci->DefineSection(0,-0.25,0.,3.50*cos(30.*deg));
79  lay_sci->DefineSection(1,0.25,0.,(3.50)*cos(30.*deg));
80  //TGeoPgon* lay_sci ;
81 
82  //aluminum capsule
83  TGeoPgon* lay_capl = new TGeoPgon("PgcapL",0., 360., 6,2);
84  //first thickness 0.5 cm,
85  //second 1.cm
86  lay_capl->DefineSection(0,-4.2,0.,3.144);
87  lay_capl->DefineSection(1,4.2,0.,(3.627));
88 
89  TGeoPgon* lay_caps = new TGeoPgon("PgcapS",0., 360., 6,2);
90  //first thickness 0.5 cm,
91  //second 1.cm
92  lay_caps->DefineSection(0,-3.95,0.,3.078);
93  lay_caps->DefineSection(1,3.95,0.,(3.553));
94 
95 
96  //###real geometry
97  //###TGeoCompositeShape* logicCrystal_test = new TGeoCompositeShape("test","Pgon*tube");
98 
99  TGeoPgon* logicCrystal_test;
100  logicCrystal_test = Crystal_Polyhedra ;
101  TGeoCompositeShape* lay_cap = new TGeoCompositeShape("test","(PgcapL)-(PgcapS)");
102 
103 
104 
105  // GeCluster* HypGe= new GeCluster();
106 
107 
108  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
109  30.0, 168.0, 90., 2.0,1000,fcrl,fclr); // L
110  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
111  30.0, 168.0,18, 2.0,1200,fcrl,fclr); // N
112 
113  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
114  30.1, 168.0, 162, 2.0,400,fcrl,fclr); // D
115 
116  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
117  30.3, 168.0, 234, 2.0,600,fcrl,fclr); // F
118 
119  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
120  30., 168.0, 306., 2.0,1400,fcrl,fclr); // Q
121 
122  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,32.2,
123  148.0, 54.0, 2.0,1100,fcrl,fclr); // M
124  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,32.9,
125  148.0, 126.0, 2.0,300,fcrl,fclr); // C
126  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,31.4,
127  148.0, 198.0, 2.0,500,fcrl,fclr); // E
128  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,32.6,
129  148.0, 270.0, 2.0,1500,fcrl,fclr); // R
130  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,30.3,
131  148.0, 342.0, 2.0,1300,fcrl,fclr); // P
132  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
133  30.0, 140.0, 90., 2.0,700,fcrl,fclr); // G
134  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
135  30.2, 140.0,18., 2.0,800,fcrl,fclr); // J
136 
137  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
138  30.15, 140.0, 162., 2.0,100,fcrl,fclr); // A
139 
140  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
141  30.15, 140.0, 234., 2.0,200,fcrl,fclr); // B
142 
143  CreateCluster(logicCrystal_test,lay_sci,lay_cap,*sphere,
144  30.2, 140.0, 306., 2.0,900,fcrl,fclr); // K
145 
146 
147 
148  top->AddNode(sphere, 0,new TGeoCombiTrans(0., 0., -78., new TGeoRotation (0)));
149 
150 
151  gGeoMan->CloseGeometry();
152  top->Write();
153  fi->Close();
154 // gGeoManager->Export(outfile);
155  top->Draw("ogl");
156 }
157 
158 void CreateCluster(TGeoPgon* logicCrystal_test,
159  TGeoPgon* lay_sci,
160  TGeoCompositeShape* lay_cap,
161  TGeoVolume& sphere,
162  Double_t distance,
163  Double_t theta,
164  Double_t phi,
165  Double_t depth_first_interaction,
166  int id_start, int fnum_crystal,int fnum_cluster)
167 {
168 char cluster_name[10];
169  if(id_start==800) sprintf(cluster_name,"J");
170  if(id_start==700) sprintf(cluster_name,"G");
171  if(id_start==100) sprintf(cluster_name,"A");
172  if(id_start==200) sprintf(cluster_name,"B");
173  if(id_start==900) sprintf(cluster_name,"K");
174 
175  if(id_start==1100) sprintf(cluster_name,"M");
176  if(id_start==300) sprintf(cluster_name,"C");
177  if(id_start==500) sprintf(cluster_name,"E");
178  if(id_start==1500) sprintf(cluster_name,"R");
179  if(id_start==1300) sprintf(cluster_name,"P");
180 
181  if(id_start== 400) sprintf(cluster_name,"D");
182  if(id_start== 600) sprintf(cluster_name,"F");
183  if(id_start==1000) sprintf(cluster_name,"L");
184  if(id_start==1200) sprintf(cluster_name,"N");
185 
186  if(id_start==1400) sprintf(cluster_name,"Q");
187 
188  Double_t deg = TMath::Pi()/180.;
189  Double_t phi1;
190 
191  fnum_crystal = 0; // each cluster consists of several crystals
192  fnum_cluster = 0; // cluster number
193 
194  distance = distance + 2.5;
195 
196  phi1=phi*deg;
197 
198  if(phi1 >360.0*deg) phi1 = phi1 - 360.0*deg;
199  if(phi1>=0.0 && phi1 <=180.0*deg) phi1 = 180.0*deg - phi1;
200  if(phi1>180.0*deg && phi1 <360.0*deg) phi1 = 540.0*deg - phi1;
201  phi1 = phi1+90.0*deg;
202 
203  if(phi >360.0) phi = phi - 360.0;
204  if(phi>=0.0 && phi <=180.0) phi = 180.0 - phi;
205  if(phi>180.0 && phi <360.0) phi = 540.0 - phi;
206  phi = phi+90.0;
207 
208  TGeoRotation* rot3D = new TGeoRotation();
209 
210  TVector3* rotatedPos;
211  //Double_t rotatedPos* = new Double_t();
212 
213  TGeoVolume* crystal_centre[60];
214  TGeoVolume* laySci[60];
215  TGeoVolume* Cap[60];
216 
217  //TGeoCombiTrans *c1;
218 
219  TGeoRotation* rot3D2 = new TGeoRotation();
220  TVector3* rotatedPos2;
221  //Double_t rotatedPos2* = new Double_t();
222  //TGeoCombiTrans *c2;
223  rot3D2->SetAngles(0, 0, 0);
224  rot3D2->RotateX(theta);
225  rot3D2->RotateZ(phi);
226 
227 
228  rot3D->SetAngles(0, 0, 0);
229  rot3D->RotateX(theta);
230  rot3D->RotateZ(phi);
231 
232 
233  rotatedPos = new TVector3(rotate(TVector3(0.,0.,62.22-8.05/2.0),rot3D)
234  - rotate(TVector3(0.,0.,62.22-8.05),rot3D2)
235  + rotate(TVector3(0.,0.,distance),rot3D2));
236 
237 
238 
239  rotatedPos2 = new TVector3(rotate(TVector3(0.,0.,62.22-8.05/2.0),rot3D)
240  - rotate(TVector3(0.,0.,62.22-8.05),rot3D2)
241  + rotate(TVector3(0.,0.,distance-5.0),rot3D2));
242 
243  TGeoMedium *Ge = gGeoManager->GetMedium("germanium");
244  //scintillator plate
245  TGeoMedium *sci = gGeoManager->GetMedium("polypropylene");
246  //capsule
247  TGeoMedium *cap = gGeoManager->GetMedium("HYPaluminium");
248 
249  std::cout<<"material "<<Ge<<" "<<sci<<" "<<cap<<std::endl;
250 
251 
252  char name_crystal[13];
253  sprintf (name_crystal,"GeCrystal%02d",id_start+6);
254 
255 
256  char name_lay[13];
257  char name_cap[13];
258  //scint.
259  sprintf (name_lay,"laySci%02d",id_start+6);
260  //capsula
261  sprintf (name_cap,"Cap%02d",id_start+6);
262 
263 
264  crystal_centre[fnum_crystal] = new TGeoVolume(name_crystal,
265  (const TGeoShape *)logicCrystal_test,
266  (const TGeoMedium *)Ge );
267 
268  laySci[fnum_crystal] = new TGeoVolume(name_lay,
269  (const TGeoShape *)lay_sci,
270  (const TGeoMedium *)sci );
271 
272  Cap[fnum_crystal] = new TGeoVolume(name_cap,(const TGeoShape *)lay_cap,
273  (const TGeoMedium *)cap );
274 
275  //DataG4 data = read.GetData(id_start+6);
276 
277 
278  // sphere.AddNode(crystal_centre[fnum_crystal],id_start+6,
279  // new TGeoCombiTrans(rotatedPos.X(),
280  // rotatedPos.Y(),
281  // rotatedPos.Z(),
282  // new TGeoRotation (rot3D)));
283 
284  sphere.AddNode(crystal_centre[fnum_crystal],id_start+6,
285  new TGeoCombiTrans(rotatedPos->x(),
286  rotatedPos->y(),
287  rotatedPos->z(),
288  new TGeoRotation (*rot3D)));
289 
290  sphere.AddNode(Cap[fnum_crystal],id_start+6,
291  new TGeoCombiTrans(rotatedPos->x(),
292  rotatedPos->y(),
293  rotatedPos->z(),
294  new TGeoRotation (*rot3D)));
295 
296  //scint
297  //DataG4 data2 = readsci.GetData(id_start+6);
299 
300 
301  sphere.AddNode(laySci[fnum_crystal],id_start+6,
302  new TGeoCombiTrans(rotatedPos2->x(),
303  rotatedPos2->y(),
304  rotatedPos2->z(),
305  new TGeoRotation (*rot3D)));
306 
307 
308 
309  //------------------------------------------------
310  // Sensitive detectors
311  //------------------------------------------------
312 
313 
314 
315 
316 
317  cout << id_start+6 << " " <<rotatedPos->x()<<" "
318  <<rotatedPos->y()<<" "
319  <<rotatedPos->z()<<endl;
320 
321 
322 
323  //------------------------------------------------
324  // Crystals surrounding the central one.
325  //------------------------------------------------
326 
327  fnum_crystal = fnum_crystal + 1;
328 
329  //Float_t v1[6],v2[6],v3[6];
330 
331 
332  char tempname[100];
333  char periname[13];
334  char cappname[13];
335  TGeoVolume* crystal[60];
336  TGeoVolume* laySci_peri[60];
337  TGeoVolume* Cap_peri[60];
338 
339  TGeoRotation* rot3D3 = new TGeoRotation();
340  TGeoRotation* rot3D4 = new TGeoRotation();
341 
342  TVector3* v1;
343  TVector3* v2;
344 
345 
346 
347  for(int iii=0;iii<=1;iii++)
348  {
349 
350 
351  rot3D3->SetAngles(0., 0., 0.);
352  rot3D3->RotateX(-6.66);
353  rot3D3->RotateZ((-1*iii+1)*60.0);
354  rot3D3->RotateX(theta);
355  rot3D3->RotateZ(phi);
356 
357  // TVector3 v7;
358  // v7 = ((((TVector3(0.,0.,62.22-8.05/2.0)).RotateX(-6.66)).RotateZ((-1*iii+1)*60.0)).RotateX(theta)).RotateZ(phi)-((TVector3(0.,0.,62.22-8.05)).RotateX(theta)).RotateZ(phi)+((TVector3(0.,0.,distance)).RotateX(theta)).RotateZ(phi);
359 
360  v1 = new TVector3(rotate(TVector3(0.,0.,62.22-8.05/2.0),rot3D3)
361  - rotate(TVector3(0.,0.,62.22-8.05),rot3D2)
362  + rotate(TVector3(0.,0.,distance),rot3D2));
363 
364  //cout<<" crystal v1 "<<iii<<endl;
365  v1->Print();
366 
367  v2 = new TVector3(rotate(TVector3(0.,0.,62.22-8.05/2.0),rot3D3)
368  - rotate(TVector3(0.,0.,62.22-8.05),rot3D2)
369  + rotate(TVector3(0.,0.,distance-5.0),rot3D2));
370 
371  //cout<<" crystal v2 "<<iii<<endl;
372  v2->Print();
373 
374  sprintf(tempname,"GeCrystal%04d",id_start+iii);
375  //TGeoVolume *crystal[fnum_crystal];
376 
377  //scint
378  sprintf (periname,"laySci%04d",id_start+iii);
379  //capsula
380  sprintf (cappname,"Cap%04d",id_start+iii);
381 
382  crystal[fnum_crystal] = new TGeoVolume(tempname,
383  (const TGeoShape *)logicCrystal_test,
384  (const TGeoMedium *)Ge);
385 
386  laySci_peri[fnum_crystal] = new TGeoVolume(periname,
387  (const TGeoShape *)lay_sci,
388  (const TGeoMedium *)sci);
389 
390  Cap_peri[fnum_crystal] = new TGeoVolume(cappname,
391  (const TGeoShape *)lay_cap,
392  (const TGeoMedium *)cap);
393 
394 
395 
396  //DataG4 data = read.GetData(id_start+iii);
397 
398  // sphere.AddNode(crystal[fnum_crystal],id_start+iii,
399  // new TGeoCombiTrans(PosOrb.X(),PosOrb.Y(),
400  // PosOrb.Z(), new TGeoRotation (rot3D)));
401 
402  sphere.AddNode(crystal[fnum_crystal],id_start+iii,
403  new TGeoCombiTrans(v1->x(),v1->y(),v1->z(),
404  new TGeoRotation (*rot3D3)));
405 
406  //capsula
407  sphere.AddNode(Cap_peri[fnum_crystal],id_start+iii,
408  new TGeoCombiTrans(v1->x(),v1->y(),v1->z(),
409  new TGeoRotation (*rot3D3)));
410 
411  //DataG4 data2 = readsci.GetData(id_start+iii);
412 
413 
414 
415  sphere.AddNode(laySci_peri[fnum_crystal],id_start+iii,
416  new TGeoCombiTrans(v2->x(),v2->y(),v2->z(),
417  new TGeoRotation (*rot3D)));
418 
419 
420  //------------------------------------------------
421  // Sensitive detectors
422  //------------------------------------------------
423 
424  //AddSensitiveVolume(crystal[fnum_crystal]);
425 
426  // cout << id_start+iii << " " << v1[0] << " "<< v1[2]<< " "<< v1[2]<< std::endl;
427  cout << id_start+iii << " " <<v1->x()<< " " <<v1->y()<< " " <<v1->z()<<endl;
428 
429  fnum_crystal = fnum_crystal + 1;
430 
431  }
432 
433  fnum_cluster = fnum_cluster + 1;
434 
435 
436 
437 
438 }
439 
440 TVector3 rotate(TVector3 vec, TGeoRotation* rotma)
441 {
442  TVector3 v;
443  v.SetXYZ(0,0,0);
444 
445 
446  const Double_t* b;
447  b = rotma->GetRotationMatrix();
448  //a = new TRotation(b[0],b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8]);
449 
450  for (Int_t i=0; i<9; i++) {
451  if (TMath::Abs(b[i])<1E-10) b[i] = 0;
452  if (TMath::Abs(b[i]-1)<1E-10) b[i] = 1;
453  if (TMath::Abs(b[i]+1)<1E-10) b[i] = -1;
454  }
455  for(i=0;i<9;i++){
456  //cout<<" matrix "<<b[i]<<endl;
457  }
458  v.SetX(b[0]*vec.X()+b[1]*vec.Y()+b[2]*vec.Z());
459  v.SetY(b[3]*vec.X()+b[4]*vec.Y()+b[5]*vec.Z());
460  v.SetZ(b[6]*vec.X()+b[7]*vec.Y()+b[8]*vec.Z());
461  //cout<<" rotate "<<v.X()<<" "<<v.Y()<<" "<<v.Z()<<endl;
462 
463  return v;
464 
465 
466 }
467 
468 
theta
double theta
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
geoLoad
FairGeoLoader * geoLoad
Definition: createRootGeoFile.C:16
i
Int_t i
Definition: run_full.C:25
Media
FairGeoMedia * Media
Definition: createRootFscGeometryFile.C:68
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
gGeoMan
TGeoManager * gGeoMan
Definition: createRootFscGeometryFile.C:85
HPGe_geo
void HPGe_geo()
Definition: create_HPGe_rootgeo.C:1
gGeoManager
TGeoManager * gGeoManager
Definition: createRootGeoFile.C:97
v
__m128 v
Definition: P4_F32vec4.h:4
top
TGeoVolume * top
Definition: createRootFscGeometryFile.C:89
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
CreateCluster
void CreateCluster(TGeoPgon *logicCrystal_test, TGeoPgon *lay_sci, TGeoCompositeShape *lay_cap, TGeoVolume &sphere, Double_t distance, Double_t theta, Double_t phi, Double_t depth_first_interaction, int id_start, int fnum_crystal, int fnum_cluster)
Definition: create_HPGe_rootgeo.C:158
geobuild
FairGeoBuilder * geobuild
Definition: createRootFscGeometryFile.C:69
fi
TFile * fi
Definition: createRootGeoFile.C:238
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
rotate
TVector3 rotate(TVector3 vec, TGeoRotation *rotma)
Definition: create_HypSTxy3Cmic_RootGeo.C:37
phi
double phi
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
v1
TVector3 v1
Definition: bump_analys.C:40
v2
TVector3 v2
Definition: bump_analys.C:40
geoFace
FairGeoInterface * geoFace
Definition: createRootGeoFile.C:17
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
vec
dble_vec_t vec[12]
Definition: ranlxd.cxx:380
outfile
TString outfile
Definition: createRootFscGeometryFile.C:57
Generated on Wed Apr 3 2019 08:43:49 for FairRoot/PandaRoot by   doxygen 1.8.5