Page MenuHomeHEPForge
Files F11222053

No OneTemporary
Actions

View Options

Index: trunk/CMS_2012_I1090064.cc
===================================================================
--- trunk/CMS_2012_I1090064.cc (revision 178)
+++ trunk/CMS_2012_I1090064.cc (revision 179)
@@ -1,267 +1,267 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/RivetAIDA.hh"
#include <boost/lexical_cast.hpp>
#include <map>
namespace Rivet {
/// @brief CMS charged particle fragmentation function in PbPb
class CMS_2012_I1090064 : public Analysis {
public:
/// Constructor
CMS_2012_I1090064() : Analysis("CMS_2012_I1090064")
{
setBeams(PROTON, PROTON);
}
/// Book projections and histograms
void init() {
FinalState fs(-5.0, 5.0, 0.*GeV);
addProjection(fs, "FS");
ChargedFinalState cfs(-5.0, 5.0, 4.*GeV);
addProjection(cfs, "CFS");
FastJets fj03(fs, FastJets::ANTIKT, 0.3);
addProjection(fj03, "Jets03");
_h_deltaphi_1 = bookHistogram1D(1, 1, 1);
_h_deltaphi_2 = bookHistogram1D(1, 1, 2);
_h_deltaphi_3 = bookHistogram1D(1, 1, 3);
_h_deltaphi_4 = bookHistogram1D(1, 1, 4);
_h_deltaphi_5 = bookHistogram1D(1, 1, 5);
_h_deltaphi_6 = bookHistogram1D(1, 1, 6);
_h_dijets_p = bookHistogram1D("N_di-jets_pp", binEdges(2, 1, 1));
_h_leadjets_p = bookHistogram1D("N_lead-jets_pp", binEdges(2, 1, 1));
_h_dijets_pb = bookHistogram1D("N_di-jets_PbPb", binEdges(2, 1, 2));
_h_leadjets_pb= bookHistogram1D("N_lead-jets_PbPb", binEdges(2, 1, 2));
_h_aj_1 = bookHistogram1D(3, 1, 1);
_h_aj_2 = bookHistogram1D(3, 1, 2);
_h_aj_3 = bookHistogram1D(3, 1, 3);
_h_aj_4 = bookHistogram1D(3, 1, 4);
_h_aj_5 = bookHistogram1D(3, 1, 5);
_h_aj_6 = bookHistogram1D(3, 1, 6);
_h_aj_all = bookHistogram1D(4, 1, 1);
_h_ptratio1_1 = bookHistogram1D(5, 1, 1);
_h_ptratio1_2 = bookHistogram1D(5, 1, 2);
_h_ptratio1_3 = bookHistogram1D(5, 1, 3);
_h_ptratio1_4 = bookHistogram1D(5, 1, 4);
_h_ptratio1_5 = bookHistogram1D(5, 1, 5);
_h_ptratio1_6 = bookHistogram1D(5, 1, 6);
_p_ptratio2_p = bookProfile1D(6, 1, 1);
_p_ptratio2_pb= bookProfile1D(6, 1, 2);
}
struct jetcomp{
bool operator() (const Jet jet1, const Jet jet2) const
{
if (jet1.momentum().pT() < jet2.momentum().pT()) return true;
else return false;
}
};
double ptsmear(const double pttrue) {
double c(0.0246);
double s(1.213);
// double n(0.001);
// double n(3.88);
// double n(5.10);
double n(5.23);
double sigma(sqrt(c*c+s*s/pttrue+n*n/(pttrue*pttrue)));
double r1(1.0*rand()/RAND_MAX), r2(1.0*rand()/RAND_MAX);
double fac(max(sqrt(-2.*log(r1))*cos(2.*M_PI*r2)*sigma+1.,0.));
// std::cout<<sigma<<std::endl;
return fac*pttrue;
}
/// Do the analysis
void analyze(const Event& event) {
const double weight = event.weight();
const double centrality = (event.genEvent().heavy_ion()?event.genEvent().heavy_ion()->impact_parameter():-1.);
const Jets jets03 = applyProjection<FastJets>(event, "Jets03").jetsByPt(3.*GeV);
const ChargedFinalState tracks = applyProjection<ChargedFinalState>(event, "CFS");
map<Jet, double, jetcomp> jetmap;
foreach (Jet jet, jets03) {
double ptreco, ptgen;
ptgen = jet.momentum().pT();
ptreco = ptsmear(ptgen);
jetmap[jet] = ptreco;
}
Jets jets;
foreach (Jet jet, jets03) {
if (fabs(jet.momentum().eta()) < 2.0 && jetmap[jet] > 30.*GeV)
jets.push_back(jet);
}
if (jets.size() < 1) vetoEvent;
Jet leadjet, partner;
double ptlead(0.), ptsublead(0.);
foreach (Jet jet, jets) {
if (jetmap[jet] > ptlead) {
ptlead = jetmap[jet];
leadjet = jet;
}
}
if (ptlead < 120.*GeV) vetoEvent;
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_leadjets_p->fill(ptlead,weight);
if (centrality<0.) _h_leadjets_pb->fill(ptlead,weight);
if (jets.size() < 2) vetoEvent;
foreach (Jet jet, jets) {
if (jetmap[jet] > ptsublead && jet.momentum()!=leadjet.momentum()) {
ptsublead = jetmap[jet];
partner = jet;
}
}
double deltaphi = fabs(deltaPhi(leadjet,partner));
if (ptlead < 150*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_deltaphi_1->fill(deltaphi,weight);
}
else if (ptlead < 180*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_deltaphi_2->fill(deltaphi,weight);
}
else if (ptlead < 220*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_deltaphi_3->fill(deltaphi,weight);
}
else if (ptlead < 260*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_deltaphi_4->fill(deltaphi,weight);
}
else if (ptlead < 300*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_deltaphi_5->fill(deltaphi,weight);
}
else if (ptlead < 500*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_deltaphi_6->fill(deltaphi,weight);
}
if (deltaphi < 2.*M_PI/3.) vetoEvent;
if (centrality<0.) _h_dijets_p->fill(ptlead,weight);
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_dijets_pb->fill(ptlead,weight);
if (centrality<0.) _p_ptratio2_p->fill(ptlead,ptsublead/ptlead,weight);
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _p_ptratio2_pb->fill(ptlead,ptsublead/ptlead,weight);
Rivet::FourMomentum leading, sublead;
leading = leadjet.momentum();
sublead = partner.momentum();
double aj = (ptlead-ptsublead)/(ptlead+ptsublead);
if (centrality<0. || (centrality>0.0 && centrality<0.1)) _h_aj_all->fill(aj,weight);
if (ptlead < 150*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_aj_1->fill(aj,weight);
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_ptratio1_1->fill(ptsublead/ptlead,weight);
}
else if (ptlead < 180*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_aj_2->fill(aj,weight);
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_ptratio1_2->fill(ptsublead/ptlead,weight);
}
else if (ptlead < 220*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_aj_3->fill(aj,weight);
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_ptratio1_3->fill(ptsublead/ptlead,weight);
}
else if (ptlead < 260*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_aj_4->fill(aj,weight);
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_ptratio1_4->fill(ptsublead/ptlead,weight);
}
else if (ptlead < 300*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_aj_5->fill(aj,weight);
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_ptratio1_5->fill(ptsublead/ptlead,weight);
}
else if (ptlead < 500*GeV) {
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_aj_6->fill(aj,weight);
if (centrality<0. || (centrality>0.0 && centrality<0.2)) _h_ptratio1_6->fill(ptsublead/ptlead,weight);
}
}
/// Finalize
void finalize() {
double Dbin_Dphi, Dbin_aj, Dbin_ptratio;
Dbin_Dphi = M_PI/30.;
- Dbin_aj = 1./20.;
+ Dbin_aj = 0.72/12.;
Dbin_ptratio = 1./10.;
normalize(_h_deltaphi_1,Dbin_Dphi);
normalize(_h_deltaphi_2,Dbin_Dphi);
normalize(_h_deltaphi_3,Dbin_Dphi);
normalize(_h_deltaphi_4,Dbin_Dphi);
normalize(_h_deltaphi_5,Dbin_Dphi);
normalize(_h_deltaphi_6,Dbin_Dphi);
normalize(_h_aj_all,Dbin_aj);
normalize(_h_aj_1,Dbin_aj);
normalize(_h_aj_2,Dbin_aj);
normalize(_h_aj_3,Dbin_aj);
normalize(_h_aj_4,Dbin_aj);
normalize(_h_aj_5,Dbin_aj);
normalize(_h_aj_6,Dbin_aj);
normalize(_h_ptratio1_1,Dbin_ptratio);
normalize(_h_ptratio1_2,Dbin_ptratio);
normalize(_h_ptratio1_3,Dbin_ptratio);
normalize(_h_ptratio1_4,Dbin_ptratio);
normalize(_h_ptratio1_5,Dbin_ptratio);
normalize(_h_ptratio1_6,Dbin_ptratio);
AIDA::IHistogramFactory& hf = histogramFactory();
const string dir = histoDir();
hf.divide(dir + "/d02-x01-y01", *_h_dijets_p, *_h_leadjets_p);
hf.destroy(_h_dijets_p);
hf.destroy(_h_leadjets_p);
hf.divide(dir + "/d02-x01-y02", *_h_dijets_pb, *_h_leadjets_pb);
hf.destroy(_h_dijets_pb);
hf.destroy(_h_leadjets_pb);
getLog() << Log::DEBUG << "sumOfWeights() = " << sumOfWeights() << std::endl;
}
private:
AIDA::IHistogram1D * _h_deltaphi_1;
AIDA::IHistogram1D * _h_deltaphi_2;
AIDA::IHistogram1D * _h_deltaphi_3;
AIDA::IHistogram1D * _h_deltaphi_4;
AIDA::IHistogram1D * _h_deltaphi_5;
AIDA::IHistogram1D * _h_deltaphi_6;
AIDA::IHistogram1D * _h_aj_all;
AIDA::IHistogram1D * _h_aj_1;
AIDA::IHistogram1D * _h_aj_2;
AIDA::IHistogram1D * _h_aj_3;
AIDA::IHistogram1D * _h_aj_4;
AIDA::IHistogram1D * _h_aj_5;
AIDA::IHistogram1D * _h_aj_6;
AIDA::IHistogram1D * _h_ptratio1_1;
AIDA::IHistogram1D * _h_ptratio1_2;
AIDA::IHistogram1D * _h_ptratio1_3;
AIDA::IHistogram1D * _h_ptratio1_4;
AIDA::IHistogram1D * _h_ptratio1_5;
AIDA::IHistogram1D * _h_ptratio1_6;
AIDA::IHistogram1D * _h_dijets_p;
AIDA::IHistogram1D * _h_leadjets_p;
AIDA::IHistogram1D * _h_dijets_pb;
AIDA::IHistogram1D * _h_leadjets_pb;
AIDA::IProfile1D * _p_ptratio2_p;
AIDA::IProfile1D * _p_ptratio2_pb;
};
// This global object acts as a hook for the plugin system
AnalysisBuilder<CMS_2012_I1090064> plugin_CMS_2012_I1090064;
}

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 11:20 AM (10 h, 57 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111400
Default Alt Text
(9 KB)

Event Timeline