|
[Rivet-svn] r2441 - in trunk: . include/Rivet/Projections plugindemo src/Analyses src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgTue May 11 02:29:19 BST 2010
Author: buckley Date: Tue May 11 11:58:03 2010 New Revision: 2441 Log: Replacing uses of TotalVisibleMomentum with MissingMomentum, with corresponding improvements of MM interface, plus some small code tidying and highlighting of bad/apparently unused code in CDF 1994 and D0 2004 analyses. Modified: trunk/ChangeLog trunk/include/Rivet/Projections/MissingMomentum.hh trunk/include/Rivet/Projections/TotalVisibleMomentum.hh trunk/plugindemo/RootAnalysis.cc trunk/src/Analyses/CDF_1994_S2952106.cc trunk/src/Analyses/CDF_2005_S6217184.cc trunk/src/Analyses/D0_2004_S5992206.cc trunk/src/Analyses/MC_SUSY.cc trunk/src/Projections/WFinder.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Fri May 7 21:11:18 2010 (r2440) +++ trunk/ChangeLog Tue May 11 11:58:03 2010 (r2441) @@ -1,3 +1,8 @@ +2010-05-11 Andy Buckley <andy at insectnation.org> + + * Replacing TotalVisibleMomentum with MissingMomentum in analyses + and WFinder. Using vector ET rather than scalar ET in some places. + 2010-05-07 Andy Buckley <andy at insectnation.org> * Revamping the AnalysisHandler constructor and data writing, with Modified: trunk/include/Rivet/Projections/MissingMomentum.hh ============================================================================== --- trunk/include/Rivet/Projections/MissingMomentum.hh Fri May 7 21:11:18 2010 (r2440) +++ trunk/include/Rivet/Projections/MissingMomentum.hh Tue May 11 11:58:03 2010 (r2441) @@ -46,13 +46,16 @@ public: - /// The projected four-momentum vector - FourMomentum& momentum() { return _momentum; } + /// The vector-summed visible four-momentum in the event. + FourMomentum& visibleMomentum() { return _momentum; } - /// The projected four-momentum vector - const FourMomentum& momentum() const { return _momentum; } + /// The vector-summed visible four-momentum in the event. + const FourMomentum& visibleMomentum() const { return _momentum; } - /// The projected scalar transverse energy + /// The vector-summed (in)visible transverse energy in the event + double vectorET() const { return _momentum.Et(); } + + /// The scalar-summed (in)visible transverse energy in the event. double scalarET() const { return _set; } Modified: trunk/include/Rivet/Projections/TotalVisibleMomentum.hh ============================================================================== --- trunk/include/Rivet/Projections/TotalVisibleMomentum.hh Fri May 7 21:11:18 2010 (r2440) +++ trunk/include/Rivet/Projections/TotalVisibleMomentum.hh Tue May 11 11:58:03 2010 (r2441) @@ -22,6 +22,8 @@ { setName("TotalVisibleMomentum"); addProjection(fsp, "FS"); + getLog() << Log::WARNING << "TotalVisibleMomentum projection is deprecated: " + << "please use the MissingMomentum projection instead." << endl; } /// Clone on the heap. Modified: trunk/plugindemo/RootAnalysis.cc ============================================================================== --- trunk/plugindemo/RootAnalysis.cc Fri May 7 21:11:18 2010 (r2440) +++ trunk/plugindemo/RootAnalysis.cc Tue May 11 11:58:03 2010 (r2441) @@ -2,7 +2,7 @@ #include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/ChargedLeptons.hh" -#include "Rivet/Projections/TotalVisibleMomentum.hh" +#include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/FastJets.hh" // ROOT stuff @@ -16,13 +16,13 @@ /// @brief Book and fill a ROOT tree with simulated data. /// /// This does some things, e.g. access parton level information, which - /// are not recommended in Rivet analyses, since the information is + /// are not recommended in Rivet analyses, since the information is /// unphysical and so cannot be compared to data, and also may be generator dependent. - /// + /// class RootAnalysis : public Analysis { public: - RootAnalysis() : Analysis("ROOTANALYSIS") { + RootAnalysis() : Analysis("ROOTANALYSIS") { // Choose cuts _jet_pt_cut = 20*GeV; _subj_pt_cut = 20*GeV; @@ -30,24 +30,15 @@ _store_partons = true; _treeFileName = "rivetTree.root"; } - - + + void init() { const FinalState fs(-4.0, 4.0, 0.0*GeV); addProjection(fs, "FS"); addProjection(ChargedLeptons(fs), "ChLeptons"); addProjection(FastJets(fs, FastJets::KT, 0.7), "Jets"); - - /// Veto neutrinos, antineutrinos and LSP - VetoedFinalState vfs(fs); - vfs - .addVetoDetail(NU_E, 10.0*GeV, 50.0*GeV) - .addVetoPairId(NU_MU) - .addVetoPairId(NU_TAU) - .addVetoId(1000022); // Assumes that neutralino_1 is the LSP - addProjection(vfs, "VFS"); - addProjection(TotalVisibleMomentum(vfs), "TotalVisMom"); - + addProjection(MissingMomentum(fs), "TotalVisMom"); + ZFinder zs(fs, ELECTRON, 80*GeV, 100*GeV, 0.2); addProjection(zs, "Zs"); @@ -55,12 +46,12 @@ _treeFile = new TFile(_treeFileName, "recreate"); _rivetTree = new TTree("Rivet Tree", "Rivet Example Tree"); _rivetTree->Branch("nevt", &_nevt, "nevt/I"); - + // Vector bosons _rivetTree->Branch("nvb", &_nvb, "nvb/I"); _rivetTree->Branch("vbtype", &_vbtype, "vbtype[nvb]/I"); _rivetTree->Branch("vbvec", &_vbvec, "vbvec[nvb][4]/F"); - // Jets + // Jets _rivetTree->Branch("njet", &_njet, "njet/I"); _rivetTree->Branch("vjet", &_vjet, "vjet[njet][4]/F"); // Subjets @@ -79,13 +70,13 @@ // Missing Et _rivetTree->Branch("esumr", &_esumr, "esumr[4]/F"); } - + // Do the analysis void analyze(const Event& event) { const GenEvent& ev = event.genEvent(); _nevt = ev.event_number(); - + // Get the vector bosons _nvb = 0; const FinalState& zs = applyProjection<FinalState>(event, "Zs"); @@ -98,13 +89,13 @@ _vbtype[_nvb] = 1; ++_nvb; } - + // Get the partons. This is generator-dependent and should not be // used in normal analyses. _npart = 0; if (_store_partons) { foreach (const HepMC::GenParticle* p, particles(event.genEvent())) { - // Only include particles which are documentation line (status >1) + // Only include particles which are documentation line (status >1) // The result/meaning will be generator dependent. if (p->status() >= 2) { const FourMomentum p4 = p->momentum(); @@ -130,8 +121,8 @@ } } } - - + + // Get the jets in decreasing pT order. const FastJets& jets = applyProjection<FastJets>(event, "Jets"); PseudoJets jetList = jets.pseudoJetsByPt(); @@ -156,12 +147,12 @@ _ysubsj[_nsub][i] = 0; } } - ++_nsub; + ++_nsub; } ++_njet; } } - + // Loop over leptons _nlep = 0; const ChargedLeptons& cl = applyProjection<ChargedLeptons>(event, "ChLeptons"); @@ -175,24 +166,24 @@ ++_nlep; } } - - // Missing Et/total energy - const TotalVisibleMomentum& tvm = applyProjection<TotalVisibleMomentum>(event, "TotalVisMom"); - _esumr[0] = tvm.momentum().E()/GeV; - _esumr[1] = tvm.momentum().px()/GeV; - _esumr[2] = tvm.momentum().py()/GeV; - _esumr[3] = tvm.momentum().pz()/GeV; - + + // Total energy vector == -1 * invisible momentum + const MissingMomentum& tvm = applyProjection<MissingMomentum>(event, "TotalVisMom"); + _esumr[0] = tvm.visibleMomentum().E()/GeV; + _esumr[1] = tvm.visibleMomentum().px()/GeV; + _esumr[2] = tvm.visibleMomentum().py()/GeV; + _esumr[3] = tvm.visibleMomentum().pz()/GeV; + // Finally fill the tree _rivetTree->Fill(); } - - - void finalize() { + + + void finalize() { // Write the tree to file. _rivetTree->Write(); } - + //@} @@ -200,7 +191,7 @@ /// The tree TTree* _rivetTree; - + /// The file for the Tree TFile* _treeFile; @@ -211,22 +202,22 @@ /// @name The ntuple variables. //@{ /// Event number - int _nevt; + int _nevt; /// Number of W bosons - int _nvb; + int _nvb; /// 4 momentum of W bosons. float _vbvec[8][4]; /// Type (i.e. decay mode) of W bosons. - int _vbtype[8]; + int _vbtype[8]; /// Number of jets - int _njet; + int _njet; /// Four momentum of the jets - float _vjet[50][4]; + float _vjet[50][4]; /// Number of jets for which the subjet analysis was performed. - int _nsub; + int _nsub; /// Four vector of jets for which we found subjets. float _sjet3[200][4]; /// y 1->2, 2->3, 3->4, 4->5 for the above jets. @@ -239,7 +230,7 @@ float _vlep[150][4]; /// Number of partons - int _npart; + int _npart; float _ppart[4000][4]; int _pid[4000]; int _mo[4000]; @@ -262,7 +253,7 @@ }; - + // This global object acts as a hook for the plugin system AnalysisBuilder<RootAnalysis> plugin_RootAnalysis; Modified: trunk/src/Analyses/CDF_1994_S2952106.cc ============================================================================== --- trunk/src/Analyses/CDF_1994_S2952106.cc Fri May 7 21:11:18 2010 (r2440) +++ trunk/src/Analyses/CDF_1994_S2952106.cc Tue May 11 11:58:03 2010 (r2441) @@ -2,9 +2,10 @@ #include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" -#include "Rivet/Projections/VetoedFinalState.hh" -#include "Rivet/Projections/TotalVisibleMomentum.hh" #include "Rivet/Projections/FastJets.hh" +#include "Rivet/Projections/VetoedFinalState.hh" +#include "Rivet/Projections/VisibleFinalState.hh" +#include "Rivet/Projections/MissingMomentum.hh" namespace Rivet { @@ -30,11 +31,11 @@ const FinalState fs(-4.2, 4.2); addProjection(fs, "FS"); addProjection(FastJets(fs, FastJets::CDFJETCLU, 0.7), "ConeJets"); - addProjection(TotalVisibleMomentum(fs), "CalMET"); + addProjection(MissingMomentum(fs), "CalMET"); // Veto (anti)neutrinos, and muons with pT above 1.0 GeV - VetoedFinalState vfs(fs); - vfs.vetoNeutrinos(); + /// @todo This doesn't seem to be being used for anything... + VetoedFinalState vfs(VisibleFinalState(fs)); vfs.addVetoPairDetail(MUON, 1.0*GeV, MAXDOUBLE); addProjection(vfs, "VFS"); @@ -79,7 +80,7 @@ /// @todo Need better title //_histR23Ideal = bookHistogram1D("R23Ideal", 50, 0., 5.); _histR23Ideal = bookHistogram1D(7,1,2); - + /// @todo Need better title //_histJet3etaCDF = bookHistogram1D("Jet3etaCDF", 42, -4., 4.); @@ -88,20 +89,20 @@ /// @todo Need better title //_histJet3etaIdeal = bookHistogram1D("Jet3etaIdeal", 42, -4., 4.); _histJet3etaIdeal = bookHistogram1D(8,1,2); - + //_histEta3Corr = bookHistogram1D(3,1,1); - const string hnameEta3 = "Eta3Corr"; + const string hnameEta3 = "Eta3Corr"; //const string htitleEta3 = "Eta3Corr"; _histEta3Corr = bookHistogram1D(hnameEta3, 40, -4., 4.); //_histR23Corr = bookHistogram1D(4,1,1); - const string hnameR23 = "R23Corr"; + const string hnameR23 = "R23Corr"; //const string htitleR23 = "R23Corr"; _histR23Corr = bookHistogram1D(hnameR23, 35, 0., 4.375); //_histAlphaCorr = bookHistogram1D(5,1,1); - const string hnameAlpha = "AlphaCorr"; + const string hnameAlpha = "AlphaCorr"; //const string htitleAlpha = "AlphaCorr"; _histAlphaCorr = bookHistogram1D(hnameAlpha, 40, -90., 90.); @@ -115,10 +116,8 @@ getLog() << Log::DEBUG << "Jet multiplicity before any cuts = " << jets.size() << endl; // Check there isn't too much missing Et - const TotalVisibleMomentum& caloMissEt = applyProjection<TotalVisibleMomentum>(event, "CalMET"); + const MissingMomentum& caloMissEt = applyProjection<MissingMomentum>(event, "CalMET"); getLog() << Log::DEBUG << "Missing pT = " << caloMissEt.momentum().pT()/GeV << " GeV" << endl; - /// @todo should this really be scalar ET here, and not caloMissEt.momentum().Et()? - //if ((caloMissEt.momentum().pT()/GeV)/sqrt(caloMissEt.scalarET()/GeV) > 6.0) vetoEvent; if ((caloMissEt.momentum().pT()/GeV)/sqrt(caloMissEt.momentum().Et()/GeV) > 6.0) vetoEvent; // Check jet requirements @@ -130,8 +129,8 @@ if (fabs(pj1.eta()) > 0.7 || fabs(pj2.eta()) > 0.7) vetoEvent; getLog() << Log::DEBUG << "Jet 1 & 2 eta, pT requirements fulfilled" << endl; - // back-to-bak within 20 degrees in phi - if (deltaPhi(pj1.phi(), pj2.phi()) < PI-PI/18.) vetoEvent; + // Require that jets are back-to-back within 20 degrees in phi + if (deltaPhi(pj1.phi(), pj2.phi()) < 17.0/18.0 * PI) vetoEvent; getLog() << Log::DEBUG << "Jet 1 & 2 phi requirement fulfilled" << endl; const double weight = event.weight(); @@ -154,127 +153,131 @@ /// Finalize void finalize() { - - const double eta3_CDF_sim[] = {0.0013, 0.0037, 0.0047, 0.0071, 0.0093, - 0.0117, 0.0151, 0.0149, 0.0197, 0.0257, - 0.0344, 0.0409, 0.0481, 0.0454, 0.0394, - 0.0409, 0.0387, 0.0387, 0.0322, 0.0313, - 0.029, 0.0309, 0.0412, 0.0417, 0.0412, - 0.0397, 0.0417, 0.0414, 0.0376, 0.0316, - 0.027, 0.0186, 0.0186, 0.0132, 0.0127, + + const double eta3_CDF_sim[] = {0.0013, 0.0037, 0.0047, 0.0071, 0.0093, + 0.0117, 0.0151, 0.0149, 0.0197, 0.0257, + 0.0344, 0.0409, 0.0481, 0.0454, 0.0394, + 0.0409, 0.0387, 0.0387, 0.0322, 0.0313, + 0.029, 0.0309, 0.0412, 0.0417, 0.0412, + 0.0397, 0.0417, 0.0414, 0.0376, 0.0316, + 0.027, 0.0186, 0.0186, 0.0132, 0.0127, 0.0106, 0.0071, 0.004, 0.002, 0.0013}; - - const double eta3_CDF_sim_err[] = {0.0009, 0.0009, 0.0007, 0.0007, 0.0007, - 0.001, 0.0012, 0.0012, 0.0013, 0.0016, - 0.0017, 0.002, 0.002, 0.0022, 0.002, - 0.002, 0.0018, 0.0018, 0.0016, 0.0017, - 0.0017, 0.0019, 0.002, 0.0021, 0.002, - 0.002, 0.0019, 0.002, 0.0018, 0.0017, - 0.0017, 0.0014, 0.0014, 0.0009, 0.001, + + const double eta3_CDF_sim_err[] = {0.0009, 0.0009, 0.0007, 0.0007, 0.0007, + 0.001, 0.0012, 0.0012, 0.0013, 0.0016, + 0.0017, 0.002, 0.002, 0.0022, 0.002, + 0.002, 0.0018, 0.0018, 0.0016, 0.0017, + 0.0017, 0.0019, 0.002, 0.0021, 0.002, + 0.002, 0.0019, 0.002, 0.0018, 0.0017, + 0.0017, 0.0014, 0.0014, 0.0009, 0.001, 0.0009, 0.0009, 0.0008, 0.0008, 0.0009}; - - const double eta3_Ideal_sim[] = {0.0017, 0.003, 0.0033, 0.0062, 0.0062, - 0.0112, 0.0177, 0.0164, 0.0196, 0.0274, - 0.0351, 0.0413, 0.052, 0.0497, 0.0448, - 0.0446, 0.0375, 0.0329, 0.0291, 0.0272, - 0.0233, 0.0288, 0.0384, 0.0396, 0.0468, - 0.0419, 0.0459, 0.0399, 0.0355, 0.0329, - 0.0274, 0.023, 0.0201, 0.012, 0.01, + + const double eta3_Ideal_sim[] = {0.0017, 0.003, 0.0033, 0.0062, 0.0062, + 0.0112, 0.0177, 0.0164, 0.0196, 0.0274, + 0.0351, 0.0413, 0.052, 0.0497, 0.0448, + 0.0446, 0.0375, 0.0329, 0.0291, 0.0272, + 0.0233, 0.0288, 0.0384, 0.0396, 0.0468, + 0.0419, 0.0459, 0.0399, 0.0355, 0.0329, + 0.0274, 0.023, 0.0201, 0.012, 0.01, 0.008, 0.0051, 0.0051, 0.001, 0.001}; - - for (int ibin = 0; ibin < 40; ++ibin) //x-value + weight=err^2/y-value - _histEta3Corr->fill(-4. + (ibin+0.5)*8./40., //x-value - eta3_CDF_sim_err[ibin]/eta3_Ideal_sim[ibin] //error - / eta3_CDF_sim[ibin] * eta3_Ideal_sim[ibin] //norm. weight to value - ); //(Ideal error identical zero) - - - - const double R23_CDF_sim[] = {0.0005, 0.0161, 0.057, 0.0762, 0.0723, - 0.0705, 0.0598, 0.0563, 0.0557, 0.0579, - 0.0538, 0.0522, 0.0486, 0.0449, 0.0418, - 0.0361, 0.0326, 0.0304, 0.0252, 0.0212, - 0.0173, 0.0176, 0.0145, 0.0127, 0.0103, - 0.0065, 0.0049, 0.0045, 0.0035, 0.0029, + + for (int ibin = 0; ibin < 40; ++ibin) { + // x-value + weight = err^2/y-value + _histEta3Corr->fill(-4. + (ibin+0.5)*8./40., //x-value + eta3_CDF_sim_err[ibin]/eta3_Ideal_sim[ibin] //error + / eta3_CDF_sim[ibin] * eta3_Ideal_sim[ibin] //norm. weight to value + ); //(Ideal error identical zero) + } + + const double R23_CDF_sim[] = {0.0005, 0.0161, 0.057, 0.0762, 0.0723, + 0.0705, 0.0598, 0.0563, 0.0557, 0.0579, + 0.0538, 0.0522, 0.0486, 0.0449, 0.0418, + 0.0361, 0.0326, 0.0304, 0.0252, 0.0212, + 0.0173, 0.0176, 0.0145, 0.0127, 0.0103, + 0.0065, 0.0049, 0.0045, 0.0035, 0.0029, 0.0024, 0.0014, 0.0011, 0.001, 0.0009}; - - const double R23_CDF_sim_err[] = {0.0013, 0.0009, 0.0022, 0.0029, 0.0026, - 0.0024, 0.0022, 0.0025, 0.0023, 0.0024, - 0.0021, 0.0021, 0.0021, 0.0021, 0.0021, - 0.0019, 0.0019, 0.0016, 0.0017, 0.0014, - 0.001, 0.0014, 0.0012, 0.0013, 0.001, - 0.0011, 0.001, 0.001, 0.001, 0.0011, + + const double R23_CDF_sim_err[] = {0.0013, 0.0009, 0.0022, 0.0029, 0.0026, + 0.0024, 0.0022, 0.0025, 0.0023, 0.0024, + 0.0021, 0.0021, 0.0021, 0.0021, 0.0021, + 0.0019, 0.0019, 0.0016, 0.0017, 0.0014, + 0.001, 0.0014, 0.0012, 0.0013, 0.001, + 0.0011, 0.001, 0.001, 0.001, 0.0011, 0.0011, 0.0009, 0.0008, 0.0008, 0.0009}; - - const double R23_Ideal_sim[] = {0.0005, 0.0176, 0.0585, 0.0862, 0.0843, - 0.0756, 0.0673, 0.0635, 0.0586, 0.0619, - 0.0565, 0.0515, 0.0466, 0.0472, 0.0349, - 0.0349, 0.0266, 0.0254, 0.0204, 0.0179, - 0.0142, 0.0134, 0.0101, 0.009, 0.008, - 0.0034, 0.003, 0.0033, 0.0027, 0.0021, + + const double R23_Ideal_sim[] = {0.0005, 0.0176, 0.0585, 0.0862, 0.0843, + 0.0756, 0.0673, 0.0635, 0.0586, 0.0619, + 0.0565, 0.0515, 0.0466, 0.0472, 0.0349, + 0.0349, 0.0266, 0.0254, 0.0204, 0.0179, + 0.0142, 0.0134, 0.0101, 0.009, 0.008, + 0.0034, 0.003, 0.0033, 0.0027, 0.0021, 0.0012, 0.0006, 0.0004, 0.0005, 0.0003}; - - for (int ibin = 0; ibin < 35; ++ibin) //x-value + weight=err^2/y-value - _histR23Corr->fill((ibin+0.5)*4.375/35., //x-value - R23_CDF_sim_err[ibin]/R23_Ideal_sim[ibin] //error - / R23_CDF_sim[ibin] * R23_Ideal_sim[ibin] //norm. weight to y-value - ); //(Ideal error identical zero) - - - - - const double alpha_CDF_sim[] = {0.0517, 0.0461, 0.049, 0.0452, 0.0451, - 0.0435, 0.0317, 0.0287, 0.0294, 0.0261, - 0.0231, 0.022, 0.0233, 0.0192, 0.0213, - 0.0166, 0.0176, 0.0146, 0.0136, 0.0156, - 0.0142, 0.0152, 0.0151, 0.0147, 0.0164, - 0.0186, 0.018, 0.021, 0.0198, 0.0189, - 0.0197, 0.0211, 0.027, 0.0236, 0.0243, + + for (int ibin = 0; ibin < 35; ++ibin) { + // x-value + weight=err^2/y-value + _histR23Corr->fill((ibin+0.5)*4.375/35., //x-value + R23_CDF_sim_err[ibin]/R23_Ideal_sim[ibin] //error + / R23_CDF_sim[ibin] * R23_Ideal_sim[ibin] //norm. weight to y-value + ); //(Ideal error identical zero) + } + + + + const double alpha_CDF_sim[] = {0.0517, 0.0461, 0.049, 0.0452, 0.0451, + 0.0435, 0.0317, 0.0287, 0.0294, 0.0261, + 0.0231, 0.022, 0.0233, 0.0192, 0.0213, + 0.0166, 0.0176, 0.0146, 0.0136, 0.0156, + 0.0142, 0.0152, 0.0151, 0.0147, 0.0164, + 0.0186, 0.018, 0.021, 0.0198, 0.0189, + 0.0197, 0.0211, 0.027, 0.0236, 0.0243, 0.0269, 0.0257, 0.0276, 0.0246, 0.0286}; - - const double alpha_CDF_sim_err[] = {0.0024, 0.0025, 0.0024, 0.0024, 0.0024, - 0.0022, 0.0019, 0.0018, 0.0019, 0.0016, - 0.0017, 0.0017, 0.0019, 0.0013, 0.0017, - 0.0014, 0.0016, 0.0013, 0.0012, 0.0009, - 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, - 0.0015, 0.0014, 0.0016, 0.0016, 0.0015, - 0.0016, 0.0016, 0.0019, 0.0017, 0.0019, + + const double alpha_CDF_sim_err[] = {0.0024, 0.0025, 0.0024, 0.0024, 0.0024, + 0.0022, 0.0019, 0.0018, 0.0019, 0.0016, + 0.0017, 0.0017, 0.0019, 0.0013, 0.0017, + 0.0014, 0.0016, 0.0013, 0.0012, 0.0009, + 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, + 0.0015, 0.0014, 0.0016, 0.0016, 0.0015, + 0.0016, 0.0016, 0.0019, 0.0017, 0.0019, 0.0018, 0.0018, 0.0018, 0.0018, 0.0019}; - - const double alpha_Ideal_sim[] = {0.0552, 0.0558, 0.0583, 0.055, 0.0495, - 0.0433, 0.0393, 0.0346, 0.0331, 0.0296, - 0.0258, 0.0196, 0.0171, 0.0179, 0.0174, - 0.0141, 0.0114, 0.0096, 0.0076, 0.0087, - 0.0099, 0.0079, 0.0102, 0.0114, 0.0124, - 0.013, 0.0165, 0.016, 0.0177, 0.019, - 0.0232, 0.0243, 0.0238, 0.0248, 0.0235, + + const double alpha_Ideal_sim[] = {0.0552, 0.0558, 0.0583, 0.055, 0.0495, + 0.0433, 0.0393, 0.0346, 0.0331, 0.0296, + 0.0258, 0.0196, 0.0171, 0.0179, 0.0174, + 0.0141, 0.0114, 0.0096, 0.0076, 0.0087, + 0.0099, 0.0079, 0.0102, 0.0114, 0.0124, + 0.013, 0.0165, 0.016, 0.0177, 0.019, + 0.0232, 0.0243, 0.0238, 0.0248, 0.0235, 0.0298, 0.0292, 0.0291, 0.0268, 0.0316}; - - for (int ibin = 0; ibin < 40; ++ibin) {//x-value + weight=err^2/y-value - _histAlphaCorr->fill(-90. + (ibin+0.5)*180./40., //y-value - alpha_CDF_sim_err[ibin]/alpha_Ideal_sim[ibin] //error - / alpha_CDF_sim[ibin] * alpha_Ideal_sim[ibin] //norm. weight to y-value - ); //(Ideal error identical zero) - cout << "bin i=" << ibin << " val=" << alpha_CDF_sim[ibin]/alpha_Ideal_sim[ibin] - << " err^2/value=" << alpha_CDF_sim_err[ibin]/alpha_Ideal_sim[ibin] - / alpha_CDF_sim[ibin] * alpha_Ideal_sim[ibin] << endl; + + for (int ibin = 0; ibin < 40; ++ibin) { + //x-value + weight=err^2/y-value + _histAlphaCorr->fill(-90. + (ibin+0.5)*180./40., //y-value + alpha_CDF_sim_err[ibin]/alpha_Ideal_sim[ibin] //error + / alpha_CDF_sim[ibin] * alpha_Ideal_sim[ibin] //norm. weight to y-value + ); //(Ideal error identical zero) + //getLog() << Log::TRACE << "bin i=" << ibin << " val=" << alpha_CDF_sim[ibin]/alpha_Ideal_sim[ibin] + // << " err^2/value=" << alpha_CDF_sim_err[ibin]/alpha_Ideal_sim[ibin] / alpha_CDF_sim[ibin] * alpha_Ideal_sim[ibin] << endl; } - + AIDA::IHistogramFactory& hf = histogramFactory(); - + + /// @todo Histo factory output paths don't work this way //hf.multiply(histoDir() + "/d03-x01-y01", *_histJet3eta, *_histEta3Corr); hf.multiply("/_histJet3eta", *_histJet3eta, *_histEta3Corr); //_histJet3eta = hf.multiply("/_histJet3eta", *_histJet3eta, *_histEta3Corr); hf.destroy(_histEta3Corr); - + + /// @todo Histo factory output paths don't work this way //hf.multiply(histoDir() + "/d04-x01-y01", *_histR23, *_histR23Corr); hf.multiply("/_histR23", *_histR23, *_histR23Corr); hf.destroy(_histR23Corr); + /// @todo Histo factory output paths don't work this way //hf.multiply(histoDir() + "/d05-x01-y01", *_histAlpha, *_histAlphaCorr); hf.multiply("/_histAlpha", *_histAlpha, *_histAlphaCorr); hf.destroy(_histAlphaCorr); - + //getLog() << Log::INFO << "Cross-section = " << crossSection()/picobarn << " pb" << endl; normalize(_histJet1Et, 12.3025); //norm to integral of Ref data Modified: trunk/src/Analyses/CDF_2005_S6217184.cc ============================================================================== --- trunk/src/Analyses/CDF_2005_S6217184.cc Fri May 7 21:11:18 2010 (r2440) +++ trunk/src/Analyses/CDF_2005_S6217184.cc Tue May 11 11:58:03 2010 (r2441) @@ -3,7 +3,8 @@ #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FastJets.hh" -#include "Rivet/Projections/TotalVisibleMomentum.hh" +#include "Rivet/Projections/VetoedFinalState.hh" +#include "Rivet/Projections/VisibleFinalState.hh" #include "Rivet/Projections/JetShape.hh" namespace Rivet { @@ -30,13 +31,11 @@ const FinalState fs(-2.0, 2.0); addProjection(fs, "FS"); addProjection(FastJets(fs, FastJets::CDFMIDPOINT, 0.7), "Jets"); - addProjection(TotalVisibleMomentum(fs), "CalMET"); // Veto (anti)neutrinos, and muons with pT above 1.0 GeV - VetoedFinalState vfs(fs); - vfs.vetoNeutrinos(); + VisibleFinalState visfs(fs); + VetoedFinalState vfs(visfs); vfs.addVetoPairDetail(MUON, 1.0*GeV, MAXDOUBLE); - addProjection(vfs, "VFS"); addProjection(JetShape(vfs, _jetaxes, 0.0, 0.7, 0.1, 0.3), "JetShape"); // Specify pT bins @@ -62,7 +61,7 @@ void analyze(const Event& event) { // Get jets and require at least one to pass pT and y cuts - const Jets jets = applyProjection<FastJets>(event, "Jets").jetsByPt(); + const Jets& jets = applyProjection<FastJets>(event, "Jets").jetsByPt(); getLog() << Log::DEBUG << "Jet multiplicity before cuts = " << jets.size() << endl; // Determine the central jet axes Modified: trunk/src/Analyses/D0_2004_S5992206.cc ============================================================================== --- trunk/src/Analyses/D0_2004_S5992206.cc Fri May 7 21:11:18 2010 (r2440) +++ trunk/src/Analyses/D0_2004_S5992206.cc Tue May 11 11:58:03 2010 (r2441) @@ -3,7 +3,9 @@ #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FastJets.hh" -#include "Rivet/Projections/TotalVisibleMomentum.hh" +#include "Rivet/Projections/VetoedFinalState.hh" +#include "Rivet/Projections/VisibleFinalState.hh" +#include "Rivet/Projections/MissingMomentum.hh" namespace Rivet { @@ -45,12 +47,13 @@ const FinalState fs(-3.0, 3.0); addProjection(fs, "FS"); addProjection(FastJets(FinalState(), FastJets::D0ILCONE, 0.7), "Jets"); - addProjection(TotalVisibleMomentum(fs), "CalMET"); + addProjection(MissingMomentum(fs), "CalMET"); // Veto neutrinos, and muons with pT above 1.0 GeV - VetoedFinalState vfs(fs); - vfs.vetoNeutrinos(); - vfs.addVetoPairDetail(MUON, 1.0, MAXDOUBLE); + /// @todo This doesn't seem to be used for anything: fix! + VisibleFinalState visfs(fs); + VetoedFinalState vfs(visfs); + vfs.addVetoPairDetail(MUON, 1.0*GeV, MAXDOUBLE); addProjection(vfs, "VFS"); // Book histograms @@ -82,9 +85,9 @@ getLog() << Log::DEBUG << "Jet eta and pT requirements fulfilled" << endl; const double pT1 = jets[0].momentum().pT(); - const TotalVisibleMomentum& caloMissEt = applyProjection<TotalVisibleMomentum>(event, "CalMET"); - getLog() << Log::DEBUG << "Missing Et = " << caloMissEt.momentum().pT()/GeV << endl; - if (caloMissEt.momentum().pT() > 0.7*pT1) { + const MissingMomentum& caloMissEt = applyProjection<MissingMomentum>(event, "CalMET"); + getLog() << Log::DEBUG << "Missing scalar Et = " << caloMissEt.scalarET()/GeV << " GeV" << endl; + if (caloMissEt.scalarET() > 0.7*pT1) { vetoEvent; } Modified: trunk/src/Analyses/MC_SUSY.cc ============================================================================== --- trunk/src/Analyses/MC_SUSY.cc Fri May 7 21:11:18 2010 (r2440) +++ trunk/src/Analyses/MC_SUSY.cc Tue May 11 11:58:03 2010 (r2441) @@ -7,7 +7,7 @@ #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" -#include "Rivet/Projections/TotalVisibleMomentum.hh" +#include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" namespace Rivet { @@ -48,11 +48,8 @@ mufs.acceptIdPair(MUON); addProjection(mufs, "Muons"); - VetoedFinalState missing(fs); - missing.vetoNeutrinos(); - missing.addVetoId(1000022); // lightest neutralino = usual LSP - missing.addVetoId(1000039); // gravitino = LSP in GMSB - addProjection(TotalVisibleMomentum(missing), "MET"); + MissingMomentum missing(fs); + addProjection(missing, "MET"); LeadingParticlesFinalState lpfs(fs); lpfs.addParticleIdPair(ELECTRON); @@ -206,9 +203,8 @@ } // Calculate and fill missing Et histos - const TotalVisibleMomentum& met = applyProjection<TotalVisibleMomentum>(evt, "MET"); - /// @todo should this really be the scalarET sum, and not met.momentum().Et()? - _hist_met->fill(met.scalarET()/GeV); + const MissingMomentum& met = applyProjection<MissingMomentum>(evt, "MET"); + _hist_met->fill(met.vectorET()/GeV); // Choose highest-pT leptons of each sign and flavour for dilepton mass edges const FinalState& lpfs = applyProjection<FinalState>(evt, "LeadingParticles"); Modified: trunk/src/Projections/WFinder.cc ============================================================================== --- trunk/src/Projections/WFinder.cc Fri May 7 21:11:18 2010 (r2440) +++ trunk/src/Projections/WFinder.cc Tue May 11 11:58:03 2010 (r2441) @@ -1,7 +1,7 @@ // -*- C++ -*- #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/InvMassFinalState.hh" -#include "Rivet/Projections/TotalVisibleMomentum.hh" +#include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/MergedFinalState.hh" #include "Rivet/Projections/ClusteredPhotons.hh" #include "Rivet/Projections/VetoedFinalState.hh" @@ -44,7 +44,7 @@ void WFinder::_init(const std::vector<std::pair<double, double> >& etaRanges, - double pTmin, + double pTmin, PdgId pid, double m2_min, double m2_max, double missingET, @@ -87,7 +87,7 @@ l_nu_ids += make_pair(-abs(pid), abs(nu_pid)); InvMassFinalState imfs(mergedFS, l_nu_ids, m2_min, m2_max); addProjection(imfs, "IMFS"); - + // A projection for clustering photons on to the charged lepton ClusteredPhotons cphotons(FinalState(), imfs, dRmax); addProjection(cphotons, "CPhotons"); @@ -95,11 +95,11 @@ // Projection for all signal constituents MergedFinalState signalFS(imfs, cphotons); addProjection(cphotons, "SignalParticles"); - - // Add TotalVisibleMomentum proj to calc MET - TotalVisibleMomentum vismom(signalFS); + + // Add MissingMomentum proj to calc MET + MissingMomentum vismom(signalFS); addProjection(vismom, "MissingET"); - + // FS for non-signal bits of the event VetoedFinalState remainingFS; remainingFS.addVetoOnThisFinalState(signalFS); @@ -171,10 +171,10 @@ msg << " = " << pW; // Check missing ET - const TotalVisibleMomentum& vismom = applyProjection<TotalVisibleMomentum>(e, "MissingET"); - /// @todo Restrict missing momentum eta range? + const MissingMomentum& vismom = applyProjection<MissingMomentum>(e, "MissingET"); + /// @todo Restrict missing momentum eta range? Use vectorET()? if (vismom.scalarET() < _etMiss) { - getLog() << Log::DEBUG << "Not enough missing ET: " << vismom.scalarET()/GeV + getLog() << Log::DEBUG << "Not enough missing ET: " << vismom.scalarET()/GeV << " GeV vs. " << _etMiss/GeV << " GeV" << endl; return; } @@ -182,7 +182,7 @@ // Check mass range again if (!inRange(pW.mass()/GeV, _m2_min, _m2_max)) return; getLog() << Log::DEBUG << msg.str() << endl; - + // Make W Particle and insert into particles list const PdgId wpid = (wcharge == 1) ? WPLUSBOSON : WMINUSBOSON; _theParticles.push_back(Particle(wpid, pW));
More information about the Rivet-svn mailing list |