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