source: ZHANGProjects/ICOSIM/CPP/trunk/source/StandardCollimator.cc @ 2

#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include "StandardCollimator.h"
using namespace std;


StandardCollimator::StandardCollimator(const double& ALFX, const double& ALFY, const double& APER_1, const double& APER_2, const double& APER_3, const double& APER_4, const string& APERTYPE, const double& BETX, const double& BETY, const double& DPX, const double& DPY, const double& DX, const double& DY, const string& KEYWORD, const double& L, const double& MUX, const double& MUY, const string& NAME, const double& PTC, const double& PXC, const double& PYC, const double& S, const double& TC, const double& XC, const double& YC, const double& K0L, const double& K0SL, const double& K1L, const double& K1SL, const double& K2L, const double& K2SL, const string& PARENT, const string& meth, const long double& hgap, const long double& hgap2, const double& collang, const long double& pdepth, const long double& pdepth2, const double& tcang, const double& nsig)
    : Collimator(ALFX, ALFY, APER_1, APER_2, APER_3, APER_4, APERTYPE, BETX, BETY, DPX, DPY, DX, DY, KEYWORD, L, MUX, MUY, NAME, PTC, PXC, PYC, S, TC, XC, YC, K0L, K0SL, K1L, K1SL, K2L, K2SL, PARENT, meth, hgap, hgap2, collang, pdepth, pdepth2, tcang, nsig)
{};

StandardCollimator::StandardCollimator(Element elt, const double& tcang, const double& nsig, const string& meth, const string& material)
    : Collimator(elt, tcang, nsig, meth, material)
{
};


StandardCollimator::StandardCollimator(const StandardCollimator& obj)
    : Collimator(obj)
{};


