Fork me on GitHub

source: svn/trunk/external/Hector/H_RectEllipticAperture.cc@ 1360

Last change on this file since 1360 was 1360, checked in by Pavel Demin, 11 years ago

add Hector module

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 4.5 KB
Line 
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * *
2 * *
3* --<--<-- A fast simulator --<--<-- *
4* / --<--<-- of particle --<--<-- *
5* ----HECTOR----< *
6* \ -->-->-- transport through -->-->-- *
7* -->-->-- generic beamlines -->-->-- *
8* *
9* JINST 2:P09005 (2007) *
10* X Rouby, J de Favereau, K Piotrzkowski (CP3) *
11* http://www.fynu.ucl.ac.be/hector.html *
12* *
13* Center for Cosmology, Particle Physics and Phenomenology *
14* Universite catholique de Louvain *
15* Louvain-la-Neuve, Belgium *
16 * *
17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */
18
19/// \file H_RectEllipticAperture.cc
20/// \brief Defines the Rect-Elliptic aperture of beamline elements.
21
22// C++ #includes
23#include <iostream>
24
25// C #includes
26#include <cmath> // needed for fabs
27
28// ROOT #includes
29#include "TPolyLine.h"
30
31// local #includes
32#include "H_RectEllipticAperture.h"
33using namespace std;
34
35H_RectEllipticAperture::H_RectEllipticAperture(const float l, const float h, const float L, const float H, const float posx, const float posy) :H_Aperture(RECTELLIPSE,((l==0)?L:l),((h==0)?H:h),L,H,posx,posy) {
36 /// @param l, h, L, H are the geometrical parameters of the rect-ellipse shape
37 /// @param posx, posy defines the (x,y) of the center of the shape
38}
39
40H_RectEllipticAperture* H_RectEllipticAperture::clone() const {
41 return new H_RectEllipticAperture(x1,x2,x3,x4,fx,fy);
42}
43
44
45TPolyLine * rectellipse(const float a_e = 2, const float b_e = 1, const float a_r = 1, const float b_r = 2, const float center_x = 0, const float center_y =0) {
46 const int n = 20; // number of points per segment
47 const int N = 4*n; // there are 4 segments
48 float x[N+1], y[N+1];
49
50 if(a_e>a_r) {
51 // a rectellipse has 4 segments
52 // 1) upper one
53 for (int i=0; i<n; i++) {
54 x[i] = -a_r + i*(2*a_r)/(float)n;
55 y[i] = b_e * sqrt(1-pow(x[i]/a_e,2));
56 }
57
58 // 2) right vertical segment
59 // upper right corner
60 const float y2 = b_e * sqrt(1-pow(a_r/a_e,2));
61 // lower right corner
62 const float y3 = -b_e * sqrt(1-pow(a_r/a_e,2));
63 for (int i=n; i<2*n; i++) {
64 x[i] = a_r;
65 y[i] = y2 - (i-n)*(2*y2)/(float)n;
66 }
67
68 // 3) lower side
69 for (int i=2*n; i<3*n; i++) {
70 x[i] = a_r - (i-2*n)*(2*a_r)/(float)n;
71 y[i] = -b_e * sqrt(1-pow(x[i]/a_e,2));
72 }
73
74 // 4) left vertical segment
75 // lower left corner
76 const float y4 = y3;
77 for (int i=3*n; i<4*n; i++) {
78 x[i] = -a_r;
79 y[i] = y4 + (i-3*n)*(2*y2)/(float)n;
80 }
81 } else {
82 // 1) upper one : flat
83 const float x1 = -a_e * sqrt(1-pow(b_r/b_e,2));
84 const float x2 = -x1;
85 for (int i=0; i<n; i++) {
86 y[i] = b_r;
87 x[i] = x1 + i * (x2-x1)/(float)n;
88 }
89
90 // 2) right curved border
91 for (int i=n; i<2*n; i++) {
92 y[i] = b_r - (i-n) * (2*b_r)/(float)n;
93 x[i] = a_e * sqrt(1-pow(y[i]/b_e,2));
94 }
95
96 // 3) lower side : flat
97 for (int i=2*n; i<3*n; i++) {
98 y[i] = -b_r;
99 x[i] = x2 - (i-2*n) * (2*x2)/(float)n;
100 }
101
102 // 4) left curved border
103 for (int i=3*n; i<4*n; i++) {
104 y[i] = -b_r + (i-3*n) * (2*b_r)/(float)n;
105 x[i] = -a_e * sqrt(1-pow(y[i]/b_e,2));
106 }
107 }
108
109 // closing the polyline
110 x[N] = x[0];
111 y[N] = y[0];
112
113 // shifting the center
114 for (int i=0; i<N+1; i++) {
115 x[i] += center_x;
116 y[i] += center_y;
117 }
118
119 return new TPolyLine(N+1,x,y);
120}
121
122
123void H_RectEllipticAperture::draw(const float scale) const {
124 TPolyLine * re = rectellipse(x3*scale, x4*scale, x1*scale, x2*scale, fx*scale, fy*scale);
125 re->SetLineColor(39);
126 re->SetLineWidth(2);
127 re->Draw("l");
128 return;
129}
130
131bool H_RectEllipticAperture::isInside(const float x, const float y) const {
132 return((fabs(fx-x)<x1)&&(fabs(fy-y)<x2)&&(((x-fx)/x3)*((x-fx)/x3)+((y-fy)/x4)*((y-fy)/x4)<1));
133}
134
135std::ostream& operator<< (std::ostream& os, const H_RectEllipticAperture& ap) {
136 os << "Aperture shape:" << ap.aptypestring << ", parameters " << ap.x1 <<", "<< ap.x2 <<", "<< ap.x3 <<", "<< ap.x4 << endl;
137 os << " \t Center : " << ap.fx <<", "<< ap.fy <<endl;
138 return os;
139}
Note: See TracBrowser for help on using the repository browser.