Collimator.cpp 17.7 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
// ------------------------------------------------------------------------
// $RCSfile: Collimator.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Collimator
//   Defines the abstract interface for a beam Collimator.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
// ------------------------------------------------------------------------

#include "AbsBeamline/Collimator.h"
kraus's avatar
kraus committed
22
#include "Physics/Physics.h"
23
#include "Algorithms/PartBunch.h"
gsell's avatar
gsell committed
24
#include "AbsBeamline/BeamlineVisitor.h"
25
#include "Fields/Fieldmap.h"
gsell's avatar
gsell committed
26
#include "Structure/LossDataSink.h"
kraus's avatar
kraus committed
27
#include "Utilities/Options.h"
kraus's avatar
kraus committed
28
#include "Solvers/SurfacePhysicsHandler.hh"
29
#include <memory>
gsell's avatar
gsell committed
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

extern Inform *gmsg;

using namespace std;

// Class Collimator
// ------------------------------------------------------------------------

Collimator::Collimator():
    Component(),
    filename_m(""),
    plane_m(OFF),
    position_m(0.0),
    PosX_m(0),
    PosY_m(0),
    PosZ_m(0),
    MomentumX_m(0),
    MomentumY_m(0),
    MomentumZ_m(0),
    time_m(0),
    id_m(0),
    informed_m(false),
    a_m(0.0),
    b_m(0.0),
    x0_m(0.0),
    y0_m(0.0),
56 57 58 59 60
    xstart_m(0.0),
    xend_m(0.0),
    ystart_m(0.0),
    yend_m(0.0),
    width_m(0.0),
gsell's avatar
gsell committed
61 62 63 64 65 66 67 68
    isAPepperPot_m(false),
    isASlit_m(false),
    isARColl_m(false),
    isACColl_m(false),
    isAWire_m(false),
    rHole_m(0.0),
    nHolesX_m(0),
    nHolesY_m(0),
69
    pitch_m(0.0),
70
    losses_m(0),
71 72
    lossDs_m(nullptr),
    sphys_m(NULL)
gsell's avatar
gsell committed
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
{}


Collimator::Collimator(const Collimator &right):
    Component(right),
    filename_m(right.filename_m),
    plane_m(right.plane_m),
    position_m(right.position_m),
    PosX_m(right.PosX_m),
    PosY_m(right.PosY_m),
    PosZ_m(right.PosZ_m),
    MomentumX_m(right.MomentumX_m),
    MomentumY_m(right.MomentumY_m),
    MomentumZ_m(right.MomentumZ_m),
    time_m(right.time_m),
    id_m(right.id_m),
    informed_m(right.informed_m),
    a_m(right.a_m),
    b_m(right.b_m),
    x0_m(right.x0_m),
    y0_m(right.y0_m),
94 95 96 97
    xstart_m(right.xstart_m),
    xend_m(right.xend_m),
    ystart_m(right.ystart_m),
    yend_m(right.yend_m),
98 99
    zstart_m(right.zstart_m),
    zend_m(right.zend_m),
100
    width_m(right.width_m),
gsell's avatar
gsell committed
101 102 103 104 105 106 107 108
    isAPepperPot_m(right.isAPepperPot_m),
    isASlit_m(right.isASlit_m),
    isARColl_m(right.isARColl_m),
    isACColl_m(right.isACColl_m),
    isAWire_m(right.isAWire_m),
    rHole_m(right.rHole_m),
    nHolesX_m(right.nHolesX_m),
    nHolesY_m(right.nHolesY_m),
109
    pitch_m(right.pitch_m),
110
    losses_m(0),
111 112
    lossDs_m(nullptr),
    sphys_m(NULL)
gsell's avatar
gsell committed
113
{
kraus's avatar
kraus committed
114
    setGeom();
gsell's avatar
gsell committed
115 116 117
}


118
Collimator::Collimator(const std::string &name):
gsell's avatar
gsell committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
    Component(name),
    filename_m(""),
    plane_m(OFF),
    position_m(0.0),
    PosX_m(0),
    PosY_m(0),
    PosZ_m(0),
    MomentumX_m(0),
    MomentumY_m(0),
    MomentumZ_m(0),
    time_m(0),
    id_m(0),
    informed_m(false),
    a_m(0.0),
    b_m(0.0),
    x0_m(0.0),
    y0_m(0.0),
136 137 138 139
    xstart_m(0.0),
    xend_m(0.0),
    ystart_m(0.0),
    yend_m(0.0),
140 141
    zstart_m(0.0),
    zend_m(0.0),
142
    width_m(0.0),
