Some Tutorial Examples

Note that the provided tutorials focus on simulating with a circular accelerator - linac tutorials will be added in future.

A basic knowledge of C++ is assumed, so includes, namespaces, defines and typedefs will not be discussed.

The tutorials are designed to done methodically in order, i.e. only new code will be discussed in each subsequent tutorial section.

Please don't take these examples as templates to be modified to match your own requirements. That's not their purpose. They're written to make things clear to a novice, which is not necessarily the 'best' code. If you understand the tutorials you should be in shape to write your own program rather than adapting someone else's

• Tutorial 1 - LatticeConfig
• Tutorial 2 - LatticeConfigMAD
• Tutorial 3 - LatticeManipulation
• Tutorial 4 - ParticleTracking
• Tutorial 5 - ParticleBunchTracking
• Tutorial 6 - TheLHClattice
• Tutorial 7 - CollimationAndScattering i

Tutorial 1 - LatticeConfig

The following snippet shows how to manually define a simple FODO lattice storage ring, made solely of quadrupole and dipole,. within Merlin++.

First, an instance of the AcceleratorModelConstructor class is opened, and a AcceleratorModelConstructor::NewModel() is defined. Subsequently, a loop appends individual components each with defined length, field and distance from the previous component. A final drift is appended at the end of the loop to complete the circle. Finally, an instance of AcceleratorModel is defined by the member function AcceleratorModelConstructor::GetModel().

Tutorial 2 - LatticeConfigMAD

The provided LatticeConfigMAD tutorial demonstrates how to build an accelerator model by importing a MAD .tfs file with the Merlin++ API.

The script constructs a stable FODO lattice storage consisting of dipoles, quadrupoles and sextupole corrector magnets. Furthermore, this tutorial goes on to show how to calculate and plot commonly sought after lattice functions.

• The following shows how to import a lattice directly from a MAD output .tfs file. The a MADInterface class imports the .tfs using C++ ifstream and the information is stored in a DataTable. The model is then constructed using the MADInterface::ConstructModel() function.

  
  // Instantiate MADInterface with lattice file input
  MADInterface MADinput(lattice_path, beamenergy);
  // Construct Model
  AcceleratorModel* theModel = MADinput.ConstructModel();
  

  ---

• Lattice functions are calculated and stored using the LatticeFunctionTable class. They can then be printed to a file, if desired.

  ---

  // Calculate beta and dispersion functions
  LatticeFunctionTable latticeFunctions(theModel,beamenergy);
  latticeFunctions.Calculate();
  // Write lattice functions to output file
  ofstream latticeFunctionLog("build/tutorial2.out");
  latticeFunctions.PrintTable(latticeFunctionLog);
  

  ---

• A corresponding python script is provided to plot the lattice function output. The result should be as shown here:

Tutorial 3 - LatticeManipulation

The LatticeManipulation tutorial shows how to alter individual component parameters within a user script directly to simulation alignment/field errors etc. This is done via the MagnetMover and MutipoleField classes.

• The following snippet shows how to manually realign magnet components using the MagnetMover class. Two examples are shown. The first offsets the 20th magnet vertically (i.e. along the y-axis) by 100 μm. The second offsets the 40th magnet horizontally (i.e. along the x-axis) by 50 μm.

  ---

  // Extract a list of magnet movers from the AcceleratorModel
  MagnetMoverList magnetMovers;
  theModel->ExtractTypedElements(magnetMovers);
  // Offset 20th magnet 100 um on y-axis
  magnetMovers[20]->SetY(100.0e-6);
  // Offset 40th magnet 50 um on x-axis
  magnetMovers[40]->SetX(50.0e-6);
  

  ---

• The following snippet shows how to manually alter component fields. Two examples are shown. The first increases the 5th quadrupole field by 20%. The second increases the 20th quadrupole by 10%.

  ---

  // Extract a list of the quadrupoles from the AcceleratorModel
  vector quadVec;
  theModel->ExtractTypedElements(quadVec);
  // Simulate a 20% gradient error on the 5th quadrupoles
  MultipoleField& field = quadVec[5]->GetField();
  Complex b1a = field.GetComponent(1);
  field.SetComponent(1, b1a.real() * 1.20, b1a.imag() * 1.20);
  // Simulate a 10% gradient error on the 20th quadrupoles
  MultipoleField& field2 = quadVec[20]->GetField();
  Complex b1b = field2.GetComponent(1);
  field2.SetComponent(1, b1b.real() * 1.10, b1b.imag() * 1.10);
  

  ---

• • A corresponding python script is provided to plot the changes in lattice functions as result of the manipulations. The result should be as follows: 20

