Herwig/Analysis/JetsPlusAnalysis.h

0001 // -*- C++ -*-
0002 #ifndef Herwig_JetsPlusAnalysis_H
0003 #define Herwig_JetsPlusAnalysis_H
0004 //
0005 // This is the declaration of the JetsPlusAnalysis class.
0006 //
0007
0008 #include "ThePEG/Handlers/AnalysisHandler.h"
0009 #include "ThePEG/Cuts/JetFinder.h"
0010 #include "ThePEG/Cuts/JetRegion.h"
0011 #include "ThePEG/Repository/EventGenerator.h"
0012 #include "Herwig/Utilities/Statistics/Histogram.h"
0013
0014 namespace Herwig {
0015
0016 using namespace ThePEG;
0017
0018 /**
0019  * Here is the documentation of the JetsPlusAnalysis class.
0020  *
0021  * @see \ref JetsPlusAnalysisInterfaces "The interfaces"
0022  * defined for JetsPlusAnalysis.
0023  */
0024 class JetsPlusAnalysis: public AnalysisHandler {
0025
0026 public:
0027
0028   /**
0029    * The default constructor.
0030    */
0031   JetsPlusAnalysis();
0032
0033 public:
0034
0035   /** @name Virtual functions required by the AnalysisHandler class. */
0036   //@{
0037   /**
0038    * Analyze a given Event. Note that a fully generated event
0039    * may be presented several times, if it has been manipulated in
0040    * between. The default version of this function will call transform
0041    * to make a lorentz transformation of the whole event, then extract
0042    * all final state particles and call analyze(tPVector) of this
0043    * analysis object and those of all associated analysis objects. The
0044    * default version will not, however, do anything on events which
0045    * have not been fully generated, or have been manipulated in any
0046    * way.
0047    * @param event pointer to the Event to be analyzed.
0048    * @param ieve the event number.
0049    * @param loop the number of times this event has been presented.
0050    * If negative the event is now fully generated.
0051    * @param state a number different from zero if the event has been
0052    * manipulated in some way since it was last presented.
0053    */
0054   virtual void analyze(tEventPtr event, long ieve, int loop, int state);
0055   //@}
0056
0057 protected:
0058
0059   /**
0060    * Analyze one subprocess, given the event number it belongs to
0061    */
0062   void analyze(ParticleVector&, long, double);
0063
0064   /**
0065    * Clear the hard objects and jets for the next event
0066    */
0067   void clear() {
0068     theHardObjects.clear();
0069     theJets.clear();
0070   }
0071
0072   /**
0073    * Reconstruct the desired electroweak objects and fill the
0074    * respective momenta. Remove the reconstructed particles from the
0075    * list.
0076    */
0077   virtual void reconstructHardObjects(ParticleVector&) {}
0078
0079   /**
0080    * Set the momentum of the indicated electroweak object.
0081    */
0082   LorentzMomentum& hardObjectMomentum(const string& id) {
0083     return theHardObjects[id];
0084   }
0085
0086   /**
0087    * Reconstruct the jets and fill the respective momenta.
0088    */
0089   virtual void reconstructJets(const ParticleVector&);
0090
0091   /**
0092    * The jet finder to use
0093    */
0094   Ptr<JetFinder>::tptr jetFinder() const {
0095     return theJetFinder;
0096   }
0097
0098   /**
0099    * The jet regions to match.
0100    */
0101   const vector<Ptr<JetRegion>::ptr>& jetRegions() const { return theJetRegions; }
0102
0103   /**
0104    * Return the number of matched jets
0105    */
0106   unsigned int nJets() const { return theJets.size(); }
0107
0108   /**
0109    * Set the momentum of the indicated jet.
0110    */
0111   LorentzMomentum& jetMomentum(const unsigned int id) {
0112     return theJets[id];
0113   }
0114
0115 protected:
0116
0117   /**
0118    * Finalize this object. Called in the run phase just after a
0119    * run has ended. Used eg. to write out statistics.
0120    */
0121   virtual void dofinish();
0122
0123 public:
0124
0125   /** @name Functions used by the persistent I/O system. */
0126   //@{
0127   /**
0128    * Function used to write out object persistently.
0129    * @param os the persistent output stream written to.
0130    */
0131   void persistentOutput(PersistentOStream & os) const;
0132
0133   /**
0134    * Function used to read in object persistently.
0135    * @param is the persistent input stream read from.
0136    * @param version the version number of the object when written.
0137    */
0138   void persistentInput(PersistentIStream & is, int version);
0139   //@}
0140
0141   /**
0142    * The standard Init function used to initialize the interfaces.
0143    * Called exactly once for each class by the class description system
0144    * before the main function starts or
0145    * when this class is dynamically loaded.
0146    */
0147   static void Init();
0148
0149 protected:
0150
0151   /** @name Clone Methods. */
0152   //@{
0153   /**
0154    * Make a simple clone of this object.
0155    * @return a pointer to the new object.
0156    */
0157   virtual IBPtr clone() const;
0158
0159   /** Make a clone of this object, possibly modifying the cloned object
0160    * to make it sane.
0161    * @return a pointer to the new object.
0162    */
0163   virtual IBPtr fullclone() const;
0164   //@}
0165
0166
0167 // If needed, insert declarations of virtual function defined in the
0168 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0169
0170 protected:
0171
0172   /**
0173    * Collection of object histograms; ranges are adjusted to the
0174    * maximum, so range constraints and rebinning can be applied later.
0175    */
0176   struct ObjectProperties {
0177
0178     /**
0179      * Transverse momentum
0180      */
0181     Statistics::Histogram pt;
0182     Statistics::Histogram pt_logx;
0183
0184     /**
0185      * Rapidity
0186      */
0187     Statistics::Histogram y;
0188
0189     /**
0190      * Azimuth
0191      */
0192     Statistics::Histogram phi;
0193
0194     /**
0195      * Mass
0196      */
0197     Statistics::Histogram mass;
0198
0199     /**
0200      * Default constructor
0201      */
0202     ObjectProperties() {}
0203
0204     /**
0205      * Construct given Ecm
0206      */
0207     ObjectProperties(const string& name, Energy)
0208       : pt(name + "Pt",Statistics::Histogram::regularBinEdges(0,1000,1000),true,false),
0209         pt_logx(name + "PtLogX",Statistics::Histogram::logBinEdges(0.1,1000,1000),true,false),
0210     y(name + "Y",Statistics::Histogram::regularBinEdges(-6,6,120),false,false),
0211     phi(name + "Phi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,32),
0212         make_pair(-Constants::pi,Constants::pi)),
0213     mass(name + "Mass",Statistics::Histogram::regularBinEdges(0,1000,1000),true,false) {}
0214
0215     /**
0216      * Count given momentum, weight and id
0217      */
0218     void count(const LorentzMomentum& p, double weight, unsigned int id) {
0219       pt.count(Statistics::EventContribution(p.perp()/GeV,weight,1.),id);
0220       pt_logx.count(Statistics::EventContribution(p.perp()/GeV,weight,1.),id);
0221       y.count(Statistics::EventContribution(p.rapidity(),weight,0.1),id);
0222       phi.count(Statistics::EventContribution(p.phi(),weight,0.1),id);
0223       mass.count(Statistics::EventContribution(p.m()/GeV,weight,1.),id);
0224     }
0225
0226     /**
0227      * Count given momentum components, weight and id
0228      */
0229     void count(Energy perp, double rapidity,
0230            double xphi, Energy m,
0231            double weight, unsigned int id) {
0232       pt.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
0233       pt_logx.count(Statistics::EventContribution(perp/GeV,weight,1.),id);
0234       y.count(Statistics::EventContribution(rapidity,weight,0.1),id);
0235       phi.count(Statistics::EventContribution(xphi,weight,0.1),id);
0236       mass.count(Statistics::EventContribution(m/GeV,weight,1.),id);
0237     }
0238
0239     /**
0240      * Convert to XML
0241      */
0242     void finalize(XML::Element& elem) {
0243       pt.finalize(); elem.append(pt.toXML());
0244       pt_logx.finalize(); elem.append(pt_logx.toXML());
0245       y.finalize(); elem.append(y.toXML());
0246       phi.finalize(); elem.append(phi.toXML());
0247       mass.finalize(); elem.append(mass.toXML());
0248     }
0249
0250   };
0251
0252   /**
0253    * Collection of pair histograms; ranges are adjusted to the
0254    * maximum, so range constraints and rebinning can be applied later.
0255    */
0256   struct PairProperties
0257     : public ObjectProperties {
0258
0259     /**
0260      * Calculate deltaPhi
0261      */
0262     static double dPhi(const LorentzMomentum& a,
0263                const LorentzMomentum& b){
0264       double phi1 = a.phi();
0265       double phi2 = b.phi();
0266       double diff=phi1-phi2;
0267       if(diff<-Constants::pi){
0268     diff+=(2.0*Constants::pi);
0269       }
0270       else if (diff>Constants::pi){
0271     diff-=(2.0*Constants::pi);
0272       }
0273       return diff;
0274     }
0275
0276     /**
0277      * Calculate deltaY
0278      */
0279     static double dY(const LorentzMomentum& a,
0280              const LorentzMomentum& b){
0281       return abs(a.rapidity()-b.rapidity());
0282     }
0283
0284     /**
0285      * Calculate deltaR
0286      */
0287     static double dR(const LorentzMomentum& a,
0288              const LorentzMomentum& b){
0289       return sqrt(sqr(dPhi(a,b))+sqr(dY(a,b)));
0290     }
0291
0292     /**
0293      * Calculate ydoty
0294      */
0295     static double yy(const LorentzMomentum& a,
0296              const LorentzMomentum& b){
0297       double ya = a.rapidity();
0298       double yb = b.rapidity();
0299       double yres = sqrt(abs(ya*yb));
0300       return ya*yb < 0. ? -yres : yres;
0301     }
0302
0303     /**
0304      * Delta y
0305      */
0306     Statistics::Histogram deltaY;
0307
0308     /**
0309      * Delta phi
0310      */
0311     Statistics::Histogram deltaPhi;
0312
0313     /**
0314      * Delta phi
0315      */
0316     Statistics::Histogram deltaR;
0317
0318     /**
0319      * Product of the rapidities
0320      */
0321     Statistics::Histogram yDotY;
0322
0323     /**
0324      * Default constructor
0325      */
0326     PairProperties()
0327       : ObjectProperties() {}
0328
0329     /**
0330      * Construct given Ecm
0331      */
0332     PairProperties(const string& name, Energy ecm)
0333       : ObjectProperties(name,ecm),
0334     deltaY(name + "DeltaY",Statistics::Histogram::regularBinEdges(0,6,60),true,false),
0335     deltaPhi(name + "DeltaPhi",Statistics::Histogram::regularBinEdges(-Constants::pi,Constants::pi,32),
0336          make_pair(-Constants::pi,Constants::pi)),
0337     deltaR(name + "DeltaR",Statistics::Histogram::regularBinEdges(0,10,100),true,false),
0338     yDotY(name + "YDotY",Statistics::Histogram::regularBinEdges(-6,6,120),false,false) {}
0339
0340     /**
0341      * Count given momentum, weight and id
0342      */
0343     void count(const LorentzMomentum& p, const LorentzMomentum& q, double weight, unsigned int id) {
0344       ObjectProperties::count(p+q,weight,id);
0345       deltaY.count(Statistics::EventContribution(dY(p,q),weight,0.1),id);
0346       deltaPhi.count(Statistics::EventContribution(dPhi(p,q),weight,0.1),id);
0347       deltaR.count(Statistics::EventContribution(dR(p,q),weight,0.1),id);
0348       yDotY.count(Statistics::EventContribution(yy(p,q),weight,0.1),id);
0349     }
0350
0351     /**
0352      * Convert to XML
0353      */
0354     void finalize(XML::Element& elem) {
0355       ObjectProperties::finalize(elem);
0356       deltaY.finalize(); elem.append(deltaY.toXML());
0357       deltaPhi.finalize(); elem.append(deltaPhi.toXML());
0358       deltaR.finalize(); elem.append(deltaR.toXML());
0359       yDotY.finalize(); elem.append(yDotY.toXML());
0360     }
0361
0362   };
0363
0364 private:
0365
0366   /**
0367    * Switch between fixed order and showered
0368    */
0369   bool theIsShowered;
0370
0371   /**
0372    * The jet finder to use
0373    */
0374   Ptr<JetFinder>::ptr theJetFinder;
0375
0376   /**
0377    * The jet regions to match.
0378    */
0379   vector<Ptr<JetRegion>::ptr> theJetRegions;
0380
0381   /**
0382    * The reconstructed hard objects.
0383    */
0384   map<string,LorentzMomentum> theHardObjects;
0385
0386   /**
0387    * The reconstructed jets
0388    */
0389   map<unsigned int,LorentzMomentum> theJets;
0390
0391   /**
0392    * Hard object properties
0393    */
0394   map<string,ObjectProperties> theHardObjectProperties;
0395
0396   /**
0397    * Jet properties
0398    */
0399   map<unsigned int,ObjectProperties> theJetProperties;
0400
0401   /**
0402    * Exclusive jet properties
0403    */
0404   map<unsigned int,ObjectProperties> theExclusiveJetProperties;
0405
0406   /**
0407    * Jet-inclusive properties
0408    */
0409   ObjectProperties theJetInclusiveProperties;
0410
0411   /**
0412    * Jet-summed properties
0413    */
0414   ObjectProperties theJetSummedProperties;
0415
0416   /**
0417    * Jet-average properties
0418    */
0419   ObjectProperties theJetAverageProperties;
0420
0421   /**
0422    * Inclusive jet multiplicities
0423    */
0424   Statistics::Histogram theNJetsInclusive;
0425
0426   /**
0427    * Exclusive jet multiplicities
0428    */
0429   Statistics::Histogram theNJetsExclusive;
0430
0431   /**
0432    * Hard object pair properties
0433    */
0434   map<pair<string,string>,PairProperties> theHardPairProperties;
0435
0436   /**
0437    * Jet pair properties
0438    */
0439   map<pair<unsigned int,unsigned int>,PairProperties> theJetPairProperties;
0440
0441   /**
0442    * Jet/hard pair properties
0443    */
0444   map<pair<unsigned int,string>,PairProperties> theJetHardPairProperties;
0445
0446   /**
0447    * Trijet properties
0448    */
0449   map<std::tuple<unsigned int,unsigned int,unsigned int>,ObjectProperties> theThreeJetProperties;
0450
0451   /**
0452    * Fourjet properties
0453    */
0454   map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties> theFourJetProperties;
0455
0456 protected:
0457
0458   /**
0459    * Hard object properties
0460    */
0461   ObjectProperties& hardObjectProperties(const string& id) {
0462     map<string,ObjectProperties>::iterator h =
0463       theHardObjectProperties.find(id);
0464     if ( h != theHardObjectProperties.end() )
0465       return h->second;
0466     return
0467       theHardObjectProperties[id] =
0468       ObjectProperties(id,generator()->maximumCMEnergy());
0469   }
0470
0471   /**
0472    * Jet properties
0473    */
0474   ObjectProperties& jetProperties(const unsigned int id) {
0475     map<unsigned int,ObjectProperties>::iterator h =
0476       theJetProperties.find(id);
0477     if ( h != theJetProperties.end() )
0478       return h->second;
0479     ostringstream ids; ids << "Jet" << id;
0480     return
0481       theJetProperties[id] =
0482       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0483   }
0484
0485   /**
0486    * Exclusive jet properties
0487    */
0488   ObjectProperties& exclusiveJetProperties(const unsigned int id) {
0489     map<unsigned int,ObjectProperties>::iterator h =
0490       theExclusiveJetProperties.find(id);
0491     if ( h != theExclusiveJetProperties.end() )
0492       return h->second;
0493     ostringstream ids; ids << "ExclusiveJet" << id;
0494     return
0495       theExclusiveJetProperties[id] =
0496       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0497   }
0498
0499   /**
0500    * Jet-inclusive properties
0501    */
0502   ObjectProperties& jetInclusiveProperties() {
0503     if ( !theJetInclusiveProperties.pt.bins().empty() )
0504       return theJetInclusiveProperties;
0505     return
0506       theJetInclusiveProperties =
0507       ObjectProperties("JetInclusive",generator()->maximumCMEnergy());
0508   }
0509
0510   /**
0511    * Jet-summed properties
0512    */
0513   ObjectProperties& jetSummedProperties() {
0514     if ( !theJetSummedProperties.pt.bins().empty() )
0515       return theJetSummedProperties;
0516     return
0517       theJetSummedProperties =
0518       ObjectProperties("JetSummed",generator()->maximumCMEnergy());
0519   }
0520
0521   /**
0522    * Jet-average properties
0523    */
0524   ObjectProperties& jetAverageProperties() {
0525     if ( !theJetAverageProperties.pt.bins().empty() )
0526       return theJetAverageProperties;
0527     return
0528       theJetAverageProperties =
0529       ObjectProperties("JetAverage",generator()->maximumCMEnergy());
0530   }
0531
0532
0533   /**
0534    * Inclusive jet multiplicities
0535    */
0536   Statistics::Histogram& nJetsInclusive() {
0537     if ( !theNJetsInclusive.bins().empty() )
0538       return theNJetsInclusive;
0539     return
0540       theNJetsInclusive =
0541       Statistics::Histogram("NJetsInclusive",
0542                 Statistics::Histogram::regularBinEdges(-0.5,theJetRegions.size()+0.5,
0543                                    theJetRegions.size()+1),
0544                 true,true);
0545   }
0546
0547   /**
0548    * Exclusive jet multiplicities
0549    */
0550   Statistics::Histogram& nJetsExclusive() {
0551     if ( !theNJetsExclusive.bins().empty() )
0552       return theNJetsExclusive;
0553     return
0554       theNJetsExclusive =
0555       Statistics::Histogram("NJetsExclusive",
0556                 Statistics::Histogram::regularBinEdges(-0.5,theJetRegions.size()+0.5,
0557                                    theJetRegions.size()+1),
0558                 true,true);
0559   }
0560
0561   /**
0562    * Hard object pair properties
0563    */
0564   PairProperties& hardPairProperties(const string& id, const string& jd) {
0565     map<pair<string,string>,PairProperties>::iterator h =
0566       theHardPairProperties.find(make_pair(id,jd));
0567     if ( h != theHardPairProperties.end() )
0568       return h->second;
0569     return theHardPairProperties[make_pair(id,jd)] =
0570       PairProperties(id+jd,generator()->maximumCMEnergy());
0571   }
0572
0573   /**
0574    * Jet pair properties
0575    */
0576   PairProperties& jetPairProperties(const unsigned int id, const unsigned int jd) {
0577     map<pair<unsigned int,unsigned int>,PairProperties>::iterator h =
0578       theJetPairProperties.find(make_pair(id,jd));
0579     if ( h != theJetPairProperties.end() )
0580       return h->second;
0581     ostringstream ids; ids << "Jet" << id << jd;
0582     return theJetPairProperties[make_pair(id,jd)] =
0583       PairProperties(ids.str(),generator()->maximumCMEnergy());
0584   }
0585
0586   /**
0587    * Jet/hard pair properties
0588    */
0589   PairProperties& jetHardPairProperties(const unsigned int id, const string& jd) {
0590     map<pair<unsigned int,string>,PairProperties>::iterator h =
0591       theJetHardPairProperties.find(make_pair(id,jd));
0592     if ( h != theJetHardPairProperties.end() )
0593       return h->second;
0594     ostringstream ids; ids << "Jet" << id << jd;
0595     return theJetHardPairProperties[make_pair(id,jd)] =
0596       PairProperties(ids.str(),generator()->maximumCMEnergy());
0597   }
0598
0599   /**
0600    * Trijet properties
0601    */
0602   ObjectProperties& threeJetProperties(const unsigned int id1, const unsigned int id2,
0603                        const unsigned int id3) {
0604     map<std::tuple<unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
0605       theThreeJetProperties.find(std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3));
0606     if ( it != theThreeJetProperties.end() )
0607       return it->second;
0608     ostringstream ids;
0609     ids << "Jet" << id1 << id2 << id3;
0610     return theThreeJetProperties[std::tuple<unsigned int,unsigned int,unsigned int>(id1,id2,id3)] =
0611       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0612   }
0613
0614   /**
0615    * Fourjet properties
0616    */
0617   ObjectProperties& fourJetProperties(const unsigned int id1, const unsigned int id2,
0618                       const unsigned int id3, const unsigned int id4) {
0619     map<std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>,ObjectProperties>::iterator it =
0620       theFourJetProperties.find(std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4));
0621     if ( it != theFourJetProperties.end() )
0622       return it->second;
0623     ostringstream ids;
0624     ids << "Jet" << id1 << id2 << id3 << id4;
0625     return theFourJetProperties[std::tuple<unsigned int,unsigned int,unsigned int,unsigned int>(id1,id2,id3,id4)] =
0626       ObjectProperties(ids.str(),generator()->maximumCMEnergy());
0627   }
0628
0629   /**
0630    * Perform any additional analysis required
0631    */
0632   virtual void analyzeSpecial(long, double) {}
0633
0634   /**
0635    * Append any additional histograms to the given histogram element
0636    */
0637   virtual void finalize(XML::Element&) {}
0638
0639 private:
0640
0641   /**
0642    * The assignment operator is private and must never be called.
0643    * In fact, it should not even be implemented.
0644    */
0645   JetsPlusAnalysis & operator=(const JetsPlusAnalysis &) = delete;
0646
0647 };
0648
0649 }
0650
0651 #endif /* Herwig_JetsPlusAnalysis_H */