gsell's avatar
gsell committed
143 144 145 146 147 148 149 150
    isAPepperPot_m(false),
    isASlit_m(false),
    isARColl_m(false),
    isACColl_m(false),
    isAWire_m(false),
    rHole_m(0.0),
    nHolesX_m(0),
    nHolesY_m(0),
151
    pitch_m(0.0),
152
    losses_m(0),
153 154
    lossDs_m(nullptr),
    sphys_m(NULL)
gsell's avatar
gsell committed
155 156 157 158
{}


Collimator::~Collimator() {
kraus's avatar
kraus committed
159 160
    if(online_m)
        goOffline();
gsell's avatar
gsell committed
161 162 163 164 165 166 167
}


void Collimator::accept(BeamlineVisitor &visitor) const {
    visitor.visitCollimator(*this);
}

168 169 170

inline bool Collimator::isInColl(Vector_t R, Vector_t P, double recpgamma) {
    /**
kraus's avatar
kraus committed
171 172
       check if we are in the longitudinal
       range of the collimator
173
    */
kraus's avatar
kraus committed
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
    const double z = R(2) + P(2) * recpgamma;

    if((z > position_m) && (z <= position_m + getElementLength())) {
        if(isAPepperPot_m) {

            /**
               ------------
               |(0)|  |(0)|
               ----   -----
               |    a)    |
               |          |
               ----   -----
               |(0)|  |(0)|
               yL------------
               xL
               |---| d
               |--| pitch
               Observation: the area in a) is much larger than the
               area(s) (0). In a) particles are lost in (0)
               particles they are not lost.

            */
            const double h  =   pitch_m;
            const double xL = - 0.5 * h * (nHolesX_m - 1);
            const double yL = - 0.5 * h * (nHolesY_m - 1);
            bool alive = false;

            for(unsigned int m = 0; (m < nHolesX_m && (!alive)); m++) {
                for(unsigned int n = 0; (n < nHolesY_m && (!alive)); n++) {
                    double x_m = xL  + (m * h);
                    double y_m = yL  + (n * h);
                    /** are we in a) ? */
                    double rr = std::pow((R(0) - x_m) / rHole_m, 2) + std::pow((R(1) - y_m) / rHole_m, 2);
                    alive = (rr < 1.0);
                }
            }
            return !alive;
        } else if(isASlit_m) {
            return (R(0) <= -getXsize() || R(1) <= -getYsize() || R(0) >= getXsize() || R(1) >= getYsize());
        } else if (isARColl_m) {
            return (R(0) <= -getXsize() || R(1) <= -getYsize() || R(0) >= getXsize() || R(1) >= getYsize());
        } else if(isAWire_m) {
            ERRORMSG("Not yet implemented");
        } else {
            // case of an elliptic collimator
            const double trm1 = ((R(0)*R(0))/(getXsize()*getXsize()));
            const double trm2 = ((R(1)*R(1))/(getYsize()*getYsize()));
            return (trm1 + trm2) > 1.0;
        }
223
    }
kraus's avatar
kraus committed
224
    return false;
225 226
}

227
bool Collimator::apply(const size_t &i, const double &t, double E[], double B[]) {
gsell's avatar
gsell committed
228 229 230 231
    Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
    return apply(i, t, Ev, Bv);
}

232
bool Collimator::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
gsell's avatar
gsell committed
233 234 235 236
    const Vector_t &R = RefPartBunch_m->R[i] - Vector_t(dx_m, dy_m, ds_m); // including the missaligment
    const Vector_t &P = RefPartBunch_m->P[i];
    const double recpgamma = Physics::c * RefPartBunch_m->getdT() / sqrt(1.0  + dot(P, P));
    bool pdead = false;
237 238 239
    pdead = isInColl(R,P,recpgamma);

    if(pdead) {
kraus's avatar
kraus committed
240 241 242
        double frac = (R(2) - position_m) / P(2) * recpgamma;
        lossDs_m->addParticle(R,P, RefPartBunch_m->ID[i], t + frac * RefPartBunch_m->getdT(), 0);
        ++losses_m;
kraus's avatar
kraus committed
243
    }
gsell's avatar
gsell committed
244 245 246 247 248 249 250
    return pdead;
}

bool Collimator::apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) {
    return false;
}

251 252
bool Collimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {

kraus's avatar
kraus committed
253 254 255 256 257 258 259 260 261 262 263
    double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
    double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
    double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
    bool isDead = false;
    if(rmax(2) >= zstart_m && rmin(2) <= zend_m) {
        if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
            if(r(2) < zend_m && r(2) > zstart_m ) {
                int pflag = checkPoint(r(0), r(1));
                isDead = (pflag != 0);
            }
        }
264
    }
