Skip to content

Commit 42bdd11

Browse files
osbornjdosbornjd
authored
Merge pull request #4288 from tkhemmick/PHGarfield
Initial version. Requires later updates to remove reference to user …
2 parents 82cee04 + 3f70256 commit 42bdd11

8 files changed

Lines changed: 846 additions & 0 deletions

File tree

‎offline/packages/PHGarfield/GasModel.cc‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
#include <cmath>
2+
#include <iostream>
3+
4+
#include <Garfield/ComponentUser.hh>
5+
#include <Garfield/DriftLineRKF.hh>
6+
#include <Garfield/MediumMagboltz.hh>
7+
#include <Garfield/Sensor.hh>
8+
9+
//------------------------------------------------------------
10+
// This standalone executable makes gas calculations for
11+
// whatever mixture you specify and range of electric,
12+
// magnetic, and angle between values that you select.
13+
// TKH 5/27/2026
14+
//
15+
//
16+
//------------------------------------------------------------
17+
18+
int main(int argc, char* argv[])
19+
{
20+
if (argc != 11)
21+
{
22+
std::cerr
23+
<< "Usage:\n"
24+
<< argv[0]
25+
<< " Emin Emax nE Bmin Bmax nB Amin Amax nA output_file_name\n\n"
26+
<< "Units:\n"
27+
<< " E: V/cm\n"
28+
<< " B: Tesla\n"
29+
<< " angle: radians\n";
30+
return 1;
31+
}
32+
33+
const double Emin = std::atof(argv[1]);
34+
const double Emax = std::atof(argv[2]);
35+
const int nE = std::atoi(argv[3]);
36+
37+
const double Bmin = std::atof(argv[4]);
38+
const double Bmax = std::atof(argv[5]);
39+
const int nB = std::atoi(argv[6]);
40+
41+
const double Amin = std::atof(argv[7]);
42+
const double Amax = std::atof(argv[8]);
43+
const int nA = std::atoi(argv[9]);
44+
45+
const std::string output_file(argv[10]);
46+
47+
std::cout << "E grid: "
48+
<< Emin << " -> " << Emax
49+
<< " with " << nE << " points\n";
50+
51+
std::cout << "B grid: "
52+
<< Bmin << " -> " << Bmax
53+
<< " with " << nB << " points\n";
54+
55+
std::cout << "Angle grid: "
56+
<< Amin << " -> " << Amax
57+
<< " with " << nA << " points\n";
58+
59+
std::cout << "Output File: " << output_file << std::endl;
60+
61+
// ------------------------------------------------------------
62+
// Gas: Ar/CF4/isobutane = 75/20/5.
63+
// ------------------------------------------------------------
64+
Garfield::MediumMagboltz gas;
65+
gas.SetComposition("ar", 75., "cf4", 20., "isobutane", 5.);
66+
gas.SetTemperature(301.65); // K from Grafana
67+
gas.SetPressure(762.); // Torr from Grafana
68+
69+
// Try to load an existing gas table.
70+
bool LogGrid = false;
71+
gas.SetFieldGrid(Emin, Emax, nE, LogGrid,
72+
Bmin, Bmax, nB,
73+
Amin, Amax, nA);
74+
75+
gas.GenerateGasTable(10);
76+
gas.WriteGasFile(output_file);
77+
}

