Collimator.cpp 17 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 114 115 116 117
{
  setGeom();
}


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() {
adelmann's avatar
adelmann 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 171 172 173 174

inline bool Collimator::isInColl(Vector_t R, Vector_t P, double recpgamma) {
  /**
     check if we are in the longitudinal
     range of the collimator
  */
  const double z = R(2) + P(2) * recpgamma;
kraus's avatar
kraus committed
175

176 177
  if((z > position_m) && (z <= position_m + getElementLength())) {
    if(isAPepperPot_m) {
kraus's avatar
kraus committed
178

179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
    /**
       ------------
       |(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.
kraus's avatar
kraus committed
194

195 196 197 198 199
    */
      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;
kraus's avatar
kraus committed
200

201 202 203 204 205 206 207 208 209 210 211
      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) {
212 213 214 215 216
        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");
kraus's avatar
kraus committed
217
    } else {
218 219
        // case of an elliptic collimator
        const double trm1 = ((R(0)*R(0))/(getXsize()*getXsize()));
kraus's avatar
kraus committed
220
        const double trm2 = ((R(1)*R(1))/(getYsize()*getYsize()));
221
        return (trm1 + trm2) > 1.0;
222 223 224 225 226
    }
  }
  return false;
}

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 240
    pdead = isInColl(R,P,recpgamma);

    if(pdead) {
      double frac = (R(2) - position_m) / P(2) * recpgamma;
adelmann's avatar
adelmann committed
241
      lossDs_m->addParticle(R,P, RefPartBunch_m->ID[i], t + frac * RefPartBunch_m->getdT(), 0);
242
      ++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 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
bool Collimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {

  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);
      }
    }
  }
  return isDead;
}

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) {
282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
      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] == 0 && 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;
	    }
	  }
	}
      }
    }
    reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
    if (flagNeedUpdate && sphys_m) {
      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 313 314 315 316 317
    if (!sphys_m) {
      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));
    }
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 328 329 330 331 332

    if (!sphys_m) {
      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));
    }
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()
{
adelmann's avatar
adelmann committed
339 340
  if(online_m)
    goOffline();
adelmann's avatar
Cleanup  
adelmann committed
341
  *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;
    }

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

    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
385 386 387
void Collimator::print() {
    if(RefPartBunch_m == NULL) {
        if(!informed_m) {
388
            std::string errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
adelmann's avatar
adelmann committed
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415
            *gmsg << errormsg << "\n"
                << endl;
            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
416
void Collimator::goOffline() {
adelmann's avatar
adelmann committed
417
  if (online_m && lossDs_m)
adelmann's avatar
adelmann committed
418
    lossDs_m->save();
adelmann's avatar
adelmann committed
419
  online_m = false;
gsell's avatar
gsell committed
420 421 422 423 424 425
}

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

426
void Collimator::setOutputFN(std::string fn) {
gsell's avatar
gsell committed
427 428 429 430
    filename_m = fn;
}

string Collimator::getOutputFN() {
adelmann's avatar
adelmann committed
431 432
    if (filename_m == std::string(""))
        return getName();
kraus's avatar
kraus committed
433
    else
adelmann's avatar
adelmann committed
434
        return filename_m.substr(0, filename_m.rfind("."));
gsell's avatar
gsell committed
435 436
}

437 438 439 440
unsigned int Collimator::getLosses() const {
    return losses_m;
}

gsell's avatar
gsell committed
441 442 443 444 445 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
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
}
475 476
void Collimator::setXStart(double xstart) {
    xstart_m = xstart;
gsell's avatar
gsell committed
477 478
}

479 480
void Collimator::setXEnd(double xend) {
    xend_m = xend;
gsell's avatar
gsell committed
481 482
}

483 484
void Collimator::setYStart(double ystart) {
    ystart_m = ystart;
gsell's avatar
gsell committed
485 486
}

487 488
void Collimator::setYEnd(double yend) {
    yend_m = yend;
gsell's avatar
gsell committed
489 490
}

491 492 493 494 495 496 497 498
void Collimator::setZStart(double zstart) {
    zstart_m = zstart;
}

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

499 500
void Collimator::setWidth(double width) {
    width_m = width;
gsell's avatar
gsell committed
501 502
}

503 504
double Collimator::getXStart() {
    return xstart_m;
gsell's avatar
gsell committed
505 506
}

507 508
double Collimator::getXEnd() {
    return xend_m;
gsell's avatar
gsell committed
509 510
}

511 512
double Collimator::getYStart() {
    return ystart_m;
gsell's avatar
gsell committed
513 514
}

515 516
double Collimator::getYEnd() {
    return yend_m;
gsell's avatar
gsell committed
517 518
}

519 520 521 522 523 524 525 526
double Collimator::getZStart() {
    return zstart_m;
}

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

gsell's avatar
gsell committed
527
double Collimator::getWidth() {
528
    return width_m;
gsell's avatar
gsell committed
529 530 531 532 533 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
}

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

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;

}

572 573
ElementBase::ElementType Collimator::getType() const {
    return COLLIMATOR;
gsell's avatar
gsell committed
574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592
}

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() {

593 594 595 596 597 598 599 600 601 602 603
    double slope;
    if (xend_m == xstart_m)
      slope = 1.0e12;
    else
      slope = (yend_m - ystart_m) / (xend_m - xstart_m);

    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
604

605 606
    geom_m[1].x = xstart_m + halfdist * coeff1;
    geom_m[1].y = ystart_m - halfdist / coeff2;
gsell's avatar
gsell committed
607

608 609
    geom_m[2].x = xend_m + halfdist * coeff1;
    geom_m[2].y = yend_m - halfdist  / coeff2;
gsell's avatar
gsell committed
610

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

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

617 618 619 620 621 622
    if (zstart_m > zend_m){
      double tempz = 0.0;
      tempz = zstart_m;
      zstart_m = zend_m;
      zend_m = tempz;
    }
gsell's avatar
gsell committed
623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639
}


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)
}