DESY Hbb Analysis Framework
JetAnalyser.cc
Go to the documentation of this file.
2 
3 // system include files
4 #include "boost/program_options.hpp"
5 #include "boost/algorithm/string.hpp"
6 #include <string>
7 #include <iostream>
8 #include <fstream>
9 #include <vector>
10 //
11 // user include files
13 
14 //
15 // class declaration
16 //
17 
18 using namespace analysis;
19 using namespace analysis::tools;
20 
22 {
23 }
24 
25 
26 JetAnalyser::JetAnalyser(int argc, char * argv[]) : BaseAnalyser(argc,argv)
27 {
28  // Physics objects
29  // Jets
30  jetsanalysis_ = ( analysis_->addTree<Jet> ("Jets",config_->jetsCollection()) != nullptr );
31 
32  applyjer_ = false;
33 
34  if ( config_->btagScaleFactors() != "" )
35  {
36  bsf_reader_["loose"] = analysis_->btagCalibration(config_->btagAlgorithm(), config_->btagScaleFactors(), "loose");
37  bsf_reader_["medium"] = analysis_->btagCalibration(config_->btagAlgorithm(), config_->btagScaleFactors(), "medium");
38  bsf_reader_["tight"] = analysis_->btagCalibration(config_->btagAlgorithm(), config_->btagScaleFactors(), "tight");
39  }
40 
41  if ( config_->jerPtRes() != "" && config_->jerSF() != "" && this->genJetsAnalysis() ) // FIXME: check if files exist
42  {
43  jerinfo_ = analysis_->jetResolutionInfo(config_->jerPtRes(),config_->jerSF());
44  applyjer_ = ( jerinfo_ != nullptr && jetsanalysis_ );
45  }
46 
47  if ( config_->isMC() )
48  {
49  flavours_ = {"udsg","c","b"};
50  if ( config_->useJetsExtendedFlavour() )
51  {
52  flavours_.push_back("cc");
53  flavours_.push_back("bb");
54  }
55  }
56 // histograms("jet",config_->nJetsMin());
57 
58 }
59 
61 {
62  // do anything here that needs to be done at desctruction time
63  // (e.g. close files, deallocate resources etc.)
64 }
65 
66 
67 //
68 // member functions
69 //
70 // ------------ method called for each event ------------
71 
73 {
74  jets_.clear();
75  selectedJets_.clear();
76 
77  // trigger emulation
78  // L1 jets
79  std::string triggerObjectsL1Jets;
80  if ( config_->triggerObjectsL1Jets() != "" )
81  {
82  triggerObjectsL1Jets = config_->triggerObjectsL1Jets();
83  if ( config_->triggerEmulateL1Jets() != "" && config_->triggerEmulateL1JetsNMin() > 0 )
84  {
85  int nmin = config_->triggerEmulateL1JetsNMin();
86  float ptmin = config_->triggerEmulateL1JetsPtMin();
87  float etamax = config_->triggerEmulateL1JetsEtaMax();
88  std::string newL1Jets = config_->triggerEmulateL1Jets();
89  triggerEmulation(triggerObjectsL1Jets,nmin,ptmin,etamax,newL1Jets);
90  triggerObjectsL1Jets = newL1Jets;
91  }
92  }
93 
94  // Calo jets
95  std::string triggerObjectsCaloJets;
96  if ( config_->triggerObjectsCaloJets() != "" )
97  {
98  triggerObjectsCaloJets = config_->triggerObjectsCaloJets();
99  if ( config_->triggerEmulateCaloJets() != "" && config_->triggerEmulateCaloJetsNMin() > 0 )
100  {
101  int nmin = config_->triggerEmulateCaloJetsNMin();
102  float ptmin = config_->triggerEmulateCaloJetsPtMin();
103  float etamax = config_->triggerEmulateCaloJetsEtaMax();
104  std::string newCaloJets = config_->triggerEmulateCaloJets();
105  triggerEmulation(triggerObjectsCaloJets,nmin,ptmin,etamax,newCaloJets);
106  triggerObjectsCaloJets = newCaloJets;
107  }
108  }
109 
110  // PF jets
111  std::string triggerObjectsPFJets;
112  if ( config_->triggerObjectsPFJets() != "" )
113  {
114  triggerObjectsPFJets = config_->triggerObjectsPFJets();
115  if ( config_->triggerEmulatePFJets() != "" && config_->triggerEmulatePFJetsNMin() > 0 )
116  {
117  int nmin = config_->triggerEmulatePFJetsNMin();
118  float ptmin = config_->triggerEmulatePFJetsPtMin();
119  float etamax = config_->triggerEmulatePFJetsEtaMax();
120  std::string newPFJets = config_->triggerEmulatePFJets();
121  triggerEmulation(triggerObjectsPFJets,nmin,ptmin,etamax,newPFJets);
122  triggerObjectsPFJets = newPFJets;
123  }
124  }
125 
126 
127  if ( ! jetsanalysis_ ) return false;
128 
129  cutflow(Form("Using Jet collection: %s",(config_->jetsCollection()).c_str()));
130 
131 
132  if ( config_->triggerObjectsL1Jets() != "" )
133  analysis_->match<Jet,TriggerObject>("Jets",triggerObjectsL1Jets,config_-> triggerMatchL1JetsDrMax());
134  if ( config_->triggerObjectsCaloJets() != "" )
135  analysis_->match<Jet,TriggerObject>("Jets",triggerObjectsCaloJets,config_-> triggerMatchCaloJetsDrMax());
136  if ( config_->triggerObjectsPFJets() != "" )
137  analysis_->match<Jet,TriggerObject>("Jets",triggerObjectsPFJets,config_-> triggerMatchPFJetsDrMax());
138  analysis_->match<Jet,TriggerObject>("Jets",config_->triggerObjectsBJets(),config_-> triggerMatchCaloBJetsDrMax());
139 
140  // std::shared_ptr< Collection<Jet> >
141  auto jets = analysis_->collection<Jet>("Jets");
142 
143  if ( this->genParticlesAnalysis() && config_ -> useJetsExtendedFlavour() )
144  {
145  auto particles = analysis_->collection<GenParticle>("GenParticles");
146  jets->associatePartons(particles,0.4,1.,config_->pythia8());
147  }
148 
149  if ( this->genJetsAnalysis() )
150  {
151  auto genjets = analysis_->collection<GenJet>("GenJets");
152  jets->addGenJets(genjets);
153  }
154 
155  for ( int j = 0 ; j < jets->size() ; ++j ) jets_.push_back(std::make_shared<Jet>(jets->at(j)));
156 
158 
159  return true;
160 }
161 
162 
163 void JetAnalyser::jets(const std::string & col)
164 {
165  analysis_->addTree<Jet> ("Jets",col);
166 }
167 
168 void JetAnalyser::jetHistograms( const std::string & label )
169 {
170  this -> jetHistograms(config_->nJetsMin(), label);
171 }
172 void JetAnalyser::jetHistograms( const int & n, const std::string & label )
173 {
174  this->output()->cd();
175  this->output()->mkdir(label.c_str());
176  this->output()->cd(label.c_str());
177 
178  n_hjets_ = n;
179 
180  h1_[Form("jet_hist_weight_%s",label.c_str())] = std::make_shared<TH1F>(Form("jet_hist_weight_%s",label.c_str()) , Form("jet_hist_weight_%s",label.c_str()) ,1 , 0. , 1. );
181 
182  // btag variable binning
183  float min1 = 0.0;
184  float max1 = 0.1;
185  float size1 = 0.0001;
186  int nbins1 = int((max1-min1)/size1);
187  float min2 = 0.1;
188  float max2 = 0.9;
189  float size2 = 0.001;
190  int nbins2 = int((max2-min2)/size2);
191  float min3 = 0.9;
192  float max3 = 1.0;
193  float size3 = 0.0001;
194  int nbins3 = int((max3-min3)/size3);
195  int nbins_btag = nbins1+nbins2+nbins3;
196 
197  std::vector<float> bins_btag;
198  bins_btag.clear();
199  int counter = 0;
200  for ( int i = 0; i<nbins1; ++i) { bins_btag.push_back(min1 + size1*i); ++counter; }
201  for ( int i = 0; i<nbins2; ++i) { bins_btag.push_back(min2 + size2*i); ++counter; }
202  for ( int i = 0; i<nbins3+1; ++i) { bins_btag.push_back(min3 + size3*i); ++counter; }
203 
204  // uniform binning for btag (comment it out if you want the variable binning above)
205  float size = 0.0002;
206  nbins_btag = int(1./size);
207  bins_btag.clear();
208  for ( int i = 0; i<nbins_btag+1; ++i) { bins_btag.push_back(size*i); ++counter; }
209 
210 
211  for ( int j = 0; j < n; ++j ) // loop over jets
212  {
213  // 1D histograms
214  h1_[Form("pt_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d" , j+1) , Form("pt_jet%d_%s" , j+1,label.c_str()) ,1500 , 0 , 1500 );
215  h1_[Form("eta_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("eta_jet%d" , j+1) , Form("eta_jet%d_%s" , j+1,label.c_str()) , 600 , -3, 3 );
216  h1_[Form("phi_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("phi_jet%d" , j+1) , Form("phi_jet%d_%s" , j+1,label.c_str()) , 360 , -180, 180 );
217  h1_[Form("btag_jet%d_%s", j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_jet%d", j+1) , Form("btag_jet%d_%s", j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
218  h1_[Form("btaglog_jet%d_%s", j+1,label.c_str())] = std::make_shared<TH1F>(Form("btaglog_jet%d", j+1) , Form("btaglog_jet%d_%s", j+1,label.c_str()) , 200 , 1.e-6, 10 );
219  h1_[Form("qglikelihood_jet%d_%s", j+1,label.c_str())] = std::make_shared<TH1F>(Form("qglikelihood_jet%d", j+1) , Form("qglikelihood_jet%d_%s", j+1,label.c_str()) , 200 , 0, 1 );
220  h1_[Form("nconstituents_jet%d_%s", j+1,label.c_str())] = std::make_shared<TH1F>(Form("nconstituents_jet%d", j+1) , Form("nconstituents_jet%d_%s", j+1,label.c_str()) , 200 , 0, 200 );
221  // pt distributions separate for barrel and endcap, minus and plus sides and overlaps
222  if ( config_-> histogramJetsRegionSplit() )
223  {
224  h1_[Form("pt_jet%d_me_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d_me" , j+1) , Form("pt_jet%d_me_%s" , j+1,label.c_str()) ,1500 , 0 , 1500 );
225  h1_[Form("pt_jet%d_pe_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d_pe" , j+1) , Form("pt_jet%d_pe_%s" , j+1,label.c_str()) ,1500 , 0 , 1500 );
226  h1_[Form("pt_jet%d_mb_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d_mb" , j+1) , Form("pt_jet%d_mb_%s" , j+1,label.c_str()) ,1500 , 0 , 1500 );
227  h1_[Form("pt_jet%d_pb_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d_pb" , j+1) , Form("pt_jet%d_pb_%s" , j+1,label.c_str()) ,1500 , 0 , 1500 );
228  h1_[Form("pt_jet%d_mbe_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d_mbe" , j+1) , Form("pt_jet%d_mbe_%s" , j+1,label.c_str()) ,1500 , 0 , 1500 );
229  h1_[Form("pt_jet%d_pbe_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d_pbe" , j+1) , Form("pt_jet%d_pbe_%s" , j+1,label.c_str()) ,1500 , 0 , 1500 );
230  }
231 
232  h1_[Form("pt_jet%d_%s" , j+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d p_{T} [GeV]",j+1));
233  h1_[Form("eta_jet%d_%s" , j+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d #eta",j+1));
234  h1_[Form("phi_jet%d_%s" , j+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d #phi",j+1));
235  h1_[Form("btag_jet%d_%s", j+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d btag discriminator",j+1));
236  h1_[Form("btaglog_jet%d_%s", j+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d -ln(1-btag discriminator)",j+1));
237  h1_[Form("qglikelihood_jet%d_%s", j+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d q-g likelihood",j+1));
238  h1_[Form("nconstituents_jet%d_%s", j+1,label.c_str())]-> GetXaxis() -> SetTitle(Form("Jet %d n constituents",j+1));
239 
240  if ( config_->btagAlgorithm() == "deepcsv")
241  {
242  h1_[Form("btag_light_jet%d_%s", j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_light_jet%d", j+1) , Form("btag_light_jet%d_%s", j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
243  h1_[Form("btag_c_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_c_jet%d" , j+1) , Form("btag_c_jet%d_%s" , j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
244  h1_[Form("btag_b_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_b_jet%d" , j+1) , Form("btag_b_jet%d_%s" , j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
245  h1_[Form("btag_bb_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_bb_jet%d" , j+1) , Form("btag_bb_jet%d_%s" , j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
246  h1_[Form("btag_cc_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_cc_jet%d" , j+1) , Form("btag_cc_jet%d_%s" , j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
247  }
248  if ( config_->btagAlgorithm() == "deepflavour" || config_->btagAlgorithm() == "deepjet" )
249  {
250  h1_[Form("btag_light_jet%d_%s", j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_light_jet%d", j+1) , Form("btag_light_jet%d_%s", j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
251  h1_[Form("btag_g_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_g_jet%d" , j+1) , Form("btag_g_jet%d_%s" , j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
252  h1_[Form("btag_c_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_c_jet%d" , j+1) , Form("btag_c_jet%d_%s" , j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
253  h1_[Form("btag_b_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_b_jet%d" , j+1) , Form("btag_b_jet%d_%s" , j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
254  h1_[Form("btag_bb_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_bb_jet%d" , j+1) , Form("btag_bb_jet%d_%s" , j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
255  h1_[Form("btag_lepb_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH1F>(Form("btag_lepb_jet%d" , j+1) , Form("btag_lepb_jet%d_%s" , j+1,label.c_str()) , nbins_btag, &bins_btag[0] );
256  }
257 
258  // 2D histograms
259  h2_[Form("pt_eta_jet%d_%s" , j+1,label.c_str())] = std::make_shared<TH2F>(Form("pt_eta_jet%d" , j+1) , Form("pt_eta_jet%d_%s" , j+1,label.c_str()) ,1500 , 0 , 1500, 600, -3, 3 );
260  h2_[Form("pt_eta_jet%d_%s" , j+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d p_{T} [GeV]",j+1));
261  h2_[Form("pt_eta_jet%d_%s" , j+1,label.c_str())] -> GetYaxis() -> SetTitle(Form("Jet %d #eta",j+1));
262 
263  if ( config_->isMC() && config_->histogramJetsPerFlavour() )
264  {
265  for ( auto & flv : flavours_ ) // flavour dependent histograms
266  {
267  // 1D histograms
268  h1_[Form("pt_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d_%s" , j+1, flv.c_str()) , Form("pt_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) ,1500 , 0 , 1500 );
269  h1_[Form("eta_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("eta_jet%d_%s" , j+1, flv.c_str()) , Form("eta_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , 600 , -3, 3 );
270  h1_[Form("phi_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("phi_jet%d_%s" , j+1, flv.c_str()) , Form("phi_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , 360 , -180, 180 );
271  h1_[Form("btag_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_jet%d_%s" , j+1, flv.c_str()) , Form("btag_jet%d_%s_%s", j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
272  h1_[Form("qglikelihood_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("qglikelihood_jet%d_%s" , j+1, flv.c_str()) , Form("qglikelihood_jet%d_%s_%s", j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
273  h1_[Form("nconstituents_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("nconstituents_jet%d_%s" , j+1, flv.c_str()) , Form("nconstituents_jet%d_%s_%s", j+1,label.c_str(),flv.c_str()) , 200 , 0, 200 );
274 
275  h1_[Form("pt_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) p_{T} [GeV]" , j+1, flv.c_str()));
276  h1_[Form("eta_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) #eta" , j+1, flv.c_str()));
277  h1_[Form("phi_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) #phi" , j+1, flv.c_str()));
278  h1_[Form("btag_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) btag discriminator" , j+1, flv.c_str()));
279  h1_[Form("qglikelihood_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) q-g likelihood" , j+1, flv.c_str()));
280  h1_[Form("nconstituents_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) n constituents" , j+1, flv.c_str()));
281 
282  if ( config_->btagAlgorithm() == "deepcsv")
283  {
284  h1_[Form("btag_light_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_light_jet%d_%s", j+1,flv.c_str()) , Form("btag_light_jet%d_%s_%s", j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
285  h1_[Form("btag_c_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_c_jet%d_%s" , j+1,flv.c_str()) , Form("btag_c_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
286  h1_[Form("btag_b_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_b_jet%d_%s" , j+1,flv.c_str()) , Form("btag_b_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
287  h1_[Form("btag_bb_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_bb_jet%d_%s" , j+1,flv.c_str()) , Form("btag_bb_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
288  h1_[Form("btag_cc_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_cc_jet%d_%s" , j+1,flv.c_str()) , Form("btag_cc_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
289  }
290  if ( config_->btagAlgorithm() == "deepflavour" || config_->btagAlgorithm() == "deepjet" )
291  {
292  h1_[Form("btag_light_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_light_jet%d_%s", j+1,flv.c_str()) , Form("btag_light_jet%d_%s_%s", j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
293  h1_[Form("btag_g_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_g_jet%d_%s" , j+1,flv.c_str()) , Form("btag_g_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
294  h1_[Form("btag_c_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_c_jet%d_%s" , j+1,flv.c_str()) , Form("btag_c_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
295  h1_[Form("btag_b_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_b_jet%d_%s" , j+1,flv.c_str()) , Form("btag_b_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
296  h1_[Form("btag_bb_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_bb_jet%d_%s" , j+1,flv.c_str()) , Form("btag_bb_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
297  h1_[Form("btag_lepb_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH1F>(Form("btag_lepb_jet%d_%s" , j+1,flv.c_str()) , Form("btag_lepb_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) , nbins_btag, &bins_btag[0] );
298  }
299  // 2D histograms
300  h2_[Form("pt_eta_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] = std::make_shared<TH2F>(Form("pt_eta_jet%d_%s" , j+1,flv.c_str()) , Form("pt_eta_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str()) ,1500 , 0 , 1500, 600, -3, 3 );
301  h2_[Form("pt_eta_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) p_{T} [GeV]" , j+1,flv.c_str()));
302  h2_[Form("pt_eta_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> GetYaxis() -> SetTitle(Form("Jet %d (%s) #eta" ,j+1,flv.c_str()));
303 
304  }
305  }
306  if ( config_->doDijet() ) // dijet histograms
307  {
308  for ( int k = j+1; k < n && j < n; ++k )
309  {
310  h1_[Form("dptrel_jet%d%d_%s" , j+1,k+1,label.c_str())] = std::make_shared<TH1F>(Form("dptrel_jet%d%d" , j+1,k+1) , Form("dptrel_jet%d%d_%s" , j+1,k+1,label.c_str()) ,1000 , 0,1 );
311  h1_[Form("dpt_jet%d%d_%s" , j+1,k+1,label.c_str())] = std::make_shared<TH1F>(Form("dpt_jet%d%d" , j+1,k+1) , Form("dpt_jet%d%d_%s" , j+1,k+1,label.c_str()) ,1000 , 0,1000 );
312  h1_[Form("dr_jet%d%d_%s" , j+1,k+1,label.c_str())] = std::make_shared<TH1F>(Form("dr_jet%d%d" , j+1,k+1) , Form("dr_jet%d%d_%s" , j+1,k+1,label.c_str()) , 100 , 0, 5 );
313  h1_[Form("deta_jet%d%d_%s" , j+1,k+1,label.c_str())] = std::make_shared<TH1F>(Form("deta_jet%d%d", j+1,k+1) , Form("deta_jet%d%d_%s" , j+1,k+1,label.c_str()) , 100 , 0,10 );
314  h1_[Form("dphi_jet%d%d_%s" , j+1,k+1,label.c_str())] = std::make_shared<TH1F>(Form("dphi_jet%d%d", j+1,k+1) , Form("dphi_jet%d%d_%s" , j+1,k+1,label.c_str()) , 315 , 0, 3.15 );
315  h1_[Form("pt_jet%d%d_%s" , j+1,k+1,label.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d%d" , j+1,k+1) , Form("pt_jet%d%d_%s" , j+1,k+1,label.c_str()) , 300 , 0,3000 );
316  h1_[Form("eta_jet%d%d_%s" , j+1,k+1,label.c_str())] = std::make_shared<TH1F>(Form("eta_jet%d%d" , j+1,k+1) , Form("eta_jet%d%d_%s" , j+1,k+1,label.c_str()) , 200 , -10,10 );
317  h1_[Form("phi_jet%d%d_%s" , j+1,k+1,label.c_str())] = std::make_shared<TH1F>(Form("phi_jet%d%d" , j+1,k+1) , Form("phi_jet%d%d_%s" , j+1,k+1,label.c_str()) , 360 , -180,180 );
318  h1_[Form("m_jet%d%d_%s" , j+1,k+1,label.c_str())] = std::make_shared<TH1F>(Form("m_jet%d%d" , j+1,k+1) , Form("m_jet%d%d_%s" , j+1,k+1,label.c_str()) ,3000 , 0,3000 );
319 
320  h1_[Form("dptrel_jet%d%d_%s" , j+1,k+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("#DeltaP_{T}(Jet %d, Jet %d)/Jet %d p_{T}",j+1,k+1,j+1));
321  h1_[Form("dpt_jet%d%d_%s" , j+1,k+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("#Delta p_{T}(Jet %d, Jet %d) [GeV]",j+1,k+1));
322  h1_[Form("dr_jet%d%d_%s" , j+1,k+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("#DeltaR(Jet %d, Jet %d)",j+1,k+1));
323  h1_[Form("deta_jet%d%d_%s" , j+1,k+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("#Delta#eta(Jet %d, Jet %d)",j+1,k+1));
324  h1_[Form("dphi_jet%d%d_%s" , j+1,k+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("#Delta#phi(Jet %d, Jet %d)",j+1,k+1));
325  h1_[Form("pt_jet%d%d_%s" , j+1,k+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d + Jet %d p_{T} [GeV]",j+1,k+1));
326  h1_[Form("eta_jet%d%d_%s" , j+1,k+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d + Jet %d #eta",j+1,k+1));
327  h1_[Form("phi_jet%d%d_%s" , j+1,k+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d + Jet %d #phi",j+1,k+1));
328  h1_[Form("m_jet%d%d_%s" , j+1,k+1,label.c_str())] -> GetXaxis() -> SetTitle(Form("M_{%d%d} [GeV]",j+1,k+1));
329 
330  if ( config_->isMC() && config_->histogramJetsPerFlavour() )
331  {
332  for ( auto & flv1 : flavours_ )
333  {
334  for ( auto & flv2 : flavours_ )
335  {
336  h1_[Form("dptrel_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] = std::make_shared<TH1F>(Form("dptrel_jet%d%d_%s_%s" , j+1,k+1,flv1.c_str(),flv2.c_str()) , Form("dptrel_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str()) , 100 , 0,1 );
337  h1_[Form("dpt_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] = std::make_shared<TH1F>(Form("dpt_jet%d%d_%s_%s" , j+1,k+1,flv1.c_str(),flv2.c_str()) , Form("dpt_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str()) , 100 , 0,1000 );
338  h1_[Form("dr_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] = std::make_shared<TH1F>(Form("dr_jet%d%d_%s_%s" , j+1,k+1,flv1.c_str(),flv2.c_str()) , Form("dr_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str()) , 100 , 0, 5 );
339  h1_[Form("deta_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] = std::make_shared<TH1F>(Form("deta_jet%d%d_%s_%s" , j+1,k+1,flv1.c_str(),flv2.c_str()) , Form("deta_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str()) , 100 , 0,10 );
340  h1_[Form("dphi_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] = std::make_shared<TH1F>(Form("dphi_jet%d%d_%s_%s" , j+1,k+1,flv1.c_str(),flv2.c_str()) , Form("dphi_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str()) , 315 , 0, 3.15 );
341  h1_[Form("pt_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] = std::make_shared<TH1F>(Form("pt_jet%d%d_%s_%s" , j+1,k+1,flv1.c_str(),flv2.c_str()) , Form("pt_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str()) , 300 , 0,300 );
342  h1_[Form("eta_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] = std::make_shared<TH1F>(Form("eta_jet%d%d_%s_%s" , j+1,k+1,flv1.c_str(),flv2.c_str()) , Form("eta_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str()) , 200 , -10,10 );
343  h1_[Form("phi_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] = std::make_shared<TH1F>(Form("phi_jet%d%d_%s_%s" , j+1,k+1,flv1.c_str(),flv2.c_str()) , Form("phi_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str()) , 360 , -180,180 );
344  h1_[Form("m_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] = std::make_shared<TH1F>(Form("m_jet%d%d_%s_%s" , j+1,k+1,flv1.c_str(),flv2.c_str()) , Form("m_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str()) ,3000 , 0,300 );
345 
346  h1_[Form("dptrel_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> GetXaxis() -> SetTitle(Form("#DeltaP_{T}(Jet %d (%s), Jet %d (%s))/Jet %d p_{T}",j+1,flv1.c_str(),k+1,flv2.c_str(),j+1));
347  h1_[Form("dpt_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> GetXaxis() -> SetTitle(Form("#Delta p_{T}(Jet %d (%s), Jet %d (%s)) [GeV]",j+1,flv1.c_str(),k+1,flv2.c_str()));
348  h1_[Form("dr_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> GetXaxis() -> SetTitle(Form("#DeltaR(Jet %d (%s), Jet %d (%s))",j+1,flv1.c_str(),k+1,flv2.c_str()));
349  h1_[Form("deta_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> GetXaxis() -> SetTitle(Form("#Delta#eta(Jet %d (%s), Jet %d (%s))",j+1,flv1.c_str(),k+1,flv2.c_str()));
350  h1_[Form("dphi_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> GetXaxis() -> SetTitle(Form("#Delta#phi(Jet %d (%s), Jet %d (%s))",j+1,flv1.c_str(),k+1,flv2.c_str()));
351  h1_[Form("pt_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) + Jet %d (%s) p_{T} [GeV]",j+1,flv1.c_str(),k+1,flv2.c_str()));
352  h1_[Form("eta_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) + Jet %d (%s) #eta",j+1,flv1.c_str(),k+1,flv2.c_str()));
353  h1_[Form("phi_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> GetXaxis() -> SetTitle(Form("Jet %d (%s) + Jet %d (%s) #phi",j+1,flv1.c_str(),k+1,flv2.c_str()));
354  h1_[Form("m_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> GetXaxis() -> SetTitle(Form("M_{%d%d} (%s)(%s) [GeV]",j+1,k+1,flv1.c_str(),flv2.c_str()));
355  }
356  }
357  }
358  }
359  }
360  }
361 
362  this->output()->cd();
363 }
364 
365 
366 float JetAnalyser::btag(const Jet & jet, const std::string & algo)
367 {
368  float btag;
369  if ( config_->btagAlgorithm() == "csvivf" || config_->btagAlgorithm() == "csv" ) { btag = jet.btag("btag_csvivf"); }
370  else if ( config_->btagAlgorithm() == "deepcsv" ) { btag = jet.btag("btag_deepb") + jet.btag("btag_deepbb"); }
371  else if ( config_->btagAlgorithm() == "deepbvsall" ) { btag = jet.btag("btag_deepbvsall"); }
372  else if ( config_->btagAlgorithm() == "deepflavour" || config_->btagAlgorithm() == "deepjet" ) { btag = jet.btag("btag_dfb") + jet.btag("btag_dfbb") + jet.btag("btag_dflepb"); }
373  else { btag = -9999; }
374 
375  return btag;
376 }
377 
378 bool JetAnalyser::selectionJet(const int & r)
379 {
380  if ( r > config_->nJetsMin() ) return true;
381  int j = r-1;
382 
383  float pt_min = config_->jetsPtMin()[j];
384  float pt_max = -1.;
385  float eta_max = config_->jetsEtaMax()[j];
386 
387  if ( config_->jetsPtMax().size() > 0 && config_->jetsPtMax()[j] > config_->jetsPtMin()[j] ) pt_max = config_->jetsPtMax()[j];
388 
389  bool isgood = this -> selectionJet(r,pt_min,eta_max,pt_max);
390 
391  return isgood;
392 }
393 
394 bool JetAnalyser::selectionJet(const int & r, const float & pt_min, const float &eta_max, const float &pt_max)
395 {
396  if ( r > config_->nJetsMin() ) return true;
397  if ( (int)selectedJets_.size() < r ) return false; // is this correct?
398 
399  bool isgood = true;
400 
401  std::string label = Form("Jet %d: pt > %5.1f GeV and |eta| < %3.1f",r,pt_min, eta_max );
402  if ( pt_max > pt_min )
403  label = Form("Jet %d: pt > %5.1f GeV and pt < %5.1f GeV and |eta| < %3.1f",r,pt_min, pt_max, eta_max );
404 
405  int j = r-1;
406 
407  // kinematic selection
408  if ( selectedJets_[j] -> pt() < pt_min && !(pt_min < 0) ) isgood = false;
409  if ( fabs(selectedJets_[j] -> eta()) > eta_max && !(eta_max < 0) ) isgood = false;
410  if ( config_->jetsPtMax().size() > 0 )
411  {
412  if ( selectedJets_[j] -> pt() > pt_max && !(pt_max < pt_min ) ) isgood = false;
413  }
414 
415  cutflow(label,isgood);
416 
417  return isgood;
418 }
419 
420 
421 bool JetAnalyser::selectionJetDeta(const int & r1, const int & r2, const float & delta)
422 {
423  if ( r1 > config_->nJetsMin() || r2 > config_->nJetsMin() || delta == 0 ) return true;
424 
425  bool isgood = true;
426 
427  std::string label = Form("Deta(jet %d, jet %d) < %4.2f",r1,r2,fabs(delta));
428  if ( delta < 0 )
429  label = Form("Deta(jet %d, jet %d) > %4.2f",r1,r2,fabs(delta));
430 
431  int j1 = r1-1;
432  int j2 = r2-1;
433 
434  if ( delta > 0 )
435  isgood = ( fabs(selectedJets_[j1]->eta() - selectedJets_[j2]->eta()) < fabs(delta) );
436  else
437  isgood = ( fabs(selectedJets_[j1]->eta() - selectedJets_[j2]->eta()) > fabs(delta) );
438 
439  cutflow(label,isgood);
440 
441  return isgood;
442 
443 }
444 
445 bool JetAnalyser::selectionJetDeta(const int & r1, const int & r2)
446 {
447  bool ok = true;
448  if (config_->jetsDetaMax() < 0 )
449  {
450  ok = ok && true;
451  }
452  else
453  {
454  ok = ok && selectionJetDeta(r1,r2,config_->jetsDetaMax());
455  }
456 
457  if (config_->jetsDetaMin() < 0 )
458  {
459  ok = ok && true;
460  }
461  else
462  {
463  ok = ok && selectionJetDeta(r1,r2,-1*config_->jetsDetaMin());
464  }
465  return ok;
466 
467 }
468 
469 bool JetAnalyser::selectionJetDphi(const int & r1, const int & r2, const float & delta)
470 {
471  if ( r1 > config_->nJetsMin() || r2 > config_->nJetsMin() || delta == 0 ) return true;
472 
473  bool isgood = true;
474 
475  std::string label = Form("Dphi(jet %d, jet %d) < %4.2f",r1,r2,fabs(delta));
476  if ( delta < 0 )
477  label = Form("Dphi(jet %d, jet %d) > %4.2f",r1,r2,fabs(delta));
478 
479  int j1 = r1-1;
480  int j2 = r2-1;
481 
482  if ( delta > 0 )
483  isgood = ( fabs(selectedJets_[j1]->deltaPhi(*selectedJets_[j2])) < fabs(delta) );
484  else
485  isgood = ( fabs(selectedJets_[j1]->deltaPhi(*selectedJets_[j2])) > fabs(delta) );
486 
487  cutflow(label,isgood);
488 
489  return isgood;
490 
491 }
492 bool JetAnalyser::selectionJetDphi(const int & r1, const int & r2)
493 {
494  bool ok = true;
495  if (config_->jetsDphiMax() < 0 )
496  {
497  ok = ok && true;
498  }
499  else
500  {
501  ok = ok && selectionJetDphi(r1,r2,config_->jetsDphiMax());
502  }
503 
504  if (config_->jetsDphiMin() < 0 )
505  {
506  ok = ok && true;
507  }
508  else
509  {
510  ok = ok && selectionJetDphi(r1,r2,-1*config_->jetsDphiMin());
511  }
512  return ok;
513 
514 }
515 
516 
517 bool JetAnalyser::selectionJetDr(const int & r1, const int & r2, const float & delta)
518 {
519  if ( r1 > config_->nJetsMin() || r2 > config_->nJetsMin() || delta == 0 ) return true;
520 
521  bool isgood = true;
522 
523  std::string label = Form("DR(jet %d, jet %d) < %4.2f",r1,r2,fabs(delta));
524  if ( delta < 0 )
525  label = Form("DR(jet %d, jet %d) > %4.2f",r1,r2,fabs(delta));
526 
527  int j1 = r1-1;
528  int j2 = r2-1;
529 
530  if ( delta > 0 )
531  isgood = ( selectedJets_[j1]->deltaR(*selectedJets_[j2]) < fabs(delta) );
532  else
533  isgood = ( selectedJets_[j1]->deltaR(*selectedJets_[j2]) > fabs(delta) );
534 
535  cutflow(label,isgood);
536 
537  return isgood;
538 
539 }
540 
541 bool JetAnalyser::selectionJetDr(const int & r1, const int & r2)
542 {
543  bool ok = true;
544  if (config_->jetsDrMax() < 0 )
545  {
546  ok = ok && true;
547  }
548  else
549  {
550  ok = ok && selectionJetDr(r1,r2,config_->jetsDrMax());
551  }
552 
553  if (config_->jetsDrMin() < 0 )
554  {
555  ok = ok && true;
556  }
557  else
558  {
559  ok = ok && selectionJetDr(r1,r2,-1*config_->jetsDrMin());
560  }
561  return ok;
562 }
563 
564 bool JetAnalyser::selectionJetPtImbalance(const int & r1, const int & r2, const float & delta)
565 {
566  if ( ! jetsanalysis_ ) return true;
567  if ( r1 > config_->nJetsMin() || r2 > config_->nJetsMin() || delta == 0 ) return true;
568 
569  bool isgood = true;
570  std::string label = Form("DpT(jet %d, jet %d)/jet %d pT < %4.2f",r1,r2,r1,fabs(delta));
571  if ( delta < 0 )
572  label = Form("DpT(jet %d, jet %d)/jet %d pT > %4.2f",r1,r2,r1,fabs(delta));
573 
574  int j1 = r1-1;
575  int j2 = r2-1;
576 
577  if ( delta > 0 )
578  isgood = ( fabs(selectedJets_[j1]->pt() - selectedJets_[j2]->pt())/selectedJets_[j1]->pt() < fabs(delta) );
579  else
580  isgood = ( fabs(selectedJets_[j1]->pt() - selectedJets_[j2]->pt())/selectedJets_[j1]->pt() > fabs(delta) );
581 
582 
583  cutflow(label,isgood);
584 
585  return isgood;
586 
587 }
588 bool JetAnalyser::selectionJetPtImbalance(const int & r1, const int & r2)
589 {
590  bool ok = true;
591  if (config_->jetsPtImbalanceMax() < 0 )
592  {
593  ok = ok && true;
594  }
595  else
596  {
597  ok = ok && selectionJetPtImbalance(r1,r2,config_->jetsPtImbalanceMax());
598  }
599 
600  if (config_->jetsPtImbalanceMin() < 0 )
601  {
602  ok = ok && true;
603  }
604  else
605  {
606  ok = ok && selectionJetPtImbalance(r1,r2,-1*config_->jetsPtImbalanceMin());
607  }
608  return ok;
609 
610 }
611 
612 
613 std::vector< std::shared_ptr<Jet> > JetAnalyser::jets()
614 {
615  return jets_;
616 }
617 
618 std::vector< std::shared_ptr<Jet> > JetAnalyser::selectedJets()
619 {
620  return selectedJets_;
621 }
622 
623 
625 {
626  if ( ! jetsanalysis_ ) return true;
627 
628  bool isgood = true;
629  std::string label = Form("JetId: %s",config_->jetsId().c_str());
630 
631  auto jet = std::begin(selectedJets_);
632  while ( jet != std::end(selectedJets_) )
633  {
634  if ( ! (*jet)->id(config_->jetsId() ) )
635  jet = selectedJets_.erase(jet);
636  else
637  ++jet;
638  }
639  isgood = ( selectedJets_.size() > 0 );
640 
641  cutflow(label,isgood);
642 
643  return isgood;
644 }
645 
647 {
648  if ( ! jetsanalysis_ ) return true;
649 
650  bool isgood = true;
651  std::string label = Form("JetPileupId: %s",config_->jetsPuId().c_str());
652 
653  auto jet = std::begin(selectedJets_);
654  while ( jet != std::end(selectedJets_) )
655  {
656  if ( ! (*jet)->pileupJetIdFullId(config_->jetsPuId()) )
657  jet = selectedJets_.erase(jet);
658  else
659  ++jet;
660  }
661 
662  isgood = ( selectedJets_.size() > 0 );
663 
664  cutflow(label,isgood);
665 
666  return isgood;
667 
668 }
669 
671 {
672  if ( config_->nJetsMin() < 0 ) return true;
673 
674  bool isgood = true;
675  std::string label;
676 
677  if ( config_->nJetsMax() <= 0 )
678  {
679  label = Form("NJets >= %d",config_->nJetsMin());
680  isgood = ((int)selectedJets_.size() >= config_->nJetsMin());
681  }
682  else if ( config_->nJets() >= 0 )
683  {
684  label = Form("NJets = %d",config_->nJets());
685  isgood = ((int)selectedJets_.size() == config_->nJets());
686  }
687  else
688  {
689  if ( config_->nJetsMin() == 0 )
690  {
691  label = Form("NJets <= %d",config_->nJetsMax());
692  isgood = ((int)selectedJets_.size() <= config_->nJetsMax());
693  }
694  else
695  {
696  label = Form("%d <= NJets <= %d",config_->nJetsMin(),config_->nJetsMax());
697  isgood = ((int)selectedJets_.size() >= config_->nJetsMin() && (int)selectedJets_.size() <= config_->nJetsMax());
698  }
699  }
700 
701  cutflow(label,isgood);
702 
703  return isgood;
704 
705 }
706 
707 
708 bool JetAnalyser::selectionBJet(const int & r )
709 {
710  if ( config_->nJetsMin() < config_->nBJetsMin() || config_->nBJetsMin() < 1 || r > config_->nBJetsMin() || (int)(config_->jetsBtagWP()).size() < config_->nBJetsMin() ) return true;
711 
712  if ( ! config_->signalRegion() && r == config_->revBtagJet() ) return this->selectionNonBJet(r);
713 
714  int j = r-1;
715  if ( config_->btagWP(config_->jetsBtagWP()[j]) < 0 ) return true; // there is no selection here, so will not update the cutflow
716 
717  bool isgood = true;
718  std::string label = Form("Jet %d: %s btag > %6.4f (%s)",r,config_->btagAlgorithm().c_str(),config_->btagWP(config_->jetsBtagWP()[j]),config_->jetsBtagWP()[j].c_str());
719 
720  isgood = ( btag(*selectedJets_[j],config_->btagAlgorithm()) > config_->btagWP(config_->jetsBtagWP()[j]) );
721 
722  cutflow(label,isgood);
723 
724  return isgood;
725 }
726 
727 
728 bool JetAnalyser::selectionNonBJet(const int & r )
729 {
730  if ( config_->btagWP(config_->revBtagWP()) < 0 ) return true; // there is no selection here, so will not update the cutflow
731 
732  bool isgood = true;
733  std::string label = Form("Jet %d: %s btag < %6.4f (%s) [reverse btag]",r,config_->btagAlgorithm().c_str(),config_->btagWP(config_->revBtagWP()),config_->revBtagWP().c_str());
734 
735  int j = r-1;
736 
737  // jet non btag
738  isgood = ( btag(*selectedJets_[j],config_->btagAlgorithm()) < config_->btagWP(config_->revBtagWP()) );
739 
740  cutflow(label,isgood);
741 
742  return isgood;
743 }
744 
745 
747 {
748  if ( config_->triggerObjectsL1Jets() == "" && config_->triggerObjectsCaloJets() == "" && config_->triggerObjectsPFJets() == "") return true;
749  if ( config_->nJetsMin() < 0 ) return true;
750 
751  bool isgood = true;
752  std::string label = Form("Jet %d: online jet match (deltaR: L1 < %4.3f, Calo < %4.3f, PF < %4.3f)",r,config_-> triggerMatchL1JetsDrMax(),config_-> triggerMatchCaloJetsDrMax(),config_-> triggerMatchPFJetsDrMax());
753 
754  int j = r-1;
755 
756  std::string triggerObjectsL1Jets = config_->triggerObjectsL1Jets();
757  if ( config_->triggerEmulateL1Jets() != "" && config_->triggerEmulateL1JetsNMin() > 0 )
758  triggerObjectsL1Jets = config_->triggerEmulateL1Jets();
759  std::string triggerObjectsCaloJets = config_->triggerObjectsCaloJets();
760  if ( config_->triggerEmulateCaloJets() != "" && config_->triggerEmulateCaloJetsNMin() > 0 )
761  triggerObjectsCaloJets = config_->triggerEmulateCaloJets();
762  std::string triggerObjectsPFJets = config_->triggerObjectsPFJets();
763  if ( config_->triggerEmulatePFJets() != "" && config_->triggerEmulatePFJetsNMin() > 0 )
764  triggerObjectsPFJets = config_->triggerEmulatePFJets();
765 
766 
767  std::shared_ptr<Jet> jet = selectedJets_[j];
768  isgood = ( jet->matched(triggerObjectsL1Jets ) );
769  isgood = ( isgood && jet->matched(triggerObjectsCaloJets ) );
770  isgood = ( isgood && jet->matched(triggerObjectsPFJets ) );
771 
772  cutflow(label,isgood);
773 
774  return isgood;
775 }
776 
777 
779 {
780  if ( config_->triggerObjectsBJets() == "" ) return true;
781 
782  bool isgood = true;
783  std::string label = Form("Jet %d: online b jet match (deltaR < %4.3f)",r,config_-> triggerMatchCaloBJetsDrMax());
784 
785  int j = r-1;
786 
787  std::shared_ptr<Jet> jet = selectedJets_[j];
788  isgood = ( jet->matched(config_->triggerObjectsBJets()) );
789 
790  cutflow(label,isgood);
791 
792  return isgood;
793 }
794 
795 void JetAnalyser::fillJetHistograms(const std::string & label)
796 {
797  this->output()->cd();
798 
799  this->output()->cd(label.c_str());
800 
801  int n = n_hjets_;
802 
803  if ( n > config_->nJetsMin() ) n = config_->nJetsMin();
804 
805  h1_[Form("jet_hist_weight_%s",label.c_str())] -> Fill(0.,weight_);
806 
807  for ( int j = 0; j < n; ++j )
808  {
809  // fill each jet with addtional SF to weight = 1
810  this->fillJetHistograms(j+1,label,1.);
811 
812  if ( config_ -> doDijet() )
813  {
814  for ( int k = j+1; k < n && j < n; ++k )
815  {
817 
818  h1_[Form("dptrel_jet%d%d_%s" , j+1,k+1,label.c_str())] -> Fill(fabs(selectedJets_[j]->pt()-selectedJets_[k]->pt())/selectedJets_[j]->pt(),weight_);
819  h1_[Form("dpt_jet%d%d_%s" , j+1,k+1,label.c_str())] -> Fill(fabs(selectedJets_[j]->pt()-selectedJets_[k]->pt()),weight_);
820  h1_[Form("dr_jet%d%d_%s",j+1,k+1,label.c_str())] -> Fill(c_ij.deltaR(),weight_);
821  h1_[Form("deta_jet%d%d_%s",j+1,k+1,label.c_str())] -> Fill(c_ij.deltaEta(),weight_);
822  h1_[Form("dphi_jet%d%d_%s",j+1,k+1,label.c_str())] -> Fill(fabs(selectedJets_[j]->deltaPhi(*selectedJets_[k])));
823  h1_[Form("pt_jet%d%d_%s",j+1,k+1,label.c_str())] -> Fill(c_ij.pt(),weight_);
824  h1_[Form("eta_jet%d%d_%s",j+1,k+1,label.c_str())] -> Fill(c_ij.eta(),weight_);
825  h1_[Form("phi_jet%d%d_%s",j+1,k+1,label.c_str())] -> Fill(c_ij.phi()*180./acos(-1.),weight_);
826  if ( config_->isMC() || !config_->signalRegion() )
827  {
828  h1_[Form("m_jet%d%d_%s",j+1,k+1,label.c_str())] -> Fill(c_ij.m(),weight_);
829  }
830  else // blind
831  {
832  h1_[Form("m_jet%d%d_%s",j+1,k+1,label.c_str())] -> Fill(0.,weight_);
833  }
834  if ( config_->isMC() && config_->histogramJetsPerFlavour() )
835  {
836  std::string flv1 = "udsg";
837  std::string flv2 = "udsg";
838  if ( ! config_ -> useJetsExtendedFlavour() )
839  {
840  if ( abs(selectedJets_[j]->flavour()) == 4 ) flv1 = "c";
841  if ( abs(selectedJets_[j]->flavour()) == 5 ) flv1 = "b";
842  if ( abs(selectedJets_[k]->flavour()) == 4 ) flv2 = "c";
843  if ( abs(selectedJets_[k]->flavour()) == 5 ) flv2 = "b";
844  }
845  else
846  {
847  flv1 = selectedJets_[j]->extendedFlavour();
848  flv2 = selectedJets_[k]->extendedFlavour();
849  }
850  h1_[Form("dptrel_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> Fill(fabs(selectedJets_[j]->pt()-selectedJets_[k]->pt())/selectedJets_[j]->pt(),weight_);
851  h1_[Form("dpt_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> Fill(fabs(selectedJets_[j]->pt()-selectedJets_[k]->pt()),weight_);
852  h1_[Form("dr_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> Fill(c_ij.deltaR(),weight_);
853  h1_[Form("deta_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> Fill(c_ij.deltaEta(),weight_);
854  h1_[Form("dphi_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> Fill(fabs(selectedJets_[j]->deltaPhi(*selectedJets_[k])));
855  h1_[Form("pt_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> Fill(c_ij.pt(),weight_);
856  h1_[Form("eta_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> Fill(c_ij.eta(),weight_);
857  h1_[Form("phi_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> Fill(c_ij.phi()*180./acos(-1.),weight_);
858  h1_[Form("m_jet%d%d_%s_%s_%s" , j+1,k+1,label.c_str(),flv1.c_str(),flv2.c_str())] -> Fill(c_ij.m(),weight_);
859 
860  }
861  }
862  }
863  }
864  this->output()->cd();
865 
866  cutflow(Form("*** Filling jets histograms - %s",label.c_str()));
867 
868 }
869 
870 void JetAnalyser::fillJetHistograms(const int & r, const std::string & label, const float & sf, const bool & workflow)
871 {
872  if ( r < 1 || r > n_hjets_ ) return;
873 
874  if ( workflow ) // BE CAREFUL with this
875  {
876  this->output()->cd();
877  ++ cutflow_;
878  if ( std::string(h1_["cutflow"] -> GetXaxis()-> GetBinLabel(cutflow_+1)) == "" )
879  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("*** Filling jet # %d histograms - %s",r, label.c_str()));
880  this->output()->cd(label.c_str());
881  }
882 
883  int j = r-1;
884 
885  // 1D histograms
886  h1_[Form("pt_jet%d_%s",j+1,label.c_str())] -> Fill(selectedJets_[j]->pt(),weight_*sf);
887  // barrel and endcap pt distributions
888  float eta = selectedJets_[j]->eta();
889  if ( config_-> histogramJetsRegionSplit() )
890  {
891  if ( eta < -1.4 ) h1_[Form("pt_jet%d_me_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->pt(),weight_*sf);
892  if ( eta > 1.4 ) h1_[Form("pt_jet%d_pe_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->pt(),weight_*sf);
893  if ( eta < 0.0 && eta > -1.0 ) h1_[Form("pt_jet%d_mb_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->pt(),weight_*sf);
894  if ( eta > 0.0 && eta < 1.0 ) h1_[Form("pt_jet%d_pb_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->pt(),weight_*sf);
895  if ( eta < -1.0 && eta > -1.4 ) h1_[Form("pt_jet%d_mbe_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->pt(),weight_*sf);
896  if ( eta > 1.0 && eta < 1.4 ) h1_[Form("pt_jet%d_pbe_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->pt(),weight_*sf);
897  }
898 
899  //
900  h1_[Form("eta_jet%d_%s",j+1,label.c_str())] -> Fill(selectedJets_[j]->eta(),weight_*sf);
901  h1_[Form("phi_jet%d_%s",j+1,label.c_str())] -> Fill(selectedJets_[j]->phi()*180./acos(-1.),weight_*sf);
902  float mybtag = btag(*selectedJets_[j],config_->btagAlgorithm());
903  float mybtaglog = 1.e-7;
904  if ( mybtag > 0 ) mybtaglog = -log(1.-mybtag);
905  h1_[Form("btag_jet%d_%s",j+1,label.c_str())] -> Fill(mybtag,weight_*sf);
906  h1_[Form("btaglog_jet%d_%s",j+1,label.c_str())] -> Fill(mybtaglog,weight_*sf);
907  h1_[Form("qglikelihood_jet%d_%s", j+1,label.c_str())] -> Fill(selectedJets_[j]->qgLikelihood(),weight_*sf);
908  h1_[Form("nconstituents_jet%d_%s", j+1,label.c_str())] -> Fill(selectedJets_[j]->constituents(),weight_*sf);
909 
910  if ( config_->btagAlgorithm() == "deepcsv")
911  {
912  h1_[Form("btag_light_jet%d_%s", j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_deeplight"),weight_*sf);
913  h1_[Form("btag_c_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_deepc"),weight_*sf);
914  h1_[Form("btag_b_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_deepb"),weight_*sf);
915  h1_[Form("btag_bb_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_deepbb"),weight_*sf);
916  h1_[Form("btag_cc_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_deepcc"),weight_*sf);
917  }
918  if ( config_->btagAlgorithm() == "deepflavour" || config_->btagAlgorithm() == "deepjet" )
919  {
920  h1_[Form("btag_light_jet%d_%s", j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_dflight"),weight_*sf);
921  h1_[Form("btag_g_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_dfg") ,weight_*sf);
922  h1_[Form("btag_c_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_dfc") ,weight_*sf);
923  h1_[Form("btag_b_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_dfb") ,weight_*sf);
924  h1_[Form("btag_bb_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_dfbb") ,weight_*sf);
925  h1_[Form("btag_lepb_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->btag("btag_dflepb") ,weight_*sf);
926  }
927  // 2D histograms
928  h2_[Form("pt_eta_jet%d_%s" , j+1,label.c_str())] -> Fill(selectedJets_[j]->pt(), selectedJets_[j]->eta(), weight_*sf);
929 
930 
931  if ( config_->isMC() && config_->histogramJetsPerFlavour() )
932  {
933  std::string flv = "udsg";
934  if ( ! config_ -> useJetsExtendedFlavour() )
935  {
936  if ( abs(selectedJets_[j]->flavour()) == 4 ) flv = "c";
937  if ( abs(selectedJets_[j]->flavour()) == 5 ) flv = "b";
938  }
939  else
940  {
941  flv = selectedJets_[j]->extendedFlavour();
942  }
943  // 1D histograms
944  h1_[Form("pt_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->pt(),weight_*sf);
945  h1_[Form("eta_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->eta(),weight_*sf);
946  h1_[Form("phi_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->phi()*180./acos(-1.),weight_*sf);
947  h1_[Form("btag_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] -> Fill(btag(*selectedJets_[j],config_->btagAlgorithm()),weight_*sf);
948  h1_[Form("qglikelihood_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->qgLikelihood(),weight_*sf);
949  h1_[Form("nconstituents_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->constituents(),weight_*sf);
950  if ( config_->btagAlgorithm() == "deepcsv")
951  {
952  h1_[Form("btag_light_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_deeplight"),weight_*sf);
953  h1_[Form("btag_c_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_deepc"),weight_*sf);
954  h1_[Form("btag_b_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_deepb"),weight_*sf);
955  h1_[Form("btag_bb_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_deepbb"),weight_*sf);
956  h1_[Form("btag_cc_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_deepcc"),weight_*sf);
957  }
958  if ( config_->btagAlgorithm() == "deepflavour" || config_->btagAlgorithm() == "deepjet")
959  {
960  h1_[Form("btag_light_jet%d_%s_%s", j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_dflight"),weight_*sf);
961  h1_[Form("btag_g_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_dfg") ,weight_*sf);
962  h1_[Form("btag_c_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_dfc") ,weight_*sf);
963  h1_[Form("btag_b_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_dfb") ,weight_*sf);
964  h1_[Form("btag_bb_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_dfbb") ,weight_*sf);
965  h1_[Form("btag_lepb_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->btag("btag_dflepb") ,weight_*sf);
966 
967 
968  }
969  // 2D histograms
970  h2_[Form("pt_eta_jet%d_%s_%s" , j+1,label.c_str(),flv.c_str())] -> Fill(selectedJets_[j]->pt(), selectedJets_[j]->eta(), weight_*sf);
971  }
972 
973  this->output()->cd();
974  if ( workflow ) h1_["cutflow"] -> Fill(cutflow_,weight_*sf);
975 
976 }
977 
978 
979 
980 ScaleFactors JetAnalyser::btagSF(const int & r, const std::string & wp)
981 {
982  ScaleFactors sf;
983  int j = r-1;
984  sf.nominal = selectedJets_[j]->btagSF(bsf_reader_[wp]);
985  sf.up = selectedJets_[j]->btagSFup(bsf_reader_[wp]);
986  sf.down = selectedJets_[j]->btagSFdown(bsf_reader_[wp]);
987 
988  return sf;
989 }
990 
991 
993 {
994  if ( ! jetsanalysis_ || ! isMC_ ) return;
995 
996  std::string label = "WARNING: NO JER smearing (*** missing JER Info and/or GenJet collection ***)";
997 
998  if ( applyjer_ )
999  {
1000  std::string bnpt = basename(config_->jerPtRes());
1001  std::string bnsf = basename(config_->jerSF());
1002  label = Form("JER smearing (%s,%s)",bnpt.c_str(),bnsf.c_str());
1003  for ( auto & j : selectedJets_ )
1004  j -> applyJER(*jerinfo_,0.2);
1005  }
1006 
1007  cutflow(label);
1008 }
1009 
1010 float JetAnalyser::actionApplyBtagSF(const int & r, const bool & global_weight)
1011 {
1012  float sf = 1.;
1013  if ( ! config_-> isMC() || config_->btagScaleFactors() == "" ) return sf; // will not apply btag SF
1014  if ( ! config_->signalRegion() && r == config_->revBtagJet() ) return sf;
1015 
1016  int j = r-1;
1017  std::string label = Form("Jet %d: btag SF applied (%s %s WP)",r,config_->btagAlgorithm().c_str(),config_->jetsBtagWP()[j].c_str());
1018 
1019  if ( config_->jetsBtagWP()[j] == "xxx" )
1020  label = Form("Jet %d: btag SF = 1 applied (%s %s WP)",r,config_->btagAlgorithm().c_str(),config_->jetsBtagWP()[j].c_str());
1021 
1022  if ( global_weight || config_->jetsBtagWP()[j] != "xxx" ) sf = this->btagSF(r,config_->jetsBtagWP()[j]).nominal;
1023 
1024  weight_ *= sf;
1025 
1026  cutflow(label);
1027 
1028  return sf;
1029 
1030 }
1031 
1032 float JetAnalyser::getBtagSF(const int & r)
1033 {
1034  float sf = 1.;
1035  int j = r-1;
1036  if ( ! config_-> isMC() || config_->btagScaleFactors() == "" ) return sf; // will not apply btag SF
1037  if ( ! config_->signalRegion() && r == config_->revBtagJet() ) return sf;
1038 
1039  if ( config_->jetsBtagWP()[j] != "xxx" ) sf = this->btagSF(r,config_->jetsBtagWP()[j]).nominal;
1040 
1041  return sf;
1042 }
1043 
1045 {
1046  if ( ! config_->bRegression() ) return;
1047 
1048  for ( auto & j : selectedJets_ )
1049  j -> applyBjetRegression();
1050 
1051  cutflow("b jet energy regression");
1052 }
1053 
1055 {
1056  if ( ! config_->bRegression() ) return;
1057 
1058  int j = r-1;
1059  selectedJets_[j]-> applyBjetRegression();
1060 
1061  cutflow(Form("Jet %d: b jet energy regression",r));
1062 }
1063 
1064 
1065 
1066 bool JetAnalyser::selectionDiJetMass(const int & r1, const int & r2)
1067 {
1068  float min = config_->massMin();
1069  float max = config_->massMax();
1070  if ( min < 0. && max < 0. ) return true;
1071 
1072  ++cutflow_;
1073  if ( std::string(h1_["cutflow"] -> GetXaxis()-> GetBinLabel(cutflow_+1)) == "" )
1074  {
1075  std::string label = Form("M(Jet %d + Jet %d)",r1,r2);
1076  if ( min > 0. )
1077  {
1078  if ( max > 0. && max > min ) label = Form("%5.1f GeV < %s < %5.1f GeV",min,label.c_str(),max);
1079  else label = Form("%s > %5.1f GeV",label.c_str(),min);
1080  }
1081  else
1082  {
1083  label = Form("%s < %5.1f GeV",label.c_str(),max);
1084  }
1085 
1086  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,label.c_str());
1087  }
1088 
1089  int j1 = r1-1;
1090  int j2 = r2-1;
1091  Composite<Jet,Jet> c_j1j2(*(selectedJets_[j1]),*(selectedJets_[j2]));
1092  float mass = c_j1j2.m();
1093 
1094  if ( min > 0. && mass < min ) return false;
1095  if ( max > 0. && mass > max ) return false;
1096 
1097  h1_["cutflow"] -> Fill(cutflow_,weight_);
1098 
1099  return true;
1100 }
1101 
1102 
1103 
1104 void JetAnalyser::jetSwap(const int & r1, const int & r2)
1105 {
1106  if ( r1 == r2 ) return;
1107  int j1 = r1-1;
1108  int j2 = r2-1;
1109 // ++ cutflow_;
1110 // if ( std::string(h1_["cutflow"] -> GetXaxis()-> GetBinLabel(cutflow_+1)) == "" )
1111 // h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d <--> Jet %d: jets ranking was swapped ",r1,r2));
1112 
1113  auto jet1 = selectedJets_[j1];
1114  auto jet2 = selectedJets_[j2];
1115 
1116  selectedJets_[j1] = jet2;
1117  selectedJets_[j2] = jet1;
1118 
1119 // h1_["cutflow"] -> Fill(cutflow_,weight_);
1120 
1121 }
1122 
1123 
1124 
1125 bool JetAnalyser::selectionJetQGlikelihood(const int & r, const float & cut)
1126 {
1127  bool isgood = true;
1128  std::string label = Form("Jet %d Q-G likelihood < %4.2f",r,fabs(cut));
1129 
1130  int j = r-1;
1131 
1132  if ( cut < 0 )
1133  label = Form("Jet %d Q-G likelihood > %4.2f",r,fabs(cut));
1134 
1135  if ( cut > 0 )
1136  isgood = ( selectedJets_[j]->qgLikelihood() < fabs(cut) );
1137  else
1138  isgood = ( selectedJets_[j]->qgLikelihood() > fabs(cut) );
1139 
1140  cutflow(label,isgood);
1141 
1142  return isgood;
1143 
1144 }
1145 
1147 {
1148  int j = r-1;
1149  bool ok = true;
1150  if ( config_->jetsQGmax().size() == 0 || config_->jetsQGmax()[j] < 0 || (int)config_->jetsQGmax().size() < r )
1151  {
1152  ok = ok && true;
1153  }
1154  else
1155  {
1156  ok = ok && selectionJetQGlikelihood(r,config_->jetsQGmax()[j]);
1157  }
1158 
1159  if (config_->jetsQGmin().size() == 0 || config_->jetsQGmin()[j] < 0 || (int)config_->jetsQGmin().size() < r )
1160  {
1161  ok = ok && true;
1162  }
1163  else
1164  {
1165  ok = ok && selectionJetQGlikelihood(r,-1*config_->jetsQGmin()[j] );
1166  }
1167  return ok;
1168 
1169 }
1170 
1171 
1173 {
1174  if ( config_->jetsBtagProbB().size() == 0 ) return true;
1175  int j = r-1;
1176  float wp = config_->jetsBtagProbB()[j];
1177  std::string algo = config_->btagAlgorithm();
1178  if ( fabs(wp) > 1 ) return true; // there is no selection here, so will not update the cutflow
1179 
1180  ++ cutflow_;
1181  if ( std::string(h1_["cutflow"] -> GetXaxis()-> GetBinLabel(cutflow_+1)) == "" )
1182  {
1183  if ( wp > 0 )
1184  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob b < %6.4f",r,algo.c_str(),fabs(wp)));
1185  else
1186  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob b > %6.4f",r,algo.c_str(),fabs(wp)));
1187  }
1188 
1189 
1190  if ( r > config_->nBJetsMin() )
1191  {
1192  std::cout << "* warning * - JetAnalyser::selectionBJetProbB(): given jet rank > nbjetsmin. Returning false! " << std::endl;
1193  return false;
1194  }
1195 
1196  // jet btag
1197  if ( algo == "deepcsv" )
1198  {
1199  if ( wp > 0 && selectedJets_[j]->btag("btag_deepb") > fabs(wp) ) return false;
1200  if ( wp < 0 && selectedJets_[j]->btag("btag_deepb") < fabs(wp) ) return false;
1201  }
1202  if ( algo == "deepflavour" || algo == "deepjet" )
1203  {
1204  if ( wp > 0 && selectedJets_[j]->btag("btag_dfb") > fabs(wp) ) return false;
1205  if ( wp < 0 && selectedJets_[j]->btag("btag_dfb") < fabs(wp) ) return false;
1206  }
1207 
1208  h1_["cutflow"] -> Fill(cutflow_,weight_);
1209 
1210  return true;
1211 }
1212 
1213 
1214 
1215 
1217 {
1218  if ( config_->jetsBtagProbBB().size() == 0 ) return true;
1219  int j = r-1;
1220  float wp = config_->jetsBtagProbBB()[j];
1221  std::string algo = config_->btagAlgorithm();
1222  if ( fabs(wp) > 1 ) return true; // there is no selection here, so will not update the cutflow
1223 
1224  ++ cutflow_;
1225  if ( std::string(h1_["cutflow"] -> GetXaxis()-> GetBinLabel(cutflow_+1)) == "" )
1226  {
1227  if ( wp > 0 )
1228  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob bb < %6.4f",r,algo.c_str(),fabs(wp)));
1229  else
1230  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob bb > %6.4f",r,algo.c_str(),fabs(wp)));
1231  }
1232 
1233 
1234  if ( r > config_->nBJetsMin() )
1235  {
1236  std::cout << "* warning * - JetAnalyser::selectionBJetProbBB(): given jet rank > nbjetsmin. Returning false! " << std::endl;
1237  return false;
1238  }
1239 
1240  // jet btag
1241  if ( algo == "deepcsv" )
1242  {
1243  if ( wp > 0 && selectedJets_[j]->btag("btag_deepbb") > fabs(wp) ) return false;
1244  if ( wp < 0 && selectedJets_[j]->btag("btag_deepbb") < fabs(wp) ) return false;
1245  }
1246  if ( algo == "deepflavour" || algo == "deepjet" )
1247  {
1248  if ( wp > 0 && selectedJets_[j]->btag("btag_dfbb") > fabs(wp) ) return false;
1249  if ( wp < 0 && selectedJets_[j]->btag("btag_dfbb") < fabs(wp) ) return false;
1250  }
1251 
1252  h1_["cutflow"] -> Fill(cutflow_,weight_);
1253 
1254  return true;
1255 }
1256 
1258 {
1259  if ( config_->jetsBtagProbLepB().size() == 0 ) return true;
1260  int j = r-1;
1261  float wp = config_->jetsBtagProbLepB()[j];
1262  std::string algo = config_->btagAlgorithm();
1263  if ( fabs(wp) > 1 || algo == "deepcsv" ) return true; // there is no selection here, so will not update the cutflow
1264 
1265  ++ cutflow_;
1266  if ( std::string(h1_["cutflow"] -> GetXaxis()-> GetBinLabel(cutflow_+1)) == "" )
1267  {
1268  if ( wp > 0 )
1269  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob lepb < %6.4f",r,algo.c_str(),fabs(wp)));
1270  else
1271  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob lepb > %6.4f",r,algo.c_str(),fabs(wp)));
1272  }
1273 
1274 
1275  if ( r > config_->nBJetsMin() )
1276  {
1277  std::cout << "* warning * - JetAnalyser::selectionBJetProbLepB(): given jet rank > nbjetsmin. Returning false! " << std::endl;
1278  return false;
1279  }
1280 
1281  // jet btag
1282  if ( algo == "deepflavour" || algo == "deepjet" )
1283  {
1284  if ( wp > 0 && selectedJets_[j]->btag("btag_dflepb") > fabs(wp) ) return false;
1285  if ( wp < 0 && selectedJets_[j]->btag("btag_dflepb") < fabs(wp) ) return false;
1286  }
1287 
1288  h1_["cutflow"] -> Fill(cutflow_,weight_);
1289 
1290  return true;
1291 }
1292 
1293 
1294 
1296 {
1297  if ( config_->jetsBtagProbC().size() == 0 ) return true;
1298  int j = r-1;
1299  float wp = config_->jetsBtagProbC()[j];
1300  std::string algo = config_->btagAlgorithm();
1301  if ( fabs(wp) > 1 ) return true; // there is no selection here, so will not update the cutflow
1302 
1303  ++ cutflow_;
1304  if ( std::string(h1_["cutflow"] -> GetXaxis()-> GetBinLabel(cutflow_+1)) == "" )
1305  {
1306  if ( wp > 0 )
1307  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob c < %6.4f",r,algo.c_str(),fabs(wp)));
1308  else
1309  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob c > %6.4f",r,algo.c_str(),fabs(wp)));
1310  }
1311 
1312 
1313  if ( r > config_->nBJetsMin() )
1314  {
1315  std::cout << "* warning * - JetAnalyser::selectionBJetProbC(): given jet rank > nbjetsmin. Returning false! " << std::endl;
1316  return false;
1317  }
1318 
1319  // jet btag
1320  if ( algo == "deepcsv" )
1321  {
1322  if ( wp > 0 && selectedJets_[j]->btag("btag_deepc") > fabs(wp) ) return false;
1323  if ( wp < 0 && selectedJets_[j]->btag("btag_deepc") < fabs(wp) ) return false;
1324  }
1325  if ( algo == "deepflavour" || algo == "deepjet" )
1326  {
1327  if ( wp > 0 && selectedJets_[j]->btag("btag_dfc") > fabs(wp) ) return false;
1328  if ( wp < 0 && selectedJets_[j]->btag("btag_dfc") < fabs(wp) ) return false;
1329  }
1330 
1331  h1_["cutflow"] -> Fill(cutflow_,weight_);
1332 
1333  return true;
1334 }
1335 
1337 {
1338  if ( config_->jetsBtagProbG().size() == 0 ) return true;
1339  int j = r-1;
1340  float wp = config_->jetsBtagProbG()[j];
1341  std::string algo = config_->btagAlgorithm();
1342  if ( fabs(wp) > 1 || algo == "deepcsv" ) return true; // there is no selection here, so will not update the cutflow
1343 
1344  ++ cutflow_;
1345  if ( std::string(h1_["cutflow"] -> GetXaxis()-> GetBinLabel(cutflow_+1)) == "" )
1346  {
1347  if ( wp > 0 )
1348  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob g < %6.4f",r,algo.c_str(),fabs(wp)));
1349  else
1350  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob g > %6.4f",r,algo.c_str(),fabs(wp)));
1351  }
1352 
1353 
1354  if ( r > config_->nBJetsMin() )
1355  {
1356  std::cout << "* warning * - JetAnalyser::selectionBJetProbG(): given jet rank > nbjetsmin. Returning false! " << std::endl;
1357  return false;
1358  }
1359 
1360  // jet btag
1361  if ( algo == "deepflavour" || algo == "deepjet" )
1362  {
1363  if ( wp > 0 && selectedJets_[j]->btag("btag_dfg") > fabs(wp) ) return false;
1364  if ( wp < 0 && selectedJets_[j]->btag("btag_dfg") < fabs(wp) ) return false;
1365  }
1366 
1367  h1_["cutflow"] -> Fill(cutflow_,weight_);
1368 
1369  return true;
1370 }
1371 
1373 {
1374  if ( config_->jetsBtagProbLight().size() == 0 ) return true;
1375  int j = r-1;
1376  float wp = config_->jetsBtagProbLight()[j];
1377  std::string algo = config_->btagAlgorithm();
1378  if ( fabs(wp) > 1 ) return true; // there is no selection here, so will not update the cutflow
1379 
1380  ++ cutflow_;
1381  if ( std::string(h1_["cutflow"] -> GetXaxis()-> GetBinLabel(cutflow_+1)) == "" )
1382  {
1383  if ( wp > 0 )
1384  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob light < %6.4f",r,algo.c_str(),fabs(wp)));
1385  else
1386  h1_["cutflow"] -> GetXaxis()-> SetBinLabel(cutflow_+1,Form("Jet %d: %s btag prob light > %6.4f",r,algo.c_str(),fabs(wp)));
1387  }
1388 
1389 
1390  if ( r > config_->nBJetsMin() )
1391  {
1392  std::cout << "* warning * - JetAnalyser::selectionBJetProbLight(): given jet rank > nbjetsmin. Returning false! " << std::endl;
1393  return false;
1394  }
1395 
1396  // jet btag
1397  if ( algo == "deepcsv" )
1398  {
1399  if ( wp > 0 && selectedJets_[j]->btag("btag_deeplight") > fabs(wp) ) return false;
1400  if ( wp < 0 && selectedJets_[j]->btag("btag_deeplight") < fabs(wp) ) return false;
1401  }
1402  if ( algo == "deepflavour" || algo == "deepjet" )
1403  {
1404  if ( wp > 0 && selectedJets_[j]->btag("btag_dflight") > fabs(wp) ) return false;
1405  if ( wp < 0 && selectedJets_[j]->btag("btag_dflight") < fabs(wp) ) return false;
1406  }
1407 
1408  h1_["cutflow"] -> Fill(cutflow_,weight_);
1409 
1410  return true;
1411 }
1412 
1413 // float JetAnalyser::btagEfficiency(const int & r)
1414 // {
1415 // float eff = 1;
1416 // auto btageff = this -> btagEfficiencies();
1417 // int j = r-1;
1418 //
1419 // return eff;
1420 //
1421 // }
1422 
1423 // void JetAnalyser::actionApplyBtagEfficiencyWeight(const int & r)
1424 // {
1425 // weight_ *= btagEfficiency(r);
1426 //
1427 // }
1428 
1430 {
1431 // CORRECTIONS
1432  // b energy regression
1433  if ( this->config()->bRegression() ) this->actionApplyBjetRegression();
1434  // Jet energy resolution smearing
1435  this->actionApplyJER();
1436 
1437  return true;
1438 
1439 }
virtual bool selectionBJetProbLight(const int &)
virtual bool selectionJetDphi(const int &r1, const int &r2)
Given the rankings r1 and r2 of two jets, it returns whether the jets passes the delta_phi selection;...
Definition: JetAnalyser.cc:492
virtual bool selectionNonBJet(const int &)
Definition: JetAnalyser.cc:728
int cutflow()
get cutflow index
float eta() const
returns the pseudorapidity
Definition: Candidate.cc:134
virtual void jetSwap(const int &, const int &)
virtual bool selectionBJetProbB(const int &)
std::shared_ptr< Config > config()
returns pointer to Config object
std::shared_ptr< TFile > output()
returns pointer to output root file
float deltaR() const
Definition: Composite.cc:46
virtual float getBtagSF(const int &)
float m() const
returns the mass
Definition: Candidate.cc:137
virtual bool analysisWithJets()
Definition: JetAnalyser.cc:72
virtual ScaleFactors btagSF(const int &, const std::string &)
Definition: JetAnalyser.cc:980
std::vector< std::shared_ptr< Jet > > jets_
Definition: JetAnalyser.h:46
bool genParticlesAnalysis() const
returns whether analysis of gen particles can be done
virtual bool selectionJet(const int &r)
Given the ranking &#39;r&#39; of a jet, it returns whether the jet passes the pt_min and |eta_max|, optionally pt_max, where the values of the thresholds pt_min and |eta_max|, pt_max are passed by configuration file.
Definition: JetAnalyser.cc:378
bool isMC_
flag for MC sample
Definition: BaseAnalyser.h:102
virtual bool selectionJetQGlikelihood(const int &, const float &)
std::map< std::string, std::shared_ptr< BTagCalibrationReader > > bsf_reader_
Definition: JetAnalyser.h:55
virtual void fillJetHistograms(const std::string &label="x")
Fill the pre-defined histograms created by the jetHistograms() method.
Definition: JetAnalyser.cc:795
bool genJetsAnalysis() const
returns whether analysis of gen jets can be done
virtual bool selectionBJetProbBB(const int &)
std::string basename(const std::string &)
returns the basename of a path
std::vector< std::shared_ptr< Jet > > jets()
vector of pointers of all jets from the "Jets" collection
Definition: JetAnalyser.cc:613
float deltaEta() const
Definition: Composite.cc:52
virtual bool onlineJetMatching(const int &)
Definition: JetAnalyser.cc:746
std::shared_ptr< Config > config_
Config objects.
Definition: BaseAnalyser.h:63
float phi() const
returns the azimuthal angle
Definition: Candidate.cc:135
std::shared_ptr< Analysis > analysis_
Analysis objects.
Definition: BaseAnalyser.h:61
float btag(const Jet &jet, const std::string &algo)
Returns the btag value of the jet for a given btag algorithm.
Definition: JetAnalyser.cc:366
virtual bool onlineBJetMatching(const int &)
Definition: JetAnalyser.cc:778
virtual bool selectionBJetProbC(const int &)
virtual bool selectionJetPtImbalance(const int &r1, const int &r2)
Given the rankings r1 and r2 of two jets, it returns whether the jets passes the pt imbalance selecti...
Definition: JetAnalyser.cc:588
float weight_
event weight
Definition: BaseAnalyser.h:72
virtual bool selectionBJet(const int &)
Definition: JetAnalyser.cc:708
virtual void jetHistograms(const int &n, const std::string &label="x")
Creates pre-defined histograms in directory &#39;label&#39; for analysis with &#39;n&#39; jets.
Definition: JetAnalyser.cc:172
virtual bool selectionJetDeta(const int &r1, const int &r2)
Given the rankings r1 and r2 of two jets, it returns whether the jets passes the delta_eta selection;...
Definition: JetAnalyser.cc:445
float pt() const
returns the transverse momentum
Definition: Candidate.cc:133
std::vector< std::shared_ptr< Jet > > selectedJets_
Definition: JetAnalyser.h:47
std::map< std::string, std::shared_ptr< TH1F > > h1_
1D histograms mapping
Definition: BaseAnalyser.h:80
int cutflow_
Cutflow index.
Definition: BaseAnalyser.h:66
virtual bool selectionJetDr(const int &r1, const int &r2)
Given the rankings r1 and r2 of two jets, it returns whether the jets passes the delta_R selection; t...
Definition: JetAnalyser.cc:541
virtual bool selectionDiJetMass(const int &, const int &)
virtual float actionApplyBtagSF(const int &, const bool &global_weight=true)
std::shared_ptr< JetResolutionInfo > jerinfo_
Definition: JetAnalyser.h:57
virtual bool jetCorrections()
multiple actions: apply JER and b-tag regression corrections
float btag(const std::string &) const
returns the btag value of btag_csvivf
Definition: Jet.cc:57
virtual void actionApplyBjetRegression()
virtual bool triggerEmulation(const std::string &, const int &, const float &, const float &, const std::string &)
emulate l1 trigger
virtual bool selectionBJetProbG(const int &)
virtual bool selectionJetPileupId()
Definition: JetAnalyser.cc:646
std::vector< std::string > flavours_
Definition: JetAnalyser.h:59
virtual bool selectionBJetProbLepB(const int &)
std::vector< std::shared_ptr< Jet > > selectedJets()
vector of pointers of the selectedJets
Definition: JetAnalyser.cc:618
std::map< std::string, std::shared_ptr< TH2F > > h2_
2D histograms mapping
Definition: BaseAnalyser.h:82