‎offline/packages/PHGarfield/Makefile.am‎

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
AUTOMAKE_OPTIONS = foreign
2+
3+
ROOT_CFLAGS = $(shell root-config --cflags)
4+
ROOT_LIBS = $(shell root-config --libs)
5+
6+
AM_CPPFLAGS = \
7+
-I$(includedir) \
8+
-isystem$(OFFLINE_MAIN)/include \
9+
-isystem$(ROOTSYS)/include
10+
11+
AM_LDFLAGS = \
12+
-L$(libdir) \
13+
-L$(OFFLINE_MAIN)/lib \
14+
-L$(OFFLINE_MAIN)/lib64
15+
16+
pkginclude_HEADERS = \
17+
PHGarfield.h
18+
19+
lib_LTLIBRARIES = \
20+
libPHGarfield.la
21+
22+
libPHGarfield_la_LIBADD = \
23+
-lffamodules \
24+
-lffarawobjects \
25+
-lcdbobjects \
26+
-ltpc \
27+
-lNoRootEvent \
28+
-lphool \
29+
-lphfield \
30+
-lGarfield \
31+
-lSubsysReco
32+
33+
ROOTDICTS = \
34+
PHGarfield_Dict.cc
35+
36+
pcmdir = $(libdir)
37+
nobase_dist_pcm_DATA = \
38+
PHGarfield_Dict_rdict.pcm
39+
40+
libPHGarfield_la_SOURCES = \
41+
$(ROOTDICTS) \
42+
PHGarfield.cc
43+
44+
bin_PROGRAMS = \
45+
GasModel \
46+
MergeGasFiles
47+
48+
MergeGasFiles_SOURCES = MergeGasFiles.cc
49+
50+
MergeGasFiles_LDFLAGS = \
51+
-L$(OFFLINE_MAIN)/lib64 \
52+
$(ROOT_LIBS)
53+
54+
MergeGasFiles_LDADD = \
55+
-lGarfield
56+
57+
GasModel_SOURCES = GasModel.cc
58+
59+
GasModel_LDFLAGS = \
60+
-L$(OFFLINE_MAIN)/lib64 \
61+
$(ROOT_LIBS)
62+
63+
GasModel_LDADD = \
64+
-lGarfield
65+
66+
# Rule for generating table CINT dictionaries.
67+
%_Dict.cc: %.h %LinkDef.h
68+
rootcint -f $@ $(DEFAULT_INCLUDES) $(AM_CPPFLAGS) $^
69+
70+
#just to get the dependency
71+
%_Dict_rdict.pcm: %_Dict.cc ;
72+
73+
################################################
74+
BUILT_SOURCES = testexternals.cc
75+
76+
noinst_PROGRAMS = \
77+
testexternals
78+
79+
testexternals_SOURCES = testexternals.cc
80+
testexternals_LDADD = libPHGarfield.la
81+
82+
testexternals.cc:
83+
echo "//*** this is a generated file. Do not commit, do not edit" > $@
84+
echo "int main()" >> $@
85+
echo "{" >> $@
86+
echo " return 0;" >> $@
87+
echo "}" >> $@
88+
89+
clean-local:
90+
rm -f $(BUILT_SOURCES)

‎offline/packages/PHGarfield/MergeGasFiles.cc‎

