|
[Rivet-svn] r3450 - in trunk: . bin data/anainfo data/plotinfo data/refdata src/Analysesblackhole at projects.hepforge.org blackhole at projects.hepforge.orgWed Oct 19 17:16:37 BST 2011
Author: buckley Date: Wed Oct 19 17:16:36 2011 New Revision: 3450 Log: Adding new version of LHCB_2010_S8758301, based on submission from Alex Grecu. There is some slightly dodgy-looking GenParticle* fiddling going on, but apparently it's necessary (and hopefully robust). Modified: trunk/ChangeLog trunk/bin/rivet-mkhtml trunk/data/anainfo/LHCB_2010_S8758301.info trunk/data/plotinfo/LHCB_2010_S8758301.plot trunk/data/refdata/LHCB_2010_S8758301.aida trunk/src/Analyses/LHCB_2010_S8758301.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Tue Oct 18 11:23:27 2011 (r3449) +++ trunk/ChangeLog Wed Oct 19 17:16:36 2011 (r3450) @@ -1,4 +1,16 @@ +2011-10-19 Andy Buckley <andy at insectnation.org> + + * Adding new version of LHCB_2010_S8758301, based on submission + from Alex Grecu. There is some slightly dodgy-looking GenParticle* + fiddling going on, but apparently it's necessary (and hopefully robust). + +2011-10-17 Andy Buckley <andy at insectnation.org> + + * bin/rivet-nopy linker line tweak to make compilation work with + GCC 4.6 (-lHepMC has to be explicitly added for some reason). + 2011-10-13 Frank Siegert <frank.siegert at cern.ch> + * Add four CMS QCD analyses kindly provided by CMS. 2011-10-12 Andy Buckley <andy at insectnation.org> Modified: trunk/bin/rivet-mkhtml ============================================================================== --- trunk/bin/rivet-mkhtml Tue Oct 18 11:23:27 2011 (r3449) +++ trunk/bin/rivet-mkhtml Wed Oct 19 17:16:36 2011 (r3450) @@ -186,7 +186,7 @@ index = open(os.path.join(opts.OUTPUTDIR, "index.html"), "w") index.write('<html>\n<head>\n<title>%s</title>\n%s</head>\n<body>' % (opts.TITLE, style)) -if opts.BOOKLET and opts.VECTORFORMAT=="PDF": +if opts.BOOKLET and opts.VECTORFORMAT == "PDF": index.write('<h2><a href="booklet.pdf">%s</a></h2>\n\n' % opts.TITLE) else: index.write('<h2>%s</h2>\n\n' % opts.TITLE) Modified: trunk/data/anainfo/LHCB_2010_S8758301.info ============================================================================== --- trunk/data/anainfo/LHCB_2010_S8758301.info Tue Oct 18 11:23:27 2011 (r3449) +++ trunk/data/anainfo/LHCB_2010_S8758301.info Wed Oct 19 17:16:36 2011 (r3450) @@ -1,13 +1,14 @@ Name: LHCB_2010_S8758301 Year: 2010 -Summary: Differential cross section measurement of $K^0_S$ production in three - rapidity windows at $\sqrt{s}=\unit{0.9}{\TeV}$ +Summary: LHCb differential cross section measurement of prompt $K^0_S$ production in three + rapidity windows at $\sqrt{s}$=0.9 TeV Experiment: LHCB -Collider: LHC +Collider: LHC 900 GeV SpiresID: 8758301 -Status: UNVALIDATED +Status: VALIDATED Authors: - - Holger Schulz holger.schulz at physik.hu-berlin.de + - Holger Schulz <holger.schulz at physik.hu-berlin.de> + - Alex Grecu <Alex.Grecu at cern.ch> References: - Phys.Lett.B693:69-80,2010 - arXiv:1008.3105[hep-ex] @@ -16,12 +17,14 @@ NumEvents: 1000000 Beams: [p+, p+] Energies: [900] -PtCuts: +PtCuts: [0, 1.6] Description: - 'Differential cross-section measurement of prompt $K^0_S$ production in pp collisions - at $\sqrt{s}=\unit{0.9}{\TeV}$ in the rapidity windows $2.5<y<3.0$, $3.0<y<3.5$ and $3.5<y<4.0$. - Unresolved issues at the moment are the normalisations. Not clear to me how to do this from the - paper so I normalise to reference data area.' + The paper presents the cross-section and double differential cross-section measurement + for prompt $K^0_S$ production in pp collisions at $\sqrt{s}$=0.9 TeV. The data were taken + during the LHCb run in December 2009 and cover a transversal momentum range from 0 to 1.6 GeV/c. + The differential production cross-section is calculated for three rapidity windows + $2.5 < y < 3.0$, $3.0 < y < 3.5$ and $3.5 < y < 4.0$ as well as the whole rapidity domain + $2.5 < y < 4.0$. BibKey: Aaij:2010nx BibTeX: '@Article{Aaij:2010nx, author = "Aaij, R and others", @@ -38,5 +41,3 @@ doi = "10.1016/j.physletb.2010.08.055", SLACcitation = "%%CITATION = 1008.3105;%%" }' -ToDo: - - Do not normalise to data! Modified: trunk/data/plotinfo/LHCB_2010_S8758301.plot ============================================================================== --- trunk/data/plotinfo/LHCB_2010_S8758301.plot Tue Oct 18 11:23:27 2011 (r3449) +++ trunk/data/plotinfo/LHCB_2010_S8758301.plot Wed Oct 19 17:16:36 2011 (r3450) @@ -1,17 +1,41 @@ # BEGIN PLOT /LHCB_2010_S8758301/d01-x01-y01 -Title=Diff. cross-section for $K_S^0$ in production ($2.5<y<3.0$) -XLabel=$p_\perp$ [GeV] -YLabel=$\text{d}^2\sigma(pp\rightarrow\text{K}_S^0 X) / (\text{d}y\text{d}p_\perp)$ +Title=Cross-section for prompt $K_S^0$ production ($2.5<y<3.0$) +XLabel=$p_\perp$ [GeV/c] +YLabel=$\sigma(pp\rightarrow\text{K}_S^0 X) [\mu b]$ # END PLOT # BEGIN PLOT /LHCB_2010_S8758301/d01-x01-y02 -Title=Diff. cross-section for $K_S^0$ in production ($3.0<y<3.5$) -XLabel=$p_\perp$ [GeV] -YLabel=$\text{d}^2\sigma(pp\rightarrow\text{K}_S^0 X) / (\text{d}y\text{d}p_\perp)$ +Title=Cross-section for prompt $K_S^0$ production ($3.0<y<3.5$) +XLabel=$p_\perp$ [GeV/c] +YLabel=$\sigma(pp\rightarrow\text{K}_S^0 X) [\mu b]$ # END PLOT # BEGIN PLOT /LHCB_2010_S8758301/d01-x01-y03 -Title=Diff. cross-section for $K_S^0$ in production ($3.5<y<4.0$) -XLabel=$p_\perp$ [GeV] -YLabel=$\text{d}^2\sigma(pp\rightarrow\text{K}_S^0 X) / (\text{d}y\text{d}p_\perp)$ +Title=Cross-section for prompt $K_S^0$ production ($3.5<y<4.0$) +XLabel=$p_\perp$ [GeV/c] +YLabel=$\sigma(pp\rightarrow\text{K}_S^0 X) [\mu b]$ +# END PLOT + +# BEGIN PLOT /LHCB_2010_S8758301/d02-x01-y01 +Title=Diff. cross-section for prompt $K_S^0$ production ($2.5<y<3.0$) +XLabel=$p_\perp$ [GeV/c] +YLabel=$\text{d}^2\sigma(pp\rightarrow\text{K}_S^0 X) / (\text{d}y\text{d}p_\perp) [mb/(GeV/c)]$ +# END PLOT + +# BEGIN PLOT /LHCB_2010_S8758301/d02-x01-y02 +Title=Diff. cross-section for prompt $K_S^0$ production ($3.0<y<3.5$) +XLabel=$p_\perp$ [GeV/c] +YLabel=$\text{d}^2\sigma(pp\rightarrow\text{K}_S^0 X) / (\text{d}y\text{d}p_\perp) [mb/(GeV/c)]$ +# END PLOT + +# BEGIN PLOT /LHCB_2010_S8758301/d02-x01-y03 +Title=Diff. cross-section for prompt $K_S^0$ production ($3.5<y<4.0$) +XLabel=$p_\perp$ [GeV/c] +YLabel=$\text{d}^2\sigma(pp\rightarrow\text{K}_S^0 X) / (\text{d}y\text{d}p_\perp) [mb/(GeV/c)]$ +# END PLOT + +# BEGIN PLOT /LHCB_2010_S8758301/d03-x01-y01 +Title=Diff. cross-section for prompt $K_S^0$ production ($2.5<y<4.0$) +XLabel=$p_\perp$ [GeV/c] +YLabel=$\text{d}^2\sigma(pp\rightarrow\text{K}_S^0 X) / (\text{d}y\text{d}p_\perp) [mb/(GeV/c)]$ # END PLOT Modified: trunk/data/refdata/LHCB_2010_S8758301.aida ============================================================================== --- trunk/data/refdata/LHCB_2010_S8758301.aida Tue Oct 18 11:23:27 2011 (r3449) +++ trunk/data/refdata/LHCB_2010_S8758301.aida Wed Oct 19 17:16:36 2011 (r3450) @@ -104,4 +104,140 @@ <measurement value="23.0" errorPlus="7.483314773547883" errorMinus="7.483314773547883"/> </dataPoint> </dataPointSet> + <dataPointSet name="d02-x01-y01" dimension="2" path="/REF/LHCB_2010_S8758301" title="D2(SIG)/DPT/DYRAP IN MB/GEV" > + <dataPoint> + <measurement value="0.1" errorPlus="0.1" errorMinus="0.1"/> + <measurement value="2.935" errorPlus="1.26387697186079" errorMinus="1.26387697186079"/> + </dataPoint> + <dataPoint> + <measurement value="0.30000000000000004" errorPlus="0.09999999999999998" errorMinus="0.10000000000000003"/> + <measurement value="6.485" errorPlus="2.6359652880870796" errorMinus="2.6359652880870796"/> + </dataPoint> + <dataPoint> + <measurement value="0.5" errorPlus="0.09999999999999998" errorMinus="0.09999999999999998"/> + <measurement value="6.178" errorPlus="1.338118455145134" errorMinus="1.338118455145134"/> + </dataPoint> + <dataPoint> + <measurement value="0.7" errorPlus="0.10000000000000009" errorMinus="0.09999999999999998"/> + <measurement value="4.009" errorPlus="0.7380243898408778" errorMinus="0.7380243898408778"/> + </dataPoint> + <dataPoint> + <measurement value="0.9" errorPlus="0.09999999999999998" errorMinus="0.09999999999999998"/> + <measurement value="2.324" errorPlus="0.4262980178232125" errorMinus="0.4262980178232125"/> + </dataPoint> + <dataPoint> + <measurement value="1.1" errorPlus="0.09999999999999987" errorMinus="0.10000000000000009"/> + <measurement value="1.146" errorPlus="0.21797476918212347" errorMinus="0.21797476918212347"/> + </dataPoint> + <dataPoint> + <measurement value="1.2999999999999998" errorPlus="0.10000000000000009" errorMinus="0.09999999999999987"/> + <measurement value="0.851" errorPlus="0.17204650534085256" errorMinus="0.17204650534085256"/> + </dataPoint> + <dataPoint> + <measurement value="1.5" errorPlus="0.10000000000000009" errorMinus="0.10000000000000009"/> + <measurement value="0.462" errorPlus="0.10217631819555842" errorMinus="0.10217631819555842"/> + </dataPoint> + </dataPointSet> + <dataPointSet name="d02-x01-y02" dimension="2" path="/REF/LHCB_2010_S8758301" title="D2(SIG)/DPT/DYRAP IN MB/GEV" > + <dataPoint> + <measurement value="0.1" errorPlus="0.1" errorMinus="0.1"/> + <measurement value="3.165" errorPlus="0.9499373663563297" errorMinus="0.9499373663563297"/> + </dataPoint> + <dataPoint> + <measurement value="0.30000000000000004" errorPlus="0.09999999999999998" errorMinus="0.10000000000000003"/> + <measurement value="5.619" errorPlus="1.1197861402964406" errorMinus="1.1197861402964406"/> + </dataPoint> + <dataPoint> + <measurement value="0.5" errorPlus="0.09999999999999998" errorMinus="0.09999999999999998"/> + <measurement value="5.335" errorPlus="0.9142740289431828" errorMinus="0.9142740289431828"/> + </dataPoint> + <dataPoint> + <measurement value="0.7" errorPlus="0.10000000000000009" errorMinus="0.09999999999999998"/> + <measurement value="3.709" errorPlus="0.6382993028352765" errorMinus="0.6382993028352765"/> + </dataPoint> + <dataPoint> + <measurement value="0.9" errorPlus="0.09999999999999998" errorMinus="0.09999999999999998"/> + <measurement value="1.832" errorPlus="0.32930988445535614" errorMinus="0.32930988445535614"/> + </dataPoint> + <dataPoint> + <measurement value="1.1" errorPlus="0.09999999999999987" errorMinus="0.10000000000000009"/> + <measurement value="1.362" errorPlus="0.24686230980042295" errorMinus="0.24686230980042295"/> + </dataPoint> + <dataPoint> + <measurement value="1.2999999999999998" errorPlus="0.10000000000000009" errorMinus="0.09999999999999987"/> + <measurement value="0.699" errorPlus="0.13813037319865606" errorMinus="0.13813037319865606"/> + </dataPoint> + <dataPoint> + <measurement value="1.5" errorPlus="0.10000000000000009" errorMinus="0.10000000000000009"/> + <measurement value="0.491" errorPlus="0.1036001930500132" errorMinus="0.1036001930500132"/> + </dataPoint> + </dataPointSet> + <dataPointSet name="d02-x01-y03" dimension="2" path="/REF/LHCB_2010_S8758301" title="D2(SIG)/DPT/DYRAP IN MB/GEV" > + <dataPoint> + <measurement value="0.1" errorPlus="0.1" errorMinus="0.1"/> + <measurement value="1.96" errorPlus="0.6719181497771883" errorMinus="0.6719181497771883"/> + </dataPoint> + <dataPoint> + <measurement value="0.30000000000000004" errorPlus="0.09999999999999998" errorMinus="0.10000000000000003"/> + <measurement value="5.71" errorPlus="1.2361035555324642" errorMinus="1.2361035555324642"/> + </dataPoint> + <dataPoint> + <measurement value="0.5" errorPlus="0.09999999999999998" errorMinus="0.09999999999999998"/> + <measurement value="4.772" errorPlus="0.8199786582588598" errorMinus="0.8199786582588598"/> + </dataPoint> + <dataPoint> + <measurement value="0.7" errorPlus="0.10000000000000009" errorMinus="0.09999999999999998"/> + <measurement value="3.226" errorPlus="0.5613777694209132" errorMinus="0.5613777694209132"/> + </dataPoint> + <dataPoint> + <measurement value="0.9" errorPlus="0.09999999999999998" errorMinus="0.09999999999999998"/> + <measurement value="2.015" errorPlus="0.3681548587211637" errorMinus="0.3681548587211637"/> + </dataPoint> + <dataPoint> + <measurement value="1.1" errorPlus="0.09999999999999987" errorMinus="0.10000000000000009"/> + <measurement value="1.078" errorPlus="0.21414014102918677" errorMinus="0.21414014102918677"/> + </dataPoint> + <dataPoint> + <measurement value="1.2999999999999998" errorPlus="0.10000000000000009" errorMinus="0.09999999999999987"/> + <measurement value="0.346" errorPlus="0.10625441167311596" errorMinus="0.10625441167311596"/> + </dataPoint> + <dataPoint> + <measurement value="1.5" errorPlus="0.10000000000000009" errorMinus="0.10000000000000009"/> + <measurement value="0.234" errorPlus="0.07300684899377592" errorMinus="0.07300684899377592"/> + </dataPoint> + </dataPointSet> + <dataPointSet name="d03-x01-y01" dimension="2" path="/REF/LHCB_2010_S8758301" title="D2(SIG)/DPT/DYRAP IN MB/GEV" > + <dataPoint> + <measurement value="0.1" errorPlus="0.1" errorMinus="0.1"/> + <measurement value="2.687" errorPlus="0.7770263830784641" errorMinus="0.7770263830784641"/> + </dataPoint> + <dataPoint> + <measurement value="0.30000000000000004" errorPlus="0.09999999999999998" errorMinus="0.10000000000000003"/> + <measurement value="5.938" errorPlus="1.4693515576607252" errorMinus="1.4693515576607252"/> + </dataPoint> + <dataPoint> + <measurement value="0.5" errorPlus="0.09999999999999998" errorMinus="0.09999999999999998"/> + <measurement value="5.428" errorPlus="0.928724393994257" errorMinus="0.928724393994257"/> + </dataPoint> + <dataPoint> + <measurement value="0.7" errorPlus="0.10000000000000009" errorMinus="0.09999999999999998"/> + <measurement value="3.648" errorPlus="0.6046436967338699" errorMinus="0.6046436967338699"/> + </dataPoint> + <dataPoint> + <measurement value="0.9" errorPlus="0.09999999999999998" errorMinus="0.09999999999999998"/> + <measurement value="2.057" errorPlus="0.3465616828213991" errorMinus="0.3465616828213991"/> + </dataPoint> + <dataPoint> + <measurement value="1.1" errorPlus="0.09999999999999987" errorMinus="0.10000000000000009"/> + <measurement value="1.195" errorPlus="0.20270421801235414" errorMinus="0.20270421801235414"/> + </dataPoint> + <dataPoint> + <measurement value="1.2999999999999998" errorPlus="0.10000000000000009" errorMinus="0.09999999999999987"/> + <measurement value="0.632" errorPlus="0.11403946685248927" errorMinus="0.11403946685248927"/> + </dataPoint> + <dataPoint> + <measurement value="1.5" errorPlus="0.10000000000000009" errorMinus="0.10000000000000009"/> + <measurement value="0.396" errorPlus="0.07430343195303969" errorMinus="0.07430343195303969"/> + </dataPoint> + </dataPointSet> </aida> Modified: trunk/src/Analyses/LHCB_2010_S8758301.cc ============================================================================== --- trunk/src/Analyses/LHCB_2010_S8758301.cc Tue Oct 18 11:23:27 2011 (r3449) +++ trunk/src/Analyses/LHCB_2010_S8758301.cc Wed Oct 19 17:16:36 2011 (r3450) @@ -3,82 +3,139 @@ #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/UnstableFinalState.hh" -#include "Rivet/Tools/ParticleIdUtils.hh" +#include "Rivet/Math/Constants.hh" +#include "Rivet/Math/Units.hh" +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" +#include "HepMC/GenVertex.h" +#include "HepMC/SimpleVector.h" namespace Rivet { +using namespace HepMC; +using namespace std; + + /// lifetime cut: longest living ancestor ctau < 10^-11 [m] + const double MAX_CTAU = 1.0E-11; // [m] + const double MIN_pT = 0.0001; // [GeV/c] class LHCB_2010_S8758301 : public Analysis { public: - /// @name Constructors etc. //@{ - /// Constructor LHCB_2010_S8758301() : Analysis("LHCB_2010_S8758301"), - sumKs0_30(0.0), sumKs0_35(0.0), sumKs0_40(0.0) - { } - + sumKs0_30(0.0), sumKs0_35(0.0), + sumKs0_40(0.0), sumKs0_badnull(0), + sumKs0_badlft(0), sumKs0_all(0), + sumKs0_outup(0), sumKs0_outdwn(0), + sum_low_pt_loss(0), sum_high_pt_loss(0) + { + setNeedsCrossSection(true); + setBeams(PROTON, PROTON); + } //@} - - - public: - /// @name Analysis methods //@{ - /// Book histograms and initialise projections before the run void init() { - + MSG_DEBUG("Initializing analysis!"); + fillMap(partLftMap); + _h_K0s_pt_30 = bookHistogram1D(1,1,1); + _h_K0s_pt_35 = bookHistogram1D(1,1,2); + _h_K0s_pt_40 = bookHistogram1D(1,1,3); + _h_K0s_pt_y_30 = bookHistogram1D(2,1,1); + _h_K0s_pt_y_35 = bookHistogram1D(2,1,2); + _h_K0s_pt_y_40 = bookHistogram1D(2,1,3); + _h_K0s_pt_y_all = bookHistogram1D(3,1,1); addProjection(UnstableFinalState(), "UFS"); - - _h_K0s_pt_y_30 = bookHistogram1D(1,1,1); - _h_K0s_pt_y_35 = bookHistogram1D(1,1,2); - _h_K0s_pt_y_40 = bookHistogram1D(1,1,3); - } - /// Perform the per-event analysis void analyze(const Event& event) { + int id; + double y, pT; const double weight = event.weight(); const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(event, "UFS"); - + double ancestor_lftime; foreach (const Particle& p, ufs.particles()) { - const PdgId id = abs(p.pdgId()); - - if (id == 310) { - double y = p.momentum().rapidity(); - double pT = p.momentum().perp(); + id = p.pdgId(); + 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, weight); if (y > 2.5 && y < 3.0) { - //_h_K0s_pt_y_30->fill(pT, weight/(0.5*0.2)); _h_K0s_pt_y_30->fill(pT, weight); + _h_K0s_pt_30->fill(pT, weight); sumKs0_30 += weight; - } - else if (y > 3.0 && y < 3.5) { - //_h_K0s_pt_y_35->fill(pT, weight/(0.5*0.2)); + } else if (y > 3.0 && y < 3.5) { _h_K0s_pt_y_35->fill(pT, weight); + _h_K0s_pt_35->fill(pT, weight); sumKs0_35 += weight; - } - else if (y > 3.5 && y < 4.0) { - //_h_K0s_pt_y_40->fill(pT, weight/(0.5*0.2)); + } else if (y > 3.5 && y < 4.0) { _h_K0s_pt_y_40->fill(pT, weight); + _h_K0s_pt_40->fill(pT, weight); sumKs0_40 += weight; } + } else if (y < 2.5) { + sumKs0_outdwn ++; + } else if (y > 4.0) { + sumKs0_outup ++; } } - - } /// Normalise histograms etc., after the run void finalize() { - /// @todo We have to normalise to reference data... for now - scale(_h_K0s_pt_y_30, 4.880000e+02/sumKs0_30); - scale(_h_K0s_pt_y_35, 4.442000e+02/sumKs0_35); - scale(_h_K0s_pt_y_40, 3.868000e+02/sumKs0_40); + MSG_DEBUG("Total number Ks0: " << sumKs0_all << endl + << "Sum of weights: " << sumOfWeights() << endl + << "Weight Ks0 (2.5 < y < 3.0): " << sumKs0_30 << endl + << "Weight Ks0 (3.0 < y < 3.5): " << sumKs0_35 << endl + << "Weight Ks0 (3.5 < y < 4.0): " << 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); } //@} @@ -86,25 +143,197 @@ private: - /// @name Histograms - //@{ - - AIDA::IHistogram1D *_h_K0s_pt_y_30; - AIDA::IHistogram1D *_h_K0s_pt_y_35; - AIDA::IHistogram1D *_h_K0s_pt_y_40; - - double sumKs0_30; - double sumKs0_35; - double sumKs0_40; + /// 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<int, double>::iterator pPartLft = partLftMap.find(pid); + // search stable particle list + if (pPartLft == partLftMap.end()) { + if (pid <= 100) return 0.0; + for (unsigned int i=0; i < sizeof(stablePDGIds)/sizeof(unsigned int); i++ ) { + if (pid == stablePDGIds[i]) { 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 = NULL; + lifeTime = 1.; + if (!p.hasGenParticle()) return NULL; + const GenParticle* pmother = &(p.genParticle()); + double longest_ctau = 0.; + double mother_ctau; + int mother_pid, n_inparts; + GenVertex* ivertex = pmother->production_vertex(); + while (ivertex) { + n_inparts = ivertex->particles_in_size(); + if (n_inparts < 1) {ret = NULL; break;}; // error: should never happen! + const HepMC::GenVertex::particles_in_const_iterator 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 == NULL) ret = pmother; + continue; + } + mother_ctau = getLifeTime(mother_pid); + if (mother_ctau < 0.) { ret= NULL; 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<int, double> &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[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[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 + //@{ + AIDA::IHistogram1D *_h_K0s_pt_y_30; // histogram for 2.5 < y < 3.0 (d2sigma) + AIDA::IHistogram1D *_h_K0s_pt_y_35; // histogram for 3.0 < y < 3.5 (d2sigma) + AIDA::IHistogram1D *_h_K0s_pt_y_40; // histogram for 3.5 < y < 4.0 (d2sigma) + AIDA::IHistogram1D *_h_K0s_pt_30; // histogram for 2.5 < y < 3.0 (sigma) + AIDA::IHistogram1D *_h_K0s_pt_35; // histogram for 3.0 < y < 3.5 (sigma) + AIDA::IHistogram1D *_h_K0s_pt_40; // histogram for 3.5 < y < 4.0 (sigma) + AIDA::IHistogram1D *_h_K0s_pt_y_all; // histogram for 2.5 < y < 4.0 (d2sigma) + double sumKs0_30; // Sum of weights 2.5 < y < 3.0 + double sumKs0_35; // Sum of weights 3.0 < y < 3.5 + double 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<int, double> partLftMap; + // Set of PDG Ids for stable particles (PDG Id <= 100 are considered stable) + static const int stablePDGIds[205]; + //@} }; + // Actual initialization according to ISO C++ requirements + const int LHCB_2010_S8758301::stablePDGIds[205] = { + 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}; + - // The hook for the plugin system + // Hook for the plugin system DECLARE_RIVET_PLUGIN(LHCB_2010_S8758301); }
More information about the Rivet-svn mailing list |