Stripper.cpp 10.1 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 22
// ------------------------------------------------------------------------
// $RCSfile: Stripper.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Stripper
//   Defines the abstract interface for a stripper
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2011/07/08 11:16:04 $
// $Author: Jianjun Yang $
//
// ------------------------------------------------------------------------

#include "AbsBeamline/Stripper.h"
#include "AbsBeamline/BeamlineVisitor.h"
23
#include "Algorithms/PartBunch.h"
gsell's avatar
gsell committed
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
#include "Physics/Physics.h"
#include "Structure/LossDataSink.h"
#include <iostream>
#include <fstream>
#include "AbstractObjects/OpalData.h"


using Physics::pi;
using Physics::q_e;

extern Inform *gmsg;

using namespace std;

// Class Stripper
// ------------------------------------------------------------------------

Stripper::Stripper():
    Component(),
    filename_m(""),
    position_m(0.0),
    xstart_m(0.0),
    xend_m(0.0),
    ystart_m(0.0),
    yend_m(0.0),
    width_m(0.0),
    opcharge_m(0.0),
    opmass_m(0.0),
    stop_m(true),
    step_m(0) {
54
    A_m = yend_m - ystart_m;
55
    B_m = xstart_m - xend_m;
56 57
    R_m = sqrt(A_m*A_m+B_m*B_m);
    C_m = ystart_m*xend_m - xstart_m*yend_m;
gsell's avatar
gsell committed
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
}


Stripper::Stripper(const Stripper &right):
    Component(right),
    filename_m(right.filename_m),
    position_m(right.position_m),
    xstart_m(right.xstart_m),
    xend_m(right.xend_m),
    ystart_m(right.ystart_m),
    yend_m(right.yend_m),
    width_m(right.width_m),
    opcharge_m(right.opcharge_m),
    opmass_m(right.opmass_m),
    stop_m(right.stop_m),
    step_m(right.step_m) {
74
    A_m = yend_m - ystart_m;
75
    B_m = xstart_m - xend_m;
76 77
    R_m = sqrt(A_m*A_m+B_m*B_m);
    C_m = ystart_m*xend_m - xstart_m*yend_m;
gsell's avatar
gsell committed
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
}


Stripper::Stripper(const string &name):
    Component(name),
    filename_m(""),
    position_m(0.0),
    xstart_m(0.0),
    xend_m(0.0),
    ystart_m(0.0),
    yend_m(0.0),
    width_m(0.0),
    opcharge_m(0.0),
    opmass_m(0.0),
    stop_m(true),
    step_m(0){
94
    A_m = yend_m - ystart_m;
95
    B_m = xstart_m - xend_m;
96 97
    R_m = sqrt(A_m*A_m+B_m*B_m);
    C_m = ystart_m*xend_m - xstart_m*yend_m;
gsell's avatar
gsell committed
98 99 100 101
}

void Stripper::setGeom(const double dist) {

102
   double slope;
103
    if (xend_m == xstart_m)
104 105 106
      slope = 1.0e12;
    else
      slope = (yend_m - ystart_m) / (xend_m - xstart_m);
107

108 109 110 111 112
    double coeff2 = sqrt(1 + slope * slope);
    double coeff1 = slope / coeff2;
    double halfdist = dist / 2.0;
    geom_m[0].x = xstart_m - halfdist * coeff1;
    geom_m[0].y = ystart_m + halfdist / coeff2;
gsell's avatar
gsell committed
113

114 115
    geom_m[1].x = xstart_m + halfdist * coeff1;
    geom_m[1].y = ystart_m - halfdist / coeff2;
gsell's avatar
gsell committed
116

117 118
    geom_m[2].x = xend_m + halfdist * coeff1;
    geom_m[2].y = yend_m - halfdist  / coeff2;
gsell's avatar
gsell committed
119

120 121
    geom_m[3].x = xend_m - halfdist * coeff1;
    geom_m[3].y = yend_m + halfdist / coeff2;
gsell's avatar
gsell committed
122 123 124

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

gsell's avatar
gsell committed
126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
}


Stripper::~Stripper() {
    idrec_m.clear();
}


void Stripper::accept(BeamlineVisitor &visitor) const {
    visitor.visitStripper(*this);
}

bool Stripper::apply(const int &i, const double &t, double E[], double B[]) {
    *gmsg << "stripper1" << endl;
    Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
    return false;
}

bool Stripper::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
    *gmsg << "stripper2" << endl;
    return false;
}

