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

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


MagneticCollimator::MagneticCollimator(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, const double& Bmax, const double& thicknessMagneticField, const double& energyPerIon, const double& mass)
    : 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, Bmax, thicknessMagneticField, energyPerIon, mass)
{};

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


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

void MagneticCollimator::affiche()
{

    cout << "Magnetic collimator: " << endl;

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

double MagneticCollimator::Bcalc(const double& x, const double& y, const double& z, const double& w)
{

    double percent, n;

    int nn;

    percent = 0.1;
    double a(log(percent));
    double b(log((y - z) / y));
    n = a / b;
    nn = (int)floor(n);
    if ((nn % 2) == 1) {
        n = floor(n);
    } else {
        n = ceil(n);
    }

    return (w * pow((x / y), n)); //in Tesla
};

void MagneticCollimator::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)
{

    double sa, ca;

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


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

    double xl, xsl, pdepthStart, pdepthHalf, pdepthEnd, alf, lcolleff, z, B, BLength;


    //Calculate effective x and x' in the rotated collimator coordinate system
    //transforming to rotated system: [xl; yl]=[cos sin; -sin cos][x;y]
    //transforming back: [x;y]=[cos -sin; sin cos][xl;yl]

    xl = p1.coordonnees[0][0] * ca + p1.coordonnees[0][2] * sa;
    xsl = p1.coordonnees[0][1] * ca + p1.coordonnees[0][3] * sa;
    pdepthStart = abs(xl) - this->hgap; //the depth in the collimator at the beginning of the collimator
    pdepthHalf = abs(xl + this->L * xsl / 2) - (this->hgap2 + this->hgap) / 2; //depth in collimator at the middle of the collimator
    pdepthEnd = abs(xl + this->L * xsl) - this->hgap2; //depth in the collimator at the end of the collimator

    if (pdepthStart > 0) { //particle impacting on the front end of the collimator

        this->first.push_back(2);


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

        lcolleff = pdepthStart / alf; //alpha is negative for particles hitting on the face

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

        //Apply a thin collimator at the entrance

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

        p1.coordonnees[0][0] = p1.coordonnees[1][0];
        p1.coordonnees[0][1] = p1.coordonnees[1][1];
        p1.coordonnees[0][2] = p1.coordonnees[1][2];
        p1.coordonnees[0][3] = p1.coordonnees[1][3];
        p1.coordonnees[0][4] = p1.coordonnees[1][4];

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

        xsl = p1.coordonnees[0][1] * ca + p1.coordonnees[0][3] * sa;
        xl = p1.coordonnees[0][0] * ca + p1.coordonnees[0][2] * sa;
        pdepthEnd = abs(xl + this->L * xsl) - this->hgap2;
        p1.coordonnees[0][0] = p1.coordonnees[0][0] + scaleorbit * this->XC;
        p1.coordonnees[0][2] = p1.coordonnees[0][2] + scaleorbit * this->YC;

        if ((p1.Ap0 == 0) || (pdepthEnd >= 0)) { //particle lost, or stays inside the collimator's 2nd half
            p1.coordonnees[1][0] = p1.coordonnees[0][0] + p1.coordonnees[0][1] * delta_s;
            p1.coordonnees[1][2] = p1.coordonnees[0][2] + p1.coordonnees[0][3] * delta_s;
            p1.coordonnees[1][1] = p1.coordonnees[0][1];
            p1.coordonnees[1][3] = p1.coordonnees[0][3];
            this->second.push_back(-1);
            return;
        } else {//particle gets out of the collimator
            if (xl < 0) {
                alf = this->phi + xsl;
            } else {
                alf = this->phi - xsl;
            }
            z = pdepthStart / alf; //where particle comes out of the collimator
            z = (z + this->L) / 2; //midpoint between where particle gets out and the end of the collimator
            p1.coordonnees[1][0] = p1.coordonnees[0][0] + p1.coordonnees[0][1] * z;
            p1.coordonnees[1][2] = p1.coordonnees[0][2] + p1.coordonnees[0][3] * z;
            p1.coordonnees[1][1] = p1.coordonnees[0][1];
            p1.coordonnees[1][3] = p1.coordonnees[0][3];

            if (xl > 0) {
                B = this->Bmax;
            } else {
                B = -this->Bmax;
            }
            BLength = this->L - pdepthStart / alf;
        }

    }//end if(pdepthStart > 0)

    else if (pdepthHalf >= 0) { //particle enters the collimator in the first half

        double sImp, xint, yint, xsint, ysint;

        this->first.push_back(1);

        //Particles impacting along the jaw (outside the gap at the end but inside in beginning)

        if ((xl + this->L * xsl / 2) < 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 = pdepthEnd / alf; //effective length of orbit inside collimator

        if ((lcolleff > this->L) || (lcolleff < 0)) {
            lcolleff = this->L;
        }
        sImp = -pdepthStart / 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];

        double p1xtemp, p1ytemp, p1x0temp, p1xs0temp, p1y0temp, p1ys0temp;;

        p1xtemp = p1.coordonnees[0][0];
        p1ytemp = p1.coordonnees[0][2];
        p1x0temp = p1.coordonnees[0][0];
        p1xs0temp = p1.coordonnees[0][1];
        p1y0temp = p1.coordonnees[0][2];
        p1ys0temp = p1.coordonnees[0][3];

        p1.coordonnees[0][0] = xint;
        p1.coordonnees[0][1] = xsint;
        p1.coordonnees[0][2] = yint;
        p1.coordonnees[0][3] = ysint;

        //Apply thin collimator at impact position

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

        xint = p1.coordonnees[1][0];
        xsint = p1.coordonnees[1][1];
        yint = p1.coordonnees[1][2];
        ysint = p1.coordonnees[1][3];
        p1.coordonnees[0][4] = p1.coordonnees[1][4];

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

        xl = (xint + xsint * (this->L - sImp)) * ca + (yint + ysint * (this->L - sImp)) * sa;
        pdepthEnd = abs(xl) - this->hgap2;

        xint = xint + scaleorbit * this->XC; //coordinates not relative to reference particle/middle point between collimators
        yint = yint + scaleorbit * this->YC;
        p1.coordonnees[0][0] = p1x0temp + scaleorbit * this->XC;
        p1.coordonnees[0][2] = p1y0temp + scaleorbit * this->YC;
        p1.coordonnees[0][1] = p1xs0temp;
        p1.coordonnees[0][3] = p1ys0temp;

        if ((pdepthEnd >= 0) || (p1.Ap0 == 0)) {
            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;
            this->second.push_back(-1);
            return;
        } else {

            //Move particle as defined in MADX.
            //Drift => angles don't change, x,y change as straight line
            p1.coordonnees[1][0] = xint + xsint * (this->L - sImp) / 2;
            p1.coordonnees[1][2] = yint + ysint * (this->L - sImp) / 2;
            p1.coordonnees[1][1] = xsint;
            p1.coordonnees[1][3] = ysint;

            if (xl > 0) {
                B = this->Bmax;
            } else {
                B = -this->Bmax;
            }
            BLength = this->L;
            z = (this->L + sImp) / 2;
        }

    }//end if(pdepthHalf >=0)
    else {//particles that don't hit the first half of the collimator: drift to the middle of the collimator

        this->first.push_back(0);
        p1.coordonnees[0][0] = p1.coordonnees[0][0] + scaleorbit * this->XC;
        p1.coordonnees[0][2] = p1.coordonnees[0][2] + scaleorbit * this->YC;
        p1.coordonnees[1][0] = p1.coordonnees[0][0] + p1.coordonnees[0][1] * this->L / 2;
        p1.coordonnees[1][2] = p1.coordonnees[0][2] + p1.coordonnees[0][3] * this->L / 2;
        xl = (p1.coordonnees[1][0] - scaleorbit * this->XC) * ca + (p1.coordonnees[1][2] - scaleorbit * this->YC) * sa;
        p1.coordonnees[1][1] = p1.coordonnees[0][1];
        p1.coordonnees[1][3] = p1.coordonnees[0][3];
        p1.coordonnees[1][4] = p1.coordonnees[0][4];
        BLength = this->L;
        B = Bcalc(xl, (this->hgap + this->hgap2) / 2, this->thicknessMagneticField, this->Bmax);
        z = this->L / 2;

    }//end dernier else


    //applying kick

    double Preference, p, delta_angle, pdepthZ, sImp, xint, xsint, yint, ysint;

    Preference = sqrt(this->energyPerIon * this->energyPerIon - this->mass * this->mass) / Apr; //momentum per nucleon for reference particle. unit: GeV/c
    p = (1 + p1.coordonnees[0][4]) * Preference * p1.Ap0; //momentum. unit: GeV/c
    delta_angle = B * p1.Zp0 * BLength * 299792458 / (p * 1000000000);
    p1.coordonnees[1][1] = p1.coordonnees[1][1] + delta_angle * ca;
    p1.coordonnees[1][3] = p1.coordonnees[1][3] + delta_angle * sa;

    p1.coordonnees[1][0] = p1.coordonnees[1][0] - scaleorbit * this->XC; //coordinates relative to reference particle/middle point between collimators
    p1.coordonnees[1][2] = p1.coordonnees[1][2] - scaleorbit * this->YC;
    xl = p1.coordonnees[1][0] * ca + p1.coordonnees[1][2] * sa; //ca is cosine of collimator's angle
    xsl = p1.coordonnees[1][1] * ca + p1.coordonnees[1][3] * sa;

    pdepthZ = abs(xl) - (this->hgap * (this->L - z) + this->hgap2 * z) / this->L; //depth at the point the kick is given
    pdepthEnd = abs(xl + (this->L - z) * xsl) - this->hgap2;
    if (pdepthEnd > 0) { //enters collimator

        this->second.push_back(1);

        //Particles impacting along the jaw (outside the gap at the end but inside in beginning)

        if (xl < 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 = pdepthEnd / alf; //effective length of orbit inside collimator

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

        sImp = -pdepthZ / alf; //distance from current position to impact


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

        double p1xtemp, p1ytemp, p1x0temp, p1xs0temp, p1y0temp, p1ys0temp;

        p1xtemp = p1.coordonnees[1][0];
        p1ytemp = p1.coordonnees[1][2];
        p1x0temp = p1.coordonnees[0][0];
        p1xs0temp = p1.coordonnees[0][1];
        p1y0temp = p1.coordonnees[0][2];
        p1ys0temp = p1.coordonnees[0][3];

        p1.coordonnees[0][0] = xint;
        p1.coordonnees[0][1] = xsint;
        p1.coordonnees[0][2] = yint;
        p1.coordonnees[0][3] = ysint;

        //Apply thin collimator at impact position

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

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

        xint = p1.coordonnees[1][0];
        xsint = p1.coordonnees[1][1];
        yint = p1.coordonnees[1][2];
        ysint = p1.coordonnees[1][3];
        p1.coordonnees[0][4] = p1.coordonnees[1][4];

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

        xint = xint + scaleorbit * this->XC;
        yint = yint + scaleorbit * this->YC;
        p1.coordonnees[1][0] = p1.coordonnees[1][0] + scaleorbit * this->XC;
        p1.coordonnees[1][2] = p1.coordonnees[1][2] + scaleorbit * this->YC;
        p1.coordonnees[0][0] = p1x0temp;
        p1.coordonnees[0][2] = p1y0temp;
        p1.coordonnees[0][1] = p1xs0temp;
        p1.coordonnees[0][3] = p1ys0temp;

        //Move particle to next element of the collimator region as defined in MADX.
        //Drift => angles don't change, x,y change as straight line

        p1.coordonnees[1][0] = xint + xsint * (delta_s - z);
        p1.coordonnees[1][2] = yint + ysint * (delta_s - z);
        p1.coordonnees[1][1] = xsint;
        p1.coordonnees[1][3] = ysint;

    }//end if(pdepthend > 0)
    else {

        this->second.push_back(0);
        p1.coordonnees[1][0] = p1.coordonnees[1][0] + scaleorbit * this->XC;
        p1.coordonnees[1][2] = p1.coordonnees[1][2] + scaleorbit * this->YC;
        p1.coordonnees[1][0] = p1.coordonnees[1][0] + p1.coordonnees[1][1] * (delta_s - z);
        p1.coordonnees[1][2] = p1.coordonnees[1][2] + p1.coordonnees[1][3] * (delta_s - z);
        p1.coordonnees[1][4] = p1.coordonnees[0][4];

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

};