[281] | 1 | /* * * * * * * * * * * * * * * * * * * * * * * * * * * *
|
---|
| 2 | * *
|
---|
| 3 | * --<--<-- A fast simulator --<--<-- *
|
---|
| 4 | * / --<--<-- of particle --<--<-- *
|
---|
| 5 | * ----HECTOR----< *
|
---|
| 6 | * \ -->-->-- transport through -->-->-- *
|
---|
| 7 | * -->-->-- generic beamlines -->-->-- *
|
---|
| 8 | * *
|
---|
| 9 | * JINST 2:P09005 (2007) *
|
---|
| 10 | * X Rouby, J de Favereau, K Piotrzkowski (CP3) *
|
---|
| 11 | * http://www.fynu.ucl.ac.be/hector.html *
|
---|
| 12 | * *
|
---|
| 13 | * Center for Cosmology, Particle Physics and Phenomenology *
|
---|
| 14 | * Universite catholique de Louvain *
|
---|
| 15 | * Louvain-la-Neuve, Belgium *
|
---|
| 16 | * *
|
---|
| 17 | * * * * * * * * * * * * * * * * * * * * * * * * * * * */
|
---|
[3] | 18 |
|
---|
| 19 | /// \file H_BeamParticle.cc
|
---|
| 20 | /// \brief Class aiming at simulating a particle in the LHC beam
|
---|
| 21 |
|
---|
| 22 | // from IP to RP, with emission of a photon of defined energy and Q.
|
---|
| 23 | // Units : angles [rad], distances [m], energies [GeV], c=[1].
|
---|
| 24 | // !!! no comment statement at the end of a #define line !!!
|
---|
| 25 |
|
---|
| 26 | // c++ #includes
|
---|
| 27 | #include <iostream>
|
---|
| 28 | #include <iomanip>
|
---|
| 29 | #include <vector>
|
---|
| 30 | #include <string>
|
---|
| 31 | #include <cmath>
|
---|
| 32 |
|
---|
| 33 | // ROOT #includes
|
---|
| 34 | #include "H_Parameters.h"
|
---|
| 35 | #include "TRandom.h"
|
---|
| 36 | #include "TVectorD.h"
|
---|
| 37 | #ifdef _include_pythia_
|
---|
| 38 | #include "TPythia6.h"
|
---|
| 39 | #endif
|
---|
| 40 |
|
---|
| 41 | // local #includes
|
---|
| 42 | #include "H_OpticalElement.h"
|
---|
| 43 | #include "H_BeamParticle.h"
|
---|
| 44 | #include "H_Drift.h"
|
---|
| 45 |
|
---|
| 46 | using namespace std;
|
---|
| 47 |
|
---|
[281] | 48 | void H_BeamParticle::init() { // for more efficiency, put the objects away from init!
|
---|
[3] | 49 | mp = MP;
|
---|
| 50 | qp = QP;
|
---|
| 51 | fx = 0;
|
---|
| 52 | fy = 0;
|
---|
| 53 | thx = 0;
|
---|
| 54 | thy = 0;
|
---|
| 55 | fs = 0;
|
---|
| 56 | energy = BE;
|
---|
| 57 | hasstopped = false;
|
---|
| 58 | hasemitted = false;
|
---|
| 59 | isphysical = true;
|
---|
| 60 | addPosition(fx,thx,fy,thy,fs);
|
---|
| 61 | stop_position = new TVectorD(LENGTH_VEC);
|
---|
| 62 | for (int i=0; i<LENGTH_VEC; i++) (*stop_position)[i] = -1;
|
---|
| 63 | stop_element = 0;
|
---|
| 64 | }
|
---|
| 65 |
|
---|
[281] | 66 | H_BeamParticle::H_BeamParticle() {
|
---|
[3] | 67 | init();
|
---|
| 68 | }
|
---|
| 69 |
|
---|
[281] | 70 | H_BeamParticle::H_BeamParticle(const H_BeamParticle& p):
|
---|
| 71 | mp(p.mp), qp(p.qp), fs(p.fs), fx(p.fx), fy(p.fy), thx(p.thx), thy(p.thy),
|
---|
| 72 | energy(p.energy), hasstopped(p.hasstopped), hasemitted(p.hasemitted),
|
---|
| 73 | isphysical(p.isphysical), stop_position(new TVectorD(*(p.stop_position))),
|
---|
| 74 | stop_element(0), positions(p.positions) {
|
---|
[3] | 75 | if(p.hasstopped) stop_element = new H_OpticalElement(*(p.stop_element));
|
---|
| 76 | }
|
---|
| 77 |
|
---|
| 78 | H_BeamParticle::H_BeamParticle(const double mass, const double charge) {
|
---|
| 79 | init();
|
---|
| 80 | mp = mass;
|
---|
| 81 | qp = (mass==0) ? 0 : charge;
|
---|
| 82 | /// rejects particles with mass = 0 and charge != 0
|
---|
| 83 | }
|
---|
| 84 |
|
---|
| 85 | H_BeamParticle& H_BeamParticle::operator=(const H_BeamParticle& p) {
|
---|
| 86 | if(this==&p) return *this;
|
---|
| 87 | mp = p.mp;
|
---|
| 88 | qp = p.qp;
|
---|
| 89 | fx = p.fx;
|
---|
| 90 | fy = p.fy;
|
---|
| 91 | thx = p.thx;
|
---|
| 92 | thy = p.thy;
|
---|
| 93 | fs = p.fs;
|
---|
| 94 | energy = p.energy;
|
---|
| 95 | hasstopped = p.hasstopped;
|
---|
| 96 | hasemitted = p.hasemitted;
|
---|
| 97 | isphysical = p.isphysical;
|
---|
| 98 | stop_position = new TVectorD(*(p.stop_position));
|
---|
| 99 | if(p.hasstopped) stop_element = new H_OpticalElement(*(p.stop_element));
|
---|
[281] | 100 | else stop_element = 0;
|
---|
[3] | 101 | positions = p.positions;
|
---|
| 102 | return *this;
|
---|
| 103 | }
|
---|
| 104 |
|
---|
[281] | 105 | const bool H_BeamParticle::stopped(const H_AbstractBeamLine * beamline) {
|
---|
[3] | 106 | vector<TVectorD>::const_iterator position_i;
|
---|
| 107 | for(position_i = positions.begin(); position_i < positions.end()-1; position_i++) {
|
---|
| 108 | const unsigned int pos = position_i-positions.begin();
|
---|
| 109 | if(beamline->getElement(pos)->getAperture()->getType()!=NONE) {
|
---|
| 110 | bool has_passed_entrance = beamline->getElement(pos)->isInside((*position_i)[INDEX_X],(*position_i)[INDEX_Y]);
|
---|
| 111 | bool has_passed_exit = beamline->getElement(pos)->isInside((*(position_i+1))[INDEX_X],(*(position_i+1))[INDEX_Y]);
|
---|
| 112 | // here we should distinguish between particles passing the input or not.
|
---|
| 113 | // - particles not passing the input are logically stopped at the input position. period.
|
---|
| 114 | // - particles passing the input but not the output are stopped somewhere in the element
|
---|
| 115 | // - particles passing both the input and the output could have been stopped somewhere inside too (less likely)
|
---|
| 116 | if(!has_passed_entrance) {
|
---|
| 117 | if(VERBOSE) cout<<"Stopped at the entrance of "<<beamline->getElement(pos)->getName();
|
---|
| 118 | hasstopped=true;
|
---|
| 119 | stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos));
|
---|
| 120 | *stop_position = *position_i;
|
---|
| 121 | if(VERBOSE) cout<<" at s = "<<(*stop_position)[4]<<endl;
|
---|
| 122 | return hasstopped;
|
---|
| 123 | } else if (!has_passed_exit) {
|
---|
| 124 | if(VERBOSE) cout<<"Stopped inside "<<beamline->getElement(pos)->getName();
|
---|
| 125 | hasstopped=true;
|
---|
| 126 | stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos));
|
---|
| 127 | // this should be computed using the element-based method "H_OpticalElement::getHitPosition" (caution : always nonlinear).
|
---|
| 128 | *stop_position = beamline->getElement(pos)->getHitPosition(*position_i,BE-energy,mp,qp);
|
---|
| 129 | if(VERBOSE) cout<<" at s = "<<(*stop_position)[4]<<endl;
|
---|
| 130 | return hasstopped;
|
---|
| 131 | }
|
---|
| 132 | /*
|
---|
| 133 | if(!(has_passed_entrance && has_passed_exit)) {
|
---|
| 134 | if(VERBOSE) cout<<"particle stopped at "<<(beamline->getElement(pos))->getName();
|
---|
| 135 | if(VERBOSE) cout<<" (s = "<<(*position_i)[4] << ")" << endl;
|
---|
| 136 | hasstopped=true;
|
---|
| 137 | stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos));
|
---|
| 138 | *stop_position = *position_i;
|
---|
| 139 | return hasstopped;
|
---|
| 140 | } // if
|
---|
| 141 | */
|
---|
| 142 | } // if
|
---|
| 143 | } // for
|
---|
| 144 | return hasstopped;
|
---|
| 145 | }
|
---|
| 146 |
|
---|
| 147 | void H_BeamParticle::addPosition(const double x, const double tx, const double y, const double ty, const double s) {
|
---|
| 148 | // [x] = [y] = m ; [tx] = [ty] = rad; [s] = m
|
---|
| 149 | double xys[LENGTH_VEC];
|
---|
| 150 | xys[INDEX_X]=x;
|
---|
| 151 | xys[INDEX_TX]=tx;
|
---|
| 152 | xys[INDEX_Y]=y;
|
---|
| 153 | xys[INDEX_TY]=ty;
|
---|
| 154 | xys[INDEX_S]=s;
|
---|
| 155 |
|
---|
| 156 | TVectorD temp_vec(LENGTH_VEC,xys);
|
---|
| 157 | positions.push_back(temp_vec);
|
---|
| 158 | }
|
---|
| 159 |
|
---|
| 160 |
|
---|
[281] | 161 |
|
---|
| 162 | void H_BeamParticle::smearPos(const double dx,const double dy, TRandom* r) {
|
---|
[3] | 163 | // the beam is centered on (fx,fy) at IP
|
---|
[281] | 164 | fx = r->Gaus(fx,dx);
|
---|
| 165 | fy = r->Gaus(fy,dy);
|
---|
| 166 | positions.clear();
|
---|
| 167 | addPosition(fx,thx,fy,thy,fs);
|
---|
| 168 | return;
|
---|
[3] | 169 | }
|
---|
| 170 |
|
---|
[281] | 171 | void H_BeamParticle::smearAng(const double tx, const double ty, TRandom* r) {
|
---|
[3] | 172 | // the beam transverse direction is centered on (thx,thy) at IP
|
---|
[281] | 173 | thx = r->Gaus(thx,tx);
|
---|
| 174 | thy = r->Gaus(thy,ty);
|
---|
[3] | 175 | positions.clear();
|
---|
| 176 | addPosition(fx,thx,fy,thy,fs);
|
---|
| 177 | return;
|
---|
| 178 | }
|
---|
| 179 |
|
---|
[281] | 180 | void H_BeamParticle::smearE(const double erre, TRandom* r) {
|
---|
| 181 | energy = r->Gaus(energy,erre);
|
---|
[3] | 182 | return;
|
---|
| 183 | }
|
---|
| 184 |
|
---|
[281] | 185 | void H_BeamParticle::smearS(const double errs, TRandom* r) {
|
---|
| 186 | fs= r->Gaus(fs,errs);
|
---|
[3] | 187 | positions.clear();
|
---|
| 188 | addPosition(fx,thx,fy,thy,fs);
|
---|
| 189 | return;
|
---|
| 190 | }
|
---|
| 191 |
|
---|
| 192 |
|
---|
| 193 | void H_BeamParticle::set4Momentum(const double px, const double py, const double pz, const double ene) {
|
---|
| 194 | /// @param px, py, pz, ene is \f$ (\vec p , E) [GeV]\f$
|
---|
| 195 | ///
|
---|
| 196 | /// Clears the H_BeamParticle::positions vector.
|
---|
| 197 | positions.clear();
|
---|
| 198 | if(pz==0) {
|
---|
| 199 | cout<<" ERROR in H_BeamParticle::set4Momentum : no momentum in the beamline direction !"<<endl;
|
---|
| 200 | return;
|
---|
| 201 | }
|
---|
| 202 | thx = thx + URAD*atan(px/pz);
|
---|
| 203 | thy = thy + URAD*atan(py/pz);
|
---|
| 204 | energy = ene;
|
---|
| 205 | positions.clear();
|
---|
| 206 | addPosition(fx,thx,fy,thy,fs);
|
---|
| 207 | return;
|
---|
| 208 | }
|
---|
| 209 |
|
---|
| 210 | void H_BeamParticle::set4Momentum(const TLorentzVector& pmu) {
|
---|
| 211 | /// @param pmu is the particle 4-momentum \f$ p^\mu \f$
|
---|
| 212 | ///
|
---|
| 213 | /// Clears the H_BeamParticle::positions vector.
|
---|
| 214 | set4Momentum(pmu.Px(), pmu.Py(), pmu.Pz(), pmu.E());
|
---|
| 215 | }
|
---|
| 216 |
|
---|
| 217 | void H_BeamParticle::setE(const double ene) {
|
---|
| 218 | energy = ene;
|
---|
| 219 | return;
|
---|
| 220 | }
|
---|
| 221 |
|
---|
| 222 | void H_BeamParticle::setPosition(const double x, const double y, const double tx, const double ty, const double s) {
|
---|
| 223 | /// @param x, y are the transverse positions in \f$ \mu \f$ m
|
---|
| 224 | /// @param tx, ty are the angles in \f$ \mu \f$ rad
|
---|
| 225 | /// @param s is the longitudinal coordinate in m
|
---|
| 226 | // clear positions and sets the initial one.
|
---|
| 227 | fx=x;
|
---|
| 228 | fy=y;
|
---|
| 229 | thx=tx;
|
---|
| 230 | thy=ty;
|
---|
| 231 | fs = s;
|
---|
| 232 | positions.clear();
|
---|
| 233 | addPosition(fx,thx,fy,thy,s);
|
---|
| 234 |
|
---|
| 235 | return;
|
---|
| 236 | }
|
---|
| 237 |
|
---|
| 238 | const H_OpticalElement* H_BeamParticle::getStoppingElement() const{
|
---|
| 239 | if(hasstopped) return stop_element;
|
---|
| 240 | else { H_OpticalElement * dummy_el = new H_Drift("",0,0); return dummy_el;}
|
---|
| 241 | }
|
---|
| 242 |
|
---|
| 243 | void H_BeamParticle::emitGamma(const double gee, const double gq2, const double phimin, const double phimax) {
|
---|
| 244 | /// @param gee = \f$ E_{\gamma} \f$ is the photon energy
|
---|
| 245 | /// @param gq2 = Q < 0 is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$
|
---|
| 246 | /// @param phimin : lower bound for \f$ phi \f$
|
---|
| 247 | /// @param phimax : higher bound for \f$ phi \f$
|
---|
| 248 |
|
---|
| 249 | if(gq2==0) {
|
---|
| 250 | if(VERBOSE) cout<<"No virtuality : only energy has changed"<<endl;
|
---|
| 251 | setE(energy-gee);
|
---|
| 252 | return;
|
---|
| 253 | }
|
---|
| 254 |
|
---|
| 255 | double m1 = mp;
|
---|
| 256 |
|
---|
| 257 | double e1 = energy , e2 = energy - gee; // particle energy : before (1) / after (2)
|
---|
| 258 |
|
---|
| 259 | double p1 = sqrt(pow(e1,2) - pow(m1,2)), p2 = sqrt(pow(e2,2) - pow(m1,2)); // particle momentum : before (1) / after (2)
|
---|
| 260 | double q2min = pow(gee,2) - pow(p1+p2,2); // lower bound from kinematics E^2 - (p1 + p2)^2
|
---|
| 261 | double q2max = -2 * pow(m1*gee/(p1+p2),2) * (1 + (pow(e1,2) + pow(e2,2) -pow(m1,2))/(e1*e2 + p1*p2) );
|
---|
| 262 | // upper bound from kinematics ; E^2 - (p1-p2)^2; is bad for numerical computations
|
---|
| 263 |
|
---|
| 264 |
|
---|
| 265 | // if q2min < q2 < q2max is NOT true, there will be mathematical problems (like cos(eta) > 1).
|
---|
| 266 | // so if q2 > q2max, we force q2 = q2max (-> cos(eta) = 1)
|
---|
| 267 | // and if q2 < q2min, we force q2 = q2min (-> cos(eta) = 1)
|
---|
| 268 | // BUT the user knows something was wrong with the value of "H_BeamParticle::isphysical"
|
---|
| 269 | const double q2 = (gq2 > q2max ) ? q2max : (gq2 < q2min) ? q2min : gq2;
|
---|
| 270 |
|
---|
| 271 | if( (gq2>q2max) || (gq2<q2min)) {
|
---|
| 272 | if(VERBOSE) cout<<"<H_BeamParticle> WARNING : Non physical particle ! Q2 (" << q2 << ") and E ("<<gee << ") are not compatible." << endl;
|
---|
| 273 | isphysical = false;
|
---|
| 274 | }
|
---|
| 275 |
|
---|
| 276 | if(hasemitted) { cout<<"<H_BeamParticle> WARNING : particle has already emitted at least one gamma !"<<endl;}
|
---|
| 277 | hasemitted = true;
|
---|
| 278 | energy = energy - gee;
|
---|
| 279 | // gkk is k
|
---|
| 280 | double gkk = sqrt(pow(gee,2)-q2);
|
---|
| 281 | // eta is the angle between gamma and initial direction of the gamma-emitting particle
|
---|
| 282 | // ceta = cos(eta) and seta = sin(eta)
|
---|
| 283 |
|
---|
| 284 | double ceta = sqrt( pow(mp/p1,2) + 1 ) * sqrt( q2/pow(gkk,2) + 1 ) - q2/(2*p1*gkk);
|
---|
| 285 | double seta = sqrt(1 - ceta*ceta);
|
---|
| 286 | // theta is the angle between particle and beam
|
---|
| 287 | double theta = URAD*atan(seta/(BE/gkk - ceta));
|
---|
| 288 | double phi = phimin + gRandom->Uniform(phimax-phimin);
|
---|
| 289 | thx = thx + theta*cos(phi);
|
---|
| 290 | thy = thy - theta*sin(phi);
|
---|
| 291 |
|
---|
| 292 | // caution : emitting a photon erases all known positions !
|
---|
| 293 | positions.clear();
|
---|
| 294 | addPosition(fx,thx,fy,thy,fs);
|
---|
| 295 | return;
|
---|
| 296 | }
|
---|
| 297 |
|
---|
| 298 | void H_BeamParticle::doInelastic() {
|
---|
| 299 | #ifdef _include_pythia_
|
---|
| 300 | // if(!gROOT->GetClass("TPythia6")) {
|
---|
| 301 | // gROOT->Reset();
|
---|
| 302 | // gSystem->Load("libPythia6");
|
---|
| 303 | // gSystem->Load("libEG");
|
---|
| 304 | // gSystem->Load("libEGPythia6");
|
---|
| 305 | // }
|
---|
| 306 |
|
---|
| 307 | TPythia6 gen;
|
---|
| 308 | // select AB -> AX process
|
---|
| 309 | gen.SetMSEL(0);
|
---|
| 310 | gen.SetMSUB(93,1);
|
---|
| 311 | // no showers/decays
|
---|
| 312 | gen.SetMSTP(111,0);
|
---|
| 313 | // no printouts
|
---|
| 314 | gen.SetMSTP(122,0);
|
---|
| 315 | gen.SetMSTU(12,0);
|
---|
| 316 | // generator initialization
|
---|
| 317 | gen.Initialize("CMS","p","p",14000);
|
---|
| 318 | // event generation
|
---|
| 319 | gen.GenerateEvent();
|
---|
| 320 | // list particles
|
---|
| 321 | // gen.Pylist(1);
|
---|
| 322 | thx = thx + URAD*atan(gen.GetP(5,1)/gen.GetP(5,3));
|
---|
| 323 | thy = thy + URAD*atan(gen.GetP(5,2)/gen.GetP(5,3));
|
---|
| 324 | energy = gen.GetP(5,4);
|
---|
| 325 | positions.clear();
|
---|
| 326 | addPosition(fx,thx,fy,thy,fs);
|
---|
| 327 | #endif
|
---|
| 328 | return;
|
---|
| 329 | }
|
---|
| 330 |
|
---|
[281] | 331 | std::ostream& operator<< (std::ostream& os, const H_BeamParticle& p) {
|
---|
| 332 | os << " M = " << p.getM() << "GeV ";
|
---|
| 333 | os << " Q = " << p.getQ() << "e";
|
---|
| 334 | os << " fx = " << p.getX() << "m ";
|
---|
| 335 | os << " fy = " << p.getY() << "m ";
|
---|
| 336 | os << " thx = " << p.getTX() << "rad ";
|
---|
| 337 | os << " thy = " << p.getTY() << "rad ";
|
---|
| 338 | os << endl;
|
---|
| 339 | return os;
|
---|
[3] | 340 | }
|
---|
| 341 |
|
---|
| 342 | /// The phase space vector is (x,x',y,y',E)
|
---|
| 343 | /// [x] = [y] = meters
|
---|
| 344 | /// [x'] = [y'] = 1 with x' = dx/ds = tan (thetaX)
|
---|
| 345 | /// [E] = GeV
|
---|
| 346 | const TMatrixD * H_BeamParticle::getV() const {
|
---|
| 347 | double vec[MDIM] = {fx/URAD, tan(thx/URAD), fy/URAD, tan(thy/URAD),energy};
|
---|
| 348 | TMatrixD * mat = new TMatrixD(1,MDIM,vec);
|
---|
| 349 | return mat;
|
---|
| 350 | }
|
---|
| 351 |
|
---|
| 352 | void H_BeamParticle::printV() const {
|
---|
| 353 | TMatrixD X(*getV());
|
---|
| 354 | cout << " x = " << (X.GetMatrixArray())[0] << "m ";
|
---|
| 355 | cout << " x' = " << (X.GetMatrixArray())[1] << " ";
|
---|
| 356 | cout << " y = " << (X.GetMatrixArray())[2] << "m ";
|
---|
| 357 | cout << " y' = " << (X.GetMatrixArray())[3] << " ";
|
---|
| 358 | cout << endl;
|
---|
| 359 | return;
|
---|
| 360 | }
|
---|
| 361 |
|
---|
| 362 | void H_BeamParticle::propagate(const double position) {
|
---|
| 363 | /// @param position is the s coordinate in m to reach
|
---|
| 364 | /// Does not propagate if position is in the middle of an otics element of the beamline.
|
---|
| 365 | vector<TVectorD>::const_iterator position_i = positions.begin();
|
---|
| 366 | double l = 0.;
|
---|
| 367 | if(position != fs) { // avoid repeating the computation if already done at this position
|
---|
| 368 | if(position == (*position_i)[INDEX_S]) {
|
---|
| 369 | fs = position;
|
---|
| 370 | fx = (*position_i)[INDEX_X];
|
---|
| 371 | fy = (*position_i)[INDEX_Y];
|
---|
| 372 | thx= (*position_i)[INDEX_TX];
|
---|
| 373 | thy= (*position_i)[INDEX_TY];
|
---|
| 374 | return;
|
---|
| 375 | } else
|
---|
| 376 | for(position_i = positions.begin(); position_i < positions.end(); position_i++) {
|
---|
| 377 | if((*position_i)[INDEX_S]>=position) {
|
---|
| 378 | if(position_i==positions.begin()) {
|
---|
| 379 | if(VERBOSE) cout<<"<H_BeamParticle> ERROR : non reachable value"<<endl;
|
---|
| 380 | return;
|
---|
| 381 | }
|
---|
| 382 | l = (*position_i)[INDEX_S] - (*(position_i-1))[INDEX_S];
|
---|
| 383 | if(l==0) {
|
---|
| 384 | if(VERBOSE) cout<<"<H_BeamParticle> WARNING : no luck in choosing position, no propagation done"<<endl;
|
---|
| 385 | return;
|
---|
| 386 | }
|
---|
| 387 | fs = position;
|
---|
| 388 | fx = (*(position_i-1))[INDEX_X] + (position-(*(position_i-1))[INDEX_S])*((*position_i)[INDEX_X] - (*(position_i-1))[INDEX_X])/l;
|
---|
| 389 | fy = (*(position_i-1))[INDEX_Y] + (position-(*(position_i-1))[INDEX_S])*((*position_i)[INDEX_Y] - (*(position_i-1))[INDEX_Y])/l;
|
---|
| 390 | thx = (*(position_i-1))[INDEX_TX];
|
---|
| 391 | thy = (*(position_i-1))[INDEX_TY];
|
---|
| 392 | return;
|
---|
| 393 | }
|
---|
| 394 | }
|
---|
| 395 | position_i = positions.begin();
|
---|
| 396 | cout << "Desired position is : " << position << " & positions.begin() is " << (*position_i)[INDEX_S] << endl;
|
---|
| 397 | cout<<"<H_BeamParticle> ERROR : position not reachable"<<endl;
|
---|
| 398 | return;
|
---|
| 399 | }
|
---|
| 400 | }
|
---|
| 401 |
|
---|
| 402 | /// Caution : do not use this method (obsolete) !!!
|
---|
| 403 | void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const H_OpticalElement * element) {
|
---|
| 404 | TMatrixD X(*getV());
|
---|
| 405 | X *= beam->getPartialMatrix(element);
|
---|
| 406 | fx = URAD*(X.GetMatrixArray())[0];
|
---|
| 407 | thx = URAD*atan((X.GetMatrixArray())[1]);
|
---|
| 408 | fy = URAD*(X.GetMatrixArray())[2];
|
---|
| 409 | thy = URAD*atan((X.GetMatrixArray())[3]);
|
---|
| 410 | return;
|
---|
| 411 | }
|
---|
| 412 |
|
---|
[216] | 413 | void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const string& el_name) {
|
---|
[3] | 414 | propagate(beam->getElement(el_name)->getS()+beam->getElement(el_name)->getLength());
|
---|
| 415 | return;
|
---|
| 416 | }
|
---|
| 417 |
|
---|
| 418 | void H_BeamParticle::propagate(const H_AbstractBeamLine * beam) {
|
---|
| 419 | TMatrixD X(*getV());
|
---|
| 420 | X *= beam->getBeamMatrix();
|
---|
| 421 | fx = URAD*(X.GetMatrixArray())[0];
|
---|
| 422 | thx = URAD*atan((X.GetMatrixArray())[1]);
|
---|
| 423 | fy = URAD*(X.GetMatrixArray())[2];
|
---|
| 424 | thy = URAD*atan((X.GetMatrixArray())[3]);
|
---|
| 425 | return;
|
---|
| 426 | }
|
---|
| 427 |
|
---|
| 428 | void H_BeamParticle::showPositions() const{
|
---|
| 429 | vector<TVectorD>::const_iterator position_i;
|
---|
| 430 | TVectorD temp_vec(LENGTH_VEC);
|
---|
| 431 |
|
---|
| 432 | for(position_i = positions.begin(); position_i < positions.end(); position_i++) {
|
---|
| 433 | cout << "Vector (x,y,s) = (" << (*position_i)[INDEX_X] << ", " << (*position_i)[INDEX_Y] << ", " << (*position_i)[INDEX_S] << "). " << endl;
|
---|
| 434 | }
|
---|
| 435 | return ;
|
---|
| 436 | }
|
---|
| 437 |
|
---|
| 438 | TGraph * H_BeamParticle::getPath(const int x_or_y, const int color) const{
|
---|
| 439 | /// @param x_or_y = 0(1) draws the x(y) component;
|
---|
| 440 |
|
---|
| 441 | const int N = (int) positions.size();
|
---|
| 442 | int mycolor = color;
|
---|
| 443 | if(N<2) cout<<"<H_BeamParticle> WARNING : particle positions not calculated : please run computePath"<<endl;
|
---|
| 444 | double * s = new double[N], * graph = new double[N];
|
---|
| 445 |
|
---|
| 446 | int index;
|
---|
| 447 | if(x_or_y==0) {index = INDEX_X;} else {index = INDEX_Y;}
|
---|
| 448 |
|
---|
| 449 | vector<TVectorD>::const_iterator position_i;
|
---|
| 450 | for(position_i = positions.begin(); position_i < positions.end(); position_i++) {
|
---|
| 451 | graph[(int)(position_i-positions.begin())] = (*position_i)[index];
|
---|
| 452 | s[(int)(position_i-positions.begin())] = (*position_i)[INDEX_S];
|
---|
| 453 | }
|
---|
| 454 |
|
---|
| 455 | TGraph * ppath = new TGraph(N,s,graph);
|
---|
| 456 | ppath->SetLineColor(mycolor);
|
---|
| 457 | delete [] s;
|
---|
| 458 | delete [] graph;
|
---|
| 459 | return ppath;
|
---|
| 460 | }
|
---|
| 461 |
|
---|
| 462 | // should be removed later, to keep only computePath(const H_AbstractBeamLine & , const bool)
|
---|
| 463 | void H_BeamParticle::computePath(const H_AbstractBeamLine * beam, const bool NonLinear) {
|
---|
| 464 | TMatrixD temp_mat(MDIM,MDIM);
|
---|
| 465 | double temp_x, temp_y, temp_s, temp_tx, temp_ty;
|
---|
| 466 |
|
---|
| 467 | temp_x = (positions.front())[INDEX_X];
|
---|
| 468 | temp_tx = (positions.front())[INDEX_TX];
|
---|
| 469 | temp_y = (positions.front())[INDEX_Y];
|
---|
| 470 | temp_ty = (positions.front())[INDEX_TY];
|
---|
| 471 | temp_s = (positions.front())[INDEX_S];
|
---|
| 472 |
|
---|
| 473 | double vec[MDIM] = {temp_x/URAD, tan(temp_tx/URAD), temp_y/URAD, tan(temp_ty/URAD),energy,1};
|
---|
| 474 |
|
---|
| 475 | extern bool relative_energy;
|
---|
| 476 | if(relative_energy) {
|
---|
| 477 | vec[4] = energy-BE;
|
---|
| 478 | } else {
|
---|
| 479 | vec[4] = energy;
|
---|
| 480 | }
|
---|
| 481 |
|
---|
| 482 | TMatrixD mat(1,MDIM,vec);
|
---|
| 483 |
|
---|
| 484 | const int N =beam->getNumberOfElements();
|
---|
| 485 | double xys[LENGTH_VEC];
|
---|
| 486 |
|
---|
| 487 | double energy_loss = NonLinear?BE-energy:0;
|
---|
| 488 |
|
---|
| 489 | for (int i=0; i<N; i++) {
|
---|
| 490 | const unsigned pos = i;
|
---|
| 491 | if(fs < beam->getElement(pos)->getS() && fs > beam->getElement(pos-1)->getS()) {
|
---|
| 492 | if(beam->getElement(pos-1)->getType()!=DRIFT) {
|
---|
| 493 | cout<<"Path starts inside element "<<beam->getElement(pos-1)->getName()<<endl;
|
---|
| 494 | } else {
|
---|
| 495 | cout<<"Path starts inside unnamed drift "<<endl;
|
---|
| 496 | }
|
---|
| 497 | H_OpticalElement* temp_el = beam->getElement(pos-1)->clone();
|
---|
| 498 | temp_el->setS(fs);
|
---|
| 499 | temp_el->setLength(beam->getElement(pos)->getS() - fs);
|
---|
| 500 | mat[0][0] = mat[0][0] - temp_el->getX();
|
---|
| 501 | mat[0][1] = mat[0][1] - tan(temp_el->getTX());
|
---|
| 502 | mat[0][2] = mat[0][2] - temp_el->getY();
|
---|
| 503 | mat[0][3] = mat[0][3] - tan(temp_el->getTY());
|
---|
| 504 | mat *= temp_el->getMatrix(energy_loss,mp,qp);
|
---|
| 505 | mat[0][0] = mat[0][0] + temp_el->getX();
|
---|
| 506 | mat[0][1] = mat[0][1] + tan(temp_el->getTX());
|
---|
| 507 | mat[0][2] = mat[0][2] + temp_el->getY();
|
---|
| 508 | mat[0][3] = mat[0][3] + tan(temp_el->getTY());
|
---|
| 509 | xys[0] = mat.GetMatrixArray()[0]*URAD;
|
---|
| 510 | xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
|
---|
| 511 | xys[2] = mat.GetMatrixArray()[2]*URAD;
|
---|
| 512 | xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
|
---|
| 513 | xys[4] = temp_el->getS()+temp_el->getLength();
|
---|
| 514 | addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
|
---|
| 515 | if (temp_el) delete temp_el;
|
---|
| 516 | }
|
---|
| 517 | if(fs <= beam->getElement(pos)->getS()) {
|
---|
| 518 | mat[0][0] = mat[0][0] - beam->getElement(pos)->getX();
|
---|
| 519 | mat[0][1] = mat[0][1] - tan(beam->getElement(pos)->getTX()/URAD)*URAD;
|
---|
| 520 | mat[0][2] = mat[0][2] - beam->getElement(pos)->getY();
|
---|
| 521 | mat[0][3] = mat[0][3] - tan(beam->getElement(pos)->getTY()/URAD)*URAD;
|
---|
| 522 | mat *= beam->getElement(pos)->getMatrix(energy_loss,mp,qp);
|
---|
| 523 | mat[0][0] = mat[0][0] + beam->getElement(pos)->getX();
|
---|
| 524 | mat[0][1] = mat[0][1] + tan(beam->getElement(pos)->getTX()/URAD)*URAD;
|
---|
| 525 | mat[0][2] = mat[0][2] + beam->getElement(pos)->getY();
|
---|
| 526 | mat[0][3] = mat[0][3] + tan(beam->getElement(pos)->getTY()/URAD)*URAD;
|
---|
| 527 | xys[0] = mat.GetMatrixArray()[0]*URAD;
|
---|
| 528 | xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
|
---|
| 529 | xys[2] = mat.GetMatrixArray()[2]*URAD;
|
---|
| 530 | xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
|
---|
| 531 | xys[4] = beam->getElement(pos)->getS()+beam->getElement(pos)->getLength();
|
---|
| 532 | addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
|
---|
| 533 | }
|
---|
| 534 | }
|
---|
| 535 | }
|
---|
| 536 |
|
---|
| 537 | // part about non-ip particle is not ready yet. use the above method in the meantime
|
---|
| 538 | void H_BeamParticle::computePath(const H_AbstractBeamLine & beam, const bool NonLinear) {
|
---|
| 539 | TMatrixD temp_mat(MDIM,MDIM);
|
---|
| 540 | double temp_x, temp_y, temp_s, temp_tx, temp_ty;
|
---|
| 541 |
|
---|
| 542 | temp_x = (positions.front())[INDEX_X];
|
---|
| 543 | temp_tx = (positions.front())[INDEX_TX];
|
---|
| 544 | temp_y = (positions.front())[INDEX_Y];
|
---|
| 545 | temp_ty = (positions.front())[INDEX_TY];
|
---|
| 546 | temp_s = (positions.front())[INDEX_S];
|
---|
| 547 |
|
---|
| 548 | double vec[MDIM] = {temp_x/URAD, tan(temp_tx/URAD), temp_y/URAD, tan(temp_ty/URAD),energy,1};
|
---|
| 549 |
|
---|
| 550 | extern bool relative_energy;
|
---|
| 551 | if(relative_energy) {
|
---|
| 552 | vec[4] = energy-BE;
|
---|
| 553 | } else {
|
---|
| 554 | vec[4] = energy;
|
---|
| 555 | }
|
---|
| 556 |
|
---|
| 557 | TMatrixD mat(1,MDIM,vec);
|
---|
| 558 |
|
---|
| 559 | const int N =beam.getNumberOfElements();
|
---|
| 560 | double xys[LENGTH_VEC];
|
---|
| 561 |
|
---|
| 562 | double energy_loss = NonLinear?BE-energy:0;
|
---|
| 563 |
|
---|
| 564 | // modify here to allow starting at non-IP positions
|
---|
| 565 | // s is distance to IP
|
---|
| 566 | // initial position is already in positions vector ?
|
---|
| 567 | for (int i=0; i<N; i++) {
|
---|
| 568 | const unsigned pos = i;
|
---|
| 569 | // if we are inside an element, we should start by adding the action
|
---|
| 570 | // of the rest of this element
|
---|
| 571 | if(pos > 0 && fs < beam.getElement(pos)->getS() && fs > beam.getElement(pos-1)->getS()) {
|
---|
| 572 | cout<<"Path starts inside element "<<beam.getElement(pos-1)->getName()<<endl;
|
---|
| 573 | H_OpticalElement* temp_el = new H_OpticalElement(*(beam.getElement(pos-1)));
|
---|
| 574 | temp_el->setS(fs);
|
---|
| 575 | temp_el->setLength(beam.getElement(pos)->getS() - fs);
|
---|
| 576 | mat[0][0] = mat[0][0] - temp_el->getX();
|
---|
| 577 | mat[0][1] = mat[0][1] - tan(temp_el->getTX());
|
---|
| 578 | mat[0][2] = mat[0][2] - temp_el->getY();
|
---|
| 579 | mat[0][3] = mat[0][3] - tan(temp_el->getTY());
|
---|
| 580 | mat *= temp_el->getMatrix(energy_loss,mp,qp);
|
---|
| 581 | mat[0][0] = mat[0][0] + temp_el->getX();
|
---|
| 582 | mat[0][1] = mat[0][1] + tan(temp_el->getTX());
|
---|
| 583 | mat[0][2] = mat[0][2] + temp_el->getY();
|
---|
| 584 | mat[0][3] = mat[0][3] + tan(temp_el->getTY());
|
---|
| 585 | } else if(fs >= beam.getElement(pos)->getS()) {
|
---|
| 586 | mat[0][0] = mat[0][0] - beam.getElement(pos)->getX();
|
---|
| 587 | mat[0][1] = mat[0][1] - tan(beam.getElement(pos)->getTX());
|
---|
| 588 | mat[0][2] = mat[0][2] - beam.getElement(pos)->getY();
|
---|
| 589 | mat[0][3] = mat[0][3] - tan(beam.getElement(pos)->getTY());
|
---|
| 590 | mat *= beam.getElement(pos)->getMatrix(energy_loss,mp,qp);
|
---|
| 591 | mat[0][0] = mat[0][0] + beam.getElement(pos)->getX();
|
---|
| 592 | mat[0][1] = mat[0][1] + tan(beam.getElement(pos)->getTX());
|
---|
| 593 | mat[0][2] = mat[0][2] + beam.getElement(pos)->getY();
|
---|
| 594 | mat[0][3] = mat[0][3] + tan(beam.getElement(pos)->getTY());
|
---|
| 595 | }
|
---|
| 596 | xys[0] = mat.GetMatrixArray()[0]*URAD;
|
---|
| 597 | xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
|
---|
| 598 | xys[2] = mat.GetMatrixArray()[2]*URAD;
|
---|
| 599 | xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
|
---|
| 600 | xys[4] = beam.getElement(pos)->getS()+beam.getElement(pos)->getLength();
|
---|
| 601 | addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
|
---|
| 602 | fx = xys[0];
|
---|
| 603 | fy = xys[2];
|
---|
| 604 | thx = xys[1];
|
---|
| 605 | thy = xys[3];
|
---|
| 606 | }
|
---|
| 607 | }
|
---|
| 608 |
|
---|
| 609 | void H_BeamParticle::resetPath() {
|
---|
| 610 | double temp_x, temp_y, temp_s, temp_tx, temp_ty;
|
---|
| 611 |
|
---|
| 612 | temp_x = (positions.front())[INDEX_X];
|
---|
| 613 | temp_tx = (positions.front())[INDEX_TX];
|
---|
| 614 | temp_y = (positions.front())[INDEX_Y];
|
---|
| 615 | temp_ty = (positions.front())[INDEX_TY];
|
---|
| 616 | temp_s = (positions.front())[INDEX_S];
|
---|
| 617 | positions.clear();
|
---|
| 618 | addPosition(temp_x,temp_tx,temp_y,temp_ty,temp_s);
|
---|
| 619 | }
|
---|
| 620 |
|
---|
| 621 | const TVectorD * H_BeamParticle::getPosition(const int element_position) const {
|
---|
| 622 | const int N = (element_position<0)?0:(( ((unsigned int) element_position)>positions.size()-1)?positions.size()-1:element_position);
|
---|
| 623 | return &(*(positions.begin()+N));// same as "return &positions[N];", but more efficient
|
---|
| 624 | }
|
---|