Skip to content

Commit bd106d5

Browse files
nathanielerowenathanielerowe
committed
add "back-up" FOM calculations
1 parent b440c1b commit bd106d5

6 files changed

Lines changed: 102 additions & 22 deletions

File tree

‎sbncode/BeamSpillInfoRetriever/ICARUSBNBRetriever/ICARUSBNBRetriever_module.cc‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -451,8 +451,10 @@ int sbn::ICARUSBNBRetriever::matchMultiWireData(
451451
}//end loop over MWR devices
452452

453453
sbn::BNBSpillInfo spillInfo = sbn::pot::makeBNBSpillInfo(eventID, times_temps[i], MWRdata, matched_MWR, bfp, offsets, vp873);
454-
double const spillFOM = sbn::getBNBqualityFOM(spillInfo);
455-
spillInfo.FOM = spillFOM;
454+
std::tuple<float, float, float> allFOM = sbn::getBNBqualityFOM(spillInfo);
455+
spillInfo.FOM = std::get<0>(allFOM);
456+
spillInfo.PreFitFOM = std::get<1>(allFOM);
457+
spillInfo.NoMultiWireFOM = std::get<2>(allFOM);
456458

457459
beamInfos.push_back(std::move(spillInfo));
458460

‎sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -222,9 +222,10 @@ int sbn::SBNDBNBRetriever::matchMultiWireData(
222222
}//end loop over MWR devices
223223

224224
sbn::BNBSpillInfo spillInfo = makeBNBSpillInfo(eventID, times_temps[i], MWRdata, matched_MWR, bfp, offsets, vp873);
225-
226-
double const spillFOM = sbn::getBNBqualityFOM(spillInfo);
227-
spillInfo.FOM = spillFOM;
225+
std::tuple<float, float, float> allFOM = sbn::getBNBqualityFOM(spillInfo);
226+
spillInfo.FOM = std::get<0>(allFOM);
227+
spillInfo.PreFitFOM = std::get<1>(allFOM);
228+
spillInfo.NoMultiWireFOM = std::get<2>(allFOM);
228229

229230
beamInfos.push_back(std::move(spillInfo));
230231

‎sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -228,8 +228,10 @@ void sbn::SBNDBNBZEROBIASRetriever::matchMultiWireData(
228228
}//end loop over MWR devices
229229

230230
sbn::BNBSpillInfo spillInfo = makeBNBSpillInfo(eventID, times_temps[i], MWRdata, matched_MWR, bfp, offsets, vp873);
231-
double spillFOM = sbn::getBNBqualityFOM(spillInfo);
232-
spillInfo.FOM = spillFOM;
231+
std::tuple<float, float, float> allFOM = sbn::getBNBqualityFOM(spillInfo);
232+
spillInfo.FOM = std::get<0>(allFOM);
233+
spillInfo.PreFitFOM = std::get<1>(allFOM);
234+
spillInfo.NoMultiWireFOM = std::get<2>(allFOM);
233235

234236
beamInfos.push_back(std::move(spillInfo));
235237

‎sbncode/BeamSpillInfoRetriever/getFOM.cpp‎