bool Stripper::apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) {
    *gmsg << "stripper3" << endl;
    return false;
}

void Stripper::initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) {

}

void Stripper::initialise(PartBunch *bunch, const double &scaleFactor) {
    // initialize DataSink with H5Part output enabled
    bool doH5 = false;
    lossDs_m = new LossDataSink(bunch->getTotalNum(), doH5);
162
    //lossDs_m->openH5(getName());
gsell's avatar
gsell committed
163 164 165 166 167 168 169 170 171 172 173 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 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
}

void Stripper::finalise() {
    *gmsg << "Finalize probe" << endl;
    if(lossDs_m)
        delete lossDs_m;
}

bool Stripper::bends() const {
    return false;
}
void Stripper::goOnline() {
    online_m = true;
}

void Stripper::goOffline() {
    online_m = false;
}

void  Stripper::setXstart(double xstart) {
    xstart_m = xstart;
}

void  Stripper::setXend(double xend) {
    xend_m = xend;
}

void  Stripper::setYstart(double ystart) {
    ystart_m = ystart;
}

void  Stripper::setYend(double yend) {
    yend_m = yend;
}
void  Stripper::setWidth(double width) {
    width_m = width;
}

void  Stripper::setOPCharge(double charge) {
    opcharge_m = charge;
}

void  Stripper::setOPMass(double mass) {
    opmass_m = mass;
}

void  Stripper::setStop(bool flag) {
    stop_m = flag;

}

double  Stripper::getXstart() const {
    return xstart_m;
}

double  Stripper::getXend() const {
    return xend_m;
}

double  Stripper::getYstart() const {
    return ystart_m;
}

double  Stripper::getYend() const {
    return yend_m;
}
double  Stripper::getWidth() const {
    return width_m;
}

double  Stripper::getOPCharge() const {
    return opcharge_m;
}

double  Stripper::getOPMass() const {
    return opmass_m;
}

bool  Stripper::getStop () const {
    return stop_m;
}