kraus's avatar
kraus committed
265
    return isDead;
266 267
}

gsell's avatar
gsell committed
268 269

// rectangle collimators in cyclotron cyclindral coordiantes
270
// without surfacephysics, the particle hitting collimator is deleted directly
gsell's avatar
gsell committed
271 272 273 274
bool Collimator::checkCollimator(PartBunch &bunch, const int turnnumber, const double t, const double tstep) {

    bool flagNeedUpdate = false;
    Vector_t rmin, rmax;
275

gsell's avatar
gsell committed
276
    bunch.get_bounds(rmin, rmax);
277 278
    double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
    double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
gsell's avatar
gsell committed
279
    double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
280 281

    if(rmax(2) >= zstart_m && rmin(2) <= zend_m) {
kraus's avatar
kraus committed
282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
        if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
            size_t tempnum = bunch.getLocalNum();
            int pflag = 0;
            for(unsigned int i = 0; i < tempnum; ++i) {
                if(bunch.PType[i] == ParticleType::REGULAR && bunch.R[i](2) < zend_m && bunch.R[i](2) > zstart_m ) {
                    pflag = checkPoint(bunch.R[i](0), bunch.R[i](1));
                    if(pflag != 0) {
                        if (!sphys_m)
                            lossDs_m->addParticle(bunch.R[i], bunch.P[i], bunch.ID[i]);
                        bunch.Bin[i] = -1;
                        flagNeedUpdate = true;
                    }
                }
            }
        }
297 298 299
    }
    reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
    if (flagNeedUpdate && sphys_m) {
kraus's avatar
kraus committed
300
        sphys_m->apply(bunch);
gsell's avatar
gsell committed
301
    }
302
    return flagNeedUpdate;
gsell's avatar
gsell committed
303 304 305 306 307 308
}

void Collimator::initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) {
    RefPartBunch_m = bunch;
    position_m = startField;
    endField = position_m + getElementLength();
309

kraus's avatar
kraus committed
310 311
    sphys_m = getSurfacePhysics();

adelmann's avatar
adelmann committed
312
    if (!sphys_m) {
kraus's avatar
kraus committed
313 314 315 316
        if (filename_m == std::string(""))
            lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
        else
            lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
adelmann's avatar
adelmann committed
317
    }
kraus's avatar
kraus committed
318

kraus's avatar
kraus committed
319
    goOnline(-1e6);
gsell's avatar
gsell committed
320 321 322 323
}

void Collimator::initialise(PartBunch *bunch, const double &scaleFactor) {
    RefPartBunch_m = bunch;
324 325

    sphys_m = getSurfacePhysics();
adelmann's avatar
adelmann committed
326 327

    if (!sphys_m) {
kraus's avatar
kraus committed
328 329 330 331
        if (filename_m == std::string(""))
            lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
        else
            lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
adelmann's avatar
adelmann committed
332
    }
kraus's avatar
kraus committed
333

kraus's avatar
kraus committed
334
    goOnline(-1e6);
gsell's avatar
gsell committed
335 336 337 338
}

void Collimator::finalise()
{
kraus's avatar
kraus committed
339 340 341
    if(online_m)
        goOffline();
    *gmsg << "* Finalize probe" << endl;
gsell's avatar
gsell committed
342 343
}

kraus's avatar
kraus committed
344
void Collimator::goOnline(const double &) {
gsell's avatar
gsell committed
345 346
    if(RefPartBunch_m == NULL) {
        if(!informed_m) {
347
            std::string errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
348
            ERRORMSG(errormsg << endl);
gsell's avatar
gsell committed
349 350 351 352 353 354 355 356 357 358
            if(Ippl::myNode() == 0) {
                ofstream omsg("errormsg.txt", ios_base::app);
                omsg << errormsg << endl;
                omsg.close();
            }
            informed_m = true;
        }
        return;
    }

kraus's avatar
kraus committed
359 360 361 362 363 364 365 366 367 368 369 370 371 372
    *gmsg << level3;
    if(isAPepperPot_m)
        *gmsg << "* Pepperpot x= " << a_m << " y= " << b_m << " r= " << rHole_m
              << " nx= " << nHolesX_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
    else if(isASlit_m)
        *gmsg << "* Slit x= " << getXsize() << " Slit y= " << getYsize()
              << " start= " << position_m << " fn= " << filename_m << endl;
    else if(isARColl_m)
        *gmsg << "* RCollimator a= " << getXsize() << " b= " << b_m << " start= " << position_m
              << " fn= " << filename_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
    else if(isACColl_m)
        *gmsg << "* CCollimator angle start " << xstart_m << " (Deg) angle end " << ystart_m << " (Deg) "
              << "R start " << xend_m << " (mm) R rend " << yend_m << " (mm)" << endl;
    else if(isAWire_m)
373
        *gmsg << "* Wire x= " << x0_m << " y= " << y0_m << endl;
kraus's avatar
kraus committed
374 375 376 377
    else
        *gmsg << "* ECollimator a= " << getXsize() << " b= " << b_m << " start= " << position_m
              << " fn= " << filename_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;

gsell's avatar
gsell committed
378 379 380 381 382 383 384 385 386 387 388 389

    PosX_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    PosY_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    PosZ_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    MomentumX_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    MomentumY_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    MomentumZ_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    time_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    id_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    online_m = true;
}

