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,  JetAlg::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 weight100 = 100.;  //to get results in %
 38
 39      //keeps jets with pt>20 geV and ordered in decreasing pt
 40      Jets jetAr = apply<FastJets>(event, "AntiKT04").jetsByPt(Cuts::pT > 20*GeV);
 41
 42      int flav[2]={1,1};
 43      vector<FourMomentum> leadjets;
 44
 45      //get b/c-hadrons
 46      vector<ConstGenParticlePtr> B_hadrons, C_hadrons;
 47      vector<ConstGenParticlePtr> allParticles = HepMCUtils::particles(event.genEvent());
 48      for (size_t i = 0; i < allParticles.size(); i++) {
 49        ConstGenParticlePtr p = allParticles.at(i);
 50        if(p->momentum().perp()*GeV < 5) continue;
 51        if ( (Rivet::PID::isHadron ( p->pdg_id() ) &&
 52              Rivet::PID::hasBottom( p->pdg_id() )    ) ) {
 53          B_hadrons.push_back(p);
 54        }
 55        if ( (Rivet::PID::isHadron( p->pdg_id() ) &&
 56              Rivet::PID::hasCharm( p->pdg_id() )    ) ) {
 57          C_hadrons.push_back(p);
 58        }
 59      }
 60
 61      //select dijet
 62      for (const Jet& jet : jetAr) {
 63
 64        const double pT   = jet.pT();
 65        const double absy = jet.absrap();
 66
 67        bool isBjet = jet.bTagged();
 68        bool isCjet = jet.cTagged();
 69
 70        int jetflav=1;
 71        if      (isBjet)jetflav=5;
 72        else if (isCjet)jetflav=4;
 73
 74        if (absy <= 2.1 && leadjets.size() < 2) {
 75
 76          if (pT > 500*GeV) continue;
 77          if ((leadjets.empty() && pT < 40*GeV) || pT < 20*GeV) continue;
 78
 79          leadjets.push_back(jet.momentum());
 80
 81          if (leadjets.size()==1) flav[0] = jetflav;
 82          if (leadjets.size()==2) flav[1] = jetflav;
 83        }
 84      }
 85
 86      if (leadjets.size() < 2) vetoEvent;
 87
 88
 89      double pBinsLJ[7] = {40.,60.,80.,120.,160.,250.,500.};
 90      int    iPBinLJ = -1;
 91
 92      for (int k = 0 ; k < 7 ; k++) {
 93        if (leadjets[0].pT() > pBinsLJ[k]*GeV) iPBinLJ=k;
 94        else break;
 95      }
 96
 97      bool c_ljpt  = (iPBinLJ != -1);
 98      bool c_nljpt = leadjets[1].pT() > 20*GeV;
 99      bool c_dphi  = fabs( deltaPhi(leadjets[0],leadjets[1]) ) > 2.1;
100      bool isDijet = c_ljpt & c_nljpt & c_dphi;
101      if (!isDijet) vetoEvent;
102
103      _h_temp[6]->fill(leadjets[0].pT());
104
105      if (flav[0]==5 && flav[1]==5)                                  // BB dijet
106        _h_temp[0]->fill(leadjets[0].pT(), weight100);
107
108      if ((flav[0]==5 && flav[1]==4) || (flav[0]==4 && flav[1]==5))  // BC dijet
109        _h_temp[1]->fill(leadjets[0].pT(), weight100);
110
111      if (flav[0]==4 && flav[1]==4)                                  // CC dijet
112        _h_temp[2]->fill(leadjets[0].pT(), weight100);
113
114      if ((flav[0]==5 && flav[1]==1) || (flav[0]==1 && flav[1]==5))  // B-light dijet
115        _h_temp[3]->fill(leadjets[0].pT(), weight100);
116
117      if ((flav[0]==4 && flav[1]==1) || (flav[0]==1 && flav[1]==4))  // C-light dijet
118        _h_temp[4]->fill(leadjets[0].pT(), weight100);
119
120      if (flav[0]==1 && flav[1]==1)                                  // light-light dijet
121        _h_temp[5]->fill(leadjets[0].pT(), weight100);
122    }
123
124
125    void finalize() {
126      divide(_h_temp[0], _h_temp[6], _h_results[0]);
127      divide(_h_temp[1], _h_temp[6], _h_results[1]);
128      divide(_h_temp[2], _h_temp[6], _h_results[2]);
129      divide(_h_temp[3], _h_temp[6], _h_results[3]);
130      divide(_h_temp[4], _h_temp[6], _h_results[4]);
131      divide(_h_temp[5], _h_temp[6], _h_results[5]);
132    }
133
134  private:
135
136    Histo1DPtr   _h_temp[7];
137    Estimate1DPtr _h_results[6];
138  };
139
140
141  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1188891);
142
143
144}