Tutorial 4 - ParticleTracking

The ParticleTracking tutorial shows how to create and track a ‘bunch’ of 2 particles at different starting real space coordinates using the Particle, ParticleBunch and ParticleTracker classes.

• First we use a loop to construct a particle bunch, f just 2 particles in this case) adding each defined particle individually.

  ---

  // Define particle bunch and add individual particles in a loop
  Particle p(0);
  ParticleBunch* theBunch = new ParticleBunch(beamenergy);
  for(int xi = 1; xi <= 2; xi++)
  {
    p.x() = xi * 0.001; // Each particle offset by 1mm on x-axis
    theParticles->AddParticle(p);
  }
  

  ---

• Then we constuct a ParticleTracker using the accelerator and particle bunch constructed earlier..

  ---

  // Construct a ParticleTracker to perform the tracking
  ParticleTracker tracker(theModel->GetBeamline(), theParticles);

  ---

• Then we use a loop to track the constructed particles around the ring for 100 turns. The phase-space coordinates of each particle are also recorded after every turn using the ParticleTracker::GetTrackedBunch() and ParticleBunch::Output() member functions.

  ---

  // Track and record phase-space coords after each turn for 100 turns
  ofstream trackingLog("build/tutorial4.out");
  for(int turn = 0; turn < 100; turn++)
  {
    cout << "Tracking... turn: " << turn+1 << endl;
    if(turn == 0)
      tracker.Run();
    else
      tracker.Continue();
    ParticleBunch& tracked = tracker.GetTrackedBunch();
    tracked.Output(trackingLog);
  }
  

  ---

• A python script is provided which plots the particle evolution recorded in the log file, in phase-space coordinates. The result should be

  

Tutorial 5 - ParticleBunchTracking

The ParticleTracking tutorial shows how to create and track a normally distributed particle bunch of 10,000 particles using the ParticleDistributionGenerator, rather than doing it 'by hand', as in the previous tutorial The tutorial continues to show a different method of recording the results, using beam position monitor (BPM) component buffers. This method records the bunch centroid in real space coordinates at every BPM component present in the constructed lattice.

• The following defines the number of particles to be generated, as well as the initial beam parameters. which are spcified as BeamData.

  ---

  // Define number of particles in bunch
  size_t npart = (size_t) 1e4;
  // Define basic beam parameters (example)
  BeamData beam;
  beam.charge = 1e8;
  beam.beta_x = 0.5500005011 * meter;
  beam.beta_y = 0.5499999849 * meter;
  beam.alpha_x = -7.115569055e-7 * meter;
  beam.alpha_y = 1.797781918e-7 * meter;
  beam.emit_x = 5.026457122e-10 * meter;
  beam.emit_y = 5.026457122e-10 * meter;
  beam.sig_z = 75.5 * millimeter;
  beam.sig_dp = 0.000113;
  beam.p0 = beamenergy;
  

  ---

• The following shows to generator a normally distributed bunch with the above beam parameters.

  ---

  // Define distribution type
  ParticleDistributionGenerator* bunchDist = new NormalParticleDistributionGenerator();
  // Generate corresponding bunch
  ParticleBunch* particleBunch = new ParticleBunch(npart, *bunchDist, beam);
  

  ---

• The following shows how to initialise the BPM buffers and track real space coordinate evolution of the bunch centroid.

  ---

  // Construct a BPMBuffer to record the bunch centroid at each BPM
  BPMVectorBuffer* bpmVecBuffer = new BPMVectorBuffer();
  BPM::SetDefaultBuffer(bpmVecBuffer);
  // While tracking: Write turn tracking results to a file
  ofstream bpmLog("build/tutorial5.out");
  vector& theBPMBuffer = bpmVecBuffer->BPMReading;
  for(vector::iterator bpm_iter = theBPMBuffer.begin(); bpm_iter != theBPMBuffer.end(); bpm_iter++)
  {
    bpmLog << std::setw(14) << (bpm_iter->x).value;
    bpmLog << std::setw(14) << (bpm_iter->x).error;
    bpmLog << std::setw(14) << (bpm_iter->y).value;
    bpmLog << std::setw(14) << (bpm_iter->y).error;
    bpmLog << endl;
  }

  ---

• A corresponding python script is provided which plots the recorded bunch centroid evolution in real space coordinates at each BPM over 2 turns. The result should be as follows:

Tutorial 6 - TheLHClattice

TheLHClattice tutorial shows how to import the LHC lattice, aperture and collimator information. A simple closed orbit calculation is performed to confirm stability and the lattice functions are subsequently calculated and recorded.