Lines changed: 86 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -17,9 +17,11 @@ namespace sbn {
1717

1818
bool onePlot = true;
1919

20-
float getBNBqualityFOM(BNBSpillInfo & spill )
20+
std::tuple<float, float, float> getBNBqualityFOM(BNBSpillInfo & spill )
2121
{
2222
double fom=0;
23+
double prefitfom=0;
24+
double noMWfom=0;
2325

2426
double hp875_offset= spill.HP875Offset;
2527
double vp875_offset= spill.VP875Offset;
@@ -57,6 +59,16 @@ namespace sbn {
5759
std::vector<double> vptg1;
5860
std::vector<double> hptg2;
5961
std::vector<double> vptg2;
62+
63+
std::vector<double> m875hs; // Multiwire station after Mag 875, Fit to Horizontal Sigma
64+
std::vector<double> m875hm; // Multiwire station after Mag 875, Fit to Horizontal Mean
65+
std::vector<double> m875vs; // Multiwire station after Mag 875, Fit to Vertical Sigma
66+
std::vector<double> m875vm; // Multiwire station after Mag 875, Fit to Vertical Mean
67+
68+
std::vector<double> m876hs; // Multiwire station after Mag 876, Fit to Horizontal Sigma
69+
std::vector<double> m876hm; // Multiwire station after Mag 876, Fit to Horizontal Mean
70+
std::vector<double> m876vs; // Multiwire station after Mag 876, Fit to Vertical Sigma
71+
std::vector<double> m876vm; // Multiwire station after Mag 876, Fit to Vertical Mean
6072

6173
std::vector<double> mw875(spill.M875BB.begin(), spill.M875BB.end());
6274
std::vector<double> mw876(spill.M876BB.begin(), spill.M876BB.end());
@@ -73,13 +85,23 @@ namespace sbn {
7385
hptg2.push_back(spill.HPTG2);
7486
vptg2.push_back(spill.VPTG2);
7587

88+
m875hs.push_back(spill.M875HS);
89+
m875hm.push_back(spill.M875HM);
90+
m875vs.push_back(spill.M875VS);
91+
m875vm.push_back(spill.M875VM);
92+
93+
m876hs.push_back(spill.M876HS);
94+
m876hm.push_back(spill.M876HM);
95+
m876vs.push_back(spill.M876VS);
96+
m876vm.push_back(spill.M876VM);
97+
7698
double tor;
7799
if (!tor860.empty())
78100
tor=tor860[0];
79101
else if (!tor875.empty())
80102
tor=tor875[0];
81103
else
82-
return -1;
104+
return {-1,-1,-1};
83105

84106
/**
85107
* @brief when creating ntuples for pot counting script the variables are filled with -999
@@ -104,8 +126,8 @@ namespace sbn {
104126
return std::pair(ang, pos);
105127
};
106128

107-
// return 2 when missing essential beam horizontal position data:
108-
if (hp875.empty() || (hptg1.empty() && hptg2.empty())) return 2;
129+
// return 2 when missing essential beam horizontal position data:
130+
if (hp875.empty() || (hptg1.empty() && hptg2.empty())) return {2,2,2};
109131
bool const doUseHTG1 = (useHTG == 1) || hptg2.empty();
110132
auto const [ Tanhorang, horpos ] = doUseHTG1?
111133
interpolate_hp875(hptg1[0] - hptg1_offset, hptg1_zpos):
@@ -121,8 +143,8 @@ namespace sbn {
121143
return std::pair(ang, pos);
122144
};
123145

124-
// return 2 when missing essential beam horizontal position data:
125-
if (vp875.empty() || (vptg1.empty() && vptg2.empty())) return 3;
146+
// return 3 when missing essential beam horizontal position data:
147+
if (vp875.empty() || (vptg1.empty() && vptg2.empty())) return {3,3,3};
126148
bool const doUseVTG1 = (useVTG == 1) || vptg2.empty();
127149
auto const [ Tanverang, verpos ] = doUseVTG1?
128150
interpolate_vp875(vp873[0] - vp873_offset, vp873_zpos):
@@ -138,6 +160,8 @@ namespace sbn {
138160
constexpr size_t FirstXMWtgt =0;
139161
constexpr size_t FirstYMWtgt =48;
140162
const double smallSigmaX =0.5, largeSigmaX = 10, smallSigmaY = 0.3, largeSigmaY =10, maxChi2X = 20, maxChi2Y = 20;
163+
164+
// First calculate try the version of the FOM where we fit the width ourselves
141165
if (mwtgt.size()>0) {
142166
processBNBprofile(&mwtgt[FirstXMWtgt], xx, sx,chi2x);
143167
processBNBprofile(&mwtgt[FirstYMWtgt], yy, sy, chi2y);
@@ -170,12 +194,53 @@ namespace sbn {
170194
}
171195
}
172196
if (!good_tgt && !good_876 && !good_875) {
173-
//failed getting multiwire data
174-
return 4;
197+
//failed getting multiwire data
198+
fom=-999;
175199
}
176-
177-
fom=1-pow(10,sbn::calcFOM(horpos,horang,verpos,verang,tor,tgtsx,tgtsy));
178-
return fom;
200+
else{
201+
fom=1-pow(10,sbn::calcFOM(horpos,horang,verpos,verang,tor,tgtsx,tgtsy));
202+
}
203+
204+
bool good_prefit876=false;
205+
bool good_prefit875=false;
206+
sx = 0;
207+
sy = 0;
208+
209+
// Now calculate the "pre-fit" FOM values
210+
if (!m876hs.empty() && !m876vs.empty() && !m876hm.empty() && !m876vm.empty()) {
211+
sx = m876hs[0];
212+
sy = m876vs[0];
213+
double tgtsx876=p876x[0]+p876x[1]*sx+p876x[2]*sx*sx;
214+
double tgtsy876=p876y[0]+p876y[1]*sy+p876y[2]*sy*sy;
215+
if (tgtsx876>smallSigmaX && tgtsx876<largeSigmaX && tgtsy876>smallSigmaY && tgtsy876<largeSigmaY && chi2x< maxChi2X && chi2y< maxChi2Y) {
216+
tgtsx=tgtsx876;
217+
tgtsy=tgtsy876;
218+
good_prefit876=true;
219+
}
220+
}
221+
if (!good_prefit876 && !m875hs.empty() && !m875vs.empty() && !m875hm.empty() && !m875vm.empty()) {
222+
sx = m875hs[0];
223+
sy = m875vs[0];
224+
double tgtsx875=p875x[0]+p875x[1]*sx+p875x[2]*sx*sx;
225+
double tgtsy875=p875y[0]+p875y[1]*sy+p875y[2]*sy*sy;
226+
if (tgtsx875>smallSigmaX && tgtsx875<largeSigmaX && tgtsy875>smallSigmaY && tgtsy875<largeSigmaY && chi2x< maxChi2X && chi2y< maxChi2Y) {
227+
tgtsx=tgtsx875;
228+
tgtsy=tgtsy875;
229+
good_prefit875=true;
230+
}
231+
}
232+
if (!good_prefit876 && !good_prefit875) {
233+
//failed getting multiwire data
234+
prefitfom=-999;
235+
}
236+
else{
237+
prefitfom=1-pow(10,sbn::calcFOM(horpos,horang,verpos,verang,tor,tgtsx,tgtsy));
238+
}
239+
240+
// Lastly calculate the "no multiwire" fom values, defaults values for tgtsx,y will be caught
241+
// internally in order to set scalex=scaley=1
242+
noMWfom=1-pow(10,sbn::calcFOM(horpos,horang,verpos,verang,tor));
243+
return {fom, prefitfom, noMWfom};
179244
}
180245

181246

@@ -318,8 +383,16 @@ namespace sbn {
318383
sbn::swimBNB(centroid1,sigma1,
319384
cnttoups, begtoups,
320385
cx, cy, sx, sy, rho);
321-
double scalex=tgtsx/sx;
322-
double scaley=tgtsy/sy;
386+
387+
double scalex, scaley;
388+
if( tgtsx !=-999 && tgtsy!=-999){
389+
scalex=tgtsx/sx;
390+
scaley=tgtsy/sy;
391+
}
392+
else{
393+
scalex=1;
394+
scaley=1;
395+
}
323396
double fom_a=sbn::func_intbivar(cx, cy, sx*scalex, sy*scaley, rho);
324397
//swim to center of target
325398
sbn::swimBNB(centroid1,sigma1,

‎sbncode/BeamSpillInfoRetriever/getFOM.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,12 +19,12 @@ namespace sbn
1919
* The figure of merit is described in [SBN DocDB 41901](https://sbn-docdb.fnal.gov/cgi-bin/sso/ShowDocument?docid=41901).
2020
* Inputs the BNBSpillInfo and returns the BNB Quality Metric called FOM, derived from MicroBooNE's FOM
2121
*/
22-
float getBNBqualityFOM(BNBSpillInfo& spill);
22+
std::tuple<float, float, float> getBNBqualityFOM(BNBSpillInfo& spill);
2323

2424
/**
2525
* @brief Inside the getFOM script, takes the positions and angles of the beam and calculates the BNB FOM
2626
*/
27-
double calcFOM(double horpos,double horang,double verpos,double verang,double tor,double tgtsx,double tgtsy);
27+
double calcFOM(double horpos,double horang,double verpos,double verang,double tor,double tgtsx=-999,double tgtsy=-999);
2828

2929
/**
3030
* @brief Takes in the centroid and sigma of the beam, along with transfer matrices, and will determine the beam's

‎sbncode/CAFMaker/FillExposure.cxx‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,8 @@ namespace caf
2424
single_store.HPTG2 = info.HPTG2;
2525
single_store.VPTG2 = info.VPTG2;
2626
single_store.FOM = info.FOM;
27+
single_store.PreFitFOM = info.PreFitFOM;
28+
single_store.NoMultiWireFOM = info.NoMultiWireFOM;
2729
single_store.VP873 = info.VP873;
2830
single_store.VP873Offset = info.VP873Offset;
2931
single_store.VP875Offset = info.VP875Offset;

0 commit comments

Comments
 (0)