rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2014_I1282447

W + charm production at 7 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1282447
Status: VALIDATED
Authors:
  • Kristin Lohwasser
References: Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • W+c and inclusive double production in pp collisions

This routine implements the measurement of $W$ boson production in association with a single charm quark with $4.6 \text{fb}^{-1}$ of $pp$ collision data at $\sqrt{s} = 7$ TeV collected with the ATLAS detector at the Large Hadron Collider. Results are quoted for the $W$ boson production decaying into either an electron or muon and the charm quark being identified as any charmed hadron above 5 GeV inside $\Delta R = 0.4$ of a jet with more than 25 GeV. Alternatively the presence of the charm quark is indicated by the presence of a $D$ or a $D^*$ meson with a $p_\text{T}$ above 8 GeV. The cross sections are quoted for the number of opposite sign minus same sign events, where the signs under consideration are the charge of the $W$ boson and the charmed hadrons tagging the event. Given are the integrated and differential cross sections as a function of the pseudorapidity of the lepton from the $W$-boson decay are measured. Additionally, the cross-section ratios $\sigma(W^+ + bar{c}) / \sigma(W^- + c)$ as well as the $p_\text{T}$ dependent cross sections of the $D$/$D^*$ mesons normalized to the $W$ inclusive cross sections are published. The measured data is unfolded to the Born level. One should therefore take care to run on samples without QED radiation off of the electrons. IMPORTANT NOTICE --- For the MC predictions in the paper, the branching fractions to $D$ and $D^*$ mesons have been adapted to be 0.2256 ($D$) and 0.2287 ($D^*$) (LEP/HERA combination). Some suggestion on how to do post-processing -- in case separate samples for $W^+ + c$ and $W^- + c$ and inclusive $W$ are used -- are included in the cc-file. This post-processing is needed to properly display the histograms for the cross section ratio plots $\sigma(W^+ + bar{c}) / \sigma(W^- + c)$ as well as for the cross section ratios of $W + D^{(*)}$ production over inclusive $W$ production ($\sigma(W^{+/-} D^{(*)}) / \sigma(W^{+/-})$) as a function of the $D^{(*)}$ meson transverse momentum.

