9 hasHepMCFilterInfo_(
false),
11 filterOnlyEffStat_(0,0,0,0,0.,0.,0.,0.),
12 hepMCFilterEffStat_(0,0,0,0,0.,0.,0.,0.)
65 if (!genLumiInfo.
isValid())
return;
69 std::vector<GenLumiInfoProduct::ProcessInfo> theProcesses = genLumiInfo->getProcessInfos();
76 edm::LogError(
"GenXSecAnalyzer::endLuminosityBlock") <<
"Pure parton shower has thisProcessInfos size!=1";
79 if(theProcesses[0].lheXSec().
value()<1
e-6){
80 edm::LogError(
"GenXSecAnalyzer::endLuminosityBlock") <<
"cross section value = " << theProcesses[0].lheXSec().value();
91 else if(!
products_[
i].samePhysics(*genLumiInfo))
93 edm::LogError(
"GenXSecAnalyzer::endLuminosityBlock") <<
"Merging samples that come from different physics processes";
128 double passw = temp_killed.
sum();
129 double passw2 = temp_killed.
sum2();
130 double totalw = temp_selected.
sum();
131 double totalw2 = temp_selected.
sum2();
134 theProcesses[ip].nPassPos(),
135 theProcesses[ip].nPassNeg(),
136 theProcesses[ip].nTotalPos(),
137 theProcesses[ip].nTotalNeg(),
145 theProcesses[ip].nPassPos(),
146 theProcesses[ip].nPassNeg(),
147 theProcesses[ip].nTotalPos(),
148 theProcesses[ip].nTotalNeg(),
159 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
161 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
163 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
164 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
165 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
166 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg())
170 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
172 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
174 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
175 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
176 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
177 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg())
201 sigSum += proc.
tried().
sum() * hepxsec_value;
205 double sigAve = totalN>1
e-6? sigSum/totalN: 0;
214 double sum_numerator_before[sizeOfInfos];
215 double sum_denominator_before[sizeOfInfos];
218 double sum_numerator_after[sizeOfInfos];
219 double sum_denominator_after[sizeOfInfos];
222 for(
unsigned int i=0;
i < sizeOfInfos;
i++)
224 sum_numerator_before[
i]=0;
225 sum_denominator_before[
i]=0;
226 sum_numerator_after[
i]=0;
227 sum_denominator_after[
i]=0;
234 double sigSelSum = 0.0;
235 double errSel2Sum = 0.0;
237 double err2Sum = 0.0;
240 unsigned int vectorSize =
products_[
i].getProcessInfos().size();
241 for(
unsigned int ip=0; ip < vectorSize; ip++){
256 fracAcc = ntotal > 1
e-6? npass/ntotal: -1;
263 if(fracAcc<1
e-6)
continue;
266 double sigmaFin = hepxsec_value * fracAcc;
275 double effp = ntotal_pos > 1
e-6?
276 (double)proc.
nPassPos()/ntotal_pos:0;
277 double effp_err2 = ntotal_pos > 1
e-6?
278 (1-effp)*effp/ntotal_pos: 0;
281 double effn = ntotal_neg > 1
e-6?
282 (double)proc.
nPassNeg()/ntotal_neg:0;
283 double effn_err2 = ntotal_neg > 1
e-6?
284 (1-effn)*effn/ntotal_neg: 0;
286 efferr2 = ntotal > 0 ?
287 (ntotal_pos*ntotal_pos*effp_err2 +
288 ntotal_neg*ntotal_neg*effn_err2)/ntotal/ntotal:0;
298 double numerator = (passw2*failw*failw + failw2*passw*passw);
300 efferr2 = denominator>1
e-6?
301 numerator/denominator:0;
305 double delta2Veto = efferr2/fracAcc/fracAcc;
309 double sigma2Sum, sigma2Err;
310 sigma2Sum = hepxsec_value * hepxsec_value;
311 sigma2Err = hepxsec_error * hepxsec_error;
314 sum_denominator_before[ip] += (sigma2Err > 0)? 1/sigma2Err: 0;
315 sum_numerator_before[ip] += (sigma2Err > 0)? hepxsec_value/sigma2Err: 0;
318 double delta2Sum = delta2Veto
319 + sigma2Err / sigma2Sum;
320 relErr = (delta2Sum > 0.0 ?
322 double deltaFin = sigmaFin * relErr;
325 sigSelSum += hepxsec_value;
326 errSel2Sum += sigma2Err;
328 err2Sum += deltaFin * deltaFin;
331 sum_denominator_after[ip] += (deltaFin > 0)? 1/(deltaFin * deltaFin): 0;
332 sum_numerator_after[ip] += (deltaFin > 0)? sigmaFin/(deltaFin * deltaFin): 0;
337 sum_denominator_before[sizeOfInfos-1] += (errSel2Sum> 0)? 1/errSel2Sum: 0;
338 sum_numerator_before[sizeOfInfos-1] += (errSel2Sum> 0)? sigSelSum/errSel2Sum: 0;
340 sum_denominator_after[sizeOfInfos-1] += (err2Sum>0)? 1/err2Sum: 0;
341 sum_numerator_after[sizeOfInfos-1] += (err2Sum>0)? sigSum/err2Sum: 0;
347 for(
unsigned int i=0;
i<sizeOfInfos;
i++)
349 double final_value = (sum_denominator_before[
i]>0) ? (sum_numerator_before[
i]/sum_denominator_before[
i]):0;
350 double final_error = (sum_denominator_before[
i]>0) ? (1/
sqrt(sum_denominator_before[i])):-1;
353 double final_value2 = (sum_denominator_after[
i]>0) ? (sum_numerator_after[i]/sum_denominator_after[i]):0;
354 double final_error2 = (sum_denominator_after[
i]>0) ? 1/
sqrt(sum_denominator_after[i]):-1;
367 <<
"------------------------------------" <<
"\n"
368 <<
"GenXsecAnalyzer:" <<
"\n"
369 <<
"------------------------------------";
372 edm::LogPrint(
"GenXSecAnalyzer") <<
"------------------------------------" <<
"\n"
373 <<
"Cross-section summary not available" <<
"\n"
374 <<
"------------------------------------";
383 <<
"-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- \n"
384 <<
"Overall cross-section summary:" <<
"\n"
385 <<
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------";
387 <<
"Process\t\txsec_before [pb]\t\tpassed\tnposw\tnnegw\ttried\tnposw\tnnegw \txsec_match [pb]\t\t\taccepted [%]\t event_eff [%]";
390 const int last = sizeOfInfos-1;
393 for(
int i=0;
i < sizeOfInfos;
i++){
395 double jetmatch_eff=0;
396 double jetmatch_err=0;
403 <<
"-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ";
412 title[
i] = Form(
"%d",
i);
419 << title[
i] <<
"\t\t"
420 << std::scientific << std::setprecision(3)
429 << std::scientific << std::setprecision(3)
432 << std::fixed << std::setprecision(1)
433 << (jetmatch_eff*100) <<
" +/- " << (jetmatch_err*100) <<
"\t"
434 << std::fixed << std::setprecision(1)
442 <<
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------";
445 <<
"Before matching: total cross section = "
446 << std::scientific << std::setprecision(3)
453 <<
"After matching: total cross section = "
454 << std::scientific << std::setprecision(3)
455 << xsec_match <<
" +- " << error_match <<
" pb";
459 double hepMCFilter_eff = 1.0;
460 double hepMCFilter_err = 0.0;
466 <<
"HepMC filter efficiency (taking into account weights)= "
471 << std::scientific << std::setprecision(3)
472 << hepMCFilter_eff <<
" +- " << hepMCFilter_err;
476 double hepMCFilter_event_eff = hepMCFilter_event_total > 0 ? hepMCFilter_event_pass/ hepMCFilter_event_total : 0;
477 double hepMCFilter_event_err = hepMCFilter_event_total > 0 ?
478 sqrt((1-hepMCFilter_event_eff)*hepMCFilter_event_eff/hepMCFilter_event_total): -1;
480 <<
"HepMC filter efficiency (event-level)= "
481 <<
"(" << hepMCFilter_event_pass <<
")"
483 <<
"(" << hepMCFilter_event_total <<
")"
485 << std::scientific << std::setprecision(3)
486 << hepMCFilter_event_eff <<
" +- " << hepMCFilter_event_err;
494 <<
"Filter efficiency (taking into account weights)= "
499 << std::scientific << std::setprecision(3)
500 << filterOnly_eff <<
" +- " << filterOnly_err;
504 double filterOnly_event_eff = filterOnly_event_total > 0 ? filterOnly_event_pass/filterOnly_event_total : 0;
505 double filterOnly_event_err = filterOnly_event_total > 0 ?
506 sqrt((1-filterOnly_event_eff)*filterOnly_event_eff/filterOnly_event_total): -1;
508 <<
"Filter efficiency (event-level)= "
509 <<
"(" << filterOnly_event_pass <<
")"
511 <<
"(" << filterOnly_event_total <<
")"
513 << std::scientific << std::setprecision(3)
514 << filterOnly_event_eff <<
" +- " << filterOnly_event_err;
517 double xsec_final = xsec_match*filterOnly_eff*hepMCFilter_eff;
518 double error_final = xsec_final*
sqrt(error_match*error_match/xsec_match/xsec_match+
519 filterOnly_err*filterOnly_err/filterOnly_eff/filterOnly_eff +
520 hepMCFilter_err*hepMCFilter_err/hepMCFilter_eff/hepMCFilter_eff
524 <<
"After filter: final cross section = "
525 << std::scientific << std::setprecision(3)
526 << xsec_final <<
" +- " << error_final <<
" pb";
std::vector< GenLumiInfoProduct::XSec > xsecBeforeMatching_
std::vector< GenFilterInfo > jetMatchEffStat_
unsigned int nTotalNeg() const
std::vector< GenLumiInfoProduct > products_
virtual void endJob() override
edm::EDGetTokenT< GenFilterInfo > hepMCFilterInfoToken_
virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) override
unsigned int numTotalPositiveEvents() const
virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) override
bool getByToken(EDGetToken token, Handle< PROD > &result) const
TrainProcessor *const proc
double filterEfficiency(int idwtup=+3) const
GenFilterInfo filterOnlyEffStat_
edm::EDGetTokenT< GenLumiInfoProduct > genLumiInfoToken_
unsigned int nPassPos() const
unsigned int theProcesses_size
GenLumiInfoProduct::XSec xsec_
GenFilterInfo hepMCFilterEffStat_
unsigned int nPassNeg() const
bool mergeProduct(GenFilterInfo const &other)
std::vector< GenFilterInfo > eventEffStat_
double sumPassWeights() const
virtual void analyze(const edm::Event &, const edm::EventSetup &) override
virtual void beginJob() override
unsigned int numTotalNegativeEvents() const
double filterEfficiencyError(int idwtup=+3) const
unsigned int numPassPositiveEvents() const
FinalStat selected() const
edm::EDGetTokenT< GenFilterInfo > genFilterInfoToken_
GenXSecAnalyzer(const edm::ParameterSet &)
std::vector< GenLumiInfoProduct::XSec > xsecAfterMatching_
volatile std::atomic< bool > shutdown_flag false
unsigned int nTotalPos() const
double sumWeights() const
Power< A, B >::type pow(const A &a, const B &b)
unsigned int numPassNegativeEvents() const