tom is hosted by Hepforge, IPPP Durham
Template Overlap Method  0.3.0
A Template Matching Tool for Jet Substructure
TemplateBuilder.hh
1 #ifndef MATCHING_H
2 #define MATCHING_H
3 
4 
5 #include <string>
6 #include <vector>
7 #include <iostream>
8 #include <fstream>
9 #include <sstream>
10 
11 #include <iomanip>
12 
13 
14 #include "fastjet/PseudoJet.hh"
15 #include "fastjet/ClusterSequence.hh"
16 
17 enum TemplateModel {TOP, HIGGS2, HIGGS3};
18 
19 std::ostream & operator<<(std::ostream & ostr, const fastjet::PseudoJet & jet) {
20 
21  ostr << " " << std::setw(7) << jet.px()
22  << " " << std::setw(7) << jet.py()
23  << " " << std::setw(7) << jet.pz()
24  << " " << std::setw(7) << jet.e() ;
25  return ostr;
26 
27 
28 }
29 
30 double operator *(const fastjet::PseudoJet & jet1, const fastjet::PseudoJet & jet2) {
31 
32  return dot_product(jet1,jet2);
33 
34 
35 }
36 
37 
38 void TemplateBuilder (std::ofstream & File,
39  const fastjet::PseudoJet & axis,
40  const double etaMax,
41  const double phiMax,
42  const double minPt,
43  const int nEta,
44  const int nPhi,
45  const int nPt,
46  const double R,
47  TemplateModel mode
48  )
49 {
50 
51  fastjet::Strategy strategy = fastjet::Best;
52  fastjet::RecombinationScheme recombScheme = fastjet::E_scheme;
53  fastjet::JetDefinition *jetDef = NULL;
54  jetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, R,
55  recombScheme, strategy);
56 
57  std::vector <fastjet::PseudoJet> fjInputs;
58 
59 
60  int nJets = 0;
61 
62  const double mW = 80.;
63  const double Mass = axis.m();
64 
65  const double pTmin = 0., pTmax = axis.perp();
66 
67 
68  switch(mode){
69  case HIGGS2:
70  for( int iEtaNow2 = 0; iEtaNow2 < nEta; ++iEtaNow2 ) {
71  for( int iPhiNow2 = 0; iPhiNow2 < nPhi; ++iPhiNow2 ) {
72  double etaCell2 = (etaMax / nEta) * (2 * iEtaNow2 - 1 - nEta);
73  double phiCell2 = (phiMax / nPhi) * (2 * iPhiNow2 - 1 - nPhi);
74  double pTcell2 = Mass * Mass / ( 2 * (axis.e() * cosh(etaCell2) - axis.perp() * cos(phiCell2)) );
75  fastjet::PseudoJet p2 = fastjet::PtYPhiM(pTcell2, etaCell2, phiCell2);
76  fastjet::PseudoJet p1 = axis - p2;
77  fjInputs.resize(0);
78  if (p1.perp() < minPt) continue;
79  if (p2.perp() < minPt) continue;
80  fjInputs.push_back(p1);
81  fjInputs.push_back(p2);
82  std::vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
83  fastjet::ClusterSequence clustSeq(fjInputs, *jetDef);
84  inclusiveJets = clustSeq.inclusive_jets();
85  sortedJets = sorted_by_pt(inclusiveJets);
86  if (sortedJets.size()>1) continue;
87  File << int(2) <<std::endl;
88  File << axis << std::endl;
89  File << p1 << std::endl;
90  File << p2 << std::endl;
91  ++nJets;}}
92  std::cout << "Found "<<nJets << std::endl;
93  break;
94  case HIGGS3:
95  for( int ipTnow = 0; ipTnow < nPt; ++ipTnow ) {
96  for( int iEtaNow1 = 0; iEtaNow1 < nEta; ++iEtaNow1 ) {
97  for( int iPhiNow1 = 0; iPhiNow1 < nPhi; ++iPhiNow1 ) {
98  for( int iEtaNow2 = 0; iEtaNow2 < nEta; ++iEtaNow2 ) {
99  for( int iPhiNow2 = 0; iPhiNow2 < nPhi; ++iPhiNow2 ) {
100  double etaCell1 = (etaMax / nEta) * (2 * iEtaNow1 - 1 - nEta);
101  double phiCell1 = (phiMax / nPhi) * (2 * iPhiNow1 - 1 - nPhi);
102  double etaCell2 = (etaMax / nEta) * (2 * iEtaNow2 - 1 - nEta);
103  double phiCell2 = (phiMax / nPhi) * (2 * iPhiNow2 - 1 - nPhi);
104  double pTcell1 = ( (pTmax - pTmin)/nPt) * ipTnow + pTmin;
105  fastjet::PseudoJet p1 = fastjet::PtYPhiM( pTcell1, etaCell1, phiCell1);
106  double p1P =axis * p1;
107  double pTcell2 = (Mass * Mass - 2 * p1P) / ( 2 * pTcell1 * (cos(phiCell1 -phiCell2) - cosh(etaCell1 - etaCell2) )
108  + 2 * (axis.e() * cosh(etaCell2) - axis.perp() * cos(phiCell2)) );
109  fastjet::PseudoJet p2 = fastjet::PtYPhiM(pTcell2, etaCell2, phiCell2);
110  fastjet::PseudoJet p3 = axis - p1 - p2;
111 
112  // Reset Fastjet input
113  fjInputs.resize(0);
114  //min pt cut on momenta
115  if (p1.perp() < minPt) continue;
116  if (p2.perp() < minPt) continue;
117  if (p3.perp() < minPt) continue;
118 
119  if (p1.m() < -0.0001) continue;
120  if (p2.m() < -0.0001) continue;
121  if (p3.m() < -0.0001) continue;
122 
123  if (p1.e() < 0.) continue;
124  if (p2.e() < 0.) continue;
125  if (p3.e() < 0.) continue;
126 
127  if (p1.perp() <p2.perp() || p1.perp() < p3.perp() ) continue;
128 
129  fjInputs.push_back(p1);
130  fjInputs.push_back(p2);
131  fjInputs.push_back(p3);
132 
133  // Run Fastjet algorithm
134  std::vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
135  fastjet::ClusterSequence clustSeq(fjInputs, *jetDef);
136  inclusiveJets = clustSeq.inclusive_jets();
137  sortedJets = sorted_by_pt(inclusiveJets);
138 
139  if (sortedJets.size()>1) continue;
140  File << int(3) <<std::endl;
141 
142  File << axis;
143  File << p1 << std::endl;
144  File << p2 << std::endl;
145  File << p3 << std::endl;
146  ++nJets;
147  }}}}}
148  std::cout << "Found "<<nJets << std::endl;
149 
150  break;
151  case TOP:
152  for( int iEtaNow1 = 0; iEtaNow1 < nEta; ++iEtaNow1 ) {
153  for( int iPhiNow1 = 0; iPhiNow1 < nPhi; ++iPhiNow1 ) {
154  for( int iEtaNow2 = 0; iEtaNow2 < nEta; ++iEtaNow2 ) {
155  for( int iPhiNow2 = 0; iPhiNow2 < nPhi; ++iPhiNow2 ) {
156  double etaCell1 = (etaMax / nEta) * (2 * iEtaNow1 - 1 - nEta);
157  double phiCell1 = (phiMax / nPhi) * (2 * iPhiNow1 - 1 - nPhi);
158  double etaCell2 = (etaMax / nEta) * (2 * iEtaNow2 - 1 - nEta);
159  double phiCell2 = (phiMax / nPhi) * (2 * iPhiNow2 - 1 - nPhi);
160  fastjet::PseudoJet e1 = fastjet::PtYPhiM(1., etaCell1, phiCell1);
161  fastjet::PseudoJet e2 = fastjet::PtYPhiM(1., etaCell2, phiCell2);
162  double p1P = axis * e1;
163  double p2P = axis * e2;
164  double p12 = e1 * e2;
165  double Delta = 1- (8* mW*mW*p1P*p2P)/((Mass*Mass+mW*mW)*(Mass*Mass+mW*mW) *p12);
166  if (Delta <0.) continue;
167  double pT1 = ( (Mass * Mass + mW * mW) / (4 * p1P ) ) *( 1+ sqrt(Delta ));
168  double pT2 = mW*mW/(2 * pT1 *p12);
169  fastjet::PseudoJet p1 = pT1 * e1;
170  fastjet::PseudoJet p2 = pT2 * e2;
171  fastjet::PseudoJet p3 = axis - p1 - p2;
172  // Reset Fastjet input
173  fjInputs.resize(0);
174  if (p1.m() < -0.0001) continue;
175  if (p2.m() < -0.0001) continue;
176  if (p3.m() < -0.0001) continue;
177  if (p1.e() < 0.) continue;
178  if (p2.e() < 0.) continue;
179  if (p3.e() < 0.) continue;
180  if (p1.perp() < minPt) continue;
181  if (p2.perp() < minPt) continue;
182  if (p3.perp() < minPt) continue;
183  fjInputs.push_back(p1);
184  fjInputs.push_back(p2);
185  fjInputs.push_back(p3);
186  // Run Fastjet algorithm
187  std::vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
188  fastjet::ClusterSequence clustSeq(fjInputs, *jetDef);
189  inclusiveJets = clustSeq.inclusive_jets();
190  sortedJets = sorted_by_pt(inclusiveJets);
191  if (sortedJets.size()>1) continue;
192  File << int(3) <<std::endl;
193  File << axis <<std::endl;
194  File << p1 << std::endl;
195  File << p2 << std::endl;
196  File << p3 << std::endl;
197  ++nJets;
198  }}}}
199  std::cout << "Found "<<nJets << std::endl;
200 
201  break;
202  default: std::cout << "This does not correspond to a valid mode"<<std::endl;
203  break;
204 
205 }
206 
207 }
208 
209 
210 
211 
212 #endif // end MATCHING_H
Generated on Fri Dec 14 2012 11:41:56 for Template Overlap Method by   doxygen 1.8.2