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