|
|
[Rivet-svn] r4273 - trunk/src/Analyses
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Mon May 13 21:19:10 BST 2013
Author: buckley
Date: Mon May 13 21:19:10 2013
New Revision: 4273
Log:
Addressing a few YODA TODOs in analyses
Modified:
trunk/src/Analyses/ATLAS_2011_I945498.cc
trunk/src/Analyses/ATLAS_2012_I1083318.cc
Modified: trunk/src/Analyses/ATLAS_2011_I945498.cc
==============================================================================
--- trunk/src/Analyses/ATLAS_2011_I945498.cc Mon May 13 21:18:12 2013 (r4272)
+++ trunk/src/Analyses/ATLAS_2011_I945498.cc Mon May 13 21:19:10 2013 (r4273)
@@ -29,7 +29,6 @@
ATLAS_2011_I945498()
: Analysis("ATLAS_2011_I945498")
{
- /// @todo Set whether your finalize method needs the generator cross section
setNeedsCrossSection(true);
for (size_t chn = 0; chn < 3; ++chn) {
weights_nj0[chn] = 0.0;
@@ -45,7 +44,6 @@
public:
-
/// Book histograms and initialise projections before the run
void init() {
@@ -54,9 +52,9 @@
addProjection(zfinder_mu, "ZFinder_mu");
std::vector<std::pair<double, double> > eta_e;
- eta_e.push_back(make_pair(-2.47,-1.52));
- eta_e.push_back(make_pair(-1.37,1.37));
- eta_e.push_back(make_pair(1.52,2.47));
+ eta_e.push_back(make_pair(-2.47, -1.52));
+ eta_e.push_back(make_pair(-1.37, 1.37));
+ eta_e.push_back(make_pair(1.52, 2.47));
ZFinder zfinder_el(eta_e, 20, PID::ELECTRON, 66.0*GeV, 116.0*GeV, 0.1, true, false);
addProjection(zfinder_el, "ZFinder_el");
@@ -86,12 +84,13 @@
}
}
+
// Jet selection criteria universal for electron and muon channel
+ /// @todo Replace with a Cut passed to jetsByPt
Jets selectJets(const ZFinder* zf, const Event& event) {
- FourMomentum l1=zf->constituents()[0].momentum();
- FourMomentum l2=zf->constituents()[1].momentum();
+ FourMomentum l1 = zf->constituents()[0].momentum();
+ FourMomentum l2 = zf->constituents()[1].momentum();
Jets jets;
-
foreach (const Jet& jet, applyProjection<FastJets>(event, "jets").jetsByPt(30.0*GeV)) {
FourMomentum jmom = jet.momentum();
if (fabs(jmom.rapidity()) < 4.4 && deltaR(l1, jmom) > 0.5 && deltaR(l2, jmom) > 0.5) {
@@ -101,6 +100,7 @@
return jets;
}
+
/// Perform the per-event analysis
void analyze(const Event& event) {
@@ -110,14 +110,12 @@
zfs.push_back(& (applyProjection<ZFinder>(event, "ZFinder_el")));
zfs.push_back(& (applyProjection<ZFinder>(event, "ZFinder_mu")));
-
// Require exactly one electronic or muonic Z-decay in the event
- if (!( (zfs[0]->bosons().size()==1 && zfs[1]->bosons().size()!=1) || (zfs[1]->bosons().size()==1 && zfs[0]->bosons().size()!=1) )) vetoEvent;
-
- int chn = (zfs[0]->bosons().size()==1 ? 0:1);
-
- Jets jets=selectJets(zfs[chn], event);
+ if (!( (zfs[0]->bosons().size() == 1 && zfs[1]->bosons().size() != 1) ||
+ (zfs[1]->bosons().size() == 1 && zfs[0]->bosons().size() != 1) )) vetoEvent;
+ int chn = zfs[0]->bosons().size() == 1 ? 0 : 1;
+ Jets jets = selectJets(zfs[chn], event);
/// TODO: Holger wrote in his commit message that the njet_ratio
/// histograms need fixing/checking, and therefore the analysis
@@ -126,25 +124,26 @@
// Some silly weight counters for the njet-ratio histo
// --- not sure about the njet=0 case, the Figure ca[ption says
// that selected events require at least one jet with 20 GeV
- if (jets.size() == 0) {
+ switch (jets.size()) {
+ case 0:
weights_nj0[chn] += 1.0;
- }
- else if (jets.size() == 1) {
+ break;
+ case 1:
weights_nj0[chn] += 1.0;
weights_nj1[chn] += 1.0;
- }
- else if (jets.size() == 2) {
+ break;
+ case 2:
weights_nj0[chn] += 1.0;
weights_nj1[chn] += 1.0;
weights_nj2[chn] += 1.0;
- }
- else if (jets.size() == 3) {
+ break;
+ case 3:
weights_nj0[chn] += 1.0;
weights_nj1[chn] += 1.0;
weights_nj2[chn] += 1.0;
weights_nj3[chn] += 1.0;
- }
- else if (jets.size() >= 4) {
+ break;
+ default: // >= 4
weights_nj0[chn] += 1.0;
weights_nj1[chn] += 1.0;
weights_nj2[chn] += 1.0;
@@ -152,108 +151,97 @@
weights_nj4[chn] += 1.0;
}
- if (jets.size() <1) vetoEvent;
+ // Require at least one jet
+ if (jets.size() < 1) vetoEvent;
- _h_njet_incl[chn] ->fill(jets.size(), weight);
- _h_njet_incl[2] ->fill(jets.size(), weight);
+ // Fill jet multiplicities
+ _h_njet_incl[chn]->fill(jets.size(), weight);
+ _h_njet_incl[2]->fill(jets.size(), weight);
// Loop over selected jets, fill inclusive jet distributions
- for (unsigned int ijet=0; ijet<jets.size(); ijet++){
- _h_ptjet[chn] ->fill(jets[ijet].momentum().pT()/GeV, weight);
- _h_ptjet[2] ->fill(jets[ijet].momentum().pT()/GeV, weight);
- _h_yjet[chn] ->fill(fabs(jets[ijet].momentum().rapidity()), weight);
- _h_yjet[2] ->fill(fabs(jets[ijet].momentum().rapidity()), weight);
+ for (size_t ijet = 0; ijet < jets.size(); ++ijet) {
+ _h_ptjet[chn]->fill(jets[ijet].momentum().pT()/GeV, weight);
+ _h_ptjet[2] ->fill(jets[ijet].momentum().pT()/GeV, weight);
+ _h_yjet[chn] ->fill(fabs(jets[ijet].momentum().rapidity()), weight);
+ _h_yjet[2] ->fill(fabs(jets[ijet].momentum().rapidity()), weight);
+ }
+
+ // Leading jet histos
+ const double ptlead = jets[0].momentum().pT()/GeV;
+ const double yabslead = fabs(jets[0].momentum().rapidity());
+ _h_ptlead[chn]->fill(ptlead, weight);
+ _h_ptlead[2] ->fill(ptlead, weight);
+ _h_ylead[chn] ->fill(yabslead, weight);
+ _h_ylead[2] ->fill(yabslead, weight);
+
+ if (jets.size() >= 2) {
+ // Second jet histos
+ const double pt2ndlead = jets[1].momentum().pT()/GeV;
+ const double yabs2ndlead = fabs(jets[1].momentum().rapidity());
+ _h_ptseclead[chn] ->fill(pt2ndlead, weight);
+ _h_ptseclead[2] ->fill(pt2ndlead, weight);
+ _h_yseclead[chn] ->fill(yabs2ndlead, weight);
+ _h_yseclead[2] ->fill(yabs2ndlead, weight);
+
+ // Dijet histos
+ const double deltaphi = fabs(deltaPhi(jets[1], jets[0]));
+ const double deltarap = fabs(jets[0].momentum().rapidity() - jets[1].momentum().rapidity()) ;
+ const double deltar = fabs(deltaR(jets[0], jets[1], RAPIDITY));
+ const double mass = (jets[0].momentum() + jets[1].momentum()).mass();
+ _h_mass[chn] ->fill(mass, weight);
+ _h_mass[2] ->fill(mass, weight);
+ _h_deltay[chn] ->fill(deltarap, weight);
+ _h_deltay[2] ->fill(deltarap, weight);
+ _h_deltaphi[chn]->fill(deltaphi, weight);
+ _h_deltaphi[2] ->fill(deltaphi, weight);
+ _h_deltaR[chn] ->fill(deltar, weight);
+ _h_deltaR[2] ->fill(deltar, weight);
}
- // The leading jet histos
- if (jets.size()>=1) {
- double ptlead = jets[0].momentum().pT()/GeV;
- double yabslead = fabs(jets[0].momentum().rapidity());
- _h_ptlead[chn] ->fill(ptlead, weight);
- _h_ptlead[2] ->fill(ptlead, weight);
+ }
- _h_ylead[chn] ->fill(yabslead, weight);
- _h_ylead[2] ->fill(yabslead, weight);
- }
- if (jets.size()>=2) {
- // The second to leading jet histos
- double pt2ndlead = jets[1].momentum().pT()/GeV;
- double yabs2ndlead = fabs(jets[1].momentum().rapidity());
-
- _h_ptseclead[chn] ->fill(pt2ndlead, weight);
- _h_ptseclead[2] ->fill(pt2ndlead, weight);
- _h_yseclead[chn] ->fill(yabs2ndlead, weight);
- _h_yseclead[2] ->fill(yabs2ndlead, weight);
+ /// @name Ratio calculator util functions
+ //@{
- // Dijet histos
- double deltaphi = fabs(deltaPhi(jets[1], jets[0]));
- double deltarap = fabs(jets[0].momentum().rapidity() - jets[1].momentum().rapidity()) ;
- double deltar = fabs(deltaR(jets[0], jets[1], RAPIDITY));
- double mass = (jets[0].momentum() + jets[1].momentum()).mass();
-
- _h_mass[chn] ->fill(mass, weight);
- _h_mass[2] ->fill(mass, weight);
- _h_deltay[chn] ->fill(deltarap, weight);
- _h_deltay[2] ->fill(deltarap, weight);
- _h_deltaphi[chn] ->fill(deltaphi, weight);
- _h_deltaphi[2] ->fill(deltaphi, weight);
- _h_deltaR[chn] ->fill(deltar, weight);
- _h_deltaR[2] ->fill(deltar, weight);
- }
+ /// Calculate the ratio, being careful about div-by-zero
+ double ratio(double a, double b) {
+ return (b != 0) ? a/b : 0;
}
- /// Normalise histograms etc., after the run
+ /// Calculate the ratio error, being careful about div-by-zero
+ double ratio_err(double a, double b) {
+ return (b != 0) ? sqrt(a/sqr(b) + sqr(a)/(b*b*b)) : 0;
+ }
+ /// Calculate combined ratio from muon and electron channels
+ double comb_ratio(double* as, double* bs) {
+ return ratio(as[0]+as[1], bs[0]+bs[1]);
+ }
- std::vector<double> ratio(double a, double b) {
- double ratio=0.0;
- double ratio_err=0.0;
- std::vector<double> temp;
- cout << "a: " << a << " b: " << b << endl;
- if (b>0. && a>0.) {
- ratio = a/b;
- ratio_err = sqrt(a/pow(b,2) + pow(a,2)*b/pow(b,4));
- }
- temp.push_back(ratio);
- temp.push_back(ratio_err);
- return temp;
+ /// Calculate combined ratio error from muon and electron channels
+ double comb_ratio_err(double* as, double* bs) {
+ return ratio_err(as[0]+as[1], bs[0]+bs[1]);
}
- void finalize() {
+ //@}
- // Fill RATIO histograms (DataPointSets)
- vector<double> yvals_ratio[3];
- vector<double> yerrs_ratio[3];
- for (size_t chn = 0; chn < 2; ++chn) {
- yvals_ratio[chn].push_back(ratio(weights_nj1[chn], weights_nj0[chn])[0]);
- yvals_ratio[chn].push_back(ratio(weights_nj2[chn], weights_nj1[chn])[0]);
- yvals_ratio[chn].push_back(ratio(weights_nj3[chn], weights_nj2[chn])[0]);
- yvals_ratio[chn].push_back(ratio(weights_nj4[chn], weights_nj3[chn])[0]);
- // Errors
- yerrs_ratio[chn].push_back(ratio(weights_nj1[chn], weights_nj0[chn])[1]);
- yerrs_ratio[chn].push_back(ratio(weights_nj2[chn], weights_nj1[chn])[1]);
- yerrs_ratio[chn].push_back(ratio(weights_nj3[chn], weights_nj2[chn])[1]);
- yerrs_ratio[chn].push_back(ratio(weights_nj4[chn], weights_nj3[chn])[1]);
- // Actually fill histo
- /// \todo YODA
- //_h_njet_ratio[chn] ->setCoordinate(1, yvals_ratio[chn], yerrs_ratio[chn]);
- // Combined histos
- yvals_ratio[2].push_back(ratio(weights_nj1[chn], weights_nj0[chn])[0]);
- yvals_ratio[2].push_back(ratio(weights_nj2[chn], weights_nj1[chn])[0]);
- yvals_ratio[2].push_back(ratio(weights_nj3[chn], weights_nj2[chn])[0]);
- yvals_ratio[2].push_back(ratio(weights_nj4[chn], weights_nj3[chn])[0]);
- // Errors
- yerrs_ratio[2].push_back(ratio(weights_nj1[chn], weights_nj0[chn])[1]);
- yerrs_ratio[2].push_back(ratio(weights_nj2[chn], weights_nj1[chn])[1]);
- yerrs_ratio[2].push_back(ratio(weights_nj3[chn], weights_nj2[chn])[1]);
- yerrs_ratio[2].push_back(ratio(weights_nj4[chn], weights_nj3[chn])[1]);
- }
- /// \todo YODA
- //_h_njet_ratio[2] ->setCoordinate(1, yvals_ratio[2], yerrs_ratio[2]);
+ void finalize() {
+ // Fill ratio histograms
+ for (size_t chn = 0; chn < 2; ++chn) {
+ _h_njet_ratio[chn]->point(0).setY( ratio(weights_nj1[chn], weights_nj0[chn]), ratio_err(weights_nj1[chn], weights_nj0[chn]) );
+ _h_njet_ratio[chn]->point(1).setY( ratio(weights_nj2[chn], weights_nj1[chn]), ratio_err(weights_nj2[chn], weights_nj1[chn]) );
+ _h_njet_ratio[chn]->point(2).setY( ratio(weights_nj3[chn], weights_nj2[chn]), ratio_err(weights_nj3[chn], weights_nj2[chn]) );
+ _h_njet_ratio[chn]->point(3).setY( ratio(weights_nj4[chn], weights_nj3[chn]), ratio_err(weights_nj4[chn], weights_nj3[chn]) );
+ }
+ _h_njet_ratio[2]->point(0).setY( comb_ratio(weights_nj1, weights_nj0), comb_ratio_err(weights_nj1, weights_nj0) );
+ _h_njet_ratio[2]->point(1).setY( comb_ratio(weights_nj2, weights_nj1), comb_ratio_err(weights_nj2, weights_nj1) );
+ _h_njet_ratio[2]->point(2).setY( comb_ratio(weights_nj3, weights_nj2), comb_ratio_err(weights_nj3, weights_nj2) );
+ _h_njet_ratio[2]->point(3).setY( comb_ratio(weights_nj4, weights_nj3), comb_ratio_err(weights_nj4, weights_nj3) );
+ // Scale other histos
const double xs = crossSectionPerEvent()/picobarn;
for (size_t chn = 0; chn < 3; ++chn) {
scale(_h_njet_incl[chn], xs);
@@ -268,6 +256,7 @@
scale(_h_deltaR[chn] , xs);
scale(_h_mass[chn] , xs);
}
+
}
//@}
@@ -281,7 +270,6 @@
double weights_nj3[3];
double weights_nj4[3];
- //
Scatter2DPtr _h_njet_ratio[3];
Histo1DPtr _h_njet_incl[3];
Histo1DPtr _h_ptjet[3];
@@ -297,6 +285,8 @@
};
- // This global object acts as a hook for the plugin system
+
DECLARE_RIVET_PLUGIN(ATLAS_2011_I945498);
+
+
}
Modified: trunk/src/Analyses/ATLAS_2012_I1083318.cc
==============================================================================
--- trunk/src/Analyses/ATLAS_2012_I1083318.cc Mon May 13 21:18:12 2013 (r4272)
+++ trunk/src/Analyses/ATLAS_2012_I1083318.cc Mon May 13 21:19:10 2013 (r4273)
@@ -38,7 +38,7 @@
IdentifiedFinalState allleptons;
allleptons.acceptIdPair(PID::ELECTRON);
allleptons.acceptIdPair(PID::MUON);
- std::vector<std::pair<double, double> > etaRanges;
+ vector< pair<double, double> > etaRanges;
etaRanges.push_back(make_pair(-2.5, 2.5));
LeptonClusters leptons(fs, allleptons,
0.1, true,
@@ -61,7 +61,7 @@
jets.useInvisibles(true);
addProjection(jets, "jets");
- for (size_t i=0; i<2; ++i) {
+ for (size_t i = 0; i < 2; ++i) {
_h_NjetIncl[i] = bookHisto1D(1, 1, i+1);
_h_RatioNjetIncl[i] = bookScatter2D(2, 1, i+1);
_h_FirstJetPt_1jet[i] = bookHisto1D(3, 1, i+1);
@@ -98,30 +98,29 @@
const vector<ClusteredLepton>& leptons = applyProjection<LeptonClusters>(event, "leptons").clusteredLeptons();
Particles neutrinos = applyProjection<FinalState>(event, "neutrinos").particlesByPt();
- if (leptons.size()!=1 || (neutrinos.size()==0)) {
+ if (leptons.size() != 1 || (neutrinos.size() == 0)) {
vetoEvent;
}
- FourMomentum lepton=leptons[0].momentum();
+ FourMomentum lepton = leptons[0].momentum();
FourMomentum p_miss = neutrinos[0].momentum();
- if (p_miss.Et()<25.0*GeV) {
+ if (p_miss.Et() < 25.0*GeV) {
vetoEvent;
}
- double mT=sqrt(2.0*lepton.pT()*p_miss.Et()*(1.0-cos(lepton.phi()-p_miss.phi())));
- if (mT<40.0*GeV) {
+ double mT = sqrt(2.0 * lepton.pT() * p_miss.Et() * (1.0 - cos( lepton.phi()-p_miss.phi()) ) );
+ if (mT < 40.0*GeV) {
vetoEvent;
}
- double jetcuts[] = {30.0*GeV, 20.0*GeV};
+ double jetcuts[] = { 30.0*GeV, 20.0*GeV };
const FastJets& jetpro = applyProjection<FastJets>(event, "jets");
- for (size_t i=0; i<2; ++i) {
+ for (size_t i = 0; i < 2; ++i) {
vector<FourMomentum> jets;
- double HT=lepton.pT()+p_miss.pT();
+ double HT = lepton.pT() + p_miss.pT();
foreach (const Jet& jet, jetpro.jetsByPt(jetcuts[i])) {
- if (fabs(jet.momentum().rapidity())<4.4 &&
- deltaR(lepton, jet.momentum())>0.5) {
+ if (fabs(jet.momentum().rapidity()) < 4.4 && deltaR(lepton, jet.momentum()) > 0.5) {
jets.push_back(jet.momentum());
HT += jet.momentum().pT();
}
@@ -130,7 +129,7 @@
_h_NjetIncl[i]->fill(0.0, weight);
// Njet>=1 observables
- if (jets.size()<1) continue;
+ if (jets.size() < 1) continue;
_h_NjetIncl[i]->fill(1.0, weight);
_h_FirstJetPt_1jet[i]->fill(jets[0].pT(), weight);
_h_JetRapidity[i]->fill(jets[0].rapidity(), weight);
@@ -139,7 +138,7 @@
_h_SumYElecJet[i]->fill(lepton.rapidity()+jets[0].rapidity(), weight);
// Njet>=2 observables
- if (jets.size()<2) continue;
+ if (jets.size() < 2) continue;
_h_NjetIncl[i]->fill(2.0, weight);
_h_FirstJetPt_2jet[i]->fill(jets[0].pT(), weight);
_h_SecondJetPt_2jet[i]->fill(jets[1].pT(), weight);
@@ -151,7 +150,7 @@
_h_DeltaPhi_2jet[i]->fill(deltaPhi(jets[0], jets[1]), weight);
// Njet>=3 observables
- if (jets.size()<3) continue;
+ if (jets.size() < 3) continue;
_h_NjetIncl[i]->fill(3.0, weight);
_h_FirstJetPt_3jet[i]->fill(jets[0].pT(), weight);
_h_SecondJetPt_3jet[i]->fill(jets[1].pT(), weight);
@@ -161,7 +160,7 @@
_h_Minv_3jet[i]->fill(m2_3jet>0.0 ? sqrt(m2_3jet) : 0.0, weight);
// Njet>=4 observables
- if (jets.size()<4) continue;
+ if (jets.size() < 4) continue;
_h_NjetIncl[i]->fill(4.0, weight);
_h_FirstJetPt_4jet[i]->fill(jets[0].pT(), weight);
_h_SecondJetPt_4jet[i]->fill(jets[1].pT(), weight);
@@ -172,7 +171,7 @@
_h_Minv_4jet[i]->fill(m2_4jet>0.0 ? sqrt(m2_4jet) : 0.0, weight);
// Njet>=5 observables
- if (jets.size()<5) continue;
+ if (jets.size() < 5) continue;
_h_NjetIncl[i]->fill(5.0, weight);
}
}
@@ -180,24 +179,19 @@
/// Normalise histograms etc., after the run
void finalize() {
- for (size_t i=0; i<2; ++i) {
- /// @todo YODA
- //// first construct jet multi ratio
- //int Nbins = _h_NjetIncl[i]->numBins();
- //std::vector<double> ratio(Nbins-1, 0.0);
- //std::vector<double> err(Nbins-1, 0.0);
- //for (int n = 0; n < Nbins-1; ++n) {
- // if (_h_NjetIncl[i]->binHeight(n) > 0.0 && _h_NjetIncl[i]->binHeight(n+1) > 0.0) {
- // ratio[n] = _h_NjetIncl[i]->binHeight(n+1)/_h_NjetIncl[i]->binHeight(n);
- // double relerr_n = _h_NjetIncl[i]->binError(n)/_h_NjetIncl[i]->binHeight(n);
- // double relerr_m = _h_NjetIncl[i]->binError(n+1)/_h_NjetIncl[i]->binHeight(n+1);
- // err[n] = ratio[n] * (relerr_n + relerr_m);
- // }
- //}
- //_h_RatioNjetIncl[i]->setCoordinate(1, ratio, err);
+ for (size_t i = 0; i < 2; ++i) {
- // scale all histos to the cross section
- double factor = crossSection()/sumOfWeights();
+ // Construct jet multiplicity ratio
+ for (int n = 0; n < _h_NjetIncl[i]->numPoints()-1; ++n) {
+ YODA::Bin& b0 = _h_NjetIncl[i]->bin(n);
+ YODA::Bin& b1 = _h_NjetIncl[i]->bin(n+1);
+ if (b0.height() == 0.0 || b1.height() == 0.0) continue;
+ _h_RatioNjetIncl[i]->point(n).setY(b1.height()/b0.height());
+ _h_RatioNjetIncl[i]->point(n).setYErr(b1.height()/b0.height() * (b0.relErr() + b1.relErr()));
+ }
+
+ // Scale all histos to the cross section
+ const double factor = crossSection()/sumOfWeights();
scale(_h_DeltaPhi_2jet[i], factor);
scale(_h_DeltaR_2jet[i], factor);
scale(_h_DeltaY_2jet[i], factor);
More information about the Rivet-svn
mailing list
|