herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
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
14namespace Herwig {
15
16using namespace ThePEG;
17
25
26public:
27
32
33public:
34
54 virtual void analyze(tEventPtr event, long ieve, int loop, int state);
56
57protected:
58
62 void analyze(ParticleVector&, long, double);
63
67 void clear() {
68 theHardObjects.clear();
69 theJets.clear();
70 }
71
78
82 LorentzMomentum& hardObjectMomentum(const string& id) {
83 return theHardObjects[id];
84 }
85
89 virtual void reconstructJets(const ParticleVector&);
90
94 Ptr<JetFinder>::tptr jetFinder() const {
95 return theJetFinder;
96 }
97
101 const vector<Ptr<JetRegion>::ptr>& jetRegions() const { return theJetRegions; }
102
106 unsigned int nJets() const { return theJets.size(); }
107
111 LorentzMomentum& jetMomentum(const unsigned int id) {
112 return theJets[id];
113 }
114
115protected:
116
121 virtual void dofinish();
122
123public:
124
132
138 void persistentInput(PersistentIStream & is, int version);
140
147 static void Init();
148
149protected:
150
157 virtual IBPtr clone() const;
158
163 virtual IBPtr fullclone() const;
165
166
167// If needed, insert declarations of virtual function defined in the
168// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
169
170protected:
171
177
182 Statistics::Histogram pt_logx;
183
188
193
198
203
207 ObjectProperties(const string& name, Energy)
208 : pt(name + "Pt",Statistics::Histogram::regularBinEdges(0,1000,1000),true,false),
209 pt_logx(name + "PtLogX",Statistics::Histogram::logBinEdges(0.1,1000,1000),true,false),
210 y(name + "Y",Statistics::Histogram::regularBinEdges(-6,6,120),false,false),
211 phi(name + "Phi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,32),
212 make_pair(-Constants::pi,Constants::pi)),
213 mass(name + "Mass",Statistics::Histogram::regularBinEdges(0,1000,1000),true,false) {}
214
218 void count(const LorentzMomentum& p, double weight, unsigned int id) {
219 pt.count(Statistics::EventContribution(p.perp()/GeV,weight,1.),id);
220 pt_logx.count(Statistics::EventContribution(p.perp()/GeV,weight,1.),id);
221 y.count(Statistics::EventContribution(p.rapidity(),weight,0.1),id);
222 phi.count(Statistics::EventContribution(p.phi(),weight,0.1),id);
223 mass.count(Statistics::EventContribution(p.m()/GeV,weight,1.),id);
224 }
225
229 void count(Energy perp, double rapidity,
230 double xphi, Energy m,
231 double weight, unsigned int id) {
232 pt.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
233 pt_logx.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
234 y.count(Statistics::EventContribution(rapidity,weight,0.1),id);
235 phi.count(Statistics::EventContribution(xphi,weight,0.1),id);
236 mass.count(Statistics::EventContribution(m/GeV,weight,1.),id);
237 }
238
242 void finalize(XML::Element& elem) {
243 pt.finalize(); elem.append(pt.toXML());
244 pt_logx.finalize(); elem.append(pt_logx.toXML());
245 y.finalize(); elem.append(y.toXML());
246 phi.finalize(); elem.append(phi.toXML());
247 mass.finalize(); elem.append(mass.toXML());
248 }
249
250 };
251
257 : public ObjectProperties {
258
262 static double dPhi(const LorentzMomentum& a,
263 const LorentzMomentum& b){
264 double phi1 = a.phi();
265 double phi2 = b.phi();
266 double diff=phi1-phi2;
267 if(diff<-Constants::pi){
268 diff+=(2.0*Constants::pi);
269 }
270 else if (diff>Constants::pi){
271 diff-=(2.0*Constants::pi);
272 }
273 return diff;
274 }
275
279 static double dY(const LorentzMomentum& a,
280 const LorentzMomentum& b){
281 return abs(a.rapidity()-b.rapidity());
282 }
283
287 static double dR(const LorentzMomentum& a,
288 const LorentzMomentum& b){
289 return sqrt(sqr(dPhi(a,b))+sqr(dY(a,b)));
290 }
291
295 static double yy(const LorentzMomentum& a,
296 const LorentzMomentum& b){
297 double ya = a.rapidity();
298 double yb = b.rapidity();
299 double yres = sqrt(abs(ya*yb));
300 return ya*yb < 0. ? -yres : yres;
301 }
302
307
312
317
322
327 : ObjectProperties() {}
328
332 PairProperties(const string& name, Energy ecm)
333 : ObjectProperties(name,ecm),
334 deltaY(name + "DeltaY",Statistics::Histogram::regularBinEdges(0,6,60),true,false),
335 deltaPhi(name + "DeltaPhi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,32),
336 make_pair(-Constants::pi,Constants::pi)),
337 deltaR(name + "DeltaR",Statistics::Histogram::regularBinEdges(0,10,100),true,false),
338 yDotY(name + "YDotY",Statistics::Histogram::regularBinEdges(-6,6,120),false,false) {}
339
343 void count(const LorentzMomentum& p, const LorentzMomentum& q, double weight, unsigned int id) {
344 ObjectProperties::count(p+q,weight,id);
345 deltaY.count(Statistics::EventContribution(dY(p,q),weight,0.1),id);
347 deltaR.count(Statistics::EventContribution(dR(p,q),weight,0.1),id);
348 yDotY.count(Statistics::EventContribution(yy(p,q),weight,0.1),id);
349 }
350
354 void finalize(XML::Element& elem) {
356 deltaY.finalize(); elem.append(deltaY.toXML());
358 deltaR.finalize(); elem.append(deltaR.toXML());
359 yDotY.finalize(); elem.append(yDotY.toXML());
360 }
361
362 };
363
364private:
365
370
374 Ptr<JetFinder>::ptr theJetFinder;
375
379 vector<Ptr<JetRegion>::ptr> theJetRegions;
380
384 map<string,LorentzMomentum> theHardObjects;
385
389 map<unsigned int,LorentzMomentum> theJets;
390
394 map<string,ObjectProperties> theHardObjectProperties;
395
399 map<unsigned int,ObjectProperties> theJetProperties;
400
404 map<unsigned int,ObjectProperties> theExclusiveJetProperties;
405
410
415
420
425
430
434 map<pair<string,string>,PairProperties> theHardPairProperties;
435
439 map<pair<unsigned int,unsigned int>,PairProperties> theJetPairProperties;
440
444 map<pair<unsigned int,string>,PairProperties> theJetHardPairProperties;
445
449 map<std::tuple<unsigned int,unsigned int,unsigned int>,ObjectProperties> theThreeJetProperties;
450
454 map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourJetProperties;
455
456protected:
457
462 map<string,ObjectProperties>::iterator h =
464 if ( h != theHardObjectProperties.end() )
465 return h->second;
466 return
468 ObjectProperties(id,generator()->maximumCMEnergy());
469 }
470
474 ObjectProperties& jetProperties(const unsigned int id) {
475 map<unsigned int,ObjectProperties>::iterator h =
476 theJetProperties.find(id);
477 if ( h != theJetProperties.end() )
478 return h->second;
479 ostringstream ids; ids << "Jet" << id;
480 return
481 theJetProperties[id] =
482 ObjectProperties(ids.str(),generator()->maximumCMEnergy());
483 }
484
488 ObjectProperties& exclusiveJetProperties(const unsigned int id) {
489 map<unsigned int,ObjectProperties>::iterator h =
491 if ( h != theExclusiveJetProperties.end() )
492 return h->second;
493 ostringstream ids; ids << "ExclusiveJet" << id;
494 return
496 ObjectProperties(ids.str(),generator()->maximumCMEnergy());
497 }
498
503 if ( !theJetInclusiveProperties.pt.bins().empty() )
505 return
507 ObjectProperties("JetInclusive",generator()->maximumCMEnergy());
508 }
509
514 if ( !theJetSummedProperties.pt.bins().empty() )
516 return
518 ObjectProperties("JetSummed",generator()->maximumCMEnergy());
519 }
520
525 if ( !theJetAverageProperties.pt.bins().empty() )
527 return
529 ObjectProperties("JetAverage",generator()->maximumCMEnergy());
530 }
531
532
537 if ( !theNJetsInclusive.bins().empty() )
538 return theNJetsInclusive;
539 return
541 Statistics::Histogram("NJetsInclusive",
543 theJetRegions.size()+1),
544 true,true);
545 }
546
551 if ( !theNJetsExclusive.bins().empty() )
552 return theNJetsExclusive;
553 return
555 Statistics::Histogram("NJetsExclusive",
557 theJetRegions.size()+1),
558 true,true);
559 }
560
564 PairProperties& hardPairProperties(const string& id, const string& jd) {
565 map<pair<string,string>,PairProperties>::iterator h =
566 theHardPairProperties.find(make_pair(id,jd));
567 if ( h != theHardPairProperties.end() )
568 return h->second;
569 return theHardPairProperties[make_pair(id,jd)] =
570 PairProperties(id+jd,generator()->maximumCMEnergy());
571 }
572
576 PairProperties& jetPairProperties(const unsigned int id, const unsigned int jd) {
577 map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
578 theJetPairProperties.find(make_pair(id,jd));
579 if ( h != theJetPairProperties.end() )
580 return h->second;
581 ostringstream ids; ids << "Jet" << id << jd;
582 return theJetPairProperties[make_pair(id,jd)] =
583 PairProperties(ids.str(),generator()->maximumCMEnergy());
584 }
585
589 PairProperties& jetHardPairProperties(const unsigned int id, const string& jd) {
590 map<pair<unsigned int,string>,PairProperties>::iterator h =
591 theJetHardPairProperties.find(make_pair(id,jd));
592 if ( h != theJetHardPairProperties.end() )
593 return h->second;
594 ostringstream ids; ids << "Jet" << id << jd;
595 return theJetHardPairProperties[make_pair(id,jd)] =
596 PairProperties(ids.str(),generator()->maximumCMEnergy());
597 }
598
602 ObjectProperties& threeJetProperties(const unsigned int id1, const unsigned int id2,
603 const unsigned int id3) {
604 map<std::tuple<unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
605 theThreeJetProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
606 if ( it != theThreeJetProperties.end() )
607 return it->second;
608 ostringstream ids;
609 ids << "Jet" << id1 << id2 << id3;
610 return theThreeJetProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
611 ObjectProperties(ids.str(),generator()->maximumCMEnergy());
612 }
613
617 ObjectProperties& fourJetProperties(const unsigned int id1, const unsigned int id2,
618 const unsigned int id3, const unsigned int id4) {
619 map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
620 theFourJetProperties.find(std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
621 if ( it != theFourJetProperties.end() )
622 return it->second;
623 ostringstream ids;
624 ids << "Jet" << id1 << id2 << id3 << id4;
625 return theFourJetProperties[std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
626 ObjectProperties(ids.str(),generator()->maximumCMEnergy());
627 }
628
632 virtual void analyzeSpecial(long, double) {}
633
637 virtual void finalize(XML::Element&) {}
638
639private:
640
646
647};
648
649}
650
651#endif /* Herwig_JetsPlusAnalysis_H */
The Histogram class is a simple histogram for the Analysis handlers.
Definition: Histogram.h:43
Here is the documentation of the JetsPlusAnalysis class.
const vector< Ptr< JetRegion >::ptr > & jetRegions() const
The jet regions to match.
virtual void dofinish()
Finalize this object.
virtual IBPtr clone() const
Make a simple clone of this object.
ObjectProperties & jetInclusiveProperties()
Jet-inclusive properties.
static void Init()
The standard Init function used to initialize the interfaces.
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.
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
map< pair< unsigned int, unsigned int >, PairProperties > theJetPairProperties
Jet pair properties.
PairProperties & jetHardPairProperties(const unsigned int id, const string &jd)
Jet/hard pair properties.
ObjectProperties & hardObjectProperties(const string &id)
Hard object properties.
ObjectProperties & threeJetProperties(const unsigned int id1, const unsigned int id2, const unsigned int id3)
Trijet properties.
LorentzMomentum & hardObjectMomentum(const string &id)
Set the momentum of the indicated electroweak object.
virtual void finalize(XML::Element &)
Append any additional histograms to the given histogram element.
virtual void analyze(tEventPtr event, long ieve, int loop, int state)
Analyze a given Event.
Ptr< JetFinder >::tptr jetFinder() const
The jet finder to use.
JetsPlusAnalysis()
The default constructor.
bool theIsShowered
Switch between fixed order and showered.
ObjectProperties & exclusiveJetProperties(const unsigned int id)
Exclusive jet properties.
virtual void reconstructJets(const ParticleVector &)
Reconstruct the jets and fill the respective momenta.
ObjectProperties & jetProperties(const unsigned int id)
Jet properties.
unsigned int nJets() const
Return the number of matched jets.
void clear()
Clear the hard objects and jets for the next event.
map< std::tuple< unsigned int, unsigned int, unsigned int, unsigned int >, ObjectProperties > theFourJetProperties
Fourjet properties.
map< std::tuple< unsigned int, unsigned int, unsigned int >, ObjectProperties > theThreeJetProperties
Trijet properties.
Statistics::Histogram & nJetsInclusive()
Inclusive jet multiplicities.
ObjectProperties theJetInclusiveProperties
Jet-inclusive properties.
Statistics::Histogram & nJetsExclusive()
Exclusive jet multiplicities.
map< string, LorentzMomentum > theHardObjects
The reconstructed hard objects.
virtual void analyzeSpecial(long, double)
Perform any additional analysis required.
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
Statistics::Histogram theNJetsExclusive
Exclusive jet multiplicities.
ObjectProperties & fourJetProperties(const unsigned int id1, const unsigned int id2, const unsigned int id3, const unsigned int id4)
Fourjet properties.
PairProperties & hardPairProperties(const string &id, const string &jd)
Hard object pair properties.
map< unsigned int, ObjectProperties > theJetProperties
Jet properties.
Ptr< JetFinder >::ptr theJetFinder
The jet finder to use.
LorentzMomentum & jetMomentum(const unsigned int id)
Set the momentum of the indicated jet.
ObjectProperties & jetSummedProperties()
Jet-summed properties.
Statistics::Histogram theNJetsInclusive
Inclusive jet multiplicities.
map< string, ObjectProperties > theHardObjectProperties
Hard object properties.
map< unsigned int, LorentzMomentum > theJets
The reconstructed jets.
map< pair< unsigned int, string >, PairProperties > theJetHardPairProperties
Jet/hard pair properties.
ObjectProperties & jetAverageProperties()
Jet-average properties.
vector< Ptr< JetRegion >::ptr > theJetRegions
The jet regions to match.
void analyze(ParticleVector &, long, double)
Analyze one subprocess, given the event number it belongs to.
ObjectProperties theJetSummedProperties
Jet-summed properties.
map< unsigned int, ObjectProperties > theExclusiveJetProperties
Exclusive jet properties.
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
map< pair< string, string >, PairProperties > theHardPairProperties
Hard object pair properties.
JetsPlusAnalysis & operator=(const JetsPlusAnalysis &)=delete
The assignment operator is private and must never be called.
ObjectProperties theJetAverageProperties
Jet-average properties.
A pointlike or boxlike eventContribution; serves to define the EventContribution concept.
A (one dimensional) histogram.
void finalize()
Finalize this histogram.
bool count(EventContribution event, size_t id)
Book a contribution to the current event.
const std::vector< Bin > & bins() const
Return the bins.
XML::Element toXML() const
Return an XML element for the data of this histogram.
static std::vector< double > regularBinEdges(double lower, double upper, size_t nBins)
Create equally spaced bin edges.
InitState state() const
string name() const
tEGPtr generator() const
Element represents a (tree of) XML elements.
Definition: Element.h:56
Element & append(const Element &)
Append a child element to this element.
-*- C++ -*-
constexpr double pi
double sqrt(int x)
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
ThePEG::Ptr< Event >::transient_pointer tEventPtr
vector< PPtr > ParticleVector
constexpr auto sqr(const T &x) -> decltype(x *x)
Collection of object histograms; ranges are adjusted to the maximum, so range constraints and rebinni...
ObjectProperties(const string &name, Energy)
Construct given Ecm.
void count(Energy perp, double rapidity, double xphi, Energy m, double weight, unsigned int id)
Count given momentum components, weight and id.
void finalize(XML::Element &elem)
Convert to XML.
Statistics::Histogram pt
Transverse momentum.
void count(const LorentzMomentum &p, double weight, unsigned int id)
Count given momentum, weight and id.
Collection of pair histograms; ranges are adjusted to the maximum, so range constraints and rebinning...
static double dR(const LorentzMomentum &a, const LorentzMomentum &b)
Calculate deltaR.
static double dY(const LorentzMomentum &a, const LorentzMomentum &b)
Calculate deltaY.
Statistics::Histogram deltaY
Delta y.
void count(const LorentzMomentum &p, const LorentzMomentum &q, double weight, unsigned int id)
Count given momentum, weight and id.
static double yy(const LorentzMomentum &a, const LorentzMomentum &b)
Calculate ydoty.
PairProperties(const string &name, Energy ecm)
Construct given Ecm.
static double dPhi(const LorentzMomentum &a, const LorentzMomentum &b)
Calculate deltaPhi.
Statistics::Histogram yDotY
Product of the rapidities.
Statistics::Histogram deltaPhi
Delta phi.
Statistics::Histogram deltaR
Delta phi.
void finalize(XML::Element &elem)
Convert to XML.