herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
JetsPlusAnalysis.h
1 // -*- C++ -*-
2 #ifndef Herwig_JetsPlusAnalysis_H
3 #define Herwig_JetsPlusAnalysis_H
4 //
5 // This is the declaration of the JetsPlusAnalysis class.
6 //
7 
8 #include "ThePEG/Handlers/AnalysisHandler.h"
9 #include "ThePEG/Cuts/JetFinder.h"
10 #include "ThePEG/Cuts/JetRegion.h"
11 #include "ThePEG/Repository/EventGenerator.h"
12 #include "Herwig/Utilities/Statistics/Histogram.h"
13 
14 namespace Herwig {
15 
16 using namespace ThePEG;
17 
25 
26 public:
27 
34 
38  virtual ~JetsPlusAnalysis();
40 
41 public:
42 
62  virtual void analyze(tEventPtr event, long ieve, int loop, int state);
64 
65 protected:
66 
70  void analyze(ParticleVector&, long, double);
71 
75  void clear() {
76  theHardObjects.clear();
77  theJets.clear();
78  }
79 
86 
90  LorentzMomentum& hardObjectMomentum(const string& id) {
91  return theHardObjects[id];
92  }
93 
97  virtual void reconstructJets(const ParticleVector&);
98 
102  Ptr<JetFinder>::tptr jetFinder() const {
103  return theJetFinder;
104  }
105 
109  const vector<Ptr<JetRegion>::ptr>& jetRegions() const { return theJetRegions; }
110 
114  unsigned int nJets() const { return theJets.size(); }
115 
119  LorentzMomentum& jetMomentum(const unsigned int id) {
120  return theJets[id];
121  }
122 
123 protected:
124 
129  virtual void dofinish();
130 
131 public:
132 
139  void persistentOutput(PersistentOStream & os) const;
140 
146  void persistentInput(PersistentIStream & is, int version);
148 
155  static void Init();
156 
157 protected:
158 
165  virtual IBPtr clone() const;
166 
171  virtual IBPtr fullclone() const;
173 
174 
175 // If needed, insert declarations of virtual function defined in the
176 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
177 
178 protected:
179 
185 
190  Statistics::Histogram pt_logx;
191 
196 
201 
206 
211 
215  ObjectProperties(const string& name, Energy)
216  : pt(name + "Pt",Statistics::Histogram::regularBinEdges(0,1000,1000),true,false),
217  pt_logx(name + "PtLogX",Statistics::Histogram::logBinEdges(0.1,1000,1000),true,false),
218  y(name + "Y",Statistics::Histogram::regularBinEdges(-6,6,120),false,false),
219  phi(name + "Phi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,32),
220  make_pair(-Constants::pi,Constants::pi)),
221  mass(name + "Mass",Statistics::Histogram::regularBinEdges(0,1000,1000),true,false) {}
222 
226  void count(const LorentzMomentum& p, double weight, unsigned int id) {
227  pt.count(Statistics::EventContribution(p.perp()/GeV,weight,1.),id);
228  pt_logx.count(Statistics::EventContribution(p.perp()/GeV,weight,1.),id);
229  y.count(Statistics::EventContribution(p.rapidity(),weight,0.1),id);
230  phi.count(Statistics::EventContribution(p.phi(),weight,0.1),id);
231  mass.count(Statistics::EventContribution(p.m()/GeV,weight,1.),id);
232  }
233 
237  void count(Energy perp, double rapidity,
238  double xphi, Energy m,
239  double weight, unsigned int id) {
240  pt.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
241  pt_logx.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
242  y.count(Statistics::EventContribution(rapidity,weight,0.1),id);
243  phi.count(Statistics::EventContribution(xphi,weight,0.1),id);
244  mass.count(Statistics::EventContribution(m/GeV,weight,1.),id);
245  }
246 
250  void finalize(XML::Element& elem) {
251  pt.finalize(); elem.append(pt.toXML());
252  pt_logx.finalize(); elem.append(pt_logx.toXML());
253  y.finalize(); elem.append(y.toXML());
254  phi.finalize(); elem.append(phi.toXML());
255  mass.finalize(); elem.append(mass.toXML());
256  }
257 
258  };
259 
265  : public ObjectProperties {
266 
270  static double dPhi(const LorentzMomentum& a,
271  const LorentzMomentum& b){
272  double phi1 = a.phi();
273  double phi2 = b.phi();
274  double diff=phi1-phi2;
275  if(diff<-Constants::pi){
276  diff+=(2.0*Constants::pi);
277  }
278  else if (diff>Constants::pi){
279  diff-=(2.0*Constants::pi);
280  }
281  return diff;
282  }
283 
287  static double dY(const LorentzMomentum& a,
288  const LorentzMomentum& b){
289  return abs(a.rapidity()-b.rapidity());
290  }
291 
295  static double dR(const LorentzMomentum& a,
296  const LorentzMomentum& b){
297  return sqrt(sqr(dPhi(a,b))+sqr(dY(a,b)));
298  }
299 
303  static double yy(const LorentzMomentum& a,
304  const LorentzMomentum& b){
305  double ya = a.rapidity();
306  double yb = b.rapidity();
307  double yres = sqrt(abs(ya*yb));
308  return ya*yb < 0. ? -yres : yres;
309  }
310 
315 
320 
325 
330 
335  : ObjectProperties() {}
336 
340  PairProperties(const string& name, Energy ecm)
341  : ObjectProperties(name,ecm),
342  deltaY(name + "DeltaY",Statistics::Histogram::regularBinEdges(0,6,60),true,false),
343  deltaPhi(name + "DeltaPhi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,32),
344  make_pair(-Constants::pi,Constants::pi)),
345  deltaR(name + "DeltaR",Statistics::Histogram::regularBinEdges(0,10,100),true,false),
346  yDotY(name + "YDotY",Statistics::Histogram::regularBinEdges(-6,6,120),false,false) {}
347 
351  void count(const LorentzMomentum& p, const LorentzMomentum& q, double weight, unsigned int id) {
352  ObjectProperties::count(p+q,weight,id);
353  deltaY.count(Statistics::EventContribution(dY(p,q),weight,0.1),id);
354  deltaPhi.count(Statistics::EventContribution(dPhi(p,q),weight,0.1),id);
355  deltaR.count(Statistics::EventContribution(dR(p,q),weight,0.1),id);
356  yDotY.count(Statistics::EventContribution(yy(p,q),weight,0.1),id);
357  }
358 
362  void finalize(XML::Element& elem) {
363  ObjectProperties::finalize(elem);
364  deltaY.finalize(); elem.append(deltaY.toXML());
365  deltaPhi.finalize(); elem.append(deltaPhi.toXML());
366  deltaR.finalize(); elem.append(deltaR.toXML());
367  yDotY.finalize(); elem.append(yDotY.toXML());
368  }
369 
370  };
371 
372 private:
373 
378 
382  Ptr<JetFinder>::ptr theJetFinder;
383 
387  vector<Ptr<JetRegion>::ptr> theJetRegions;
388 
392  map<string,LorentzMomentum> theHardObjects;
393 
397  map<unsigned int,LorentzMomentum> theJets;
398 
402  map<string,ObjectProperties> theHardObjectProperties;
403 
407  map<unsigned int,ObjectProperties> theJetProperties;
408 
412  map<unsigned int,ObjectProperties> theExclusiveJetProperties;
413 
418 
423 
428 
433 
438 
442  map<pair<string,string>,PairProperties> theHardPairProperties;
443 
447  map<pair<unsigned int,unsigned int>,PairProperties> theJetPairProperties;
448 
452  map<pair<unsigned int,string>,PairProperties> theJetHardPairProperties;
453 
457  map<std::tuple<unsigned int,unsigned int,unsigned int>,ObjectProperties> theThreeJetProperties;
458 
462  map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourJetProperties;
463 
464 protected:
465 
470  map<string,ObjectProperties>::iterator h =
471  theHardObjectProperties.find(id);
472  if ( h != theHardObjectProperties.end() )
473  return h->second;
474  return
475  theHardObjectProperties[id] =
476  ObjectProperties(id,generator()->maximumCMEnergy());
477  }
478 
482  ObjectProperties& jetProperties(const unsigned int id) {
483  map<unsigned int,ObjectProperties>::iterator h =
484  theJetProperties.find(id);
485  if ( h != theJetProperties.end() )
486  return h->second;
487  ostringstream ids; ids << "Jet" << id;
488  return
489  theJetProperties[id] =
490  ObjectProperties(ids.str(),generator()->maximumCMEnergy());
491  }
492 
496  ObjectProperties& exclusiveJetProperties(const unsigned int id) {
497  map<unsigned int,ObjectProperties>::iterator h =
498  theExclusiveJetProperties.find(id);
499  if ( h != theExclusiveJetProperties.end() )
500  return h->second;
501  ostringstream ids; ids << "ExclusiveJet" << id;
502  return
503  theExclusiveJetProperties[id] =
504  ObjectProperties(ids.str(),generator()->maximumCMEnergy());
505  }
506 
511  if ( !theJetInclusiveProperties.pt.bins().empty() )
512  return theJetInclusiveProperties;
513  return
514  theJetInclusiveProperties =
515  ObjectProperties("JetInclusive",generator()->maximumCMEnergy());
516  }
517 
522  if ( !theJetSummedProperties.pt.bins().empty() )
523  return theJetSummedProperties;
524  return
525  theJetSummedProperties =
526  ObjectProperties("JetSummed",generator()->maximumCMEnergy());
527  }
528 
533  if ( !theJetAverageProperties.pt.bins().empty() )
534  return theJetAverageProperties;
535  return
536  theJetAverageProperties =
537  ObjectProperties("JetAverage",generator()->maximumCMEnergy());
538  }
539 
540 
545  if ( !theNJetsInclusive.bins().empty() )
546  return theNJetsInclusive;
547  return
548  theNJetsInclusive =
549  Statistics::Histogram("NJetsInclusive",
550  Statistics::Histogram::regularBinEdges(-0.5,theJetRegions.size()+0.5,
551  theJetRegions.size()+1),
552  true,true);
553  }
554 
559  if ( !theNJetsExclusive.bins().empty() )
560  return theNJetsExclusive;
561  return
562  theNJetsExclusive =
563  Statistics::Histogram("NJetsExclusive",
564  Statistics::Histogram::regularBinEdges(-0.5,theJetRegions.size()+0.5,
565  theJetRegions.size()+1),
566  true,true);
567  }
568 
572  PairProperties& hardPairProperties(const string& id, const string& jd) {
573  map<pair<string,string>,PairProperties>::iterator h =
574  theHardPairProperties.find(make_pair(id,jd));
575  if ( h != theHardPairProperties.end() )
576  return h->second;
577  return theHardPairProperties[make_pair(id,jd)] =
578  PairProperties(id+jd,generator()->maximumCMEnergy());
579  }
580 
584  PairProperties& jetPairProperties(const unsigned int id, const unsigned int jd) {
585  map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
586  theJetPairProperties.find(make_pair(id,jd));
587  if ( h != theJetPairProperties.end() )
588  return h->second;
589  ostringstream ids; ids << "Jet" << id << jd;
590  return theJetPairProperties[make_pair(id,jd)] =
591  PairProperties(ids.str(),generator()->maximumCMEnergy());
592  }
593 
597  PairProperties& jetHardPairProperties(const unsigned int id, const string& jd) {
598  map<pair<unsigned int,string>,PairProperties>::iterator h =
599  theJetHardPairProperties.find(make_pair(id,jd));
600  if ( h != theJetHardPairProperties.end() )
601  return h->second;
602  ostringstream ids; ids << "Jet" << id << jd;
603  return theJetHardPairProperties[make_pair(id,jd)] =
604  PairProperties(ids.str(),generator()->maximumCMEnergy());
605  }
606 
610  ObjectProperties& threeJetProperties(const unsigned int id1, const unsigned int id2,
611  const unsigned int id3) {
612  map<std::tuple<unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
613  theThreeJetProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
614  if ( it != theThreeJetProperties.end() )
615  return it->second;
616  ostringstream ids;
617  ids << "Jet" << id1 << id2 << id3;
618  return theThreeJetProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
619  ObjectProperties(ids.str(),generator()->maximumCMEnergy());
620  }
621 
625  ObjectProperties& fourJetProperties(const unsigned int id1, const unsigned int id2,
626  const unsigned int id3, const unsigned int id4) {
627  map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
628  theFourJetProperties.find(std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
629  if ( it != theFourJetProperties.end() )
630  return it->second;
631  ostringstream ids;
632  ids << "Jet" << id1 << id2 << id3 << id4;
633  return theFourJetProperties[std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
634  ObjectProperties(ids.str(),generator()->maximumCMEnergy());
635  }
636 
640  virtual void analyzeSpecial(long, double) {}
641 
645  virtual void finalize(XML::Element&) {}
646 
647 private:
648 
653  JetsPlusAnalysis & operator=(const JetsPlusAnalysis &) = delete;
654 
655 };
656 
657 }
658 
659 #endif /* Herwig_JetsPlusAnalysis_H */
A (one dimensional) histogram.
double sqrt(int x)
void clear()
Clear the hard objects and jets for the next event.
map< pair< string, string >, PairProperties > theHardPairProperties
Hard object pair properties.
Collection of object histograms; ranges are adjusted to the maximum, so range constraints and rebinni...
Collection of pair histograms; ranges are adjusted to the maximum, so range constraints and rebinning...
Element & append(const Element &)
Append a child element to this element.
void finalize()
Finalize this histogram.
static double dY(const LorentzMomentum &a, const LorentzMomentum &b)
Calculate deltaY.
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
constexpr double pi
The Histogram class is a simple histogram for the Analysis handlers.
Definition: Histogram.h:43
constexpr auto sqr(const T &x) -> decltype(x *x)
Statistics::Histogram yDotY
Product of the rapidities.
static double dPhi(const LorentzMomentum &a, const LorentzMomentum &b)
Calculate deltaPhi.
ObjectProperties & jetAverageProperties()
Jet-average properties.
Statistics::Histogram deltaPhi
Delta phi.
Definition: Bin.h:16
Statistics::Histogram y
Rapidity.
Ptr< JetFinder >::tptr jetFinder() const
The jet finder to use.
ObjectProperties & threeJetProperties(const unsigned int id1, const unsigned int id2, const unsigned int id3)
Trijet properties.
ObjectProperties(const string &name, Energy)
Construct given Ecm.
Statistics::Histogram theNJetsInclusive
Inclusive jet multiplicities.
PairProperties & hardPairProperties(const string &id, const string &jd)
Hard object pair properties.
ObjectProperties theJetInclusiveProperties
Jet-inclusive properties.
map< std::tuple< unsigned int, unsigned int, unsigned int, unsigned int >, ObjectProperties > theFourJetProperties
Fourjet properties.
ObjectProperties & hardObjectProperties(const string &id)
Hard object properties.
PairProperties(const string &name, Energy ecm)
Construct given Ecm.
ObjectProperties & jetProperties(const unsigned int id)
Jet properties.
map< pair< unsigned int, string >, PairProperties > theJetHardPairProperties
Jet/hard pair properties.
void finalize(XML::Element &elem)
Convert to XML.
Statistics::Histogram deltaR
Delta phi.
map< string, LorentzMomentum > theHardObjects
The reconstructed hard objects.
Statistics::Histogram deltaY
Delta y.
Statistics::Histogram theNJetsExclusive
Exclusive jet multiplicities.
void count(const LorentzMomentum &p, const LorentzMomentum &q, double weight, unsigned int id)
Count given momentum, weight and id.
map< unsigned int, ObjectProperties > theExclusiveJetProperties
Exclusive jet properties.
ObjectProperties & jetInclusiveProperties()
Jet-inclusive properties.
Statistics::Histogram & nJetsInclusive()
Inclusive jet multiplicities.
Ptr< JetFinder >::ptr theJetFinder
The jet finder to use.
const std::vector< Bin > & bins() const
Return the bins.
bool theIsShowered
Switch between fixed order and showered.
ObjectProperties & jetSummedProperties()
Jet-summed properties.
PairProperties & jetPairProperties(const unsigned int id, const unsigned int jd)
Jet pair properties.
virtual void reconstructHardObjects(ParticleVector &)
Reconstruct the desired electroweak objects and fill the respective momenta.
XML::Element toXML() const
Return an XML element for the data of this histogram.
ObjectProperties theJetSummedProperties
Jet-summed properties.
void count(Energy perp, double rapidity, double xphi, Energy m, double weight, unsigned int id)
Count given momentum components, weight and id.
static std::vector< double > regularBinEdges(double lower, double upper, size_t nBins)
Create equally spaced bin edges.
static double yy(const LorentzMomentum &a, const LorentzMomentum &b)
Calculate ydoty.
A pointlike or boxlike eventContribution; serves to define the EventContribution concept.
Statistics::Histogram phi
Azimuth.
Element represents a (tree of) XML elements.
Definition: Element.h:56
PairProperties & jetHardPairProperties(const unsigned int id, const string &jd)
Jet/hard pair properties.
map< unsigned int, ObjectProperties > theJetProperties
Jet properties.
map< string, ObjectProperties > theHardObjectProperties
Hard object properties.
ObjectProperties & exclusiveJetProperties(const unsigned int id)
Exclusive jet properties.
ThePEG::Ptr< Event >::transient_pointer tEventPtr
LorentzMomentum & hardObjectMomentum(const string &id)
Set the momentum of the indicated electroweak object.
bool count(EventContribution event, size_t id)
Book a contribution to the current event.
Statistics::Histogram & nJetsExclusive()
Exclusive jet multiplicities.
map< std::tuple< unsigned int, unsigned int, unsigned int >, ObjectProperties > theThreeJetProperties
Trijet properties.
virtual void analyzeSpecial(long, double)
Perform any additional analysis required.
void finalize(XML::Element &elem)
Convert to XML.
-*- C++ -*-
virtual void finalize(XML::Element &)
Append any additional histograms to the given histogram element.
ObjectProperties & fourJetProperties(const unsigned int id1, const unsigned int id2, const unsigned int id3, const unsigned int id4)
Fourjet properties.
vector< Ptr< JetRegion >::ptr > theJetRegions
The jet regions to match.
ObjectProperties theJetAverageProperties
Jet-average properties.
unsigned int nJets() const
Return the number of matched jets.
void count(const LorentzMomentum &p, double weight, unsigned int id)
Count given momentum, weight and id.
vector< PPtr > ParticleVector
static double dR(const LorentzMomentum &a, const LorentzMomentum &b)
Calculate deltaR.
Here is the documentation of the JetsPlusAnalysis class.
Statistics::Histogram pt
Transverse momentum.
LorentzMomentum & jetMomentum(const unsigned int id)
Set the momentum of the indicated jet.
map< pair< unsigned int, unsigned int >, PairProperties > theJetPairProperties
Jet pair properties.
map< unsigned int, LorentzMomentum > theJets
The reconstructed jets.
const vector< Ptr< JetRegion >::ptr > & jetRegions() const
The jet regions to match.