Lines changed: 220 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,220 @@
1+
#include <TMath.h>
2+
#include <TRandom3.h>
3+
#include <ctime>
4+
#include <filesystem>
5+
#include <iostream>
6+
#include <string>
7+
8+
#include <TFile.h>
9+
#include <TNtuple.h>
10+
11+
#include <Garfield/MediumMagboltz.hh>
12+
13+
14+
int main()
15+
{
16+
// This is a utility to test whether the process of "merging" files is actually different from a single file.
17+
// It may be of no further use afterthe development was complete.
18+
// TKH 6/2/2026
19+
int nValid = 10000;
20+
TNtuple *Validity = new TNtuple("Validity", "Validity", "Valid:e:b:a:Vxerr:Vyerr:Vzerr");
21+
22+
// New version chooses to NOT write output to a file (which seems broken),
23+
// but to instead just tries to merge the files and validate the copy in memory.
24+
const std::string dir = "gasfiles";
25+
26+
auto filename = [&](const int i)
27+
{
28+
return dir + "/PART_" + std::to_string(i) + ".gas";
29+
};
30+
31+
Garfield::MediumMagboltz gas;
32+
Garfield::MediumMagboltz gas0;
33+
34+
const std::string first = filename(0);
35+
if (!std::filesystem::exists(first))
36+
{
37+
std::cerr << "Missing first gas file: " << first << std::endl;
38+
return 1;
39+
}
40+
41+
if (!gas.LoadGasFile(first))
42+
{
43+
std::cerr << "Failed to load " << first << std::endl;
44+
return 1;
45+
}
46+
47+
// Gas 0 only loads the FIRST file. This will test the memory validity of the merge...
48+
if (!gas0.LoadGasFile(first))
49+
{
50+
std::cerr << "Failed to load " << first << std::endl;
51+
return 1;
52+
}
53+
54+
for (int i = 1;; ++i)
55+
{
56+
const std::string file = filename(i);
57+
58+
if (!std::filesystem::exists(file))
59+
{
60+
std::cout << "Stopping at first missing file: " << file << std::endl;
61+
break;
62+
}
63+
64+
std::cout << "Merging " << file << std::endl;
65+
66+
if (!gas.MergeGasFile(file, true))
67+
{
68+
std::cerr << "Failed to merge " << file << std::endl;
69+
return 1;
70+
}
71+
}
72+
73+
// Don't write out since it crashes?
74+
// gas.WriteGasFile("test.gas");
75+
76+
// Now perform the validation test...
77+
double emin = 400;
78+
double emax = 400;
79+
// double ne=1;
80+
81+
double bmin = 1.15;
82+
double bmax = 1.45;
83+
double nb = 50;
84+
85+
double amin = 0.0;
86+
double amax = 0.2;
87+
// double na=50;
88+
89+
// Initialize using the current system time
90+
TRandom3 Randy(time(nullptr)); // new initialization each run
91+
std::cout << std::endl
92+
<< std::endl
93+
<< "Valid Calls: " << std::endl;
94+
for (int i = 0; i < nValid; i++)
95+
{
96+
double eMag = Randy.Uniform(emin, emax);
97+
double bMag = Randy.Uniform(bmin, bmin + 0.2 * (bmax - bmin) / nb); // Comes from file0...
98+
double a = Randy.Uniform(amin, amax);
99+
double PHI = Randy.Uniform(0.0, 2.0 * TMath::Pi());
100+
101+
double ex = 0;
102+
double ey = 0;
103+
double ez = eMag;
104+
105+
double bx = bMag * sin(a) * cos(PHI);
106+
double by = bMag * sin(a) * sin(PHI);
107+
double bz = bMag * cos(a);
108+
109+
double vx;
110+
double vy;
111+
double vz;
112+
113+
double vx0;
114+
double vy0;
115+
double vz0;
116+
117+
gas.ElectronVelocity(ex, ey, ez, bx, by, bz, vx, vy, vz);
118+
gas0.ElectronVelocity(ex, ey, ez, bx, by, bz, vx0, vy0, vz0);
119+
120+
/*
121+
std::cout << " i:" <<i;
122+
std::cout << " Ex:"<<ex;
123+
std::cout << " Ey:"<<ey;
124+
std::cout << " Ez:"<<ez;
125+
std::cout << " Bx:"<<bx;
126+
std::cout << " By:"<<by;
127+
std::cout << " Bz:"<<bz;
128+
std::cout << " Vx:"<<vx;
129+
std::cout << " Vy:"<<vy;
130+
std::cout << " Vz:"<<vz;
131+
std::cout << std::endl;
132+
std::cout << " i:" <<i;
133+
std::cout << " Ex:"<<ex;
134+
std::cout << " Ey:"<<ey;
135+
std::cout << " Ez:"<<ez;
136+
std::cout << " Bx:"<<bx;
137+
std::cout << " By:"<<by;
138+
std::cout << " Bz:"<<bz;
139+
std::cout << " Vx0:"<<vx0;
140+
std::cout << " Vy0:"<<vy0;
141+
std::cout << " Vz0:"<<vz0;
142+
std::cout << " Vx%:"<<(vx-vx0)/vx;
143+
std::cout << " Vy%:"<<(vy-vy0)/vy;
144+
std::cout << " Vy%:"<<(vz-vz0)/vz;
145+
std::cout << std::endl;
146+
std::cout << std::endl;
147+
*/
148+
149+
double DelVx = (vx - vx0) / ((vx + vx0) / 2.);
150+
double DelVy = (vy - vy0) / ((vy + vy0) / 2.);
151+
double DelVz = (vz - vz0) / ((vz + vz0) / 2.);
152+
Validity->Fill(1, sqrt(ex * ex + ey * ey + ez * ez), sqrt(bx * bx + by * by + bz * bz), a, DelVx, DelVy, DelVz);
153+
}
154+
155+
std::cout << std::endl
156+
<< std::endl
157+
<< "Invalid Calls: " << std::endl;
158+
for (int i = 0; i < nValid; i++)
159+
{
160+
double eMag = Randy.Uniform(emin, emax);
161+
double bMag = Randy.Uniform(bmin + 5.0 * (bmax - bmin) / nb, bmax); // Comes from beyond file0...
162+
double a = Randy.Uniform(amin, amax);
163+
double PHI = Randy.Uniform(0.0, 2.0 * TMath::Pi());
164+
165+
double ex = 0;
166+
double ey = 0;
167+
double ez = eMag;
168+
169+
double bx = bMag * sin(a) * cos(PHI);
170+
double by = bMag * sin(a) * sin(PHI);
171+
double bz = bMag * cos(a);
172+
173+
double vx;
174+
double vy;
175+
double vz;
176+
double vx0;
177+
double vy0;
178+
double vz0;
179+
gas.ElectronVelocity(ex, ey, ez, bx, by, bz, vx, vy, vz);
180+
gas0.ElectronVelocity(ex, ey, ez, bx, by, bz, vx0, vy0, vz0);
181+
/*
182+
std::cout << " i:" <<i;
183+
std::cout << " Ex:"<<ex;
184+
std::cout << " Ey:"<<ey;
185+
std::cout << " Ez:"<<ez;
186+
std::cout << " Bx:"<<bx;
187+
std::cout << " By:"<<by;
188+
std::cout << " Bz:"<<bz;
189+
std::cout << " Vx:"<<vx;
190+
std::cout << " Vy:"<<vy;
191+
std::cout << " Vz:"<<vz;
192+
std::cout << std::endl;
193+
std::cout << " i:" <<i;
194+
std::cout << " Ex:"<<ex;
195+
std::cout << " Ey:"<<ey;
196+
std::cout << " Ez:"<<ez;
197+
std::cout << " Bx:"<<bx;
198+
std::cout << " By:"<<by;
199+
std::cout << " Bz:"<<bz;
200+
std::cout << " Vx0:"<<vx0;
201+
std::cout << " Vy0:"<<vy0;
202+
std::cout << " Vz0:"<<vz0;
203+
std::cout << " Vx%:"<<(vx-vx0)/vx;
204+
std::cout << " Vy%:"<<(vy-vy0)/vy;
205+
std::cout << " Vy%:"<<(vz-vz0)/vz;
206+
std::cout << std::endl;
207+
std::cout << std::endl;
208+
*/
209+
double DelVx = (vx - vx0) / ((vx + vx0) / 2.);
210+
double DelVy = (vy - vy0) / ((vy + vy0) / 2.);
211+
double DelVz = (vz - vz0) / ((vz + vz0) / 2.);
212+
Validity->Fill(0, sqrt(ex * ex + ey * ey + ez * ez), sqrt(bx * bx + by * by + bz * bz), a, DelVx, DelVy, DelVz);
213+
}
214+
215+
TFile *output = new TFile("GarfieldValidity.root", "RECREATE");
216+
Validity->Write();
217+
output->Close();
218+
219+
return 0;
220+
}

0 commit comments

Comments
 (0)