14 xsecPreFilter_(-1,-1),
17 filterOnlyEffRun_(0,0,0,0,0.,0.,0.,0.),
18 hepMCFilterEffRun_(0,0,0,0,0.,0.,0.,0.),
19 filterOnlyEffStat_(0,0,0,0,0.,0.,0.,0.),
20 hepMCFilterEffStat_(0,0,0,0,0.,0.,0.,0.)
96 if (!genLumiInfo.
isValid())
return;
99 std::vector<GenLumiInfoProduct::ProcessInfo> theProcesses = genLumiInfo->
getProcessInfos();
101 unsigned int theProcesses_size = theProcesses.size();
107 if(theProcesses_size!=1){
108 edm::LogError(
"GenXSecAnalyzer::endLuminosityBlock") <<
"Pure parton shower has thisProcessInfos size!=1";
114 for(
unsigned int ip=0; ip < theProcesses_size; ip++)
117 if(theProcesses[ip].lheXSec().
value()<0){
118 edm::LogError(
"GenXSecAnalyzer::endLuminosityBlock") <<
"cross section of process " << ip <<
" value = " << theProcesses[ip].lheXSec().value();
147 for(
unsigned int ip=0; ip < theProcesses_size; ip++)
149 int id = theProcesses[ip].process();
154 double passw = temp_killed.
sum();
155 double passw2 = temp_killed.
sum2();
156 double totalw = temp_selected.
sum();
157 double totalw2 = temp_selected.
sum2();
159 theProcesses[ip].nPassPos(),
160 theProcesses[ip].nPassNeg(),
161 theProcesses[ip].nTotalPos(),
162 theProcesses[ip].nTotalNeg(),
170 double currentValue = theProcesses[ip].lheXSec().value();
171 double currentError = theProcesses[ip].lheXSec().error();
180 if(currentValue != previousValue)
183 double xsec = y.
value();
184 double err = y.
error();
196 y = theProcesses[ip].lheXSec();
219 for (
unsigned int iSize = 0 ; iSize < thisHeprup_.
XSECUP.size() ; iSize++ ) {
234 std::vector<GenLumiInfoProduct::ProcessInfo> newInfos;
239 temp.
setLheXSec(iter->second.value(),iter->second.error());
240 newInfos.push_back(temp);
248 double thisHepFilterEff = 1;
249 double thisHepFilterErr = 0;
255 if(thisHepFilterEff<0)
257 thisHepFilterEff = 1;
258 thisHepFilterErr = 0;
263 double thisGenFilterEff = 1;
264 double thisGenFilterErr = 0;
270 if(thisGenFilterEff<0)
272 thisGenFilterEff = 1;
273 thisGenFilterErr = 0;
277 double thisXsec = thisRunXSecPre.
value() > 0 ? thisHepFilterEff*thisGenFilterEff*thisRunXSecPre.
value() : 0;
278 double thisErr = thisRunXSecPre.
value() > 0 ? thisXsec*
280 pow(thisHepFilterErr/thisHepFilterEff,2)+
281 pow(thisGenFilterErr/thisGenFilterEff,2)) : 0;
290 GenXSecAnalyzer::combine(
double& finalValue,
double& finalError,
double& finalWeight,
const double& currentValue,
const double& currentError,
const double & currentWeight)
295 finalValue = currentValue;
296 finalError = currentError;
297 finalWeight += currentWeight;
301 double wgt1 = (finalError <=0 || currentError <=0)?
303 1/(finalError*finalError);
304 double wgt2 = (finalError <=0 || currentError <=0)?
306 1/(currentError*currentError);
307 double xsec = (wgt1 * finalValue + wgt2 * currentValue) /(wgt1 + wgt2);
308 double err = (finalError <=0 || currentError <=0)? 0 :
312 finalWeight += currentWeight;
323 double thisValue = thisRunXSec.
value();
324 double thisError = thisRunXSec.
error();
325 combine(value,error,totalw,thisValue,thisError,thisw);
335 double sigSelSum = 0.0;
336 double err2SelSum = 0.0;
338 double err2Sum = 0.0;
340 std::vector<GenLumiInfoProduct::XSec> tempVector_before;
341 std::vector<GenLumiInfoProduct::XSec> tempVector_after;
345 for(
unsigned int ip=0; ip < vectorSize; ip++){
351 sigSelSum += hepxsec_value;
352 err2SelSum += hepxsec_error*hepxsec_error;
367 fracAcc = ntotal > 0? npass/ntotal: -1;
381 double sigmaFin = hepxsec_value * fracAcc;
390 double effp = ntotal_pos > 0?
391 (double)proc.
nPassPos()/ntotal_pos:0;
392 double effp_err2 = ntotal_pos > 0?
393 (1-effp)*effp/ntotal_pos: 0;
396 double effn = ntotal_neg > 0?
397 (double)proc.
nPassNeg()/ntotal_neg:0;
398 double effn_err2 = ntotal_neg > 0?
399 (1-effn)*effn/ntotal_neg: 0;
401 efferr2 = ntotal > 0 ?
402 (ntotal_pos*ntotal_pos*effp_err2 +
403 ntotal_neg*ntotal_neg*effn_err2)/ntotal/ntotal:0;
413 double numerator = (passw2*failw*failw + failw2*passw*passw);
415 efferr2 = denominator>0?
416 numerator/denominator:0;
420 double delta2Veto = efferr2/fracAcc/fracAcc;
424 double sigma2Sum, sigma2Err;
425 sigma2Sum = hepxsec_value * hepxsec_value;
426 sigma2Err = hepxsec_error * hepxsec_error;
429 double delta2Sum = delta2Veto
430 + sigma2Err / sigma2Sum;
431 relErr = (delta2Sum > 0.0 ?
433 double deltaFin = sigmaFin * relErr;
439 err2Sum += deltaFin * deltaFin;
449 double xsec_after = sigSelSum*total_matcheff;
450 double xsecerr_after = (total_matcheff > 0 && sigSelSum > 0)? xsec_after*
sqrt(err2SelSum/sigSelSum/sigSelSum +
451 total_matcherr*total_matcherr/total_matcheff/total_matcheff):0;
454 tempVector_after.push_back(result);
467 <<
"------------------------------------" <<
"\n" 468 <<
"GenXsecAnalyzer:" <<
"\n" 469 <<
"------------------------------------";
472 edm::LogPrint(
"GenXSecAnalyzer") <<
"------------------------------------" <<
"\n" 473 <<
"Cross-section summary not available" <<
"\n" 474 <<
"------------------------------------";
484 <<
"-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- \n" 485 <<
"Overall cross-section summary \n" 486 <<
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------";
488 <<
"Process\t\txsec_before [pb]\t\tpassed\tnposw\tnnegw\ttried\tnposw\tnnegw \txsec_match [pb]\t\t\taccepted [%]\t event_eff [%]";
491 const unsigned last = sizeOfInfos-1;
494 double jetmatch_eff=0;
495 double jetmatch_err=0;
497 for(std::map<int, GenFilterInfo>::const_iterator iter =
jetMatchEffStat_.begin();
521 <<
"-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ";
525 title[
i] = Form(
"%d",i);
531 << title[
i] <<
"\t\t" 532 << std::scientific << std::setprecision(3)
541 << std::scientific << std::setprecision(3)
545 << (jetmatch_eff*100) <<
" +/- " << (jetmatch_err*100) <<
"\t" 554 <<
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------";
557 <<
"Before matching: total cross section = " 558 << std::scientific << std::setprecision(3)
562 <<
"After matching: total cross section = " 563 << std::scientific << std::setprecision(3)
568 <<
"Before Filtrer: total cross section = " 569 << std::scientific << std::setprecision(3)
573 double hepMCFilter_eff = 1.0;
574 double hepMCFilter_err = 0.0;
579 <<
"HepMC filter efficiency (taking into account weights)= " 584 << std::scientific << std::setprecision(3)
585 << hepMCFilter_eff <<
" +- " << hepMCFilter_err;
589 double hepMCFilter_event_eff = hepMCFilter_event_total > 0 ? hepMCFilter_event_pass/ hepMCFilter_event_total : 0;
590 double hepMCFilter_event_err = hepMCFilter_event_total > 0 ?
591 sqrt((1-hepMCFilter_event_eff)*hepMCFilter_event_eff/hepMCFilter_event_total): -1;
593 <<
"HepMC filter efficiency (event-level)= " 594 <<
"(" << hepMCFilter_event_pass <<
")" 596 <<
"(" << hepMCFilter_event_total <<
")" 598 << std::scientific << std::setprecision(3)
599 << hepMCFilter_event_eff <<
" +- " << hepMCFilter_event_err;
608 <<
"Filter efficiency (taking into account weights)= " 613 << std::scientific << std::setprecision(3)
614 << filterOnly_eff <<
" +- " << filterOnly_err;
618 double filterOnly_event_eff = filterOnly_event_total > 0 ? filterOnly_event_pass/filterOnly_event_total : 0;
619 double filterOnly_event_err = filterOnly_event_total > 0 ?
620 sqrt((1-filterOnly_event_eff)*filterOnly_event_eff/filterOnly_event_total): -1;
622 <<
"Filter efficiency (event-level)= " 623 <<
"(" << filterOnly_event_pass <<
")" 625 <<
"(" << filterOnly_event_total <<
")" 627 << std::scientific << std::setprecision(3)
628 << filterOnly_event_eff <<
" +- " << filterOnly_event_err;
633 <<
"After filter: final cross section = " 634 << std::scientific << std::setprecision(3)
std::vector< GenLumiInfoProduct::XSec > xsecBeforeMatching_
const std::vector< ProcessInfo > & getProcessInfos() const
std::map< int, GenFilterInfo > jetMatchEffStat_
unsigned int nTotalNeg() const
virtual void endJob() override
void combine(GenLumiInfoProduct::XSec &, double &, const GenLumiInfoProduct::XSec &, const double &)
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
std::map< int, GenLumiInfoProduct::XSec > currentLumiBlockLHEXSec_
const lhef::HEPRUP & heprup() const
GenLumiInfoProduct::XSec xsec_
GenFilterInfo hepMCFilterEffRun_
GenFilterInfo hepMCFilterEffStat_
unsigned int nPassNeg() const
bool mergeProduct(GenFilterInfo const &other)
double sumPassWeights() const
unsigned int numEventsPassed() const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
virtual void analyze(const edm::Event &, const edm::EventSetup &) override
virtual void beginJob() override
virtual void beginRun(edm::Run const &, edm::EventSetup const &) override
GenLumiInfoProduct::XSec compute(const GenLumiInfoProduct &)
unsigned int numTotalNegativeEvents() const
std::vector< double > XERRUP
unsigned int numEventsTotal() const
std::map< int, GenLumiInfoProduct::XSec > previousLumiBlockLHEXSec_
std::vector< double > XMAXUP
double filterEfficiencyError(int idwtup=+3) const
virtual bool mergeProduct(const GenLumiInfoProduct &other)
unsigned int numPassPositiveEvents() const
void setLheXSec(double value, double err)
double sumWeights2() const
GenFilterInfo filterOnlyEffRun_
FinalStat selected() const
edm::EDGetTokenT< GenFilterInfo > genFilterInfoToken_
GenXSecAnalyzer(const edm::ParameterSet &)
std::vector< GenLumiInfoProduct::XSec > xsecAfterMatching_
virtual void endRun(edm::Run const &, edm::EventSetup const &) override
edm::EDGetTokenT< LHERunInfoProduct > lheRunInfoToken_
std::vector< double > XSECUP
unsigned int nTotalPos() const
double sumWeights() const
Power< A, B >::type pow(const A &a, const B &b)
void setProcessInfo(const std::vector< ProcessInfo > &processes)
GenLumiInfoProduct::XSec xsecPreFilter_
GenLumiInfoProduct product_
const int getHEPIDWTUP() const
unsigned int numPassNegativeEvents() const