| 1 |
dnisson |
1.2 |
// TODO find why energies suddenly bigger
|
| 2 |
dnisson |
1.1 |
/* A simple jet-finding analyzer */
|
| 3 |
|
|
|
| 4 |
|
|
#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
|
| 5 |
|
|
|
| 6 |
|
|
#include "fastjet/ClusterSequence.hh"
|
| 7 |
|
|
#include "FWCore/ServiceRegistry/interface/Service.h"
|
| 8 |
dnisson |
1.2 |
#include "FWCore/MessageLogger/interface/MessageLogger.h"
|
| 9 |
dnisson |
1.1 |
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
|
| 10 |
|
|
|
| 11 |
dnisson |
1.2 |
#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
|
| 12 |
|
|
#include "DataFormats/Candidate/interface/Particle.h"
|
| 13 |
|
|
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
|
| 14 |
|
|
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
|
| 15 |
|
|
#include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
|
| 16 |
|
|
|
| 17 |
dnisson |
1.1 |
#include <map>
|
| 18 |
|
|
#include <vector>
|
| 19 |
|
|
#include <limits>
|
| 20 |
|
|
#include <cmath>
|
| 21 |
|
|
#include <cstdlib>
|
| 22 |
|
|
#include <fstream>
|
| 23 |
|
|
#include <sstream>
|
| 24 |
|
|
|
| 25 |
|
|
#include "TFormula.h"
|
| 26 |
|
|
#include "TF2.h"
|
| 27 |
|
|
|
| 28 |
|
|
#define PI 3.141593
|
| 29 |
|
|
|
| 30 |
|
|
using namespace std;
|
| 31 |
|
|
using namespace fastjet;
|
| 32 |
|
|
|
| 33 |
dnisson |
1.2 |
// Class to represent a "generic" particle, whether raw or reconstructed
|
| 34 |
|
|
class GenericParticle {
|
| 35 |
|
|
public:
|
| 36 |
|
|
GenericParticle(double __px, double __py, double __pz, double __e,
|
| 37 |
|
|
double __charge = 0.0)
|
| 38 |
|
|
: _px(__px), _py(__py), _pz(__pz), _e(__e), _charge(__charge) {
|
| 39 |
|
|
}
|
| 40 |
|
|
GenericParticle(const HepMC::GenParticle &genParticle)
|
| 41 |
|
|
: _charge(0.0) {
|
| 42 |
|
|
_px = genParticle.momentum().px();
|
| 43 |
|
|
_py = genParticle.momentum().py();
|
| 44 |
|
|
_pz = genParticle.momentum().pz();
|
| 45 |
|
|
_e = genParticle.momentum().e();
|
| 46 |
|
|
}
|
| 47 |
|
|
GenericParticle(const reco::PFCandidate &pfCandidate)
|
| 48 |
|
|
: _charge(0.0) {
|
| 49 |
|
|
_px = pfCandidate.px();
|
| 50 |
|
|
_py = pfCandidate.py();
|
| 51 |
|
|
_pz = pfCandidate.pz();
|
| 52 |
|
|
_e = pfCandidate.energy();
|
| 53 |
|
|
}
|
| 54 |
|
|
double px() {
|
| 55 |
|
|
return _px;
|
| 56 |
|
|
}
|
| 57 |
|
|
double py() {
|
| 58 |
|
|
return _py;
|
| 59 |
|
|
}
|
| 60 |
|
|
double pz() {
|
| 61 |
|
|
return _pz;
|
| 62 |
|
|
}
|
| 63 |
|
|
double e() {
|
| 64 |
|
|
return _e;
|
| 65 |
|
|
}
|
| 66 |
|
|
double charge() {
|
| 67 |
|
|
return _charge;
|
| 68 |
|
|
}
|
| 69 |
|
|
double eta() {
|
| 70 |
|
|
double theta = acos(_pz/sqrt(_px*_px + _py*_py + _pz*_pz));
|
| 71 |
|
|
return -log(tan(theta*0.5));
|
| 72 |
|
|
}
|
| 73 |
|
|
double phi() {
|
| 74 |
|
|
double phi0 = acos(_px/sqrt(_px*_px + _py*_py));
|
| 75 |
|
|
return _py > 0.0 ? phi0 : -phi0;
|
| 76 |
|
|
}
|
| 77 |
|
|
double m() {
|
| 78 |
|
|
return sqrt(_e*_e - _px*_px - _py*_py - _pz*_pz);
|
| 79 |
|
|
}
|
| 80 |
|
|
private:
|
| 81 |
|
|
double _px;
|
| 82 |
|
|
double _py;
|
| 83 |
|
|
double _pz;
|
| 84 |
|
|
double _e;
|
| 85 |
|
|
double _charge;
|
| 86 |
|
|
};
|
| 87 |
|
|
|
| 88 |
dnisson |
1.1 |
class JetFinderAnalyzer : public JetFitAnalyzer {
|
| 89 |
|
|
public:
|
| 90 |
|
|
struct jet {
|
| 91 |
|
|
double energy;
|
| 92 |
|
|
double eta;
|
| 93 |
|
|
double phi;
|
| 94 |
|
|
};
|
| 95 |
|
|
|
| 96 |
|
|
explicit JetFinderAnalyzer( const edm::ParameterSet&);
|
| 97 |
|
|
~JetFinderAnalyzer() {}
|
| 98 |
|
|
|
| 99 |
|
|
private:
|
| 100 |
|
|
static map<TH2 *, vector< vector<jet> > > unique_jets;
|
| 101 |
|
|
|
| 102 |
|
|
static double phi_cutoff_;
|
| 103 |
|
|
|
| 104 |
|
|
static double g2int(double xlo, double xhi, double ylo, double yhi,
|
| 105 |
|
|
double *pval) {
|
| 106 |
|
|
double sum1 = 0.0;
|
| 107 |
|
|
double sum2 = 0.0;
|
| 108 |
|
|
double xmid = 0.5 * (xlo + xhi);
|
| 109 |
|
|
double ymid = 0.5 * (ylo + yhi);
|
| 110 |
|
|
double xstep = (xhi - xlo) / 50.0;
|
| 111 |
|
|
double ystep = (yhi - ylo) / 50.0;
|
| 112 |
|
|
for (int i = 0; i < 50; i++) {
|
| 113 |
|
|
double x = (static_cast<double>(i) + 0.5) * xstep + xlo;
|
| 114 |
|
|
sum1 += xstep * jetfit::fit_fcn(x, ymid, pval);
|
| 115 |
|
|
}
|
| 116 |
|
|
for (int i = 0; i < 50; i++) {
|
| 117 |
|
|
double y = (static_cast<double>(i) + 0.5) * ystep + ylo;
|
| 118 |
|
|
sum2 += ystep * jetfit::fit_fcn(xmid, y, pval);
|
| 119 |
|
|
}
|
| 120 |
|
|
return sum1 * sum2;
|
| 121 |
|
|
}
|
| 122 |
|
|
|
| 123 |
|
|
static void jetfinder(TMinuit *gMinuit, TH2 *hist, int ngauss) {
|
| 124 |
|
|
double dist_sq = numeric_limits<double>::infinity();
|
| 125 |
|
|
unique_jets[hist].resize(ngauss);
|
| 126 |
|
|
int nbinsX = hist->GetXaxis()->GetNbins();
|
| 127 |
|
|
int nbinsY = hist->GetYaxis()->GetNbins();
|
| 128 |
|
|
double XbinSize = (hist->GetXaxis()->GetXmax()
|
| 129 |
|
|
- hist->GetXaxis()->GetXmin())
|
| 130 |
|
|
/ static_cast<double>(nbinsX);
|
| 131 |
|
|
double YbinSize = (hist->GetYaxis()->GetXmax()
|
| 132 |
|
|
- hist->GetYaxis()->GetXmin())
|
| 133 |
|
|
/ static_cast<double>(nbinsY);
|
| 134 |
|
|
for (int i = 0; i < ngauss; i++) {
|
| 135 |
|
|
double N, mu_x, mu_y, sig, err, lo, hi;
|
| 136 |
|
|
int iuint;
|
| 137 |
|
|
TString name;
|
| 138 |
|
|
gMinuit->mnpout(4*i, name, N, err, lo, hi, iuint);
|
| 139 |
|
|
gMinuit->mnpout(4*i + 1, name, mu_x, err, lo, hi, iuint);
|
| 140 |
|
|
gMinuit->mnpout(4*i + 2, name, mu_y, err, lo, hi, iuint);
|
| 141 |
|
|
gMinuit->mnpout(4*i + 3, name, sig, err, lo, hi, iuint);
|
| 142 |
|
|
for (int j = 0; j < i; j++) {
|
| 143 |
|
|
double N2, mu_x2, mu_y2, sig2;
|
| 144 |
|
|
gMinuit->mnpout(4*j, name, N2, err, lo, hi, iuint);
|
| 145 |
|
|
gMinuit->mnpout(4*j + 1, name, mu_x2, err, lo, hi, iuint);
|
| 146 |
|
|
gMinuit->mnpout(4*j + 2, name, mu_y2, err, lo, hi, iuint);
|
| 147 |
|
|
gMinuit->mnpout(4*j + 3, name, sig2, err, lo, hi, iuint);
|
| 148 |
|
|
double _dist_sq = (mu_x2 - mu_x)*(mu_x2 - mu_x)
|
| 149 |
|
|
+ (mu_y2 - mu_y)*(mu_y2 - mu_y);
|
| 150 |
|
|
if (_dist_sq < dist_sq)
|
| 151 |
|
|
dist_sq = _dist_sq;
|
| 152 |
|
|
}
|
| 153 |
|
|
|
| 154 |
|
|
jet j;
|
| 155 |
|
|
j.energy = N;
|
| 156 |
|
|
j.eta = mu_x; j.phi = mu_y;
|
| 157 |
|
|
unique_jets[hist][ngauss-1].push_back(j);
|
| 158 |
|
|
}
|
| 159 |
|
|
}
|
| 160 |
|
|
|
| 161 |
|
|
virtual void beginJob(const edm::EventSetup&);
|
| 162 |
|
|
virtual TH2D* make_histo(const edm::Event&, const edm::EventSetup&);
|
| 163 |
|
|
virtual jetfit::model_def& make_model_def(const edm::Event&,
|
| 164 |
|
|
const edm::EventSetup&,
|
| 165 |
|
|
TH2 *);
|
| 166 |
|
|
virtual void analyze_results(jetfit::results r,
|
| 167 |
|
|
std::vector<jetfit::trouble> t, TH2 *);
|
| 168 |
dnisson |
1.2 |
vector<GenericParticle *> get_particles(const edm::Event&);
|
| 169 |
|
|
void fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>&,
|
| 170 |
|
|
const edm::InputTag&, const edm::Event&) const;
|
| 171 |
dnisson |
1.1 |
|
| 172 |
|
|
fstream ofs;
|
| 173 |
dnisson |
1.2 |
edm::InputTag inputTagPFCandidates_;
|
| 174 |
|
|
int info_type_;
|
| 175 |
dnisson |
1.1 |
double smear_;
|
| 176 |
|
|
int smear_coord_;
|
| 177 |
|
|
};
|
| 178 |
|
|
|
| 179 |
|
|
map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
|
| 180 |
|
|
JetFinderAnalyzer::unique_jets;
|
| 181 |
|
|
|
| 182 |
|
|
JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
|
| 183 |
|
|
: JetFitAnalyzer(pSet) // this is important!
|
| 184 |
|
|
{
|
| 185 |
dnisson |
1.2 |
info_type_ = pSet.getUntrackedParameter("info_type", 0);
|
| 186 |
|
|
|
| 187 |
|
|
if (info_type_ == 1) {
|
| 188 |
|
|
inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
|
| 189 |
|
|
}
|
| 190 |
|
|
|
| 191 |
dnisson |
1.1 |
smear_ = pSet.getUntrackedParameter("smear", 0.02);
|
| 192 |
|
|
smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0);
|
| 193 |
|
|
// 0 = eta-phi smear
|
| 194 |
|
|
// 1 = proper angle smear
|
| 195 |
|
|
set_user_minuit(jetfinder);
|
| 196 |
|
|
}
|
| 197 |
|
|
|
| 198 |
dnisson |
1.2 |
void
|
| 199 |
|
|
JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>& c,
|
| 200 |
|
|
const edm::InputTag& tag,
|
| 201 |
|
|
const edm::Event& iEvent) const {
|
| 202 |
|
|
|
| 203 |
|
|
bool found = iEvent.getByLabel(tag, c);
|
| 204 |
|
|
|
| 205 |
|
|
if(!found ) {
|
| 206 |
|
|
ostringstream err;
|
| 207 |
|
|
err<<" cannot get PFCandidates: "
|
| 208 |
|
|
<<tag<<endl;
|
| 209 |
|
|
edm::LogError("PFCandidates")<<err.str();
|
| 210 |
|
|
throw cms::Exception( "MissingProduct", err.str());
|
| 211 |
|
|
}
|
| 212 |
|
|
|
| 213 |
|
|
}
|
| 214 |
|
|
|
| 215 |
|
|
vector<GenericParticle *> JetFinderAnalyzer::get_particles(const edm::Event &evt) {
|
| 216 |
|
|
// fill unreduced histo
|
| 217 |
|
|
edm::Handle<edm::HepMCProduct> hRaw;
|
| 218 |
|
|
edm::Handle<reco::PFCandidateCollection> hPFlow;
|
| 219 |
|
|
if (info_type_ == 0) {
|
| 220 |
|
|
evt.getByLabel("source", hRaw);
|
| 221 |
|
|
}
|
| 222 |
|
|
if (info_type_ == 1) {
|
| 223 |
|
|
fetchCandidateCollection(hPFlow,
|
| 224 |
|
|
inputTagPFCandidates_,
|
| 225 |
|
|
evt);
|
| 226 |
|
|
}
|
| 227 |
|
|
|
| 228 |
|
|
vector<GenericParticle *> particles;
|
| 229 |
|
|
|
| 230 |
|
|
switch (info_type_) {
|
| 231 |
|
|
case 0:
|
| 232 |
|
|
const HepMC::GenEvent *hmcEvt = hRaw->GetEvent();
|
| 233 |
|
|
for (HepMC::GenEvent::particle_const_iterator
|
| 234 |
|
|
pit = hmcEvt->particles_begin(); pit != hmcEvt->particles_end();
|
| 235 |
|
|
pit++) {
|
| 236 |
|
|
if ((*pit)->status() == 1) {
|
| 237 |
|
|
particles.push_back(new GenericParticle(**pit));
|
| 238 |
|
|
}
|
| 239 |
|
|
}
|
| 240 |
|
|
|
| 241 |
|
|
break;
|
| 242 |
|
|
case 1:
|
| 243 |
|
|
for (unsigned i = 0; i < hPFlow->size(); i++) {
|
| 244 |
|
|
particles.push_back(new GenericParticle((*hPFlow)[i]));
|
| 245 |
|
|
}
|
| 246 |
|
|
break;
|
| 247 |
|
|
default:
|
| 248 |
|
|
cerr << "Unknown event type" << endl; // TODO use MessageLogger
|
| 249 |
|
|
}
|
| 250 |
|
|
|
| 251 |
|
|
return particles;
|
| 252 |
|
|
}
|
| 253 |
|
|
|
| 254 |
dnisson |
1.1 |
TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
|
| 255 |
|
|
ostringstream oss;
|
| 256 |
|
|
oss << "eta_phi_energy_unred"<<evt.id().event() << flush;
|
| 257 |
|
|
TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(),
|
| 258 |
|
|
600, -2.5, 2.5, 600, -PI, PI);
|
| 259 |
|
|
|
| 260 |
dnisson |
1.2 |
vector<GenericParticle *> particles = get_particles(evt);
|
| 261 |
dnisson |
1.1 |
for (int i = 0; i < particles.size(); i++) {
|
| 262 |
dnisson |
1.2 |
unred_histo->Fill(particles[i]->eta(),
|
| 263 |
|
|
particles[i]->phi(),
|
| 264 |
|
|
particles[i]->e());
|
| 265 |
dnisson |
1.1 |
}
|
| 266 |
|
|
|
| 267 |
|
|
// reduce histo
|
| 268 |
|
|
ostringstream oss2;
|
| 269 |
|
|
oss2 << "eta_phi_energy_red"<<evt.id().event() << flush;
|
| 270 |
|
|
edm::Service<TFileService> fs;
|
| 271 |
|
|
// draw cone of radius 0.5 around highest energy bin, reduce
|
| 272 |
|
|
double maxE = 0.0;
|
| 273 |
|
|
int max_i = 29, max_j = 29;
|
| 274 |
|
|
for (int i = 0; i < unred_histo->GetNbinsX(); i++) {
|
| 275 |
|
|
for (int j = 0; j < unred_histo->GetNbinsY(); j++) {
|
| 276 |
|
|
double E = unred_histo->GetBinContent(i+1, j+1);
|
| 277 |
|
|
if (E > maxE) {
|
| 278 |
|
|
maxE = E;
|
| 279 |
|
|
max_i = i;
|
| 280 |
|
|
max_j = j;
|
| 281 |
|
|
}
|
| 282 |
|
|
}
|
| 283 |
|
|
}
|
| 284 |
|
|
|
| 285 |
|
|
double rcone = 0.5;
|
| 286 |
|
|
double Xlo = unred_histo->GetXaxis()->GetXmin();
|
| 287 |
|
|
double Xhi = unred_histo->GetXaxis()->GetXmax();
|
| 288 |
|
|
double Ylo = unred_histo->GetYaxis()->GetXmin();
|
| 289 |
|
|
double Yhi = unred_histo->GetYaxis()->GetXmax();
|
| 290 |
|
|
double XbinSize = (Xhi - Xlo) /
|
| 291 |
|
|
static_cast<double>(unred_histo->GetXaxis()->GetNbins());
|
| 292 |
|
|
double YbinSize = (Yhi - Ylo) /
|
| 293 |
|
|
static_cast<double>(unred_histo->GetYaxis()->GetNbins());
|
| 294 |
|
|
double max_x = (static_cast<double>(max_i) + 0.5) * XbinSize + Xlo;
|
| 295 |
|
|
double max_y = (static_cast<double>(max_j) + 0.5) * YbinSize + Ylo;
|
| 296 |
|
|
TH2D *histo = fs->make<TH2D>(oss2.str().c_str(), oss2.str().c_str(),
|
| 297 |
|
|
60, max_x-rcone, max_x+rcone,
|
| 298 |
|
|
60, max_y-rcone, max_y+rcone);
|
| 299 |
|
|
|
| 300 |
|
|
// create an unsmeared reduced histo
|
| 301 |
|
|
TH2D *histo_unsmeared = fs->make<TH2D>((oss2.str()+"_unsmeared").c_str(),
|
| 302 |
|
|
(oss2.str()+"_unsmeared").c_str(),
|
| 303 |
|
|
60, max_x-rcone, max_x+rcone,
|
| 304 |
|
|
60, max_y-rcone, max_y+rcone);
|
| 305 |
|
|
for (int i = 0; i < particles.size(); i++) {
|
| 306 |
dnisson |
1.2 |
double N = particles[i]->e();
|
| 307 |
|
|
double x = particles[i]->eta();
|
| 308 |
|
|
double y = particles[i]->phi();
|
| 309 |
dnisson |
1.1 |
histo_unsmeared->Fill(x, y, N);
|
| 310 |
|
|
}
|
| 311 |
|
|
|
| 312 |
|
|
// create a smeared reduced histo
|
| 313 |
|
|
// create a temporary 2D vector for smeared energies
|
| 314 |
|
|
XbinSize = (histo->GetXaxis()->GetXmax()
|
| 315 |
|
|
- histo->GetXaxis()->GetXmin()) /
|
| 316 |
|
|
static_cast<double>(histo->GetXaxis()->GetNbins());
|
| 317 |
|
|
YbinSize = (histo->GetYaxis()->GetXmax()
|
| 318 |
|
|
- histo->GetYaxis()->GetXmin()) /
|
| 319 |
|
|
static_cast<double>(histo->GetYaxis()->GetNbins());
|
| 320 |
|
|
vector< vector<double> > smeared(60, vector<double>(60, 0.0) );
|
| 321 |
|
|
switch (smear_coord_) {
|
| 322 |
|
|
case 1:
|
| 323 |
|
|
for (int i = 0; i < particles.size(); i++) {
|
| 324 |
dnisson |
1.2 |
double N = particles[i]->e();
|
| 325 |
|
|
double x = particles[i]->eta();
|
| 326 |
|
|
double y = particles[i]->phi();
|
| 327 |
dnisson |
1.1 |
// loop over bins and add Gaussian in proper angle to smeared
|
| 328 |
|
|
for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
|
| 329 |
|
|
for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
|
| 330 |
|
|
double eta = static_cast<double>((signed int)i2) * XbinSize +
|
| 331 |
|
|
max_x - rcone - x;
|
| 332 |
|
|
double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize +
|
| 333 |
|
|
max_y - rcone - y));
|
| 334 |
|
|
phi = sin(phi) > 0 ? phi : -phi;
|
| 335 |
|
|
|
| 336 |
|
|
// transform eta, phi to proper angle
|
| 337 |
|
|
double theta = 2.0*atan(exp(-eta));
|
| 338 |
|
|
double iota = asin(sin(theta)*sin(phi));
|
| 339 |
|
|
|
| 340 |
|
|
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
|
| 341 |
|
|
* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
|
| 342 |
|
|
}
|
| 343 |
|
|
}
|
| 344 |
|
|
}
|
| 345 |
|
|
break;
|
| 346 |
|
|
case 0:
|
| 347 |
|
|
default:
|
| 348 |
|
|
for (int i = 0; i < particles.size(); i++) {
|
| 349 |
dnisson |
1.2 |
double N = particles[i]->e();
|
| 350 |
|
|
double x = particles[i]->eta();
|
| 351 |
|
|
double y = particles[i]->phi();
|
| 352 |
dnisson |
1.1 |
// loop over bins and add Gaussian to smeared
|
| 353 |
|
|
for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) {
|
| 354 |
|
|
for (vector< double >::size_type j2 = 0; j2 < 60; j2++) {
|
| 355 |
|
|
double eta = static_cast<double>((signed int)i2) * XbinSize
|
| 356 |
|
|
+ max_x - rcone - x;
|
| 357 |
|
|
double phi = acos(cos(static_cast<double>((signed int)j2) * YbinSize
|
| 358 |
|
|
+ max_y - rcone - y));
|
| 359 |
|
|
phi = sin(phi) > 0 ? phi : -phi;
|
| 360 |
|
|
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
|
| 361 |
|
|
* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
|
| 362 |
|
|
}
|
| 363 |
|
|
}
|
| 364 |
|
|
}
|
| 365 |
|
|
}
|
| 366 |
|
|
// set histogram to match smear vector
|
| 367 |
|
|
for (int i = 1; i <= 60; i++) {
|
| 368 |
|
|
for (int j = 1; j <= 60; j++) {
|
| 369 |
|
|
histo->SetBinContent(i, j, smeared[i-1][j-1]);
|
| 370 |
|
|
}
|
| 371 |
|
|
}
|
| 372 |
|
|
|
| 373 |
|
|
return histo;
|
| 374 |
|
|
}
|
| 375 |
|
|
|
| 376 |
dnisson |
1.2 |
void seed_with_CA(vector<GenericParticle *> gParticles, TH2 *histo,
|
| 377 |
dnisson |
1.1 |
jetfit::model_def &_mdef) {
|
| 378 |
|
|
// create a PseudoJet vector
|
| 379 |
|
|
vector<PseudoJet> particles;
|
| 380 |
dnisson |
1.2 |
for (unsigned i = 0; i < gParticles.size(); i++) {
|
| 381 |
|
|
double x_max = (histo->GetXaxis()->GetXmax()
|
| 382 |
|
|
+ histo->GetXaxis()->GetXmin()) / 2.0;
|
| 383 |
|
|
double y_max = (histo->GetYaxis()->GetXmax()
|
| 384 |
|
|
+ histo->GetYaxis()->GetXmin()) / 2.0;
|
| 385 |
|
|
valarray<double> pmom(4);
|
| 386 |
|
|
pmom[0] = gParticles[i]->px();
|
| 387 |
|
|
pmom[1] = gParticles[i]->py();
|
| 388 |
|
|
pmom[2] = gParticles[i]->pz();
|
| 389 |
|
|
pmom[3] = gParticles[i]->e();
|
| 390 |
|
|
double eta = gParticles[i]->eta();
|
| 391 |
|
|
double phi = gParticles[i]->phi();
|
| 392 |
|
|
if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) {
|
| 393 |
|
|
PseudoJet j(pmom);
|
| 394 |
|
|
particles.push_back(j);
|
| 395 |
dnisson |
1.1 |
}
|
| 396 |
|
|
}
|
| 397 |
|
|
|
| 398 |
|
|
// choose a jet definition
|
| 399 |
|
|
double R = 0.2;
|
| 400 |
|
|
JetDefinition jet_def(cambridge_algorithm, R);
|
| 401 |
|
|
|
| 402 |
|
|
// run clustering and extract the jets
|
| 403 |
|
|
ClusterSequence cs(particles, jet_def);
|
| 404 |
|
|
vector<PseudoJet> jets = cs.inclusive_jets();
|
| 405 |
|
|
|
| 406 |
|
|
double XbinSize = (histo->GetXaxis()->GetXmax()
|
| 407 |
|
|
- histo->GetXaxis()->GetXmin()) /
|
| 408 |
|
|
static_cast<double>(histo->GetXaxis()->GetNbins());
|
| 409 |
|
|
double YbinSize = (histo->GetYaxis()->GetXmax()
|
| 410 |
|
|
- histo->GetYaxis()->GetXmin()) /
|
| 411 |
|
|
static_cast<double>(histo->GetYaxis()->GetNbins());
|
| 412 |
|
|
|
| 413 |
|
|
// seed with C-A jets
|
| 414 |
|
|
int ijset = 0;
|
| 415 |
|
|
for (unsigned ij = 0; ij < jets.size(); ij++) {
|
| 416 |
|
|
double N = jets[ij].e();
|
| 417 |
dnisson |
1.2 |
if (N > 50.0) {
|
| 418 |
dnisson |
1.1 |
_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
|
| 419 |
|
|
0.0, 1.0e6);
|
| 420 |
|
|
_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01,
|
| 421 |
|
|
0.0, 0.0);
|
| 422 |
|
|
double mdef_phi = jets[ij].phi() > PI ? jets[ij].phi() - 2*PI
|
| 423 |
|
|
: jets[ij].phi();
|
| 424 |
|
|
_mdef.set_special_par(ijset, 2, mdef_phi, 0.01,
|
| 425 |
|
|
0.0, 0.0);
|
| 426 |
|
|
_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
|
| 427 |
|
|
ijset++;
|
| 428 |
|
|
}
|
| 429 |
|
|
}
|
| 430 |
|
|
}
|
| 431 |
|
|
|
| 432 |
|
|
jetfit::model_def& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
|
| 433 |
|
|
const edm::EventSetup&,
|
| 434 |
|
|
TH2 *histo) {
|
| 435 |
|
|
class jf_model_def : public jetfit::model_def {
|
| 436 |
|
|
public:
|
| 437 |
|
|
virtual double chisquare_error(double E) {
|
| 438 |
|
|
return 0.97*E + 14.0;
|
| 439 |
|
|
// study from 08-27-09
|
| 440 |
|
|
}
|
| 441 |
|
|
};
|
| 442 |
|
|
|
| 443 |
|
|
jf_model_def *_mdef = new jf_model_def();
|
| 444 |
|
|
TFormula *formula = new TFormula("gaus2d",
|
| 445 |
|
|
"[0]*exp(-0.5*((x-[1])**2 + (y-[2])**2)/([3]**2))/(2*pi*[3]**2)");
|
| 446 |
|
|
_mdef->set_formula(formula);
|
| 447 |
|
|
_mdef->set_indiv_max_E(0);
|
| 448 |
|
|
_mdef->set_indiv_max_x(1);
|
| 449 |
|
|
_mdef->set_indiv_max_y(2);
|
| 450 |
|
|
_mdef->set_indiv_par(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
|
| 451 |
|
|
_mdef->set_indiv_par(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
|
| 452 |
|
|
_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
|
| 453 |
|
|
_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
|
| 454 |
|
|
|
| 455 |
dnisson |
1.2 |
seed_with_CA(get_particles(evt), histo, *_mdef);
|
| 456 |
dnisson |
1.1 |
|
| 457 |
|
|
jetfit::set_model_def(_mdef);
|
| 458 |
|
|
|
| 459 |
|
|
// generate initial fit histogram
|
| 460 |
|
|
edm::Service<TFileService> fs;
|
| 461 |
|
|
TH2D *init_fit_histo = fs->make<TH2D>(("init_fit_"+string(histo->GetName()))
|
| 462 |
|
|
.c_str(),
|
| 463 |
|
|
("Initial fit for "
|
| 464 |
|
|
+string(histo->GetName())).c_str(),
|
| 465 |
|
|
histo->GetXaxis()->GetNbins(),
|
| 466 |
|
|
histo->GetXaxis()->GetXmin(),
|
| 467 |
|
|
histo->GetXaxis()->GetXmax(),
|
| 468 |
|
|
histo->GetXaxis()->GetNbins(),
|
| 469 |
|
|
histo->GetXaxis()->GetXmin(),
|
| 470 |
|
|
histo->GetXaxis()->GetXmax());
|
| 471 |
|
|
double XbinSize = (histo->GetXaxis()->GetXmax()
|
| 472 |
|
|
- histo->GetXaxis()->GetXmin()) /
|
| 473 |
|
|
static_cast<double>(histo->GetXaxis()->GetNbins());
|
| 474 |
|
|
double YbinSize = (histo->GetYaxis()->GetXmax()
|
| 475 |
|
|
- histo->GetYaxis()->GetXmin()) /
|
| 476 |
|
|
static_cast<double>(histo->GetYaxis()->GetNbins());
|
| 477 |
|
|
double Xlo = histo->GetXaxis()->GetXmin();
|
| 478 |
|
|
double Xhi = histo->GetXaxis()->GetXmax();
|
| 479 |
|
|
double Ylo = histo->GetYaxis()->GetXmin();
|
| 480 |
|
|
double Yhi = histo->GetYaxis()->GetXmax();
|
| 481 |
|
|
|
| 482 |
|
|
for (int i = 0; i < 60; i++) {
|
| 483 |
|
|
for (int j = 0; j < 60; j++) {
|
| 484 |
|
|
double x = (static_cast<double>(i) + 0.5)*XbinSize + Xlo;
|
| 485 |
|
|
double y = (static_cast<double>(j) + 0.5)*YbinSize + Ylo;
|
| 486 |
|
|
double pval[256];
|
| 487 |
|
|
if (_mdef->get_n_special_par_sets() > 64) {
|
| 488 |
|
|
cerr << "Parameter overload" << endl;
|
| 489 |
|
|
return *_mdef;
|
| 490 |
|
|
}
|
| 491 |
|
|
else {
|
| 492 |
|
|
for (int is = 0; is < _mdef->get_n_special_par_sets(); is++) {
|
| 493 |
|
|
for (int ii = 0; ii < 4; ii++) {
|
| 494 |
|
|
double spval, sperr, splo, sphi;
|
| 495 |
|
|
_mdef->get_special_par(is, ii, spval, sperr, splo, sphi);
|
| 496 |
|
|
pval[4*is + ii] = spval;
|
| 497 |
|
|
}
|
| 498 |
|
|
}
|
| 499 |
|
|
}
|
| 500 |
|
|
jetfit::set_ngauss(_mdef->get_n_special_par_sets());
|
| 501 |
|
|
init_fit_histo->SetBinContent(i+1, j+1,
|
| 502 |
|
|
jetfit::fit_fcn(x, y, pval));
|
| 503 |
|
|
}
|
| 504 |
|
|
}
|
| 505 |
|
|
|
| 506 |
|
|
return *_mdef;
|
| 507 |
|
|
}
|
| 508 |
|
|
|
| 509 |
|
|
void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
|
| 510 |
|
|
ofs.open("jetfindlog.txt", ios::out);
|
| 511 |
|
|
if (ofs.fail()) {
|
| 512 |
|
|
cerr << "Opening jetfindlog.txt FAILED" << endl;
|
| 513 |
|
|
}
|
| 514 |
|
|
ofs << "Jetfinder log" << endl
|
| 515 |
|
|
<< "=============" << endl << endl;
|
| 516 |
|
|
}
|
| 517 |
|
|
|
| 518 |
|
|
ostream& operator<<(ostream &out, jetfit::trouble t) {
|
| 519 |
|
|
string action, error_string;
|
| 520 |
|
|
|
| 521 |
|
|
if (t.istat != 3) {
|
| 522 |
|
|
switch(t.occ) {
|
| 523 |
|
|
case jetfit::T_NULL:
|
| 524 |
|
|
action = "Program"; break;
|
| 525 |
|
|
case jetfit::T_SIMPLEX:
|
| 526 |
|
|
action = "SIMPLEX"; break;
|
| 527 |
|
|
case jetfit::T_MIGRAD:
|
| 528 |
|
|
action = "MIGRAD"; break;
|
| 529 |
|
|
case jetfit::T_MINOS:
|
| 530 |
|
|
action = "MINOS"; break;
|
| 531 |
|
|
default:
|
| 532 |
|
|
action = "Program"; break;
|
| 533 |
|
|
}
|
| 534 |
|
|
|
| 535 |
|
|
switch (t.istat) {
|
| 536 |
|
|
case 0:
|
| 537 |
|
|
error_string = "Unable to calculate error matrix"; break;
|
| 538 |
|
|
case 1:
|
| 539 |
|
|
error_string = "Error matrix a diagonal approximation"; break;
|
| 540 |
|
|
case 2:
|
| 541 |
|
|
error_string = "Error matrix not positive definite"; break;
|
| 542 |
|
|
case 3:
|
| 543 |
|
|
error_string = "Converged successfully"; break;
|
| 544 |
|
|
default:
|
| 545 |
|
|
ostringstream oss;
|
| 546 |
|
|
oss<<"Unknown status code "<<t.istat << flush;
|
| 547 |
|
|
error_string = oss.str(); break;
|
| 548 |
|
|
}
|
| 549 |
|
|
|
| 550 |
|
|
if (t.occ != jetfit::T_NULL)
|
| 551 |
|
|
out << action<<" trouble: "<<error_string;
|
| 552 |
|
|
else
|
| 553 |
|
|
out << "Not calculated" << endl;
|
| 554 |
|
|
}
|
| 555 |
|
|
|
| 556 |
|
|
return out;
|
| 557 |
|
|
}
|
| 558 |
|
|
|
| 559 |
|
|
void JetFinderAnalyzer::analyze_results(jetfit::results r,
|
| 560 |
|
|
std::vector<jetfit::trouble> t,
|
| 561 |
|
|
TH2 *hist_orig) {
|
| 562 |
|
|
ofs << "Histogram "<<hist_orig->GetName() << endl;
|
| 563 |
|
|
for (int i = 0; i < unique_jets[hist_orig].size(); i++) {
|
| 564 |
|
|
ofs << "For "<<i+1<<" gaussians: " << endl
|
| 565 |
|
|
<< t.at(i) << endl
|
| 566 |
|
|
<< unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
|
| 567 |
|
|
for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
|
| 568 |
|
|
jet _jet = unique_jets[hist_orig][i][j];
|
| 569 |
|
|
ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta
|
| 570 |
|
|
<< ", phi = "<<_jet.phi << endl;
|
| 571 |
|
|
}
|
| 572 |
|
|
ofs << endl;
|
| 573 |
|
|
}
|
| 574 |
|
|
ofs << endl;
|
| 575 |
|
|
|
| 576 |
|
|
// save fit function histograms to root file
|
| 577 |
|
|
edm::Service<TFileService> fs;
|
| 578 |
|
|
for (vector< vector<double> >::size_type i = 0;
|
| 579 |
|
|
i < r.pval.size(); i++) {
|
| 580 |
|
|
jetfit::set_ngauss(r.pval[i].size() / 4);
|
| 581 |
|
|
TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2,
|
| 582 |
|
|
hist_orig->GetXaxis()->GetXmin(),
|
| 583 |
|
|
hist_orig->GetXaxis()->GetXmax(),
|
| 584 |
|
|
hist_orig->GetYaxis()->GetXmin(),
|
| 585 |
|
|
hist_orig->GetYaxis()->GetXmax(),
|
| 586 |
|
|
r.pval[i].size());
|
| 587 |
|
|
for (vector<double>::size_type j = 0; j < r.pval[i].size(); j++) {
|
| 588 |
|
|
tf2->SetParameter(j, r.pval[i][j]);
|
| 589 |
|
|
}
|
| 590 |
|
|
ostringstream fit_histo_oss;
|
| 591 |
|
|
fit_histo_oss << hist_orig->GetName()<<"_fit_"<<i << flush;
|
| 592 |
|
|
tf2->SetNpx(hist_orig->GetXaxis()->GetNbins());
|
| 593 |
|
|
tf2->SetNpy(hist_orig->GetYaxis()->GetNbins());
|
| 594 |
|
|
TH2D *fit_histo = fs->make<TH2D>(fit_histo_oss.str().c_str(),
|
| 595 |
|
|
fit_histo_oss.str().c_str(),
|
| 596 |
|
|
hist_orig->GetXaxis()->GetNbins(),
|
| 597 |
|
|
hist_orig->GetXaxis()->GetXmin(),
|
| 598 |
|
|
hist_orig->GetXaxis()->GetXmax(),
|
| 599 |
|
|
hist_orig->GetYaxis()->GetNbins(),
|
| 600 |
|
|
hist_orig->GetYaxis()->GetXmin(),
|
| 601 |
|
|
hist_orig->GetYaxis()->GetXmax());
|
| 602 |
|
|
TH1 *tf2_histo = tf2->CreateHistogram();
|
| 603 |
|
|
double XbinSize = (fit_histo->GetXaxis()->GetXmax()
|
| 604 |
|
|
- fit_histo->GetXaxis()->GetXmin())
|
| 605 |
|
|
/ static_cast<double>(fit_histo->GetXaxis()->GetNbins());
|
| 606 |
|
|
double YbinSize = (fit_histo->GetYaxis()->GetXmax()
|
| 607 |
|
|
- fit_histo->GetYaxis()->GetXmin())
|
| 608 |
|
|
/ static_cast<double>(fit_histo->GetYaxis()->GetNbins());
|
| 609 |
|
|
for (int ih = 0; ih < tf2->GetNpx(); ih++) {
|
| 610 |
|
|
for (int jh = 0; jh < tf2->GetNpy(); jh++) {
|
| 611 |
|
|
fit_histo->SetBinContent(ih+1, jh+1,
|
| 612 |
|
|
tf2_histo->GetBinContent(ih+1, jh+1)
|
| 613 |
|
|
* XbinSize * YbinSize);
|
| 614 |
|
|
}
|
| 615 |
|
|
}
|
| 616 |
|
|
}
|
| 617 |
|
|
|
| 618 |
|
|
// save results to file
|
| 619 |
|
|
ostringstream res_tree_oss, rt_title_oss;
|
| 620 |
|
|
res_tree_oss << hist_orig->GetName()<<"_results" << flush;
|
| 621 |
|
|
rt_title_oss << "Fit results for "<<hist_orig->GetName() << flush;
|
| 622 |
|
|
}
|
| 623 |
|
|
|
| 624 |
|
|
DEFINE_FWK_MODULE(JetFinderAnalyzer);
|
