Skip to content

Commit

Permalink
Ran format.sh (#177)
Browse files Browse the repository at this point in the history 
Co-authored-by: Greg Rischbieter <grischbieter@albany.edu>
  • Loading branch information
@grischbieter
grischbieter and Greg Rischbieter authored Aug 4, 2023
1 parent 9da152f commit 2418fef
Show file tree
Hide file tree
Showing 11 changed files with 197 additions and 184 deletions.
15 changes: 8 additions & 7 deletions GarfieldppIntegration/GenerateGarfieldGasTableForLiquidNoble.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,18 +127,19 @@ void PassHeaderInformation(std::ofstream& outFile, std::string element,
"-7----.----8----.----9----.---10----.---11----.---12----.---13--"
<< std::endl;
char tempString[200]; // Hopefully not initializing this is okay...
snprintf(tempString, 200,
"%% Created 07/01/20 at 09.04.28 < none > GAS \"none "
" \"");
snprintf(
tempString, 200,
"%% Created 07/01/20 at 09.04.28 < none > GAS \"none "
" \"");
outFile << tempString << std::endl;
outFile << " Version : 12" << std::endl;
outFile << " GASOK bits: TFTFFFFTFFFFFFFFFFFF" << std::endl;
outFile << " Identifier: " << element << " 100%, T=293.15 K, p=1.77 atm"
<< std::endl;
outFile << " Clusters :" << std::endl;
snprintf(tempString, 200,
" Dimension : F %d 1 1 0 0",
nFields);
" Dimension : F %d 1 1 0 0",
nFields);
outFile << tempString << std::endl;
outFile << " E fields " << std::endl;
outFile << " ";
Expand Down Expand Up @@ -356,8 +357,8 @@ int main(int argc, char** argv) {
// Define a file to which we print the "gas" table.
std::ofstream outFile;
char oFileName[100];
snprintf(oFileName, 100, "GasTable_%s_%dK.gas", std::get<0>(inputArgs).c_str(),
(int)std::get<5>(inputArgs));
snprintf(oFileName, 100, "GasTable_%s_%dK.gas",
std::get<0>(inputArgs).c_str(), (int)std::get<5>(inputArgs));
outFile.open(oFileName);

// Pass the header information into the file
Expand Down
4 changes: 2 additions & 2 deletions examples/LArNEST/LArNESTNeutronCapture.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,8 +51,8 @@ int main(int argc, char* argv[]) {
for (size_t v = 0; v < electric_field.size(); v++) {
// iterate over energy values
for (size_t i = 0; i < energy_vals.size(); i++) {
result = larnest.FullCalculation(NEST::LArInteraction::ER, energy_vals[i], 0,
electric_field[v], density, false);
result = larnest.FullCalculation(NEST::LArInteraction::ER, energy_vals[i],
0, electric_field[v], density, false);
output_file << "ER,";
output_file << energy_vals[i] << ",";
output_file << electric_field[v] << ",";
Expand Down
65 changes: 33 additions & 32 deletions examples/rootNEST.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,14 +36,14 @@
#define VSTEP 1e-3 // step size in keV for convolving WIMP recoil E with eff
#define SKIP 0 // number of lines of energy to skip in efficiency file (lim)
#define NUMSIGABV -0. // # of std dev to offset NR band central line, up\dn
#define UL_MIN 0.5 // minimum upper limit on # WIMP events for limit setting
#define UL_MIN 0.5 // minimum upper limit on # WIMP events for limit setting

using namespace std;
double dayNumber = 0.;

vector<vector<double> > GetBand_Gaussian(vector<vector<double> > signals);
vector<vector<double> > GetBand(vector<double> S1s, vector<double> S2s,
bool resol);
bool resol);
TRandom3 r; // r.SetSeed(10);

double band[NUMBINS_MAX][17], band2[NUMBINS_MAX][17], medians[NUMBINS_MAX];
Expand Down Expand Up @@ -284,11 +284,11 @@ int main(int argc, char** argv) {
} else {
TFeldmanCousins fc(CL);
if (numBGeventsExp != numBGeventsObs ||
numBGeventsObs == int(numBGeventsObs))
Ul = fc.CalculateUpperLimit(numBGeventsObs, numBGeventsExp);
numBGeventsObs == int(numBGeventsObs))
Ul = fc.CalculateUpperLimit(numBGeventsObs, numBGeventsExp);
else {
Ul = fc.CalculateUpperLimit(numBGeventsObs, numBGeventsExp);
//Ul = expectedUlFc(numBGeventsExp,fc); // if you want MEAN not median
Ul = fc.CalculateUpperLimit(numBGeventsObs, numBGeventsExp);
// Ul = expectedUlFc(numBGeventsExp,fc); // if you want MEAN not median
}
double powCon =
fc.CalculateUpperLimit(0., 0.); // can't do better than 2.44 ever!
Expand Down Expand Up @@ -328,7 +328,9 @@ int main(int argc, char** argv) {
for (double j = VSTEP; j < hiE;
j +=
VSTEP) { // Iterate across energies within each sample wimp mass
double eff = pow(10., 2. - aa * exp(-bb * pow(j, cc)) - dd * exp(-ee * pow(j, ff))) / 100.;
double eff = pow(10., 2. - aa * exp(-bb * pow(j, cc)) -
dd * exp(-ee * pow(j, ff))) /
100.;
// cerr << j << " " << eff << endl;
if (eff > 1. || eff < 0.) {
if (verbosity > 0)
Expand All @@ -341,12 +343,12 @@ int main(int argc, char** argv) {
xEff * NRacc; // integrating (Riemann, left sum)
// mass-dependent differential rate with
// effxacc and step size
}
}
}
// CUSTOMIZE: your upper limit here (to mimic a PLR for example)
if (Ul < UL_MIN)
Ul = UL_MIN; // maintain what is physically possible,
// mathematically/statistically
// mathematically/statistically
xSect[i] = 1e-36 * Ul /
(sigAboveThr[i] * fidMass *
time); // derive the cross-section based upon the
Expand All @@ -365,9 +367,8 @@ int main(int argc, char** argv) {
} // end spin-dep. code block
++i;
if (xSect[i - 1] < DBL_MAX && xSect[i - 1] > 0. &&
!std::isnan(
xSect[i - 1])) { // Print the results, skipping any weirdness
// (low WIMP masses prone)
!std::isnan(xSect[i - 1])) { // Print the results, skipping any
// weirdness (low WIMP masses prone)
if (verbosity > 1) cout << Ul << "\t\t";
cout << mass[i - 1] << "\t\t\t" << xSect[i - 1]
<< endl; // final answer
Expand Down Expand Up @@ -428,17 +429,15 @@ int main(int argc, char** argv) {
<< endl;
return 1;
}
if ((useS2 == 0 && std::abs(band[i][2]) < 5.) || useS2 != 0) {
error = sqrt(pow(band[i][5], 2.) + pow(band2[i][5], 2.));
chi2[0] += pow((band2[i][2] - band[i][2]) / error, 2.);
error = sqrt(pow(band[i][6], 2.) + pow(band2[i][6], 2.));
chi2[1] += pow((band2[i][3] - band[i][3]) / error, 2.);
chi2[2] += 100. * (band2[i][2] - band[i][2]) / band2[i][2];
chi2[3] += 100. * (band2[i][3] - band[i][3]) / band2[i][3];
}
else if (verbosity > 0)
cerr << "Ignoring outliers for chi^2 calculation"
<< endl;
if ((useS2 == 0 && std::abs(band[i][2]) < 5.) || useS2 != 0) {
error = sqrt(pow(band[i][5], 2.) + pow(band2[i][5], 2.));
chi2[0] += pow((band2[i][2] - band[i][2]) / error, 2.);
error = sqrt(pow(band[i][6], 2.) + pow(band2[i][6], 2.));
chi2[1] += pow((band2[i][3] - band[i][3]) / error, 2.);
chi2[2] += 100. * (band2[i][2] - band[i][2]) / band2[i][2];
chi2[3] += 100. * (band2[i][3] - band[i][3]) / band2[i][3];
} else if (verbosity > 0)
cerr << "Ignoring outliers for chi^2 calculation" << endl;
}
chi2[0] /= double(DoF - 1);
chi2[2] /= numBins;
Expand Down Expand Up @@ -503,7 +502,7 @@ int main(int argc, char** argv) {
for (i = 0; i < numBins; ++i) {
if (skewness <= 1) {
if (ERis2nd) {
band[i][2] += NUMSIGABV * band[i][3];
band[i][2] += NUMSIGABV * band[i][3];
numSigma[i] = (band2[i][2] - band[i][2]) / band2[i][3];
errorBars[i][0] =
(band2[i][2] - band2[i][5] - band[i][2] - band[i][5]) /
Expand All @@ -519,7 +518,7 @@ int main(int argc, char** argv) {
2. * owens_t((band[i][2] - band2[i][9]) / band2[i][11],
band2[i][4]);
} else {
band2[i][2] += NUMSIGABV * band2[i][3];
band2[i][2] += NUMSIGABV * band2[i][3];
numSigma[i] = (band[i][2] - band2[i][2]) / band[i][3];
errorBars[i][0] =
(band[i][2] - band[i][5] - band2[i][2] - band2[i][5]) /
Expand Down Expand Up @@ -579,7 +578,7 @@ int main(int argc, char** argv) {
if (ERis2nd) {
numSigma[i] = (band2[i][2] - band[i][2]) / band2[i][3];
NRbandX[i] = band[i][0];
band[i][2] += NUMSIGABV * band[i][3];
band[i][2] += NUMSIGABV * band[i][3];
NRbandY[i] = band[i][2];
leakage[i] =
0.5 +
Expand Down Expand Up @@ -614,7 +613,7 @@ int main(int argc, char** argv) {
} else {
numSigma[i] = (band[i][2] - band2[i][2]) / band[i][3];
NRbandX[i] = band2[i][0];
band2[i][2] += NUMSIGABV * band2[i][3];
band2[i][2] += NUMSIGABV * band2[i][3];
NRbandY[i] = band2[i][2];
leakage[i] =
0.5 +
Expand Down Expand Up @@ -653,7 +652,7 @@ int main(int argc, char** argv) {
}
if (NRbandCenter < 0 && !ERis2nd) {
for (int nb; nb < numBins; ++nb) {
NRbandY[nb] = medians[nb] + NUMSIGABV * band2[nb][3];
NRbandY[nb] = medians[nb] + NUMSIGABV * band2[nb][3];
}
}
TGraph* gr1 = new TGraph(numBins, NRbandX, NRbandY);
Expand Down Expand Up @@ -724,7 +723,8 @@ int main(int argc, char** argv) {
finalSums[1] += (double)outputs[i].size();
}
fprintf(stderr,
"\nOVERALL DISCRIMINATION (ER) or NON-ACCEPTANCE (NR) between min and maxS1 = "
"\nOVERALL DISCRIMINATION (ER) or NON-ACCEPTANCE (NR) between min "
"and maxS1 = "
"%.12f%%, total: Gaussian & non-Gaussian (tot=counting) Leakage "
"Fraction = %.12e\n",
(1. - finalSums[0] / finalSums[1]) * 100.,
Expand All @@ -743,7 +743,8 @@ int main(int argc, char** argv) {
"with corresponding leakage of %.12e\n",
(1. - LowValue) * 100., LowValue);
fprintf(stderr,
"OVERALL DISCRIMINATION between min and maxS1 = "
"OVERALL DISCRIMINATION between min "
"and maxS1 = "
"%.12f%%, Gauss, or skew-normal fits (whatever you ran) Leakage "
"Fraction = %.12e\n",
(1. - finalSums[2] / numBins) * 100., finalSums[2] / numBins);
Expand Down Expand Up @@ -824,7 +825,7 @@ void GetFile(char* fileName) {
}
}
fclose(ifp);
if (E_keV.size() < 20000 && numBins > 1 && skewness != 0) { //used to be 1e5
if (E_keV.size() < 20000 && numBins > 1 && skewness != 0) { // used to be 1e5
skewness = 0;
cerr << "WARNING: Not enough stats (at least 10^5 events) for skew fits so "
"doing Gaussian"
Expand Down Expand Up @@ -1080,7 +1081,7 @@ vector<vector<double> > GetBand(vector<double> S1s, vector<double> S2s,
if (NRbandCenter < 0 && !save && medians[j] == 0.) {
std::sort(signals[j].begin(), signals[j].end());
medians[j] = signals[j][int(floor(double(numPts) / 2. + 0.5))];
//medians[j] += NUMSIGABV * band[j][3];
// medians[j] += NUMSIGABV * band[j][3];
// band[j][2] = medians[j];
}
for (i = 0; i < (int)numPts; ++i) {
Expand Down
16 changes: 8 additions & 8 deletions include/Detectors/XENONnT_SR0.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,8 +56,8 @@ class XENONnT : public VDetector {
// Ionization and Secondary Scintillation (S2) parameters
g1_gas = 0.1533; // phd per S2 photon in gas, used to get SE size
// Modified to have SE as 31.1518/1.21
s2Fano = 1.0; // Fano-like fudge factor for SE width It will scale
// up s2 resolution. set as 1 for now.
s2Fano = 1.0; // Fano-like fudge factor for SE width It
// will scale up s2 resolution. set as 1 for now.
s2_thr = 100.; // the S2 threshold in phe or PE, *not* phd. Affects NR most
E_gas = 6.8903; // field in kV/cm between liquid/gas border and anode
// Extraction efficiency: 52.8%
Expand All @@ -82,19 +82,19 @@ class XENONnT : public VDetector {

radius = 607.3; // millimeters (fiducial rad)
// SR0 fiducial volume cut
radmax = 664.; // actual physical geo. limit TPC
// radius
radmax = 664.; // actual physical geo. limit
// TPC radius

// Electrodes positions in mm above 0 at bottom, + above
// a z-axis value of 0 means the bottom of the detector (cathode OR bottom
// PMTs) In 2-phase, TopDrift=liquid/gas border. nT dimensions are from
// https://xe1t-wiki.lngs.infn.it/doku.php?id=xenon:xenonnt:analysis:coordinate_system
TopDrift = 1550.8; // mm not cm or us (but, this *is* where dt=0)
// 5.0 mm liquid level
anode = 1553.8; // the level of the anode grid-wire plane in mm 3.0 mm
// above the interface
gate = 1545.8; // mm. This is where the E-field changes (higher)
cathode = 60.2; // mm. Defines point below which events are gamma-X
anode = 1553.8; // the level of the anode grid-wire plane in mm 3.0 mm
// above the interface
gate = 1545.8; // mm. This is where the E-field changes (higher)
cathode = 60.2; // mm. Defines point below which events are gamma-X

// 2-D (X & Y) Position Reconstruction
// Set these to zero to implement "perfect" position corrections
Expand Down
Loading

0 comments on commit 2418fef

Please sign in to comment.