Hello,
I’m hoping that this is just a simple thing that I am just not seeing because I’ve been staring at this code all day.
Basic situation, I have code that simulates an ionization chamber with two kinds of ion beams entering from the bottom. One is for 26-Al, the other for 29-Si. Everything is working perfectly fine except for one detail.
I want 26-Al and the 29-Si tracks to be different colors. (The exact colors don’t really matter, I’ve been aiming for 26-Al being Green and 29-Si being Purple but really I just need the tracks to be different enough that colorblindness isn’t an issue). I can control the colors of the tracks just fine - except I have yet to find a way to code it that results in different colors. They always end up the same color.
Following is my code, if needed.
#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
0.72 MeV/u, 1.12 um -> el = 0.162
1.441 MeV/u, 0.71 um -> el = 0.046
1.441 MeV/u, 1.12 um -> el = 0.073
ekin = (1 - el)*T
This is a crude approach, but good enough for now.
*/
double ekin = 34.725e6;
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 driftView; // To visualise electron drift lines
driftView.SetCanvas(cD); // Set canvas to cD
driftView.SetArea(xmin, ymin, xmax, ymax); // Restrict area to be plotted
drift.EnablePlotting(&driftView); // Also need to specifically enable plotting drift
Al26tr.EnablePlotting(&driftView); // TRIM plotting
Si29tr.EnablePlotting(&driftView);
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";
///////////////////////////////////////////////////////////////////////////////////////////////
// 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+1; ++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;
Si29tr.NewTrack(xt + 1, yt, 0., 0., 0., 1., 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, ekin);
Si29tr.GetCluster(xc, yc, zc, tc, nc, ec, ekin);
}
}
}
output << "\n \n";
}
//driftView.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;
driftView.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 assistance would be appreciated.