Howdy,
So my program works pretty great now thanks to previous help here, thank you. I have discovered an issue that’s more of a logic error on my end than anything in my code but I’d appreciate another set of eyes on it. I think what I want to do is doable, I just don’t know the syntax.
Basically what I want my code to do is have the initial Kinetic Energy of Silicon-29 Ions in my system dependent on the Recoil Angle. I already have the relationships all calculated, and I already have a Si-29 gas file in place with an Energy that covers everything I need, I just always get the full energy when I try to tell it to start paying attention when the initial energy should be at a different, lower amount.
Here’s my code if it’ll help.
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <math.h>
#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/TrackSrim.hh" //Probably not needed
#include "Garfield/TrackTrim.hh"
using namespace Garfield;
int main(int argc, char* argv[]) {
TApplication app("app", &argc, argv);
srand(time(NULL));
///////////////////////////////////////////////////////////////////////////////////////////////
// Make a gas medium.
MediumMagboltz gas;
gas.LoadGasFile("isobutane.gas");
///////////////////////////////////////////////////////////////////////////////////////////////
// Make a component with analytic electric field.
ComponentAnalyticField cmp;
cmp.SetMedium(&gas); // Attach gas to cmp
// Set parameters
const int nElec = 19; // Number of electrodes - assumes CACACACA.... //(19)
const int nWires = 50; // Wires per electrode //(50)
const double dWire = 0.0007*2.54; // Diameter of each wire [cm]
const double elecspac = 0.497 * 2.54; // Spacing between electrodes [cm]
const double wirespac = 0.2; // Spacing between wires [cm]
//Only have one of each Cathode and Anode not commented out at a time
//const double vCathode = -250.; // Potential of cathodes [V]
//const double vAnode = 250.; // Potential of anodes [V]
const double vCathode = 0.; // Potential of cathodes [V]
const double vAnode = 500.; // Potential of anodes [V] //!! 500ish?
const double chamberRad = (nWires * wirespac / 2.) +2.5; // Radius of chamber [cm] //!!
double vWire;
std::string label;
// Create electrodes
bool cathode = true; // First electrode is a cathode
for (int j = 0; j < nElec; j++) {
// Set y positions
double ywire = elecspac * j;
// Set potential and label
if (cathode) {
vWire = vCathode;
label = "C"; // Elements with same label are grouped together
}
else if (!cathode) {
vWire = vAnode;
label = "A";
}
// Create wires
for (int i = -nWires / 2; i < nWires / 2; i++) {
double xwire = wirespac * i + wirespac/2.;
cmp.AddWire(xwire, ywire, dWire, vWire, label);
}
// Set 'cathode' for next electrode
cathode = !cathode; // Alternate CACACA...
}
// Define area (for sensor and outer casing)
const double xmin = -chamberRad;
const double xmax = chamberRad;
const double ymin = -elecspac; // !! Should actually be size of dead region
const double ymax = nElec * elecspac; // !!?
const double zmin = -chamberRad;
const double zmax = chamberRad;
const unsigned int nx = 1000; //Should give point density, Default is 200
const unsigned int ny = 1000; //Should give point density, Default is 200
// Create outer casing - model w 4 planes
cmp.AddPlaneX(xmin, 0., "P1");
cmp.AddPlaneY(ymin, 0., "P2");
cmp.AddPlaneX(xmax, 0., "P3");
cmp.AddPlaneY(ymax, 0., "P4");
// Request calculation of the weighting field.
cmp.AddReadout("A"); // Only anodes read out
///////////////////////////////////////////////////////////////////////////////////////////////
// Make a sensor.
Sensor sensor;
sensor.AddComponent(&cmp); // Signal from cmp
// Set area // Area the sensor considers
sensor.SetArea(xmin, ymin, zmin, xmax, ymax, zmax);
sensor.ClearSignal();
///////////////////////////////////////////////////////////////////////////////////////////////
// Set up TRIM
TrackTrim Al26tr, Si29tr;
Al26tr.SetSensor(&sensor);
Si29tr.SetSensor(&sensor);
//const std::string file = "Aluminum-26 in No. 340 Isobutane (gas).txt";
//tr.ReadFile(file);
const unsigned int nIons = 10000;
const unsigned int nSkip = 0;
if ((!Al26tr.ReadFile("Al26EXYZ.txt", nIons, nSkip)) || (!Si29tr.ReadFile("Si29EXYZ.txt", nIons, nSkip))) {
std::cerr << "Reading TRIM EXYZ files failed.\n";
return 1;
}
/*
Kinetic Energy -
Ecm = 2.5 MeV , 0.72 MeV/u -> T = 18.73e6
Ecm = 5 MeV , 1.441 MeV/u -> T = 37.46e6
Energy losses relative to a blocker of mylar either 1.12 or 0.71 um thick
0.72 MeV/u, 0.71 um -> el = 0.100 -> AlEkin = 16.857e6
0.72 MeV/u, 1.12 um -> el = 0.162 -> AlEkin = 15.696e6
1.441 MeV/u, 0.71 um -> el = 0.046 -> AlEkin = 35.737e6
1.441 MeV/u, 1.12 um -> el = 0.073 -> AlEkin = 34.725e6
ekin = (1 - el)*T
This is a crude approach, but good enough for now.
*/
double AlEkin = 16.857e6;
double SiEkin;
Al26tr.SetWorkFunction(30); // Work Function, just set to 30 eV
Si29tr.SetWorkFunction(30);
//tr.SetAtomicMassNumbers(58.12,34); // Set (A,Z) of gas of medium in SRIM
///////////////////////////////////////////////////////////////////////////////////////////////
// Set up drift lines
std::cout << "RKF integration\n";
DriftLineRKF drift;
drift.SetMaximumStepSize(0.1 * wirespac); // Setting Max Step Size; Experimental to get rid of "Hooks"
drift.SetSensor(&sensor); // Attach to the sensor
///////////////////////////////////////////////////////////////////////////////////////////////
// Set canvas
std::cout << "Set canvas for drift\n";
TCanvas* cD = nullptr; // Canvas cD will be used to plot the electron drift lines
cD = new TCanvas("cD", "", 600, 600); // Sets stuff up but won't be plotted yet
//SetNumberOfSamples2d(nx, ny); // Set number of points, default 200
// Defined outside loop to allow it to be used after loop
ViewCell cellView; // To visualise outline of cell i.e. tube, wire etc
cellView.SetComponent(&cmp); // Viewing cmp
ViewDrift driftViewAl, driftViewSi; // To visualise electron drift lines
driftViewAl.SetCanvas(cD); // Set canvas to cD
driftViewAl.SetArea(xmin, ymin, xmax, ymax); // Restrict area to be plotted
driftViewSi.SetCanvas(cD); // Set canvas to cD
driftViewSi.SetArea(xmin, ymin, xmax, ymax); // Restrict area to be plotted
drift.EnablePlotting(&driftViewAl); // Also need to specifically enable plotting drift
drift.EnablePlotting(&driftViewSi);
Al26tr.EnablePlotting(&driftViewAl); // TRIM plotting
Si29tr.EnablePlotting(&driftViewSi);
cellView.SetCanvas(cD); // cellView will be plotted on same canvas
cellView.SetArea(xmin, ymin, xmax, ymax);
///////////////////////////////////////////////////////////////////////////////////////////////
// Open output file
std::ofstream output;
output.open( "out.txt" );
output << "# xi yi wire time \n";
///////////////////////////////////////////////////////////////////////////////////////////////
float const a = 0.090, b = 0.185, c = 17.68, d = 7.02; //explanation in loop below
double theta;
// Main loop
int n1 = 100; // 10 Electrons per row
//int n1 = 100; // 100 Electrons per row
int n2 = 100;//100; // Number of rows - "10"
//int n2 = 100; // Number of rows - "100"
for (int j = 0; j < n2; ++j) { // Each row //~20mins per row of 200 electrons
for (int k = 0; k < n1; ++k) { // Each electron in the row
// Simulate a track.
//const double xt = xmin;
const double xt = ((xmin + xmax) / 2) + (100*j + k) * (xmax - xmin) / 2500000; //y-track; 2500000 goes to 1 mm width
//const double yt = ymin + (10*j + k) * (ymax - ymin) / (n1 * n2);
const double yt = ymin; //y-track
//track.NewTrack(xt, yt, 0., 0., 1., 0., 0.);
Al26tr.NewTrack(xt, yt, 0., 0., 0., 1., 0.); //y-track;
///////////////////////////////////////////////////////////////////
//Silicon-29 Ion Angular Spread and Chamber Energy
//SiEkin = c + b*theta - a*theta*theta, theta_limit = d
//(a, b, c, d) calculated using LISE++ and Catkin
//(1.441 MeV/u, 1.12 um, 0 Excitation) -> (0.214, 0.342, 34.94, 5.69)
//(1.441 MeV/u, 1.12 um, 8 Excitation) -> (0.314, 0.072, 31.08, 1.86)
//(0.72 MeV/u, 1.12 um, 0 Excitation) -> (0.089, 0.181, 16.485, 7.02)
//(1.441 MeV/u, 0.71 um, 0 Excitation) -> (0.212, 0.341, 36.06, 5.69)
//(1.441 MeV/u, 0.71 um, 8 Excitation) -> (0.311, 0.070, 32.22, 1.86)
//(0.72 MeV/u, 0.71 um, 0 Excitation) -> (0.090, 0.185, 17.68, 7.02)
//////////////////////////////////////////////////////////////////
theta = -d + (100*j + k) * d / 5000;
SiEkin = (c + b*abs(theta) - a*theta*theta) * 1000000;
Si29tr.NewTrack(xt + 16*tan( theta * M_PI / 180), yt, 0., 0., sin( theta * M_PI / 180),cos( theta * M_PI / 180), 0.); //
////////////////////////////////////////////////////////////////////////////////
//Test of Track Output
// Loop over the clusters.
double xc, yc, zc, tc, ec, extra;
int nc;
while (Al26tr.GetCluster(xc, yc, zc, tc, nc, ec, extra) && Si29tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
for (int j = 0; j < nc; ++j) {
Al26tr.GetCluster(xc, yc, zc, tc, nc, ec, AlEkin);
Si29tr.GetCluster(xc, yc, zc, tc, nc, ec, SiEkin);
}
}
}
output << "\n \n";
}
driftViewSi.SetColourTracks(kViolet);
///////////////////////////////////////////////////////////////////////////////
// Close output file
output.close();
// Plot drift
std::cout << "Plot drift\n";
cD->Clear();
cellView.Plot2d(); // Plot2d to plot cell (set up earlier)
constexpr bool twod = true;
constexpr bool drawaxis = false;
driftViewAl.Plot(twod, drawaxis); // Plot to plot drift lines (set up earlier) - NOTE: when drawaxis is true, cell gets overwritten... FIX(?)
driftViewSi.Plot(twod, drawaxis); // Plot to plot drift lines (set up earlier) - NOTE: when drawaxis is true, cell gets overwritten... FIX(?)\
///////////////////////////////////////////////////////////////////////////////////////////////
// Set up fieldView to plot the electric potential
std::cout << "Set up fieldView\n";
TCanvas* cF = nullptr; // Canvas cF will be used to plot potential/field
ViewField fieldView; // ViewField object manages this
cF = new TCanvas("cF", "", 700, 600); // Doesn't rely on the main simulation loop so we can set up then plot immediately
fieldView.SetCanvas(cF);
fieldView.SetSensor(&sensor); // ViewField gets the field/potentials from sensor (sensor is an interface between a lot of things)
fieldView.SetPlane(0., 0., 1., 0., 0., 0.); // Set plane to plot - xy plane
cellView.SetCanvas(cF); // Also plot cell on same plot (this is why cellView was set up BEFORE main loop)
cellView.SetPlane(0., 0., 1., 0., 0., 0.); // Set plane to plot - xy plane
cellView.SetArea(xmin, ymin, xmax, ymax); // Restrict area (so that it lines up with area set by sensor)
// Plot field
std::cout << "Plot field\n";
//fieldView.Plot("e","CONT2"); // ARR looks promising but too many cells
fieldView.PlotContour("v"); // Plot "e" for field magnitude, "v" for potential (see table 4.2, pg55)
cellView.Plot2d(); // Also plot cell
///////////////////////////////////////////////////////////////////////////////////////////////
app.Run(kTRUE);
}
Any help would be appreciated.
Edit:- So some experimentation and it looks like I can send in energy values all I want, the program just goes with the max value I gave the gas file. Which I did make a point to make sure those values weren’t crazy, but doesn’t help with what I’m trying to do.