Addes script to compare Naiive, Poisson to Hits, Primaries
[u/mrichter/AliRoot.git] / FMD / scripts / SimpleGeometry.C
1 //____________________________________________________________________
2 //
3 // Script I used for rapid prototyping of the FMD3 geometry - in
4 // particular the support cone 
5 //
6 //____________________________________________________________________
7 TObjArray*
8 waferParameters(double dl, double dh, double theta, double r, 
9                 Bool_t verbose=kFALSE)
10 {      
11   double tan_theta  = tan(theta * TMath::Pi() / 180.);
12   double tan_theta2 = pow(tan_theta,2);
13   double r2         = pow(r,2);
14   double y_A        = tan_theta * dl;
15   double x_D        = dl + sqrt(r2 - tan_theta2 * pow(dl,2));
16   double x_D_2      = dl - sqrt(r2 - tan_theta2 * pow(dl,2));
17   
18   double y_B       = sqrt(r2 - pow(dh,2) + 2 * dh * x_D - pow(x_D,2));
19   double x_C       = (x_D + sqrt(-tan_theta2 * pow(x_D,2) + r2 
20                                  + r2 * tan_theta2)) / (1 + tan_theta2);
21   double y_C       = tan_theta * x_C;
22
23   if (verbose) {
24     cout << "A: (" << dl << "," << y_A << ")" << endl;
25     cout << "B: (" << dh << "," << y_B << ")" << endl;
26     cout << "C: (" << x_C << "," << y_C << ")" << endl;
27     cout << "D: (" << x_D << ",0)" << endl;
28     
29     cout << "Recentred at D:"  << endl;
30     cout << "A: (" << dl - x_D  << "," << y_A << ")" << endl;
31     cout << "B: (" << dh - x_D  << "," << y_B << ")" << endl;
32     cout << "C: (" << x_C - x_D << "," << y_C << ")" << endl;
33   }
34   
35   TObjArray* verticies = new TObjArray(6);
36   verticies->AddAt(new TVector2(dl,   y_A), 0);
37   verticies->AddAt(new TVector2(x_C,  y_C), 1);
38   verticies->AddAt(new TVector2(dh,   y_B), 2);
39   verticies->AddAt(new TVector2(dh,  -y_B), 3);
40   verticies->AddAt(new TVector2(x_C, -y_C), 4);
41   verticies->AddAt(new TVector2(dl,  -y_A), 5);
42   
43   return verticies;
44 }
45
46 //____________________________________________________________________
47 TShape* 
48 createModuleShape(const Char_t* name, double rl, double rh, double th, 
49                   double r, double dz) 
50 {
51   std::cout << "Creating Module shape for " << name << std::flush;
52   // TShape* virtualShape = new TTUBS(Form("%sVirtual", name), 
53   //                                  Form("Virtual %s", name), 
54   //                                  "", rl, rh, 1, -th, th);
55   
56   TObjArray* v = waferParameters(rl, rh, th, r);
57   TXTRU* moduleShape = new TXTRU(Form("%sModule", name),
58                                  Form("Module %s", name), 
59                                  "", 6, 2);
60   for (Int_t i = 0; i  < 6; i++) {
61     std::cout << "." << std::flush;
62     TVector2* vv = static_cast<TVector2*>(v->At(i));
63     moduleShape->DefineVertex(i, vv->X(), vv->Y());
64   }
65   moduleShape->DefineSection(0, -dz, 1, 0, 0);
66   moduleShape->DefineSection(1,  dz, 1, 0, 0);
67   std::cout << std::endl;
68
69   return (TShape*)moduleShape;
70 }
71
72 //____________________________________________________________________
73 TNode* 
74 createRing(const char* name, double rl, double rh, double th, 
75            double siThick, double waferR, double staggering, double z) 
76 {
77   std::cout << "Creating Ring node for " << name << std::flush;
78   TShape* ringShape   = new TTUBE(Form("%sShape", name), "Ring Shape", 
79                                   "", rl, rh, staggering + siThick);
80   TNode*  ringNode    = new TNode(Form("%sNode", name), "Ring Node", 
81                                   ringShape, 0, 0, z, 0);
82   TShape* moduleShape = createModuleShape(name, rl, rh, th, waferR, siThick);
83   Int_t n = 360 / 2 / th;
84   for (Int_t i = 0; i < n; i++) {
85     std::cout << "." << std::flush;
86     ringNode->cd();
87     Double_t theta = 2  * th * i;
88     Double_t z     = (i % 2 ? 0 : staggering);
89     TRotMatrix* rot = new TRotMatrix(Form("%sRotation%02d", name, i), 
90                                      "Rotation", 90, theta, 90, 
91                                      fmod(90 + theta, 360), 0, 0);
92     TNode* moduleNode = new TNode(Form("%sModule%02d", name, i), 
93                                   "Module", moduleShape, 0, 0, z,
94                                   rot);
95     moduleNode->SetFillColor(5);
96     moduleNode->SetLineColor(5);
97     moduleNode->SetLineWidth(2);
98   }
99   std::cout << std::endl;
100   ringNode->SetVisibility(0);
101   return ringNode;
102 }
103
104 //____________________________________________________________________
105 TNode*
106 createSupport(double noseRl, double noseRh, double noseDz, double noseZ, 
107               double backRl, double backRh, double backDz, double coneL, 
108               double beamW, double beamDz,  double flangeR) 
109 {
110   TShape* noseShape = new TTUBE("noseShape", "Nose Shape", "", 
111                                 noseRl, noseRh, noseDz);
112   TNode*  noseNode  = new TNode("noseNode", "noseNode", noseShape, 
113                                 0, 0, noseZ - noseDz, 0);
114   noseNode->SetLineColor(0);
115   
116   Double_t zdist = coneL - 2 * backDz - 2 * noseDz;
117   Double_t tdist = backRh - noseRh;
118   Double_t beamL = TMath::Sqrt(zdist * zdist + tdist * tdist);
119   Double_t theta = -TMath::ATan2(tdist, zdist);
120   
121
122   TShape* backShape = new TTUBE("backShape", "Back Shape", "", 
123                                  backRl, backRh, backDz);
124   TNode*  backNode  = new TNode("backNode", "backNode", backShape, 
125                                 0, 0, noseZ - 2 * noseDz - zdist - backDz, 0);
126   backNode->SetLineColor(0);
127
128
129   TShape* beamShape = new TBRIK("beamShape", "beamShape", "", 
130                                 beamDz, beamW / 2 , beamL / 2);
131   Int_t    n = 8;
132   Double_t r = noseRl + tdist / 2;
133   for (Int_t i = 0; i < n; i++) {
134     Double_t phi   = 360. / n * i;
135     Double_t t     = 180. * theta / TMath::Pi();
136     TRotMatrix* beamRotation = new TRotMatrix(Form("beamRotation%d", i), 
137                                                 Form("beamRotation%d", i),
138                                                 180-t, phi, 90, 90+phi, 
139                                                 t, phi);
140     TNode* beamNode = new TNode(Form("beamNode%d", i), 
141                                 Form("beamNode%d", i), beamShape, 
142                                 r * TMath::Cos(phi / 180 * TMath::Pi()), 
143                                 r * TMath::Sin(phi / 180 * TMath::Pi()), 
144                                 noseZ - 2 * noseDz - zdist / 2,  beamRotation);
145     beamNode->SetLineColor(0);
146   }
147   
148   Double_t flangel = (flangeR - backRh) / 2;
149   TShape* flangeShape = new TBRIK("flangeShape", "FlangeShape", "", 
150                                   flangel, beamW / 2, backDz);
151   n = 4;
152   r = backRh + flangel;
153   for (Int_t i = 0; i < n; i++) {
154     Double_t phi = 360. / n * i + 180. / n;
155     TRotMatrix* flangeRotation = new TRotMatrix(Form("flangeRotation%d", i),
156                                                 Form("Flange Rotation %d", i),
157                                                 90, phi, 90, 90+phi, 0, 0);
158     TNode* flangeNode = new TNode(Form("flangeNode%d", i), 
159                                   Form("flangeNode%d", i), 
160                                   flangeShape,
161                                   r * TMath::Cos(phi / 180 * TMath::Pi()), 
162                                   r * TMath::Sin(phi / 180 * TMath::Pi()),
163                                   noseZ - 2 * noseDz - zdist - backDz, 
164                                   flangeRotation);
165     flangeNode->SetLineColor(0);
166                                   
167   }
168   return 0;
169 }
170
171                                  
172
173
174 //____________________________________________________________________
175 void
176 SimpleGeometry() 
177 {
178   TGeometry* geometry = new TGeometry("geometry","geometry");
179   TShape* topShape = new TBRIK("topShape", "topShape", "", 40, 40, 150);
180   TNode* topNode = new TNode("topNode", "topNode", topShape, 0, 0, 0, 0);
181   topNode->SetVisibility(0);
182   topNode->cd();
183   
184   Double_t waferR     = 13.4 / 2;
185   Double_t siThick    = .03;
186   Double_t staggering = 1;
187   Double_t innerRh    = 17.2;
188   Double_t innerRl    = 4.3;
189   Double_t innerTh    = 18;
190   Double_t innerZ     = -62.8;
191   Double_t outerRh    = 28;
192   Double_t outerRl    = 15.6;
193   Double_t outerTh    = 9;
194   Double_t outerZ     = -75.2;
195   Double_t noseRl     = 5.5;
196   Double_t noseRh     = 6.7;
197   Double_t noseDz     = 2.8 / 2;
198   Double_t noseZ      = -46;
199   Double_t coneL      = 30.9;
200   Double_t backRl     = 61 / 2;
201   Double_t backRh     = 66.8 /2;
202   Double_t backDz     = 1.4 / 2;
203   Double_t beamDz     = .5 / 2;
204   Double_t beamW      = 6;
205   Double_t flangeR    = 49.25;
206
207   Double_t zdist      = coneL - 2 * backDz - 2 * noseDz;
208   Double_t tdist      = backRh - noseRh;
209   Double_t alpha      = tdist / zdist;
210
211   Double_t x, rl, rh, z;
212   z  = noseZ - coneL / 2;
213   TPCON* fmd3Shape = new TPCON("fmd3Shape", "FMD 3 Shape", "", 0, 360, 7);
214   
215   x  = noseZ;
216   rl = noseRl;
217   rh = noseRh;
218   fmd3Shape->DefineSection(0, x - z, rl, rh);
219
220   x  = noseZ-2*noseDz;
221   fmd3Shape->DefineSection(1, x - z, rl, rh);
222
223   x  = innerZ - staggering - siThick; 
224   rl = innerRl;
225   rh = noseRh + alpha * TMath::Abs(x-noseZ + 2 * noseDz);
226   fmd3Shape->DefineSection(2, x - z, rl, rh);
227
228   x  = outerZ;
229   rl = outerRl;
230   rh = backRh;
231   fmd3Shape->DefineSection(3, x - z, rl, rh);
232
233   x  = noseZ - zdist - 2 * noseDz;
234   rl = outerRl;
235   rh = backRh;
236   fmd3Shape->DefineSection(4, x - z, rl, rh);
237
238   x  = noseZ - zdist - 2 * noseDz;
239   rl = outerRl;
240   rh = flangeR;
241   fmd3Shape->DefineSection(5, x - z, rl, rh);
242
243   x  = noseZ - coneL;
244   rl = outerRl;
245   rh = flangeR;
246   fmd3Shape->DefineSection(6, x - z, rl, rh);
247
248   TNode* fmd3Node = new TNode("fmd3Node", "FMD3 Node", fmd3Shape, 
249                               0, 0, z, 0);
250   fmd3Node->SetLineColor(3);
251   fmd3Node->SetVisibility(1);
252
253   fmd3Node->cd();
254   TNode* innerNode = createRing("inner", innerRl, innerRh, innerTh, 
255                                 siThick, waferR, staggering, innerZ-z);
256
257
258   fmd3Node->cd();
259   TNode* outerNode = createRing("outer", outerRl, outerRh, outerTh, 
260                                 siThick, waferR, staggering, outerZ-z);
261   
262
263   fmd3Node->cd();
264   TNode* supportNode = createSupport(noseRl, noseRh, noseDz, noseZ-z, 
265                                      backRl, backRh, backDz, coneL,
266                                      beamW, beamDz, flangeR);
267   
268   TCanvas* c = new TCanvas("c", "c", 800, 800);
269   c->SetFillColor(1);
270   geometry->Draw();
271   // c->x3d("ogl");
272 }