adelmann's avatar
adelmann committed
390 391 392
void Collimator::print() {
    if(RefPartBunch_m == NULL) {
        if(!informed_m) {
393
            std::string errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
adelmann's avatar
adelmann committed
394
            *gmsg << errormsg << "\n"
kraus's avatar
kraus committed
395
                  << endl;
adelmann's avatar
adelmann committed
396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420
            if(Ippl::myNode() == 0) {
                ofstream omsg("errormsg.txt", ios_base::app);
                omsg << errormsg << endl;
                omsg.close();
            }
            informed_m = true;
        }
        return;
    }

    if(isAPepperPot_m)
        *gmsg << "Pepperpot x= " << a_m << " y= " << b_m << " r= " << rHole_m << " nx= " << nHolesX_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
    else if(isASlit_m)
        *gmsg << "Slit x= " << getXsize() << " Slit y= " << getYsize() << " start= " << position_m << " fn= " << filename_m << endl;
    else if(isARColl_m)
        *gmsg << "RCollimator a= " << getXsize() << " b= " << b_m << " start= " << position_m << " fn= " << filename_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
    else if(isACColl_m)
        *gmsg << "CCollimator angstart= " << xstart_m << " angend " << ystart_m << " rstart " << xend_m << " rend " << yend_m << endl;
    else if(isAWire_m)
        *gmsg << "Wire x= " << x0_m << " y= " << y0_m << endl;
    else
        *gmsg << "ECollimator a= " << getXsize() << " b= " << b_m << " start= " << position_m << " fn= " << filename_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;

}

gsell's avatar
gsell committed
421
void Collimator::goOffline() {
kraus's avatar
kraus committed
422 423 424
    if (online_m && lossDs_m)
        lossDs_m->save();
    online_m = false;
gsell's avatar
gsell committed
425 426 427 428 429 430
}

bool Collimator::bends() const {
    return false;
}

431
void Collimator::setOutputFN(std::string fn) {
gsell's avatar
gsell committed
432 433 434 435
    filename_m = fn;
}

string Collimator::getOutputFN() {
adelmann's avatar
adelmann committed
436 437
    if (filename_m == std::string(""))
        return getName();
kraus's avatar
kraus committed
438
    else
adelmann's avatar
adelmann committed
439
        return filename_m.substr(0, filename_m.rfind("."));
gsell's avatar
gsell committed
440 441
}

442 443 444 445
unsigned int Collimator::getLosses() const {
    return losses_m;
}

gsell's avatar
gsell committed
446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479
void Collimator::setXsize(double a) {
    a_m = a;
}

void Collimator::setYsize(double b) {
    b_m = b;
}

void Collimator::setXpos(double x0) {
    x0_m = x0;
}

void Collimator::setYpos(double y0) {
    y0_m = y0;
}


double Collimator::getXsize(double a) {
    return a_m;
}

double Collimator::getYsize(double b) {
    return b_m;
}

double Collimator::getXpos() {
    return x0_m;
}

double Collimator::getYpos() {
    return y0_m;

    // --------Cyclotron collimator
}
480 481
void Collimator::setXStart(double xstart) {
    xstart_m = xstart;
gsell's avatar
gsell committed
482 483
}

484 485
void Collimator::setXEnd(double xend) {
    xend_m = xend;
gsell's avatar
gsell committed
486 487
}

488 489
void Collimator::setYStart(double ystart) {
    ystart_m = ystart;
gsell's avatar
gsell committed
490 491
}

492 493
void Collimator::setYEnd(double yend) {
    yend_m = yend;
gsell's avatar
gsell committed
494 495
}

496 497 498 499 500 501 502 503
void Collimator::setZStart(double zstart) {
    zstart_m = zstart;
}

