MC_PHOTONJETUE.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "Rivet/Analysis.hh"
00003 #include "Rivet/Tools/Logging.hh"
00004 #include "Rivet/RivetAIDA.hh"
00005 #include "Rivet/Projections/ChargedFinalState.hh"
00006 #include "Rivet/Projections/IdentifiedFinalState.hh"
00007 #include "Rivet/Projections/FastJets.hh"
00008 
00009 namespace Rivet {
00010 
00011 
00012   /* Underlying event in jet + isolated photon events
00013    * @author Andy Buckley
00014    */
00015   class MC_PHOTONJETUE : public Analysis {
00016   public:
00017  
00018     /// Constructor
00019     MC_PHOTONJETUE()
00020       : Analysis("MC_PHOTONJETUE")
00021     {    }
00022  
00023  
00024     /// @name Analysis methods
00025     //@{
00026  
00027     // Book histograms and projections
00028     void init() {
00029       // Final state for the jet finding
00030       const FinalState fsj(-4.0, 4.0, 0.1*GeV);
00031       addProjection(fsj, "FSJ");
00032       addProjection(FastJets(fsj, FastJets::ANTIKT, 0.7), "Jets");
00033    
00034       // Charged final state for the distributions
00035       const ChargedFinalState cfs(-2.0, 2.0, 0.2*GeV);
00036       addProjection(cfs, "Tracks");
00037 
00038       // Photons (for isolation)
00039       const FinalState fsp(-2.0, 2.0, 20.0*GeV);
00040       IdentifiedFinalState photonfs(fsp);
00041       photonfs.acceptId(PHOTON);
00042       addProjection(photonfs, "Photons");
00043 
00044       // Histograms
00045       _hist_jetgamma_dR   = bookHistogram1D("gammajet-dR", 52, 0.0, 5.2);
00046       _hist_jetgamma_dphi = bookHistogram1D("gammajet-dphi", 50, 0.0, PI);
00047       //
00048       const double MAXPT1 = 50.0;
00049       _hist_pnchg_jet      = bookProfile1D("trans-nchg-jet",     50, 0.0, MAXPT1);
00050       _hist_pmaxnchg_jet   = bookProfile1D("trans-maxnchg-jet",  50, 0.0, MAXPT1);
00051       _hist_pminnchg_jet   = bookProfile1D("trans-minnchg-jet",  50, 0.0, MAXPT1);
00052       _hist_pcptsum_jet    = bookProfile1D("trans-ptsum-jet",    50, 0.0, MAXPT1);
00053       _hist_pmaxcptsum_jet = bookProfile1D("trans-maxptsum-jet", 50, 0.0, MAXPT1);
00054       _hist_pmincptsum_jet = bookProfile1D("trans-minptsum-jet", 50, 0.0, MAXPT1);
00055       _hist_pcptave_jet    = bookProfile1D("trans-ptavg-jet",    50, 0.0, MAXPT1);
00056       //
00057       _hist_pnchg_gamma      = bookProfile1D("trans-nchg-gamma",     50, 0.0, MAXPT1);
00058       _hist_pmaxnchg_gamma   = bookProfile1D("trans-maxnchg-gamma",  50, 0.0, MAXPT1);
00059       _hist_pminnchg_gamma   = bookProfile1D("trans-minnchg-gamma",  50, 0.0, MAXPT1);
00060       _hist_pcptsum_gamma    = bookProfile1D("trans-ptsum-gamma",    50, 0.0, MAXPT1);
00061       _hist_pmaxcptsum_gamma = bookProfile1D("trans-maxptsum-gamma", 50, 0.0, MAXPT1);
00062       _hist_pmincptsum_gamma = bookProfile1D("trans-minptsum-gamma", 50, 0.0, MAXPT1);
00063       _hist_pcptave_gamma    = bookProfile1D("trans-ptavg-gamma",    50, 0.0, MAXPT1);
00064     }
00065 
00066 
00067     // Do the analysis
00068     void analyze(const Event& evt) {
00069 
00070       // Get jets
00071       const Jets jets = applyProjection<FastJets>(evt, "Jets").jetsByPt();
00072       getLog() << Log::DEBUG << "Jet multiplicity = " << jets.size() << endl;
00073       if (jets.size() < 1) {
00074         getLog() << Log::DEBUG << "No jets found" << endl;
00075         vetoEvent;
00076       }
00077 
00078       // Get leading jet properties
00079       const FourMomentum pjet = jets.front().momentum();
00080       const double jeteta = pjet.eta();
00081       const double jetphi = pjet.phi();
00082       const double jetpT  = pjet.pT();
00083       getLog() << Log::DEBUG << "Leading jet: pT = " << jetpT/GeV << " GeV"
00084                << ", eta = " << jeteta
00085                << ", phi = " << jetphi << endl;
00086 
00087       // Require the leading jet to be within |eta| < 2
00088       if (fabs(jeteta) > 2) {
00089         getLog() << Log::DEBUG << "Failed leading jet eta cut" << endl;
00090         vetoEvent;
00091       }
00092 
00093       // Get the leading photon
00094       const FinalState& photonfs = applyProjection<FinalState>(evt, "Photons");
00095       if (photonfs.size() < 1) {
00096         getLog() << Log::DEBUG << "No hard photons found" << endl;
00097         vetoEvent;
00098       }
00099       const FourMomentum pgamma = photonfs.particlesByPt().front().momentum();
00100 
00101       // Check that leading photon is isolated from jets
00102       bool isolated = true;
00103       foreach (const Jet& j, jets) {
00104         if (deltaR(j.momentum(), pgamma) < 0.1) {
00105           isolated = false;
00106           break;
00107         }
00108       }
00109       if (!isolated) {
00110         getLog() << Log::DEBUG << "Leading photon is not isolated from jets" << endl;
00111         vetoEvent;
00112       }
00113 
00114       // Get leading photon properties
00115       const double gammaeta = pgamma.eta();
00116       const double gammaphi = pgamma.phi();
00117       const double gammapT  = pgamma.pT();
00118       getLog() << Log::DEBUG << "Leading photon: pT = " << gammapT/GeV << " GeV"
00119                << ", eta = " << gammaeta
00120                << ", phi = " << gammaphi << endl;
00121 
00122       // Get the event weight
00123       const double weight = evt.weight();
00124 
00125       // Fill jet1-photon separation histos
00126       const double dR_jetgamma = deltaR(pgamma, pjet);
00127       _hist_jetgamma_dR->fill(dR_jetgamma, weight);
00128       const double dPhi_jetgamma = deltaPhi(gammaphi, jetphi);
00129       _hist_jetgamma_dphi->fill(dPhi_jetgamma, weight);
00130 
00131       /// Cut on back-to-backness of jet-photon
00132       if (dPhi_jetgamma < 0.5) vetoEvent;
00133 
00134       /// @todo Plot evolution of UE as a function of jet-photon angle
00135       /// @todo Plot evolution of UE as a function of photon pT
00136 
00137       // Get the final states to work with for filling the distributions
00138       const FinalState& cfs = applyProjection<ChargedFinalState>(evt, "Tracks");
00139 
00140       // Whole-event counters
00141       unsigned int   numOverall(0);
00142       double ptSumOverall(0.0), ptMaxOverall(0.0);
00143       // Jet-oriented counters
00144       unsigned int   numToward_jet(0),     numTrans1_jet(0),     numTrans2_jet(0),     numAway_jet(0);
00145       double ptSumToward_jet(0.0), ptSumTrans1_jet(0.0), ptSumTrans2_jet(0.0), ptSumAway_jet(0.0);
00146       double ptMaxToward_jet(0.0), ptMaxTrans1_jet(0.0), ptMaxTrans2_jet(0.0), ptMaxAway_jet(0.0);
00147       // Photon-oriented counters
00148       unsigned int   numToward_gamma(0),     numTrans1_gamma(0),     numTrans2_gamma(0),     numAway_gamma(0);
00149       double ptSumToward_gamma(0.0), ptSumTrans1_gamma(0.0), ptSumTrans2_gamma(0.0), ptSumAway_gamma(0.0);
00150       double ptMaxToward_gamma(0.0), ptMaxTrans1_gamma(0.0), ptMaxTrans2_gamma(0.0), ptMaxAway_gamma(0.0);
00151 
00152       // Calculate all the charged stuff
00153       foreach (const Particle& p, cfs.particles()) {
00154         ++numOverall;
00155         const double pT = p.momentum().pT();
00156         ptSumOverall += pT;
00157         if (pT > ptMaxOverall) ptMaxOverall = pT;
00158 
00159         // Increment jet-oriented variables
00160         const double dPhi_jet = jetphi - p.momentum().phi();
00161         if (fabs(dPhi_jet) < PI/3.0) {
00162           ptSumToward_jet += pT;
00163           ++numToward_jet;
00164           if (pT > ptMaxToward_jet) ptMaxToward_jet = pT;
00165         }
00166         else if (fabs(dPhi_jet) < 2*PI/3.0) {
00167           if (sign(dPhi_jet) == MINUS) {
00168             ptSumTrans1_jet += pT;
00169             ++numTrans1_jet;
00170             if (pT > ptMaxTrans1_jet) ptMaxTrans1_jet = pT;
00171           } else {
00172             ptSumTrans2_jet += pT;
00173             ++numTrans2_jet;
00174             if (pT > ptMaxTrans2_jet) ptMaxTrans2_jet = pT;
00175           }
00176         }
00177         else {
00178           ptSumAway_jet += pT;
00179           ++numAway_jet;
00180           if (pT > ptMaxAway_jet) ptMaxAway_jet = pT;
00181         }
00182 
00183 
00184         // Increment photon-oriented variables
00185         const double dPhi_gamma = gammaphi - p.momentum().phi();
00186         if (fabs(dPhi_gamma) < PI/3.0) {
00187           ptSumToward_gamma += pT;
00188           ++numToward_gamma;
00189           if (pT > ptMaxToward_gamma) ptMaxToward_gamma = pT;
00190         }
00191         else if (fabs(dPhi_gamma) < 2*PI/3.0) {
00192           if (sign(dPhi_gamma) == MINUS) {
00193             ptSumTrans1_gamma += pT;
00194             ++numTrans1_gamma;
00195             if (pT > ptMaxTrans1_gamma) ptMaxTrans1_gamma = pT;
00196           } else {
00197             ptSumTrans2_gamma += pT;
00198             ++numTrans2_gamma;
00199             if (pT > ptMaxTrans2_gamma) ptMaxTrans2_gamma = pT;
00200           }
00201         }
00202         else {
00203           ptSumAway_gamma += pT;
00204           ++numAway_gamma;
00205           if (pT > ptMaxAway_gamma) ptMaxAway_gamma = pT;
00206         }
00207 
00208 
00209       }
00210    
00211    
00212       // Fill the histograms
00213       _hist_pnchg_jet->fill(jetpT/GeV, (numTrans1_jet+numTrans2_jet)/(4*PI/3), weight);
00214       _hist_pmaxnchg_jet->fill(jetpT/GeV, (numTrans1_jet>numTrans2_jet ? numTrans1_jet : numTrans2_jet)/(2*PI/3), weight);
00215       _hist_pminnchg_jet->fill(jetpT/GeV, (numTrans1_jet<numTrans2_jet ? numTrans1_jet : numTrans2_jet)/(2*PI/3), weight);
00216       _hist_pcptsum_jet->fill(jetpT/GeV, (ptSumTrans1_jet+ptSumTrans2_jet)/GeV/(4*PI/3), weight);
00217       _hist_pmaxcptsum_jet->fill(jetpT/GeV, (ptSumTrans1_jet>ptSumTrans2_jet ? ptSumTrans1_jet : ptSumTrans2_jet)/GeV/(2*PI/3), weight);
00218       _hist_pmincptsum_jet->fill(jetpT/GeV, (ptSumTrans1_jet<ptSumTrans2_jet ? ptSumTrans1_jet : ptSumTrans2_jet)/GeV/(2*PI/3), weight);
00219       if ((numTrans1_jet+numTrans2_jet) > 0) {
00220         _hist_pcptave_jet->fill(jetpT/GeV, (ptSumTrans1_jet+ptSumTrans2_jet)/GeV/(numTrans1_jet+numTrans2_jet), weight);
00221       }
00222       //
00223       _hist_pnchg_gamma->fill(gammapT/GeV, (numTrans1_gamma+numTrans2_gamma)/(4*PI/3), weight);
00224       _hist_pmaxnchg_gamma->fill(gammapT/GeV, (numTrans1_gamma>numTrans2_gamma ? numTrans1_gamma : numTrans2_gamma)/(2*PI/3), weight);
00225       _hist_pminnchg_gamma->fill(gammapT/GeV, (numTrans1_gamma<numTrans2_gamma ? numTrans1_gamma : numTrans2_gamma)/(2*PI/3), weight);
00226       _hist_pcptsum_gamma->fill(gammapT/GeV, (ptSumTrans1_gamma+ptSumTrans2_gamma)/GeV/(4*PI/3), weight);
00227       _hist_pmaxcptsum_gamma->fill(gammapT/GeV, (ptSumTrans1_gamma>ptSumTrans2_gamma ? ptSumTrans1_gamma : ptSumTrans2_gamma)/GeV/(2*PI/3), weight);
00228       _hist_pmincptsum_gamma->fill(gammapT/GeV, (ptSumTrans1_gamma<ptSumTrans2_gamma ? ptSumTrans1_gamma : ptSumTrans2_gamma)/GeV/(2*PI/3), weight);
00229       if ((numTrans1_gamma+numTrans2_gamma) > 0) {
00230         _hist_pcptave_gamma->fill(gammapT/GeV, (ptSumTrans1_gamma+ptSumTrans2_gamma)/GeV/(numTrans1_gamma+numTrans2_gamma), weight);
00231       }
00232 
00233     }
00234  
00235  
00236     void finalize() {
00237       //
00238     }
00239  
00240  
00241   private:
00242 
00243     AIDA::IHistogram1D* _hist_jetgamma_dR;
00244     AIDA::IHistogram1D* _hist_jetgamma_dphi;
00245  
00246     AIDA::IProfile1D *_hist_pnchg_jet, *_hist_pnchg_gamma;
00247     AIDA::IProfile1D *_hist_pmaxnchg_jet, *_hist_pmaxnchg_gamma;
00248     AIDA::IProfile1D *_hist_pminnchg_jet, *_hist_pminnchg_gamma;
00249     AIDA::IProfile1D *_hist_pcptsum_jet, *_hist_pcptsum_gamma;
00250     AIDA::IProfile1D *_hist_pmaxcptsum_jet, *_hist_pmaxcptsum_gamma;
00251     AIDA::IProfile1D *_hist_pmincptsum_jet, *_hist_pmincptsum_gamma;
00252     AIDA::IProfile1D *_hist_pcptave_jet, *_hist_pcptave_gamma;
00253 
00254   };
00255 
00256 
00257 
00258   // This global object acts as a hook for the plugin system
00259   AnalysisBuilder<MC_PHOTONJETUE> plugin_MC_PHOTONJETUE;
00260 
00261 }