1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Particle.hh"
  6
  7namespace Rivet {
  8
  9
 10  class ATLAS_2012_I1188891 : public Analysis {
 11  public:
 12
 13
 14    ATLAS_2012_I1188891()
 15      : Analysis("ATLAS_2012_I1188891")
 16    {    }
 17
 18
 19    void init() {
 20
 21      const FinalState fs;
 22      FastJets fj04(fs,  FastJets::ANTIKT, 0.4);
 23      declare(fj04, "AntiKT04");
 24
 25      string histotitle[7]={"BBfraction","BCfraction","CCfraction","BUfraction","CUfraction","UUfraction","Total"};
 26      for (int i = 0 ; i < 7 ; i++){
 27        book(_h_temp[i] ,"TMP/"+histotitle[i],refData(1,1,1));
 28        if (i < 6) {
 29          book(_h_results[i], i+1, 1, 1);
 30        }
 31      }
 32    }
 33
 34
 35    void analyze(const Event& event) {
 36
 37      double weight    = 1.0;
 38      double weight100 = 1.0 * 100.;  //to get results in %
 39
 40      //keeps jets with pt>20 geV and ordered in decreasing pt
 41      Jets jetAr = apply<FastJets>(event, "AntiKT04").jetsByPt(20*GeV);
 42
 43      int flav[2]={1,1};
 44      vector<FourMomentum> leadjets;
 45
 46      //get b/c-hadrons
 47      vector<ConstGenParticlePtr> B_hadrons, C_hadrons;
 48      vector<ConstGenParticlePtr> allParticles = HepMCUtils::particles(event.genEvent());
 49      for (size_t i = 0; i < allParticles.size(); i++) {
 50        ConstGenParticlePtr p = allParticles.at(i);
 51        if(p->momentum().perp()*GeV < 5) continue;
 52        if ( (Rivet::PID::isHadron ( p->pdg_id() ) &&
 53              Rivet::PID::hasBottom( p->pdg_id() )    ) ) {
 54          B_hadrons.push_back(p);
 55        }
 56        if ( (Rivet::PID::isHadron( p->pdg_id() ) &&
 57              Rivet::PID::hasCharm( p->pdg_id() )    ) ) {
 58          C_hadrons.push_back(p);
 59        }
 60      }
 61
 62      //select dijet
 63      for (const Jet& jet : jetAr) {
 64
 65        const double pT   = jet.pT();
 66        const double absy = jet.absrap();
 67
 68        bool isBjet = jet.bTagged();
 69        bool isCjet = jet.cTagged();
 70
 71        int jetflav=1;
 72        if      (isBjet)jetflav=5;
 73        else if (isCjet)jetflav=4;
 74
 75        if (absy <= 2.1 && leadjets.size() < 2) {
 76
 77          if (pT > 500*GeV) continue;
 78          if ((leadjets.empty() && pT < 40*GeV) || pT < 20*GeV) continue;
 79
 80          leadjets.push_back(jet.momentum());
 81
 82          if (leadjets.size()==1) flav[0] = jetflav;
 83          if (leadjets.size()==2) flav[1] = jetflav;
 84        }
 85      }
 86
 87      if (leadjets.size() < 2) vetoEvent;
 88
 89
 90      double pBinsLJ[7] = {40.,60.,80.,120.,160.,250.,500.};
 91      int    iPBinLJ = -1;
 92
 93      for (int k = 0 ; k < 7 ; k++) {
 94        if (leadjets[0].pT() > pBinsLJ[k]*GeV) iPBinLJ=k;
 95        else break;
 96      }
 97
 98      bool c_ljpt  = (iPBinLJ != -1);
 99      bool c_nljpt = leadjets[1].pT() > 20*GeV;
100      bool c_dphi  = fabs( deltaPhi(leadjets[0],leadjets[1]) ) > 2.1;
101      bool isDijet = c_ljpt & c_nljpt & c_dphi;
102      if (!isDijet) vetoEvent;
103
104      _h_temp[6]->fill(leadjets[0].pT(), weight);
105
106      if (flav[0]==5 && flav[1]==5)                                  // BB dijet
107        _h_temp[0]->fill(leadjets[0].pT(), weight100);
108
109      if ((flav[0]==5 && flav[1]==4) || (flav[0]==4 && flav[1]==5))  // BC dijet
110        _h_temp[1]->fill(leadjets[0].pT(), weight100);
111
112      if (flav[0]==4 && flav[1]==4)                                  // CC dijet
113        _h_temp[2]->fill(leadjets[0].pT(), weight100);
114
115      if ((flav[0]==5 && flav[1]==1) || (flav[0]==1 && flav[1]==5))  // B-light dijet
116        _h_temp[3]->fill(leadjets[0].pT(), weight100);
117
118      if ((flav[0]==4 && flav[1]==1) || (flav[0]==1 && flav[1]==4))  // C-light dijet
119        _h_temp[4]->fill(leadjets[0].pT(), weight100);
120
121      if (flav[0]==1 && flav[1]==1)                                  // light-light dijet
122        _h_temp[5]->fill(leadjets[0].pT(), weight100);
123    }
124
125
126    void finalize() {
127      divide(_h_temp[0], _h_temp[6], _h_results[0]);
128      divide(_h_temp[1], _h_temp[6], _h_results[1]);
129      divide(_h_temp[2], _h_temp[6], _h_results[2]);
130      divide(_h_temp[3], _h_temp[6], _h_results[3]);
131      divide(_h_temp[4], _h_temp[6], _h_results[4]);
132      divide(_h_temp[5], _h_temp[6], _h_results[5]);
133    }
134
135  private:
136
137    Histo1DPtr   _h_temp[7];
138    Scatter2DPtr _h_results[6];
139  };
140
141
142  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1188891);
143
144
145}