~max-mcginley/third-order-lie-map/generator : revision 14

16

17

18

namespace MAUS {

19

BASEparticle::BASEparticle() : TransferMatrixGenerator("BASEparticle"), usingLattice(false) {

19

20

Matrix<double> TransferMatrixGenerator::getKineticMatrix() const {

21

Matrix<double> shift(4,4,0.0);

22

shift(1,1) = 1.0/sigma;

23

shift(2,2) = 1.0/sigma;

24

shift(3,3) = 1.0/refMomentum;

25

shift(4,4) = 1.0/refMomentum;

26

shift(3,2) = -getIntegrator()->getBField(getZIn())*lightSpeed*muonCharge/(2.e6*refMomentum);

27

shift(4,1) = getIntegrator()->getBField(getZIn())*lightSpeed*muonCharge/(2.e6*refMomentum);

28

29

Matrix<double> back(4,4,0.0);

30

back(1,1) = sigma;

31

back(2,2) = sigma;

32

back(3,3) = refMomentum;

33

back(4,4) = refMomentum;

34

back(3,2) = getIntegrator()->getBField(getZOut())*lightSpeed*muonCharge/2.e6;

35

back(4,1) = -getIntegrator()->getBField(getZOut())*lightSpeed*muonCharge/2.e6;

36

37

Matrix<double> tran = getTransferMatrix();

38

39

return back * tran * shift;

40

}

41

42

BASEparticle::BASEparticle(double zIn, double zOut, double refMom) : TransferMatrixGenerator("BASEparticle",zIn,zOut,refMom), usingLattice(false) {

20

43

getCurrentStream() << "Constructing BASEparticle" << std::endl;

21

44

matrixGen = false;

22

45

particlesRun = false;

31

54

_numberOfParticles = 0;

32

55

}

33

56

34

BASEparticle::BASEparticle(AnalyticLattice* lattice) : TransferMatrixGenerator("BASEparticle"), usingLattice(true), _lattice(lattice) {

57

BASEparticle::BASEparticle(AnalyticLattice* lattice, double zIn, double zOut, double refMom) : TransferMatrixGenerator("BASEparticle",zIn,zOut,refMom), usingLattice(true), _lattice(lattice) {

35

58

getCurrentStream() << "Constructing BASEparticle using lattice" << std::endl;

36

59

matrixGen = false;

37

60

particlesRun = false;

76

99

return v;

77

100

}

78

101

79

void BASEparticle::runParticles(double zOut) {

102

void BASEparticle::runParticles() {

80

103

if (usingLattice) {

81

104

getCurrentStream() << "Running reference particle" << std::endl;

82

105

trackedReference = new double[8];

83

106

for (int i = 0; i < 8; i++)

84

107

trackedReference[i] = inReference[i];

85

BTTracker::integrate(zOut,trackedReference,getLattice()->getFieldGroup(),BTTracker::z, stepSize, 1.0);

108

BTTracker::integrate(getZOut(),trackedReference,getLattice()->getFieldGroup(),BTTracker::z, stepSize, -1.0);

86

109

getCurrentStream() << "Running all particles" << std::endl;

87

110

trackedParticles = new double[8*numberOfParticles];

88

111

for (int i = 0; i < 8*numberOfParticles; i++)

89

112

trackedParticles[i] = inParticles[i];

90

113

for (int i = 0; i < numberOfParticles; i++)

91

BTTracker::integrate(zOut,trackedParticles + (i*8),getLattice()->getFieldGroup(),BTTracker::z, stepSize, 1.0);

114

BTTracker::integrate(getZOut(),trackedParticles + (i*8),getLattice()->getFieldGroup(),BTTracker::z, stepSize, -1.0);

92

115

particlesRun = true;

93

116

/*

94

117

for (int i = 0; i < numberOfParticles; i++) {

220

243

return retval;

221

244

}

222

245

223

void BASEparticle::generateParticles(double zIn) {

246

void BASEparticle::generateParticles() {

224

247

getCurrentStream() << "Generating " << numberOfParticles << " particles..." << std::endl;

225

generateReferenceParticle(zIn);

248

generateReferenceParticle();

226

249

if (!usingLattice) {

227

250

Json::Value retval(Json::arrayValue);

228

251

for(int i = 0; i < numberOfParticles; i++) {

235

258

double pz = std::sqrt(val["energy"].asDouble()*val["energy"].asDouble() - px*px - py*py);

236

259

val["position"]["x"] = x;

237

260

val["position"]["y"] = y;

238

val["position"]["z"] = zIn;

261

val["position"]["z"] = getZIn();

239

262

val["momentum"]["x"] = px;

240

263

val["momentum"]["y"] = py;

241

264

val["momentum"]["z"] = pz;

259

282

parts[i*8] = 0.0;

260

283

parts[i*8 + 1] = x;

261

284

parts[i*8 + 2] = y;

262

parts[i*8 + 3] = zIn;

285

parts[i*8 + 3] = getZIn();

263

286

parts[i*8 + 4] = refEnergy;

264

287

parts[i*8 + 5] = px;

265

288

parts[i*8 + 6] = py;

280

303

return *polynomialMatrix;

281

304

}

282

305

283

void BASEparticle::generateReferenceParticle(double zIn) {

306

void BASEparticle::generateReferenceParticle() {

284

307

if (usingLattice) {

285

308

double* ret = new double[8];

286

309

ret[0] = 0.0; ret[1] = 0.0; ret[2] = 0.0; ret[5] = 0.0; ret[6] = 0.0;

287

ret[3] = zIn;

310

ret[3] = getZIn();

288

311

ret[4] = refEnergy;

289

312

ret[7] = refMomentum;

290

313

inReference = ret;

292

315

else setRefPart((MAUSGeant4Manager::GetInstance())->GetReferenceParticle());

293

316

}

294

317

295

double* BASEparticle::generateRandomCoordinates() {

318

double* BASEparticle::generateRandomCoordinates(double scale) {

296

319

double* ret = (double*) malloc(sizeof(double)*4);

297

for(int i = 0; i < 4; i++) ret[i] = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)/10.0;

320

for(int i = 0; i < 4; i++) ret[i] = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*scale;

298

321

return ret;

299

322

}

300

323

301

void BASEparticle::advanceCell(double* coordsIn, double* coordsOut, double zIn, double zOut) {

324

void BASEparticle::advanceCell(double* coordsIn, double* coordsOut) {

302

325

for (int i = 0; i < 4; i++)

303

326

coordsOut[i] = 0.0;

304

327

if (!usingLattice) return;

306

329

properCoords[0] = 0.;

307

330

properCoords[1] = coordsIn[0]*sigma;

308

331

properCoords[2] = coordsIn[1]*sigma;

309

std::vector<double> A1 = integrator->getVectorPotentialAnalytic(properCoords[0],properCoords[1],zIn);

310

properCoords[3] = zIn;

332

std::vector<double> A1 = integrator->getVectorPotentialAnalytic(properCoords[0],properCoords[1],getZIn());

333

properCoords[3] = getZIn();

311

334

properCoords[4] = refEnergy;

312

335

properCoords[5] = coordsIn[2]*refMomentum - A1[0]*lightSpeed/1.0e-6;

313

336

properCoords[6] = coordsIn[3]*refMomentum - A1[1]*lightSpeed/1.0e-6;

314

337

properCoords[7] = std::sqrt(refEnergy*refEnergy - properCoords[5]*properCoords[5] - properCoords[6]*properCoords[6]);

315

BTTracker::integrate(zOut,properCoords,getLattice()->getFieldGroup(),BTTracker::z, stepSize, 1.0);

316

std::vector<double> A2 = integrator->getVectorPotentialAnalytic(properCoords[0],properCoords[1],zOut);

338

BTTracker::integrate(getZOut(),properCoords,getLattice()->getFieldGroup(),BTTracker::z, stepSize, 1.0);

339

std::vector<double> A2 = integrator->getVectorPotentialAnalytic(properCoords[0],properCoords[1],getZOut());

317

340

coordsOut[0] = properCoords[1]/sigma;

318

341

coordsOut[1] = properCoords[2]/sigma;

319

342

coordsOut[2] = (properCoords[5] + A2[0]*lightSpeed/1.0e-6)/sigma;

396

419

double py2 = v_hit2["momentum"]["y"].asDouble();

397

420

double z2 = v_hit2["position"]["z"].asDouble();

398

421

399

//conversion to canonical momenta

400

std::vector<double> A1 = integrator->getVectorPotentialAnalytic(x1,y1,z1);

401

std::vector<double> A2 = integrator->getVectorPotentialAnalytic(x2,y2,z2);

402

403

//from T m * {fundamental charge unit} -> MeV/c is lightSpeed (m s^-1) / 1,000,000

404

px1 += A1[0] * lightSpeed / 1.0e6;

405

py1 += A1[1] * lightSpeed / 1.0e6;

406

px2 += A2[0] * lightSpeed / 1.0e6;

407

py2 += A2[1] * lightSpeed / 1.0e6;

408

409

//rescaling to dimensionless MARYLIE variables

410

x1 /= sigma;

411

y1 /= sigma;

412

x2 /= sigma;

413

y2 /= sigma;

414

px1 /= refMom;

415

py1 /= refMom;

416

px2 /= refMom;

417

py2 /= refMom;

418

419

/*getCurrentStream() << "Particle IN: z = " << step1["position"]["z"] << std::endl;

420

getCurrentStream() << "x = " << x1 << " y = " << y1 << " px = " << px1 << " py = " << py1 << std::endl;

421

getCurrentStream() << "Particle OUT: z = " << step2["position"]["z"] << std::endl;

422

getCurrentStream() << "x = " << x2 << " y = " << y2 << " px = " << px2 << " py = " << py2 << "\n" << std::endl;*/

423

424

std::vector<double> phaseCoordIn(4);

425

std::vector<double> phaseCoordOut(4);

426

427

phaseCoordIn[0] = x1;

428

phaseCoordIn[1] = y1;

429

phaseCoordIn[2] = px1;

430

phaseCoordIn[3] = py1;

431

432

phaseCoordOut[0] = x2;

433

phaseCoordOut[1] = y2;

434

phaseCoordOut[2] = px2;

435

phaseCoordOut[3] = py2;

422

std::vector<double> phaseCoordIn = convertToPhaseSpaceCoords(x1,y1,px1,py1,z1,refMom,integrator);

423

std::vector<double> phaseCoordOut = convertToPhaseSpaceCoords(x2,y2,px2,py2,z2,refMom,integrator);

436

424

437

425

std::pair<std::vector<double>, std::vector<double> > retval(phaseCoordIn, phaseCoordOut);

438

426

441

429

442

430

std::pair<std::vector<double>, std::vector<double> > BASEparticle::getVirtualCoordsFromTracked(double* particlesIn, double* particlesOut, double refMom, bool thirdOrder) {

443

431

//initial kinetic variables

432

/*

444

433

double x1 = particlesIn[1];

445

434

double y1 = particlesIn[2];

446

435

double px1 = particlesIn[5];

454

443

double py2 = particlesOut[6];

455

444

double z2 = particlesOut[3];

456

445

457

//conversion to canonical momenta

458

std::vector<double> A1 = integrator->getVectorPotentialAnalytic(x1,y1,z1);

459

std::vector<double> A2 = integrator->getVectorPotentialAnalytic(x2,y2,z2);

460

461

//from T m * {fundamental charge unit} -> MeV/c is lightSpeed (m s^-1) / 1,000,000

462

px1 += A1[0] * lightSpeed / 1.0e6;

463

py1 += A1[1] * lightSpeed / 1.0e6;

464

px2 += A2[0] * lightSpeed / 1.0e6;

465

py2 += A2[1] * lightSpeed / 1.0e6;

466

467

//rescaling to dimensionless MARYLIE variables

468

x1 /= sigma;

469

y1 /= sigma;

470

x2 /= sigma;

471

y2 /= sigma;

472

px1 /= refMom;

473

py1 /= refMom;

474

px2 /= refMom;

475

py2 /= refMom;

476

477

/*getCurrentStream() << "Particle IN: z = " << step1["position"]["z"] << std::endl;

478

getCurrentStream() << "x = " << x1 << " y = " << y1 << " px = " << px1 << " py = " << py1 << std::endl;

479

getCurrentStream() << "Particle OUT: z = " << step2["position"]["z"] << std::endl;

480

getCurrentStream() << "x = " << x2 << " y = " << y2 << " px = " << px2 << " py = " << py2 << "\n" << std::endl;*/

446

447

std::vector<double> phaseCoordIn = convertToPhaseSpaceCoords(x1,y1,px1,py1,z1,refMom,integrator);

448

std::vector<double> phaseCoordOut = convertToPhaseSpaceCoords(x2,y2,px2,py2,z2,refMom,integrator);

449

450

std::pair<std::vector<double>, std::vector<double> > retval(phaseCoordIn, phaseCoordOut);

451

*/

481

452

482

453

std::vector<double> phaseCoordIn(4);

483

454

std::vector<double> phaseCoordOut(4);

455

456

phaseCoordIn[0] = particlesIn[1];

457

phaseCoordIn[1] = particlesIn[2];

458

phaseCoordIn[2] = particlesIn[5];

459

phaseCoordIn[3] = particlesIn[6];

460

461

//final kinetic variables

462

phaseCoordOut[0] = particlesOut[1];

463

phaseCoordOut[1] = particlesOut[2];

464

phaseCoordOut[2] = particlesOut[5];

465

phaseCoordOut[3] = particlesOut[6];

466

467

std::pair<std::vector<double>, std::vector<double> > retval(phaseCoordIn, phaseCoordOut);

468

469

470

return retval;

471

}

472

473

std::vector<double> BASEparticle::convertToPhaseSpaceCoords(double x, double y, double px, double py, double z, double refMom, BfieldIntegrator* integ) {

474

std::vector<double> phaseCoord(4);

475

//conversion to canonical momenta

476

std::vector<double> A = integ->getVectorPotential(x,y,z);

477

478

//from T m * {fundamental charge unit} -> MeV/c is lightSpeed (m s^-1) / 1,000,000

479

px += muonCharge * A[0] * lightSpeed / 1.0e6;

480

py += muonCharge * A[1] * lightSpeed / 1.0e6;

481

482

//rescaling to dimensionless MARYLIE variables

483

x /= sigma;

484

y /= sigma;

485

px /= refMom;

486

py /= refMom;

487

488

phaseCoord[0] = x;

489

phaseCoord[1] = y;

490

phaseCoord[2] = px;

491

phaseCoord[3] = py;

492

493

return phaseCoord;

494

}

495

496

std::vector<double> BASEparticle::convertToPhaseSpaceCoords(std::vector<double> v, double z, double refMom, BfieldIntegrator* integ) {

497

return convertToPhaseSpaceCoords(v[0],v[1],v[2],v[3],z,refMom,integ);

498

}

499

500

std::vector<double> BASEparticle::convertToKineticCoords(double x, double y, double px, double py, double z, double refMom, BfieldIntegrator* integ) {

501

std::vector<double> kinCoord(4);

502

//conversion to canonical momenta

503

std::vector<double> A = integ->getVectorPotential(x,y,z);

504

505

x *= sigma;

506

y *= sigma;

507

px *= refMom;

508

py *= refMom;

484

509

485

phaseCoordIn[0] = x1;

486

phaseCoordIn[1] = y1;

487

phaseCoordIn[2] = px1;

488

phaseCoordIn[3] = py1;

510

//from T m * {fundamental charge unit} -> MeV/c is lightSpeed (m s^-1) / 1,000,000

511

px -= muonCharge * A[0] * lightSpeed / 1.0e6;

512

py -= muonCharge * A[1] * lightSpeed / 1.0e6;

513

514

//rescaling to dimensionless MARYLIE variables

515

516

517

kinCoord[0] = x;

518

kinCoord[1] = y;

519

kinCoord[2] = px;

520

kinCoord[3] = py;

521

522

return kinCoord;

523

}

489

524

490

phaseCoordOut[0] = x2;

491

phaseCoordOut[1] = y2;

492

phaseCoordOut[2] = px2;

493

phaseCoordOut[3] = py2;

494

495

std::pair<std::vector<double>, std::vector<double> > retval(phaseCoordIn, phaseCoordOut);

496

497

return retval;

525

std::vector<double> BASEparticle::convertToKineticCoords(std::vector<double> v, double z, double refMom, BfieldIntegrator* integ) {

526

return convertToKineticCoords(v[0],v[1],v[2],v[3],z,refMom,integ);

498

527

}

499

528

500

529

/*

508

537

addTest(test);

509

538

}*/

510

539

511

Matrix<double> TransferMatrixGenerator::getFirstOrderEllipse() {

512

return EmittanceCalculator::getEllipseFromMap(getTransferMatrix());

513

}

514

515

std::pair<double,double> TransferMatrixGenerator::getActionAngleCoordinates(double x, double px) {

516

Matrix<double> ellipse = getFirstOrderEllipse();

517

std::pair<Vector<complex>, Matrix<complex> > eig = eigensystem(ellipse);

518

Vector<double> realVal = real(eig.first);

519

Vector<double> imagVal = real(eig.first);

520

Matrix<double> realOrth = real(eig.second);

521

Matrix<double> imagOrth = real(eig.second);

522

for (int i = 1; i <= 4; i++) {

523

if (imagVal(i) != 0.0) throw std::runtime_error("LieMap::getActionAngleCoordinates() - eigenvalues are not pure real");

524

for (int j = 1; j <= 4; j++) {

525

if (imagOrth(i,j) != 0.0) throw std::runtime_error("LieMap::getActionAngleCoordinates() - eigenvectors are not pure real");

526

}

540

double TransferMatrixGenerator::getEmittance(double x, double y, double px, double py) {

541

Matrix<double> mat = getKineticMatrix();

542

double phi = std::acos((mat(1,1) + mat(3,3))/2.);

543

if (mat(1,3) < 0.0) phi = 2.*M_PI - phi;

544

double alpha = (mat(1,1) - mat(3,3))/(2.*std::sin(phi));

545

double beta = mat(1,3)/std::sin(phi);

546

double gamma = -mat(3,1)/std::sin(phi);

547

if (std::abs(beta*gamma - alpha*alpha - 1.) > 1.e-9) {

548

std::cerr << "beta*gamma - alpha*alpha = " << beta*gamma - alpha*alpha << std::endl;

549

throw std::runtime_error("Symplectic condition beta*gamma = 1 + alpha^2 not satisfied - was there an overal Larmor rotation (method only currently implemented for cells which have a null integral of Bfield)");

527

550

}

528

Vector<double> coords(2,0.0);

529

coords(1) = x;

530

coords(2) = px;

531

Vector<double> ellipseCoords = realOrth * coords;

532

Vector<double> circleCoords(2,0.0);

533

circleCoords(1) = ellipseCoords(1) * std::sqrt(realVal(1));

534

circleCoords(2) = ellipseCoords(2) * std::sqrt(realVal(2));

535

std::pair<double, double> ret;

536

ret.first = circleCoords(1)*circleCoords(1) + circleCoords(2)*circleCoords(2);

537

ret.second = std::atan(circleCoords(2)/circleCoords(1));

538

return ret;

551

double emx = gamma*x*x + 2.*alpha*x*px + beta*px*px;

552

double emy = gamma*y*y + 2.*alpha*y*py + beta*py*py;

553

double em = std::sqrt(emx + emy);

554

//std::cerr << "alpha = " << alpha << " beta = " << beta << " gamma = " << gamma << " emittance = " << em << std::endl;

555

return em;

556

}

557

558

double TransferMatrixGenerator::getEmittance(std::vector<double> v) {

559

return this->getEmittance(v[0],v[1],v[2],v[3]);

560

}

561

562

double TransferMatrixGenerator::getPhaseAdvance() {

563

Matrix<double> mat = getKineticMatrix();

564

double phi = std::acos((mat(1,1) + mat(3,3))/2.);

565

return (mat(1,3) > 0.) ? phi : (2.*M_PI - phi);

539

566

}

540

567

541

568

void saveEllipseDataToFile(TransferMatrixGenerator* map, double zIn, double zOut, std::vector<double> initialScales, const char* filename) {

545

572

if (!stream.is_open()) {

546

573

throw InvalidFileException("LieMap::saveStabilityDataToFile",filename);

547

574

}

548

static const int timesRoundEllipse = 70;

575

static const int timesRoundEllipse = 400;

549

576

550

577

551

578

for (int i = 0; i < initialScales.size(); i++) {

553

580

double* coords = new double[4];

554

581

coords[0] = initialScales[i];

555

582

coords[1] = 0.0;

556

coords[2] = initialScales[i];

557

coords[3] = 0.0;

583

coords[2] = initialScales[i]/1.414;

584

coords[3] = -initialScales[i]/1.414;

558

585

double* coordsOut = new double[4];

559

586

for (int j = 0; j <= timesRoundEllipse; j++) {

560

map->advanceCell(coords, coordsOut,zIn,zOut);

561

562

if (coordsOut[0] > 1. || coordsOut[2] > 1.) break;

563

std::pair<double, double> actAng = map->getActionAngleCoordinates(coordsOut[0],coordsOut[2]);

564

stream << coordsOut[0]*std::cos(coordsOut[2]) << "\t" << coordsOut[0]*std::sin(coordsOut[2]) << "\n";

587

map->advanceCell(coords, coordsOut);

588

589

std::vector<double> kinCoords = BASEparticle::convertToKineticCoords(std::vector<double>(coordsOut,coordsOut+4),zOut,refMomentum,map->getIntegrator());

590

double emittance = map->getEmittance(kinCoords);

591

if (std::isnan(emittance)) break;

592

double phase = atan2(kinCoords[2], kinCoords[0]);

593

if (j == 0) std::cerr << "Phase advance = " << phase << std::endl;

594

stream << emittance*std::cos(phase) << "\t" << emittance*std::sin(phase) << "\n";

565

595

566

596

delete[] coords;

567

597

coords = coordsOut;