H1_1994_S2919893.cc
Go to the documentation of this file.
00001 // -*- C++ -*- 00002 #include "Rivet/Analysis.hh" 00003 #include "Rivet/Math/Constants.hh" 00004 #include "Rivet/Projections/FinalState.hh" 00005 #include "Rivet/Projections/DISKinematics.hh" 00006 00007 namespace Rivet { 00008 00009 00010 /// @brief H1 energy flow and charged particle spectra 00011 /// @author Peter Richardson 00012 /// Based on the equivalent HZTool analysis 00013 class H1_1994_S2919893 : public Analysis { 00014 public: 00015 00016 /// Constructor 00017 H1_1994_S2919893() 00018 : Analysis("H1_1994_S2919893") 00019 { 00020 00021 // Initialise member variables 00022 _w77 = make_pair(0.0, 0.0); 00023 _w122 = make_pair(0.0, 0.0); 00024 _w169 = make_pair(0.0, 0.0); 00025 _w117 = make_pair(0.0, 0.0); 00026 _wEnergy = make_pair(0.0, 0.0); 00027 } 00028 00029 00030 00031 /// @name Analysis methods 00032 //@{ 00033 00034 void analyze(const Event& event) { 00035 const FinalState& fs = applyProjection<FinalState>(event, "FS"); 00036 const DISKinematics& dk = applyProjection<DISKinematics>(event, "Kinematics"); 00037 const DISLepton& dl = applyProjection<DISLepton>(event,"Lepton"); 00038 00039 // Get the DIS kinematics 00040 double x = dk.x(); 00041 double w2 = dk.W2(); 00042 double w = sqrt(w2); 00043 00044 // Momentum of the scattered lepton 00045 FourMomentum leptonMom = dl.out().momentum(); 00046 double ptel = pT(leptonMom); 00047 double enel = leptonMom.E(); 00048 double thel = leptonMom.angle(dk.beamHadron().momentum())/degree; 00049 00050 // Extract the particles other than the lepton 00051 Particles particles; 00052 particles.reserve(fs.particles().size()); 00053 const GenParticle* dislepGP = dl.out().genParticle(); 00054 foreach (const Particle& p, fs.particles()) { 00055 const GenParticle* loopGP = p.genParticle(); 00056 if (loopGP == dislepGP) continue; 00057 particles.push_back(p); 00058 } 00059 00060 // Cut on the forward energy 00061 double efwd = 0.0; 00062 foreach (const Particle& p, particles) { 00063 double th = p.momentum().angle(dk.beamHadron().momentum())/degree; 00064 if (th > 4.4 && th < 15.) { 00065 efwd += p.momentum().E(); 00066 } 00067 } 00068 00069 // Apply the cuts 00070 // Lepton energy and angle, w2 and forward energy 00071 MSG_DEBUG("enel/GeV = " << enel/GeV << ", thel = " << thel 00072 << ", w2 = " << w2 << ", efwd/GeV = " << efwd/GeV); 00073 bool cut = enel/GeV > 14. && thel > 157. && thel < 172.5 && w2 >= 3000. && efwd/GeV > 0.5; 00074 if (!cut) vetoEvent; 00075 00076 // Weight of the event 00077 const double weight = event.weight(); 00078 // weights for x<1e-3 and x>1e-3 00079 if (x < 1e-3) { 00080 _wEnergy.first += weight; 00081 } else { 00082 _wEnergy.second += weight; 00083 } 00084 00085 // Boost to hadronic CM 00086 const LorentzTransform hcmboost = dk.boostHCM(); 00087 // Loop over the particles 00088 long ncharged(0); 00089 for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) { 00090 const Particle& p = particles[ip1]; 00091 00092 double th = p.momentum().angle(dk.beamHadron().momentum()) / degree; 00093 // Boost momentum to lab 00094 const FourMomentum hcmMom = hcmboost.transform(p.momentum()); 00095 // Angular cut 00096 if (th <= 4.4) continue; 00097 00098 // Energy flow histogram 00099 double et = fabs(Et(hcmMom)); 00100 double eta = hcmMom.pseudorapidity(); 00101 if (x < 1e-3) { 00102 _histEnergyFlowLowX ->fill(eta, et*weight); 00103 } else { 00104 _histEnergyFlowHighX->fill(eta, et*weight); 00105 } 00106 if (PID::threeCharge(p.pdgId()) != 0) { 00107 /// @todo Use units in w comparisons... what are the units? 00108 if (w > 50. && w <= 200.) { 00109 double xf= 2 * hcmMom.z() / w; 00110 double pt2 = pT2(hcmMom); 00111 if (w > 50. && w <= 100.) { 00112 _histSpectraW77 ->fill(xf, weight); 00113 } else if (w > 100. && w <= 150.) { 00114 _histSpectraW122->fill(xf, weight); 00115 } else if (w > 150. && w <= 200.) { 00116 _histSpectraW169->fill(xf, weight); 00117 } 00118 _histSpectraW117->fill(xf, weight); 00119 /// @todo Is this profile meant to be filled with 2 weight factors? 00120 _histPT2->fill(xf, pt2*weight/GeV2, weight); 00121 ++ncharged; 00122 } 00123 } 00124 00125 00126 // Energy-energy correlation 00127 if (th <= 8.) continue; 00128 double phi1 = p.momentum().azimuthalAngle(ZERO_2PI); 00129 double eta1 = p.momentum().pseudorapidity(); 00130 double et1 = fabs(Et(p.momentum())); 00131 for (size_t ip2 = ip1+1; ip2 < particles.size(); ++ip2) { 00132 const Particle& p2 = particles[ip2]; 00133 00134 //double th2 = beamAngle(p2.momentum(), order); 00135 double th2 = p2.momentum().angle(dk.beamHadron().momentum()) / degree; 00136 if (th2 <= 8.) continue; 00137 double phi2 = p2.momentum().azimuthalAngle(ZERO_2PI); 00138 00139 /// @todo Use angle function 00140 double deltaphi = phi1 - phi2; 00141 if (fabs(deltaphi) > PI) 00142 deltaphi = fabs(fabs(deltaphi) - TWOPI); 00143 double eta2 = p2.momentum().pseudorapidity(); 00144 double omega = sqrt(sqr(eta1-eta2) + sqr(deltaphi)); 00145 double et2 = fabs(Et(p2.momentum())); 00146 double wt = et1*et2 / sqr(ptel) * weight; 00147 if(x < 1e-3) { 00148 _histEECLowX ->fill(omega, wt); 00149 } else { 00150 _histEECHighX->fill(omega,wt); 00151 } 00152 } 00153 } 00154 00155 // Factors for normalization 00156 if (w > 50. && w <= 200.) { 00157 if (w <= 100.) { 00158 _w77.first += ncharged*weight; 00159 _w77.second += weight; 00160 } else if (w <= 150.) { 00161 _w122.first += ncharged*weight; 00162 _w122.second += weight; 00163 } else { 00164 _w169.first += ncharged*weight; 00165 _w169.second += weight; 00166 } 00167 _w117.first += ncharged*weight; 00168 _w117.second += weight; 00169 } 00170 } 00171 00172 00173 00174 void init() { 00175 // Projections 00176 addProjection(DISLepton(), "Lepton"); 00177 addProjection(DISKinematics(), "Kinematics"); 00178 addProjection(FinalState(), "FS"); 00179 00180 // Histos 00181 _histEnergyFlowLowX = bookHisto1D(1, 1, 1); 00182 _histEnergyFlowHighX = bookHisto1D(1, 1, 2); 00183 00184 _histEECLowX = bookHisto1D(2, 1, 1); 00185 _histEECHighX = bookHisto1D(2, 1, 2); 00186 00187 _histSpectraW77 = bookHisto1D(3, 1, 1); 00188 _histSpectraW122 = bookHisto1D(3, 1, 2); 00189 _histSpectraW169 = bookHisto1D(3, 1, 3); 00190 _histSpectraW117 = bookHisto1D(3, 1, 4); 00191 00192 _histPT2 = bookProfile1D(4, 1, 1); 00193 } 00194 00195 00196 /// Finalize 00197 void finalize() { 00198 // Normalize inclusive single particle distributions to the average number 00199 // of charged particles per event. 00200 double avgNumParts = _w77.first/_w77.second; 00201 normalize(_histSpectraW77, avgNumParts); 00202 00203 avgNumParts = _w122.first/_w122.second; 00204 normalize(_histSpectraW122, avgNumParts); 00205 00206 avgNumParts = _w169.first/_w169.second; 00207 normalize(_histSpectraW169, avgNumParts); 00208 00209 avgNumParts = _w117.first/_w117.second; 00210 normalize(_histSpectraW117, avgNumParts); 00211 00212 scale(_histEnergyFlowLowX , 1./_wEnergy.first ); 00213 scale(_histEnergyFlowHighX, 1./_wEnergy.second); 00214 00215 scale(_histEECLowX , 1./_wEnergy.first ); 00216 scale(_histEECHighX, 1./_wEnergy.second); 00217 } 00218 00219 00220 //@} 00221 00222 00223 private: 00224 00225 /** 00226 * Polar angle with right direction of the beam 00227 */ 00228 inline double beamAngle(const FourVector& v, const bool & order) { 00229 double thel = v.polarAngle()/degree; 00230 if(thel<0.) thel+=180.; 00231 if(!order) thel = 180.-thel; 00232 return thel; 00233 } 00234 00235 /// @name Histograms 00236 //@{ 00237 Histo1DPtr _histEnergyFlowLowX; 00238 Histo1DPtr _histEnergyFlowHighX; 00239 Histo1DPtr _histEECLowX; 00240 Histo1DPtr _histEECHighX; 00241 Histo1DPtr _histSpectraW77; 00242 Histo1DPtr _histSpectraW122; 00243 Histo1DPtr _histSpectraW169; 00244 Histo1DPtr _histSpectraW117; 00245 Profile1DPtr _histPT2; 00246 //@} 00247 00248 /// @name storage of weight to calculate averages for normalisation 00249 //@{ 00250 pair<double,double> _w77,_w122,_w169,_w117,_wEnergy; 00251 //@} 00252 }; 00253 00254 00255 00256 // The hook for the plugin system 00257 DECLARE_RIVET_PLUGIN(H1_1994_S2919893); 00258 00259 } Generated on Fri Oct 25 2013 12:41:45 for The Rivet MC analysis system by ![]() |