00001 // -*- C++ -*-
00002 
00003 // Field & Stuart underlying event analysis at CDF.
00004 // Phys.Rev.D65:092002,2002 - no arXiv code.
00005 // FNAL-PUB 01/211-E
00006 
00007 #include "Rivet/Analysis.hh"
00008 #include "Rivet/RivetAIDA.hh"
00009 #include "Rivet/Tools/Logging.hh"
00010 #include "Rivet/Projections/ChargedFinalState.hh"
00011 #include "Rivet/Projections/ConstLossyFinalState.hh"
00012 #include "Rivet/Projections/FastJets.hh"
00013 #include "Rivet/Projections/TriggerCDFRun0Run1.hh"
00014 #include "LWH/Profile1D.h"
00015 
00016 namespace Rivet {
00017 
00018 
00019   /* @brief "Field-Stuart" CDF Run I track-jet underlying event analysis
00020    * @author Andy Buckley
00021    *
00022    * The "original" underlying event analysis, using a non-standard track-jet algorithm.
00023    *
00024    * @par Run conditions
00025    *
00026    * @arg \f$ \sqrt{s} = \f$ 1800 GeV
00027    * @arg Run with generic QCD events.
00028    * @arg Several \f$ p_\perp^\text{min} \f$ cutoffs are probably required to fill the profile histograms:
00029    *   @arg \f$ p_\perp^\text{min} = \f$ 0 (min bias), 10, 20 GeV
00030    *
00031    */
00032   class CDF_2001_S4751469 : public Analysis {
00033   public:
00034 
00035     /// Constructor: cuts on final state are \f$ -1 < \eta < 1 \f$
00036     /// and \f$ p_T > 0.5 \f$ GeV.
00037     CDF_2001_S4751469()
00038       : Analysis("CDF_2001_S4751469"),
00039         _totalNumTrans2(0), _totalNumTrans5(0), _totalNumTrans30(0),
00040         _sumWeightsPtLead2(0),_sumWeightsPtLead5(0), _sumWeightsPtLead30(0)
00041     {
00042       setBeams(PROTON, ANTIPROTON);
00043     }
00044 
00045 
00046     /// @name Analysis methods
00047     //@{
00048 
00049     // Book histograms
00050     void init() {
00051       addProjection(TriggerCDFRun0Run1(), "Trigger");
00052       // Randomly discard 8% of charged particles as a kind of hacky detector correction.
00053       const ChargedFinalState cfs(-1.0, 1.0, 0.5*GeV);
00054       const ConstLossyFinalState lfs(cfs, 0.08);
00055 
00056       // const LossyFinalState<StarRandomFilter> lfs(cfs, StarRandomFilter());
00057       addProjection(lfs, "FS");
00058       addProjection(FastJets(lfs, FastJets::TRACKJET, 0.7), "TrackJet");
00059 
00060       _numvsDeltaPhi2 =  bookProfile1D(1, 1, 1);
00061       _numvsDeltaPhi5 =  bookProfile1D(1, 1, 2);
00062       _numvsDeltaPhi30 = bookProfile1D(1, 1, 3);
00063       _pTvsDeltaPhi2 =   bookProfile1D(2, 1, 1);
00064       _pTvsDeltaPhi5 =   bookProfile1D(2, 1, 2);
00065       _pTvsDeltaPhi30 =  bookProfile1D(2, 1, 3);
00066 
00067       _numTowardMB = bookProfile1D(3, 1, 1);
00068       _numTransMB = bookProfile1D(3, 1, 2);
00069       _numAwayMB = bookProfile1D(3, 1, 3);
00070       _numTowardJ20 = bookProfile1D(4, 1, 1);
00071       _numTransJ20 = bookProfile1D(4, 1, 2);
00072       _numAwayJ20 = bookProfile1D(4, 1, 3);
00073 
00074       _ptsumTowardMB = bookProfile1D(5, 1, 1);
00075       _ptsumTransMB = bookProfile1D(5, 1, 2);
00076       _ptsumAwayMB = bookProfile1D(5, 1, 3);
00077       _ptsumTowardJ20 = bookProfile1D(6, 1, 1);
00078       _ptsumTransJ20 = bookProfile1D(6, 1, 2);
00079       _ptsumAwayJ20 = bookProfile1D(6, 1, 3);
00080 
00081       _ptTrans2 = bookHistogram1D(7, 1, 1);
00082       _ptTrans5 = bookHistogram1D(7, 1, 2);
00083       _ptTrans30 = bookHistogram1D(7, 1, 3);
00084     }
00085 
00086 
00087     /// Do the analysis
00088     void analyze(const Event& event) {
00089       // Trigger
00090       const bool trigger = applyProjection<TriggerCDFRun0Run1>(event, "Trigger").minBiasDecision();
00091       if (!trigger) vetoEvent;
00092 
00093       // Analyse, with pT > 0.5 GeV AND |eta| < 1
00094       const JetAlg& tj = applyProjection<JetAlg>(event, "TrackJet");
00095 
00096       // Final state (lossy) charged particles
00097       const FinalState& fs = applyProjection<FinalState>(event, "FS");
00098 
00099       // Get jets, sorted by pT
00100       const Jets jets = tj.jetsByPt();
00101       if (jets.empty()) {
00102         vetoEvent;
00103       }
00104 
00105       Jet leadingJet = jets.front();
00106       const double phiLead = leadingJet.ptWeightedPhi();
00107       const double ptLead = leadingJet.ptSum();
00108 
00109       // Cut on highest pT jet: combined 0.5 GeV < pT(lead) < 50 GeV
00110       if (ptLead/GeV < 0.5) vetoEvent;
00111       if (ptLead/GeV > 50.0) vetoEvent;
00112 
00113       // Get the event weight
00114       const double weight = event.weight();
00115 
00116       // Count sum of all event weights in three pT_lead regions
00117       if (ptLead/GeV > 2.0) {
00118         _sumWeightsPtLead2 += weight;
00119       }
00120       if (ptLead/GeV > 5.0) {
00121         _sumWeightsPtLead5 += weight;
00122       }
00123       if (ptLead/GeV > 30.0) {
00124         _sumWeightsPtLead30 += weight;
00125       }
00126 
00127       // Run over tracks
00128       double ptSumToward(0.0), ptSumAway(0.0), ptSumTrans(0.0);
00129       size_t numToward(0), numTrans(0), numAway(0);
00130 
00131       // Temporary histos that bin N and pT in dphi
00132       /// @todo Copy the permanent histos to get the binnings more robustly
00133       LWH::Profile1D hist_num_dphi_2(50, 0, 180), hist_num_dphi_5(50, 0, 180), hist_num_dphi_30(50, 0, 180);
00134       LWH::Profile1D hist_pt_dphi_2(50, 0, 180), hist_pt_dphi_5(50, 0, 180), hist_pt_dphi_30(50, 0, 180);
00135 
00136       foreach (const Particle& p, fs.particles()) {
00137         // Calculate DeltaPhi(p,leadingJet)
00138         const double dPhi = deltaPhi(p.momentum().phi(), phiLead);
00139         const double pT = p.momentum().pT();
00140 
00141         if (dPhi < PI/3.0) {
00142           ptSumToward += pT;
00143           ++numToward;
00144         }
00145         else if (dPhi < 2*PI/3.0) {
00146           ptSumTrans += pT;
00147           ++numTrans;
00148           // Fill transverse pT distributions
00149           if (ptLead/GeV > 2.0) {
00150             _ptTrans2->fill(pT/GeV, weight);
00151             _totalNumTrans2 += weight;
00152           }
00153           if (ptLead/GeV > 5.0) {
00154             _ptTrans5->fill(pT/GeV, weight);
00155             _totalNumTrans5 += weight;
00156           }
00157           if (ptLead/GeV > 30.0) {
00158             _ptTrans30->fill(pT/GeV, weight);
00159             _totalNumTrans30 += weight;
00160           }
00161         }
00162         else {
00163           ptSumAway += pT;
00164           ++numAway;
00165         }
00166 
00167         // Fill tmp histos to bin event's track Nch & pT in dphi
00168         const double dPhideg = 180*dPhi/PI;
00169         if (ptLead/GeV > 2.0) {
00170           hist_num_dphi_2.fill(dPhideg, 1);
00171           hist_pt_dphi_2.fill (dPhideg, pT/GeV);
00172         }
00173         if (ptLead/GeV > 5.0) {
00174           hist_num_dphi_5.fill(dPhideg, 1);
00175           hist_pt_dphi_5.fill (dPhideg, pT/GeV);
00176         }
00177         if (ptLead/GeV > 30.0) {
00178           hist_num_dphi_30.fill(dPhideg, 1);
00179           hist_pt_dphi_30.fill (dPhideg, pT/GeV);
00180         }
00181       }
00182 
00183       // Update the "proper" dphi profile histograms
00184       for (int i = 0; i < 50; i++) {
00185         if (ptLead/GeV > 2.0) {
00186           _numvsDeltaPhi2->fill(hist_num_dphi_2.binMean(i), hist_num_dphi_2.binHeight(i), weight);
00187           _pTvsDeltaPhi2->fill(hist_pt_dphi_2.binMean(i), hist_pt_dphi_2.binHeight(i), weight);
00188         }
00189         if (ptLead/GeV > 5.0) {
00190           _numvsDeltaPhi5->fill(hist_num_dphi_5.binMean(i), hist_num_dphi_5.binHeight(i), weight);
00191           _pTvsDeltaPhi5->fill(hist_pt_dphi_5.binMean(i), hist_pt_dphi_5.binHeight(i), weight);
00192         }
00193         if (ptLead/GeV > 30.0) {
00194           _numvsDeltaPhi30->fill(hist_num_dphi_30.binMean(i), hist_num_dphi_30.binHeight(i), weight);
00195           _pTvsDeltaPhi30->fill(hist_pt_dphi_30.binMean(i), hist_pt_dphi_30.binHeight(i), weight);
00196         }
00197       }
00198 
00199       // Log some event details about pT
00200       getLog() << Log::DEBUG
00201                << "pT [lead; twd, away, trans] = ["
00202                << ptLead << "; "
00203                << ptSumToward << ", "
00204                << ptSumAway << ", "
00205                << ptSumTrans << "]"
00206                << endl;
00207 
00208       // Update the pT profile histograms
00209       _ptsumTowardMB->fill(ptLead/GeV, ptSumToward/GeV, weight);
00210       _ptsumTowardJ20->fill(ptLead/GeV, ptSumToward/GeV, weight);
00211 
00212       _ptsumTransMB->fill(ptLead/GeV, ptSumTrans/GeV, weight);
00213       _ptsumTransJ20->fill(ptLead/GeV, ptSumTrans/GeV, weight);
00214 
00215       _ptsumAwayMB->fill(ptLead/GeV, ptSumAway/GeV, weight);
00216       _ptsumAwayJ20->fill(ptLead/GeV, ptSumAway/GeV, weight);
00217 
00218       // Log some event details about Nch
00219       getLog() << Log::DEBUG
00220                << "N [twd, away, trans] = ["
00221                << numToward << ", "
00222                << numTrans << ", "
00223                << numAway << "]"
00224                << endl;
00225 
00226       // Update the N_track profile histograms
00227       _numTowardMB->fill(ptLead/GeV, numToward, weight);
00228       _numTowardJ20->fill(ptLead/GeV, numToward, weight);
00229 
00230       _numTransMB->fill(ptLead/GeV, numTrans, weight);
00231       _numTransJ20->fill(ptLead/GeV, numTrans, weight);
00232 
00233       _numAwayMB->fill(ptLead/GeV, numAway, weight);
00234       _numAwayJ20->fill(ptLead/GeV, numAway, weight);
00235     }
00236 
00237 
00238     /// Normalize histos
00239     void finalize() {
00240       normalize(_ptTrans2, _totalNumTrans2 / _sumWeightsPtLead2);
00241       normalize(_ptTrans5, _totalNumTrans5 / _sumWeightsPtLead5);
00242       normalize(_ptTrans30, _totalNumTrans30 / _sumWeightsPtLead30);
00243     }
00244 
00245     //@}
00246 
00247 
00248   private:
00249 
00250     /// Sum total number of charged particles in the trans region, in 3 \f$ p_\perp^\text{lead} \f$ bins.
00251     double _totalNumTrans2, _totalNumTrans5, _totalNumTrans30;
00252 
00253     /// Sum the total number of events in 3 \f$ p_\perp^\text{lead} \f$ bins.
00254     double _sumWeightsPtLead2,_sumWeightsPtLead5, _sumWeightsPtLead30;
00255 
00256 
00257     /// @name Histogram collections
00258     //@{
00259     // These histos (binned in dphi) are filled per event and then reset
00260     // TODO: use LWH
00261     AIDA::IProfile1D *_hist_num_dphi_2, *_hist_num_dphi_5, *_hist_num_dphi_30;
00262     AIDA::IProfile1D *_hist_pt_dphi_2, *_hist_pt_dphi_5, *_hist_pt_dphi_30;
00263 
00264     // The sumpt vs. dphi and Nch vs. dphi histos
00265     AIDA::IProfile1D *_numvsDeltaPhi2, *_numvsDeltaPhi5, *_numvsDeltaPhi30;
00266     AIDA::IProfile1D *_pTvsDeltaPhi2, *_pTvsDeltaPhi5, *_pTvsDeltaPhi30;
00267 
00268 
00269     /// Profile histograms, binned in the \f$ p_T \f$ of the leading jet, for
00270     /// the \f$ p_T \f$ sum in the toward, transverse and away regions.
00271     AIDA::IProfile1D *_ptsumTowardMB,  *_ptsumTransMB,  *_ptsumAwayMB;
00272     AIDA::IProfile1D *_ptsumTowardJ20, *_ptsumTransJ20, *_ptsumAwayJ20;
00273 
00274     /// Profile histograms, binned in the \f$ p_T \f$ of the leading jet, for
00275     /// the number of charged particles per jet in the toward, transverse and
00276     /// away regions.
00277     AIDA::IProfile1D *_numTowardMB,  *_numTransMB,  *_numAwayMB;
00278     AIDA::IProfile1D *_numTowardJ20, *_numTransJ20, *_numAwayJ20;
00279 
00280     /// Histogram of \f$ p_T \f$ distribution for 3 different \f$ p_{T1} \f$ IR cutoffs.
00281     AIDA::IHistogram1D *_ptTrans2, *_ptTrans5, *_ptTrans30;
00282     //@}
00283 
00284   };
00285 
00286 
00287 
00288   // This global object acts as a hook for the plugin system
00289   AnalysisBuilder<CDF_2001_S4751469> plugin_CDF_2001_S4751469;
00290 
00291 }
CDF_2001_S4751469.cc