//change the stripped particles to outcome particles
bool  Stripper::checkStripper(PartBunch &bunch, const int turnnumber, const double t, const double tstep) {

    bool flagNeedUpdate = false;
    bool flagresetMQ = false;
    Vector_t rmin, rmax, strippoint;
    bunch.get_bounds(rmin, rmax);
    double r_ref = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
    double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));

    if(r1 > r_ref - 10.0 ){

        size_t count = 0;
        size_t tempnum = bunch.getLocalNum();
        int pflag = 0;
261

gsell's avatar
gsell committed
262 263 264 265 266 267 268 269 270 271
	Vector_t meanP(0.0, 0.0, 0.0);
	for(unsigned int i = 0; i < bunch.getLocalNum(); ++i) {
	  for(int d = 0; d < 3; ++d) {
            meanP(d) += bunch.P[i](d);
	  }
	}
	reduce(meanP, meanP, OpAddAssign());
	meanP = meanP / Vector_t(bunch.getTotalNum());

	double sk1, sk2, stangle = 0.0;
272
	if ( B_m == 0.0 ){
gsell's avatar
gsell committed
273
	  sk1 = meanP(1)/meanP(0);
274 275 276 277
 	  if(sk1 == 0.0)
	    stangle =1.0e12;
	  else
	    stangle = abs(1/sk1);
gsell's avatar
gsell committed
278
	}else if (meanP(0) == 0.0 ){
279 280 281 282 283
	  sk2 = - A_m/B_m;
 	  if(sk2 == 0.0)
	    stangle =1.0e12;
	  else
	    stangle = abs(1/sk2);
gsell's avatar
gsell committed
284 285
	}else {
	  sk1 = meanP(1)/meanP(0);
286
	  sk2 = - A_m/B_m;
287
	  stangle = abs(( sk1-sk2 )/(1 + sk1*sk2));
gsell's avatar
gsell committed
288 289
	}
	double lstep = (sqrt(1.0-1.0/(1.0+dot(meanP, meanP))) * Physics::c) * tstep*1.0e-6; // [mm]
290
	double Swidth = lstep /  sqrt( 1+1/stangle/stangle ) * 1.2;
gsell's avatar
gsell committed
291 292 293 294 295 296 297
	setGeom(Swidth);

        for(unsigned int i = 0; i < tempnum; ++i) {
            if(bunch.PType[i] == 0) {
	      pflag = checkPoint(bunch.R[i](0), bunch.R[i](1));
	      if(pflag != 0) {
		  // dist1 > 0, right hand, dt > 0; dist1 < 0, left hand, dt < 0
298
		  double dist1 = (A_m*bunch.R[i](0)+B_m*bunch.R[i](1)+C_m)/R_m/1000.0;
gsell's avatar
gsell committed
299
		  double k1, k2, tangle = 0.0;
300
		  if ( B_m == 0.0 ){
gsell's avatar
gsell committed
301
		    k1 = bunch.P[i](1)/bunch.P[i](0);
302 303 304 305
		    if (k1 == 0.0)
		      tangle = 1.0e12;
		    else
		      tangle = abs(1/k1);
gsell's avatar
gsell committed
306
		  }else if (bunch.P[i](0) == 0.0 ){
307 308 309 310 311
		    k2 = -A_m/B_m;
		    if (k2 == 0.0)
		      tangle = 1.0e12;
		    else
		      tangle = abs(1/k2);
gsell's avatar
gsell committed
312 313
		  }else {
		    k1 = bunch.P[i](1)/bunch.P[i](0);
314
		    k2 = -A_m/B_m;
315
		    tangle = abs(( k1-k2 )/(1 + k1*k2));
gsell's avatar
gsell committed
316
		  }
317
		  double dist2 = dist1 * sqrt( 1+1/tangle/tangle );
gsell's avatar
gsell committed
318 319 320 321
		  double dt = dist2/(sqrt(1.0-1.0/(1.0 + dot(bunch.P[i], bunch.P[i]))) * Physics::c)*1.0e9;
		  strippoint(0) = (B_m*B_m*bunch.R[i](0) - A_m*B_m*bunch.R[i](1)-A_m*C_m)/(R_m*R_m);
		  strippoint(1) = (A_m*A_m*bunch.R[i](1) - A_m*B_m*bunch.R[i](0)-B_m*C_m)/(R_m*R_m);
		  strippoint(2) = bunch.R[i](2);
322 323
		  lossDs_m->addParticle_time(strippoint, bunch.P[i], bunch.ID[i], t+dt, turnnumber);

gsell's avatar
gsell committed
324 325
		  if (stop_m) {
		    bunch.Bin[i] = -1;
326
		    flagNeedUpdate = true;
gsell's avatar
gsell committed
327 328 329 330 331 332 333 334 335 336 337 338
		  }else{
                    bunch.M[i] = opmass_m;
                    bunch.Q[i] = opcharge_m * q_e;
                    bunch.PType[i] = 1;

                    //create a new particle
                    bunch.create(1);
                    bunch.R[tempnum+count] = bunch.R[i];
                    bunch.P[tempnum+count] = bunch.P[i];
                    bunch.Q[tempnum+count] = bunch.Q[i];
                    bunch.M[tempnum+count] = bunch.M[i];

339
                    // once the particle is stripped, change PType from 0 to 1 as a flag so as to avoid repetitive stripping.
gsell's avatar
gsell committed
340 341 342 343 344 345
                    bunch.PType[tempnum+count] = 1;

                    if(bunch.weHaveBins())
                        bunch.Bin[bunch.getLocalNum()-1] = bunch.Bin[i];

                    // change charge and mass of PartData when the reference particle hits the stripper.
346
                    if(bunch.ID[i] == 0)
gsell's avatar
gsell committed
347
		      flagresetMQ = true;
348

gsell's avatar
gsell committed
349 350 351 352 353 354 355 356
		    count++;
                    flagNeedUpdate = true;
		  }
	      }
            }
        }
    }
    reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
357

gsell's avatar
gsell committed
358 359 360 361 362 363 364 365
    if(!stop_m){
      reduce(&flagresetMQ, &flagresetMQ + 1, &flagresetMQ, OpBitwiseOrAssign());
      if(flagresetMQ){
	bunch.resetM(opmass_m * 1.0e9); // GeV -> eV
	bunch.resetQ(opcharge_m);     // elementary charge
      }
    }

366
    if(flagNeedUpdate) lossDs_m->save_time(getName());
gsell's avatar
gsell committed
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396
    return flagNeedUpdate;
}


void Stripper::getDimensions(double &zBegin, double &zEnd) const {
    zBegin = position_m - 0.005;
    zEnd = position_m + 0.005;
}

const string &Stripper::getType() const {
    static const string type("Stripper");
    return type;
}


int Stripper::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)

}