9 hasHepMCFilterInfo_(
false),
11 filterOnlyEffStat_(0,0,0,0,0.,0.,0.,0.),
12 hepMCFilterEffStat_(0,0,0,0,0.,0.,0.,0.)
45 if(iLumi.
getByLabel(
"genFilterEfficiencyProducer", genFilter))
54 bool foundLumiInfo = iLumi.
getByLabel(
"generator",genLumiInfo);
55 if (!foundLumiInfo)
return;
59 std::vector<GenLumiInfoProduct::ProcessInfo> theProcesses = genLumiInfo->getProcessInfos();
66 edm::LogError(
"GenXSecAnalyzer::endLuminosityBlock") <<
"Pure parton shower has thisProcessInfos size!=1";
69 if(theProcesses[0].lheXSec().
value()<1
e-6){
70 edm::LogError(
"GenXSecAnalyzer::endLuminosityBlock") <<
"cross section value = " << theProcesses[0].lheXSec().value();
81 else if(!
products_[
i].samePhysics(*genLumiInfo))
83 edm::LogError(
"GenXSecAnalyzer::endLuminosityBlock") <<
"Merging samples that come from different physics processes";
118 double passw = temp_killed.
sum();
119 double passw2 = temp_killed.
sum2();
120 double totalw = temp_selected.
sum();
121 double totalw2 = temp_selected.
sum2();
124 theProcesses[ip].nPassPos(),
125 theProcesses[ip].nPassNeg(),
126 theProcesses[ip].nTotalPos(),
127 theProcesses[ip].nTotalNeg(),
135 theProcesses[ip].nPassPos(),
136 theProcesses[ip].nPassNeg(),
137 theProcesses[ip].nTotalPos(),
138 theProcesses[ip].nTotalNeg(),
149 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
151 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
153 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
154 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
155 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
156 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg())
160 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
162 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
164 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
165 theProcesses[ip].nPassPos()+theProcesses[ip].nPassNeg(),
166 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg(),
167 theProcesses[ip].nTotalPos()+theProcesses[ip].nTotalNeg())
191 sigSum += proc.
tried().
sum() * hepxsec_value;
195 double sigAve = totalN>1
e-6? sigSum/totalN: 0;
204 double sum_numerator_before[sizeOfInfos];
205 double sum_denominator_before[sizeOfInfos];
208 double sum_numerator_after[sizeOfInfos];
209 double sum_denominator_after[sizeOfInfos];
212 for(
unsigned int i=0;
i < sizeOfInfos;
i++)
214 sum_numerator_before[
i]=0;
215 sum_denominator_before[
i]=0;
216 sum_numerator_after[
i]=0;
217 sum_denominator_after[
i]=0;
224 double sigSelSum = 0.0;
225 double errSel2Sum = 0.0;
227 double err2Sum = 0.0;
230 unsigned int vectorSize =
products_[
i].getProcessInfos().size();
231 for(
unsigned int ip=0; ip < vectorSize; ip++){
246 fracAcc = ntotal > 1
e-6? npass/ntotal: -1;
253 if(fracAcc<1
e-6)
continue;
256 double sigmaFin = hepxsec_value * fracAcc;
265 double effp = ntotal_pos > 1
e-6?
266 (double)proc.
nPassPos()/ntotal_pos:0;
267 double effp_err2 = ntotal_pos > 1
e-6?
268 (1-effp)*effp/ntotal_pos: 0;
271 double effn = ntotal_neg > 1
e-6?
272 (double)proc.
nPassNeg()/ntotal_neg:0;
273 double effn_err2 = ntotal_neg > 1
e-6?
274 (1-effn)*effn/ntotal_neg: 0;
276 efferr2 = ntotal > 0 ?
277 (ntotal_pos*ntotal_pos*effp_err2 +
278 ntotal_neg*ntotal_neg*effn_err2)/ntotal/ntotal:0;
288 double numerator = (passw2*failw*failw + failw2*passw*passw);
290 efferr2 = denominator>1
e-6?
291 numerator/denominator:0;
295 double delta2Veto = efferr2/fracAcc/fracAcc;
299 double sigma2Sum, sigma2Err;
300 sigma2Sum = hepxsec_value * hepxsec_value;
301 sigma2Err = hepxsec_error * hepxsec_error;
304 sum_denominator_before[ip] += (sigma2Err > 0)? 1/sigma2Err: 0;
305 sum_numerator_before[ip] += (sigma2Err > 0)? hepxsec_value/sigma2Err: 0;
308 double delta2Sum = delta2Veto
309 + sigma2Err / sigma2Sum;
310 relErr = (delta2Sum > 0.0 ?
312 double deltaFin = sigmaFin * relErr;
315 sigSelSum += hepxsec_value;
316 errSel2Sum += sigma2Err;
318 err2Sum += deltaFin * deltaFin;
321 sum_denominator_after[ip] += (deltaFin > 0)? 1/(deltaFin * deltaFin): 0;
322 sum_numerator_after[ip] += (deltaFin > 0)? sigmaFin/(deltaFin * deltaFin): 0;
327 sum_denominator_before[sizeOfInfos-1] += (errSel2Sum> 0)? 1/errSel2Sum: 0;
328 sum_numerator_before[sizeOfInfos-1] += (errSel2Sum> 0)? sigSelSum/errSel2Sum: 0;
330 sum_denominator_after[sizeOfInfos-1] += (err2Sum>0)? 1/err2Sum: 0;
331 sum_numerator_after[sizeOfInfos-1] += (err2Sum>0)? sigSum/err2Sum: 0;
337 for(
unsigned int i=0;
i<sizeOfInfos;
i++)
339 double final_value = (sum_denominator_before[
i]>0) ? (sum_numerator_before[
i]/sum_denominator_before[
i]):0;
340 double final_error = (sum_denominator_before[
i]>0) ? (1/
sqrt(sum_denominator_before[i])):-1;
343 double final_value2 = (sum_denominator_after[
i]>0) ? (sum_numerator_after[i]/sum_denominator_after[i]):0;
344 double final_error2 = (sum_denominator_after[
i]>0) ? 1/
sqrt(sum_denominator_after[i]):-1;
357 <<
"------------------------------------" <<
"\n"
358 <<
"GenXsecAnalyzer:" <<
"\n"
359 <<
"------------------------------------";
362 edm::LogPrint(
"GenXSecAnalyzer") <<
"------------------------------------" <<
"\n"
363 <<
"Cross-section summary not available" <<
"\n"
364 <<
"------------------------------------";
373 <<
"-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- \n"
374 <<
"Overall cross-section summary:" <<
"\n"
375 <<
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------";
377 <<
"Process\t\txsec_before [pb]\t\tpassed\tnposw\tnnegw\ttried\tnposw\tnnegw \txsec_match [pb]\t\t\taccepted [%]\t event_eff [%]";
380 const int last = sizeOfInfos-1;
381 std::string
title[sizeOfInfos];
383 for(
int i=0;
i < sizeOfInfos;
i++){
385 double jetmatch_eff=0;
386 double jetmatch_err=0;
393 <<
"-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ";
402 std::ostringstream titless;
404 title[
i] = titless.str();
411 << title[
i] <<
"\t\t"
412 << std::scientific << std::setprecision(3)
421 << std::scientific << std::setprecision(3)
424 << std::fixed << std::setprecision(1)
425 << (jetmatch_eff*100) <<
" +/- " << (jetmatch_err*100) <<
"\t"
426 << std::fixed << std::setprecision(1)
433 <<
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------";
436 <<
"Before matching: total cross section = "
437 << std::scientific << std::setprecision(3)
444 <<
"After matching: total cross section = "
445 << std::scientific << std::setprecision(3)
446 << xsec_match <<
" +- " << error_match <<
" pb";
450 double hepMCFilter_eff = 1.0;
451 double hepMCFilter_err = 0.0;
457 <<
"HepMC filter efficiency (taking into account weights)= "
462 << std::scientific << std::setprecision(3)
463 << hepMCFilter_eff <<
" +- " << hepMCFilter_err;
467 double hepMCFilter_event_eff = hepMCFilter_event_total > 0 ? hepMCFilter_event_pass/ hepMCFilter_event_total : 0;
468 double hepMCFilter_event_err = hepMCFilter_event_total > 0 ?
469 sqrt((1-hepMCFilter_event_eff)*hepMCFilter_event_eff/hepMCFilter_event_total): -1;
471 <<
"HepMC filter efficiency (event-level)= "
472 <<
"(" << hepMCFilter_event_pass <<
")"
474 <<
"(" << hepMCFilter_event_total <<
")"
476 << std::scientific << std::setprecision(3)
477 << hepMCFilter_event_eff <<
" +- " << hepMCFilter_event_err;
485 <<
"Filter efficiency (taking into account weights)= "
490 << std::scientific << std::setprecision(3)
491 << filterOnly_eff <<
" +- " << filterOnly_err;
495 double filterOnly_event_eff = filterOnly_event_total > 0 ? filterOnly_event_pass/filterOnly_event_total : 0;
496 double filterOnly_event_err = filterOnly_event_total > 0 ?
497 sqrt((1-filterOnly_event_eff)*filterOnly_event_eff/filterOnly_event_total): -1;
499 <<
"Filter efficiency (event-level)= "
500 <<
"(" << filterOnly_event_pass <<
")"
502 <<
"(" << filterOnly_event_total <<
")"
504 << std::scientific << std::setprecision(3)
505 << filterOnly_event_eff <<
" +- " << filterOnly_event_err;
508 double xsec_final = xsec_match*filterOnly_eff*hepMCFilter_eff;
509 double error_final = xsec_final*
sqrt(error_match*error_match/xsec_match/xsec_match+
510 filterOnly_err*filterOnly_err/filterOnly_eff/filterOnly_eff +
511 hepMCFilter_err*hepMCFilter_err/hepMCFilter_eff/hepMCFilter_eff
515 <<
"After filter: final cross section = "
516 << std::scientific << std::setprecision(3)
517 << xsec_final <<
" +- " << error_final <<
" pb";
std::vector< GenLumiInfoProduct::XSec > xsecBeforeMatching_
std::vector< GenFilterInfo > jetMatchEffStat_
unsigned int nTotalNeg() const
std::vector< GenLumiInfoProduct > products_
unsigned int numTotalPositiveEvents() const
TrainProcessor *const proc
double filterEfficiency(int idwtup=+3) const
GenFilterInfo filterOnlyEffStat_
unsigned int nPassPos() const
bool getByLabel(std::string const &label, Handle< PROD > &result) 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 endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &)
unsigned int numTotalNegativeEvents() const
double filterEfficiencyError(int idwtup=+3) const
unsigned int numPassPositiveEvents() const
FinalStat selected() const
GenXSecAnalyzer(const edm::ParameterSet &)
std::vector< GenLumiInfoProduct::XSec > xsecAfterMatching_
virtual void analyze(const edm::Event &, const edm::EventSetup &)
unsigned int nTotalPos() const
double sumWeights() const
Power< A, B >::type pow(const A &a, const B &b)
unsigned int numPassNegativeEvents() const