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