MC_DIPHOTON.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "Rivet/Analysis.hh"
00003 #include "Rivet/RivetAIDA.hh"
00004 #include "Rivet/Projections/FinalState.hh"
00005 #include "Rivet/Projections/IdentifiedFinalState.hh"
00006 #include "Rivet/Tools/Logging.hh"
00007 
00008 namespace Rivet {
00009 
00010 
00011   class MC_DIPHOTON : public Analysis {
00012   public:
00013 
00014     /// Constructor
00015     MC_DIPHOTON() : Analysis("MC_DIPHOTON") {
00016       setNeedsCrossSection(true);
00017     }
00018 
00019 
00020     /// @name Analysis methods
00021     //@{
00022 
00023     void init() {
00024       FinalState fs;
00025       addProjection(fs, "FS");
00026    
00027       IdentifiedFinalState ifs(-2.0, 2.0, 20.0*GeV);
00028       ifs.acceptId(PHOTON);
00029       addProjection(ifs, "IFS");
00030 
00031       _h_m_PP = bookHistogram1D("m_PP", logBinEdges(50, 1.0, 0.25*sqrtS()));
00032       _h_pT_PP = bookHistogram1D("pT_PP", logBinEdges(50, 1.0, 0.25*sqrtS()));
00033       _h_dphi_PP = bookHistogram1D("dphi_PP", 20, 0.0, M_PI);
00034     }
00035 
00036 
00037     void analyze(const Event& event) {
00038       const double weight = event.weight();
00039    
00040       ParticleVector photons = applyProjection<IdentifiedFinalState>(event, "IFS").particles();
00041    
00042       if (photons.size() < 2) {
00043         vetoEvent;
00044       }
00045    
00046       // Isolate photons with ET_sum in cone
00047       ParticleVector isolated_photons;
00048       ParticleVector fs = applyProjection<FinalState>(event, "FS").particles();
00049       foreach (const Particle& photon, photons) {
00050         FourMomentum mom_in_cone;
00051         double eta_P = photon.momentum().eta();
00052         double phi_P = photon.momentum().phi();
00053         foreach (const Particle& p, fs) {
00054           if (deltaR(eta_P, phi_P, p.momentum().eta(), p.momentum().phi()) < 0.4) {
00055             mom_in_cone += p.momentum();
00056           }
00057         }
00058         if (mom_in_cone.Et()-photon.momentum().Et() < 1.0*GeV) {
00059           isolated_photons.push_back(photon);
00060         }
00061       }
00062    
00063       if (isolated_photons.size() != 2) {
00064         vetoEvent;
00065       }
00066    
00067       FourMomentum mom_PP = isolated_photons[0].momentum() + isolated_photons[1].momentum();
00068       _h_m_PP->fill(mom_PP.mass(), weight);
00069       _h_pT_PP->fill(mom_PP.pT(), weight);
00070       _h_dphi_PP->fill(mapAngle0ToPi(isolated_photons[0].momentum().phi()-
00071                                      isolated_photons[1].momentum().phi())/M_PI, weight);
00072     }
00073  
00074  
00075     void finalize() {
00076       scale(_h_m_PP, crossSection()/sumOfWeights());
00077       scale(_h_pT_PP, crossSection()/sumOfWeights());
00078       scale(_h_dphi_PP, crossSection()/sumOfWeights());
00079     }
00080 
00081     //@}
00082 
00083 
00084   private:
00085 
00086     /// @name Histograms
00087     //@{
00088     AIDA::IHistogram1D* _h_m_PP;
00089     AIDA::IHistogram1D* _h_pT_PP;
00090     AIDA::IHistogram1D* _h_dphi_PP;
00091     //@}
00092  
00093  
00094   };
00095 
00096 
00097 
00098   // This global object acts as a hook for the plugin system
00099   AnalysisBuilder<MC_DIPHOTON> plugin_MC_DIPHOTON;
00100 
00101 }
00102