CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
JetCombinatorics.h
Go to the documentation of this file.
1 #ifndef JetCombinatorics_h
2 #define JetCombinatorics_h
3 
15 #include "TLorentzVector.h"
16 #include "TString.h"
17 #include "TH1F.h"
18 #include "TFile.h"
19 #include "TMath.h"
20 #include <map>
21 #include <vector>
22 #include <iostream>
23 
24 class Combo {
25 
26 
27  public:
28 
29  Combo() {
30 
31  MW = 84.2;//79.8;
32  Mtop_h = 180.7;//175.;
33  Mtop_l = 174.9;
34  sigmaHadW = 10.5;//2.*7.6;
35  sigmaHadt = 19.2;//2.*12.5;
36  sigmaLept = 24.2;//2.*15.6;
37 
38  SumEt_ = 0.;
39  usebtag_ = false;
40  useMtop_ = true;
41 
42  useFlv_ = false;
44  }
45  ~Combo(){};
46 
47  void SetWp(const TLorentzVector& Wp) { Wp_ = Wp; }
48  void SetWq(const TLorentzVector& Wq) { Wq_ = Wq; }
49  void SetHadb(const TLorentzVector& Hadb) { Hadb_ = Hadb; }
50  void SetLepW(const TLorentzVector& LepW) { LepW_ = LepW; }
51  void SetLepb(const TLorentzVector& Lepb) { Lepb_ = Lepb; }
52  // flavor corrections
53  void ApplyFlavorCorrections(bool option=true){ useFlv_ = option;}
54  void SetFlvCorrWp( double corr ) { Wp_flv_ = corr; }
55  void SetFlvCorrWq( double corr ) { Wq_flv_ = corr; }
56  void SetFlvCorrHadb( double corr ) { Hadb_flv_ = corr; }
57  void SetFlvCorrLepb( double corr ) { Lepb_flv_ = corr; }
58  // b tagging
59  void SetWp_disc(double disc) { Wp_disc_ = disc;}
60  void SetWq_disc(double disc) { Wq_disc_= disc;}
61  void SetHadb_disc(double disc) { Hadb_disc_= disc;}
62  void SetLepb_disc(double disc) { Lepb_disc_= disc;}
63  void SetbDiscPdf(const TString& filename) {
64  pdffile_ = TFile::Open(filename);
65  hdisc_b_ = (TH1F*) gDirectory->Get("hdiscNorm_b");
66  hdisc_cl_ = (TH1F*) gDirectory->Get("hdiscNorm_cl");
67  }
68  void SetSigmas(int type=0) {
69 
70  // type == 0 take defaults
71  if (type==1) {
72  // JES +10%
73  MW = 87.2;
74  Mtop_h = 193.2;
75  Mtop_l = 179.0;
76  sigmaHadW = 13.0;
77  sigmaHadt = 22.8;
78  sigmaLept = 26.3;
79  }
80  if (type==-1) {
81  // JES -10%
82  MW = 81.6;
83  Mtop_h = 169.3;
84  Mtop_l = 171.4;
85  sigmaHadW =8.9;
86  sigmaHadt =17.9;
87  sigmaLept =22.6;
88  }
89 
90  }
91  void Usebtagging(bool option = true) { usebtag_ = option;}
92  void SetMinMassLepW( double mass ) { minMassLepW_ = mass; }
93  void SetMaxMassLepW( double mass ) { maxMassLepW_ = mass; }
94  void SetMinMassHadW( double mass ) { minMassHadW_ = mass; }
95  void SetMaxMassHadW( double mass ) { maxMassHadW_ = mass; }
96  void SetMinMassLepTop( double mass ) { minMassLepTop_ = mass; }
97  void SetMaxMassLepTop( double mass ) { maxMassLepTop_ = mass; }
98  void UseMtopConstraint(bool option=true) { useMtop_ = option; }
99 
100  void analyze() {
101 
102  if ( useFlv_ ) {
103  Wp_ = Wp_flv_ * Wp_;
104  Wq_ = Wq_flv_ * Wq_;
105  Hadb_ = Hadb_flv_ * Hadb_;
106  Lepb_ = Lepb_flv_ * Lepb_;
107  }
108 
109  HadW_ = Wp_ + Wq_;
110  HadTop_ = HadW_ + Hadb_;
111  LepTop_ = LepW_ + Lepb_;
113 
114  //double sigmaHadW = 10.5;//2.*7.6;
115  //double sigmaHadt = 19.2;//2.*12.5;
116  //double sigmaLept = 24.2;//2.*15.6;
117 
118  double chiHadW = (HadW_.M() - MW)/sigmaHadW;
119  double chiHadt = (HadTop_.M() - Mtop_h)/sigmaHadt;
120  double chiLept = (LepTop_.M() - Mtop_l)/sigmaLept;
121 
122  if ( useMtop_ ) {
123  chi2_ = chiHadW*chiHadW + chiHadt*chiHadt + chiLept*chiLept;
124  Ndof_ = 3;
125  } else {
126  chi2_ = chiHadW*chiHadW + (HadTop_.M() - LepTop_.M())*(HadTop_.M() - LepTop_.M())/(sigmaHadt*sigmaHadt+sigmaLept*sigmaLept);
127  Ndof_ = 2;
128  }
129 
130  SumEt_ = HadTop_.Pt();
131 
132  if ( usebtag_ ) {
133 
134  double gauss_norm = (2.)*TMath::Log(sigmaHadW*TMath::Sqrt(2*TMath::Pi())) +
135  (2.)*TMath::Log(sigmaHadt*TMath::Sqrt(2*TMath::Pi())) + (2.)*TMath::Log(sigmaLept*TMath::Sqrt(2*TMath::Pi()));
136 
137  double LR_Wp; double LR_Wq;
138  double LR_Hadb; double LR_Lepb;
139 
140  double LR_den = 0;
141  LR_den = ( getPdfValue("cl", Wp_disc_) + getPdfValue("b", Wp_disc_));
142  if (LR_den == 0 ) LR_Wp = 1e-5;
143  else LR_Wp = getPdfValue( "cl", Wp_disc_ )/ LR_den;
144 
145  LR_den = ( getPdfValue("cl", Wq_disc_) + getPdfValue("b", Wq_disc_));
146  if (LR_den == 0 ) LR_Wq = 1e-5;
147  else LR_Wq = getPdfValue( "cl", Wq_disc_ )/ LR_den;
148 
149  LR_den = ( getPdfValue("cl", Hadb_disc_) + getPdfValue("b", Hadb_disc_));
150  if (LR_den == 0 ) LR_Hadb = 1e-5;
151  else LR_Hadb = getPdfValue( "b", Hadb_disc_ )/ LR_den;
152 
153  LR_den = ( getPdfValue("cl", Lepb_disc_) + getPdfValue("b", Lepb_disc_));
154  if (LR_den == 0 ) LR_Lepb = 1e-5;
155  else LR_Lepb = getPdfValue( "b", Lepb_disc_ )/ LR_den;
156 
157  double btag_norm = (-0.25-TMath::Log(4)/2);
158  double btag_N2LL = btag_norm*4.*( LR_Wp * TMath::Log(LR_Wp/4) + LR_Wq*TMath::Log(LR_Wq/4) + LR_Hadb*TMath::Log(LR_Hadb/4) + LR_Lepb*TMath::Log(LR_Lepb/4) );
159 
160  chi2_ += btag_N2LL + gauss_norm;
161  Ndof_ += 3;
162  pdffile_->Close();
163  }
164  }
165 
166  TLorentzVector GetWp() { return Wp_; }
167  TLorentzVector GetWq() { return Wq_; }
168  TLorentzVector GetHadW() { return HadW_; }
169  TLorentzVector GetLepW() { return LepW_; }
170  TLorentzVector GetHadb() { return Hadb_; }
171  TLorentzVector GetLepb() { return Lepb_; }
172  TLorentzVector GetHadTop() { return HadTop_; }
173  TLorentzVector GetLepTop() { return LepTop_; }
174  TLorentzVector GetTopPair() { return TopPair_; }
175  double GetChi2() const { return chi2_; }
176  double GetNdof() { return Ndof_; }
177  double GetSumEt() const { return SumEt_; }
178  int GetIdHadb() { return IdHadb_;}
179  int GetIdWp() { return IdWp_; }
180  int GetIdWq() { return IdWq_; }
181  int GetIdLepb() { return IdLepb_;}
182  void SetIdHadb(int id) { IdHadb_ = id;}
183  void SetIdWp(int id) { IdWp_ = id; }
184  void SetIdWq(int id) { IdWq_ = id; }
185  void SetIdLepb(int id) { IdLepb_ = id;}
186  void Print() {
187  std::cout << " jet Wp : px = " << Wp_.Px() << " py = " << Wp_.Py() << " pz = " << Wp_.Pz() << " e = " << Wp_.E() << std::endl;
188  std::cout << " jet Wq : px = " << Wq_.Px() << " py = " << Wq_.Py() << " pz = " << Wq_.Pz() << " e = "<< Wq_.E() << std::endl;
189  std::cout << " jet Hadb: px = " << Hadb_.Px() << " py = " << Hadb_.Py() <<" pz = " << Hadb_.Pz() <<" e = "<< Hadb_.E() << std::endl;
190  std::cout << " jet Lepb: px = " << Lepb_.Px() << " py = " << Lepb_.Py() <<" pz = " << Lepb_.Pz() <<" e = "<< Lepb_.E() << std::endl;
191  std::cout << " chi-squared = " << chi2_ << " sumEt = " << SumEt_ << std::endl;
192  }
193  double getPdfValue(std::string flavor, double disc) {
194  double pdf= 0;
195  TH1F *hpdf;
196  if ( flavor == "b" ) hpdf = hdisc_b_;
197  else hpdf = hdisc_cl_;
198  int bin = hpdf->GetXaxis()->FindBin( disc );
199  pdf = hpdf->GetBinContent( bin );
200  if ( disc < -10 || disc >50 ) return 0;
201  //if ( pdf == 0 ) return 1.e-7;
202  return pdf;
203  }
204 
205  private:
206 
207  TLorentzVector Wp_;
208  TLorentzVector Wq_;
209  TLorentzVector HadW_;
210  TLorentzVector Hadb_;
211  TLorentzVector HadTop_;
212  TLorentzVector LepW_;
213  TLorentzVector Lepb_;
214  TLorentzVector LepTop_;
215  TLorentzVector TopPair_;
216 
217  bool usebtag_;
218  bool useMtop_;
219  double Wp_disc_;
220  double Wq_disc_;
221  double Hadb_disc_;
222  double Lepb_disc_;
223  TFile *pdffile_;
224  TH1F *hdisc_b_;
225  TH1F *hdisc_cl_;
226 
228  bool useFlv_;
229  double chi2_;
230  double Ndof_;
231  double SumEt_;
232  double minMassLepW_;
233  double maxMassLepW_;
234  double minMassHadW_;
235  double maxMassHadW_;
236 
239 
240  double MW;
241  double Mtop_h;
242  double Mtop_l;
243  double sigmaHadW;
244  double sigmaHadt;
245  double sigmaLept;
246 
247 
248  int IdHadb_;
249  int IdWp_;
250  int IdWq_;
251  int IdLepb_;
252 
253 };
254 
255 struct minChi2
256 {
257  bool operator()(const Combo& s1, const Combo& s2) const
258  {
259  return s1.GetChi2() <= s2.GetChi2();
260  }
261 };
262 
263 struct maxSumEt
264 {
265  bool operator()(const Combo& s1,const Combo& s2) const
266  {
267  return s1.GetSumEt() >= s2.GetSumEt();
268  }
269 };
270 
271 
273 
274  public:
275 
278 
279  void Verbose() {
280  verbosef = true;
281  }
282 
283  std::map< int, std::string > Combinatorics(int k, int max = 6);
284  std::map< int, std::string > NestedCombinatorics();
285 
286  void FourJetsCombinations(const std::vector<TLorentzVector>& jets, const std::vector<double>& bdiscriminators );
287  void SetFlavorCorrections(const std::vector<double >& vector ) { flavorCorrections_ = vector; }
288  void SetMaxNJets(int n) { maxNJets_ = n; }
289  Combo GetCombination(int n=0);
290  Combo GetCombinationSumEt(int n=0);
291  int GetNumberOfCombos() { return ( (int)allCombos_.size() ); }
292  //void SetCandidate( std::vector< TLorentzVector > JetCandidates );
293 
294  void SetSigmas(int type = 0) {
295  SigmasTypef = type;
296  }
297  void SetLeptonicW( const TLorentzVector& LepW ) { theLepW_ = LepW; }
298 
299  void SetMinMassLepW( double mass ) { minMassLepW_ = mass; }
300  void SetMaxMassLepW( double mass ) { maxMassLepW_ = mass; }
301  void SetMinMassHadW( double mass ) { minMassHadW_ = mass; }
302  void SetMaxMassHadW( double mass ) { maxMassHadW_ = mass; }
303  void SetMinMassLepTop( double mass ) { minMassLepTop_ = mass; }
304  void SetMaxMassLepTop( double mass ) { maxMassLepTop_ = mass; }
305 
306  void UsebTagging( bool option = true ) { UsebTagging_ = option; }
307  void ApplyFlavorCorrection( bool option = true ) { UseFlv_ = option; }
308  void UseMtopConstraint( bool option = true) { UseMtop_ = option; }
309  void SetbTagPdf(const TString& name ) { bTagPdffilename_ = name; }
310  void Clear();
311 
312  std::vector< TLorentzVector > TwoCombos();
313  std::vector< TLorentzVector > ThreeCombos();
314 
315  void RemoveDuplicates( bool option) { removeDuplicates_ = option; }
316 
317  std::vector< TLorentzVector > GetComposites();
318  void AnalyzeCombos();
319 
320 
321  private:
322 
323  //int kcombos_;
324  //int maxcombos_;
326  bool verbosef;
327  std::map< int, std::string > Template4jCombos_;
328  std::map< int, std::string > Template5jCombos_;
329  std::map< int, std::string > Template6jCombos_;
330  std::map< int, std::string > Template7jCombos_;
331 
332  std::vector< double > flavorCorrections_;
335  bool UseMtop_;
337  bool UseFlv_;
338 
339  TLorentzVector theLepW_;
340 
341  double minMassLepW_;
342  double maxMassLepW_;
343  double minMassHadW_;
344  double maxMassHadW_;
347 
348  std::map< Combo, int, minChi2 > allCombos_;
349  std::map< Combo, int, maxSumEt > allCombosSumEt_;
350 
351  Double_t minPhi_;
352  double chi2_;
353  int ndf_;
355 
356  std::vector< TLorentzVector > cand1_;
357  std::vector< TLorentzVector > cand2_;
358  std::vector< TLorentzVector > cand3_;
359 
360  //int nLists_;
361 
362  //std::vector< TLorentzVector > composites_;
363 
364 };
365 
366 #endif
double Wq_flv_
const double Pi
type
Definition: HCALResponse.h:21
double maxMassHadW_
void SetMinMassHadW(double mass)
void UseMtopConstraint(bool option=true)
void SetIdHadb(int id)
bool operator()(const Combo &s1, const Combo &s2) const
TLorentzVector Wq_
void SetLeptonicW(const TLorentzVector &LepW)
double Wp_disc_
TH1F * hdisc_b_
static const std::string LepW
TFile * pdffile_
double GetSumEt() const
void UsebTagging(bool option=true)
std::vector< double > flavorCorrections_
void FourJetsCombinations(const std::vector< TLorentzVector > &jets, const std::vector< double > &bdiscriminators)
TLorentzVector TopPair_
TLorentzVector GetHadb()
double minMassLepW_
double Wp_flv_
void SetLepb_disc(double disc)
TLorentzVector Wp_
void SetMaxMassLepW(double mass)
std::map< int, std::string > Template5jCombos_
double sigmaHadt
void SetbDiscPdf(const TString &filename)
void SetIdWp(int id)
void SetMinMassLepTop(double mass)
double GetChi2() const
std::map< int, std::string > Combinatorics(int k, int max=6)
bool operator()(const Combo &s1, const Combo &s2) const
TLorentzVector GetTopPair()
void ApplyFlavorCorrections(bool option=true)
double Hadb_flv_
std::vector< TLorentzVector > GetComposites()
std::map< Combo, int, minChi2 > allCombos_
std::vector< TLorentzVector > ThreeCombos()
TLorentzVector GetWp()
tuple s2
Definition: indexGen.py:106
void SetFlvCorrLepb(double corr)
void UseMtopConstraint(bool option=true)
std::vector< TLorentzVector > TwoCombos()
void SetFlavorCorrections(const std::vector< double > &vector)
std::vector< TLorentzVector > cand3_
double maxMassLepW_
double minMassHadW_
void SetWp_disc(double disc)
TLorentzVector Lepb_
int GetIdHadb()
void ApplyFlavorCorrection(bool option=true)
double Ndof_
double Wq_disc_
void SetMaxMassLepW(double mass)
void SetHadb_disc(double disc)
void SetFlvCorrWq(double corr)
double minMassLepTop_
std::vector< TLorentzVector > cand1_
std::map< int, std::string > NestedCombinatorics()
void SetLepb(const TLorentzVector &Lepb)
bool useFlv_
void SetWq(const TLorentzVector &Wq)
double GetNdof()
TLorentzVector HadW_
void SetMaxMassHadW(double mass)
vector< PseudoJet > jets
TLorentzVector GetHadTop()
Combo GetCombination(int n=0)
int GetIdLepb()
double Lepb_disc_
TLorentzVector LepW_
void RemoveDuplicates(bool option)
void Print()
double Mtop_h
void SetIdWq(int id)
void SetHadb(const TLorentzVector &Hadb)
void SetIdLepb(int id)
void SetMinMassLepW(double mass)
double Mtop_l
bool usebtag_
JetCorrectorParameters corr
Definition: classes.h:5
TLorentzVector theLepW_
void SetMinMassLepTop(double mass)
std::map< Combo, int, maxSumEt > allCombosSumEt_
std::map< int, std::string > Template4jCombos_
void SetFlvCorrWp(double corr)
double Lepb_flv_
TLorentzVector GetHadW()
Combo GetCombinationSumEt(int n=0)
void SetMaxMassLepTop(double mass)
void SetMaxMassLepTop(double mass)
double sigmaHadW
TH1F * hdisc_cl_
TLorentzVector GetLepb()
void SetMaxNJets(int n)
double getPdfValue(std::string flavor, double disc)
void SetWq_disc(double disc)
void SetbTagPdf(const TString &name)
int GetIdWq()
double sigmaLept
TLorentzVector HadTop_
TLorentzVector Hadb_
TLorentzVector GetWq()
TLorentzVector GetLepW()
void SetFlvCorrHadb(double corr)
void SetMinMassLepW(double mass)
void SetLepW(const TLorentzVector &LepW)
TLorentzVector GetLepTop()
void analyze()
bool useMtop_
tuple filename
Definition: lut2db_cfg.py:20
void Usebtagging(bool option=true)
void SetMinMassHadW(double mass)
double maxMassLepTop_
double SumEt_
std::map< int, std::string > Template7jCombos_
TLorentzVector LepTop_
tuple cout
Definition: gather_cfg.py:121
double MW
void SetSigmas(int type=0)
int GetIdWp()
std::vector< TLorentzVector > cand2_
std::map< int, std::string > Template6jCombos_
void SetSigmas(int type=0)
void SetWp(const TLorentzVector &Wp)
void SetMaxMassHadW(double mass)
double chi2_
double Hadb_disc_