diff --git a/analyses/pluginATLAS/ATLAS_2014_I1288706.cc b/analyses/pluginATLAS/ATLAS_2014_I1288706.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1288706.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1288706.cc @@ -1,100 +1,101 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Particle.fhh" namespace Rivet { class ATLAS_2014_I1288706 : public Analysis { public: /// Constructor ATLAS_2014_I1288706() : Analysis("ATLAS_2014_I1288706") { } public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Set up projections FinalState fs; ZFinder zfinder_ext_dressed_mu(fs, Cuts::abseta<2.4 && Cuts::pT>6.0*GeV, PID::MUON, 12.0*GeV, 66.0*GeV, 0.1); declare(zfinder_ext_dressed_mu, "ZFinder_ext_dressed_mu"); ZFinder zfinder_dressed_mu(fs, Cuts::abseta<2.4 && Cuts::pT>12*GeV, PID::MUON, 26.0*GeV, 66.0*GeV, 0.1); declare(zfinder_dressed_mu, "ZFinder_dressed_mu"); ZFinder zfinder_dressed_el(fs, Cuts::abseta<2.4 && Cuts::pT>12*GeV, PID::ELECTRON, 26.0*GeV, 66.0*GeV, 0.1); declare(zfinder_dressed_el, "ZFinder_dressed_el"); book(_hist_ext_mu_dressed ,1, 1, 1); // muon, dressed level, extended phase space book(_hist_mu_dressed ,2, 1, 1); // muon, dressed level, nominal phase space book(_hist_el_dressed ,2, 1, 2); // electron, dressed level, nominal phase space } /// Perform the per-event analysis void analyze(const Event& event) { const ZFinder& zfinder_ext_dressed_mu = apply(event, "ZFinder_ext_dressed_mu"); const ZFinder& zfinder_dressed_mu = apply(event, "ZFinder_dressed_mu" ); const ZFinder& zfinder_dressed_el = apply(event, "ZFinder_dressed_el" ); FillPlots(zfinder_ext_dressed_mu, _hist_ext_mu_dressed, 9.0); FillPlots(zfinder_dressed_mu, _hist_mu_dressed, 15.0); FillPlots(zfinder_dressed_el, _hist_el_dressed, 15.0); } void FillPlots(const ZFinder& zfinder, Histo1DPtr hist, double leading_pT) { if(zfinder.bosons().size() != 1) return; + if(zfinder.particles().size() < 2) return; const FourMomentum el1 = zfinder.particles()[0].momentum(); const FourMomentum el2 = zfinder.particles()[1].momentum(); if (el1.pT() > leading_pT*GeV || el2.pT() > leading_pT*GeV) { double mass = zfinder.bosons()[0].mass()/GeV; hist->fill(mass); } } /// Normalise histograms etc., after the run void finalize() { scale(_hist_ext_mu_dressed, crossSection()/sumOfWeights()); scale(_hist_mu_dressed, crossSection()/sumOfWeights()); scale(_hist_el_dressed, crossSection()/sumOfWeights()); } //@} private: /// @name Histograms //@{ Histo1DPtr _hist_ext_mu_dressed; Histo1DPtr _hist_mu_dressed; Histo1DPtr _hist_el_dressed; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1288706); } diff --git a/analyses/pluginATLAS/ATLAS_2016_I1419070.cc b/analyses/pluginATLAS/ATLAS_2016_I1419070.cc --- a/analyses/pluginATLAS/ATLAS_2016_I1419070.cc +++ b/analyses/pluginATLAS/ATLAS_2016_I1419070.cc @@ -1,139 +1,139 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { class ATLAS_2016_I1419070 : public Analysis { public: /// Constructor ATLAS_2016_I1419070() : Analysis("ATLAS_2016_I1419070") { } public: void init() { addProjection(FastJets(FinalState(), FastJets::ANTIKT, 0.4), "Jets"); book(forward_500MeV ,1, 1, 1); book(forward_2GeV ,2, 1, 1); book(forward_5GeV ,3, 1, 1); book(central_500MeV ,4, 1, 1); book(central_2GeV ,5, 1, 1); book(central_5GeV ,6, 1, 1); book(diff_500MeV, "d07-x01-y01", true); book(diff_2GeV , "d08-x01-y01", true); book(diff_5GeV , "d09-x01-y01", true); book(sum_500MeV, "d10-x01-y01", true); book(sum_2GeV , "d11-x01-y01", true); book(sum_5GeV , "d12-x01-y01", true); } /// Perform the per-event analysis void analyze(const Event& event) { Jets m_goodJets = applyProjection(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.1); if (m_goodJets.size() < 2) vetoEvent; if (m_goodJets[0].pT() < 50*GeV) vetoEvent; if (m_goodJets[1].pT() < 50*GeV) vetoEvent; if (fabs(1.0 - m_goodJets[0].pT()/m_goodJets[1].pT()) > 0.5) vetoEvent; bool check = m_goodJets[0].abseta() < m_goodJets[1].abseta(); int pos_f = int(check); int pos_c = int(!check); double pt500MeV_f = CalculateNCharge(m_goodJets[pos_f], 0.5); double pt2GeV_f = CalculateNCharge(m_goodJets[pos_f], 2.0); double pt5GeV_f = CalculateNCharge(m_goodJets[pos_f], 5.0); double pT_f = m_goodJets[pos_f].pT(); double pt500MeV_c = CalculateNCharge(m_goodJets[pos_c], 0.5); double pt2GeV_c = CalculateNCharge(m_goodJets[pos_c], 2.0); double pt5GeV_c = CalculateNCharge(m_goodJets[pos_c], 5.0); double pT_c = m_goodJets[pos_c].pT(); forward_500MeV->fill(pT_f, pt500MeV_f); forward_2GeV->fill( pT_f, pt2GeV_f); forward_5GeV->fill( pT_f, pt5GeV_f); central_500MeV->fill(pT_c, pt500MeV_c); central_2GeV->fill( pT_c, pt2GeV_c); central_5GeV->fill( pT_c, pt5GeV_c); } double CalculateNCharge(Jet& jet, double pTcut=0.5) { unsigned int ncharge = 0; foreach (const Particle& p, jet.particles()) { if (p.pT() < pTcut) continue; if (p.threeCharge()) ++ncharge; } if (ncharge > 60) ncharge = 60; return double(ncharge); } /// Normalise histograms etc., after the run void finalize() { if (numEvents() > 2) { for (unsigned int i = 0; i < forward_500MeV->numBins(); ++i) { ProfileBin1D bsum = central_500MeV->bin(i) + forward_500MeV->bin(i); ProfileBin1D bsum2 = central_2GeV->bin(i) + forward_2GeV->bin(i); ProfileBin1D bsum5 = central_5GeV->bin(i) + forward_5GeV->bin(i); ProfileBin1D bdiff = central_500MeV->bin(i) - forward_500MeV->bin(i); ProfileBin1D bdiff2 = central_2GeV->bin(i) - forward_2GeV->bin(i); ProfileBin1D bdiff5 = central_5GeV->bin(i) - forward_5GeV->bin(i); double ydiff = central_500MeV->bin(i).numEntries()? central_500MeV->bin(i).mean() : 0.0; double ydiff2 = central_2GeV->bin(i).numEntries()? central_2GeV->bin(i).mean() : 0.0; double ydiff5 = central_5GeV->bin(i).numEntries()? central_5GeV->bin(i).mean() : 0.0; ydiff -= forward_500MeV->bin(i).numEntries()? forward_500MeV->bin(i).mean() : 0.0; ydiff2 -= forward_2GeV->bin(i).numEntries()? forward_2GeV->bin(i).mean() : 0.0; ydiff5 -= forward_5GeV->bin(i).numEntries()? forward_5GeV->bin(i).mean() : 0.0; - double yerr = bsum.numEntries()? bsum.stdErr() : 0.0; - double yerr2 = bsum2.numEntries()? bsum2.stdErr() : 0.0; - double yerr5 = bsum5.numEntries()? bsum5.stdErr() : 0.0; + double yerr = bsum.numEntries() > 1.0 ? bsum.stdErr() : 0.0; + double yerr2 = bsum2.numEntries() > 1.0 ? bsum2.stdErr() : 0.0; + double yerr5 = bsum5.numEntries() > 1.0 ? bsum5.stdErr() : 0.0; diff_500MeV->point(i).setY(ydiff, yerr); diff_2GeV->point(i).setY(ydiff2, yerr2); diff_5GeV->point(i).setY(ydiff5, yerr5); sum_500MeV->point(i).setY(bsum.numEntries()? bsum.mean() : 0.0, yerr); sum_2GeV->point(i).setY(bsum2.numEntries()? bsum2.mean() : 0.0, yerr2); sum_5GeV->point(i).setY(bsum5.numEntries()? bsum5.mean() : 0.0, yerr5); } } } private: Profile1DPtr forward_500MeV; Profile1DPtr forward_2GeV; Profile1DPtr forward_5GeV; Profile1DPtr central_500MeV; Profile1DPtr central_2GeV; Profile1DPtr central_5GeV; Scatter2DPtr sum_500MeV; Scatter2DPtr sum_2GeV; Scatter2DPtr sum_5GeV; Scatter2DPtr diff_500MeV; Scatter2DPtr diff_2GeV; Scatter2DPtr diff_5GeV; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2016_I1419070); } diff --git a/analyses/pluginCMS/CMS_2014_I1305624.cc b/analyses/pluginCMS/CMS_2014_I1305624.cc --- a/analyses/pluginCMS/CMS_2014_I1305624.cc +++ b/analyses/pluginCMS/CMS_2014_I1305624.cc @@ -1,910 +1,919 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { namespace { /// Number of event shape variables /// @todo Move into the EventShape class const int NEVTVAR = 5; /// Number of leading jet pT thresholds /// @todo Move into the analysis class const int NJETPTMN = 5; /// Leading jet pT thresholds /// @todo Move into the analysis class const double LEADINGPTTHRESHOLD[NJETPTMN] = { 110.0, 170.0, 250.0, 320.0, 390.0 }; // Helpers for event shape calculations in hidden namespace; implementation at bottom of file /// @todo Why a class? Improve/remove this junk class EventShape { public: /// Constructor from vectors of four-vectors as input objects in the event to calculate the event shapes EventShape(const vector& px_vector, const vector& py_vector, const vector& pz_vector, const vector& e_vector, double eta_central, int irap, int nmn) : _object_px(px_vector), _object_py(py_vector), _object_pz(pz_vector), _object_e(e_vector), _eta_c(eta_central), _irap(irap), _nmnjet(nmn) { } /// @brief Returns the values of the five event shapes /// /// Event shape indices: /// 0. central transverse thrust /// 1. central total jet broadening /// 2. central total jet mass /// 3. central total transverse jet mass /// 4. central three-jet resolution threshold vector getEventShapes() { _calculate(); ///< @todo There should be some test for success/failure!! return _event_shapes; } /// Returns the global thrust axis Nx, Ny, Nz=0 vector getThrustAxis() { _calculate(); ///< @todo There should be some test for success/failure!! return _thrust_axis; } /// Returns the central thrust axis Nx, Ny, Nz=0 vector getThrustAxisC() { _calculate(); ///< @todo There should be some test for success/failure!! return _thrust_axis_c; } // /// @brief Choice of the central region // void setEtaC(double eta_central) { _eta_c = eta_central; } // // Whether to use the rapidity y (rap==1) or the pseudorapidity eta (rap==0) // void setRapType(int irap) { _irap = irap; } private: /// Calculate everything int _calculate(); /// Returns the difference in phi between two vectors double _delta_phi(double, double); /// The Lorentz scalar product double _lorentz_sp(const vector&, const vector&); // Calculates the three-jet resolutions double _three_jet_res(const vector&, const vector&, const vector&, const vector&, int); // Calculates the thrust axis and the tau values vector _thrust(const vector&, const vector&); vector _object_px, _object_py, _object_pz, _object_p; vector _object_pt, _object_e, _object_phi, _object_eta; vector _event_shapes; vector _thrust_axis, _thrust_axis_c; double _eta_c; int _irap; size_t _nmnjet; }; } class CMS_2014_I1305624 : public Analysis { public: /// Constructor CMS_2014_I1305624() : Analysis("CMS_2014_I1305624") { } /// @name Analysis methods /// Book histograms and initialise projections before the run void init() { const FastJets jets(FinalState(Cuts::abseta < 2.6), FastJets::ANTIKT, 0.5); declare(jets, "Jets"); for (int ij=0; ij < NJETPTMN; ij++) { book(_h_thrustc[ij] ,1, 1, ij+1); book(_h_broadt[ij] ,1, 2, ij+1); book(_h_tot3dmass[ij] ,1, 3, ij+1); book(_h_tottrnsmass[ij] ,1, 4, ij+1); book(_h_y23c[ij] ,1, 5, ij+1); // + } + _needBinInit = true; + } + + + /// Perform the per-event analysis + void analyze(const Event& event) { + + if (_needBinInit) { + for (int ij=0; ij < NJETPTMN; ij++) { _alow1[ij] = _h_thrustc[ij]->xMin(); _alow2[ij] = _h_broadt[ij]->xMin(); _alow3[ij] = _h_tot3dmass[ij]->xMin(); _alow4[ij] = _h_tottrnsmass[ij]->xMin(); _alow5[ij] = _h_y23c[ij]->xMin(); // _ahgh1[ij] = _h_thrustc[ij]->xMax(); _ahgh2[ij] = _h_broadt[ij]->xMax(); _ahgh3[ij] = _h_tot3dmass[ij]->xMax(); _ahgh4[ij] = _h_tottrnsmass[ij]->xMax(); _ahgh5[ij] = _h_y23c[ij]->xMax(); + + _needBinInit = false; + } } - } - - - /// Perform the per-event analysis - void analyze(const Event& event) { const Jets& jets = apply(event, "Jets").jetsByPt(30.0*GeV); if (jets.size() < 2) vetoEvent; if (jets[0].abseta() > 2.4 || jets[1].abseta() > 2.4) vetoEvent; const double leadingpt = jets[0].pT(); if (leadingpt < 110*GeV) vetoEvent; vector jtpx, jtpy, jtpz, jten; foreach (const Jet& j, jets) { if (j.abseta() < 2.4) { jtpx.push_back(j.px()); jtpy.push_back(j.py()); jtpz.push_back(j.pz()); jten.push_back(j.E()); } } EventShape eventshape(jtpx, jtpy, jtpz, jten, 2.4, 0, 2); const vector eventvar = eventshape.getEventShapes(); if (eventvar[NEVTVAR] < 0) vetoEvent; // Jets are not only one hemisphere const double weight = 1.0; for (int ij = NJETPTMN-1; ij >= 0; --ij) { if (leadingpt/GeV > LEADINGPTTHRESHOLD[ij]) { if (inRange(eventvar[0], _alow1[ij], _ahgh1[ij])) _h_thrustc[ij]->fill(eventvar[0], weight); if (inRange(eventvar[2], _alow3[ij], _ahgh3[ij])) _h_tot3dmass[ij]->fill(eventvar[2], weight); if (inRange(eventvar[3], _alow4[ij], _ahgh4[ij])) _h_tottrnsmass[ij]->fill(eventvar[3], weight); if (eventvar[NEVTVAR] >= 3) { if (inRange(eventvar[1], _alow2[ij], _ahgh2[ij])) _h_broadt[ij]->fill(eventvar[1], weight); if (inRange(eventvar[4], _alow5[ij], _ahgh5[ij])) _h_y23c[ij]->fill(eventvar[4], weight); } break; } } } /// Normalise histograms etc., after the run void finalize() { for (int ij = 0; ij < NJETPTMN; ij++) { normalize(_h_thrustc[ij]); normalize(_h_broadt[ij]); normalize(_h_tot3dmass[ij]); normalize(_h_tottrnsmass[ij]); normalize(_h_y23c[ij]); } } private: /// @name Histograms //@{ Histo1DPtr _h_thrustc[NJETPTMN]; Histo1DPtr _h_broadt[NJETPTMN]; Histo1DPtr _h_tot3dmass[NJETPTMN]; Histo1DPtr _h_tottrnsmass[NJETPTMN]; Histo1DPtr _h_y23c[NJETPTMN]; //@} // Data members + bool _needBinInit; double _alow1[NJETPTMN], _alow2[NJETPTMN], _alow3[NJETPTMN], _alow4[NJETPTMN], _alow5[NJETPTMN]; double _ahgh1[NJETPTMN], _ahgh2[NJETPTMN], _ahgh3[NJETPTMN], _ahgh4[NJETPTMN], _ahgh5[NJETPTMN]; }; DECLARE_RIVET_PLUGIN(CMS_2014_I1305624); ///////////////////// namespace { // EventShape helper class method implementations: int EventShape::_calculate() { if (!_event_shapes.empty() && !_thrust_axis.empty() && !_thrust_axis_c.empty()) return 1; //< return success if this appears to already have been run const size_t length = (size_t) _object_px.size(); if (((size_t) _object_py.size() != length) || ((size_t) _object_pz.size() != length) || ((size_t) _object_e.size() != length)) { /// @todo Change to exception or assert // cout << "ERROR!!!! Input vectors differ in size! Change that please!" << endl; // cout<<"py_size: "<<_object_py.size()<<" ,pz_size: "<<_object_pz.size() // <<" ,px_size: "<<_object_px.size()<<" ,E_size: "<<_object_e.size()< _object_e[k] + 1e-4) { /// @todo Change to exception or assert // cout << "ERROR!!! object " << k <<" has P = " << _object_p[k] // << " which is bigger than E = " << _object_e[k] <<" " // << _object_px[k] <<" "<< _object_py[k] <<" " // << _object_pz[k] <<" of total length "<< length // << endl; return 0; } //to prevent a division by zero if (_irap == 0) { if (fabs(_object_pz[k]) > 1e-5) { theta = atan(_object_pt[k]/(_object_pz[k])); } else { theta = M_PI/2; } if (theta < 0.) theta = theta + M_PI; _object_eta[k] = -log(tan(0.5*theta)); } if (_irap == 1) { if (_object_pz[k] == _object_e[k]) { /// @todo Change to exception or assert // cout << "ERROR!!! object "< object_px_in, object_py_in, object_pz_in, object_pt_in, object_e_in, object_et_in, object_eta_in; vector object_px_out, object_py_out, object_pz_out, object_e_out, object_pt_out, object_eta_out; if (!object_px_in.empty()) { //< FFS, this is impossible: it's only just been created! object_px_in.clear(); object_py_in.clear(); object_pz_in.clear(); object_pt_in.clear(); object_e_in.clear(); object_et_in.clear(); object_eta_in.clear(); object_px_out.clear(); object_py_out.clear(); object_pz_out.clear(); object_pt_out.clear(); object_e_out.clear(); object_eta_out.clear(); } size_t nin = 0; for (size_t j = 0; j < length; j++) { if (fabs(_object_eta[j]) < _eta_c) { object_px_in.push_back(_object_px[j]); object_py_in.push_back(_object_py[j]); object_pz_in.push_back(_object_pz[j]); object_e_in.push_back(_object_e[j]); object_pt_in.push_back(sqrt(pow(_object_px[j],2)+pow(_object_py[j],2))); object_et_in.push_back(sqrt((pow(_object_e[j],2)*pow(_object_pt[j],2))/(pow(_object_pt[j],2)+pow(_object_pz[j],2)))); object_eta_in.push_back(_object_eta[j]); nin += 1; } else { object_px_out.push_back(_object_px[j]); object_py_out.push_back(_object_py[j]); object_pz_out.push_back(_object_pz[j]); object_e_out.push_back(_object_e[j]); object_pt_out.push_back(sqrt(pow(_object_px[j],2)+pow(_object_py[j],2))); object_eta_out.push_back(_object_eta[j]); } } if (object_px_in.size() != nin) { /// @todo Change to exception or assert cout<<"ERROR!!! wrong dimension of 'in' momenta"<= _nmnjet) { double p_sum_c = 0; //GMA double pt_sum_c = 0; double eta_cw=0; double px_sum_in = 0; double py_sum_in = 0; for (size_t j = 0; j < nin; j++) { pt_sum_c += object_pt_in[j]; p_sum_c += sqrt(pow(object_pt_in[j],2.) + pow(object_pz_in[j], 2.0)); //GMA eta_cw += object_pt_in[j]*object_eta_in[j]; px_sum_in += object_px_in[j]; py_sum_in += object_py_in[j]; } eta_cw /= pt_sum_c; double expTerm = 0; for (size_t j = 0; j < nout; j++) { expTerm += object_pt_out[j] * exp(-fabs(object_eta_out[j]-eta_cw)); } expTerm /= pt_sum_c; //the central global transverse thrust centrthr is calculated double centrthr = 0; vector thrust_central = _thrust(object_px_in, object_py_in); for (size_t l=0; l<3; l++) _thrust_axis_c[l] = thrust_central[l]; //the variable which gets resummed is not thrust //but tau = 1 - thrust - see calculation centrthr = thrust_central[3]; _event_shapes[0] = centrthr; double alpha_c = atan2(_thrust_axis_c[1], _thrust_axis_c[0]); //central jet masses //define two jet masses in region U and D double cenjm_up = 0; double cenjm_down= 0; double dot_product = 0; vector up_sum; vector down_sum; for (size_t j=0; j<4;j++) { up_sum.push_back(0.); down_sum.push_back(0.); } for (size_t i=0;i= 0) { up_sum[0]+=object_px_in[i]; up_sum[1]+=object_py_in[i]; up_sum[2]+=object_pz_in[i]; up_sum[3]+=object_e_in[i]; } else { down_sum[0]+=object_px_in[i]; down_sum[1]+=object_py_in[i]; down_sum[2]+=object_pz_in[i]; down_sum[3]+=object_e_in[i]; } } cenjm_up = _lorentz_sp(up_sum, up_sum) / pow(p_sum_c, 2.); //GMA pow(pt_sum_c,2); cenjm_down = _lorentz_sp(down_sum, down_sum) / pow(p_sum_c, 2.); //GMA pow(pt_sum_c,2); //central total jet mass centotjm double centotjm=0; centotjm = cenjm_up + cenjm_down; _event_shapes[2]=centotjm; double centrjm_up=0, centrjm_down=0; vector upsum; vector downsum; for (size_t j = 0; j < 3; j++) { upsum.push_back(0.); downsum.push_back(0.); } for (size_t i = 0; i < nin; i++) { dot_product = object_px_in[i]*_thrust_axis_c[0]+object_py_in[i]*_thrust_axis_c[1]; if (dot_product >= 0) { upsum[0] += object_px_in[i]; upsum[1] += object_py_in[i]; upsum[2] += object_et_in[i]; } else { downsum[0] += object_px_in[i]; downsum[1] += object_py_in[i]; downsum[2] += object_et_in[i]; } } centrjm_up = _lorentz_sp(upsum, upsum) / pow(pt_sum_c, 2); centrjm_down = _lorentz_sp(downsum, downsum) / pow(pt_sum_c, 2); double centottrjm = centrjm_up + centrjm_down; _event_shapes[3] = centottrjm; //central three-jet resolution threshold double ceny3=0; if (nin < 3) { ceny3 = -1.0; } else { ceny3 = _three_jet_res(object_px_in, object_py_in, object_pz_in, object_e_in, _irap); } _event_shapes[4] = ceny3; //the central jet broadenings in the up and down region double cenbroad_up=0; double cenbroad_down=0; double eta_up=0; size_t num_up=0; double eta_down =0; size_t num_down =0; double phi_temp =0; double phi_up_aver =0; double phi_down_aver =0; double pt_sum_up =0; double pt_sum_down =0; double dot_product_b =0; vector phi_up; vector phi_down; double py_rot =0; double px_rot =0; for (size_t j = 0; j < 4; j++) { up_sum.push_back(0.); down_sum.push_back(0.); } for (size_t i=0;i=0){ pt_sum_up += object_pt_in[i]; //rotate the coordinate system so that //the central thrust axis is e_x px_rot = cos(alpha_c)*object_px_in[i]+sin(alpha_c)*object_py_in[i]; py_rot = - sin(alpha_c)*object_px_in[i]+cos(alpha_c)*object_py_in[i]; //calculate the eta and phi in the rotated system eta_up += object_pt_in[i]*object_eta_in[i]; phi_temp = atan2(py_rot,px_rot); if(phi_temp > M_PI/2){ phi_temp = phi_temp - M_PI/2; } if (phi_temp < -M_PI/2){ phi_temp = phi_temp + M_PI/2; } phi_up.push_back(phi_temp); phi_up_aver += object_pt_in[i]*phi_temp; num_up += 1; } else { eta_down += object_pt_in[i]*object_eta_in[i]; pt_sum_down += object_pt_in[i]; px_rot = cos(alpha_c)*object_px_in[i]+sin(alpha_c)*object_py_in[i]; py_rot = - sin(alpha_c)*object_px_in[i]+cos(alpha_c)*object_py_in[i]; phi_temp = atan2(py_rot,px_rot); if (phi_temp > M_PI/2) { //if phi is bigger than pi/2 in the new system calculate //the difference to the thrust axis phi_temp = M_PI -phi_temp; } if (phi_temp<-M_PI/2) { //if phi is smaller than phi_temp = -M_PI-phi_temp; } phi_down.push_back(phi_temp); //calculate the pt-weighted phi phi_down_aver += object_pt_in[i]*phi_temp; num_down += 1; } } if (num_up!=0){ eta_up = eta_up/pt_sum_up; phi_up_aver = phi_up_aver/pt_sum_up; } if (num_down!=0) { eta_down = eta_down/pt_sum_down; phi_down_aver = phi_down_aver/pt_sum_down; } size_t index_up=0, index_down=0; for (size_t i = 0; i < nin; i++) { dot_product_b = object_px_in[i]*_thrust_axis_c[0] + object_py_in[i]*_thrust_axis_c[1]; if (dot_product_b >= 0) { //calculate the broadenings of the regions with the rotated system //and the pt-weighted average of phi in the rotated system cenbroad_up += object_pt_in[i]*sqrt(pow(object_eta_in[i]-eta_up, 2) + pow(_delta_phi(phi_up[index_up], phi_up_aver), 2)); index_up += 1; } else { cenbroad_down += object_pt_in[i]*sqrt(pow(object_eta_in[i]-eta_down, 2)+ pow(_delta_phi(phi_down[index_down], phi_down_aver), 2)); index_down += 1; } } if (index_up == 0 || index_down ==0) _event_shapes[NEVTVAR] *= -1; cenbroad_up=cenbroad_up/(2*pt_sum_c); cenbroad_down=cenbroad_down/(2*pt_sum_c); //central total jet broadening double centotbroad = 0; centotbroad = cenbroad_up + cenbroad_down; _event_shapes[1] = centotbroad; for (int ij = 0; ij < 5; ij++) { if (_event_shapes[ij] < 1.e-20) _event_shapes[ij] = 1.e-20; _event_shapes[ij] = log(_event_shapes[ij]); } } return 1; } double EventShape::_three_jet_res(const vector& in_object_px, const vector& in_object_py, const vector& in_object_pz, const vector& in_object_e, int irap) { size_t y3_length = (size_t)in_object_px.size(); if (((size_t) in_object_py.size()!=y3_length) || ((size_t) in_object_pz.size()!=y3_length) || (in_object_e.size()!=y3_length)) { // cout << "ERROR!!!! Input vectors differ in size! Change that please!" << endl; // cout<<"py_size: "< in_object_p, in_object_pt, in_object_eta, in_object_phi; if (!in_object_p.empty()) { in_object_p.clear(); in_object_pt.clear(); in_object_eta.clear(); in_object_phi.clear(); } for (size_t j = 0; j < y3_length; j++) { in_object_p.push_back(0.); in_object_pt.push_back(0.); in_object_eta.push_back(0.); in_object_phi.push_back(0.); } double theta_y3_1st = 0; for (size_t k =0; k 1E-5) { theta_y3_1st = atan(in_object_pt[k]/(in_object_pz[k])); } else { theta_y3_1st = M_PI/2; } if (theta_y3_1st<0.) theta_y3_1st = theta_y3_1st + M_PI; in_object_eta[k] = - log(tan(0.5*theta_y3_1st)); } //calculates the real rapidity if (irap == 1) { in_object_eta[k]=0.5*log((in_object_e[k]+in_object_pz[k])/(in_object_e[k]-in_object_pz[k])); } in_object_phi[k] = atan2(in_object_py[k], in_object_px[k]); } //the three-jet resolution //threshold y3 double y3 = 0; //vector which will be filled with the //minimum of the distances double max_dmin_temp=0; double max_dmin = 0; //distance input object k, beam double distance_jB = 0; double distance_jB_min = 0; //distance of input object k to l double distance_jk = 0; double distance_jk_min = 0; //as we search the minimum of the distances //give them values which are for sure higher //than those we evaluate first in the for-loups size_t index_jB = 0; size_t index_j_jk = 0; size_t index_k_jk = 0; //to decide later if the minmum is a jB or jk int decide_jB = -1; vector input_pt, input_px, input_py, input_pz; vector input_p, input_e, input_phi, input_eta; if (!input_pt.empty()) { input_pt.clear(); input_px.clear(); input_px.clear(); input_pz.clear(); input_p.clear(); input_e.clear(); input_phi.clear(); input_eta.clear(); } for (size_t j = 0; j < y3_length; j++){ input_pt.push_back(in_object_pt[j]); input_px.push_back(in_object_px[j]); input_py.push_back(in_object_py[j]); input_pz.push_back(in_object_pz[j]); input_p.push_back(in_object_p[j]); input_e.push_back(in_object_e[j]); input_phi.push_back(in_object_phi[j]); input_eta.push_back(in_object_eta[j]); } if (y3_length<3) { return -1; } else { size_t rest = y3_length; for (size_t i = 0; i 2) { for (size_t j=0; j 0) { for(size_t k=0; k 1E-5){ theta_new = atan(input_pt[index_k_jk]/(input_pz[index_k_jk])); } else { theta_new = M_PI/2; } if (theta_new < 0) { theta_new = theta_new + M_PI; } input_eta[index_k_jk] = - log(tan(0.5*theta_new)); } //in the real rapidity y is wanted if (irap == 1) { input_eta[index_k_jk] = 0.5 * log((input_e[index_k_jk]+ input_pz[index_k_jk]) / (input_e[index_k_jk] - input_pz[index_k_jk])); } input_phi[index_k_jk] = atan2(input_py[index_k_jk], input_px[index_k_jk]); if (index_j_jk != rest-1) { for (size_t i = index_j_jk; i EventShape::_thrust(const vector& input_px, const vector& input_py) { double thrustmax_calc = 0; double temp_calc = 0; size_t length_thrust_calc = 0; vector thrust_values, thrust_axis_calc; vector p_thrust_max_calc, p_dec_1_calc, p_dec_2_calc, p_pt_beam_calc; if (!thrust_values.empty()){ thrust_values.clear(); thrust_axis_calc.clear(); p_thrust_max_calc.clear(); p_dec_1_calc.clear(); p_dec_2_calc.clear(); p_pt_beam_calc.clear(); } for (size_t j = 0; j < 3; j++){ p_pt_beam_calc.push_back(0.); p_dec_1_calc.push_back(0.); p_dec_2_calc.push_back(0.); p_thrust_max_calc.push_back(0.); thrust_axis_calc.push_back(0.); } for (size_t j = 0; j < 4; j++) { thrust_values.push_back(0.); } length_thrust_calc = input_px.size(); if (input_py.size() != length_thrust_calc) { /// @todo Change to exception or assert cout<<"ERROR in thrust calculation!!! Size of input vectors differs. Change that please!"<=0){ p_thrust_max_calc[0]= p_thrust_max_calc[0] + input_px[i]; p_thrust_max_calc[1]= p_thrust_max_calc[1] + input_py[i]; }else{ p_thrust_max_calc[0]= p_thrust_max_calc[0] - input_px[i]; p_thrust_max_calc[1]= p_thrust_max_calc[1] - input_py[i]; } } } p_dec_1_calc[0] = p_thrust_max_calc[0] + input_px[k]; p_dec_1_calc[1] = p_thrust_max_calc[1] + input_py[k]; p_dec_1_calc[2] = 0; p_dec_2_calc[0] = p_thrust_max_calc[0] - input_px[k]; p_dec_2_calc[1] = p_thrust_max_calc[1] - input_py[k]; p_dec_2_calc[2] = 0; temp_calc = pow(p_dec_1_calc[0], 2) + pow(p_dec_1_calc[1], 2); if (temp_calc>thrustmax_calc) { thrustmax_calc =temp_calc; for (size_t i=0; i<3; i++) { thrust_axis_calc[i] = p_dec_1_calc[i]/sqrt(thrustmax_calc); } } temp_calc = pow(p_dec_2_calc[0], 2)+pow(p_dec_2_calc[1], 2); if (temp_calc > thrustmax_calc) { thrustmax_calc =temp_calc; for (size_t i=0; i<3; i++) { thrust_axis_calc[i] = p_dec_2_calc[i]/sqrt(thrustmax_calc); } } } for (size_t j = 0; j < 3; j++) thrust_values[j] = thrust_axis_calc[j]; const double thrust_calc = sqrt(thrustmax_calc)/pt_sum_calc; // the variable which gets returned is not the thrust but tau=1-thrust thrust_values[3] = 1 - thrust_calc; if (thrust_values[3] < 1e-20) thrust_values[3] = 1e-20; return thrust_values; } double EventShape::_delta_phi(double phi1, double phi2) { double dphi = fabs(phi2 - phi1); if (dphi > M_PI) dphi = 2*M_PI - dphi; return dphi; } // Returns the scalar product between two 4 momenta double EventShape::_lorentz_sp(const vector& a, const vector& b) { size_t dim = (size_t) a.size(); if (a.size()!=b.size()) { cout<<"ERROR!!! Dimension of input vectors are different! Change that please!"<(event, "UFS"); double ancestor_lftime; foreach (const Particle& p, ufs.particles()) { id = p.pid(); if ((id != 310) && (id != -310)) continue; sumKs0_all ++; ancestor_lftime = 0.; const GenParticle* long_ancestor = getLongestLivedAncestor(p, ancestor_lftime); if ( !(long_ancestor) ) { sumKs0_badnull ++; continue; } if ( ancestor_lftime > MAX_CTAU ) { sumKs0_badlft ++; MSG_DEBUG("Ancestor " << long_ancestor->pdg_id() << ", ctau: " << ancestor_lftime << " [m]"); continue; } const FourMomentum& qmom = p.momentum(); y = 0.5 * log((qmom.E() + qmom.pz())/(qmom.E() - qmom.pz())); pT = sqrt((qmom.px() * qmom.px()) + (qmom.py() * qmom.py())); if (pT < MIN_PT) { sum_low_pt_loss ++; MSG_DEBUG("Small pT K^0_S: " << pT << " GeV/c."); } if (pT > 1.6) { sum_high_pt_loss ++; } if (y > 2.5 && y < 4.0) { _h_K0s_pt_y_all->fill(pT); if (y > 2.5 && y < 3.0) { _h_K0s_pt_y_30->fill(pT); _h_K0s_pt_30->fill(pT); sumKs0_30->fill(); } else if (y > 3.0 && y < 3.5) { _h_K0s_pt_y_35->fill(pT); _h_K0s_pt_35->fill(pT); sumKs0_35->fill(); } else if (y > 3.5 && y < 4.0) { _h_K0s_pt_y_40->fill(pT); _h_K0s_pt_40->fill(pT); sumKs0_40->fill(); } } else if (y < 2.5) { sumKs0_outdwn ++; } else if (y > 4.0) { sumKs0_outup ++; } } } /// Normalise histograms etc., after the run void finalize() { MSG_DEBUG("Total number Ks0: " << sumKs0_all << endl << "Sum of weights: " << sumOfWeights() << endl << "Weight Ks0 (2.5 < y < 3.0): " << double(sumKs0_30) << endl << "Weight Ks0 (3.0 < y < 3.5): " << double(sumKs0_35) << endl << "Weight Ks0 (3.5 < y < 4.0): " << double(sumKs0_40) << endl << "Nb. unprompt Ks0 [null mother]: " << sumKs0_badnull << endl << "Nb. unprompt Ks0 [mother lifetime exceeded]: " << sumKs0_badlft << endl << "Nb. Ks0 (y > 4.0): " << sumKs0_outup << endl << "Nb. Ks0 (y < 2.5): " << sumKs0_outdwn << endl << "Nb. Ks0 (pT < " << (MIN_PT/MeV) << " MeV/c): " << sum_low_pt_loss << endl << "Nb. Ks0 (pT > 1.6 GeV/c): " << sum_high_pt_loss << endl << "Cross-section [mb]: " << crossSection()/millibarn << endl << "Nb. events: " << numEvents()); // Compute cross-section; multiply by bin width for correct scaling // cross-section given by Rivet in pb double xsection_factor = crossSection()/sumOfWeights(); // Multiply bin width for correct scaling, xsection in mub scale(_h_K0s_pt_30, 0.2*xsection_factor/microbarn); scale(_h_K0s_pt_35, 0.2*xsection_factor/microbarn); scale(_h_K0s_pt_40, 0.2*xsection_factor/microbarn); // Divide by dy (rapidity window width), xsection in mb scale(_h_K0s_pt_y_30, xsection_factor/0.5/millibarn); scale(_h_K0s_pt_y_35, xsection_factor/0.5/millibarn); scale(_h_K0s_pt_y_40, xsection_factor/0.5/millibarn); scale(_h_K0s_pt_y_all, xsection_factor/1.5/millibarn); } //@} private: /// Get particle lifetime from hardcoded data double getLifeTime(int pid) { double lft = -1.0; if (pid < 0) pid = - pid; // Correct Pythia6 PIDs for f0(980), f0(1370) mesons if (pid == 10331) pid = 30221; if (pid == 10221) pid = 9010221; map::iterator pPartLft = partLftMap.find(pid); // search stable particle list if (pPartLft == partLftMap.end()) { - if (pid <= 100) return 0.0; + if (pid <= 100 || pid == 990) return 0.0; for ( auto id : stablePDGIds ) { if (pid == id) { lft = 0.0; break; } } } else { lft = (*pPartLft).second; } if (lft < 0.0) MSG_ERROR("Could not determine lifetime for particle with PID " << pid << "... This K_s^0 will be considered unprompt!"); return lft; } const GenParticle* getLongestLivedAncestor(const Particle& p, double& lifeTime) { const GenParticle* ret = nullptr; lifeTime = 1.; if (p.genParticle() == nullptr) return nullptr; const GenParticle* pmother = p.genParticle(); double longest_ctau = 0.; double mother_ctau; int mother_pid, n_inparts; const GenVertex* ivertex = pmother->production_vertex(); while (ivertex) { n_inparts = ivertex->particles_in_size(); if (n_inparts < 1) {ret = nullptr; break;} // error: should never happen! const auto iPart_invtx = ivertex->particles_in_const_begin(); pmother = (*iPart_invtx); // first mother particle mother_pid = pmother->pdg_id(); ivertex = pmother->production_vertex(); // get next vertex if ( (mother_pid == 2212) || (mother_pid <= 100) ) { if (ret == nullptr) ret = pmother; continue; } mother_ctau = getLifeTime(mother_pid); if (mother_ctau < 0.) { ret= nullptr; break; } // error:should never happen! if (mother_ctau > longest_ctau) { longest_ctau = mother_ctau; ret = pmother; } } if (ret) lifeTime = longest_ctau * c_light; return ret; } // Fill the PDG Id to Lifetime[seconds] map // Data was extract from LHCb Particle Table using ParticleSvc bool fillMap(map &m) { m[6] = 4.707703E-25; m[11] = 1.E+16; m[12] = 1.E+16; m[13] = 2.197019E-06; m[14] = 1.E+16; m[15] = 2.906E-13; m[16] = 1.E+16; m[22] = 1.E+16; m[23] = 2.637914E-25; m[24] = 3.075758E-25; m[25] = 9.4E-26; m[35] = 9.4E-26; m[36] = 9.4E-26; m[37] = 9.4E-26; m[84] = 3.335641E-13; m[85] = 1.290893E-12; m[111] = 8.4E-17; m[113] = 4.405704E-24; m[115] = 6.151516E-24; m[117] = 4.088275E-24; m[119] = 2.102914E-24; m[130] = 5.116E-08; m[150] = 1.525E-12; m[211] = 2.6033E-08; m[213] = 4.405704E-24; m[215] = 6.151516E-24; m[217] = 4.088275E-24; m[219] = 2.102914E-24; m[221] = 5.063171E-19; m[223] = 7.752794E-23; m[225] = 3.555982E-24; m[227] = 3.91793E-24; m[229] = 2.777267E-24; m[310] = 8.953E-11; m[313] = 1.308573E-23; m[315] = 6.038644E-24; m[317] = 4.139699E-24; m[319] = 3.324304E-24; m[321] = 1.238E-08; m[323] = 1.295693E-23; m[325] = 6.682357E-24; m[327] = 4.139699E-24; m[329] = 3.324304E-24; m[331] = 3.210791E-21; m[333] = 1.545099E-22; m[335] = 9.016605E-24; m[337] = 7.565657E-24; m[350] = 1.407125E-12; m[411] = 1.04E-12; m[413] = 6.856377E-21; m[415] = 1.778952E-23; m[421] = 4.101E-13; m[423] = 1.000003E-19; m[425] = 1.530726E-23; m[431] = 5.E-13; m[433] = 1.000003E-19; m[435] = 3.291061E-23; m[441] = 2.465214E-23; m[443] = 7.062363E-21; m[445] = 3.242425E-22; m[510] = 1.525E-12; m[511] = 1.525E-12; m[513] = 1.000019E-19; m[515] = 1.31E-23; m[521] = 1.638E-12; m[523] = 1.000019E-19; m[525] = 1.31E-23; m[530] = 1.536875E-12; m[531] = 1.472E-12; m[533] = 1.E-19; m[535] = 1.31E-23; m[541] = 4.5E-13; m[553] = 1.218911E-20; m[1112] = 4.539394E-24; m[1114] = 5.578069E-24; m[1116] = 1.994582E-24; m[1118] = 2.269697E-24; m[1212] = 4.539394E-24; m[1214] = 5.723584E-24; m[1216] = 1.994582E-24; m[1218] = 1.316424E-24; m[2112] = 8.857E+02; m[2114] = 5.578069E-24; m[2116] = 4.388081E-24; m[2118] = 2.269697E-24; m[2122] = 4.539394E-24; m[2124] = 5.723584E-24; m[2126] = 1.994582E-24; m[2128] = 1.316424E-24; m[2212] = 1.E+16; m[2214] = 5.578069E-24; m[2216] = 4.388081E-24; m[2218] = 2.269697E-24; m[2222] = 4.539394E-24; m[2224] = 5.578069E-24; m[2226] = 1.994582E-24; m[2228] = 2.269697E-24; m[3112] = 1.479E-10; m[3114] = 1.670589E-23; m[3116] = 5.485102E-24; m[3118] = 3.656734E-24; m[3122] = 2.631E-10; m[3124] = 4.219309E-23; m[3126] = 8.227653E-24; m[3128] = 3.291061E-24; m[3212] = 7.4E-20; m[3214] = 1.828367E-23; m[3216] = 5.485102E-24; m[3218] = 3.656734E-24; m[3222] = 8.018E-11; m[3224] = 1.838582E-23; m[3226] = 5.485102E-24; m[3228] = 3.656734E-24; m[3312] = 1.639E-10; m[3314] = 6.648608E-23; m[3322] = 2.9E-10; m[3324] = 7.233101E-23; m[3334] = 8.21E-11; m[4112] = 2.991874E-22; m[4114] = 4.088274E-23; m[4122] = 2.E-13; m[4132] = 1.12E-13; m[4212] = 3.999999E-22; m[4214] = 3.291061E-22; m[4222] = 2.951624E-22; m[4224] = 4.417531E-23; m[4232] = 4.42E-13; m[4332] = 6.9E-14; m[4412] = 3.335641E-13; m[4422] = 3.335641E-13; m[4432] = 3.335641E-13; m[5112] = 1.E-19; m[5122] = 1.38E-12; m[5132] = 1.42E-12; m[5142] = 1.290893E-12; m[5212] = 1.E-19; m[5222] = 1.E-19; m[5232] = 1.42E-12; m[5242] = 1.290893E-12; m[5312] = 1.E-19; m[5322] = 1.E-19; m[5332] = 1.55E-12; m[5342] = 1.290893E-12; m[5442] = 1.290893E-12; m[5512] = 1.290893E-12; m[5522] = 1.290893E-12; m[5532] = 1.290893E-12; m[5542] = 1.290893E-12; m[10111] = 2.48382E-24; m[10113] = 4.635297E-24; m[10115] = 2.54136E-24; m[10211] = 2.48382E-24; m[10213] = 4.635297E-24; m[10215] = 2.54136E-24; m[10223] = 1.828367E-24; m[10225] = 3.636531E-24; m[10311] = 2.437823E-24; m[10313] = 7.313469E-24; m[10315] = 3.538775E-24; m[10321] = 2.437823E-24; m[10323] = 7.313469E-24; m[10325] = 3.538775E-24; m[10331] = 4.804469E-24; m[10411] = 4.38E-24; m[10413] = 3.29E-23; m[10421] = 4.38E-24; m[10423] = 3.22653E-23; m[10431] = 6.5821E-22; m[10433] = 6.5821E-22; m[10441] = 6.453061E-23; m[10511] = 4.39E-24; m[10513] = 1.65E-23; m[10521] = 4.39E-24; m[10523] = 1.65E-23; m[10531] = 4.39E-24; m[10533] = 1.65E-23; m[11114] = 2.194041E-24; m[11116] = 1.828367E-24; m[11212] = 1.880606E-24; m[11216] = 1.828367E-24; m[12112] = 2.194041E-24; m[12114] = 2.194041E-24; m[12116] = 5.063171E-24; m[12126] = 1.828367E-24; m[12212] = 2.194041E-24; m[12214] = 2.194041E-24; m[12216] = 5.063171E-24; m[12224] = 2.194041E-24; m[12226] = 1.828367E-24; m[13112] = 6.582122E-24; m[13114] = 1.09702E-23; m[13116] = 5.485102E-24; m[13122] = 1.316424E-23; m[13124] = 1.09702E-23; m[13126] = 6.928549E-24; m[13212] = 6.582122E-24; m[13214] = 1.09702E-23; m[13216] = 5.485102E-24; m[13222] = 6.582122E-24; m[13224] = 1.09702E-23; m[13226] = 5.485102E-24; m[13312] = 4.135667E-22; m[13314] = 2.742551E-23; m[13324] = 2.742551E-23; m[14122] = 1.828367E-22; m[20022] = 1.E+16; m[20113] = 1.567172E-24; m[20213] = 1.567172E-24; m[20223] = 2.708692E-23; m[20313] = 3.782829E-24; m[20315] = 2.384827E-24; m[20323] = 3.782829E-24; m[20325] = 2.384827E-24; m[20333] = 1.198929E-23; m[20413] = 2.63E-24; m[20423] = 2.63E-24; m[20433] = 6.5821E-22; m[20443] = 7.395643E-22; m[20513] = 2.63E-24; m[20523] = 2.63E-24; m[20533] = 2.63E-24; m[21112] = 2.632849E-24; m[21114] = 3.291061E-24; m[21212] = 2.632849E-24; m[21214] = 6.582122E-24; m[22112] = 4.388081E-24; m[22114] = 3.291061E-24; m[22122] = 2.632849E-24; m[22124] = 6.582122E-24; m[22212] = 4.388081E-24; m[22214] = 3.291061E-24; m[22222] = 2.632849E-24; m[22224] = 3.291061E-24; m[23112] = 7.313469E-24; m[23114] = 2.991874E-24; m[23122] = 4.388081E-24; m[23124] = 6.582122E-24; m[23126] = 3.291061E-24; m[23212] = 7.313469E-24; m[23214] = 2.991874E-24; m[23222] = 7.313469E-24; m[23224] = 2.991874E-24; m[30113] = 2.632849E-24; m[30213] = 2.632849E-24; m[30221] = 1.880606E-24; m[30223] = 2.089563E-24; m[30313] = 2.056913E-24; m[30323] = 2.056913E-24; m[30443] = 2.419898E-23; m[31114] = 1.880606E-24; m[31214] = 3.291061E-24; m[32112] = 3.989164E-24; m[32114] = 1.880606E-24; m[32124] = 3.291061E-24; m[32212] = 3.989164E-24; m[32214] = 1.880606E-24; m[32224] = 1.880606E-24; m[33122] = 1.880606E-23; m[42112] = 6.582122E-24; m[42212] = 6.582122E-24; m[43122] = 2.194041E-24; m[53122] = 4.388081E-24; m[100111] = 1.645531E-24; m[100113] = 1.64553E-24; m[100211] = 1.645531E-24; m[100213] = 1.64553E-24; m[100221] = 1.196749E-23; m[100223] = 3.061452E-24; m[100313] = 2.837122E-24; m[100323] = 2.837122E-24; m[100331] = 4.459432E-25; m[100333] = 4.388081E-24; m[100441] = 4.701516E-23; m[100443] = 2.076379E-21; m[100553] = 2.056913E-20; m[200553] = 3.242425E-20; m[300553] = 3.210791E-23; m[9000111] = 8.776163E-24; m[9000211] = 8.776163E-24; m[9000443] = 8.227652E-24; m[9000553] = 5.983747E-24; m[9010111] = 3.164482E-24; m[9010211] = 3.164482E-24; m[9010221] = 9.403031E-24; m[9010443] = 8.438618E-24; m[9010553] = 8.3318E-24; m[9020221] = 8.093281E-23; m[9020443] = 1.061633E-23; m[9030221] = 6.038644E-24; m[9042413] = 2.07634E-21; m[9050225] = 1.394517E-24; m[9060225] = 3.291061E-24; m[9080225] = 4.388081E-24; m[9090225] = 2.056913E-24; m[9910445] = 2.07634E-21; m[9920443] = 2.07634E-21; return true; } /// @name Histograms //@{ Histo1DPtr _h_K0s_pt_y_30; // histogram for 2.5 < y < 3.0 (d2sigma) Histo1DPtr _h_K0s_pt_y_35; // histogram for 3.0 < y < 3.5 (d2sigma) Histo1DPtr _h_K0s_pt_y_40; // histogram for 3.5 < y < 4.0 (d2sigma) Histo1DPtr _h_K0s_pt_30; // histogram for 2.5 < y < 3.0 (sigma) Histo1DPtr _h_K0s_pt_35; // histogram for 3.0 < y < 3.5 (sigma) Histo1DPtr _h_K0s_pt_40; // histogram for 3.5 < y < 4.0 (sigma) Histo1DPtr _h_K0s_pt_y_all; // histogram for 2.5 < y < 4.0 (d2sigma) CounterPtr sumKs0_30; // Sum of weights 2.5 < y < 3.0 CounterPtr sumKs0_35; // Sum of weights 3.0 < y < 3.5 CounterPtr sumKs0_40; // Sum of weights 3.5 < y < 4.0 // Various counters mainly for debugging and comparisons between different generators size_t sumKs0_badnull; // Nb of particles for which mother could not be identified size_t sumKs0_badlft; // Nb of mesons with long lived mothers size_t sumKs0_all; // Nb of all Ks0 generated size_t sumKs0_outup; // Nb of mesons with y > 4.0 size_t sumKs0_outdwn; // Nb of mesons with y < 2.5 size_t sum_low_pt_loss; // Nb of mesons with very low pT (indicates when units are mixed-up) size_t sum_high_pt_loss; // Nb of mesons with pT > 1.6 GeV/c // Map between PDG id and particle lifetimes in seconds std::map partLftMap; // Set of PDG Ids for stable particles (PDG Id <= 100 are considered stable) static const array stablePDGIds; //@} }; // Actual initialization according to ISO C++ requirements const array LHCB_2010_S8758301::stablePDGIds{{ 311, 543, 545, 551, 555, 557, 1103, 2101, 2103, 2203, 3101, 3103, 3201, 3203, 3303, 4101, 4103, 4124, 4201, 4203, 4301, 4303, 4312, 4314, 4322, 4324, 4334, 4403, 4414, 4424, 4434, 4444, 5101, 5103, 5114, 5201, 5203, 5214, 5224, 5301, 5303, 5314, 5324, 5334, 5401, 5403, 5412, 5414, 5422, 5424, 5432, 5434, 5444, 5503, 5514, 5524, 5534, 5544, 5554, 10022, 10333, 10335, 10443, 10541, 10543, 10551, 10553, 10555, 11112, 12118, 12122, 12218, 12222, 13316, 13326, 20543, 20553, 20555, 23314, 23324, 30343, 30353, 30363, 30553, 33314, 33324, 41214, 42124, 52114, 52214, 100311, 100315, 100321, 100325, 100411, 100413, 100421, 100423, 100551, 100555, 100557, 110551, 110553, 110555, 120553, 120555, 130553, 200551, 200555, 210551, 210553, 220553, 1000001, 1000002, 1000003, 1000004, 1000005, 1000006, 1000011, 1000012, 1000013, 1000014, 1000015, 1000016, 1000021, 1000022, 1000023, 1000024, 1000025, 1000035, 1000037, 1000039, 2000001, 2000002, 2000003, 2000004, 2000005, 2000006, 2000011, 2000012, 2000013, 2000014, 2000015, 2000016, 3000111, 3000113, 3000211, 3000213, 3000221, 3000223, 3000331, 3100021, 3100111, 3100113, 3200111, 3200113, 3300113, 3400113, 4000001, 4000002, 4000011, 4000012, 5000039, 9000221, 9900012, 9900014, 9900016, 9900023, 9900024, 9900041, 9900042}}; // Hook for the plugin system DECLARE_RIVET_PLUGIN(LHCB_2010_S8758301); }