Source code: ATLAS_2014_I1282447.cc
  1// -*- C++ -*-
  2// ATLAS W+c analysis
  3
  4//////////////////////////////////////////////////////////////////////////
  5/*
  6  Description of rivet analysis ATLAS_2014_I1282447 W+c production
  7
  8  This rivet routine implements the ATLAS W+c analysis.
  9  Apart from those histograms, described and published on HEP Data, here
 10  are some helper histograms defined, these are:
 11
 12  d02-x01-y01, d02-x01-y02 and d08-x01-y01 are ratios, the nominator ("_plus")
 13  and denominator ("_minus") histograms are also given, so that the ratios can
 14  be reconstructed if need be (e.g. when running on separate samples).
 15
 16  d05 and d06 are ratios over inclusive W production.
 17  The routine has to be run on a sample for inclusive W production in order to
 18  make sure the denominator ("_winc") is correctly filled.
 19
 20  The ratios can be constructed using the following sample code:
 21  python divideWCharm.py
 22
 23  import yoda
 24  hists_wc   = yoda.read("Rivet_Wc.yoda")
 25  hists_winc = yoda.read("Rivet_Winc.yoda")
 26
 27  ## division histograms --> ONLY for different plus minus runs
 28  # (merge before using yodamerge Rivet_plus.yoda Rivet_minus.yoda > Rivet_Wc.yoda)
 29
 30  d02y01_plus = hists_wc["/ATLAS_2014_I1282447/d02-x01-y01_plus"]
 31  d02y01_minus = hists_wc["/ATLAS_2014_I1282447/d02-x01-y01_minus"]
 32  ratio_d02y01 =  d02y01_plus.divide(d02y01_minus)
 33  ratio_d02y01.path = "/ATLAS_2014_I1282447/d02-x01-y01"
 34
 35  d02y02_plus = hists_wc["/ATLAS_2014_I1282447/d02-x01-y02_plus"]
 36  d02y02_minus = hists_wc["/ATLAS_2014_I1282447/d02-x01-y02_minus"]
 37  ratio_d02y02=  d02y02_plus.divide(d02y02_minus)
 38  ratio_d02y02.path = "/ATLAS_2014_I1282447/d02-x01-y02"
 39
 40  d08y01_plus = hists_wc["/ATLAS_2014_I1282447/d08-x01-y01_plus"]
 41  d08y01_minus = hists_wc["/ATLAS_2014_I1282447/d08-x01-y01_minus"]
 42  ratio_d08y01=  d08y01_plus.divide(d08y01_minus)
 43  ratio_d08y01.path = "/ATLAS_2014_I1282447/d08-x01-y01"
 44
 45  # inclusive cross section
 46  h_winc = hists_winc["/ATLAS_2014_I1282447/d05-x01-y01"]
 47  h_d    = hists_wc["/ATLAS_2014_I1282447/d01-x01-y02"]
 48  h_dstar= hists_wc["/ATLAS_2014_I1282447/d01-x01-y03"]
 49
 50  ratio_wd      =  h_d.divide(h_winc)
 51  ratio_wd.path = "/ATLAS_2014_I1282447/d05-x01-y02"
 52
 53  ratio_wdstar      =  h_d.divide(h_winc)
 54  ratio_wdstar.path = "/ATLAS_2014_I1282447/d05-x01-y03"
 55
 56  # pT differential
 57  h_winc_plus  = hists_winc["/ATLAS_2014_I1282447/d06-x01-y01_winc"]
 58  h_winc_minus = hists_winc["/ATLAS_2014_I1282447/d06-x01-y02_winc"]
 59
 60  h_wd_plus      = hists_wc["/ATLAS_2014_I1282447/d06-x01-y01_wplus"]
 61  h_wd_minus     = hists_wc["/ATLAS_2014_I1282447/d06-x01-y02_wminus"]
 62  h_wdstar_plus  = hists_wc["/ATLAS_2014_I1282447/d06-x01-y03_wplus"]
 63  h_wdstar_minus = hists_wc["/ATLAS_2014_I1282447/d06-x01-y04_wminus"]
 64
 65  ratio_wd_plus       =  h_wd_plus.divide(h_winc_plus)
 66  ratio_wd_plus.path  = "/ATLAS_2014_I1282447/d06-x01-y01"
 67  ratio_wd_minus      =  h_wd_plus.divide(h_winc_minus)
 68  ratio_wd_minus.path = "/ATLAS_2014_I1282447/d06-x01-y02"
 69
 70  ratio_wdstar_plus       =  h_wdstar_plus.divide(h_winc_plus)
 71  ratio_wdstar_plus.path  = "/ATLAS_2014_I1282447/d06-x01-y03"
 72  ratio_wdstar_minus      =  h_wdstar_plus.divide(h_winc_minus)
 73  ratio_wdstar_minus.path = "/ATLAS_2014_I1282447/d06-x01-y04"
 74
 75  ratio_wd_plus =  h_wd_plus.divide(h_winc_plus)
 76  ratio_wd_plus.path = "/ATLAS_2014_I1282447/d06-x01-y01"
 77  ratio_wd_minus =  h_wd_plus.divide(h_winc_minus)
 78  ratio_wd_minus.path = "/ATLAS_2014_I1282447/d06-x01-y02"
 79
 80  h_winc_plus= hists_winc["/ATLAS_2014_I1282447/d06-x01-y01_winc"]
 81  h_winc_minus= hists_winc["/ATLAS_2014_I1282447/d06-x01-y02_winc"]
 82
 83  ## copy other histograms for plotting
 84
 85  d01x01y01= hists_wc["/ATLAS_2014_I1282447/d01-x01-y01"]
 86  d01x01y01.path = "/ATLAS_2014_I1282447/d01-x01-y01"
 87
 88  d01x01y02= hists_wc["/ATLAS_2014_I1282447/d01-x01-y02"]
 89  d01x01y02.path = "/ATLAS_2014_I1282447/d01-x01-y02"
 90
 91  d01x01y03= hists_wc["/ATLAS_2014_I1282447/d01-x01-y03"]
 92  d01x01y03.path = "/ATLAS_2014_I1282447/d01-x01-y03"
 93
 94  d03x01y01= hists_wc["/ATLAS_2014_I1282447/d03-x01-y01"]
 95  d03x01y01.path = "/ATLAS_2014_I1282447/d03-x01-y01"
 96
 97  d03x01y02= hists_wc["/ATLAS_2014_I1282447/d03-x01-y02"]
 98  d03x01y02.path = "/ATLAS_2014_I1282447/d03-x01-y02"
 99
100  d04x01y01= hists_wc["/ATLAS_2014_I1282447/d04-x01-y01"]
101  d04x01y01.path = "/ATLAS_2014_I1282447/d04-x01-y01"
102
103  d04x01y02= hists_wc["/ATLAS_2014_I1282447/d04-x01-y02"]
104  d04x01y02.path = "/ATLAS_2014_I1282447/d04-x01-y02"
105
106  d04x01y03= hists_wc["/ATLAS_2014_I1282447/d04-x01-y03"]
107  d04x01y03.path = "/ATLAS_2014_I1282447/d04-x01-y03"
108
109  d04x01y04= hists_wc["/ATLAS_2014_I1282447/d04-x01-y04"]
110  d04x01y04.path = "/ATLAS_2014_I1282447/d04-x01-y04"
111
112  d07x01y01= hists_wc["/ATLAS_2014_I1282447/d07-x01-y01"]
113  d07x01y01.path = "/ATLAS_2014_I1282447/d07-x01-y01"
114
115  yoda.write([ratio_d02y01,ratio_d02y02,ratio_d08y01, ratio_wd ,ratio_wdstar,ratio_wd_plus,ratio_wd_minus ,ratio_wdstar_plus,ratio_wdstar_minus,d01x01y01,d01x01y02,d01x01y03,d03x01y01,d03x01y02,d04x01y01,d04x01y02,d04x01y03,d04x01y04,d07x01y01],"validation.yoda")
116
117*/
118//////////////////////////////////////////////////////////////////////////
119
120#include "Rivet/Analysis.hh"
121#include "Rivet/Projections/UnstableParticles.hh"
122#include "Rivet/Projections/WFinder.hh"
123#include "Rivet/Projections/FastJets.hh"
124#include "Rivet/Projections/ChargedFinalState.hh"
125#include "Rivet/Projections/VetoedFinalState.hh"
126
127namespace Rivet {
128
129
130
131  class ATLAS_2014_I1282447 : public Analysis {
132  public:
133
134    /// Constructor
135    ATLAS_2014_I1282447() : Analysis("ATLAS_2014_I1282447")
136    {
137
138    }
139
140
141    /// @name Analysis methods
142    //@{
143
144    /// Book histograms and initialise projections before the run
145    void init() {
146
147      /// @todo Initialise and register projections here
148      UnstableParticles fs;
149
150      Cut cuts = Cuts::etaIn(-2.5, 2.5) & (Cuts::pT > 20*GeV);
151
152      /// should use sample WITHOUT QED radiation off the electron
153      WFinder wfinder_born_el(fs, cuts, PID::ELECTRON, 25*GeV, 8000*GeV, 15*GeV, 0.1, WFinder::ChargedLeptons::PROMPT, WFinder::ClusterPhotons::ALL, WFinder::AddPhotons::YES);
154      declare(wfinder_born_el, "WFinder_born_el");
155
156      WFinder wfinder_born_mu(fs, cuts, PID::MUON    , 25*GeV, 8000*GeV, 15*GeV, 0.1, WFinder::ChargedLeptons::PROMPT, WFinder::ClusterPhotons::ALL, WFinder::AddPhotons::YES);
157      declare(wfinder_born_mu, "WFinder_born_mu");
158
159      // all hadrons that could be coming from a charm decay --
160      // -- for safety, use region -3.5 - 3.5
161      declare(UnstableParticles(Cuts::abseta <3.5), "hadrons");
162
163      // Input for the jets: no neutrinos, no muons, and no electron which passed the electron cuts
164      // also: NO electron, muon or tau (needed due to ATLAS jet truth reconstruction feature)
165      VetoedFinalState veto;
166
167      veto.addVetoOnThisFinalState(wfinder_born_el);
168      veto.addVetoOnThisFinalState(wfinder_born_mu);
169      veto.addVetoPairId(PID::ELECTRON);
170      veto.addVetoPairId(PID::MUON);
171      veto.addVetoPairId(PID::TAU);
172
173      FastJets jets(veto, FastJets::ANTIKT, 0.4);
174      declare(jets, "jets");
175
176      // Book histograms
177
178      // charge separated integrated cross sections
179      book(_hist_wcjet_charge  ,"d01-x01-y01");
180      book(_hist_wd_charge     ,"d01-x01-y02");
181      book(_hist_wdstar_charge ,"d01-x01-y03");
182
183      // charge integrated total cross sections
184      book(_hist_wcjet_ratio,"d02-x01-y01");
185      book(_hist_wd_ratio   ,"d02-x01-y02");
186
187      book(_hist_wcjet_minus ,"d02-x01-y01_minus");
188      book(_hist_wd_minus    ,"d02-x01-y02_minus");
189
190      book(_hist_wcjet_plus  ,"d02-x01-y01_plus");
191      book(_hist_wd_plus     ,"d02-x01-y02_plus");
192
193      // eta distributions
194      book(_hist_wplus_wcjet_eta_lep   ,"d03-x01-y01");
195      book(_hist_wminus_wcjet_eta_lep  ,"d03-x01-y02");
196
197      book(_hist_wplus_wdminus_eta_lep ,"d04-x01-y01");
198      book(_hist_wminus_wdplus_eta_lep ,"d04-x01-y02");
199      book(_hist_wplus_wdstar_eta_lep  ,"d04-x01-y03");
200      book(_hist_wminus_wdstar_eta_lep ,"d04-x01-y04");
201
202      // ratio of cross section (WD over W inclusive) // postprocess!
203      book(_hist_w_inc             ,"d05-x01-y01");
204      book(_hist_wd_winc_ratio    ,"d05-x01-y02");
205      book(_hist_wdstar_winc_ratio,"d05-x01-y03");
206
207      // ratio of cross section (WD over W inclusive -- function of pT of D meson)
208      book(_hist_wplusd_wplusinc_pt_ratio      ,"d06-x01-y01");
209      book(_hist_wminusd_wminusinc_pt_ratio    ,"d06-x01-y02");
210      book(_hist_wplusdstar_wplusinc_pt_ratio  ,"d06-x01-y03");
211      book(_hist_wminusdstar_wminusinc_pt_ratio,"d06-x01-y04");
212
213      // could use for postprocessing!
214      book(_hist_wplusd_wplusinc_pt       ,"d06-x01-y01_wplus");
215      book(_hist_wminusd_wminusinc_pt     ,"d06-x01-y02_wminus");
216      book(_hist_wplusdstar_wplusinc_pt   ,"d06-x01-y03_wplus");
217      book(_hist_wminusdstar_wminusinc_pt ,"d06-x01-y04_wminus");
218
219      book(_hist_wplus_winc  ,"d06-x01-y01_winc");
220      book(_hist_wminus_winc ,"d06-x01-y02_winc");
221
222      // jet multiplicity of charge integrated W+cjet cross section (+0 or +1 jet in addition to the charm jet)
223      book(_hist_wcjet_jets  ,"d07-x01-y01");
224
225      // jet multiplicity of W+cjet cross section ratio (+0 or +1 jet in addition to the charm jet)
226      book(_hist_wcjet_jets_ratio ,"d08-x01-y01");
227      book(_hist_wcjet_jets_plus   ,"d08-x01-y01_plus");
228      book(_hist_wcjet_jets_minus  ,"d08-x01-y01_minus");
229
230    }
231
232
233    /// Perform the per-event analysis
234    void analyze(const Event& event) {
235
236      double charge_weight = 0; // account for OS/SS events
237
238      int    lepton_charge = 0;
239      double lepton_eta    = 0.;
240
241      /// Find leptons
242      const WFinder& wfinder_born_el = apply<WFinder>(event, "WFinder_born_el");
243      const WFinder& wfinder_born_mu = apply<WFinder>(event, "WFinder_born_mu");
244
245      if (wfinder_born_el.empty() && wfinder_born_mu.empty()) {
246        MSG_DEBUG("No W bosons found");
247        vetoEvent;
248      }
249
250      bool keepevent = false;
251
252      //check electrons
253      if (!wfinder_born_el.empty()) {
254        const FourMomentum nu = wfinder_born_el.constituentNeutrinos()[0];
255        if (wfinder_born_el.mT() > 40*GeV && nu.pT() > 25*GeV) {
256          keepevent = true;
257          lepton_charge = wfinder_born_el.constituentLeptons()[0].charge();
258          lepton_eta = fabs(wfinder_born_el.constituentLeptons()[0].pseudorapidity());
259        }
260      }
261
262      //check muons
263      if (!wfinder_born_mu.empty()) {
264        const FourMomentum nu = wfinder_born_mu.constituentNeutrinos()[0];
265        if (wfinder_born_mu.mT() > 40*GeV && nu.pT() > 25*GeV) {
266          keepevent = true;
267          lepton_charge = wfinder_born_mu.constituentLeptons()[0].charge();
268          lepton_eta = fabs(wfinder_born_mu.constituentLeptons()[0].pseudorapidity());
269        }
270      }
271
272      if (!keepevent) {
273        MSG_DEBUG("Event does not pass mT and MET cuts");
274        vetoEvent;
275      }
276
277      if (lepton_charge > 0) {
278        _hist_wplus_winc->fill(10.);
279        _hist_wplus_winc->fill(16.);
280        _hist_wplus_winc->fill(30.);
281        _hist_wplus_winc->fill(60.);
282        _hist_w_inc->fill(+1);
283      }
284      else if (lepton_charge < 0) {
285        _hist_wminus_winc->fill(10.);
286        _hist_wminus_winc->fill(16.);
287        _hist_wminus_winc->fill(30.);
288        _hist_wminus_winc->fill(60.);
289        _hist_w_inc->fill(-1);
290      }
291
292      // Find hadrons in the event
293      const UnstableParticles& fs = apply<UnstableParticles>(event, "hadrons");
294
295      /// FIND Different channels
296      // 1: wcjet
297      // get jets
298      const Jets& jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT>25.0*GeV && Cuts::abseta<2.5);
299      // loop over jets to select jets used to match to charm
300      Jets js;
301      int    matched_charmHadron = 0;
302      double charm_charge = 0.;
303      int    njets = 0;
304      int    nj = 0;
305      bool   mat_jet = false;
306
307      double ptcharm = 0;
308      if (matched_charmHadron > -1) {
309        for (const Jet& j : jets) {
310          mat_jet = false;
311
312          njets += 1;
313          for (const Particle& p : fs.particles()) {
314            /// @todo Avoid touching HepMC!
315            ConstGenParticlePtr part = p.genParticle();
316            if (p.hasCharm()) {
317              //if(isFromBDecay(p)) continue;
318              if (p.fromBottom()) continue;
319              if (p.pT() < 5*GeV ) continue;
320              if (hasCharmedChildren(part)) continue;
321              if (deltaR(p, j) < 0.3) {
322                mat_jet = true;
323                if (p.pT() > ptcharm) {
324                  charm_charge = part->pdg_id();
325                  ptcharm = p.pT();
326                }
327              }
328            }
329          }
330          if (mat_jet) nj++;
331        }
332
333        if (charm_charge * lepton_charge > 0)  charge_weight = -1;
334        else charge_weight = +1;
335
336        if (nj == 1)  {
337          if (lepton_charge > 0) {
338            _hist_wcjet_charge        ->fill(         1, charge_weight);
339            _hist_wcjet_plus          ->fill(         0, charge_weight);
340            _hist_wplus_wcjet_eta_lep ->fill(lepton_eta, charge_weight);
341            _hist_wcjet_jets_plus     ->fill(njets-1   , charge_weight);
342          }
343          else if (lepton_charge < 0) {
344            _hist_wcjet_charge        ->fill(        -1, charge_weight);
345            _hist_wcjet_minus         ->fill(         0, charge_weight);
346            _hist_wminus_wcjet_eta_lep->fill(lepton_eta, charge_weight);
347            _hist_wcjet_jets_minus    ->fill(njets-1   , charge_weight);
348          }
349
350          _hist_wcjet_jets->fill(njets-1, charge_weight);
351        }
352      }
353
354      // // 1/2: w+d(*) meson
355
356      for (const Particle& p : fs.particles()) {
357
358        /// @todo Avoid touching HepMC!
359        ConstGenParticlePtr part = p.genParticle();
360        if (p.pT() < 8*GeV)       continue;
361        if (fabs(p.eta()) > 2.2)  continue;
362
363        // W+D
364        if (abs(part->pdg_id()) == 411) {
365          if (lepton_charge * part->pdg_id() > 0)  charge_weight = -1;
366          else                                     charge_weight = +1;
367
368          // fill histos
369          if (lepton_charge > 0) {
370            _hist_wd_charge            ->fill(         1, charge_weight);
371            _hist_wd_plus              ->fill(         0, charge_weight);
372            _hist_wplus_wdminus_eta_lep->fill(lepton_eta, charge_weight);
373            _hist_wplusd_wplusinc_pt   ->fill(    p.pT(), charge_weight);
374          }
375          else if (lepton_charge < 0) {
376            _hist_wd_charge            ->fill(        -1, charge_weight);
377            _hist_wd_minus             ->fill(         0, charge_weight);
378            _hist_wminus_wdplus_eta_lep->fill(lepton_eta, charge_weight);
379            _hist_wminusd_wminusinc_pt ->fill(p.pT()    , charge_weight);
380          }
381        }
382
383        // W+Dstar
384        if ( abs(part->pdg_id()) == 413 ) {
385          if (lepton_charge*part->pdg_id() > 0) charge_weight = -1;
386          else charge_weight = +1;
387
388          if (lepton_charge > 0) {
389            _hist_wdstar_charge->fill(+1, charge_weight);
390            _hist_wd_plus->fill( 0, charge_weight);
391            _hist_wplus_wdstar_eta_lep->fill( lepton_eta, charge_weight);
392            _hist_wplusdstar_wplusinc_pt->fill(  p.pT(), charge_weight);
393          }
394          else if (lepton_charge < 0) {
395            _hist_wdstar_charge->fill(-1, charge_weight);
396            _hist_wd_minus->fill(0, charge_weight);
397            _hist_wminus_wdstar_eta_lep->fill(lepton_eta, charge_weight);
398            _hist_wminusdstar_wminusinc_pt->fill(p.pT(), charge_weight);
399          }
400        }
401      }
402
403    }
404
405
406    /// Normalise histograms etc., after the run
407    void finalize() {
408
409      const double sf = crossSection() / sumOfWeights();
410
411      // norm to cross section
412      // d01
413      scale(_hist_wcjet_charge,  sf);
414      scale(_hist_wd_charge,     sf);
415      scale(_hist_wdstar_charge, sf);
416
417      //d02
418      scale(_hist_wcjet_plus,  sf);
419      scale(_hist_wcjet_minus, sf);
420      scale(_hist_wd_plus,     sf);
421      scale(_hist_wd_minus,    sf);
422
423      divide(_hist_wcjet_plus, _hist_wcjet_minus, _hist_wcjet_ratio);
424      divide(_hist_wd_plus,    _hist_wd_minus,    _hist_wd_ratio   );
425
426      //d03
427      scale(_hist_wplus_wcjet_eta_lep,  sf);
428      scale(_hist_wminus_wcjet_eta_lep, sf);
429
430      //d04
431      scale(_hist_wplus_wdminus_eta_lep, crossSection()/sumOfWeights());
432      scale(_hist_wminus_wdplus_eta_lep, crossSection()/sumOfWeights());
433
434      scale(_hist_wplus_wdstar_eta_lep , crossSection()/sumOfWeights());
435      scale(_hist_wminus_wdstar_eta_lep, crossSection()/sumOfWeights());
436
437      //d05
438      scale(_hist_w_inc, 0.01 * sf); // in percent --> /100
439      divide(_hist_wd_charge,     _hist_w_inc, _hist_wd_winc_ratio    );
440      divide(_hist_wdstar_charge, _hist_w_inc, _hist_wdstar_winc_ratio);
441
442      //d06, in percentage!
443      scale(_hist_wplusd_wplusinc_pt,       sf);
444      scale(_hist_wminusd_wminusinc_pt,     sf);
445      scale(_hist_wplusdstar_wplusinc_pt,   sf);
446      scale(_hist_wminusdstar_wminusinc_pt, sf);
447
448      scale(_hist_wplus_winc,  0.01 * sf); // in percent --> /100
449      scale(_hist_wminus_winc, 0.01 * sf); // in percent --> /100
450
451      divide(_hist_wplusd_wplusinc_pt,       _hist_wplus_winc , _hist_wplusd_wplusinc_pt_ratio      );
452      divide(_hist_wminusd_wminusinc_pt,     _hist_wminus_winc, _hist_wminusd_wminusinc_pt_ratio    );
453      divide(_hist_wplusdstar_wplusinc_pt,   _hist_wplus_winc , _hist_wplusdstar_wplusinc_pt_ratio  );
454      divide(_hist_wminusdstar_wminusinc_pt, _hist_wminus_winc, _hist_wminusdstar_wminusinc_pt_ratio);
455
456
457      //d07
458      scale(_hist_wcjet_jets, sf);
459
460      //d08
461      scale(_hist_wcjet_jets_minus, sf);
462      scale(_hist_wcjet_jets_plus,  sf);
463      divide(_hist_wcjet_jets_plus, _hist_wcjet_jets_minus , _hist_wcjet_jets_ratio);
464    }
465
466    //@}
467
468
469  private:
470
471    // Data members like post-cuts event weight counters go here
472
473    // Check whether particle comes from b-decay
474    bool isFromBDecay(const Particle& p) {
475      
476      /// @todo I think we can just replicated the original behaviour with this call
477      /// Note slight difference to Rivet's native Particle::fromBottom method!
478      return p.hasAncestorWith([](const Particle &p)->bool{return p.hasBottom();});
479      /*
480      bool isfromB = false;
481
482      if (p.genParticle() == nullptr)  return false;
483
484      ConstGenParticlePtr part = p.genParticle();
485      ConstGenVertexPtr ivtx = part->production_vertex();
486      while (ivtx) {
487        if (ivtx->particles_in().size() < 1) {
488          isfromB = false;
489          break;
490        }
491        const HepMC::GenVertex::particles_in_const_iterator iPart_invtx = ivtx->particles_in_const_begin();
492        part = (*iPart_invtx);
493        if (!part) {
494          isfromB = false;
495          break;
496        }
497        isfromB = PID::hasBottom(part->pdg_id());
498
499        if (isfromB == true)  break;
500        ivtx = part->production_vertex();
501        if ( part->pdg_id() == 2212 || !ivtx )  break; // reached beam
502      }
503      return isfromB;
504       */
505    }
506
507
508    // Check whether particle has charmed children
509    /// @todo Use built-in method and avoid HepMC!
510    bool hasCharmedChildren(ConstGenParticlePtr part) {
511
512      bool hasCharmedChild = false;
513      if (part == nullptr)  return false;
514
515      ConstGenVertexPtr ivtx = part->end_vertex();
516      if (ivtx == nullptr)  return false;
517
518      // if (ivtx->particles_out_size() < 2) return false;
519      //HepMC::GenVertex::particles_out_const_iterator iPart_invtx = ivtx->particles_out_const_begin();
520      //HepMC::GenVertex::particles_out_const_iterator end_invtx = ivtx->particles_out_const_end();
521
522      for(ConstGenParticlePtr p2: HepMCUtils::particles(ivtx, Relatives::CHILDREN)){
523        if (p2 == part)  continue;
524        hasCharmedChild = PID::hasCharm(p2->pdg_id());
525        if (hasCharmedChild == true)  break;
526        hasCharmedChild = hasCharmedChildren(p2);
527        if (hasCharmedChild == true)  break;
528      }
529      return hasCharmedChild;
530    }
531
532
533  private:
534
535    /// @name Histograms
536    //@{
537
538    //d01-x01-
539    Histo1DPtr   _hist_wcjet_charge;
540    Histo1DPtr   _hist_wd_charge;
541    Histo1DPtr   _hist_wdstar_charge;
542
543    //d02-x01-
544    Scatter2DPtr _hist_wcjet_ratio;
545    Scatter2DPtr _hist_wd_ratio;
546    Histo1DPtr _hist_wcjet_plus;
547    Histo1DPtr _hist_wd_plus;
548
549    Histo1DPtr _hist_wcjet_minus;
550    Histo1DPtr _hist_wd_minus;
551
552    //d03-x01-
553    Histo1DPtr _hist_wplus_wcjet_eta_lep;
554    Histo1DPtr _hist_wminus_wcjet_eta_lep;
555
556    //d04-x01-
557    Histo1DPtr _hist_wplus_wdminus_eta_lep;
558    Histo1DPtr _hist_wminus_wdplus_eta_lep;
559
560    //d05-x01-
561    Histo1DPtr _hist_wplus_wdstar_eta_lep;
562    Histo1DPtr _hist_wminus_wdstar_eta_lep;
563
564
565    // postprocessing histos
566    //d05-x01
567    Histo1DPtr _hist_w_inc;
568    Scatter2DPtr _hist_wd_winc_ratio;
569    Scatter2DPtr _hist_wdstar_winc_ratio;
570
571    //d06-x01
572    Histo1DPtr _hist_wplus_winc;
573    Histo1DPtr _hist_wminus_winc;
574
575    Scatter2DPtr _hist_wplusd_wplusinc_pt_ratio;
576    Scatter2DPtr _hist_wminusd_wminusinc_pt_ratio;
577    Scatter2DPtr _hist_wplusdstar_wplusinc_pt_ratio;
578    Scatter2DPtr _hist_wminusdstar_wminusinc_pt_ratio;
579
580    Histo1DPtr _hist_wplusd_wplusinc_pt ;
581    Histo1DPtr _hist_wminusd_wminusinc_pt;
582    Histo1DPtr _hist_wplusdstar_wplusinc_pt;
583    Histo1DPtr _hist_wminusdstar_wminusinc_pt;
584
585    // d07-x01
586    Histo1DPtr _hist_wcjet_jets ;
587
588    //d08-x01
589    Scatter2DPtr  _hist_wcjet_jets_ratio ;
590    Histo1DPtr    _hist_wcjet_jets_plus ;
591    Histo1DPtr    _hist_wcjet_jets_minus;
592    //@}
593
594  };
595
596
597  // The hook for the plugin system
598  RIVET_DECLARE_PLUGIN(ATLAS_2014_I1282447);
599
600}