[3] | 1 | /*
|
---|
| 2 | ---- Hector the simulator ----
|
---|
| 3 | A fast simulator of particles through generic beamlines.
|
---|
| 4 | J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be
|
---|
| 5 |
|
---|
| 6 | http://www.fynu.ucl.ac.be/hector.html
|
---|
| 7 |
|
---|
| 8 | Centre de Physique des Particules et de Phénoménologie (CP3)
|
---|
| 9 | Université Catholique de Louvain (UCL)
|
---|
| 10 | */
|
---|
| 11 |
|
---|
| 12 | /// \file H_OpticalElement.cc
|
---|
| 13 | /// \brief Class aiming at describing any beam optical element.
|
---|
| 14 |
|
---|
| 15 | // c++ #includes
|
---|
| 16 | #include <iostream>
|
---|
| 17 | #include <string>
|
---|
| 18 |
|
---|
| 19 | // ROOT #includes
|
---|
| 20 | #include "TPaveLabel.h"
|
---|
| 21 |
|
---|
| 22 | // local #includes
|
---|
| 23 | #include "H_Parameters.h"
|
---|
| 24 | #include "H_TransportMatrices.h"
|
---|
| 25 | #include "H_Aperture.h"
|
---|
| 26 | #include "H_OpticalElement.h"
|
---|
| 27 | using namespace std;
|
---|
| 28 |
|
---|
| 29 | /// called by the constructors
|
---|
| 30 | void H_OpticalElement::init(const string nameE, const int typeE, const double s, const double k, const double l) {
|
---|
| 31 | // this is called by the constructors
|
---|
| 32 | // must be in public section !
|
---|
| 33 | // xpos and ypos are vectors with n point. They define the aperture shape of the optical element.
|
---|
| 34 | name = nameE;
|
---|
| 35 | fs = s;
|
---|
| 36 | fk = k;
|
---|
| 37 | xpos = 0;
|
---|
| 38 | ypos = 0;
|
---|
| 39 | txpos = 0;
|
---|
| 40 | typos = 0;
|
---|
| 41 | element_length = l;
|
---|
| 42 | type = typeE;
|
---|
| 43 | element_mat.ResizeTo(MDIM,MDIM);
|
---|
| 44 | element_mat = driftmat(l);
|
---|
| 45 |
|
---|
| 46 | // do NOT use setAperture for the initialisation ! there are protections there
|
---|
| 47 |
|
---|
| 48 | if(element_length<0) { if(VERBOSE) cout<<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl; }
|
---|
| 49 | if(element_length==0) { if(VERBOSE) cout<<"<H_OpticalElement> WARNING : 0-length element ! (" << name << ") " << " at " << fs << endl; }
|
---|
| 50 | betax =0;
|
---|
| 51 | betay =0;
|
---|
| 52 | }
|
---|
| 53 |
|
---|
| 54 | H_OpticalElement::H_OpticalElement() : element_aperture(new H_Aperture()) {
|
---|
| 55 | init("",DRIFT,0.,0.,0.1);
|
---|
| 56 | }
|
---|
| 57 |
|
---|
| 58 | H_OpticalElement::H_OpticalElement(const string nameE, const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) {
|
---|
| 59 | init(nameE,typeE,s,k,l);
|
---|
| 60 | }
|
---|
| 61 |
|
---|
| 62 | H_OpticalElement::H_OpticalElement(const string nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) : element_aperture(the_app->clone()) {
|
---|
| 63 | init(nameE,typeE,s,k,l);
|
---|
| 64 | }
|
---|
| 65 |
|
---|
| 66 | H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l, H_Aperture* the_app) : element_aperture(the_app->clone()) {
|
---|
| 67 | init("",typeE,s,k,l);
|
---|
| 68 | }
|
---|
| 69 |
|
---|
| 70 | H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) {
|
---|
| 71 | init("",typeE,s,k,l);
|
---|
| 72 | }
|
---|
| 73 |
|
---|
| 74 | H_OpticalElement::H_OpticalElement(const H_OpticalElement& el) {
|
---|
| 75 | fs = el.fs;
|
---|
| 76 | element_length = el.element_length;
|
---|
| 77 | fk = el.fk;
|
---|
| 78 | xpos = el.xpos;
|
---|
| 79 | ypos = el.ypos;
|
---|
| 80 | txpos = el.txpos;
|
---|
| 81 | typos = el.typos;
|
---|
| 82 | betax = el.betax;
|
---|
| 83 | betay = el.betay;
|
---|
| 84 | type = el.type;
|
---|
| 85 | name = el.name;
|
---|
| 86 | typestring = el.typestring;
|
---|
| 87 | element_mat.ResizeTo(MDIM,MDIM);
|
---|
| 88 | element_mat = el.element_mat;
|
---|
| 89 | element_aperture = el.element_aperture->clone();
|
---|
| 90 | }
|
---|
| 91 |
|
---|
| 92 | H_OpticalElement& H_OpticalElement::operator=(const H_OpticalElement& el) {
|
---|
| 93 | if(this==&el) return *this;
|
---|
| 94 | fs = el.fs;
|
---|
| 95 | element_length = el.element_length;
|
---|
| 96 | fk = el.fk;
|
---|
| 97 | xpos = el.xpos;
|
---|
| 98 | ypos = el.ypos;
|
---|
| 99 | txpos = el.txpos;
|
---|
| 100 | typos = el.typos;
|
---|
| 101 | betax = el.betax;
|
---|
| 102 | betay = el.betay;
|
---|
| 103 | type = el.type;
|
---|
| 104 | name = el.name;
|
---|
| 105 | typestring = el.typestring;
|
---|
| 106 | element_mat.ResizeTo(MDIM,MDIM);
|
---|
| 107 | element_mat = el.element_mat;
|
---|
| 108 | element_aperture = el.element_aperture->clone();
|
---|
| 109 | return *this;
|
---|
| 110 | }
|
---|
| 111 |
|
---|
| 112 | void H_OpticalElement::setAperture(const H_Aperture* ap) {
|
---|
| 113 | // do NOT use setAperture in your constructor, as element_aperture is not initialized
|
---|
| 114 | // this function do not take into account ap if ap=0
|
---|
| 115 | // do nothing if element_mat = ap
|
---|
| 116 | if (!ap) {cout << "<H_OpticalElement> Trying to set an empty pointer for the aperture ! Nothing done.\n"; return;}
|
---|
| 117 | if (element_aperture != ap) {
|
---|
| 118 | delete element_aperture;
|
---|
| 119 | element_aperture = ap->clone();
|
---|
| 120 | }
|
---|
| 121 | return;
|
---|
| 122 | }
|
---|
| 123 |
|
---|
| 124 | void H_OpticalElement::printProperties() const {
|
---|
| 125 | cout << typestring;
|
---|
| 126 | cout << name;
|
---|
| 127 | cout <<"\t at s = " << fs;
|
---|
| 128 | cout <<"\t length = "<< element_length;
|
---|
| 129 | if(fk!=0) cout <<"\t strength = " << fk;
|
---|
| 130 | if(element_aperture->getType()!=NONE) {
|
---|
| 131 | cout <<"\t aperture type = " << element_aperture->getTypeString();
|
---|
| 132 | element_aperture->printProperties();
|
---|
| 133 | }
|
---|
| 134 |
|
---|
| 135 | cout<<endl;
|
---|
| 136 | if(element_length<0) { if(VERBOSE) cout<<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl; }
|
---|
| 137 | if(element_length==0) { if(VERBOSE) cout<<"<H_OpticalElement> WARNING : 0-length "<< typestring << " !" << endl; }
|
---|
| 138 |
|
---|
| 139 | return;
|
---|
| 140 | }
|
---|
| 141 |
|
---|
| 142 | void H_OpticalElement::showMatrix() const {
|
---|
| 143 | printMatrix(element_mat);
|
---|
| 144 | return;
|
---|
| 145 | }
|
---|
| 146 |
|
---|
| 147 | void H_OpticalElement::drawAperture() const {
|
---|
| 148 | element_aperture->draw();
|
---|
| 149 | return;
|
---|
| 150 | }
|
---|
| 151 |
|
---|
| 152 | TMatrix H_OpticalElement::getMatrix() {
|
---|
| 153 | setMatrix(0,MP,1);
|
---|
| 154 | return element_mat;
|
---|
| 155 | }
|
---|
| 156 |
|
---|
| 157 | TMatrix H_OpticalElement::getMatrix(const float eloss, const float p_mass, const float p_charge) {
|
---|
| 158 | setMatrix(eloss,p_mass,p_charge);
|
---|
| 159 | return element_mat;
|
---|
| 160 | }
|
---|
| 161 |
|
---|
| 162 | void H_OpticalElement::draw(const float meight, const float height) const{
|
---|
| 163 | /// @param meight is the minimal extend of the graph
|
---|
| 164 | /// @param height is the maximal extend of the graph
|
---|
| 165 | float x1 = getS();
|
---|
| 166 | float x2 = getS() + getLength();
|
---|
| 167 | float y1 = meight;
|
---|
| 168 | float y2 = height;
|
---|
| 169 | TPaveLabel* cur_box = new TPaveLabel(x1,y1,x2,y2,"");
|
---|
| 170 | cur_box->SetBorderSize(1);
|
---|
| 171 | cur_box->SetFillStyle(1001);
|
---|
| 172 | cur_box->SetFillColor((int)getType());
|
---|
| 173 | cur_box->Draw();
|
---|
| 174 | }
|
---|
| 175 |
|
---|
| 176 | TVectorD H_OpticalElement::getHitPosition(TVectorD init_pos, double energy_loss, double mp, double qp) {
|
---|
| 177 | if(!element_length) {
|
---|
| 178 | // cout<<"O-length element ("<<getName()<<"), should not appear here !"<<endl;
|
---|
| 179 | return init_pos;
|
---|
| 180 | }
|
---|
| 181 | // some declarations
|
---|
| 182 | bool inside = false;
|
---|
| 183 | double vec[MDIM] = {init_pos[INDEX_X]/URAD,tan(init_pos[INDEX_TX]/URAD),init_pos[INDEX_Y]/URAD,tan(init_pos[INDEX_TY]/URAD),-energy_loss,1};
|
---|
| 184 | TMatrixD mat_init(1,MDIM,vec);
|
---|
| 185 | TMatrixD mat_min(1,MDIM,vec);
|
---|
| 186 | TMatrixD mat_max(1,MDIM,vec);
|
---|
| 187 | TMatrixD mat_stop(1,MDIM,vec);
|
---|
| 188 | H_OpticalElement* temp_el = clone();
|
---|
| 189 | // initialasing boundaries
|
---|
| 190 | double min_pos = 0;
|
---|
| 191 | double max_pos = element_length/2.;
|
---|
| 192 | double max_old = element_length/2.;
|
---|
| 193 | // number of iterations
|
---|
| 194 | // (idea : fix precision instead of number of iterations + add security)
|
---|
| 195 | // (idea : interpolate between max and min and give the error)
|
---|
| 196 | const int N = 10;
|
---|
| 197 | // starting search loop
|
---|
| 198 | for(int i = 0; i < N; i++) {
|
---|
| 199 | // fixing position to be investigated
|
---|
| 200 | temp_el->setLength(max_pos);
|
---|
| 201 | // initialising the vector at the initial vector + possible shift/tilt
|
---|
| 202 | mat_max[0][0] = mat_init[0][0] - temp_el->getX();
|
---|
| 203 | mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
|
---|
| 204 | mat_max[0][2] = mat_init[0][2] - temp_el->getY();
|
---|
| 205 | mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
|
---|
| 206 | // propagating
|
---|
| 207 | mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
|
---|
| 208 | // compensating for the previously stated shifts/tilts
|
---|
| 209 | mat_max[0][0] = mat_max[0][0] + temp_el->getX();
|
---|
| 210 | mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
|
---|
| 211 | mat_max[0][2] = mat_max[0][2] + temp_el->getY();
|
---|
| 212 | mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
|
---|
| 213 | // fixing new boundaries
|
---|
| 214 | if(temp_el->isInside(mat_max.GetMatrixArray()[0]*URAD,mat_max.GetMatrixArray()[2]*URAD)) {
|
---|
| 215 | max_old = max_pos;
|
---|
| 216 | max_pos = max_pos + (max_pos - min_pos)/2.;
|
---|
| 217 | min_pos = max_old;
|
---|
| 218 | inside = true;
|
---|
| 219 | } else {
|
---|
| 220 | max_pos = min_pos + (max_pos - min_pos)/2.;
|
---|
| 221 | inside = false;
|
---|
| 222 | }
|
---|
| 223 | // end of loop
|
---|
| 224 | }
|
---|
| 225 | // if it passes at the last iteration, choosing the other range²
|
---|
| 226 | if(inside) min_pos = max_old;
|
---|
| 227 | // here the interpolation method : now we are sure that the intercept is between min_pos and max_pos
|
---|
| 228 | // getting vector at min_pos (for first boundary) :
|
---|
| 229 | bool precision_estimate = false;
|
---|
| 230 | if(precision_estimate) {
|
---|
| 231 | temp_el->setLength(min_pos);
|
---|
| 232 | mat_min[0][0] = mat_init[0][0] - temp_el->getX();
|
---|
| 233 | mat_min[0][1] = mat_init[0][1] - tan(temp_el->getTX());
|
---|
| 234 | mat_min[0][2] = mat_init[0][2] - temp_el->getY();
|
---|
| 235 | mat_min[0][3] = mat_init[0][3] - tan(temp_el->getTY());
|
---|
| 236 | mat_min *= temp_el->getMatrix(energy_loss,mp,qp);
|
---|
| 237 | mat_min[0][0] = mat_min[0][0] + temp_el->getX();
|
---|
| 238 | mat_min[0][1] = mat_min[0][1] + tan(temp_el->getTX());
|
---|
| 239 | mat_min[0][2] = mat_min[0][2] + temp_el->getY();
|
---|
| 240 | mat_min[0][3] = mat_min[0][3] + tan(temp_el->getTY());
|
---|
| 241 | mat_min[0][4] = min_pos + init_pos[4];
|
---|
| 242 | // getting vector at max_pos (for second boundary) :
|
---|
| 243 | temp_el->setLength(max_pos);
|
---|
| 244 | mat_max[0][0] = mat_init[0][0] - temp_el->getX();
|
---|
| 245 | mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
|
---|
| 246 | mat_max[0][2] = mat_init[0][2] - temp_el->getY();
|
---|
| 247 | mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
|
---|
| 248 | mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
|
---|
| 249 | mat_max[0][0] = mat_max[0][0] + temp_el->getX();
|
---|
| 250 | mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
|
---|
| 251 | mat_max[0][2] = mat_max[0][2] + temp_el->getY();
|
---|
| 252 | mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
|
---|
| 253 | mat_max[0][4] = max_pos + init_pos[4];
|
---|
| 254 | }
|
---|
| 255 | // getting vector in the middle (for estimate) :
|
---|
| 256 | temp_el->setLength((max_pos+min_pos)/2.);
|
---|
| 257 | mat_stop[0][0] = mat_init[0][0] - temp_el->getX();
|
---|
| 258 | mat_stop[0][1] = mat_init[0][1] - tan(temp_el->getTX());
|
---|
| 259 | mat_stop[0][2] = mat_init[0][2] - temp_el->getY();
|
---|
| 260 | mat_stop[0][3] = mat_init[0][3] - tan(temp_el->getTY());
|
---|
| 261 | mat_stop *= temp_el->getMatrix(energy_loss,mp,qp);
|
---|
| 262 | mat_stop[0][0] = mat_stop[0][0] + temp_el->getX();
|
---|
| 263 | mat_stop[0][1] = mat_stop[0][1] + tan(temp_el->getTX());
|
---|
| 264 | mat_stop[0][2] = mat_stop[0][2] + temp_el->getY();
|
---|
| 265 | mat_stop[0][3] = mat_stop[0][3] + tan(temp_el->getTY());
|
---|
| 266 | mat_stop[0][4] = (max_pos+min_pos)/2. + init_pos[4];
|
---|
| 267 |
|
---|
| 268 | double xys[LENGTH_VEC];
|
---|
| 269 | xys[INDEX_X]= mat_stop[0][0]*URAD;
|
---|
| 270 | xys[INDEX_TX]= atan(mat_stop[0][1])*URAD;
|
---|
| 271 | xys[INDEX_Y]= mat_stop[0][2]*URAD;
|
---|
| 272 | xys[INDEX_TY]= atan(mat_stop[0][3])*URAD;
|
---|
| 273 | xys[INDEX_S]= mat_stop[0][4] ;
|
---|
| 274 | TVectorD temp_vec(LENGTH_VEC,xys);
|
---|
| 275 |
|
---|
| 276 | if(precision_estimate) {
|
---|
| 277 | cout<<"--- Results and precision estimates ---"<<endl;
|
---|
| 278 | cout<<"\t Stopping element : "<<getName()<<endl;
|
---|
| 279 | cout<<"\t hit point s : "<<mat_stop[0][4]<<" m +- "<<(mat_max[0][4]-mat_stop[0][4])*1000.<<" mm"<<endl;
|
---|
| 280 | cout<<"\t hit point x : "<<mat_stop[0][0]*URAD;
|
---|
| 281 | cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][0])?(mat_max[0][0]-mat_stop[0][0])*URAD:(mat_min[0][0]-mat_stop[0][0])*URAD));
|
---|
| 282 | cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][0])?(mat_stop[0][0]-mat_min[0][0])*URAD:(mat_stop[0][0]-mat_max[0][0])*URAD))<<" µm"<<endl;
|
---|
| 283 | cout<<"\t hit point y : "<<mat_stop[0][2]*URAD;
|
---|
| 284 | cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][2])?(mat_max[0][2]-mat_stop[0][2])*URAD:(mat_min[0][2]-mat_stop[0][2])*URAD));
|
---|
| 285 | cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][2])?(mat_stop[0][2]-mat_min[0][2])*URAD:(mat_stop[0][2]-mat_max[0][2])*URAD))<<" µm"<<endl;
|
---|
| 286 | cout<<"\t hit point tx : "<<mat_stop[0][1]*URAD;
|
---|
| 287 | cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][1])?(mat_max[0][1]-mat_stop[0][1])*URAD:(mat_min[0][1]-mat_stop[0][1])*URAD));
|
---|
| 288 | cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][1])?(mat_stop[0][1]-mat_min[0][1])*URAD:(mat_stop[0][1]-mat_max[0][1])*URAD))<<" µrad"<<endl;
|
---|
| 289 | cout<<"\t hit point ty : "<<mat_stop[0][3]*URAD;
|
---|
| 290 | cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][3])?(mat_max[0][3]-mat_stop[0][3])*URAD:(mat_min[0][3]-mat_stop[0][3])*URAD));
|
---|
| 291 | cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][3])?(mat_stop[0][3]-mat_min[0][3])*URAD:(mat_stop[0][3]-mat_max[0][3])*URAD))<<" µrad"<<endl;
|
---|
| 292 |
|
---|
| 293 | }
|
---|
| 294 |
|
---|
| 295 | delete temp_el;
|
---|
| 296 | return temp_vec;
|
---|
| 297 | }
|
---|