void Collimator::setZEnd(double zend) {
    zend_m = zend;
}

504 505
void Collimator::setWidth(double width) {
    width_m = width;
gsell's avatar
gsell committed
506 507
}

508 509
double Collimator::getXStart() {
    return xstart_m;
gsell's avatar
gsell committed
510 511
}

512 513
double Collimator::getXEnd() {
    return xend_m;
gsell's avatar
gsell committed
514 515
}

516 517
double Collimator::getYStart() {
    return ystart_m;
gsell's avatar
gsell committed
518 519
}

520 521
double Collimator::getYEnd() {
    return yend_m;
gsell's avatar
gsell committed
522 523
}

524 525 526 527 528 529 530 531
double Collimator::getZStart() {
    return zstart_m;
}

double Collimator::getZEnd() {
    return zend_m;
}

gsell's avatar
gsell committed
532
double Collimator::getWidth() {
533
    return width_m;
gsell's avatar
gsell committed
534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576
}

//-------------------------------

void Collimator::setRHole(double r) {
    rHole_m = r;
}
void Collimator::setNHoles(unsigned int nx, unsigned int ny) {
    nHolesX_m = nx;
    nHolesY_m = ny;
}
void Collimator::setPitch(double p) {
    pitch_m = p;
}


void Collimator::setPepperPot() {
    isAPepperPot_m = true;
}
void Collimator::setSlit() {
    isASlit_m = true;
}

void Collimator::setRColl() {
    isARColl_m = true;
}

void Collimator::setCColl() {
    isACColl_m = true;
}

void Collimator::setWire() {
    isAWire_m = true;
}
void Collimator::getDimensions(double &zBegin, double &zEnd) const {
    zBegin = position_m;
    zEnd = position_m + getElementLength();

    // zBegin = position_m - 0.005;
    //  zEnd = position_m + 0.005;

}

577 578
ElementBase::ElementType Collimator::getType() const {
    return COLLIMATOR;
gsell's avatar
gsell committed
579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597
}

string Collimator::getCollimatorShape() {
    if(isAPepperPot_m)
        return "PeperPot";
    else if(isASlit_m)
        return "Slit";
    else if(isARColl_m)
        return "RCollimator";
    else if(isACColl_m)
        return "CCollimator";
    else if(isAWire_m)
        return "Wire";
    else
        return "ECollimator";
}

void Collimator::setGeom() {

598 599
    double slope;
    if (xend_m == xstart_m)
kraus's avatar
kraus committed
600
        slope = 1.0e12;
601
    else
kraus's avatar
kraus committed
602
        slope = (yend_m - ystart_m) / (xend_m - xstart_m);
603 604 605 606 607 608

    double coeff2 = sqrt(1 + slope * slope);
    double coeff1 = slope / coeff2;
    double halfdist = width_m / 2.0;
    geom_m[0].x = xstart_m - halfdist * coeff1;
    geom_m[0].y = ystart_m + halfdist / coeff2;
gsell's avatar
gsell committed
609

610 611
    geom_m[1].x = xstart_m + halfdist * coeff1;
    geom_m[1].y = ystart_m - halfdist / coeff2;
gsell's avatar
gsell committed
612

613 614
    geom_m[2].x = xend_m + halfdist * coeff1;
    geom_m[2].y = yend_m - halfdist  / coeff2;
gsell's avatar
gsell committed
615

616 617
    geom_m[3].x = xend_m - halfdist * coeff1;
    geom_m[3].y = yend_m + halfdist / coeff2;
gsell's avatar
gsell committed
618 619 620 621

    geom_m[4].x = geom_m[0].x;
    geom_m[4].y = geom_m[0].y;

622
    if (zstart_m > zend_m){
kraus's avatar
kraus committed
623 624 625 626
        double tempz = 0.0;
        tempz = zstart_m;
        zstart_m = zend_m;
        zend_m = tempz;
627
    }
gsell's avatar
gsell committed
628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643
}


int Collimator::checkPoint(const double &x, const double &y) {
    int    cn = 0;

    for(int i = 0; i < 4; i++) {
        if(((geom_m[i].y <= y) && (geom_m[i+1].y > y))
           || ((geom_m[i].y > y) && (geom_m[i+1].y <= y))) {

            float vt = (float)(y - geom_m[i].y) / (geom_m[i+1].y - geom_m[i].y);
            if(x < geom_m[i].x + vt * (geom_m[i+1].x - geom_m[i].x))
                ++cn;
        }
    }
    return (cn & 1);  // 0 if even (out), and 1 if odd (in)
kraus's avatar
kraus committed
644
}