A lot of changes after detector review:
[u/mrichter/AliRoot.git] / FMD / scripts / NodeGeometry.C
1 //____________________________________________________________________
2 //
3 //
4 // $Id$
5 //
6 // Script I used for rapid prototyping of the FMD3 geometry - in
7 // particular the support cone 
8 //
9 /** @defgroup node_geom Simple geometry
10     @ingroup FMD_script
11 */
12 #include <TGeometry.h>
13 #include <TNode.h>
14 #include <TXTRU.h>
15 #include <TTUBE.h>
16 #include <TTUBS.h>
17 #include <TPCON.h>
18 #include <TBRIK.h>
19 #include <TCanvas.h>
20 #include <vector>
21 #include <algorithm>
22 #include <cmath>
23 #include <iostream>
24
25 //____________________________________________________________________
26 /** @brief A 2D point
27     @ingroup node_geom
28  */
29 struct point_t
30 {
31   point_t(double x=0, double y=0) : first(x), second(y) {}
32   double first;
33   double second;
34 };
35
36 //____________________________________________________________________
37 /** @brief Shape of a ring
38     @ingroup node_geom
39  */
40 struct Ring 
41 {
42   // typedef std::pair<double,double> point_t;
43   typedef std::vector<point_t>     points_t;
44   /** Constructor 
45       @param rL         Lower radius
46       @param rH         Higer radius
47       @param theta      Opening angle
48       @param waferR     Wafer radius
49       @param siThick    Silicon thickness 
50       @param staggering Staggering of modules */
51   Ring(double rL, double rH, double theta,  double waferR, 
52        double siThick, double staggering)
53     : fStaggering(staggering),
54       fInnerRadius(rL), 
55       fOuterRadius(rH), 
56       fAngle(theta), 
57       fRadius(waferR),
58       fThickness(siThick), 
59       fVerticies(6)
60   {
61     double tan_theta  = tan(fAngle * TMath::Pi() / 180.);
62     double tan_theta2 = pow(tan_theta,2);
63     double r2         = pow(fRadius,2);
64     double ir2        = pow(fInnerRadius,2);
65     double or2        = pow(fOuterRadius,2);
66     double y_A        = tan_theta * fInnerRadius;
67     double x_D        = fInnerRadius + sqrt(r2 - tan_theta2 * ir2);
68     double x_D2       = pow(x_D,2);
69     double y_B        = sqrt(r2 - or2 + 2 * fOuterRadius * x_D - x_D2);
70     double x_C        = ((x_D + sqrt(-tan_theta2 * x_D2 
71                                      + r2 * (1 + tan_theta2))) 
72                          / (1 + tan_theta2));
73     double y_C        = tan_theta * x_C;
74     
75     fVerticies[0] = point_t(fInnerRadius,  y_A);
76     fVerticies[1] = point_t(x_C,           y_C);
77     fVerticies[2] = point_t(fOuterRadius,  y_B);
78     fVerticies[3] = point_t(fOuterRadius, -y_B);
79     fVerticies[4] = point_t(x_C,          -y_C);
80     fVerticies[5] = point_t(fInnerRadius, -y_A);
81   }
82   /** Destructor */
83   virtual ~Ring() 
84   {
85     fVerticies.clear();
86   }
87   /** Create a shape 
88       @return pointer to new shape  */
89   TShape* CreateShape()
90   {
91     std::cout << "Creating Module shape" << std::flush;
92     TXTRU* moduleShape = new TXTRU("Module","Module", "", 6, 2);
93     for (Int_t i = 0; i  < 6; i++) {
94       std::cout << "." << std::flush;
95       point_t& p = fVerticies[i];
96       moduleShape->DefineVertex(i, p.first, p.second);
97     }
98     moduleShape->DefineSection(0, -fThickness/2, 1, 0, 0);
99     moduleShape->DefineSection(1,  fThickness/2, 1, 0, 0);
100     std::cout << std::endl;
101     return (TShape*)moduleShape;
102   }
103   /** Create a node that represents a ring. 
104       @return Node */
105   TNode* CreateRing(const char* name, double z) 
106   {
107     std::cout << "Creating Ring node for " << name << std::flush;
108     double bredth = fStaggering + fThickness;
109     TShape* ringShape   = new TTUBE(Form("%sShape", name), "Ring Shape", 
110                                     "", fInnerRadius, 
111                                     fOuterRadius,bredth/2);
112     TNode*  ringNode    = new TNode(Form("%sNode", name), "Ring Node", 
113                                     ringShape, 0, 0, z+bredth/2, 0);
114     TShape* moduleShape = CreateShape();
115     Int_t n = Int_t(360 / 2 / fAngle);
116     for (Int_t i = 0; i < n; i++) {
117       std::cout << "." << std::flush;
118       ringNode->cd();
119       Double_t theta  = 2  * fAngle * i;
120       Double_t z      = -(bredth+fThickness)/2+(i%2?0:fStaggering);
121       TRotMatrix* rot = new TRotMatrix(Form("%sRotation%02d", name, i), 
122                                        "Rotation", 90, theta, 90, 
123                                        fmod(90 + theta, 360), 0, 0);
124       TNode* moduleNode = new TNode(Form("%sModule%02d", name, i), 
125                                     "Module", moduleShape, 0, 0, z,
126                                     rot);
127       moduleNode->SetFillColor(2);
128       moduleNode->SetLineColor(2);
129       moduleNode->SetLineWidth(2);
130     }
131     std::cout << std::endl;
132     ringNode->SetVisibility(0);
133     return ringNode;
134   }
135   double fStaggering;
136   /** Inner radius */
137   double fInnerRadius;
138   /** Outer radius */
139   double fOuterRadius;
140   /** Opening angle (in degrees) */
141   double fAngle;
142   /** Radius (in centimeters) */
143   double fRadius;
144   /** Thickness */
145   double fThickness;
146   /** List of verticies */
147   points_t fVerticies;
148 };
149
150 //____________________________________________________________________
151 /** @brief Shape of a detector
152     @ingroup node_geom
153  */
154 struct Detector
155 {
156   /** Constructor 
157       @param id 
158       @param inner 
159       @param outer */
160   Detector(Ring* inner, double iZ, Ring* outer=0, double oZ=0) 
161     : fInner(inner), fInnerZ(iZ), fOuter(outer), fOuterZ(oZ)
162   {}
163   /** Destructor */
164   virtual ~Detector() {}
165   /** Create rings */
166   virtual void CreateRings() 
167   {
168     if (fInner) fInner->CreateRing("inner", fInnerZ);
169     if (fOuter) fOuter->CreateRing("outer", fOuterZ);
170   }
171   /** Create a node that represents the support */
172   virtual void CreateSupport(double) { }
173   /** Pointer to inner ring */
174   Ring* fInner;
175   /** Position in z of inner ring */ 
176   double fInnerZ;
177   /** Pointer to outer ring */
178   Ring* fOuter;
179   /** Position in z of inner ring */ 
180   double fOuterZ;
181 };
182
183 //____________________________________________________________________
184 /** @brief FMD3 simple node geometry 
185     @ingroup node_geom
186  */
187 struct FMD3 : public Detector
188 {
189   /** Constructor 
190       @param inner Inner ring representation  
191       @param outer Outer ring representation */
192   FMD3(Ring* inner, Ring* outer)
193     : Detector(inner, -62.8,outer, -75.2)
194   {
195     fNoseRl     = 5.5;
196     fNoseRh     = 6.7;
197     fNoseDz     = 2.8 / 2;
198     fNoseZ      = -46;
199     fConeL      = 30.9;
200     fBackRl     = 61 / 2;
201     fBackRh     = 66.8 /2;
202     fBackDz     = 1.4 / 2;
203     fBeamDz     = .5 / 2;
204     fBeamW      = 6;
205     fFlangeR    = 49.25;
206   }
207   virtual ~FMD3() {}
208   void CreateRings()
209   {
210     double zdist      = fConeL - 2 * fBackDz - 2 * fNoseDz;
211     double tdist      = fBackRh - fNoseRh;
212     double alpha      = tdist / zdist;
213     double x, rl, rh, z;
214     z  = fNoseZ - fConeL / 2;
215     TPCON* fmd3Shape = new TPCON("fmd3Shape", "FMD 3 Shape", "", 0, 360, 7);
216     x  = fNoseZ;
217     rl = fNoseRl;
218     rh = fNoseRh;
219     fmd3Shape->DefineSection(0, x - z, rl, rh);
220     x  = fNoseZ-2*fNoseDz;
221     fmd3Shape->DefineSection(1, x - z, rl, rh);
222     x  = fInnerZ - fInner->fStaggering - fInner->fThickness; 
223     rl = fInner->fInnerRadius;
224     rh = fNoseRh + alpha * TMath::Abs(x-fNoseZ + 2 * fNoseDz);
225     fmd3Shape->DefineSection(2, x - z, rl, rh);
226     x  = fOuterZ;
227     rl = fOuter->fInnerRadius;
228     rh = fBackRh;
229     fmd3Shape->DefineSection(3, x - z, rl, rh);
230     x  = fNoseZ - zdist - 2 * fNoseDz;
231     rl = fOuter->fInnerRadius;
232     rh = fBackRh;
233     fmd3Shape->DefineSection(4, x - z, rl, rh);
234     x  = fNoseZ - zdist - 2 * fNoseDz;
235     rl = fOuter->fInnerRadius;
236     rh = fFlangeR;
237     fmd3Shape->DefineSection(5, x - z, rl, rh);
238     x  = fNoseZ - fConeL;
239     rl = fOuter->fInnerRadius;
240     rh = fFlangeR;
241     fmd3Shape->DefineSection(6, x - z, rl, rh);
242
243     TNode* fmd3Node = new TNode("fmd3Node", "FMD3 Node", fmd3Shape,0,0,z,0);
244     fmd3Node->SetLineColor(11);
245     fmd3Node->SetFillColor(11);
246     fmd3Node->SetVisibility(1);
247     fmd3Node->cd();
248     if (fInner) fInner->CreateRing("inner", fInnerZ-z);
249     fmd3Node->cd();
250     if (fOuter) fOuter->CreateRing("outer", fOuterZ-z);
251     fmd3Node->cd();
252     CreateSupport(fNoseZ - z);
253   }
254   
255   /** Create support volumes  */
256   void CreateSupport(double noseZ) 
257   {
258     TShape* noseShape = new TTUBE("noseShape", "Nose Shape", "", 
259                                   fNoseRl, fNoseRh, fNoseDz);
260     TNode*  noseNode  = new TNode("noseNode", "noseNode", noseShape, 
261                                   0, 0, noseZ - fNoseDz, 0);
262     noseNode->SetLineColor(0);
263     double  zdist = fConeL - 2 * fBackDz - 2 * fNoseDz;
264     double  tdist = fBackRh - fNoseRh;
265     double  beamL = TMath::Sqrt(zdist * zdist + tdist * tdist);
266     double  theta = -TMath::ATan2(tdist, zdist);
267     TShape* backShape = new TTUBE("backShape", "Back Shape", "", 
268                                   fBackRl, fBackRh, fBackDz);
269     TNode*  backNode  = new TNode("backNode", "backNode", backShape, 
270                                   0, 0, noseZ-2*fNoseDz-zdist-fBackDz, 0);
271     backNode->SetLineColor(0);
272     TShape* beamShape = new TBRIK("beamShape", "beamShape", "", 
273                                   fBeamDz, fBeamW / 2 , beamL / 2);
274     Int_t    n = 8;
275     Double_t r = fNoseRl + tdist / 2;
276     for (Int_t i = 0; i < n; i++) {
277       Double_t phi   = 360. / n * i;
278       Double_t t     = 180. * theta / TMath::Pi();
279       TRotMatrix* beamRotation = new TRotMatrix(Form("beamRotation%d", i), 
280                                                 Form("beamRotation%d", i),
281                                                 180-t,phi,90,90+phi,t,phi);
282       TNode* beamNode = new TNode(Form("beamNode%d", i), 
283                                   Form("beamNode%d", i), beamShape, 
284                                   r * TMath::Cos(phi / 180 * TMath::Pi()), 
285                                   r * TMath::Sin(phi / 180 * TMath::Pi()), 
286                                   noseZ-2*fNoseDz-zdist/2, beamRotation);
287       beamNode->SetLineColor(0);
288     }
289     Double_t flangel    = (fFlangeR - fBackRh) / 2;
290     TShape* flangeShape = new TBRIK("flangeShape", "FlangeShape", "", 
291                                     flangel, fBeamW / 2, fBackDz);
292     n = 4;
293     r = fBackRh + flangel;
294     for (Int_t i = 0; i < n; i++) {
295       Double_t phi = 360. / n * i + 180. / n;
296       TRotMatrix* flangeRotation = new TRotMatrix(Form("flangeRotation%d", i),
297                                                   Form("Flange Rotation %d",i),
298                                                   90,phi,90,90+phi,0,0);
299       TNode* flangeNode = new TNode(Form("flangeNode%d", i), 
300                                     Form("flangeNode%d", i), 
301                                     flangeShape,
302                                     r * TMath::Cos(phi / 180 * TMath::Pi()), 
303                                     r * TMath::Sin(phi / 180 * TMath::Pi()),
304                                     noseZ-2*fNoseDz-zdist-fBackDz, 
305                                     flangeRotation);
306       flangeNode->SetLineColor(0);
307                                   
308     }
309   }
310   /** Nose inner radius */
311   double  fNoseRl;
312   /** Nose outer radius */
313   double  fNoseRh;
314   /** Nose depth */
315   double  fNoseDz;
316   /** Nose start position */
317   double  fNoseZ;
318   /** Length of whole support structure */
319   double  fConeL;
320   /** Inner radius of back ring */
321   double  fBackRl;
322   /** Outer radius of back ring */
323   double  fBackRh;
324   /** Thickness of back ring */
325   double  fBackDz;
326   /** Thickness of beams */
327   double  fBeamDz;
328   /** Width of beams */
329   double  fBeamW;
330   /** Ending radius of flanges */
331   double  fFlangeR;
332 };
333
334 //____________________________________________________________________
335 /** @brief Create a node geometry 
336     @ingroup node_geom
337     @code 
338     .x NodeGeometry.C++
339     @endcode 
340  */
341 void
342 NodeGeometry() 
343 {
344   TGeometry* geometry = new TGeometry("geometry","geometry");
345   TShape* topShape = new TBRIK("topShape", "topShape", "", 40, 40, 150);
346   TNode* topNode = new TNode("topNode", "topNode", topShape, 0, 0, 0, 0);
347   topNode->SetVisibility(0);
348   topNode->cd();
349   
350   Ring inner( 4.3, 17.2, 18, 13.4 / 2, .03, 1);
351   Ring outer(15.6, 28.0,  9, 13.4 / 2, .03, 1);
352   FMD3 fmd3(&inner, &outer);
353   fmd3.CreateRings();
354   
355   TCanvas* c = new TCanvas("c", "c", 800, 800);
356   c->SetFillColor(1);
357   geometry->Draw();
358   // c->x3d("ogl");
359 }
360 //____________________________________________________________________
361 //
362 // EOF
363 //