FairRoot/PandaRoot
createSTT.C
Go to the documentation of this file.
1 // *******************************************************************************
2 // OLD (!!!) MACRO TO PRODUCE THE GEOMETRY ASCII FILE
3 // *******************************************************************************
4 // if you want to reduce the tube length:
5 // - change tubeLength
6 // - change sttCenterZ
7 // - length and number of the short skewed straws is changed automatically
8 // *******************************************************************************
9 // to do: get media right
10 // add targetpipe support
11 // writeout function
12 // names of volumes
13 // *******************************************************************************
14 // to reproduce the file sstraws_skewed_blocks_pipe_120cm.geo use rev 7746
15 // *******************************************************************************
16 
17 #include <iostream>
18 #include <iomanip>
19 #include <fstream>
20 #include <math.h>
21 #include <string>
22 #include <sstream>
23 #include <vector>
24 
25 using namespace std;
26 
27 static int counter1;
28 static int counter2;
29 static int counter3;
30 static int counter4;
31 static int counter5;
32 static int counter6;
33 
34 //--------------------------------------------geometry definitions
35 #define outerDiam 84.4000 // 84.0000
36 #define innerDiam 30.0000 // 32.0000
37 #define tubeInnerDiam 1.0000
38 #define tubeOuterDiam 1.0060
39 // #define tubeSeperation 1.0100 // <<=== this HAS to be the tubeOuterDiam!! CHECK
41 #define wireDiam 0.0020
42 #define tubeLength 120.0000 // 150.0000
43 
44 #define safety 0.2000
45 
46 #define sttCenterX 0.
47 #define sttCenterY 0.
48 // tubeLength/2 - 40
49 #define sttCenterZ 20. // 35.
50 
51 #define innerCoverThickness 1.0000 // 1.200 // 0.1000 cm
52 #define outerCoverThickness 1.2000 // 1.200 // 0.1000 cm
53 
54 #define panelthickness 0.1000
55 #define pipeDiam 4.0800
56 #define noSupport 0
57 
58 #define skewangle 3.
59 
60 #define pi 3.141592653589793238512808959406186204433
61 
62 #define cave "cave"
63 #define sttassembly "stt01assembly"
64 #define innerCylinder "stt01innerCylinder"
65 #define outerCylinder "stt01outerCylinder"
66 #define panel1 "stt01box#1"
67 #define panel2 "stt01box#2"
68 #define panel3 "stt01box#3"
69 #define panel4 "stt01box#4"
70 
71 // materials
72 #define air "air"
73 #define AlBe "carbon"
74 #define mylar "mylar"
75 #define HeMixture "argon"
76 #define W "copper"
77 
78 static int tubeteller;
79 static int axialtubeteller;
80 static int skewedtubeteller;
82 static double maximumradius = 0;
83 static double minimumradius = 100000;
84 
85 // fuctions to write part of the declaration
86 void writename(char const *name, bool support = false, bool leftside = false)
87 {
88  if(support && leftside) cout << name << "#1" << endl;
89  else if(support && !leftside) cout << name << "#2" << endl;
90  else cout << name << endl;
91  counter1++;
92 }
93 
94 void writemother(char const *name, bool original = true, bool support = false)
95 {
96  cout << name << endl;
97  if (original && !support)
98  cout << "TUBE" << endl;
99  if(original && support)
100  cout << "TUBS" << endl;
101  counter2++;
102 }
103 
104 void writemedium(char *name)
105 {
106  cout << name << endl;
107  counter3++;
108 }
109 
110 void writetube(double inner, double outer, double length)
111 {
112  cout << 0. << " " << 0. << " " << -1. * length * 10. << endl;
113  cout << inner * 10. << " " << outer * 10. << endl;
114  cout << 0. << " " << 0. << " " << length * 10. << endl;
115  counter4++;
116 }
117 
118 void writehalftube(double inner, double outer, double length)
119 {
120  cout << 0. << " " << 0. << " " << -1. * length * 10. << endl;
121  cout << inner * 10. << " " << outer * 10. << endl;
122  cout << 0. << " " << 0. << " " << length * 10. << endl;
123  cout << atan((pipeDiam/2.)/outer) * (180. / pi) << " " << 180. - atan((pipeDiam/2.)/outer) * (180. / pi) << endl;
124  counter4++;
125 }
126 
127 void writepanel(char const *name, bool firstone, double xthick, double ythick, double length, int side)
128 {
129  xthick *= 10.;
130  ythick *= 10.;
131  length *= 10.;
132 
133  if(firstone) {
134  cout << name << endl;
135  cout << sttassembly << endl;
136  cout << "BOX" << endl;
137  cout << AlBe << endl;
138  cout << xthick/2. << " " << ythick/2. << " " << length / 2. << endl;
139  cout << -xthick/2. << " " << ythick/2. << " " << length / 2. << endl;
140  cout << xthick/2. << " " << ythick/2. << " " << -length / 2. << endl;
141  cout << -xthick/2. << " " << ythick/2. << " " << -length / 2. << endl;
142  cout << xthick/2. << " " << -ythick/2. << " " << length / 2. << endl;
143  cout << -xthick/2. << " " << -ythick/2. << " " << length / 2. << endl;
144  cout << xthick/2. << " " << -ythick/2. << " " << -length / 2. << endl;
145  cout << -xthick/2. << " " << -ythick/2. << " " << -length / 2. << endl;
146  }
147  else {
148  cout << name << endl;
149  cout << sttassembly << endl;
150  }
151 
152 }
153 
154 
155 void writetrans(double x, double y, double z)
156 {
157  cout << x * 10. << " " << y * 10. << " " << z * 10. << endl;
158  counter5++;
159 }
160 
161 void writerot(double x00, double x01, double x02, double x10, double x11, double x12, double x20, double x21, double x22)
162 {
163  cout << x00 << " " << x01 << " " << x02 << " "
164  << x10 << " " << x11 << " " << x12 << " "
165  << x20 << " " << x21 << " " << x22 << endl;
166  cout << "//----------------------------------------------------------" << endl;
167  counter6++;
168 }
169 
170 // write the declaration of a complete straw
171 bool putStraw(double posX, double posY, double posZ)
172 {
173 
174  // CHECK PIPE
175  if(posX > -(pipeDiam/2. + panelthickness + tubeOuterDiam/2.) && posX < (pipeDiam/2. + panelthickness + tubeOuterDiam/2.)) return false;
176 
177  // check if the straws fit wihin the inner and outer diameter specified
178  if (
179  (sqrt(posX * posX + posY * posY) < ((outerDiam / 2.) - outerCoverThickness - (tubeOuterDiam / 2.))) &&
180  (sqrt(posX * posX + posY * posY) > ((innerDiam / 2.) + innerCoverThickness + (tubeOuterDiam / 2.)))
181  )
182  {
183  // convert int to string
184  stringstream
185  conv;
186 
187  string
188  tubetellerStr;
189 
190  conv << tubeteller + 1;
191  conv >> tubetellerStr;
192 
193  string
194  nameItube = "stt01tube#" + tubetellerStr,
195  nameIgas = "stt01gas#" + tubetellerStr,
196  nameIwire = "stt01wire#" + tubetellerStr;
197 
198  // keep track of minimum and maximum extend of the straw package
199  if (sqrt(posX * posX + posY * posY) - (tubeOuterDiam / 2.) < minimumradius)
200  {
201  minimumradius = sqrt(posX * posX + posY * posY) - (tubeOuterDiam / 2.);
202  }
203  if (sqrt(posX * posX + posY * posY) + (tubeOuterDiam / 2.) > maximumradius)
204  {
205  maximumradius = sqrt(posX * posX + posY * posY) + (tubeOuterDiam / 2.);
206  }
207 
208  // place the tube
209  if (tubeteller == 0)
210  {
211  // original tube
212 
213  // Mylar tubes
214  writename(nameItube.c_str());
217  writetube(0, tubeOuterDiam / 2., tubeLength / 2.);
218  writetrans(posX, posY, posZ);
219  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
220 
221  // Gas filling
222  writename(nameIgas.c_str());
223  writemother(nameItube.c_str());
225  writetube(0., tubeInnerDiam / 2., tubeLength / 2.);
226  writetrans(0., 0., 0.);
227  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
228 
229  // Anode wire
230  writename(nameIwire.c_str());
231  writemother(nameIgas.c_str());
232  writemedium(W);
233  writetube(0., wireDiam / 2., tubeLength / 2.);
234  writetrans(0., 0., 0.);
235  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
236  }
237  else
238  {
239  // copy
240  writename(nameItube.c_str());
241  writemother(sttassembly, false);
242  writetrans(posX, posY, posZ);
243  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
244  }
245  tubeteller++;
246  axialtubeteller++;
247  }
248  else
249  {
250  return false;
251  }
252 
253  return true;
254 }
255 
256 // write the rotation matrix for a right skewed straw
257 void plotrotright(double x, double y, double z)
258 {
259  double
260  newskewangle = skewangle * -1.;
261 
262  cout << 1 + (1 - cos(newskewangle * (pi / 180.))) * (x * x - 1) << " "
263  << -z * sin(newskewangle * (pi / 180.)) + (1 - cos(newskewangle * (pi / 180.))) * x * y << " "
264  << y * sin(newskewangle * (pi / 180.)) + (1 - cos(newskewangle * (pi / 180.))) * x * z << " "
265 
266  << z * sin(newskewangle * (pi / 180.)) + (1 - cos(newskewangle * (pi / 180.))) * x * y << " "
267  << 1 + (1 - cos(newskewangle * (pi / 180.))) * (y * y - 1) << " "
268  << -x * sin(newskewangle * (pi / 180.)) + (1 - cos(newskewangle * (pi / 180.))) * y * z << " "
269 
270  << -y * sin(newskewangle * (pi / 180.)) + (1 - cos(newskewangle * (pi / 180.))) * x * z << " "
271  << x * sin(newskewangle * (pi / 180.)) + (1 - cos(newskewangle * (pi / 180.))) * y * z << " "
272  << 1 + (1 - cos(newskewangle * (pi / 180.))) * (z * z - 1) << endl;
273  cout << "//----------------------------------------------------------" << endl;
274  counter6++;
275 }
276 
277 // write the rotation matrix for a left skewed straw
278 void plotrotleft(double x, double y, double z)
279 {
280  cout << 1 + (1 - cos(skewangle * (pi / 180.))) * (x * x - 1) << " "
281  << -z * sin(skewangle * (pi / 180.)) + (1 - cos(skewangle * (pi / 180.))) * x * y << " "
282  << y * sin(skewangle * (pi / 180.)) + (1 - cos(skewangle * (pi / 180.))) * x * z << " "
283 
284  << z * sin(skewangle * (pi / 180.)) + (1 - cos(skewangle * (pi / 180.))) * x * y << " "
285  << 1 + (1 - cos(skewangle * (pi / 180.))) * (y * y - 1) << " "
286  << -x * sin(skewangle * (pi / 180.)) + (1 - cos(skewangle * (pi / 180.))) * y * z << " "
287 
288  << -y * sin(skewangle * (pi / 180.)) + (1 - cos(skewangle * (pi / 180.))) * x * z << " "
289  << x * sin(skewangle * (pi / 180.)) + (1 - cos(skewangle * (pi / 180.))) * y * z << " "
290  << 1 + (1 - cos(skewangle * (pi / 180.))) * (z * z - 1) << endl;
291  cout << "//----------------------------------------------------------" << endl;
292  counter6++;
293 }
294 
295 bool putStrawRotatedShortLeft(double posX, double posY, double posZ, double xvector, double yvector, double zvector, double length)
296 {
297 
298  double extr1 = posX + (length/2.) * sin(skewangle);
299  double extr2 = posX - (length/2.) * sin(skewangle);
300 
301  if (length > 0.)
302  {
303  // convert int to string
304  stringstream
305  conv;
306 
307  string
308  tubetellerStr;
309 
310  conv << tubeteller + 1;
311  conv >> tubetellerStr;
312 
313  string
314  nameItube = "stt01tube" + tubetellerStr,
315  nameIgas = "stt01gas" + tubetellerStr,
316  nameIwire = "stt01wire" + tubetellerStr;
317 
318  // N.B. we don't use copies here since the short straws are all different in length
319 
320  // Mylar tubes
321  writename(nameItube.c_str());
324  writetube(0, tubeOuterDiam / 2., length / 2.);
325  writetrans(posX, posY, posZ);
326  plotrotleft(xvector, yvector, zvector);
327 
328  // gas filling
329  writename(nameIgas.c_str());
330  writemother(nameItube.c_str());
332  writetube(0., tubeInnerDiam / 2., length / 2.);
333  writetrans(0., 0., 0.);
334  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
335 
336  // anode wires
337  writename(nameIwire.c_str());
338  writemother(nameIgas.c_str());
339  writemedium(W);
340  writetube(0., wireDiam / 2., length / 2.);
341  writetrans(0., 0., 0.);
342  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
343 
344  tubeteller++;
346 
347 
348  return true;
349  }
350 }
351 
352 bool putStrawRotatedShortRight(double posX, double posY, double posZ, double xvector, double yvector, double zvector, double length)
353 {
354  double extr1 = posX + (length/2.) * sin(skewangle);
355  double extr2 = posX - (length/2.) * sin(skewangle);
356 
357  if (length > 0.)
358  {
359 
360 
361 
362  // convert int to string
363  stringstream
364  conv;
365 
366  string
367  tubetellerStr;
368 
369  conv << tubeteller + 1;
370  conv >> tubetellerStr;
371 
372  string
373  nameItube = "stt01tube" + tubetellerStr,
374  nameIgas = "stt01gas" + tubetellerStr,
375  nameIwire = "stt01wire" + tubetellerStr;
376 
377  // N.B. don;t use copies here since all short straws have different length
378 
379  // Mylar tubes
380  writename(nameItube.c_str());
383  writetube(0, tubeOuterDiam / 2., length / 2.);
384  writetrans(posX, posY, posZ);
385  plotrotright(xvector, yvector, zvector);
386 
387  // gas filling
388  writename(nameIgas.c_str());
389  writemother(nameItube.c_str());
391  writetube(0., tubeInnerDiam / 2., length / 2.);
392  writetrans(0., 0., 0.);
393  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
394 
395  // anode wire
396  writename(nameIwire.c_str());
397  writemother(nameIgas.c_str());
398  writemedium(W);
399  writetube(0., wireDiam / 2., length / 2.);
400  writetrans(0., 0., 0.);
401  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
402 
403  tubeteller++;
405 
406  return true;
407  }
408 }
409 
410 bool putStrawRotatedLeft(double posX, double posY, double posZ, double xvector, double yvector, double zvector, double length)
411 {
412  double extr1 = posX + (length/2.) * sin(skewangle);
413  double extr2 = posX - (length/2.) * sin(skewangle);
414 
415 
416  // convert int to string
417  stringstream
418  conv;
419 
420  string
421  tubetellerStr;
422 
423  conv << tubeteller + 1;
424  conv >> tubetellerStr;
425 
426  string
427  nameItube = "stt01tube#" + tubetellerStr,
428  nameIgas = "stt01gas#" + tubetellerStr,
429  nameIwire = "stt01wire#" + tubetellerStr;
430 
431  if (tubeteller == 0)
432  {
433  // original straw
434  // Mylar tubes
435  writename(nameItube.c_str());
438  writetube(0, tubeOuterDiam / 2., length / 2.);
439  writetrans(posX, posY, posZ);
440  plotrotleft(xvector, yvector, zvector);
441 
442  // Gas filling
443  writename(nameIgas.c_str());
444  writemother(nameItube.c_str());
446  writetube(0., tubeInnerDiam / 2., length / 2.);
447  writetrans(0., 0., 0.);
448  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
449 
450  // Anode wire
451  writename(nameIwire.c_str());
452  writemother(nameIgas.c_str());
453  writemedium(W);
454  writetube(0., wireDiam / 2., length / 2.);
455  writetrans(0., 0., 0.);
456  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
457  }
458  else
459  {
460  // copy volume
461 
462  writename(nameItube.c_str());
463  writemother(sttassembly, false);
464  writetrans(posX, posY, posZ);
465  plotrotleft(xvector, yvector, zvector);
466  }
467 
468  tubeteller++;
470 
471  return true;
472 }
473 
474 bool putStrawRotatedRight(double posX, double posY, double posZ, double xvector, double yvector, double zvector, double length)
475 {
476 
477 
478  double extr1 = posX + (length/2.) * sin(skewangle);
479  double extr2 = posX - (length/2.) * sin(skewangle);
480 
481 
482  // convert int to string
483  stringstream
484  conv;
485 
486  string
487  tubetellerStr;
488 
489  conv << tubeteller + 1;
490  conv >> tubetellerStr;
491 
492  string
493  nameItube = "stt01tube#" + tubetellerStr,
494  nameIgas = "stt01gas#" + tubetellerStr,
495  nameIwire = "stt01wire#" + tubetellerStr;
496 
497  if (tubeteller == 0)
498  {
499  // original volume
500 
501  // Mylar tubes
502  writename(nameItube.c_str());
505  writetube(0, tubeOuterDiam / 2., length / 2.);
506  writetrans(posX, posY, posZ);
507  plotrotright(xvector, yvector, zvector);
508 
509  // gas filling
510  writename(nameIgas.c_str());
511  writemother(nameItube.c_str());
513  writetube(0., tubeInnerDiam / 2., length / 2.);
514  writetrans(0., 0., 0.);
515  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
516 
517  // anode wire
518  writename(nameIwire.c_str());
519  writemother(nameIgas.c_str());
520  writemedium(W);
521  writetube(0., wireDiam / 2., length / 2.);
522  writetrans(0., 0., 0.);
523  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
524  }
525  else
526  {
527  // copy volume
528  writename(nameItube.c_str());
529  writemother(sttassembly, false);
530  writetrans(posX, posY, posZ);
531  plotrotright(xvector, yvector, zvector);
532  }
533 
534  tubeteller++;
536 
537  return true;
538 }
539 
540 
541 
542 void placeSingleLayerStraightExact(double ringteller)
543 {
544  double
545  radius = tubeSeperation / 2.;
546 
547  // place a straight double layer
548  double
549  xpos = 0.,
550  ypos = ((ringteller) * 2 * radius),
551  zpos = 0.;
552 
553  // 6 loops, one for each side of the hexagon
554  // 1 / \ 4
555  // 2 | | 5
556  // 3 \ / 6
557  for(int i = 0; i < ringteller; i++)
558  {
559  xpos -= sqrt(3.) * radius;
560  ypos -= radius;
561  putStraw(xpos, ypos, zpos);
562  }
563  for(int i = 0; i < ringteller; i++)
564  {
565  ypos -= 2 * radius;
566  putStraw(xpos, ypos, zpos);
567  }
568  for(int i = 0; i < ringteller; i++)
569  {
570  xpos += sqrt(3.) * radius;
571  ypos -= radius;
572  putStraw(xpos, ypos, zpos);
573  }
574  for(int i = 0; i < ringteller; i++)
575  {
576  xpos += sqrt(3.) * radius;
577  ypos += radius;
578  putStraw(xpos, ypos, zpos);
579  }
580  for(int i = 0; i < ringteller; i++)
581  {
582  ypos += 2 * radius;
583  putStraw(xpos, ypos, zpos);
584  }
585  for(int i = 0; i < ringteller; i++)
586  {
587  xpos -= sqrt(3.) * radius;
588  ypos += radius;
589  putStraw(xpos, ypos, zpos);
590  }
591 }
592 
593 void placeSingleLayerSkewedRight(double ringposition)
594 {
595  // see void placeSingleLayerSkewedLeft(double ringposition)
596 
597  double radius = tubeSeperation / 2.;
598 
599  double
600  newskewangle = skewangle * -1.;
601 
602  double
603  newradius = radius / cos((newskewangle / 180.) * pi);
604 
605  double
606  xpos = 0.,
607  ypos = ringposition / cos((30. / 180.) * pi),
608  zpos = 0.,
609  tmpxpos,
610  tmpypos,
611  vectorx,
612  vectory,
613  vectorz;
614 
615  double
616  availableSpace = 2. * ringposition * tan((30. / 180.) * pi);
617 
618  int
619  possibleStraws = int((availableSpace - fabs(tubeLength * sin((newskewangle / 180.) * pi))) / (newradius * 2)),
620  extraStraws = -1 * possibleStraws + int(availableSpace / (newradius * 2));
621 
622  double
623  extraspace = (availableSpace - (possibleStraws + extraStraws) * (newradius * 2)) / 2.;
624 
625  cerr << "extra space: " << extraspace * 2 << endl;
626 
627  double
628  skewanglerad = (newskewangle / 180.) * pi,
629  sixtyrad = (60. / 180.) * pi,
630  translationToLeft = - 0.5 * tubeLength * tan(skewanglerad);
631 
632 
633  vector<double>
634  lengthsShort;
635 
636  double
637  xpos2 = xpos - availableSpace * sin((60. / 180.) * pi),
638  ypos2 = ypos - availableSpace * cos((60. / 180.) * pi),
639  tmpxpos2 = xpos2,
640  tmpypos2 = ypos2 - extraspace - newradius,
641  limit = ypos2 - availableSpace + newradius + safety + extraspace;
642 
643  for(int i = 0; i < possibleStraws; i++)
644  {
645  tmpypos2 -= 2 * newradius;
646  }
647 
648  int extracounter = 0; // how many short tubes can be really placed
649  for (int i = 0; i < extraStraws; i++)
650  {
651  if((tmpypos2 - limit) < 0) break;
652  extracounter++;
653  double
654  lengthShort = fabs((tmpypos2 - newradius - limit) / sin(newskewangle * (pi / 180.)));
655 
656  if (lengthShort > tubeLength)
657  lengthShort = tubeLength;
658 
659  // cerr << "length: " << lengthShort << endl;
660  lengthsShort.push_back(fabs(lengthShort));
661 
662  tmpypos2 -= 2 * newradius;
663  }
664  extraStraws = extracounter; // update the number
665 
666  double
667  additionalShift = fabs(tubeLength * sin((newskewangle / 180.) * pi));
668 
669  // =================================================
670  // planes intersecting pipe ========================
671  // pipe
672  double pipespace = (panelthickness + pipeDiam/2.) * 1./(sqrt(3)/2.);
673 
674  // the available space for a layer at the radius ringposition
675  double
676  availableSpacePipe = 2. * ringposition * tan((30. / 180.) * pi) - pipespace;
677 
678  // the number of skewed full-length straws that fit in the available space
679  int
680  possibleStrawsPipe = int((availableSpacePipe - fabs(tubeLength * sin((skewangle / 180.) * pi))) / (newradius * 2));
681 
682  // the number of short straws that fit on either side of the pack of full-length straws
683  int
684  extraStrawsPipe = -1 * possibleStrawsPipe + int(availableSpacePipe / (newradius * 2));
685  // the amount of space that is left in the layer when considering the sum of all
686  // full length and short straws together
687  double
688  extraspacePipe = (availableSpacePipe - (possibleStrawsPipe + extraStrawsPipe) * (newradius * 2)) / 2.;
689 
690  // cerr << "layer with pipe " << possibleStraws << " " << extraStraws << endl;
691 
692  // calculate the lenghts of the different short straws in this layer
693  vector<double>
694  lengthsShortPipe;
695 
696  double
697  // start at the right edge of the layer
698  xpos2Pipe = xpos - availableSpacePipe * sin((60. / 180.) * pi),
699  ypos2Pipe = ypos - availableSpacePipe * cos((60. / 180.) * pi),
700  // shift to the center of the rightmost straw
701  tmpxpos2Pipe = xpos2Pipe,
702  tmpypos2Pipe = ypos2Pipe - extraspacePipe - newradius,
703  // the limit on the left side of the plane beyond which no center of the short straws may be
704  limitPipe = ypos2Pipe - availableSpacePipe + newradius + safety + extraspacePipe;
705 
706 
707  // shift to the center of of the leftmost full-length straw
708  for(int i = 0; i < possibleStrawsPipe; i++)
709  {
710  tmpypos2Pipe -= 2 * newradius;
711  }
712 
713  // calculate the maximum allowable length for this short straw and shift to the next
714  extracounter = 0; // how many short tubes can be really placed
715  for (int i = 0; i < extraStrawsPipe; i++)
716  {
717  if((tmpypos2Pipe - limitPipe) < 0) break;
718  extracounter++;
719  double
720  lengthShort = fabs((tmpypos2Pipe - newradius - limitPipe) / sin(skewangle * (pi / 180.)));
721 
722  if (lengthShort > tubeLength)
723  lengthShort = tubeLength;
724  lengthsShortPipe.push_back(lengthShort);
725  tmpypos2Pipe -= 2 * newradius;
726  }
727  extraStrawsPipe = extracounter; // update the number
728 
729  // =================================================
730 
731  // plane 1
732  vectorx = -cos((60. / 180.) * pi);
733  vectory = sin((60. / 180.) * pi);
734  vectorz = 0;
735  tmpxpos = xpos - (extraspacePipe + newradius + pipespace) * sin((60. / 180.) * pi);
736  tmpypos = ypos - (extraspacePipe + newradius + pipespace) * cos((60. / 180.) * pi);
737 
738  for(int i = 0; i < possibleStrawsPipe; i++)
739  {
740 
741  putStrawRotatedRight(tmpxpos + (translationToLeft - additionalShift) * sin(sixtyrad),
742  tmpypos + (translationToLeft - additionalShift) * cos(sixtyrad),
743  zpos, vectorx, vectory, vectorz, tubeLength);
744 
745  tmpxpos -= 2 * newradius * sin((60. / 180.) * pi);
746  tmpypos -= 2 * newradius * cos((60. / 180.) * pi);
747  }
748 
749  for (int i = 0; i < extraStrawsPipe; i++)
750  {
751  double
752  lengthShort = lengthsShortPipe[i];
753 
754  double
755  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
756  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
757  switchSides = 2 * newradius * (possibleStrawsPipe + 1 + (2 * i)) * cos(skewanglerad);
758 
759  putStrawRotatedShortRight(tmpxpos + (translationToLeft - additionalShift - translationToFrontPer) * sin(sixtyrad),
760  tmpypos + (translationToLeft - additionalShift - translationToFrontPer) * cos(sixtyrad),
761  zpos - translationToFrontPar,
762  vectorx, vectory, vectorz, lengthShort);
763 
764  putStrawRotatedShortRight(tmpxpos + (translationToLeft - additionalShift + translationToFrontPer + switchSides) * sin(sixtyrad),
765  tmpypos + (translationToLeft - additionalShift + translationToFrontPer + switchSides) * cos(sixtyrad),
766  zpos + translationToFrontPar,
767  vectorx, vectory, vectorz, lengthShort);
768 
769 
770  tmpxpos -= 2 * newradius * sin((60. / 180.) * pi);
771  tmpypos -= 2 * newradius * cos((60. / 180.) * pi);
772  }
773 
774 
775  // plane 2
776  vectorx = -1.;
777  vectory = 0.;
778  vectorz = 0.;
779  xpos -= availableSpace * sin((60. / 180.) * pi);
780  ypos -= availableSpace * cos((60. / 180.) * pi);
781  tmpxpos = xpos;
782  tmpypos = ypos - extraspace - newradius;
783 
784  for(int i = 0; i < possibleStraws; i++)
785  {
786  putStrawRotatedRight(tmpxpos, tmpypos - additionalShift + translationToLeft, zpos, vectorx, vectory, vectorz, tubeLength);
787  tmpypos -= 2 * newradius;
788  }
789 
790  for (int i = 0; i < extraStraws; i++)
791  {
792  double
793  lengthShort = lengthsShort[i];
794 
795  double
796  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
797  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
798  switchSides = 2 * newradius * (possibleStraws + 1 + (2 * i)) * cos(skewanglerad);
799 
800  putStrawRotatedShortRight(tmpxpos,
801  tmpypos - additionalShift + (translationToLeft - translationToFrontPer),
802  zpos - translationToFrontPar,
803  vectorx, vectory, vectorz, lengthShort);
804 
805  putStrawRotatedShortRight(tmpxpos,
806  tmpypos - additionalShift + (translationToLeft + translationToFrontPer + switchSides),
807  zpos + translationToFrontPar,
808  vectorx, vectory, vectorz, lengthShort);
809 
810  tmpypos -= 2 * newradius;
811  }
812 
813  // plane 3
814  vectorx = -1 * cos((60. / 180.) * pi);
815  vectory = -1 * sin((60. / 180.) * pi);
816  vectorz = 0;
817  ypos -= availableSpace;
818  tmpxpos = xpos + (extraspacePipe + newradius) * sin((60. / 180.) * pi);
819  tmpypos = ypos - (extraspacePipe + newradius) * cos((60. / 180.) * pi);
820 
821  for(int i = 0; i < possibleStrawsPipe; i++)
822  {
823 
824  putStrawRotatedRight(tmpxpos - (translationToLeft - additionalShift) * sin(sixtyrad),
825  tmpypos + (translationToLeft - additionalShift) * cos(sixtyrad),
826  zpos, vectorx, vectory, vectorz, tubeLength);
827 
828  tmpxpos += 2 * newradius * sin((60. / 180.) * pi);
829  tmpypos -= 2 * newradius * cos((60. / 180.) * pi);
830  }
831 
832  for (int i = 0; i < extraStrawsPipe; i++)
833  {
834  double
835  lengthShort = lengthsShortPipe[i];
836 
837  double
838  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
839  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
840  switchSides = 2 * newradius * (possibleStrawsPipe + 1 + (2 * i)) * cos(skewanglerad);
841 
842  putStrawRotatedShortRight(tmpxpos - (translationToLeft - translationToFrontPer - additionalShift) * sin(sixtyrad),
843  tmpypos + (translationToLeft - translationToFrontPer - additionalShift) * cos(sixtyrad),
844  zpos - translationToFrontPar,
845  vectorx, vectory, vectorz, lengthShort);
846 
847  putStrawRotatedShortRight(tmpxpos - (translationToLeft + translationToFrontPer - additionalShift + switchSides) * sin(sixtyrad),
848  tmpypos + (translationToLeft + translationToFrontPer - additionalShift + switchSides) * cos(sixtyrad),
849  zpos + translationToFrontPar,
850  vectorx, vectory, vectorz, lengthShort);
851 
852  tmpxpos += 2 * newradius * sin((60. / 180.) * pi);
853  tmpypos -= 2 * newradius * cos((60. / 180.) * pi);
854  }
855 
856  // plane 4
857  vectorx = 1 * cos((60. / 180.) * pi);
858  vectory = -1 * sin((60. / 180.) * pi);
859  vectorz = 0;
860  xpos += availableSpace * sin((60. / 180.) * pi);
861  ypos -= availableSpace * cos((60. / 180.) * pi);
862  tmpxpos = xpos + (extraspacePipe + newradius + pipespace) * sin((60. / 180.) * pi);
863  tmpypos = ypos + (extraspacePipe + newradius + pipespace) * cos((60. / 180.) * pi);
864 
865  for(int i = 0; i < possibleStrawsPipe; i++)
866  {
867 
868  putStrawRotatedRight(tmpxpos - (translationToLeft - additionalShift) * sin(sixtyrad),
869  tmpypos - (translationToLeft - additionalShift) * cos(sixtyrad),
870  zpos, vectorx, vectory, vectorz, tubeLength);
871 
872  tmpxpos += 2 * newradius * sin((60. / 180.) * pi);
873  tmpypos += 2 * newradius * cos((60. / 180.) * pi);
874  }
875 
876  for (int i = 0; i < extraStrawsPipe; i++)
877  {
878  double
879  lengthShort = lengthsShortPipe[i];
880 
881  double
882  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
883  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
884  switchSides = 2 * newradius * (possibleStrawsPipe + 1 + (2 * i)) * cos(skewanglerad);
885 
886  putStrawRotatedShortRight(tmpxpos - (translationToLeft - additionalShift - translationToFrontPer) * sin(sixtyrad),
887  tmpypos - (translationToLeft - additionalShift - translationToFrontPer) * cos(sixtyrad),
888  zpos - translationToFrontPar,
889  vectorx, vectory, vectorz, lengthShort);
890 
891  putStrawRotatedShortRight(tmpxpos - (translationToLeft - additionalShift + translationToFrontPer + switchSides) * sin(sixtyrad),
892  tmpypos - (translationToLeft - additionalShift + translationToFrontPer + switchSides) * cos(sixtyrad),
893  zpos + translationToFrontPar,
894  vectorx, vectory, vectorz, lengthShort);
895 
896  tmpxpos += 2 * newradius * sin((60. / 180.) * pi);
897  tmpypos += 2 * newradius * cos((60. / 180.) * pi);
898  }
899 
900  // plane 5
901  vectorx = 1.;
902  vectory = 0.;
903  vectorz = 0.;
904  xpos += availableSpace * sin((60. / 180.) * pi);
905  ypos += availableSpace * cos((60. / 180.) * pi);
906  tmpxpos = xpos;
907  tmpypos = ypos + extraspace + newradius;
908 
909  for(int i = 0; i < possibleStraws; i++)
910  {
911  putStrawRotatedRight(tmpxpos, tmpypos - (translationToLeft - additionalShift), zpos, vectorx, vectory, vectorz, tubeLength);
912  tmpypos += 2 * newradius;
913  }
914 
915  for (int i = 0; i < extraStraws; i++)
916  {
917  double
918  lengthShort = lengthsShort[i];
919 
920  double
921  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
922  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
923  switchSides = 2 * newradius * (possibleStraws + 1 + (2 * i)) * cos(skewanglerad);
924 
925  putStrawRotatedShortRight(tmpxpos, tmpypos - (translationToLeft - additionalShift - translationToFrontPer),
926  zpos - translationToFrontPar,
927  vectorx, vectory, vectorz, lengthShort);
928 
929  putStrawRotatedShortRight(tmpxpos, tmpypos - (translationToLeft - additionalShift + translationToFrontPer + switchSides),
930  zpos + translationToFrontPar,
931  vectorx, vectory, vectorz, lengthShort);
932 
933  tmpypos += 2 * newradius;
934  }
935 
936  // plane 6
937  vectorx = cos((60. / 180.) * pi);
938  vectory = sin((60. / 180.) * pi);
939  vectorz = 0;
940  ypos += availableSpace;
941  tmpxpos = xpos - (extraspacePipe + newradius) * sin((60. / 180.) * pi);
942  tmpypos = ypos + (extraspacePipe + newradius) * cos((60. / 180.) * pi);
943 
944  for(int i = 0; i < possibleStrawsPipe; i++)
945  {
946 
947  putStrawRotatedRight(tmpxpos + (translationToLeft - additionalShift) * sin(sixtyrad),
948  tmpypos - (translationToLeft - additionalShift) * cos(sixtyrad),
949  zpos, vectorx, vectory, vectorz, tubeLength);
950 
951  tmpxpos -= 2 * newradius * sin((60. / 180.) * pi);
952  tmpypos += 2 * newradius * cos((60. / 180.) * pi);
953  }
954 
955  for (int i = 0; i < extraStrawsPipe; i++)
956  {
957  double
958  lengthShort = lengthsShortPipe[i];
959 
960  double
961  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
962  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
963  switchSides = 2 * newradius * (possibleStrawsPipe + 1 + (2 * i)) * cos(skewanglerad);
964 
965  putStrawRotatedShortRight(tmpxpos + (translationToLeft - additionalShift - translationToFrontPer) * sin(sixtyrad),
966  tmpypos - (translationToLeft - additionalShift - translationToFrontPer) * cos(sixtyrad),
967  zpos - translationToFrontPar,
968  vectorx, vectory, vectorz, lengthShort);
969 
970  putStrawRotatedShortRight(tmpxpos + (translationToLeft - additionalShift + translationToFrontPer + switchSides) * sin(sixtyrad),
971  tmpypos - (translationToLeft - additionalShift + translationToFrontPer + switchSides) * cos(sixtyrad),
972  zpos + translationToFrontPar,
973  vectorx, vectory, vectorz, lengthShort);
974 
975  tmpxpos -= 2 * newradius * sin((60. / 180.) * pi);
976  tmpypos += 2 * newradius * cos((60. / 180.) * pi);
977  }
978 
979 }
980 
981 
982 
983 
984 void placeSingleLayerSkewedLeft(double ringposition)
985 {
986  // calculation of some general properties of the layer layout
987  // see layerfront.ps and layertop.ps if description is not clear
988 
989  // skewed straws become ellipses in the x-y projection
990  // the short ellipse radius is still the straw radius, the
991  // long ellipse radius (within the straw plane) is given by:
992 
993  double
994  radius = tubeSeperation / 2.;
995 
996  double
997  newradius = radius / cos((skewangle / 180.) * pi);
998 
999 
1000  double
1001  xpos = 0.,
1002  ypos = ringposition / cos((30. / 180.) * pi),
1003  zpos = 0.,
1004  tmpxpos,
1005  tmpypos,
1006  vectorx,
1007  vectory,
1008  vectorz;
1009 
1010 
1011  // the available space for a layer at the radius ringposition
1012  double
1013  availableSpace = 2. * ringposition * tan((30. / 180.) * pi);
1014 
1015  // the number of skewed full-length straws that fit in the available space
1016  int
1017  possibleStraws = int((availableSpace - tubeLength * sin((skewangle / 180.) * pi)) / (newradius * 2));
1018 
1019  // the number of short straws that fit on either side of the pack of full-length straws
1020  int
1021  extraStraws = int(availableSpace / (newradius * 2)) - possibleStraws;
1022 
1023  // the amount of space that is left in the layer when considering the sum of all
1024  // full length and short straws together
1025  double
1026  extraspace = (availableSpace - (possibleStraws + extraStraws) * (newradius * 2)) / 2.;
1027 
1028  // cerr << possibleStraws << " " << extraStraws << endl;
1029 
1030  double
1031  skewanglerad = (skewangle / 180.) * pi,
1032  sixtyrad = (60. / 180.) * pi,
1033  // translation between the front end of the straw and the center
1034  translationToLeft = - 0.5 * tubeLength * tan(skewanglerad);
1035 
1036  // calculate the lenghts of the different short straws in this layer
1037  vector<double>
1038  lengthsShort;
1039 
1040  double
1041  // start at the right edge of the layer
1042  xpos2 = xpos - availableSpace * sin((60. / 180.) * pi),
1043  ypos2 = ypos - availableSpace * cos((60. / 180.) * pi),
1044  // shift to the center of the rightmost straw
1045  tmpxpos2 = xpos2,
1046  tmpypos2 = ypos2 - extraspace - newradius,
1047  // the limit on the left side of the plane beyond which no center of the short straws may be
1048  limit = ypos2 - availableSpace + newradius + safety + extraspace;
1049 
1050  // cerr << "xpos2/ypos2 " << xpos2 << " " << ypos2 << endl;
1051  // cerr << "tmpxpos2/tmpypos2 " << tmpxpos2 << " " << tmpypos2 << endl;
1052 
1053 
1054  // shift to the center of of the leftmost full-length straw
1055  for(int i = 0; i < possibleStraws; i++)
1056  {
1057  tmpypos2 -= 2 * newradius;
1058  }
1059 
1060  // calculate the maximum allowable length for this short straw and shift to the next
1061  int extracounter = 0; // how many short tubes can be really placed
1062  for (int i = 0; i < extraStraws; i++)
1063  {
1064  if((tmpypos2 - limit) < 0) break;
1065  extracounter++;
1066  double
1067  lengthShort = fabs((tmpypos2 - newradius - limit) / sin(skewangle * (pi / 180.)));
1068 
1069  if (lengthShort > tubeLength)
1070  lengthShort = tubeLength;
1071 
1072  // cerr << "length: " << lengthShort << endl;
1073  lengthsShort.push_back(lengthShort);
1074 
1075  tmpypos2 -= 2 * newradius;
1076  }
1077  extraStraws = extracounter; // update the number
1078  // =================================================
1079  // planes intersecting pipe ========================
1080  // pipe
1081  double pipespace = (panelthickness + pipeDiam/2.) * 1./(sqrt(3)/2.);
1082 
1083  // the available space for a layer at the radius ringposition
1084  double
1085  availableSpacePipe = 2. * ringposition * tan((30. / 180.) * pi) - pipespace;
1086 
1087  // the number of skewed full-length straws that fit in the available space
1088  int
1089  possibleStrawsPipe = int((availableSpacePipe - tubeLength * sin((skewangle / 180.) * pi)) / (newradius * 2));
1090 
1091  // the number of short straws that fit on either side of the pack of full-length straws
1092  int
1093  extraStrawsPipe = int(availableSpacePipe / (newradius * 2)) - possibleStrawsPipe;
1094 
1095  // the amount of space that is left in the layer when considering the sum of all
1096  // full length and short straws together
1097  double
1098  extraspacePipe = (availableSpacePipe - (possibleStrawsPipe + extraStrawsPipe) * (newradius * 2)) / 2.;
1099 
1100  // cerr << "layer with pipe " << possibleStraws << " " << extraStraws << endl;
1101 
1102  // calculate the lenghts of the different short straws in this layer
1103  vector<double>
1104  lengthsShortPipe;
1105 
1106  double
1107  // start at the right edge of the layer
1108  xpos2Pipe = xpos - availableSpacePipe * sin((60. / 180.) * pi),
1109  ypos2Pipe = ypos - availableSpacePipe * cos((60. / 180.) * pi),
1110  // shift to the center of the rightmost straw
1111  tmpxpos2Pipe = xpos2Pipe,
1112  tmpypos2Pipe = ypos2Pipe - extraspacePipe - newradius,
1113  // the limit on the left side of the plane beyond which no center of the short straws may be
1114  limitPipe = ypos2Pipe - availableSpacePipe + newradius + safety + extraspacePipe;
1115 
1116  // shift to the center of of the leftmost full-length straw
1117  for(int i = 0; i < possibleStrawsPipe; i++)
1118  {
1119  tmpypos2Pipe -= 2 * newradius;
1120  }
1121 
1122  // calculate the maximum allowable length for this short straw and shift to the next
1123  extracounter = 0; // how many short tubes can be really placed
1124  for (int i = 0; i < extraStrawsPipe; i++)
1125  {
1126  if((tmpypos2Pipe - limitPipe) < 0) break;
1127  extracounter++;
1128  double
1129  lengthShort = fabs((tmpypos2Pipe - newradius - limitPipe) / sin(skewangle * (pi / 180.)));
1130 
1131  if (lengthShort > tubeLength)
1132  lengthShort = tubeLength;
1133 
1134  lengthsShortPipe.push_back(lengthShort);
1135 
1136  tmpypos2Pipe -= 2 * newradius;
1137  }
1138  extraStrawsPipe = extracounter; // update the number
1139  // =================================================
1140  // plane 1
1141  // place full straws:
1142  vectorx = -cos((60. / 180.) * pi);
1143  vectory = sin((60. / 180.) * pi);
1144  vectorz = 0;
1145  tmpxpos = xpos - (extraspacePipe + newradius + pipespace) * sin((60. / 180.) * pi);
1146  tmpypos = ypos - (extraspacePipe + newradius + pipespace) * cos((60. / 180.) * pi);
1147 
1148  for(int i = 0; i < possibleStrawsPipe; i++)
1149  {
1150  // place a straw
1151  putStrawRotatedLeft(tmpxpos + translationToLeft * sin(sixtyrad),
1152  tmpypos + translationToLeft * cos(sixtyrad),
1153  zpos, vectorx, vectory, vectorz, tubeLength);
1154 
1155  tmpxpos -= 2 * newradius * sin((60. / 180.) * pi);
1156  tmpypos -= 2 * newradius * cos((60. / 180.) * pi);
1157  }
1158 
1159  for (int i = 0; i < extraStrawsPipe; i++)
1160  {
1161  double
1162  lengthShort = lengthsShortPipe[i];
1163 
1164  double
1165  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
1166  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
1167  // translation to the other side of the straw layer
1168  switchSides = 2 * newradius * (possibleStrawsPipe + 1 + (2 * i)) * cos(skewanglerad);
1169 
1170 
1171  putStrawRotatedShortLeft(tmpxpos + (translationToLeft + translationToFrontPer) * sin(sixtyrad),
1172  tmpypos + (translationToLeft + translationToFrontPer) * cos(sixtyrad),
1173  zpos + translationToFrontPar,
1174  vectorx, vectory, vectorz, lengthShort);
1175 
1176  putStrawRotatedShortLeft(tmpxpos + (translationToLeft - translationToFrontPer + switchSides) * sin(sixtyrad),
1177  tmpypos + (translationToLeft - translationToFrontPer + switchSides) * cos(sixtyrad),
1178  zpos - translationToFrontPar,
1179  vectorx, vectory, vectorz, lengthShort);
1180 
1181 
1182  tmpxpos -= 2 * newradius * sin((60. / 180.) * pi);
1183  tmpypos -= 2 * newradius * cos((60. / 180.) * pi);
1184  }
1185 
1186  // plane 2
1187  vectorx = -1.;
1188  vectory = 0.;
1189  vectorz = 0.;
1190  xpos -= availableSpace * sin((60. / 180.) * pi);
1191  ypos -= availableSpace * cos((60. / 180.) * pi);
1192  tmpxpos = xpos;
1193  tmpypos = ypos - extraspace - newradius;
1194 
1195  for(int i = 0; i < possibleStraws; i++)
1196  {
1197  putStrawRotatedLeft(tmpxpos, tmpypos + translationToLeft, zpos, vectorx, vectory, vectorz, tubeLength);
1198  tmpypos -= 2 * newradius;
1199  }
1200 
1201  for (int i = 0; i < extraStraws; i++)
1202  {
1203  double
1204  lengthShort = lengthsShort[i];
1205 
1206  double
1207  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
1208  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
1209  switchSides = 2 * newradius * (possibleStraws + 1 + (2 * i)) * cos(skewanglerad);
1210 
1211  putStrawRotatedShortLeft(tmpxpos,
1212  tmpypos + (translationToLeft + translationToFrontPer),
1213  zpos + translationToFrontPar,
1214  vectorx, vectory, vectorz, lengthShort);
1215 
1216  putStrawRotatedShortLeft(tmpxpos,
1217  tmpypos + (translationToLeft - translationToFrontPer + switchSides),
1218  zpos - translationToFrontPar,
1219  vectorx, vectory, vectorz, lengthShort);
1220 
1221  tmpypos -= 2 * newradius;
1222  }
1223 
1224  // plane 3
1225  vectorx = -1 * cos((60. / 180.) * pi);
1226  vectory = -1 * sin((60. / 180.) * pi);
1227  vectorz = 0;
1228  ypos -= availableSpace;
1229  tmpxpos = xpos + (extraspacePipe + newradius) * sin((60. / 180.) * pi);
1230  tmpypos = ypos - (extraspacePipe + newradius) * cos((60. / 180.) * pi);
1231 
1232  for(int i = 0; i < possibleStrawsPipe; i++)
1233  {
1234 
1235  putStrawRotatedLeft(tmpxpos - translationToLeft * sin(sixtyrad),
1236  tmpypos + translationToLeft * cos(sixtyrad),
1237  zpos, vectorx, vectory, vectorz, tubeLength);
1238 
1239  tmpxpos += 2 * newradius * sin((60. / 180.) * pi);
1240  tmpypos -= 2 * newradius * cos((60. / 180.) * pi);
1241  }
1242 
1243  for (int i = 0; i < extraStrawsPipe; i++)
1244  {
1245  double
1246  lengthShort = lengthsShortPipe[i];
1247 
1248  double
1249  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
1250  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
1251  switchSides = 2 * newradius * (possibleStrawsPipe + 1 + (2 * i)) * cos(skewanglerad);
1252 
1253  putStrawRotatedShortLeft(tmpxpos - (translationToLeft + translationToFrontPer) * sin(sixtyrad),
1254  tmpypos + (translationToLeft + translationToFrontPer) * cos(sixtyrad),
1255  zpos + translationToFrontPar,
1256  vectorx, vectory, vectorz, lengthShort);
1257 
1258  putStrawRotatedShortLeft(tmpxpos - (translationToLeft - translationToFrontPer + switchSides) * sin(sixtyrad),
1259  tmpypos + (translationToLeft - translationToFrontPer + switchSides) * cos(sixtyrad),
1260  zpos - translationToFrontPar,
1261  vectorx, vectory, vectorz, lengthShort);
1262 
1263  tmpxpos += 2 * newradius * sin((60. / 180.) * pi);
1264  tmpypos -= 2 * newradius * cos((60. / 180.) * pi);
1265  }
1266 
1267  // plane 4
1268  vectorx = 1 * cos((60. / 180.) * pi);
1269  vectory = -1 * sin((60. / 180.) * pi);
1270  vectorz = 0;
1271  xpos += availableSpace * sin((60. / 180.) * pi);
1272  ypos -= availableSpace * cos((60. / 180.) * pi);
1273  tmpxpos = xpos + (extraspacePipe + newradius + pipespace) * sin((60. / 180.) * pi);
1274  tmpypos = ypos + (extraspacePipe + newradius + pipespace) * cos((60. / 180.) * pi);
1275 
1276  for(int i = 0; i < possibleStrawsPipe; i++)
1277  {
1278 
1279  putStrawRotatedLeft(tmpxpos - translationToLeft * sin(sixtyrad),
1280  tmpypos - translationToLeft * cos(sixtyrad),
1281  zpos, vectorx, vectory, vectorz, tubeLength);
1282 
1283  tmpxpos += 2 * newradius * sin((60. / 180.) * pi);
1284  tmpypos += 2 * newradius * cos((60. / 180.) * pi);
1285  }
1286 
1287  for (int i = 0; i < extraStrawsPipe; i++)
1288  {
1289  double
1290  lengthShort = lengthsShortPipe[i];
1291 
1292  double
1293  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
1294  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
1295  switchSides = 2 * newradius * (possibleStrawsPipe + 1 + (2 * i)) * cos(skewanglerad);
1296 
1297 
1298  putStrawRotatedShortLeft(tmpxpos - (translationToLeft + translationToFrontPer) * sin(sixtyrad),
1299  tmpypos - (translationToLeft + translationToFrontPer) * cos(sixtyrad),
1300  zpos + translationToFrontPar,
1301  vectorx, vectory, vectorz, lengthShort);
1302 
1303  putStrawRotatedShortLeft(tmpxpos - (translationToLeft - translationToFrontPer + switchSides) * sin(sixtyrad),
1304  tmpypos - (translationToLeft - translationToFrontPer + switchSides) * cos(sixtyrad),
1305  zpos - translationToFrontPar,
1306  vectorx, vectory, vectorz, lengthShort);
1307 
1308  tmpxpos += 2 * newradius * sin((60. / 180.) * pi);
1309  tmpypos += 2 * newradius * cos((60. / 180.) * pi);
1310  }
1311 
1312  // plane 5
1313  vectorx = 1.;
1314  vectory = 0.;
1315  vectorz = 0.;
1316  xpos += availableSpace * sin((60. / 180.) * pi);
1317  ypos += availableSpace * cos((60. / 180.) * pi);
1318  tmpxpos = xpos;
1319  tmpypos = ypos + extraspace + newradius;
1320 
1321  for(int i = 0; i < possibleStraws; i++)
1322  {
1323 
1324  putStrawRotatedLeft(tmpxpos, tmpypos - translationToLeft, zpos, vectorx, vectory, vectorz, tubeLength);
1325 
1326  tmpypos += 2 * newradius;
1327  }
1328 
1329  for (int i = 0; i < extraStraws; i++)
1330  {
1331  double
1332  lengthShort = lengthsShort[i];
1333 
1334  double
1335  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
1336  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
1337  switchSides = 2 * newradius * (possibleStraws + 1 + (2 * i)) * cos(skewanglerad);
1338 
1339  putStrawRotatedShortLeft(tmpxpos,
1340  tmpypos - (translationToLeft + translationToFrontPer),
1341  zpos + translationToFrontPar,
1342  vectorx, vectory, vectorz, lengthShort);
1343 
1344  putStrawRotatedShortLeft(tmpxpos,
1345  tmpypos - (translationToLeft - translationToFrontPer + switchSides),
1346  zpos - translationToFrontPar,
1347  vectorx, vectory, vectorz, lengthShort);
1348 
1349  tmpypos += 2 * newradius;
1350  }
1351 
1352  // plane 6
1353  vectorx = cos((60. / 180.) * pi);
1354  vectory = sin((60. / 180.) * pi);
1355  vectorz = 0;
1356  ypos += availableSpace;
1357  tmpxpos = xpos - (extraspacePipe + newradius) * sin((60. / 180.) * pi);
1358  tmpypos = ypos + (extraspacePipe + newradius) * cos((60. / 180.) * pi);
1359 
1360  for(int i = 0; i < possibleStrawsPipe; i++)
1361  {
1362 
1363  putStrawRotatedLeft(tmpxpos + translationToLeft * sin(sixtyrad),
1364  tmpypos - translationToLeft * cos(sixtyrad),
1365  zpos, vectorx, vectory, vectorz, tubeLength);
1366 
1367  tmpxpos -= 2 * newradius * sin((60. / 180.) * pi);
1368  tmpypos += 2 * newradius * cos((60. / 180.) * pi);
1369  }
1370 
1371  for (int i = 0; i < extraStrawsPipe; i++)
1372  {
1373  double
1374  lengthShort = lengthsShortPipe[i];
1375 
1376  double
1377  translationToFrontPar = ((tubeLength - lengthShort) / 2.) * cos(skewanglerad),
1378  translationToFrontPer = ((tubeLength - lengthShort) / 2.) * sin(skewanglerad),
1379  switchSides = 2 * newradius * (possibleStrawsPipe + 1 + (2 * i)) * cos(skewanglerad);
1380 
1381  putStrawRotatedShortLeft(tmpxpos + (translationToLeft + translationToFrontPer) * sin(sixtyrad),
1382  tmpypos - (translationToLeft + translationToFrontPer) * cos(sixtyrad),
1383  zpos + translationToFrontPar,
1384  vectorx, vectory, vectorz, lengthShort);
1385 
1386  putStrawRotatedShortLeft(tmpxpos + (translationToLeft - translationToFrontPer + switchSides) * sin(sixtyrad),
1387  tmpypos - (translationToLeft - translationToFrontPer + switchSides) * cos(sixtyrad),
1388  zpos - translationToFrontPar,
1389  vectorx, vectory, vectorz, lengthShort);
1390 
1391  tmpxpos -= 2 * newradius * sin((60. / 180.) * pi);
1392  tmpypos += 2 * newradius * cos((60. / 180.) * pi);
1393  }
1394 }
1395 
1396 
1397 void makeDoubleLayerStraightExact(double &startRadius)
1398 {
1399  // close packed straight layers are only possible at specific positions
1400  // these are indicated by the integer ringteller. ringteller = 0 would be
1401  // a single straw at the origin. ringteller = 1 would be the first layer of
1402  // 6 straws closepacked around that single straw. And so on ...
1403 
1404  // look for the first ring number which is completely outside the startradius
1405  // given:
1406  double
1407  radius = tubeSeperation / 2.;
1408 
1409  int
1410  ringteller = 0;
1411 
1412  while (ringteller * sqrt(3.) * radius - radius <= startRadius)
1413  {
1414  ringteller++;
1415  }
1416  cerr << "startRadius " << startRadius << endl;
1417  // move the startRadius to the next double layer
1418  startRadius = (ringteller + 1) * sqrt(3.) * radius;
1419 
1420  // place the first layer of this double layer
1421  cerr << "positioning ring at: " << ringteller * sqrt(3.) * radius << endl;
1422  placeSingleLayerStraightExact(ringteller);
1423  // place the second layer
1424  cerr << "positioning ring at: " << (ringteller + 1) * sqrt(3.) * radius << endl;
1425  placeSingleLayerStraightExact(ringteller + 1);
1426 
1427  cerr << "ringteller " << ringteller << endl;
1428 }
1429 
1430 void makeSingleLayerStraightExact(double &startRadius)
1431 {
1432  // see comments at makeDoubleLayerStraightExact()
1433 
1434  double
1435  radius = tubeSeperation / 2.;
1436 
1437  int
1438  ringteller = 0;
1439 
1440  while (ringteller * sqrt(3.) * radius - radius + innerCoverThickness <= startRadius) // CHECK
1441  {
1442  ringteller++;
1443  }
1444 
1445  startRadius = (ringteller + 1) * sqrt(3.) * radius;
1446 
1447  cerr << "positioning ring at: " << ringteller * sqrt(3.) * radius << endl;
1448  placeSingleLayerStraightExact(ringteller);
1449  cerr << "ringteller " << ringteller << " radius " << radius << endl;
1450 
1451 }
1452 
1453 void makeDoubleLayerSkewedRight(double startRadius)
1454 {
1455  // see comments at makeDoubleLayerSkewedLeft()
1456  placeSingleLayerSkewedRight(startRadius);
1457  placeSingleLayerSkewedRight(startRadius + sqrt(3.) * (tubeSeperation / 2.));
1458 }
1459 
1460 void makeDoubleLayerSkewedLeft(double startRadius)
1461 {
1462  // place first skewed left layer at the specified radius
1463  placeSingleLayerSkewedLeft(startRadius);
1464 
1465  // the straw separation is the separation between the centers of two straws next to each other within one layer.
1466  // Close packing then dictates that separation between the layers is sqrt(3) * (separation / 2)
1467  placeSingleLayerSkewedLeft(startRadius + sqrt(3.) * (tubeSeperation / 2.));
1468 }
1469 
1470 int main()
1471 {
1472  // reset counters for the numbers of straws
1473  counter1 = 0;
1474  counter2 = 0;
1475  counter3 = 0;
1476  counter4 = 0;
1477  counter5 = 0;
1478  counter6 = 0;
1479  tubeteller = 0;
1480  axialtubeteller = 0;
1481  skewedtubeteller = 0;
1482 
1483  cout << setiosflags(ios::fixed) << setprecision(6);
1484 
1485  //------------------------------------------------------ volumes
1486 
1488  // //
1489  // STT Support: //
1490  // //
1492 
1493  if(!noSupport)
1494  {
1495  // container for all tubes!
1496  // for X > 0
1497  cout << sttassembly << endl;
1498  cout << cave << endl;
1499  cout << "ASSEMBLY" << endl;
1500  cout << air << endl;
1501  counter4++;
1502  writetrans(0., 0., sttCenterZ);
1503  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
1504 
1505  // cylinder inner x > 0
1506  writename(innerCylinder, 1, 1);
1507  writemother(sttassembly, 1, 1);
1508  writemedium(AlBe);
1510  writetrans(0., 0., 0.);
1511  writerot(0., 1., 0., -1., 0., 0., 0., 0., 1.);
1512  // writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
1513 
1514  // cylinder outer x > 0
1515  writename(outerCylinder, 1, 1);
1516  writemother(sttassembly, 1, 1);
1517  writemedium(AlBe);
1519  writetrans(0., 0., 0.);
1520  writerot(0., 1., 0., -1., 0., 0., 0., 0., 1.);
1521  // writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
1522 
1523  // cylinder inner x < 0
1524  writename(innerCylinder, 1, 0);
1525  writemother(sttassembly, 1, 1);
1526  writemedium(AlBe);
1528  writetrans(0., 0., 0.);
1529  writerot(0., -1., 0., 1., 0., 0., 0., 0., 1.);
1530  // writerot(1., 0., 0., 0., -1., 0., 0., 0., 1.);
1531 
1532  // cylinder outer x < 0
1533  writename(outerCylinder, 1, 0);
1534  writemother(sttassembly, 1, 1);
1535  writemedium(AlBe);
1537  writetrans(0., 0., 0.);
1538  writerot(0., -1., 0., 1., 0., 0., 0., 0., 1.);
1539  // writerot(1., 0., 0., 0., -1., 0., 0., 0., 1.);
1540 
1541  // around the pipe CHECK why (pipeDiam/2.) "-" panelthickness/2?? -----------------
1542  // panel up x > 0
1544  writetrans((pipeDiam/2.) + panelthickness/2., (outerDiam/2. + innerDiam/2.)/2., 0.);
1545  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
1546  // panel up x < 0
1548  writetrans(-((pipeDiam/2.) + panelthickness/2.), (outerDiam/2. + innerDiam/2.)/2., 0.);
1549  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
1550  // panel down x > 0
1551  writepanel(panel3, false, panelthickness, (outerDiam/2. - innerDiam/2.), tubeLength, 1);
1552  writetrans((pipeDiam/2.) + panelthickness/2., -(outerDiam/2. + innerDiam/2.)/2., 0.);
1553  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
1554  // panel down x < 0
1555  writepanel(panel4, false, panelthickness, (outerDiam/2. - innerDiam/2.), tubeLength, 2);
1556  writetrans(-((pipeDiam/2.) + panelthickness/2.), -(outerDiam/2. + innerDiam/2.)/2., 0.);
1557  writerot(1., 0., 0., 0., 1., 0., 0., 0., 1.);
1558 
1559  }
1560 
1562  // //
1563  // Layers Geometry and //
1564  // Placement of the logical tubes : //
1565  // //
1567 
1568  double
1569  startRadius = innerCoverThickness + innerDiam / 2.; // CHECK added CoverThickness // CHECK
1570 
1571  bool
1572  wintzdesign = true;
1573 
1574  if (wintzdesign)
1575  {
1576 
1577  // NEW ONE!
1578  // 4 double layers of straight straws
1579  makeDoubleLayerStraightExact(startRadius);
1580  makeDoubleLayerStraightExact(startRadius);
1581  makeDoubleLayerStraightExact(startRadius);
1582  makeDoubleLayerStraightExact(startRadius);
1583 
1584  cerr << "start radius before skewed " << startRadius << endl;
1585  // startRadius = (ringteller + 1) * sqrt(3.) * radius;
1586 
1587  makeDoubleLayerSkewedLeft((26 * sqrt(3.) + 2) * (tubeSeperation / 2.));
1588  makeDoubleLayerSkewedRight((27 * sqrt(3.) + 4) * (tubeSeperation / 2.) + 0.244);
1589  makeDoubleLayerSkewedLeft((28 * sqrt(3.) + 6) * (tubeSeperation / 2.) + 2*0.244);
1590  makeDoubleLayerSkewedRight((29 * sqrt(3.) + 8) * (tubeSeperation / 2.) + 3*0.244);
1591 
1592 
1593  cerr << "start radius after " << startRadius << endl;
1594  startRadius += 10 * tubeSeperation * sqrt(3.) /2.;
1595 
1596  cerr << "start radius " << startRadius << endl;
1597 
1598  // then two double layers of straight straws
1599  makeDoubleLayerStraightExact(startRadius);
1600  makeDoubleLayerStraightExact(startRadius);
1601 
1602  }
1603  else
1604  {
1605  makeDoubleLayerStraightExact(startRadius);
1606  startRadius += tubeSeperation;
1607  makeDoubleLayerSkewedLeft(startRadius);
1608  startRadius += tubeSeperation + sqrt(3.) * (tubeSeperation / 2.);
1609  makeDoubleLayerSkewedRight(startRadius);
1610  startRadius += (tubeSeperation / 2.) + sqrt(3.) * (tubeSeperation / 2.);
1611  makeDoubleLayerStraightExact(startRadius);
1612  startRadius += tubeSeperation;
1613  makeDoubleLayerSkewedLeft(startRadius);
1614  startRadius += tubeSeperation + sqrt(3.) * (tubeSeperation / 2.);
1615  makeDoubleLayerSkewedRight(startRadius);
1616  startRadius += (tubeSeperation / 2.) + sqrt(3.) * (tubeSeperation / 2.);
1617  makeDoubleLayerStraightExact(startRadius);
1618  startRadius += tubeSeperation;
1619  makeDoubleLayerSkewedLeft(startRadius);
1620  startRadius += tubeSeperation + sqrt(3.) * (tubeSeperation / 2.);
1621  makeDoubleLayerSkewedRight(startRadius);
1622  startRadius += (tubeSeperation / 2.) + sqrt(3.) * (tubeSeperation / 2.);
1623  makeDoubleLayerStraightExact(startRadius);
1624  }
1625 
1626  // fill the rest of the volume with straight straws
1627  // straws outside the outerDiam will not be placed
1628  bool
1629  fillup = true;
1630 
1631  if (fillup)
1632  {
1633  makeDoubleLayerStraightExact(startRadius);
1634  makeDoubleLayerStraightExact(startRadius);
1635  makeDoubleLayerStraightExact(startRadius);
1636  makeDoubleLayerStraightExact(startRadius);
1637  makeDoubleLayerStraightExact(startRadius);
1638  makeDoubleLayerStraightExact(startRadius);
1639  makeDoubleLayerStraightExact(startRadius);
1640  makeDoubleLayerStraightExact(startRadius);
1641  makeDoubleLayerStraightExact(startRadius);
1642  makeDoubleLayerStraightExact(startRadius);
1643  }
1644 
1645  // output of all counters
1646  cerr << "tubes: " << tubeteller << endl;
1647  cerr << "axial tubes: " << axialtubeteller << endl;
1648  cerr << "skewed tubes full length: " << skewedtubeteller << endl;
1649  cerr << "skewed tubes short: " << shortskewedtubeteller << endl;
1650  cerr << "minimum radius: " << minimumradius << endl;
1651  cerr << "maximum radius: " << maximumradius << endl;
1652 
1653  return 0;
1654 }
1655 
1656 
1657 
#define panel2
Definition: createSTT.C:67
static int axialtubeteller
Definition: createSTT.C:79
int main(int argc, char **argv)
#define tubeLength
Definition: createSTT.C:42
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
bool putStraw(double posX, double posY, double posZ)
Definition: createSTT.C:171
static int counter4
Definition: createSTT.C:30
#define innerCylinder
Definition: createSTT.C:64
#define HeMixture
Definition: createSTT.C:75
static int tubeteller
Definition: createSTT.C:78
Int_t i
Definition: run_full.C:25
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
#define panelthickness
Definition: createSTT.C:54
#define tubeInnerDiam
Definition: createSTT.C:37
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
bool putStrawRotatedShortLeft(double posX, double posY, double posZ, double xvector, double yvector, double zvector, double length)
Definition: createSTT.C:295
#define panel3
Definition: createSTT.C:68
#define sttCenterZ
Definition: createSTT.C:49
#define pi
Definition: createSTT.C:60
#define innerDiam
Definition: createSTT.C:36
static double maximumradius
Definition: createSTT.C:82
void writetrans(double x, double y, double z)
Definition: createSTT.C:155
void makeDoubleLayerStraightExact(double &startRadius)
Definition: createSTT.C:1397
static int counter2
Definition: createSTT.C:28
void writerot(double x00, double x01, double x02, double x10, double x11, double x12, double x20, double x21, double x22)
Definition: createSTT.C:161
void writename(char const *name, bool support=false, bool leftside=false)
Definition: createSTT.C:86
#define W
Definition: createSTT.C:76
void writemedium(char *name)
Definition: createSTT.C:104
#define cave
Definition: createSTT.C:62
bool putStrawRotatedRight(double posX, double posY, double posZ, double xvector, double yvector, double zvector, double length)
Definition: createSTT.C:474
bool putStrawRotatedShortRight(double posX, double posY, double posZ, double xvector, double yvector, double zvector, double length)
Definition: createSTT.C:352
#define outerCoverThickness
Definition: createSTT.C:52
static int counter5
Definition: createSTT.C:31
#define noSupport
Definition: createSTT.C:56
const Double_t zpos
#define innerCoverThickness
Definition: createSTT.C:51
#define outerDiam
Definition: createSTT.C:35
void plotrotright(double x, double y, double z)
Definition: createSTT.C:257
#define pipeDiam
Definition: createSTT.C:55
void makeSingleLayerStraightExact(double &startRadius)
Definition: createSTT.C:1430
void placeSingleLayerStraightExact(double ringteller)
Definition: createSTT.C:542
Double_t z
static int skewedtubeteller
Definition: createSTT.C:80
#define wireDiam
Definition: createSTT.C:41
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
void makeDoubleLayerSkewedRight(double startRadius)
Definition: createSTT.C:1453
#define mylar
Definition: createSTT.C:74
void placeSingleLayerSkewedRight(double ringposition)
Definition: createSTT.C:593
#define tubeSeperation
Definition: createSTT_150.C:38
bool putStrawRotatedLeft(double posX, double posY, double posZ, double xvector, double yvector, double zvector, double length)
Definition: createSTT.C:410
TString name
void writehalftube(double inner, double outer, double length)
Definition: createSTT.C:118
static double minimumradius
Definition: createSTT.C:83
#define skewangle
Definition: createSTT.C:58
void placeSingleLayerSkewedLeft(double ringposition)
Definition: createSTT.C:984
void makeDoubleLayerSkewedLeft(double startRadius)
Definition: createSTT.C:1460
Double_t x
#define panel1
Definition: createSTT.C:66
void writetube(double inner, double outer, double length)
Definition: createSTT.C:110
static int shortskewedtubeteller
Definition: createSTT.C:81
void plotrotleft(double x, double y, double z)
Definition: createSTT.C:278
#define safety
Definition: createSTT.C:44
#define sttassembly
Definition: createSTT.C:63
static int counter3
Definition: createSTT.C:29
void writepanel(char const *name, bool firstone, double xthick, double ythick, double length, int side)
Definition: createSTT.C:127
Double_t y
#define panel4
Definition: createSTT.C:69
double limit
Definition: dedx_bands.C:18
void writemother(char const *name, bool original=true, bool support=false)
Definition: createSTT.C:94
#define AlBe
Definition: createSTT.C:73
static int counter1
Definition: createSTT.C:27
#define air
Definition: createSTT.C:72
#define tubeOuterDiam
Definition: createSTT.C:38
#define outerCylinder
Definition: createSTT.C:65
static int counter6
Definition: createSTT.C:32