|
[Rivet] Rivet analysis for CMS UE in DYHendrik Hoeth hendrik.hoeth at cern.chThu Oct 4 15:11:16 BST 2012
Hi Frank, Thus spake Frank Siegert (frank.siegert at cern.ch): > All of these are related to an API change of the ZFinder when going > from Rivet 1.5 to 1.6 a bit over a year ago. I can deal with the > transition of that if it's the only problem that is left in the > analysis. sounds good. I've attached a half cleaned up version of the code for you. If you can fix the ZFinder stuff, then I'll take care of all the division-by-zero and NaN problems, as well as the manual bin division. Cheers, Hendrik -- If your dreams don't scare you, then you are not dreaming big enough. -------------- next part -------------- // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/ParticleName.hh" namespace Rivet { class CMS_2012_I1107658 : public Analysis { public: /// Constructor CMS_2012_I1107658() : Analysis("CMS_2012_I1107658") {} void init() { ZFinder zfinder(-MAXRAPIDITY, MAXRAPIDITY, 0.0*GeV, MUON, 4.0*GeV, 140.0*GeV, 0.2, false, false); addProjection(zfinder, "ZFinder"); ChargedFinalState cfs(-2.0, 2.0, 500*MeV); // For charged particles addProjection(cfs, "CFS"); _h_Nchg_towards_pTmumu = bookProfile1D(1, 1, 1); _h_Nchg_transverse_pTmumu = bookProfile1D(2, 1, 1); _h_Nchg_away_pTmumu = bookProfile1D(3, 1, 1); _h_pTsum_towards_pTmumu = bookProfile1D(4, 1, 1); _h_pTsum_transverse_pTmumu = bookProfile1D(5, 1, 1); _h_pTsum_away_pTmumu = bookProfile1D(6, 1, 1); _h_avgpT_towards_pTmumu = bookDataPointSet(7, 1, 1); _h_avgpT_transverse_pTmumu = bookDataPointSet(8, 1, 1); _h_avgpT_away_pTmumu = bookDataPointSet(9, 1, 1); _h_Nchg_towards_plus_transverse_Mmumu = bookProfile1D(10, 1, 1); _h_pTsum_towards_plus_transverse_Mmumu = bookProfile1D(11, 1, 1); _h_avgpT_towards_plus_transverse_Mmumu = bookDataPointSet(12, 1, 1); _h_Nchg_towards_zmass_81_101 = bookHistogram1D(13, 1, 1); _h_Nchg_transverse_zmass_81_101 = bookHistogram1D(14, 1, 1); _h_Nchg_away_zmass_81_101 = bookHistogram1D(15, 1, 1); _h_pT_towards_zmass_81_101 = bookHistogram1D(16, 1, 1); _h_pT_transverse_zmass_81_101 = bookHistogram1D(17, 1, 1); _h_pT_away_zmass_81_101 = bookHistogram1D(18, 1, 1); _h_Nchg_transverse_zpt_5 = bookHistogram1D(19, 1, 1); _h_pT_transverse_zpt_5 = bookHistogram1D(20, 1, 1); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); const ZFinder& zfinder = applyProjection<ZFinder>(event, "ZFinder"); const FinalState& fs = zfinder.constituentsFinalState(); if (zfinder.particles().size() != 1) vetoEvent; Particle lepton0 = fs.particles().at(0); Particle lepton1 = fs.particles().at(1); if (lepton0.pdgId() != -lepton1.pdgId()) vetoEvent; double pt0 = lepton0.momentum().pT()/GeV; double pt1 = lepton1.momentum().pT()/GeV; double eta0 = lepton0.momentum().eta(); double eta1 = lepton1.momentum().eta(); if (pt0 > 20. && pt1 > 20. && fabs(eta1) < 2.4 && fabs(eta0) < 2.4) { double Zpt = zfinder.particles()[0].momentum().pT()/GeV; double Zphi = zfinder.particles()[0].momentum().phi(); double Zmass = zfinder.particles()[0].momentum().mass()/GeV; ParticleVector particles = applyProjection<ChargedFinalState>(event, "CFS").particlesByPt(); double nTowards(0.0), ptSumTowards(0.0), nTransverse(0.0), ptSumTransverse(0.0), nAway(0.0), ptSumAway(0.0); foreach (const Particle& p, particles) { if (abs(p.pdgId()) != 13) { if (fabs(p.momentum().eta()) > 2.) continue; double _dphi = fabs(deltaPhi(Zphi, p.momentum().phi())); if ( _dphi < M_PI/3.0 ) { // Towards region nTowards += 1.0; ptSumTowards += p.momentum().pT(); double _pT = p.momentum().pT(); if (Zmass < 101. && Zmass > 81.) _h_pT_towards_zmass_81_101->fill(_pT, weight); } else if ( _dphi < 2.*M_PI/3.0 ) { // Transverse region nTransverse += 1.0; ptSumTransverse += p.momentum().pT(); double _pT = p.momentum().pT(); if (Zmass < 101. && Zmass > 81.) _h_pT_transverse_zmass_81_101->fill(_pT, weight); if (Zpt < 5.) _h_pT_transverse_zpt_5->fill(_pT, weight); } else { nAway += 1.0; ptSumAway += p.momentum().pT(); double _pT = p.momentum().pT(); if (Zmass < 101. && Zmass > 81.) _h_pT_away_zmass_81_101->fill(_pT, weight); } } // not muons } // Loop over particles const double area = 8./3.*M_PI; if (Zmass < 101. && Zmass > 81.) { _h_Nchg_towards_pTmumu-> fill(Zpt/GeV, 1./area * nTowards, weight); _h_Nchg_transverse_pTmumu-> fill(Zpt/GeV, 1./area * nTransverse, weight); _h_Nchg_away_pTmumu-> fill(Zpt/GeV, 1./area * nAway, weight); _h_pTsum_towards_pTmumu-> fill(Zpt/GeV, 1./area * ptSumTowards, weight); _h_pTsum_transverse_pTmumu-> fill(Zpt/GeV, 1./area * ptSumTransverse, weight); _h_pTsum_away_pTmumu-> fill(Zpt/GeV, 1./area * ptSumAway, weight); _h_Nchg_towards_zmass_81_101-> fill(nTowards, weight); _h_Nchg_transverse_zmass_81_101->fill(nTransverse, weight); _h_Nchg_away_zmass_81_101-> fill(nAway, weight); } if (Zpt < 5.) { _h_Nchg_towards_plus_transverse_Mmumu->fill(Zmass/GeV, (nTowards + nTransverse)/(2*area), weight); _h_pTsum_towards_plus_transverse_Mmumu->fill(Zmass/GeV, (ptSumTowards + ptSumTransverse)/(2*area), weight); _h_Nchg_transverse_zpt_5->fill(nTransverse, weight); } } // cut on leptons } /// Normalise histograms etc., after the run void finalize() { vector<double> yval_avgpt_towards, yerr_avgpt_towards, yval_avgpt_transverse, yerr_avgpt_transverse, yval_avgpt_away, yerr_avgpt_away, yval_avgpt_tt, yerr_avgpt_tt; for (size_t i = 0; i < 20; ++i) { double yval_tmp = _h_pTsum_towards_pTmumu->binHeight(i) / _h_Nchg_towards_pTmumu->binHeight(i); yval_avgpt_towards.push_back(yval_tmp); double yerr_tmp = sqrt( _h_pTsum_towards_pTmumu->binError(i) * _h_pTsum_towards_pTmumu->binError(i)/ (_h_pTsum_towards_pTmumu->binHeight(i) * _h_pTsum_towards_pTmumu->binHeight(i)) + _h_Nchg_towards_pTmumu->binError(i) * _h_Nchg_towards_pTmumu->binError(i)/ (_h_Nchg_towards_pTmumu->binHeight(i) * _h_Nchg_towards_pTmumu->binHeight(i))) * yval_tmp; yerr_avgpt_towards.push_back(yerr_tmp); yval_tmp = _h_pTsum_transverse_pTmumu->binHeight(i) / _h_Nchg_transverse_pTmumu->binHeight(i); yval_avgpt_transverse.push_back(yval_tmp); yerr_tmp = sqrt(_h_pTsum_transverse_pTmumu->binError(i) * _h_pTsum_transverse_pTmumu->binError(i)/ (_h_pTsum_transverse_pTmumu->binHeight(i) * _h_pTsum_transverse_pTmumu->binHeight(i)) + _h_Nchg_transverse_pTmumu->binError(i) * _h_Nchg_transverse_pTmumu->binError(i)/ (_h_Nchg_transverse_pTmumu->binHeight(i) * _h_Nchg_transverse_pTmumu->binHeight(i))) * yval_tmp; yerr_avgpt_transverse.push_back(yerr_tmp); yval_tmp = _h_pTsum_away_pTmumu->binHeight(i) / _h_Nchg_away_pTmumu->binHeight(i); yval_avgpt_away.push_back(yval_tmp); yerr_tmp = sqrt(_h_pTsum_away_pTmumu->binError(i) * _h_pTsum_away_pTmumu->binError(i)/ (_h_pTsum_away_pTmumu->binHeight(i) * _h_pTsum_away_pTmumu->binHeight(i)) + _h_Nchg_away_pTmumu->binError(i) * _h_Nchg_away_pTmumu->binError(i)/ (_h_Nchg_away_pTmumu->binHeight(i) * _h_Nchg_away_pTmumu->binHeight(i))) * yval_tmp; yerr_avgpt_away.push_back(yerr_tmp); } for (size_t i = 0; i < 10; ++i) { double yval_tmp = _h_pTsum_towards_plus_transverse_Mmumu->binHeight(i) / _h_Nchg_towards_plus_transverse_Mmumu->binHeight(i); yval_avgpt_tt.push_back(yval_tmp); double yerr_tmp = sqrt(_h_pTsum_towards_plus_transverse_Mmumu->binError(i) * _h_pTsum_towards_plus_transverse_Mmumu->binError(i)/ (_h_pTsum_towards_plus_transverse_Mmumu->binHeight(i) * _h_pTsum_towards_plus_transverse_Mmumu->binHeight(i)) + _h_Nchg_towards_plus_transverse_Mmumu->binError(i) * _h_Nchg_towards_plus_transverse_Mmumu->binError(i)/ (_h_Nchg_towards_plus_transverse_Mmumu->binHeight(i) * _h_Nchg_towards_plus_transverse_Mmumu->binHeight(i))) * yval_tmp; yerr_avgpt_tt.push_back(yerr_tmp); } _h_avgpT_towards_pTmumu->setCoordinate(1, yval_avgpt_towards, yerr_avgpt_towards); _h_avgpT_transverse_pTmumu->setCoordinate(1, yval_avgpt_transverse, yerr_avgpt_transverse); _h_avgpT_away_pTmumu->setCoordinate(1, yval_avgpt_away, yerr_avgpt_away); _h_avgpT_towards_plus_transverse_Mmumu->setCoordinate(1, yval_avgpt_tt, yerr_avgpt_tt); // scale(_h_pT_towards_zmass_81_101, 1.0/integral(_h_Nchg_towards_zmass_81_101)); // scale(_h_pT_transverse_zmass_81_101, 1.0/integral(_h_Nchg_transverse_zmass_81_101)); // scale(_h_pT_away_zmass_81_101, 1.0/integral(_h_Nchg_away_zmass_81_101)); // scale(_h_pT_transverse_zpt_5, 1.0/integral(_h_Nchg_transverse_zpt_5)); // scale(_h_Nchg_towards_zmass_81_101, 1.0/integral(_h_Nchg_towards_zmass_81_101)); // scale(_h_Nchg_transverse_zmass_81_101, 1.0/integral(_h_Nchg_transverse_zmass_81_101)); // scale(_h_Nchg_away_zmass_81_101, 1.0/integral(_h_Nchg_away_zmass_81_101)); // scale(_h_Nchg_transverse_zpt_5, 1.0/integral(_h_Nchg_transverse_zpt_5)); } private: AIDA::IProfile1D *_h_Nchg_towards_pTmumu; AIDA::IProfile1D *_h_Nchg_transverse_pTmumu; AIDA::IProfile1D *_h_Nchg_away_pTmumu; AIDA::IProfile1D *_h_pTsum_towards_pTmumu; AIDA::IProfile1D *_h_pTsum_transverse_pTmumu; AIDA::IProfile1D *_h_pTsum_away_pTmumu; AIDA::IDataPointSet* _h_avgpT_towards_pTmumu; AIDA::IDataPointSet* _h_avgpT_transverse_pTmumu; AIDA::IDataPointSet* _h_avgpT_away_pTmumu; AIDA::IProfile1D *_h_Nchg_towards_plus_transverse_Mmumu; AIDA::IProfile1D *_h_pTsum_towards_plus_transverse_Mmumu; AIDA::IDataPointSet*_h_avgpT_towards_plus_transverse_Mmumu; AIDA::IHistogram1D *_h_Nchg_towards_zmass_81_101; AIDA::IHistogram1D *_h_Nchg_transverse_zmass_81_101; AIDA::IHistogram1D *_h_Nchg_away_zmass_81_101; AIDA::IHistogram1D *_h_pT_towards_zmass_81_101; AIDA::IHistogram1D *_h_pT_transverse_zmass_81_101; AIDA::IHistogram1D *_h_pT_away_zmass_81_101; AIDA::IHistogram1D *_h_Nchg_transverse_zpt_5; AIDA::IHistogram1D *_h_pT_transverse_zpt_5; //@} }; // This global object acts as a hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2012_I1107658); }
More information about the Rivet mailing list |