void StandardCollimator::collipass(Particle& p1, double& dpopeff, const double& scaleorbit, const double& R11X, const double& R12X, const double& R21X, const double& R22X, const double& R11Y, const double& R12Y, const double& R21Y, const double& R22Y, const double& dx1, const double& dpx1, const double& dy1, const double& dpy1, const double& delta_s, const double& Apr, const double& Zpr, const double& betgam)
{

    long double ca, sa;

    sa = sin(this->tcang);
    ca = cos(this->tcang);

    //account for positioning of collmators on orbit
    p1.coordonnees[0][0] = p1.coordonnees[0][0] - scaleorbit * this->XC;
    p1.coordonnees[0][2] = p1.coordonnees[0][2] - scaleorbit * this->YC;

    long double xl, xsl, alf, lcolleff;


    //Calculate effective x and x' in the rotated collimator coordinate system
    xl = p1.coordonnees[0][0] * ca + p1.coordonnees[0][2] * sa;
    xsl = p1.coordonnees[0][1] * ca + p1.coordonnees[0][3] * sa;
    this->pdepth = abs(xl) - this->hgap;//the depth in the collimator at the entrance
    this->pdepth2 = abs(xl + this->L * xsl) - this->hgap2; //the depth in the collimator at the exit


    if (this->pdepth > 0) { //the particle is impacting on the front end of the collimator


        //impact angle on collimator
        if (xl < 0) {
            alf = this->phi + xsl;//impact angle on collimator
        } else {
            alf = this->phi - xsl;
        }

        //alpha is negative for particles hitting on the face

        lcolleff = this->pdepth / alf; //effective length of orbit inside collimator

        if ((lcolleff > this->L) || (lcolleff < 0)) {
            lcolleff = this->L;
        }

        //Apply a thin collimator at the entrance

        collipassInteraction(p1, Apr, Zpr, betgam, lcolleff);

        if (p1.Ap0 == 0) {
            return;
        }

        p1.coordonnees[0][0] = p1.coordonnees[1][0] + scaleorbit * this->XC;
        p1.coordonnees[0][2] = p1.coordonnees[1][2] + scaleorbit * this->YC;
        p1.coordonnees[0][4] = p1.coordonnees[1][4];

        //Twiss transform the particle to the end of the collimator

        dpopeff = (p1.Ap0 * Zpr) / (p1.Zp0 * Apr) * (1 + p1.coordonnees[0][4]) - 1;

        p1.coordonnees[1][0] = R11X * p1.coordonnees[0][0] + R12X * p1.coordonnees[1][1] + (this->DX - R11X * dx1 - R12X * dpx1) * dpopeff;
        p1.coordonnees[1][1] = R21X * p1.coordonnees[0][0] + R22X * p1.coordonnees[1][1] + (this->DPX - R21X * dx1 - R22X * dpx1) * dpopeff;
        p1.coordonnees[1][2] = R11Y * p1.coordonnees[0][2] + R12Y * p1.coordonnees[1][3] + (this->DY - R11Y * dy1 - R12Y * dpy1) * dpopeff;
        p1.coordonnees[1][3] = R21Y * p1.coordonnees[0][2] + R22Y * p1.coordonnees[1][3] + (this->DPY - R21Y * dy1 - R22Y * dpy1) * dpopeff;

    } else if (this->pdepth2 > 0) { //Particles impacting along the jaw (outside the gap at the end but inside in beginning)

        double sImp, xint, yint, xsint, ysint;

        if ((xl + L * xsl) < 0) {
            alf = -this->phi - xsl;
        } else {
            alf = -this->phi + xsl;
        }

        //impact angle on collimator. The sign determines which jaw is hit and therefore if the jaw angle should be added or subtracted.
        lcolleff = this->pdepth2 / alf;

        if ((lcolleff > this->L) || (lcolleff < 0)) {
            lcolleff = this->L;
        }

        sImp = -this->pdepth / alf; //distance from beginning of collimator to impact

        //Move particle to this point. Drift => angles don't change, x,y change as straight line
        xint = p1.coordonnees[0][0] + p1.coordonnees[0][1] * sImp;
        yint = p1.coordonnees[0][2] + p1.coordonnees[0][3] * sImp;
        xsint = p1.coordonnees[0][1];
        ysint = p1.coordonnees[0][3];

        //Apply thin collimator at impact position

        Particle ptemp(xint, xsint, yint, ysint, p1.coordonnees[0][4], 0, p1.Ap0, p1.Zp0, p1.coordonnees[0][4]);

        collipassInteraction(ptemp, Apr, Zpr, betgam, lcolleff);

        if (ptemp.Ap0 == 0) {
            return;
        }


        //Move particle to the end of collimator region as defined in MADX. Drift => angles don't change, x,y change as straight line
        xint = ptemp.coordonnees[1][0];
        yint = ptemp.coordonnees[1][2];
        xsint = ptemp.coordonnees[1][1];
        ysint = ptemp.coordonnees[1][3];
        p1.Ap0 = ptemp.Ap0;
        p1.Zp0 = ptemp.Zp0;
        p1.coordonnees[0][4] = ptemp.coordonnees[1][4];

        xint = xint + scaleorbit * this->XC;
        yint = yint + scaleorbit * this->YC;
        p1.coordonnees[0][0] = p1.coordonnees[0][0] + scaleorbit * this->XC;
        p1.coordonnees[0][2] = p1.coordonnees[0][2] + scaleorbit * this->YC;

        dpopeff = (p1.Ap0 * Zpr) / (p1.Zp0 * Apr) * (1 + p1.coordonnees[0][4]) - 1;

        p1.coordonnees[1][0] = xint + xsint * (delta_s - sImp);
        p1.coordonnees[1][2] = yint + ysint * (delta_s - sImp);
        p1.coordonnees[1][1] = xsint;
        p1.coordonnees[1][3] = ysint;
        p1.coordonnees[1][4] = p1.coordonnees[0][4];
    } else {//Particle don't hit the collimator

        p1.coordonnees[0][0] = p1.coordonnees[0][0] + scaleorbit * this->XC;
        p1.coordonnees[0][2] = p1.coordonnees[0][2] + scaleorbit * this->YC;

        //Twiss transform the particle to the end of the collimator
        dpopeff = (p1.Ap0 * Zpr) / (p1.Zp0 * Apr) * (1 + p1.coordonnees[0][4]) - 1;

        p1.coordonnees[1][0] = R11X * p1.coordonnees[0][0] + R12X * p1.coordonnees[0][1] + (this->DX - R11X * dx1 - R12X * dpx1) * dpopeff;
        p1.coordonnees[1][1] = R21X * p1.coordonnees[0][0] + R22X * p1.coordonnees[0][1] + (this->DPX - R21X * dx1 - R22X * dpx1) * dpopeff;
        p1.coordonnees[1][2] = R11Y * p1.coordonnees[0][2] + R12Y * p1.coordonnees[0][3] + (this->DY - R11Y * dy1 - R12Y * dpy1) * dpopeff;
        p1.coordonnees[1][3] = R21Y * p1.coordonnees[0][2] + R22Y * p1.coordonnees[0][3] + (this->DPY - R21Y * dy1 - R22Y * dpy1) * dpopeff;
        p1.coordonnees[1][4] = p1.coordonnees[0][4];

    }
};

void StandardCollimator::affiche()
{

    cout << "Standard collimator:" << endl;

    this->Collimator::affiche();
};