|
[Rivet-svn] r2873 - in trunk: . src/Analysesblackhole at projects.hepforge.org blackhole at projects.hepforge.orgTue Jan 11 12:16:27 GMT 2011
Author: buckley Date: Tue Jan 11 12:16:26 2011 New Revision: 2873 Log: Much-tidied ATLAS dijet analysis Modified: trunk/ChangeLog trunk/src/Analyses/ATLAS_2010_S8817804.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Tue Jan 11 11:47:53 2011 (r2872) +++ trunk/ChangeLog Tue Jan 11 12:16:26 2011 (r2873) @@ -1,3 +1,7 @@ +2011-01-11 Andy Buckley <andy at insectnation.org> + + * Clean up ATLAS dijet analysis. + 2010-12-30 Andy Buckley <andy at insectnation.org> * Adding a run timeout option, and small bug-fixes to the event Modified: trunk/src/Analyses/ATLAS_2010_S8817804.cc ============================================================================== --- trunk/src/Analyses/ATLAS_2010_S8817804.cc Tue Jan 11 11:47:53 2011 (r2872) +++ trunk/src/Analyses/ATLAS_2010_S8817804.cc Tue Jan 11 12:16:26 2011 (r2873) @@ -6,167 +6,124 @@ #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FastJets.hh" -namespace Rivet{ +namespace Rivet { + - class ATLAS_2010_S8817804; - - AnalysisBuilder<ATLAS_2010_S8817804> plugin_ATLAS_2010_S8817804; - /// @brief ATLAS inclusive jet pT spectrum, di-jet Mass and di-jet chi - class ATLAS_2010_S8817804: public Analysis{ - public: - - ATLAS_2010_S8817804(): Analysis("ATLAS_2010_S8817804"), _rapMax(2.8), _leadPTCut(60. * GeV), _asymmetryCut(log(30.)), _ySumCut(2.2){ + + ATLAS_2010_S8817804() + : Analysis("ATLAS_2010_S8817804") + { setBeams(PROTON, PROTON); setNeedsCrossSection(true); } - - void init(){ + + + private: + + enum Alg { AKT4=0, AKT6=1 }; + + + public: + + void init() { FinalState fs; addProjection(fs, "FinalState"); addProjection(FastJets(fs, FastJets::ANTIKT, 0.6), "AntiKT06"); - addProjection(FastJets(fs, FastJets::ANTIKT, 0.4), "AntiKT04"); + addProjection(FastJets(fs, FastJets::ANTIKT, 0.4), "AntiKT04"); - double ybins[] = {0.0, 0.3, 0.8, 1.2, 2.1, 2.8}; + double ybins[] = { 0.0, 0.3, 0.8, 1.2, 2.1, 2.8 }; + double massBinsForChi[] = { 340, 520, 800, 1200 }; - double massBinsForChi[] = { 340. * GeV, 520. * GeV, 800. * GeV, 1200. * GeV};//, 1700. * GeV, 2500. * GeV}; - - int ptDsOffset = 1; - int massDsOffset = 11; - int chiDsOffset = 21; - - for(Jet_Alg alg = AKT4; alg != ALG_END; alg = Jet_Alg(alg + 1)){ - for(int highbin = 0, lowbin = highbin++; highbin != sizeof(ybins) / sizeof(double); ++lowbin, ++highbin){ - - _pTHistos[alg].addHistogram(ybins[lowbin], ybins[highbin], - bookHistogram1D(lowbin + ptDsOffset, 1, 1)); - } - - ptDsOffset = 6; - - for(int highbin = 0, lowbin = highbin++; highbin != sizeof(ybins) / sizeof(double); ++lowbin, ++highbin){ - - _massVsY[alg].addHistogram(ybins[lowbin], ybins[highbin], - bookHistogram1D(massDsOffset + lowbin, 1, 1)); - } - - massDsOffset = 16; - - for(int highbin=0, lowbin = highbin++; highbin != sizeof(massBinsForChi) / sizeof(double); ++lowbin, ++highbin){ - - _chiVsMass[alg].addHistogram(massBinsForChi[lowbin], massBinsForChi[highbin], - bookHistogram1D(chiDsOffset + lowbin, 1, 1)); - - } - - chiDsOffset = 24; - + for (size_t alg = 0; alg < 2; ++alg) { + int ptDsOffset = 1; + for (size_t i = 0; i < 5; ++i) { + _pTHistos[alg].addHistogram(ybins[i], ybins[i+1], bookHistogram1D(i + ptDsOffset, 1, 1)); + } + ptDsOffset += 5; + + int massDsOffset = 11; + for (size_t i = 0; i < 5; ++i) { + _massVsY[alg].addHistogram(ybins[i], ybins[i+1], bookHistogram1D(i + massDsOffset, 1, 1)); + } + massDsOffset += 5; + + int chiDsOffset = 21; + for (size_t i = 0; i < 3; ++i) { + _chiVsMass[alg].addHistogram(massBinsForChi[i], massBinsForChi[i+1], bookHistogram1D(i + chiDsOffset, 1, 1)); + } + chiDsOffset += 3; } - - return; } - - void analyze(const Event &evt){ - - const Jets &kt6Jets = applyProjection<FastJets>(evt, "AntiKT06").jets(30.*GeV); - const Jets &kt4Jets = applyProjection<FastJets>(evt, "AntiKT04").jets(30.*GeV); - - const Jets* jetAr[2]; - jetAr[AKT4] = &kt4Jets; - jetAr[AKT6] = &kt6Jets; - - for(Jet_Alg alg = AKT4; alg != ALG_END; alg = Jet_Alg(alg + 1)){ - - const Jet *jet1 = 0, *jet2 = 0; - double pt1 = -1.; - double pt2 = -1.; - - for(Jets::const_iterator jet = jetAr[alg]->begin(); - jet != jetAr[alg]->end(); ++jet){ - double pT = jet->momentum().pT(); - double fabsy = fabs(jet->momentum().rapidity()); - _pTHistos[alg].fill(fabsy, pT, evt.weight()); - - if(fabsy < _rapMax){ - if(pT > pt1){ - if(jet1 != 0){ - jet2 = jet1; - pt2 = jet2->momentum().pT(); - } - jet1 = &(*jet); - pt1 = pT; - }else if(pT > pt2){ - jet2 = &(*jet); - pt2 = pT; - } + + + void analyze(const Event& evt) { + Jets jetAr[2]; + jetAr[AKT4] = applyProjection<FastJets>(evt, "AntiKT06").jetsByPt(30*GeV); + jetAr[AKT6] = applyProjection<FastJets>(evt, "AntiKT04").jetsByPt(30*GeV); + + // Cut values + const double rapMax(2.8), leadPTCut(60*GeV), asymmetryCut(log(30.)), ySumCut(2.2); + + for (size_t alg = 0; alg < 2; ++alg) { + vector<FourMomentum> leadjets; + foreach (const Jet& jet, jetAr[alg]) { + double pT = jet.momentum().pT(); + double fabsy = fabs(jet.momentum().rapidity()); + _pTHistos[alg].fill(fabsy, pT/GeV, evt.weight()); + + if (fabsy < rapMax && leadjets.size() < 2) { + if (leadjets.empty() && pT < leadPTCut) continue; + leadjets.push_back(jet.momentum()); } } - - if( jet1 == 0 || jet2 == 0 || pt1 < _leadPTCut){ + + if (leadjets.size() < 2) { + MSG_DEBUG("Could not find two suitable leading jets"); continue; } - - FourMomentum mom1 = jet1->momentum(); - FourMomentum mom2 = jet2->momentum(); - FourMomentum total = mom1 + mom2; - double rap1 = mom1.rapidity(); - - double absRap1 = fabs(rap1); - - double rap2 = mom2.rapidity(); - double absRap2 = fabs(rap2); - double ymax = (absRap1 > absRap2)? absRap1: absRap2; - double chi = exp(fabs(rap1 - rap2)); - double mass = total.mass(); - - if(fabs(rap1 + rap2) < _ySumCut && fabs(rap1 - rap2) < _asymmetryCut){ - _chiVsMass[alg].fill(mass, chi, evt.weight()); - } - - _massVsY[alg].fill(ymax, mass, evt.weight()); - + + const double rap1 = leadjets[0].rapidity(); + const double rap2 = leadjets[1].rapidity(); + const double mass = (leadjets[0] + leadjets[1]).mass(); + const double ymax = (fabs(rap1) > fabs(rap2)) ? fabs(rap1) : fabs(rap2); + const double chi = exp(fabs(rap1 - rap2)); + + if (fabs(rap1 + rap2) < ySumCut && fabs(rap1 - rap2) < asymmetryCut) { + _chiVsMass[alg].fill(mass/GeV, chi, evt.weight()); + } + _massVsY[alg].fill(ymax, mass/GeV, evt.weight()); + } - - return; } - - void finalize(){ - + + + void finalize() { const double xs = crossSectionPerEvent() / picobarn; - - for(Jet_Alg alg = AKT4; alg != ALG_END; alg = Jet_Alg(alg + 1)){ + for (size_t alg = 0; alg < 2; ++alg) { _pTHistos[alg].scale(xs, this); _massVsY[alg].scale(xs, this); _chiVsMass[alg].scale(xs, this); } - - return; } - + + private: - - // The max rapidity of a jet used in this analysis - double _rapMax; - // The pT cut on the leading jet - double _leadPTCut; - // maximum rapidity asymmetry between two jets in the chi plots - double _asymmetryCut; - // maximum sum of rapidities of two jets in the chi plot - double _ySumCut; - - enum Jet_Alg{AKT4=0, AKT6=1, ALG_END}; - - // The inclusive pT spectrum for akt6 and akt4 jets (array index is jet type from enum above) + + /// The inclusive pT spectrum for akt6 and akt4 jets (array index is jet type from enum above) BinnedHistogram<double> _pTHistos[2]; - - // The di-jet mass spectrum binned in rapidity for akt6 and akt4 jets (array index is jet type from enum above) + + /// The di-jet mass spectrum binned in rapidity for akt6 and akt4 jets (array index is jet type from enum above) BinnedHistogram<double> _massVsY[2]; - - // The di-jet chi distribution binned in mass for akt6 and akt4 jets (array index is jet type from enum above) + + /// The di-jet chi distribution binned in mass for akt6 and akt4 jets (array index is jet type from enum above) BinnedHistogram<double> _chiVsMass[2]; - - + }; + + + AnalysisBuilder<ATLAS_2010_S8817804> plugin_ATLAS_2010_S8817804; + }
More information about the Rivet-svn mailing list |