hpstr
The Heavy Photon Search Toolkit for Reconstruction (hpstr) provides an interface to physics data from the HPS experiment saved in the LCIO format and converts it into an ROOT based format.
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
SimpZBiOptimizationProcessor.cxx
Go to the documentation of this file.
1#include "SimpZBiOptimizationProcessor.h"
2#include <string>
3#include <cstdlib>
4#include <iostream>
5#include <fstream>
6
7SimpZBiOptimizationProcessor::SimpZBiOptimizationProcessor(const std::string& name, Process& process)
8 : Processor(name,process) {
9 std::cout << "[SimpZBiOptimizationProcessor] Constructor()" << std::endl;
10}
11
12SimpZBiOptimizationProcessor::~SimpZBiOptimizationProcessor(){}
13
14void SimpZBiOptimizationProcessor::configure(const ParameterSet& parameters) {
15 std::cout << "[SimpZBiOptimizationProcessor] configure()" << std::endl;
16 try
17 {
18 //Basic config
19 debug_ = parameters.getInteger("debug",debug_);
20 max_iteration_ = parameters.getInteger("max_iteration", max_iteration_);
21 year_ = parameters.getInteger("year",year_);
22 cuts_cfgFile_ = parameters.getString("cuts_cfgFile",cuts_cfgFile_);
23 outFileName_ = parameters.getString("outFileName",outFileName_);
24 cutVariables_ = parameters.getVString("cutVariables", cutVariables_);
25 min_ztail_events_ = parameters.getDouble("ztail_events",min_ztail_events_);
26 scan_zcut_ = parameters.getInteger("scan_zcut",scan_zcut_);
27 step_size_ = parameters.getDouble("step_size",step_size_);
28 eq_cfgFile_ = parameters.getString("eq_cfgFile", eq_cfgFile_);
29
30 //Background
31 bkgVtxAnaFilename_ = parameters.getString("bkgVtxAnaFilename", bkgVtxAnaFilename_);
32 bkgVtxAnaTreename_ = parameters.getString("bkgVtxAnaTreename",bkgVtxAnaTreename_);
33 background_sf_ = parameters.getDouble("background_sf",background_sf_);
34
35 //MC Signal
36 variableHistCfgFilename_ =
37 parameters.getString("variableHistCfgFilename",variableHistCfgFilename_);
38 signalVtxAnaFilename_ =
39 parameters.getString("signalVtxAnaFilename", signalVtxAnaFilename_);
40 signalVtxAnaTreename_ =
41 parameters.getString("signalVtxAnaTreename",signalVtxAnaTreename_);
42 signalVtxMCSelection_ =
43 parameters.getString("signalVtxMCSelection",signalVtxMCSelection_);
44 signalMCAnaFilename_ =
45 parameters.getString("signalMCAnaFilename",signalMCAnaFilename_);
46 signal_pdgid_ =
47 parameters.getString("signal_pdgid", signal_pdgid_);
48 signal_sf_ = parameters.getDouble("signal_sf", signal_sf_);
49 signal_mass_ = parameters.getDouble("signal_mass",signal_mass_);
50 mass_window_nsigma_ = parameters.getDouble("mass_window_nsigma", mass_window_nsigma_);
51 logEps2_ = parameters.getDouble("logEps2",logEps2_);
52
53 //Expected Signal Calculation
54 radFrac_ = parameters.getDouble("radFrac", radFrac_);
55 radAcc_ = parameters.getDouble("radAcc", radAcc_);
56 bkgControlRegionFilename_ = parameters.getString("bkgControlRegionFilename", bkgControlRegionFilename_);
57 bkgControlRegionTreename_ = parameters.getString("bkgControlRegionTreename", bkgControlRegionTreename_);
58 dNdm_sf_ = parameters.getDouble("dNdm_sf", dNdm_sf_);
59
60 // New Variables //
61 new_variables_ = parameters.getVString("add_new_variables", new_variables_);
62 new_variable_params_ = parameters.getVDouble("new_variable_params", new_variable_params_);
63
64 }
65 catch (std::runtime_error& error)
66 {
67 std::cout << error.what() << std::endl;
68 }
69}
70
71//USE JSON FILE TO LOAD IN NEW VARIABLES AND VARIABLE CONFIGURATIONS
72void SimpZBiOptimizationProcessor::addNewVariables(SimpAnaTTree* MTT, std::string variable, double param){
73
74 std::cout << "[SimpZBiOptimizationProcessor]::addNewVariable " << variable << " with param " << param << std::endl;
75 if(variable == "unc_vtx_ele_zalpha")
76 MTT->addVariable_unc_vtx_ele_zalpha(param);
77 if(variable == "unc_vtx_pos_zalpha")
78 MTT->addVariable_unc_vtx_pos_zalpha(param);
79
80 if(variable == "unc_vtx_zalpha_max")
81 MTT->addVariable_unc_vtx_zalpha_max(param);
82 if(variable == "unc_vtx_zalpha_min")
83 MTT->addVariable_unc_vtx_zalpha_min(param);
84
85 if(variable == "unc_vtx_ele_iso_z0err")
86 MTT->addVariable_unc_vtx_ele_iso_z0err();
87 if(variable == "unc_vtx_pos_iso_z0err")
88 MTT->addVariable_unc_vtx_pos_iso_z0err();
89
90 if(variable == "unc_vtx_ele_z0_z0err")
91 MTT->addVariable_unc_vtx_ele_z0_z0err();
92 if(variable == "unc_vtx_pos_z0_z0err")
93 MTT->addVariable_unc_vtx_pos_z0_z0err();
94
95 if(variable == "unc_vtx_ele_isolation_cut")
96 MTT->addVariable_unc_vtx_ele_isolation_cut();
97 if(variable == "unc_vtx_pos_isolation_cut")
98 MTT->addVariable_unc_vtx_pos_isolation_cut();
99
100 if(variable == "unc_vtx_ele_z0tanlambda")
101 MTT->addVariable_unc_vtx_ele_z0tanlambda();
102 if(variable == "unc_vtx_pos_z0tanlambda")
103 MTT->addVariable_unc_vtx_pos_z0tanlambda();
104
105 if(variable == "unc_vtx_ele_z0tanlambda_right")
106 MTT->addVariable_unc_vtx_ele_z0tanlambda_right(param);
107 if(variable == "unc_vtx_pos_z0tanlambda_right")
108 MTT->addVariable_unc_vtx_pos_z0tanlambda_right(param);
109
110 if(variable == "unc_vtx_ele_z0tanlambda_left")
111 MTT->addVariable_unc_vtx_ele_z0tanlambda_left(param);
112 if(variable == "unc_vtx_pos_z0tanlambda_left")
113 MTT->addVariable_unc_vtx_pos_z0tanlambda_left(param);
114
115 if(variable == "unc_vtx_abs_delta_z0tanlambda")
116 MTT->addVariable_unc_vtx_abs_delta_z0tanlambda();
117
118 //if(variable == "unc_vtx_proj_significance")
119 // MTT->addVariable_unc_vtx_proj_significance();
120
121 //else
122 // std::cout << "[SimpZBiOptimization]::ERROR::NEW VARIABLE " << variable << " IS NOT DEFINED IN SimpAnaTTree.cxx"
123 // <<std::endl;
124}
125
126void SimpZBiOptimizationProcessor::fillEventHistograms(std::shared_ptr<ZBiHistos> histos, SimpAnaTTree* MTT){
127
128 //histos->Fill histograms for each variable defined in tree
129 std::vector<std::string> variables = MTT->getAllVariables();
130 for(std::vector<std::string>::iterator it=variables.begin(); it != variables.end(); it++) {
131 std::string var = *it;
132 histos->Fill1DHisto(var+"_h", MTT->getValue(var));
133 }
134
135 //Impact Parameter
136 histos->Fill2DHisto("z0_v_recon_z_hh",MTT->getValue("unc_vtx_z"),MTT->getValue("unc_vtx_ele_track_z0"));
137 histos->Fill2DHisto("z0_v_recon_z_hh",MTT->getValue("unc_vtx_z"),MTT->getValue("unc_vtx_pos_track_z0"));
138
139 //Inv mass
140 histos->Fill2DHisto("vtx_InvM_vtx_z_hh",MTT->getValue("unc_vtx_mass")*1000.0, MTT->getValue("unc_vtx_z"));
141 //track z0
142 histos->Fill2DHisto("ele_track_z0_v_pos_track_z0_hh", MTT->getValue("unc_vtx_ele_track_z0"),
143 MTT->getValue("unc_vtx_pos_track_z0"));
144
145 //Zalpha
146 if(MTT->variableExists("unc_vtx_ele_zalpha")){
147 histos->Fill2DHisto("z0_v_unc_vtx_zalpha_hh",MTT->getValue("unc_vtx_ele_zalpha"),
148 MTT->getValue("unc_vtx_ele_track_z0"));
149 }
150 if(MTT->variableExists("unc_vtx_pos_zalpha")){
151 histos->Fill2DHisto("z0_v_unc_vtx_zalpha_hh",MTT->getValue("unc_vtx_pos_zalpha"),
152 MTT->getValue("unc_vtx_pos_track_z0"));
153 }
154 if(MTT->variableExists("unc_vtx_zalpha_max")){
155 histos->Fill2DHisto("recon_z_v_unc_vtx_zalpha_max_hh",MTT->getValue("unc_vtx_zalpha_max"),
156 MTT->getValue("unc_vtx_z"));
157 }
158
159 //v0 projection
160 if(MTT->variableExists("unc_vtx_proj_x")){
161 histos->Fill2DHisto("unc_vtx_proj_x_v_unc_vtx_proj_y_hh",
162 MTT->getValue("unc_vtx_proj_x"), MTT->getValue("unc_vtx_proj_y"));
163 }
164 if(MTT->variableExists("unc_vtx_proj_x_sig")){
165 histos->Fill2DHisto("unc_vtx_proj_x_y_significance_hh",
166 MTT->getValue("unc_vtx_proj_x_sig"), MTT->getValue("unc_vtx_proj_y_sig"));
167 }
168 if(MTT->variableExists("unc_vtx_proj_sig")){
169 histos->Fill2DHisto("recon_z_v_proj_sig_hh",
170 MTT->getValue("unc_vtx_proj_sig"), MTT->getValue("unc_vtx_z"));
171 }
172
173 //Vertex Errors
174 if(MTT->variableExists("unc_vtx_cxx")){
175 histos->Fill2DHisto("recon_z_v_cxx_hh", MTT->getValue("unc_vtx_cxx"),MTT->getValue("unc_vtx_z"));
176 histos->Fill2DHisto("recon_z_v_cyy_hh", MTT->getValue("unc_vtx_cyy"),MTT->getValue("unc_vtx_z"));
177 histos->Fill2DHisto("recon_z_v_czz_hh", MTT->getValue("unc_vtx_czz"),MTT->getValue("unc_vtx_z"));
178 histos->Fill2DHisto("recon_z_v_czx_hh", MTT->getValue("unc_vtx_czx"),MTT->getValue("unc_vtx_z"));
179 histos->Fill2DHisto("recon_z_v_czy_hh", MTT->getValue("unc_vtx_czy"),MTT->getValue("unc_vtx_z"));
180 histos->Fill2DHisto("recon_z_v_cyx_hh", MTT->getValue("unc_vtx_cyx"),MTT->getValue("unc_vtx_z"));
181 }
182 //Z0TanLambda
183 histos->Fill2DHisto("recon_z_v_z0tanlambda_hh",
184 MTT->getValue("unc_vtx_ele_track_z0")/MTT->getValue("unc_vtx_ele_track_tanLambda"),MTT->getValue("unc_vtx_z"));
185 histos->Fill2DHisto("recon_z_v_z0tanlambda_hh",
186 MTT->getValue("unc_vtx_pos_track_z0")/MTT->getValue("unc_vtx_pos_track_tanLambda"),MTT->getValue("unc_vtx_z"));
187 if(MTT->variableExists("unc_vtx_deltaZ")){
188 histos->Fill2DHisto("recon_z_v_unc_vtx_deltaZ_hh", MTT->getValue("unc_vtx_deltaZ"),MTT->getValue("unc_vtx_z"));
189 }
190 //z0 error
191 histos->Fill2DHisto("recon_z_v_Z0err_hh",
192 MTT->getValue("unc_vtx_ele_track_z0Err"),MTT->getValue("unc_vtx_z"));
193 histos->Fill2DHisto("recon_z_v_Z0err_hh",
194 MTT->getValue("unc_vtx_pos_track_z0Err"),MTT->getValue("unc_vtx_z"));
195 //track time
196 histos->Fill2DHisto("recon_z_v_track_t_hh", MTT->getValue("unc_vtx_ele_track_t"),MTT->getValue("unc_vtx_z"));
197 histos->Fill2DHisto("recon_z_v_track_t_hh", MTT->getValue("unc_vtx_pos_track_t"),MTT->getValue("unc_vtx_z"));
198
199 //Track parameters
200 histos->Fill2DHisto("recon_z_v_ele_track_d0_hh",MTT->getValue("unc_vtx_ele_track_d0"), MTT->getValue("unc_vtx_z"));
201 histos->Fill2DHisto("recon_z_v_pos_track_d0_hh",MTT->getValue("unc_vtx_pos_track_d0"), MTT->getValue("unc_vtx_z"));
202 histos->Fill2DHisto("recon_z_v_ele_track_phi0_hh",MTT->getValue("unc_vtx_ele_track_phi0"), MTT->getValue("unc_vtx_z"));
203 histos->Fill2DHisto("recon_z_v_pos_track_phi0_hh",MTT->getValue("unc_vtx_pos_track_phi0"), MTT->getValue("unc_vtx_z"));
204 histos->Fill2DHisto("recon_z_v_ele_track_px_hh",MTT->getValue("unc_vtx_ele_track_px"), MTT->getValue("unc_vtx_z"));
205 histos->Fill2DHisto("recon_z_v_pos_track_px_hh",MTT->getValue("unc_vtx_pos_track_px"), MTT->getValue("unc_vtx_z"));
206 histos->Fill2DHisto("recon_z_v_ele_track_py_hh",MTT->getValue("unc_vtx_ele_track_py"), MTT->getValue("unc_vtx_z"));
207 histos->Fill2DHisto("recon_z_v_pos_track_py_hh",MTT->getValue("unc_vtx_pos_track_py"), MTT->getValue("unc_vtx_z"));
208 histos->Fill2DHisto("recon_z_v_ele_track_pz_hh",MTT->getValue("unc_vtx_ele_track_pz"), MTT->getValue("unc_vtx_z"));
209 histos->Fill2DHisto("recon_z_v_pos_track_pz_hh",MTT->getValue("unc_vtx_pos_track_pz"), MTT->getValue("unc_vtx_z"));
210 histos->Fill2DHisto("recon_z_v_ele_track_nhits_hh",MTT->getValue("unc_vtx_ele_track_nhits"), MTT->getValue("unc_vtx_z"));
211 histos->Fill2DHisto("recon_z_v_pos_track_nhits_hh",MTT->getValue("unc_vtx_pos_track_nhits"), MTT->getValue("unc_vtx_z"));
212
213 //track params vs params
214 histos->Fill2DHisto("ele_tanlambda_vs_phi0_hh",MTT->getValue("unc_vtx_ele_track_phi0"),
215 MTT->getValue("unc_vtx_ele_track_tanLambda"));
216 histos->Fill2DHisto("pos_tanlambda_vs_phi0_hh",MTT->getValue("unc_vtx_pos_track_phi0"),
217 MTT->getValue("unc_vtx_pos_track_tanLambda"));
218 histos->Fill2DHisto("ele_cluster_energy_v_track_p_hh",MTT->getValue("unc_vtx_ele_track_p"),
219 MTT->getValue("unc_vtx_ele_clust_E"));
220 histos->Fill2DHisto("pos_cluster_energy_v_track_p_hh",MTT->getValue("unc_vtx_pos_track_p"),
221 MTT->getValue("unc_vtx_pos_clust_E"));
222 histos->Fill2DHisto("ele_z0_vs_tanlambda_hh",MTT->getValue("unc_vtx_ele_track_tanLambda"),
223 MTT->getValue("unc_vtx_ele_track_z0"));
224 histos->Fill2DHisto("pos_z0_vs_tanlambda_hh",MTT->getValue("unc_vtx_pos_track_tanLambda"),
225 MTT->getValue("unc_vtx_pos_track_z0"));
226
227}
228
229double SimpZBiOptimizationProcessor::countControlRegionBackgroundRate(std::string inFilename, std::string tree_name,
230 double m_Ap, double Mbin, double dNdm_sf){
231 double dNdm = 0.0;
232 TFile inFile(inFilename.c_str(), "READ");
233 TTree* tree = (TTree*)inFile.Get((tree_name+"/"+tree_name+"_tree").c_str());
234 double mass;
235 std::cout << "Counting: Ap mass is " << m_Ap << std::endl;
236 tree->SetBranchAddress("unc_vtx_mass", &mass);
237 std::cout << "N entries: " << tree->GetEntries() << std::endl;
238 for(int e=0; e < tree->GetEntries(); e++){
239 tree->GetEntry(e);
240 if(mass*1000.0 > m_Ap + (Mbin/2)) continue;
241 if(mass*1000.0 < m_Ap - (Mbin/2)) continue;
242 dNdm = dNdm + 1.0;
243 }
244 dNdm = dNdm/Mbin;
245 dNdm = dNdm * dNdm_sf;
246 std::cout << "Background Rate in CR: " << dNdm << std::endl;
247 inFile.Close();
248
249 return dNdm;
250}
251
252void SimpZBiOptimizationProcessor::initialize(std::string inFilename, std::string outFilename){
253 std::cout << "[SimpZBiOptimizationProcessor] Initialize " << inFilename << std::endl;
254
255 //Load Simp Equations
256 simpEqs_ = new SimpEquations(year_, eq_cfgFile_);
257 massResolution_ = simpEqs_->massResolution(signal_mass_);
258
259 //Define Mass window
260 lowMass_ = signal_mass_ - mass_window_nsigma_*massResolution_/2.0;
261 highMass_ = signal_mass_ + mass_window_nsigma_*massResolution_/2.0;
262 std::cout << "[SimpZBiOptimization]::Mass Window: " << lowMass_ << " - " << highMass_ << std::endl;
263
264 //Initialize output file
265 std::cout << "[SimpZBiOptimization]::Output File: " << outFileName_.c_str() << std::endl;
266 outFile_ = new TFile(outFileName_.c_str(),"RECREATE");
267
268 //Get the signal pretrigger simulated vertex z distribution
269 std::cout << "[SimpZBiOptimization]::Getting MC Signal pre-trigger vertex z distribution from file "
270 << signalMCAnaFilename_ << std::endl;
271 getSignalMCAnaVtxZ_h(signalMCAnaFilename_, signal_pdgid_);
272 outFile_->cd();
273 signalSimZ_h_->Write();
274
275 //Read signal ana vertex tuple, and convert to mutable tuple
276 std::cout << "[SimpZBiOptimization]::Reading Signal AnaVertex Tuple from file "
277 << signalVtxAnaFilename_.c_str() << std::endl;
278 TFile* signalVtxAnaFile = new TFile(signalVtxAnaFilename_.c_str(),"READ");
279 signalMTT_ = new SimpAnaTTree(signalVtxAnaFile,(signalVtxAnaTreename_+"/"+signalVtxAnaTreename_+"_tree").c_str());
280 signalMTT_->defineMassWindow(lowMass_, highMass_);
281
282 //Get Simp Mean Truth Energy
283 std::cout << "Get Mean Truth Energy" << std::endl;
284 TH1F* signalEnergy_h = (TH1F*)signalVtxAnaFile->Get((signalVtxMCSelection_+"/"+signalVtxMCSelection_+"_mc"+signal_pdgid_+"Energy_h").c_str());
285 std::cout << "Got the histo " << std::endl;
286 if(signalEnergy_h == nullptr) std::cout << "HISTO IS NULL" << std::endl;
287 E_Vd_ = (double)signalEnergy_h->GetMean();
288 std::cout << "Mean Energy is " << E_Vd_ << std::endl;
289
290 //Read background ana vertex tuple, and convert to mutable tuple
291 std::cout << "[SimpZBiOptimization]::Reading Background AnaVertex Tuple from file "
292 << signalVtxAnaFilename_.c_str() << std::endl;
293 TFile* bkgVtxAnaFile = new TFile(bkgVtxAnaFilename_.c_str(),"READ");
294 bkgMTT_ = new SimpAnaTTree(bkgVtxAnaFile, (bkgVtxAnaTreename_+"/"+bkgVtxAnaTreename_+"_tree").c_str());
295 bkgMTT_->defineMassWindow(lowMass_, highMass_);
296
297 //Add new variables, as defined in SimpAnaTTree.cxx, from the processor configuration script
298 for(std::vector<std::string>::iterator it = new_variables_.begin(); it != new_variables_.end(); it++){
299 int param_idx = std::distance(new_variables_.begin(), it);
300 std::cout << "[SimpZBiOptimization]::Attempting to add new variable " << *it <<
301 " with parameter " << new_variable_params_.at(param_idx) << std::endl;
302 //signalMTT_->addVariable(*it, new_variable_params_.at(param_idx));
303 //bkgMTT_->addVariable(*it, new_variable_params_.at(param_idx));
304 addNewVariables(signalMTT_, *it, new_variable_params_.at(param_idx));
305 addNewVariables(bkgMTT_, *it, new_variable_params_.at(param_idx));
306 }
307
308 //Finalize Initialization of New Mutable Tuples
309 std::cout << "[SimpZBiOptimization]::Finalizing Initialization of New Mutable Tuples" << std::endl;
310 signalMTT_->shiftVariable("unc_vtx_ele_track_z0", -0.07);
311 signalMTT_->shiftVariable("unc_vtx_pos_track_z0", -0.07);
312 signalMTT_->Fill();
313 bkgMTT_->Fill();
314
315 //Initialize Persistent Cut Selector. These cuts are applied to all events.
316 //Persistent Cut values are updated each iteration with the value of the best performing Test Cut in
317 //that iteration.
318 std::cout << "[SimpZBiOptimization]::Initializing Set of Persistent Cuts" << std::endl;
319 persistentCutsSelector_ = new IterativeCutSelector("persistentCuts", cuts_cfgFile_);
320 persistentCutsSelector_->LoadSelection();
321 persistentCutsPtr_ = persistentCutsSelector_->getPointerToCuts();
322 std::cout << "Persistent Cuts: " << std::endl;
323 persistentCutsSelector_->printCuts();
324
325 //initalize Test Cuts
326 std::cout << "[SimpZBiOptimization]::Initializing Set of Test Cuts" << std::endl;
327 testCutsSelector_ = new IterativeCutSelector("testCuts",cuts_cfgFile_);
328 testCutsSelector_->LoadSelection();
329 testCutsPtr_ = testCutsSelector_->getPointerToCuts();
330 testCutsSelector_->filterCuts(cutVariables_);
331 std::cout << "Test Cuts: " << std::endl;
332 testCutsSelector_->printCuts();
333
334 //Initialize signal histograms
335 std::cout << "[SimpZBiOptimization]::Initializing Signal Variable Histograms" << std::endl;
336 signalHistos_= std::make_shared<ZBiHistos>("signal");
337 signalHistos_->debugMode(debug_);
338 signalHistos_->loadHistoConfig(variableHistCfgFilename_);
339 signalHistos_->DefineHistos();
340
342 std::cout << "[SimpZBiOptimization]::Initializing Background Variable Histograms" << std::endl;
343 bkgHistos_ = std::make_shared<ZBiHistos>("background");
344 bkgHistos_->debugMode(debug_);
345 bkgHistos_->loadHistoConfig(variableHistCfgFilename_);
346 bkgHistos_->DefineHistos();
347
348 //Initialize Test Cut histograms
349 std::cout << "[SimpZBiOptimization]::Initializing Test Cut Variable Histograms" << std::endl;
350 testCutHistos_ = std::make_shared<ZBiHistos>("testCutHistos");
351 testCutHistos_->debugMode(debug_);
352 testCutHistos_->DefineHistos();
353
354 //Add Test Cut Analysis Histograms necessary for calculating background and signal
355 std::cout << "[SimpZBiOptimization]::Initializing Test Cut Analysis Histograms" << std::endl;
356 for(cut_iter_ it=testCutsPtr_->begin(); it!=testCutsPtr_->end(); it++){
357 std::string name = it->first;
358 //Used to select true z vertex distribution given a cut in unc_vtx_z
359 testCutHistos_->addHisto2d("unc_vtx_z_vs_true_vtx_z_"+name+"_hh","unc z_{vtx} [mm]",
360 1500, -50.0, 100.0,"true z_{vtx} [mm]",200,-50.3,149.7);
361 //signalSelZ models signal selection efficiency
362 testCutHistos_->addHisto1d("signal_SelZ_"+name+"_h","true z_{vtx} [mm]",
363 200, -50.3, 149.7);
364 //background_zVtx provides basis for background model, used to estimate nbackground in signal region
365 testCutHistos_->addHisto1d("background_zVtx_"+name+"_h","unc z_{vtx} [mm]",
366 150, -50.0, 100.0);
367 }
368
369 //Initialize processor histograms that summarize iterative results
370 std::cout << "[SimpZBiOptimization]::Initializing Itearitve Result Histograms" << std::endl;
371 processorHistos_ = std::make_shared<ZBiHistos>("zbi_processor");
372 processorHistos_->defineZBiCutflowProcessorHistograms();
373
374 //Fill Initial Signal histograms
375 std::cout << "[SimpZBiOptimization]::Filling initial signal histograms" << std::endl;
376 for(int e=0; e < signalMTT_->GetEntries(); e++){
377 signalMTT_->GetEntry(e);
378 fillEventHistograms(signalHistos_, signalMTT_);
379 }
380
381 std::cout << "[SimpZBiOptimization]::Filling initial background histograms" << std::endl;
382 //Fill Initial Background Histograms
383 for(int e=0; e < bkgMTT_->GetEntries(); e++){
384 bkgMTT_->GetEntry(e);
385 fillEventHistograms(bkgHistos_, bkgMTT_);
386 }
387
388 //Count background rate in the Control Region (used to calculate total A' Rate)
389 double m_Ap = simpEqs_->getAprimeMassFromVectorMass(signal_mass_);
390 dNdm_ = countControlRegionBackgroundRate(bkgControlRegionFilename_, bkgControlRegionTreename_, m_Ap, 30.0, dNdm_sf_);
391
392 //Write initial variable histograms for signal and background
393 std::cout << "[SimpZBiOptimization]::Writing Initial Histograms" << std::endl;
394 signalHistos_->writeHistos(outFile_,"initial_signal");
395 bkgHistos_->writeHistos(outFile_,"initial_background");
396}
397
398bool SimpZBiOptimizationProcessor::process(){
399 std::cout << "[SimpZBiOptimizationProcessor]::process()" << std::endl;
400
401 //step_size defines n% of signal distribution to cut in a given variable
402 std::cout << "step_size_ " << step_size_ << std::endl;
403 double cutFraction = step_size_;
404 if(max_iteration_ > (int)1.0/cutFraction)
405 max_iteration_ = (int)1.0/cutFraction;
406
407 std::cout << "max iteration: " << max_iteration_ << std::endl;
408 //Iteratively cut n% of the signal distribution for a given Test Cut variable
409 for(int iteration = 0; iteration < max_iteration_+1; iteration ++){
410 double cutSignal = (double)iteration*step_size_*100.0;
411 cutSignal = round(cutSignal);
412 if(debug_) std::cout << "## ITERATION " << iteration << " ##" << std::endl;
413
414 //Reset histograms at the start of each iteration
415 //Histograms will change with each iteration.
416 if(debug_) std::cout << "[SimpZBiOptimization]::Resetting Histograms" << std::endl;
417 signalHistos_->resetHistograms1d();
418 signalHistos_->resetHistograms2d();
419 bkgHistos_->resetHistograms1d();
420 bkgHistos_->resetHistograms2d();
421 testCutHistos_->resetHistograms1d();
422 testCutHistos_->resetHistograms2d();
423
424 //At the start of each iteration, save the persistent cut values
425 for(cut_iter_ it=persistentCutsPtr_->begin(); it!=persistentCutsPtr_->end(); it++){
426 std::string cutname = it->first;
427 double cutvalue = persistentCutsSelector_->getCut(cutname);
428 if(debug_){
429 std::cout << "Saving persistent cut: " << cutname << " " << cutvalue << std::endl;
430 }
431 int cutid = persistentCutsSelector_->getCutID(cutname);
432 processorHistos_->Fill2DHisto("persistent_cuts_hh",(double)cutSignal,
433 (double)cutid,cutvalue);
434 processorHistos_->set2DHistoYlabel("persistent_cuts_hh",cutid,cutname);
435 }
436
437 //Fill signal variable distributions
438 for(int e=0; e < signalMTT_->GetEntries(); e++){
439 signalMTT_->GetEntry(e);
440
441 //Apply current set of persistent cuts to all events
442 if(failPersistentCuts(signalMTT_))
443 continue;
444
445 //Fill Signal variable distributions
446 fillEventHistograms(signalHistos_, signalMTT_);
447 }
448 if(iteration == max_iteration_){
449 //Write iteration histos
450 signalHistos_->writeHistos(outFile_,"signal_pct_sig_cut_"+std::to_string(cutSignal));
451 bkgHistos_->writeHistos(outFile_,"background_pct_sig_cut_"+std::to_string(cutSignal));
452 break;
453 }
454
455 //Initialize Signal Integrals
456 if(iteration == 0){
457 //Integrate each signal variable distribution.
458 //When iterating Test Cuts, we reference these intial integrals, cutting n% of the original distribution.
459 std::cout << "[SimpZBiOptimization]::Integrating initial Signal distributions" << std::endl;
460 for(cut_iter_ it=testCutsPtr_->begin(); it!=testCutsPtr_->end(); it++){
461 std::string cutname = it->first;
462 std::string var = testCutsSelector_->getCutVar(cutname);
463 initialIntegrals_[var] = signalHistos_->integrateHistogram1D("signal_"+var+"_h");
464 if(debug_)
465 std::cout << "Initial Integral for " << var << " is " << initialIntegrals_[var] << std::endl;
466 }
467 }
468
469 //Loop over each Test Cut. Cut n% of signal distribution in Test Cut variable
470 if(debug_) std::cout << "Looping over Signal Test Cuts" << std::endl;
471 for(cut_iter_ it=testCutsPtr_->begin(); it!=testCutsPtr_->end(); it++){
472 std::string cutname = it->first;
473 std::string cutvar = testCutsSelector_->getCutVar(cutname);
474 bool isCutGT = testCutsSelector_->isCutGreaterThan(cutname);
475 double cutvalue = signalHistos_->cutFractionOfSignalVariable(cutvar, isCutGT,
476 iteration*cutFraction, initialIntegrals_[cutvar]);
477 testCutsSelector_->setCutValue(cutname, cutvalue);
478 if(debug_){
479 std::cout << "Test Cut " << cutname << " " << testCutsSelector_->getCut(cutname) <<
480 " cuts " << cutSignal << "% of signal distribution in this variable " << std::endl;
481 }
482 }
483
484 //Fill Background Histograms corresponding to each Test Cut
485 if(debug_) std::cout << "Filling Background Variables for each Test Cut" << std::endl;
486 for(int e=0; e < bkgMTT_->GetEntries(); e++){
487 bkgMTT_->GetEntry(e);
488
489 //Apply persistent cuts
490 if(failPersistentCuts(bkgMTT_))
491 continue;
492
493 fillEventHistograms(bkgHistos_, bkgMTT_);
494
495 //Loop over each Test Cut
496 for(cut_iter_ it=testCutsPtr_->begin(); it!=testCutsPtr_->end(); it++){
497 std::string cutname = it->first;
498 std::string cutvar = testCutsSelector_->getCutVar(cutname);
499
500 //apply Test Cut
501 if(failTestCut(cutname, bkgMTT_))
502 continue;
503
504 //If event passes Test Cut, fill vertex z distribution.
505 //This distribution is used to build the Background Model corresponding to each Test Cut
506 testCutHistos_->Fill1DHisto("background_zVtx_"+cutname+"_h",
507 bkgMTT_->getValue("unc_vtx_z"),background_sf_);
508 }
509 }
510
511 //For each Test Cut, build relationship between Signal truth z_vtx, and reconstructed z_vtx
512 //This is used to get the truth Signal Selection Efficiency F(z), given a Zcut in reconstructed z_vtx
513 if(debug_) std::cout << "Build Signal truth z vs recon z" << std::endl;
514 for(int e=0; e < signalMTT_->GetEntries(); e++){
515 signalMTT_->GetEntry(e);
516
517 //Apply persistent cuts
518 if(failPersistentCuts(signalMTT_))
519 continue;
520
521 //Loop over each Test Cut and plot unc_vtx_z vs true_vtx_z
522 for(cut_iter_ it=testCutsPtr_->begin(); it!=testCutsPtr_->end(); it++){
523 std::string cutname = it->first;
524 std::string cutvar = testCutsSelector_->getCutVar(cutname);
525 double cutvalue = testCutsSelector_->getCut(cutname);
526 //Apply Test Cut
527 if(!testCutsSelector_->passCutGTorLT(cutname, signalMTT_->getValue(cutvar)))
528 continue;
529 testCutHistos_->Fill2DHisto("unc_vtx_z_vs_true_vtx_z_"+cutname+"_hh",
530 signalMTT_->getValue("unc_vtx_z"),signalMTT_->getValue("vd_true_vtx_z"),1.0);
531 }
532 }
533
534 //Calcuate the Binomial Significance ZBi corresponding to each Test Cut
535 //Test Cut with maximum ZBi after cutting n% of signal distribution in a given variable is selected
536 //Test Cut value is added to list of Persistent Cuts, and is applied to all events in following iterations.
537
538 //Create Directory for this iteration
539 std::string iteration_subdir = "testCuts_pct_sig_cut_"+std::to_string(cutSignal);
540 TDirectory* dir{nullptr};
541 std::cout<<iteration_subdir.c_str()<<std::endl;
542 if (!iteration_subdir.empty()) {
543 dir = outFile_->mkdir(iteration_subdir.c_str(),"",true);
544 }
545 double best_zbi = -9999.9;
546 double best_zcut = -9999.9;
547 double best_nsig = -9999.9;
548 double best_nbkg = -9999.9;
549 std::string best_cutname;
550 double best_cutvalue;
551
552 //Loop over Test Cuts
553 if(debug_) std::cout << "Calculate ZBi for each Test Cut " << std::endl;
554 for(cut_iter_ it=testCutsPtr_->begin(); it!=testCutsPtr_->end(); it++){
555 std::string cutname = it->first;
556 double cutvalue = testCutsSelector_->getCut(cutname);
557 int cutid = testCutsSelector_->getCutID(cutname);
558 if(debug_){
559 std::cout << "Calculating ZBi for Test Cut " << cutname << std::endl;
560 std::cout << "Test Cut ID: " << cutid << " | Test Cut Value: " << cutvalue << std::endl;
561 }
562
563 //Build Background Model, used to estimate nbkg in Signal Region
564 if(debug_) std::cout << "Build Background Model" << std::endl;
565 double start_fit = 500.0*background_sf_;
566 TF1* bkg_model = (TF1*)testCutHistos_->fitExponentialTail("background_zVtx_"+cutname, start_fit);
567 //TF1* bkg_model = (TF1*)testCutHistos_->fitExponentialPlusConst("background_zVtx_"+cutname, start_fit);
568 if(debug_) std::cout << "END Build Background Model" << std::endl;
569
570 //Get signal unc_vtx_z vs true_vtx_z
571 TH2F* vtx_z_hh = (TH2F*)testCutHistos_->get2dHisto("testCutHistos_unc_vtx_z_vs_true_vtx_z_"+cutname+"_hh");
572
573 //CD to output file to save resulting plots
574 outFile_->cd();
575
576 //Graphs track the performance of a Test Cut as a function of Zcut position
577 TGraph* zcutscan_zbi_g = new TGraph();
578 zcutscan_zbi_g->SetName(("zcut_vs_zbi_"+cutname+"_g").c_str());
579 zcutscan_zbi_g->SetTitle(("zcut_vs_zbi_"+cutname+"_g;zcut [mm];zbi").c_str());
580 zcutscan_zbi_g->SetMarkerStyle(8);
581 zcutscan_zbi_g->SetMarkerSize(2.0);
582 zcutscan_zbi_g->SetMarkerColor(2);
583
584 TGraph* zcutscan_nsig_g = new TGraph();
585 zcutscan_nsig_g->SetName(("zcut_vs_nsig_"+cutname+"_g").c_str());
586 zcutscan_nsig_g->SetTitle(("zcut_vs_nsig_"+cutname+"_g;zcut [mm];nsig").c_str());
587 zcutscan_nsig_g->SetMarkerStyle(23);
588 zcutscan_nsig_g->SetMarkerSize(2.0);
589 zcutscan_nsig_g->SetMarkerColor(57);
590
591 TGraph* zcutscan_nbkg_g = new TGraph();
592 zcutscan_nbkg_g->SetName(("zcut_vs_nbkg_"+cutname+"_g").c_str());
593 zcutscan_nbkg_g->SetTitle(("zcut_vs_nbkg_"+cutname+"_g;zcut [mm];nbkg").c_str());
594 zcutscan_nbkg_g->SetMarkerStyle(45);
595 zcutscan_nbkg_g->SetMarkerSize(2.0);
596 zcutscan_nbkg_g->SetMarkerColor(49);
597
598 TGraph* nbkg_zbi_g = new TGraph();
599 nbkg_zbi_g->SetName(("nbkg_vs_zbi_"+cutname+"_g").c_str());
600 nbkg_zbi_g->SetTitle(("nbkg_vs_zbi_"+cutname+"_g;nbkg;zbi").c_str());
601
602 TGraph* nsig_zbi_g = new TGraph();
603 nsig_zbi_g->SetName(("nsig_vs_zbi_"+cutname+"_g").c_str());
604 nsig_zbi_g->SetTitle(("nsig_vs_zbi_"+cutname+"_g;nsig;zbi").c_str());
605
606 TH2F* nsig_zcut_hh = new TH2F(("nsig_v_zcut_zbi_"+cutname+"_hh").c_str(),
607 ("nsig_v_zcut_zbi_"+cutname+"_hh; zcut [mm]; Nbkg").c_str(),
608 200,-50.3,149.7,3000,0.0,300.0);
609
610 TH2F* nbkg_zcut_hh = new TH2F(("nbkg_v_zcut_zbi_"+cutname+"_hh").c_str(),
611 ("nbkg_v_zcut_zbi_"+cutname+"_hh; zcut [mm]; Nbkg").c_str(),
612 200,-50.3,149.7,3000,0.0,300.0);
613
614
615 //Find maximum position of Zcut --> ZBi calculation requires non-zero background
616 //Start the Zcut position at the target
617 double zcut_step = 0.1;
618 TH1F* bkg_zVtx_h = (TH1F*)testCutHistos_->get1dHisto("testCutHistos_background_zVtx_"+cutname+"_h");
619 double max_zcut = bkg_model->GetXmin();
620 double endIntegral = bkg_zVtx_h->GetBinLowEdge(bkg_zVtx_h->FindLastBinAbove(0.0)) + bkg_zVtx_h->GetBinWidth(1);
621 //double endIntegral = 100.0
622 double testIntegral = bkg_model->Integral(max_zcut, endIntegral);
623 if(debug_) std::cout << "Background between " << max_zcut << "and end of histo is " << testIntegral << std::endl;
624 while(testIntegral > min_ztail_events_){
625 max_zcut = max_zcut+zcut_step;
626 testIntegral = bkg_model->Integral(max_zcut, endIntegral);
627 if(testIntegral < min_ztail_events_){
628 max_zcut = max_zcut-zcut_step;
629 testIntegral = bkg_model->Integral(max_zcut, endIntegral);
630 break;
631 }
632 }
633 if(debug_) std::cout << "Maximum Zcut: " << max_zcut << " gives " << testIntegral << " background events" << std::endl;
634
635 //If configuration does not specify scanning Zcut values, use single Zcut position at maximum position.
636 double min_zcut = bkg_model->GetXmin();
637 if(!scan_zcut_)
638 min_zcut = max_zcut;
639 std::cout << "Minimum Zcut position: " << min_zcut << std::endl;
640
641 //Get the signal vtx z selection efficiency *before* zcut is applied
642 if(debug_) std::cout << "Get signal vtx z selection efficiency before Zcut" << std::endl;
643 TH1F* true_vtx_NoZ_h = (TH1F*)vtx_z_hh->ProjectionY((cutname+"_"+"true_vtx_z_projy").c_str(),1,vtx_z_hh->GetXaxis()->GetNbins(),"");
644
645 //Convert the truth vertex z distribution beyond Zcut into the appropriately binned Selection.
646 //Binning must match Signal pre-trigger distribution, in order to take Efficiency between them.
647 TH1F* signalSelNoZ_h =
648 (TH1F*)signalSimZ_h_->Clone(("testCutHistos_signal_SelNoZ_"+cutname+"_h").c_str());
649 for(int i=0; i<201; i++){
650 signalSelNoZ_h->SetBinContent(i,true_vtx_NoZ_h->GetBinContent(i));
651 }
652 TEfficiency* effCalcNoZ_h = new TEfficiency(*signalSelNoZ_h, *signalSimZ_h_);
653
654 outFile_->cd(("testCuts_pct_sig_cut_"+std::to_string(cutSignal)).c_str());
655 signalSelNoZ_h->Write();
656 effCalcNoZ_h->Write();
657 delete effCalcNoZ_h;
658
659 //Scan Zcut position and calculate ZBi
660 double best_scan_zbi = -999.9;
661 double best_scan_zcut;
662 double best_scan_nsig;
663 double best_scan_nbkg;
664 if(debug_) std::cout << "Scanning zcut position" << std::endl;
665 for(double zcut = min_zcut; zcut < (max_zcut+zcut_step); zcut = zcut+zcut_step){
666 double Nbkg = bkg_model->Integral(zcut,endIntegral);
667 //std::cout << "zcut: " << zcut << std::endl;
668 //std::cout << "Nbkg: " << Nbkg << std::endl;
669
670 //Get the Signal truth vertex z distribution beyond the reconstructed vertex Zcut
671 TH1F* true_vtx_z_h = (TH1F*)vtx_z_hh->ProjectionY((std::to_string(zcut)+"_"+cutname+"_"+"true_vtx_z_projy").c_str(),vtx_z_hh->GetXaxis()->FindBin(zcut)+1,vtx_z_hh->GetXaxis()->GetNbins(),"");
672
673 //Convert the truth vertex z distribution beyond Zcut into the appropriately binned Selection.
674 //Binning must match Signal pre-trigger distribution, in order to take Efficiency between them.
675 TH1F* signalSelZ_h =
676 (TH1F*)signalSimZ_h_->Clone(("testCutHistos_signal_SelZ_"+cutname+"_h").c_str());
677 for(int i=0; i<201; i++){
678 signalSelZ_h->SetBinContent(i,true_vtx_z_h->GetBinContent(i));
679 }
680
681 //Get Signal Selection Efficiency, as a function of truth vertex Z, F(z)
682 //if(debug_) std::cout << "Get Signal Selection Efficiency" << std::endl;
683 TEfficiency* effCalc_h = new TEfficiency(*signalSelZ_h, *signalSimZ_h_);
684 //if(zcut == min_zcut){
685 // outFile_->cd(("testCuts_pct_sig_cut_"+std::to_string(cutSignal)).c_str());
686 // effCalc_h->Write();
687 //}
688
689 double eps2 = std::pow(10, logEps2_);
690 double eps = std::sqrt(eps2);
691
692 //Calculate expected signal for Neutral Dark Vector "rho"
693 double nSigRho = simpEqs_->expectedSignalCalculation(signal_mass_,
694 eps, true, E_Vd_, effCalc_h, dNdm_, radFrac_, radAcc_, -4.3, zcut);
695
696 //Calculate expected signal for Neutral Dark Vector "rho"
697 double nSigPhi = simpEqs_->expectedSignalCalculation(signal_mass_,
698 eps, false, E_Vd_, effCalc_h, dNdm_, radFrac_, radAcc_, -4.3, zcut);
699
700 /*
701 //Calculate expected signal for Neutral Dark Vector "rho"
702 double nSigRho = background_sf_*simpEqs_->expectedSignalCalculation(signal_mass_,
703 eps, true, false, E_Vd_, effCalc_h, -4.3, zcut);
704
705 //Calculate expected signal for Neutral Dark Vector "phi"
706 double nSigPhi = background_sf_*simpEqs_->expectedSignalCalculation(signal_mass_,
707 eps, false, true, E_Vd_, effCalc_h, -4.3, zcut);
708
709 */
710 double Nsig = nSigRho + nSigPhi;
711
712 //if(debug_){
713 // std::cout << "nSigRho: " << nSigRho << std::endl;
714 // std::cout << "nSigPhi: " << nSigPhi << std::endl;
715 // std::cout << "Nsig: " << Nsig << std::endl;
716 //}
717
718 Nsig = Nsig*signal_sf_;
719 //if(debug_) std::cout << "Nsig after scale factor: " << Nsig << std::endl;
720
721
722 //CLEAR POINTERS
723 delete effCalc_h;
724
725 //Round Nsig, Nbkg, and then ZBi later
726 Nsig = round(Nsig);
727 Nbkg = round(Nbkg);
728
729 //Calculate ZBi for this Test Cut using this zcut value
730 double n_on = Nsig + Nbkg;
731 double tau = 1.0;
732 double n_off = Nbkg;
733 double ZBi = calculateZBi(n_on, n_off, tau);
734 std::cout << "ZBi before rounding: " << ZBi << std::endl;
735 ZBi = round(ZBi);
736
737 /*
738 std::cout << "[SimpZBiOptimization]::Iteration Results:" << std::endl;
739 std::cout << "Zcut = " << zcut << std::endl;
740 std::cout << "Nsig = " << Nsig << std::endl;
741 std::cout << "n_bkg: " << Nbkg << std::endl;
742 std::cout << "n_on: " << n_on << std::endl;
743 std::cout << "n_off: " << n_off << std::endl;
744 std::cout << "ZBi: " << ZBi << std::endl;
745 */
746
747 //Update Test Cut with best scan values
748 if(ZBi > best_scan_zbi){
749 best_scan_zbi = ZBi;
750 best_scan_zcut = zcut;
751 best_scan_nsig = Nsig;
752 best_scan_nbkg = Nbkg;
753 }
754
755 //Fill TGraphs
756 zcutscan_zbi_g->SetPoint(zcutscan_zbi_g->GetN(),zcut, ZBi);
757 zcutscan_nbkg_g->SetPoint(zcutscan_nbkg_g->GetN(),zcut, Nbkg);
758 zcutscan_nsig_g->SetPoint(zcutscan_nsig_g->GetN(),zcut, Nsig);
759 nbkg_zbi_g->SetPoint(nbkg_zbi_g->GetN(),Nbkg, ZBi);
760 nsig_zbi_g->SetPoint(nsig_zbi_g->GetN(),Nsig, ZBi);
761
762 nsig_zcut_hh->Fill(zcut,Nsig,ZBi);
763 nbkg_zcut_hh->Fill(zcut,Nbkg,ZBi);
764 }
765
766 //Write graph of zcut vs zbi for the Test Cut
767 writeGraph(outFile_, "testCuts_pct_sig_cut_"+std::to_string(cutSignal),
768 zcutscan_zbi_g);
769 writeGraph(outFile_, "testCuts_pct_sig_cut_"+std::to_string(cutSignal),
770 zcutscan_nbkg_g);
771 writeGraph(outFile_, "testCuts_pct_sig_cut_"+std::to_string(cutSignal),
772 zcutscan_nsig_g);
773
774 //delete TGraph pointers
775 delete zcutscan_zbi_g;
776 delete zcutscan_nsig_g;
777 delete zcutscan_nbkg_g;
778 delete nbkg_zbi_g;
779 delete nsig_zbi_g;
780 delete nsig_zcut_hh;
781 delete nbkg_zcut_hh;
782
783 //Fill Summary Histograms Test Cuts at best zcutscan value
784 processorHistos_->Fill2DHisto("test_cuts_values_hh",(double)cutSignal, (double)cutid,cutvalue);
785 processorHistos_->set2DHistoYlabel("test_cuts_values_hh",cutid,cutname);
786
787 processorHistos_->Fill2DHisto("test_cuts_ZBi_hh",(double)cutSignal, (double)cutid,best_scan_zbi);
788 processorHistos_->set2DHistoYlabel("test_cuts_ZBi_hh",cutid,cutname);
789
790 processorHistos_->Fill2DHisto("test_cuts_zcut_hh",(double)cutSignal, (double)cutid,best_scan_zcut);
791 processorHistos_->set2DHistoYlabel("test_cuts_zcut_hh",cutid,cutname);
792
793 processorHistos_->Fill2DHisto("test_cuts_nsig_hh",(double)cutSignal, (double)cutid,best_scan_nsig);
794 processorHistos_->set2DHistoYlabel("test_cuts_nsig_hh",cutid,cutname);
795
796 processorHistos_->Fill2DHisto("test_cuts_nbkg_hh",(double)cutSignal, (double)cutid,best_scan_nbkg);
797 processorHistos_->set2DHistoYlabel("test_cuts_nbkg_hh",cutid,cutname);
798
799 //Check if the best cutscan zbi for this Test Cut is the best overall ZBi for all Test Cuts
800 if(best_scan_zbi > best_zbi){
801 best_zbi = best_scan_zbi;
802 best_cutname = cutname;
803 best_cutvalue = cutvalue;
804 best_zcut = best_scan_zcut;
805 best_nsig = best_scan_nsig;
806 best_nbkg = best_scan_nbkg;
807 }
808
809 } //END LOOP OVER TEST CUTS
810
811 if(debug_){
812 std::cout << "Iteration " << iteration << " Best Test Cut is " << best_cutname
813 << " " << best_cutvalue << " with ZBi=" << best_zbi << std::endl;
814 std::cout << "Update persistent cuts list with this best test cut..." << std::endl;
815 std::cout << "[Persistent Cuts] Before update:" << std::endl;
816 persistentCutsSelector_->printCuts();
817 }
818
819 //Find best Test Cut for iteration. Add Test Cut value to Persistent Cuts list
820 processorHistos_->Fill2DHisto("best_test_cut_ZBi_hh",(double)cutSignal, best_zbi,
821 (double)testCutsSelector_->getCutID(best_cutname));
822 processorHistos_->Fill1DHisto("best_test_cut_ZBi_h",(double)cutSignal, best_zbi);
823 processorHistos_->Fill1DHisto("best_test_cut_zcut_h",(double)cutSignal, best_zcut);
824 processorHistos_->Fill1DHisto("best_test_cut_nsig_h",(double)cutSignal, best_nsig);
825 processorHistos_->Fill1DHisto("best_test_cut_nbkg_h",(double)cutSignal, best_nbkg);
826 processorHistos_->Fill1DHisto("best_test_cut_id_h",(double)cutSignal,
827 (double)testCutsSelector_->getCutID(best_cutname));
828
829
830 persistentCutsSelector_->setCutValue(best_cutname, best_cutvalue);
831 if(debug_){
832 std::cout << "[Persistent Cuts] After update:" << std::endl;
833 persistentCutsSelector_->printCuts();
834 }
835
836 //Write iteration histos
837 signalHistos_->writeHistos(outFile_,"signal_pct_sig_cut_"+std::to_string(cutSignal));
838 bkgHistos_->writeHistos(outFile_,"background_pct_sig_cut_"+std::to_string(cutSignal));
839 testCutHistos_->writeHistos(outFile_, "testCuts_pct_sig_cut_"+std::to_string(cutSignal));
840 }
841}
842
843void SimpZBiOptimizationProcessor::finalize() {
844 std::cout << "[SimpZBiOptimizationProcessor] finalize()" << std::endl;
845
846 if(debug_){
847 std::cout << "FINAL LIST OF PERSISTENT CUTS " << std::endl;
848 persistentCutsSelector_->printCuts();
849 }
850
851 processorHistos_->saveHistos(outFile_);
852 testCutHistos_->writeGraphs(outFile_,"");
853}
854
855double SimpZBiOptimizationProcessor::calculateZBi(double n_on, double n_off, double tau){
856 double P_Bi = TMath::BetaIncomplete(1./(1.+tau),n_on,n_off+1);
857 std::cout << "P_Bi: " << P_Bi << std::endl;
858 double Z_Bi = std::pow(2,0.5)*TMath::ErfInverse(1-2*P_Bi);
859 std::cout << "Z_Bi: " << Z_Bi << std::endl;
860 return Z_Bi;
861}
862
863bool SimpZBiOptimizationProcessor::failPersistentCuts(SimpAnaTTree* MTT){
864 bool failCuts = false;
865 for(cut_iter_ it=persistentCutsPtr_->begin(); it!=persistentCutsPtr_->end(); it++){
866 std::string cutname = it->first;
867 std::string cutvar = persistentCutsSelector_->getCutVar(cutname);
868 //If no value inside the tuple exists for this cut, do not apply the cut.
869 if(!MTT->variableExists(cutvar))
870 continue;
871 if(!persistentCutsSelector_->passCutGTorLT(cutname, MTT->getValue(cutvar))){
872 failCuts = true;
873 break;
874 }
875 }
876
877 return failCuts;
878}
879
880bool SimpZBiOptimizationProcessor::failTestCut(std::string cutname, SimpAnaTTree* MTT){
881
882 std::string cutvar = testCutsSelector_->getCutVar(cutname);
883 double cutvalue = testCutsSelector_->getCut(cutname);
884 //If cut variable is not found in the list of tuples, do not apply cut
885 if(!MTT->variableExists(cutvar))
886 return false;
887 if(!testCutsSelector_->passCutGTorLT(cutname, MTT->getValue(cutvar)))
888 return true;
889 else
890 return false;
891}
892
893void SimpZBiOptimizationProcessor::getSignalMCAnaVtxZ_h(std::string signalMCAnaFilename,
894 std::string signal_pdgid){
895 //Read pre-trigger Signal MCAna vertex z distribution
896 signalSimZ_h_ = new TH1F("signal_SimZ_h_","signal_SimZ;true z_{vtx} [mm];events", 200, -50.3, 149.7);
897 TFile* signalMCAnaFile = new TFile(signalMCAnaFilename_.c_str(), "READ");
898 TH1F* mcAnaSimZ_h = (TH1F*)signalMCAnaFile->Get(("mcAna/mcAna_mc"+signal_pdgid+"Z_h").c_str());
899 for(int i=0; i < 201; i++){
900 signalSimZ_h_->SetBinContent(i,mcAnaSimZ_h->GetBinContent(i));
901 }
902 signalMCAnaFile->Close();
903 delete signalMCAnaFile;
904}
905
906void SimpZBiOptimizationProcessor::writeTH1F(TFile* outF, std::string folder, TH1F* h){
907 if (outF) outF->cd();
908 TDirectory* dir{nullptr};
909 std::cout<<folder.c_str()<<std::endl;
910 if (!folder.empty()) {
911 dir = outF->mkdir(folder.c_str(),"",true);
912 dir->cd();
913 }
914 h->Write();
915}
916
917void SimpZBiOptimizationProcessor::writeGraph(TFile* outF, std::string folder, TGraph* g){
918 if (outF) outF->cd();
919 TDirectory* dir{nullptr};
920 std::cout<<folder.c_str()<<std::endl;
921 if (!folder.empty()) {
922 dir = outF->mkdir(folder.c_str(),"",true);
923 dir->cd();
924 }
925 g->Write();
926}
927
928double SimpZBiOptimizationProcessor::round(double var){
929 //float value = (int)(var * 100 + .5);
930 //return (double)value / 100;
931 return var;
932}
933
934DECLARE_PROCESSOR(SimpZBiOptimizationProcessor);
DECLARE_PROCESSOR
#define DECLARE_PROCESSOR(CLASS)
Macro which allows the framework to construct a producer given its name during configuration.
Definition Processor.h:139
BaseSelector::getCut
float getCut(const std::string &cutname)
Get cut values ?
Definition BaseSelector.h:68
BaseSelector::LoadSelection
bool LoadSelection()
description
Definition BaseSelector.cxx:36
IterativeCutSelector::isCutGreaterThan
bool isCutGreaterThan(std::string cutname)
is cut of type 'greater than'
Definition IterativeCutSelector.cxx:17
IterativeCutSelector::getCutVar
std::string getCutVar(std::string cutname)
get cut variable from name
Definition IterativeCutSelector.cxx:25
IterativeCutSelector::getPointerToCuts
std::map< std::string, std::pair< double, int > > * getPointerToCuts()
get pointer to the base class cuts
Definition IterativeCutSelector.h:93
IterativeCutSelector::filterCuts
void filterCuts(std::vector< std::string > cut_variable_list)
remove cuts that aren't specified in the list of cut variables
Definition IterativeCutSelector.cxx:86
IterativeCutSelector::getCutID
int getCutID(std::string cutname)
get cut ID
Definition IterativeCutSelector.h:79
IterativeCutSelector::setCutValue
void setCutValue(std::string cutname, double value)
set cut value
Definition IterativeCutSelector.cxx:44
IterativeCutSelector::printCuts
void printCuts()
prints cuts and values
Definition IterativeCutSelector.cxx:81
IterativeCutSelector::passCutGTorLT
bool passCutGTorLT(std::string cutname, double val)
does value pass cut
Definition IterativeCutSelector.cxx:55
MutableTTree::GetEntry
void GetEntry(int entry)
get tree entry
Definition MutableTTree.h:31
MutableTTree::defineMassWindow
void defineMassWindow(double lowMass, double highMass)
Set the mass window within which to read the input ttree.
Definition MutableTTree.cxx:31
MutableTTree::getValue
double getValue(std::string branch_name)
Get the value of a flat tuple variable.
Definition MutableTTree.cxx:13
MutableTTree::shiftVariable
void shiftVariable(std::string variable, double shift)
Apply any corrections to specified variable.
Definition MutableTTree.cxx:81
MutableTTree::GetEntries
int GetEntries()
return number of entries in tree
Definition MutableTTree.h:25
MutableTTree::getAllVariables
std::vector< std::string > getAllVariables()
Get list of all variables defined in ttree.
Definition MutableTTree.cxx:43
MutableTTree::variableExists
bool variableExists(std::string variable)
Check if a variable exists in the ttree.
Definition MutableTTree.cxx:36
MutableTTree::Fill
void Fill()
Fill ttree with new variables included.
Definition MutableTTree.cxx:52
ParameterSet
description
Definition ParameterSet.h:22
Processor
Base class for all event processing components.
Definition Processor.h:34
SimpAnaTTree
Read in flat TTree and create new cut variables for SIMP analysis.
Definition SimpAnaTTree.h:17
SimpAnaTTree::addVariable_unc_vtx_pos_zalpha
void addVariable_unc_vtx_pos_zalpha(double slope)
Definition SimpAnaTTree.cxx:125
SimpAnaTTree::addVariable_unc_vtx_pos_z0tanlambda_right
void addVariable_unc_vtx_pos_z0tanlambda_right(double slope)
Definition SimpAnaTTree.cxx:320
SimpAnaTTree::addVariable_unc_vtx_ele_iso_z0err
void addVariable_unc_vtx_ele_iso_z0err()
Definition SimpAnaTTree.cxx:204
SimpAnaTTree::addVariable_unc_vtx_pos_isolation_cut
void addVariable_unc_vtx_pos_isolation_cut()
Definition SimpAnaTTree.cxx:270
SimpAnaTTree::addVariable_unc_vtx_zalpha_max
void addVariable_unc_vtx_zalpha_max(double slope)
Definition SimpAnaTTree.cxx:144
SimpAnaTTree::addVariable_unc_vtx_pos_z0tanlambda_left
void addVariable_unc_vtx_pos_z0tanlambda_left(double slope)
Definition SimpAnaTTree.cxx:344
SimpAnaTTree::addVariable_unc_vtx_ele_z0tanlambda_left
void addVariable_unc_vtx_ele_z0tanlambda_left(double slope)
Definition SimpAnaTTree.cxx:332
SimpAnaTTree::addVariable_unc_vtx_pos_z0tanlambda
void addVariable_unc_vtx_pos_z0tanlambda()
Definition SimpAnaTTree.cxx:296
SimpAnaTTree::addVariable_unc_vtx_ele_z0_z0err
void addVariable_unc_vtx_ele_z0_z0err()
Definition SimpAnaTTree.cxx:230
SimpAnaTTree::addVariable_unc_vtx_pos_iso_z0err
void addVariable_unc_vtx_pos_iso_z0err()
Definition SimpAnaTTree.cxx:217
SimpAnaTTree::addVariable_unc_vtx_ele_z0tanlambda_right
void addVariable_unc_vtx_ele_z0tanlambda_right(double slope)
Definition SimpAnaTTree.cxx:308
SimpAnaTTree::addVariable_unc_vtx_ele_isolation_cut
void addVariable_unc_vtx_ele_isolation_cut()
Definition SimpAnaTTree.cxx:256
SimpAnaTTree::addVariable_unc_vtx_abs_delta_z0tanlambda
void addVariable_unc_vtx_abs_delta_z0tanlambda()
Definition SimpAnaTTree.cxx:356
SimpAnaTTree::addVariable_unc_vtx_zalpha_min
void addVariable_unc_vtx_zalpha_min(double slope)
Definition SimpAnaTTree.cxx:174
SimpAnaTTree::addVariable_unc_vtx_ele_z0tanlambda
void addVariable_unc_vtx_ele_z0tanlambda()
Definition SimpAnaTTree.cxx:284
SimpAnaTTree::addVariable_unc_vtx_ele_zalpha
void addVariable_unc_vtx_ele_zalpha(double slope)
New Variables.
Definition SimpAnaTTree.cxx:106
SimpAnaTTree::addVariable_unc_vtx_pos_z0_z0err
void addVariable_unc_vtx_pos_z0_z0err()
Definition SimpAnaTTree.cxx:243
SimpEquations
All SIMP Equations for calculating expected signal.
Definition SimpEquations.h:13
SimpEquations::expectedSignalCalculation
double expectedSignalCalculation(double m_V, double eps, bool rho, double E_V, TEfficiency *effCalc_h, double target_pos, double zcut)
Definition SimpEquations.cxx:240
SimpEquations::getAprimeMassFromVectorMass
double getAprimeMassFromVectorMass(double m_V)
Definition SimpEquations.h:76
SimpEquations::massResolution
double massResolution(double m_V)
Definition SimpEquations.cxx:339
SimpZBiOptimizationProcessor
Cutflow optimization tool for SIMPS.
Definition SimpZBiOptimizationProcessor.h:31
SimpZBiOptimizationProcessor::calculateZBi
double calculateZBi(double n_on, double n_off, double tau)
description
Definition SimpZBiOptimizationProcessor.cxx:855
SimpZBiOptimizationProcessor::bkgHistos_
std::shared_ptr< ZBiHistos > bkgHistos_
Definition SimpZBiOptimizationProcessor.h:161
SimpZBiOptimizationProcessor::fillEventHistograms
void fillEventHistograms(std::shared_ptr< ZBiHistos > histos, SimpAnaTTree *MTT)
description
Definition SimpZBiOptimizationProcessor.cxx:126
SimpZBiOptimizationProcessor::cut_iter_
std::map< std::string, std::pair< double, int > >::iterator cut_iter_
Definition SimpZBiOptimizationProcessor.h:166
SimpZBiOptimizationProcessor::configure
virtual void configure(const ParameterSet &parameters)
description
Definition SimpZBiOptimizationProcessor.cxx:14
SimpZBiOptimizationProcessor::countControlRegionBackgroundRate
double countControlRegionBackgroundRate(std::string inFilename, std::string tree_name, double m_Ap, double Mbin=30.0, double dNdm_sf=1.0)
description
Definition SimpZBiOptimizationProcessor.cxx:229
SimpZBiOptimizationProcessor::getSignalMCAnaVtxZ_h
void getSignalMCAnaVtxZ_h(std::string signalMCAnaFilename, std::string signal_pdgid)
description
Definition SimpZBiOptimizationProcessor.cxx:893
SimpZBiOptimizationProcessor::failTestCut
bool failTestCut(std::string cutname, SimpAnaTTree *MTT)
description
Definition SimpZBiOptimizationProcessor.cxx:880
SimpZBiOptimizationProcessor::initialize
virtual void initialize(TTree *tree)
description
Definition SimpZBiOptimizationProcessor.h:47
SimpZBiOptimizationProcessor::new_variables_
std::vector< std::string > new_variables_
Definition SimpZBiOptimizationProcessor.h:129
SimpZBiOptimizationProcessor::failPersistentCuts
bool failPersistentCuts(SimpAnaTTree *MTT)
description
Definition SimpZBiOptimizationProcessor.cxx:863
SimpZBiOptimizationProcessor::writeTH1F
void writeTH1F(TFile *outF, std::string folder, TH1F *h)
description
Definition SimpZBiOptimizationProcessor.cxx:906
SimpZBiOptimizationProcessor::writeGraph
void writeGraph(TFile *outF, std::string folder, TGraph *g)
description
Definition SimpZBiOptimizationProcessor.cxx:917
SimpZBiOptimizationProcessor::testCutsPtr_
std::map< std::string, std::pair< double, int > > * testCutsPtr_
Definition SimpZBiOptimizationProcessor.h:171
SimpZBiOptimizationProcessor::processorHistos_
std::shared_ptr< ZBiHistos > processorHistos_
Definition SimpZBiOptimizationProcessor.h:163
SimpZBiOptimizationProcessor::addNewVariables
void addNewVariables(SimpAnaTTree *MTT, std::string variable, double param)
description
Definition SimpZBiOptimizationProcessor.cxx:72
SimpZBiOptimizationProcessor::persistentCutsPtr_
std::map< std::string, std::pair< double, int > > * persistentCutsPtr_
Definition SimpZBiOptimizationProcessor.h:173
SimpZBiOptimizationProcessor::signalHistos_
std::shared_ptr< ZBiHistos > signalHistos_
Definition SimpZBiOptimizationProcessor.h:160
SimpZBiOptimizationProcessor::initialIntegrals_
std::map< std::string, double > initialIntegrals_
Definition SimpZBiOptimizationProcessor.h:167
SimpZBiOptimizationProcessor::testCutHistos_
std::shared_ptr< ZBiHistos > testCutHistos_
Definition SimpZBiOptimizationProcessor.h:162
SimpZBiOptimizationProcessor::SimpZBiOptimizationProcessor
SimpZBiOptimizationProcessor(const std::string &name, Process &process)
Definition SimpZBiOptimizationProcessor.cxx:7
Generated by doxygen 1.9.8