Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Change the way total cross section is computed in GenXSecAnalzyer and fix the print out in LHERunInfo.cc #7788

Merged
merged 1 commit into from
Mar 13, 2015
Merged

Change the way total cross section is computed in GenXSecAnalzyer and fix the print out in LHERunInfo.cc #7788

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
51 changes: 25 additions & 26 deletions GeneratorInterface/Core/plugins/GenXSecAnalyzer.cc
View file
Open in desktop
Expand Up @@ -290,22 +290,22 @@ void
GenXSecAnalyzer::combine(double& finalValue, double& finalError, double& finalWeight, const double& currentValue, const double& currentError, const double & currentWeight)
{

if(finalValue<1e-10)
if(finalValue<=0)
{
finalValue = currentValue;
finalError = currentError;
finalWeight += currentWeight;
}
else
{
double wgt1 = (finalError < 1e-10 || currentError<1e-10)?
double wgt1 = (finalError <=0 || currentError <=0)?
finalWeight :
1/(finalError*finalError);
double wgt2 = (finalError < 1e-10 || currentError<1e-10)?
double wgt2 = (finalError <=0 || currentError <=0)?
currentWeight:
1/(currentError*currentError);
double xsec = (wgt1 * finalValue + wgt2 * currentValue) /(wgt1 + wgt2);
double err = (finalError < 1e-10 || currentError<1e-10)? 0 :
double err = (finalError <=0 || currentError <=0)? 0 :
1.0 / std::sqrt(wgt1 + wgt2);
finalValue = xsec;
finalError = err;
Expand Down Expand Up @@ -345,7 +345,7 @@ GenXSecAnalyzer::compute(const GenLumiInfoProduct& iLumiInfo)
for(unsigned int ip=0; ip < vectorSize; ip++){
GenLumiInfoProduct::ProcessInfo proc = iLumiInfo.getProcessInfos()[ip];
double hepxsec_value = proc.lheXSec().value();
double hepxsec_error = proc.lheXSec().error() < 1e-10? 0:proc.lheXSec().error();
double hepxsec_error = proc.lheXSec().error() <= 0? 0:proc.lheXSec().error();
tempVector_before.push_back(GenLumiInfoProduct::XSec(hepxsec_value,hepxsec_error));

sigSelSum += hepxsec_value;
Expand All @@ -364,14 +364,14 @@ GenXSecAnalyzer::compute(const GenLumiInfoProduct& iLumiInfo)
double npass = proc.nPassPos() -proc.nPassNeg();
switch(hepidwtup_){
case 3: case -3:
fracAcc = ntotal > 1e-6? npass/ntotal: -1;
fracAcc = ntotal > 0? npass/ntotal: -1;
break;
default:
fracAcc = proc.selected().sum() > 1e-6? proc.killed().sum() / proc.selected().sum():-1;
fracAcc = proc.selected().sum() > 0? proc.killed().sum() / proc.selected().sum():-1;
break;
}

if(fracAcc<1e-6)
if(fracAcc<=0)
{
tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0,0.0));
continue;
Expand All @@ -387,15 +387,15 @@ GenXSecAnalyzer::compute(const GenLumiInfoProduct& iLumiInfo)
case 3: case -3:
{
double ntotal_pos = proc.nTotalPos();
double effp = ntotal_pos > 1e-6?
double effp = ntotal_pos > 0?
(double)proc.nPassPos()/ntotal_pos:0;
double effp_err2 = ntotal_pos > 1e-6?
double effp_err2 = ntotal_pos > 0?
(1-effp)*effp/ntotal_pos: 0;

double ntotal_neg = proc.nTotalNeg();
double effn = ntotal_neg > 1e-6?
double effn = ntotal_neg > 0?
(double)proc.nPassNeg()/ntotal_neg:0;
double effn_err2 = ntotal_neg > 1e-6?
double effn_err2 = ntotal_neg > 0?
(1-effn)*effn/ntotal_neg: 0;

efferr2 = ntotal > 0 ?
Expand Down Expand Up @@ -442,7 +442,15 @@ GenXSecAnalyzer::compute(const GenLumiInfoProduct& iLumiInfo)

} // end of loop over different processes
tempVector_before.push_back(GenLumiInfoProduct::XSec(sigSelSum, sqrt(err2SelSum)));
GenLumiInfoProduct::XSec result(sigSum,std::sqrt(err2Sum));

double total_matcheff = jetMatchEffStat_[10000].filterEfficiency(hepidwtup_);
double total_matcherr = jetMatchEffStat_[10000].filterEfficiencyError(hepidwtup_);

double xsec_after = sigSelSum*total_matcheff;
double xsecerr_after = (total_matcheff > 0 && sigSelSum > 0)? xsec_after*sqrt(err2SelSum/sigSelSum/sigSelSum +
total_matcherr*total_matcherr/total_matcheff/total_matcheff):0;

GenLumiInfoProduct::XSec result(xsec_after,xsecerr_after);
tempVector_after.push_back(result);

xsecBeforeMatching_ =tempVector_before;
Expand Down Expand Up @@ -502,28 +510,19 @@ GenXSecAnalyzer::endJob() {
);


jetmatch_eff = thisJetMatchStat.filterEfficiency(hepidwtup_);
jetmatch_err = thisJetMatchStat.filterEfficiencyError(hepidwtup_);

if(i==last)
{
title[i] = "Total";

edm::LogPrint("GenXSecAnalyzer")
<< "-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ";

double n1 = xsecBeforeMatching_[i].value();
double e1 = xsecBeforeMatching_[i].error();
double n2 = xsecAfterMatching_[i].value();
double e2 = xsecAfterMatching_[i].error();

jetmatch_eff = n1>0? n2/n1 : 0;
jetmatch_err = (n1>0 && n2>0 && pow(e2/n2,2)>pow(e1/n1,2))?
jetmatch_eff*sqrt( pow(e2/n2,2) - pow(e1/n1,2)):-1;

<< "-------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ";
}
else
{
title[i] = Form("%d",i);
jetmatch_eff = thisJetMatchStat.filterEfficiency(hepidwtup_);
jetmatch_err = thisJetMatchStat.filterEfficiencyError(hepidwtup_);

}

Expand Down
48 changes: 23 additions & 25 deletions GeneratorInterface/LHEInterface/src/LHERunInfo.cc
View file
Open in desktop
Expand Up @@ -222,37 +222,36 @@ LHERunInfo::XSec LHERunInfo::xsec() const
double npass = proc->nPassPos() -proc->nPassNeg();
switch(idwtup){
case 3: case -3:
fracAcc = ntotal > 1e-6? npass/ntotal: -1;
fracAcc = ntotal > 0? npass/ntotal: -1;
break;
default:
fracAcc = proc->selected().sum() > 1e-6? proc->killed().sum() / proc->selected().sum():-1;
fracAcc = proc->selected().sum() > 0? proc->killed().sum() / proc->selected().sum():-1;
break;
}

if(fracAcc<1e-6)continue;
if(fracAcc<=0)continue;

double fracBr = proc->accepted().sum() > 0.0 ?
proc->acceptedBr().sum() / proc->accepted().sum() : 1;
double sigmaFin = sigmaAvg * fracAcc * fracBr;
double sigmaFin = sigmaAvg * fracAcc ;
double sigmaFinBr = sigmaFin * fracBr;

double relErr = 1.0;
if (proc->killed().n() > 1) {

double efferr2=0;
switch(idwtup) {
case 3: case -3:
{
double ntotal_pos = proc->nTotalPos();
double effp = ntotal_pos > 1e-6?
double effp = ntotal_pos > 0?
(double)proc->nPassPos()/ntotal_pos:0;
double effp_err2 = ntotal_pos > 1e-6?
double effp_err2 = ntotal_pos > 0?
(1-effp)*effp/ntotal_pos: 0;

double ntotal_neg = proc->nTotalNeg();
double effn = ntotal_neg > 1e-6?
double effn = ntotal_neg > 0?
(double)proc->nPassNeg()/ntotal_neg:0;
double effn_err2 = ntotal_neg > 1e-6?
double effn_err2 = ntotal_neg > 0?
(1-effn)*effn/ntotal_neg: 0;

efferr2 = ntotal > 0 ?
Expand All @@ -269,7 +268,7 @@ LHERunInfo::XSec LHERunInfo::xsec() const
double failw2 = proc->selected().sum2() - passw2;
double numerator = (passw2*failw*failw + failw2*passw*passw);

efferr2 = denominator>1e-6?
efferr2 = denominator > 0?
numerator/denominator:0;
break;
}
Expand All @@ -279,7 +278,7 @@ LHERunInfo::XSec LHERunInfo::xsec() const
+ sigma2Err / sigma2Sum;
relErr = (delta2Sum > 0.0 ?
std::sqrt(delta2Sum) : 0.0);
}

double deltaFin = sigmaFin * relErr;
double deltaFinBr = sigmaFinBr * relErr;

Expand Down Expand Up @@ -339,38 +338,37 @@ void LHERunInfo::statistics() const
double npass = proc->nPassPos() -proc->nPassNeg();
switch(idwtup){
case 3: case -3:
fracAcc = ntotal > 1e-6? npass/ntotal: -1;
fracAcc = ntotal > 0? npass/ntotal: -1;
break;
default:
fracAcc = proc->selected().sum() > 1e-6? proc->killed().sum() / proc->selected().sum():-1;
fracAcc = proc->selected().sum() > 0? proc->killed().sum() / proc->selected().sum():-1;
break;
}

if(fracAcc<1e-6)continue;

double fracBr = proc->accepted().sum() > 0.0 ?
proc->acceptedBr().sum() / proc->accepted().sum() : 1;
double sigmaFin = sigmaAvg * fracAcc;
double sigmaFin = fracAcc >0? sigmaAvg * fracAcc : 0;
double sigmaFinBr = sigmaFin * fracBr;

double relErr = 1.0;
double relAccErr = 1.0;
double efferr2=0;

if (proc->killed().n() > 1) {
if (proc->killed().n() > 0 && fracAcc > 0) {
switch(idwtup) {
case 3: case -3:
{
double ntotal_pos = proc->nTotalPos();
double effp = ntotal_pos > 1e-6?
double effp = ntotal_pos > 0?
(double)proc->nPassPos()/ntotal_pos:0;
double effp_err2 = ntotal_pos > 1e-6?
double effp_err2 = ntotal_pos > 0?
(1-effp)*effp/ntotal_pos: 0;

double ntotal_neg = proc->nTotalNeg();
double effn = ntotal_neg > 1e-6?
double effn = ntotal_neg > 0?
(double)proc->nPassNeg()/ntotal_neg:0;
double effn_err2 = ntotal_neg > 1e-6?
double effn_err2 = ntotal_neg > 0?
(1-effn)*effn/ntotal_neg: 0;

efferr2 = ntotal > 0 ?
Expand All @@ -387,7 +385,7 @@ void LHERunInfo::statistics() const
double failw2 = proc->selected().sum2() - passw2;
double numerator = (passw2*failw*failw + failw2*passw*passw);

efferr2 = denominator>1e-6?
efferr2 = denominator>0?
numerator/denominator:0;
break;
}
Expand All @@ -404,8 +402,8 @@ void LHERunInfo::statistics() const
double deltaFinBr = sigmaFinBr * relErr;

double ntotal_proc = proc->nTotalPos()+proc->nTotalNeg();
double event_eff_proc = ntotal_proc>1e-6? (double)(proc->nPassPos()+ proc->nPassNeg())/ntotal_proc: -1;
double event_eff_err_proc = ntotal_proc>1e-6? std::sqrt((1-event_eff_proc)*event_eff_proc/ntotal_proc): -1;
double event_eff_proc = ntotal_proc>0? (double)(proc->nPassPos()+ proc->nPassNeg())/ntotal_proc: -1;
double event_eff_err_proc = ntotal_proc>0? std::sqrt((1-event_eff_proc)*event_eff_proc/ntotal_proc): -1;

std::cout << proc->process() << "\t\t"
<< std::scientific << std::setprecision(3)
Expand Down Expand Up @@ -442,8 +440,8 @@ void LHERunInfo::statistics() const
}

double ntotal_all = (nTried_pos+nTried_neg);
double event_eff_all = ntotal_all>1e-6? (double)(nAccepted_pos+nAccepted_neg)/ntotal_all: -1;
double event_eff_err_all = ntotal_all>1e-6? std::sqrt((1-event_eff_all)*event_eff_all/ntotal_all): -1;
double event_eff_all = ntotal_all>0? (double)(nAccepted_pos+nAccepted_neg)/ntotal_all: -1;
double event_eff_err_all = ntotal_all>0? std::sqrt((1-event_eff_all)*event_eff_all/ntotal_all): -1;

std::cout << "Total\t\t"
<< std::scientific << std::setprecision(3)
Expand Down