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