• The following snippet shows how to import and initialize the LHC lattice, aperture and collimator information.

  ---

  // Import and construct LHC lattice
  cout << "Loading MAD lattice file..." << endl;
  MADInterface MADinput("Tutorials/input/LHC.tfs", beamenergy);
  AcceleratorModel* theModel = MADinput.ConstructModel();
  // Import and define component aperture information
  cout << "Loading aperture information..." << endl;
  ApertureConfiguration* apertures = new ApertureConfiguration("Tutorials/input/LHCbeam1apertureinfo.tfs");
  apertures->ConfigureElementApertures(theModel);
  // Define material database
  cout << "Loading materials database..." << endl;
  MaterialDatabase* material_db = new MaterialDatabase();
  // Import and define collimator information
  cout << "Loading collimators database..." << endl;
  CollimatorDatabase* collimator_db = new
  CollimatorDatabase("Tutorials/input/LHCcollimatorinfo.dat",
  material_db, true);

  ---

• Note the use of a dynamic convergence loop for such a complex lattice.

  ---

  double bscale1 = 1e-22;
  while(true)
  {
    latticefunctions->ScaleBendPathLength(bscale1);
    latticefunctions->Calculate();
    if(!std::isnan(latticefunctions->Value(1, 1, 1, 0)))
    {
      break;

    }
  bscale1 *= 2;  
  }

  ---

• A corresponding python script is provided which plots the LHC lattice transverse beta and horizontal dispersion functions. The result should be as follows:

  >

Tutorial 7 - CollimationAndScattering

The CollimationAndScattering tutorial shows how to initialize a simulation such that aperture checks are carried out and particles are lost where they exceed aperture limits. Moreover, if particles are lost in a collimator component, they follow a collimation process and are scattered following a defined scatter process. The scattering process used is defined by R. Appleby et al. in [ref]. Each scattered and subsequently lost particle has its loss location record in a histogram. The following simulation uses the LHC accelerator lattice and tracks 10,000 particles for 20 turns.

• The following initializes collimator settings and aperture information, using the CollimatorDatabase and CollimatorAperture member functions.

  ---

  // Initialize collimator database
  collimator_db->MatchBeamEnvelope(false);
  collimator_db->EnableJawAlignmentErrors(false);
  collimator_db->SetJawPositionError(0.0 * nanometer);
  collimator_db->SetJawAngleError(0.0 * microradian);
  collimator_db->SelectImpactFactor(start_element, 1.0e-6);
  double impact = collimator_db->ConfigureCollimators(theModel, emittance, emittance, latticefunctions);
  // Initialize collimator aperture info
  vector TCP;
  int siz = theModel->ExtractTypedElements(TCP, start_element);
  Aperture *ap = (TCP[0])->GetAperture();
  CollimatorAperture* CollimatorJaw = dynamic_cast(ap);
  double h_offset = latticefunctions->Value(1, 0, 0, start_element_number);
  double JawPosition = CollimatorJaw->GetFullWidth() / 2.0;

  ---

• For this simulation we show the use of a HorizontalHalo2ParticleDistributioan, which is also filtered in accordance with the collimator jaw parameters. It is important to set the particle charge and enable scattering physics.

  ---

  // Define horizontal bunch filter
  HorizontalHaloParticleBunchFilter* hFilter = new HorizontalHaloParticleBunchFilter();
  hFilter->SetHorizontalLimit(JawPosition);
  hFilter->SetHorizontalOrbit(h_offset);
  // Construct corresponding bunch
  ProtonBunch* particleBunch = new ProtonBunch(npart, HorizontalHalo2ParticleDistributionGenerator(), beamData, hFilter);
  particleBunch->SetMacroParticleCharge(beamData.charge);

  ---

• The following shows how to enable collimation and scattering. Particle loss locations are stored and loss maps can be produced using the LossMapCollimationOutput class.

  ---

  // Enable scattering physics
  particleBunch->EnableScatteringPhysics(ProtonBunch::Merlin);
  LossMapCollimationOutput* lossOutput = new LossMapCollimationOutput(tencm);
  ScatteringModel* scatterModel = new ScatteringModelMerlin;
  CollimateProtonProcess* collimateProcess = new CollimateProtonProcess(2, 4);
  collimateProcess->SetScatteringModel(scatterModel);
  collimateProcess->ScatterAtCollimator(true);
  collimateProcess->SetLossThreshold(200.0);
  collimateProcess->SetOutputBinSize(0.1);
  collimateProcess->SetCollimationOutput(lossOutput);
  tracker->AddProcess(collimateProcess);
  

  ---

• A corresponding python script is provided which plots the LHC loss map for 10,000 over 20 turns with the above process properties. The result should be as follows: