| 13 |
|
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h" |
| 14 |
|
#include "DataFormats/ParticleFlowReco/interface/PFBlock.h" |
| 15 |
|
|
| 16 |
+ |
/* Jet reco stuff */ |
| 17 |
+ |
#include "DataFormats/JetReco/interface/PFJetCollection.h" |
| 18 |
+ |
#include "DataFormats/JetReco/interface/GenJetCollection.h" |
| 19 |
+ |
#include "DataFormats/JetReco/interface/PFJet.h" |
| 20 |
+ |
#include "DataFormats/JetReco/interface/GenJet.h" |
| 21 |
+ |
|
| 22 |
|
#include <map> |
| 23 |
|
#include <vector> |
| 24 |
|
#include <limits> |
| 70 |
|
return sum1 * sum2; |
| 71 |
|
} |
| 72 |
|
|
| 73 |
< |
static void jetfinder(TMinuit *gMinuit, TH2 *hist, int ngauss) { |
| 74 |
< |
double dist_sq = numeric_limits<double>::infinity(); |
| 69 |
< |
unique_jets[hist].resize(ngauss); |
| 70 |
< |
int nbinsX = hist->GetXaxis()->GetNbins(); |
| 71 |
< |
int nbinsY = hist->GetYaxis()->GetNbins(); |
| 72 |
< |
double XbinSize = (hist->GetXaxis()->GetXmax() |
| 73 |
< |
- hist->GetXaxis()->GetXmin()) |
| 74 |
< |
/ static_cast<double>(nbinsX); |
| 75 |
< |
double YbinSize = (hist->GetYaxis()->GetXmax() |
| 76 |
< |
- hist->GetYaxis()->GetXmin()) |
| 77 |
< |
/ static_cast<double>(nbinsY); |
| 78 |
< |
for (int i = 0; i < ngauss; i++) { |
| 79 |
< |
double N, mu_x, mu_y, sig, err, lo, hi; |
| 80 |
< |
int iuint; |
| 81 |
< |
TString name; |
| 82 |
< |
gMinuit->mnpout(4*i, name, N, err, lo, hi, iuint); |
| 83 |
< |
gMinuit->mnpout(4*i + 1, name, mu_x, err, lo, hi, iuint); |
| 84 |
< |
gMinuit->mnpout(4*i + 2, name, mu_y, err, lo, hi, iuint); |
| 85 |
< |
gMinuit->mnpout(4*i + 3, name, sig, err, lo, hi, iuint); |
| 86 |
< |
for (int j = 0; j < i; j++) { |
| 87 |
< |
double N2, mu_x2, mu_y2, sig2; |
| 88 |
< |
gMinuit->mnpout(4*j, name, N2, err, lo, hi, iuint); |
| 89 |
< |
gMinuit->mnpout(4*j + 1, name, mu_x2, err, lo, hi, iuint); |
| 90 |
< |
gMinuit->mnpout(4*j + 2, name, mu_y2, err, lo, hi, iuint); |
| 91 |
< |
gMinuit->mnpout(4*j + 3, name, sig2, err, lo, hi, iuint); |
| 92 |
< |
double _dist_sq = (mu_x2 - mu_x)*(mu_x2 - mu_x) |
| 93 |
< |
+ (mu_y2 - mu_y)*(mu_y2 - mu_y); |
| 94 |
< |
if (_dist_sq < dist_sq) |
| 95 |
< |
dist_sq = _dist_sq; |
| 96 |
< |
} |
| 97 |
< |
|
| 98 |
< |
jet j; |
| 99 |
< |
j.energy = N; |
| 100 |
< |
j.eta = mu_x; j.phi = mu_y; |
| 101 |
< |
unique_jets[hist][ngauss-1].push_back(j); |
| 102 |
< |
} |
| 103 |
< |
} |
| 104 |
< |
|
| 105 |
< |
virtual void beginJob(const edm::EventSetup&); |
| 106 |
< |
virtual TH2D* make_histo(const edm::Event&, const edm::EventSetup&); |
| 73 |
> |
virtual void beginJob(const edm::EventSetup &es); |
| 74 |
> |
virtual TH2D * make_histo(const edm::Event &, const edm::EventSetup &); |
| 75 |
|
virtual jetfit::model_def& make_model_def(const edm::Event&, |
| 76 |
|
const edm::EventSetup&, |
| 77 |
|
TH2 *); |
| 89 |
|
int info_type_; |
| 90 |
|
double smear_; |
| 91 |
|
int smear_coord_; |
| 92 |
+ |
string jet_algo_; |
| 93 |
+ |
string reclustered_jet_algo_; |
| 94 |
|
}; |
| 95 |
|
|
| 96 |
|
map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > > |
| 112 |
|
smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0); |
| 113 |
|
// 0 = eta-phi smear |
| 114 |
|
// 1 = proper angle smear |
| 115 |
< |
set_user_minuit(jetfinder); |
| 115 |
> |
jet_algo_ = pSet.getParameter<string>("jet_algo"); |
| 116 |
> |
set_user_minuit(0); // we don't need this feature--it's deprecated |
| 117 |
|
} |
| 118 |
|
|
| 119 |
|
void |
| 150 |
|
|
| 151 |
|
} |
| 152 |
|
|
| 153 |
< |
vector<reco::Candidate *> JetFinderAnalyzer::get_particles(const edm::Event &evt) { |
| 154 |
< |
// fill unreduced histo |
| 184 |
< |
edm::Handle<reco::GenParticleCollection> hRaw; |
| 185 |
< |
edm::Handle<reco::PFCandidateCollection> hPFlow; |
| 186 |
< |
if (info_type_ == 0) { |
| 187 |
< |
fetchCandidateCollection(hRaw, |
| 188 |
< |
inputTagGenParticles_, |
| 189 |
< |
evt); |
| 190 |
< |
} |
| 191 |
< |
if (info_type_ == 1) { |
| 192 |
< |
fetchCandidateCollection(hPFlow, |
| 193 |
< |
inputTagPFCandidates_, |
| 194 |
< |
evt); |
| 195 |
< |
} |
| 196 |
< |
|
| 197 |
< |
vector<reco::Candidate *> particles; |
| 153 |
> |
TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) { |
| 154 |
> |
const reco::Jet * highest_e_jet; |
| 155 |
|
|
| 156 |
< |
switch (info_type_) { |
| 157 |
< |
case 0: |
| 158 |
< |
for (unsigned i = 0; i < hRaw->size(); i++) { |
| 159 |
< |
if ((*hRaw)[i].status() == 1) { |
| 160 |
< |
particles.push_back((reco::Candidate *)&((*hRaw)[i])); |
| 156 |
> |
if (info_type_ == 0) { |
| 157 |
> |
edm::Handle< vector<reco::GenJet> > evtJets; |
| 158 |
> |
evt.getByLabel(jet_algo_, evtJets); |
| 159 |
> |
for (unsigned i = 0; i < evtJets->size(); i++) { |
| 160 |
> |
if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) { |
| 161 |
> |
highest_e_jet = &((*evtJets)[i]); |
| 162 |
|
} |
| 163 |
|
} |
| 164 |
< |
break; |
| 165 |
< |
case 1: |
| 166 |
< |
for (unsigned i = 0; i < hPFlow->size(); i++) { |
| 167 |
< |
particles.push_back((reco::Candidate *)&((*hPFlow)[i])); |
| 164 |
> |
} |
| 165 |
> |
else if (info_type_ == 1) { |
| 166 |
> |
edm::Handle< vector<reco::PFJet> > evtJets; |
| 167 |
> |
evt.getByLabel(jet_algo_, evtJets); |
| 168 |
> |
for (unsigned i = 0; i < evtJets->size(); i++) { |
| 169 |
> |
if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) { |
| 170 |
> |
highest_e_jet = &((*evtJets)[i]); |
| 171 |
> |
} |
| 172 |
|
} |
| 173 |
< |
break; |
| 174 |
< |
default: |
| 175 |
< |
cerr << "Unknown event type" << endl; // TODO use MessageLogger |
| 173 |
> |
} |
| 174 |
> |
cout << "found highest e jet" << endl; |
| 175 |
> |
|
| 176 |
> |
if (highest_e_jet == 0) { |
| 177 |
> |
cerr << "No fatjets found!" << endl; exit(1); |
| 178 |
|
} |
| 179 |
|
|
| 180 |
< |
return particles; |
| 181 |
< |
} |
| 180 |
> |
vector<const reco::Candidate *> particles = |
| 181 |
> |
highest_e_jet->getJetConstituentsQuick(); |
| 182 |
> |
cout << "got cands list of highest e jet" << endl; |
| 183 |
|
|
| 219 |
– |
TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) { |
| 184 |
|
ostringstream oss; |
| 185 |
< |
oss << "eta_phi_energy_unred"<<evt.id().event() << flush; |
| 222 |
< |
TH2D *unred_histo = new TH2D(oss.str().c_str(), oss.str().c_str(), |
| 223 |
< |
600, -2.5, 2.5, 600, -PI, PI); |
| 224 |
< |
|
| 225 |
< |
vector<reco::Candidate *> particles = get_particles(evt); |
| 226 |
< |
for (int i = 0; i < particles.size(); i++) { |
| 227 |
< |
unred_histo->Fill(particles[i]->eta(), |
| 228 |
< |
particles[i]->phi(), |
| 229 |
< |
particles[i]->energy()); |
| 230 |
< |
} |
| 231 |
< |
|
| 232 |
< |
// reduce histo |
| 233 |
< |
ostringstream oss2; |
| 234 |
< |
oss2 << "eta_phi_energy_red"<<evt.id().event() << flush; |
| 185 |
> |
oss << "eta_phi_energy"<<evt.id().event() << flush; |
| 186 |
|
edm::Service<TFileService> fs; |
| 187 |
< |
// draw cone of radius 0.5 around highest energy bin, reduce |
| 188 |
< |
double maxE = 0.0; |
| 189 |
< |
int max_i = 29, max_j = 29; |
| 190 |
< |
for (int i = 0; i < unred_histo->GetNbinsX(); i++) { |
| 191 |
< |
for (int j = 0; j < unred_histo->GetNbinsY(); j++) { |
| 192 |
< |
double E = unred_histo->GetBinContent(i+1, j+1); |
| 193 |
< |
if (E > maxE) { |
| 243 |
< |
maxE = E; |
| 244 |
< |
max_i = i; |
| 245 |
< |
max_j = j; |
| 246 |
< |
} |
| 247 |
< |
} |
| 248 |
< |
} |
| 249 |
< |
|
| 250 |
< |
double rcone = 0.5; |
| 251 |
< |
double Xlo = unred_histo->GetXaxis()->GetXmin(); |
| 252 |
< |
double Xhi = unred_histo->GetXaxis()->GetXmax(); |
| 253 |
< |
double Ylo = unred_histo->GetYaxis()->GetXmin(); |
| 254 |
< |
double Yhi = unred_histo->GetYaxis()->GetXmax(); |
| 255 |
< |
double XbinSize = (Xhi - Xlo) / |
| 256 |
< |
static_cast<double>(unred_histo->GetXaxis()->GetNbins()); |
| 257 |
< |
double YbinSize = (Yhi - Ylo) / |
| 258 |
< |
static_cast<double>(unred_histo->GetYaxis()->GetNbins()); |
| 259 |
< |
double max_x = (static_cast<double>(max_i) + 0.5) * XbinSize + Xlo; |
| 260 |
< |
double max_y = (static_cast<double>(max_j) + 0.5) * YbinSize + Ylo; |
| 261 |
< |
TH2D *histo = fs->make<TH2D>(oss2.str().c_str(), oss2.str().c_str(), |
| 262 |
< |
60, max_x-rcone, max_x+rcone, |
| 263 |
< |
60, max_y-rcone, max_y+rcone); |
| 264 |
< |
|
| 265 |
< |
// create an unsmeared reduced histo |
| 266 |
< |
TH2D *histo_unsmeared = fs->make<TH2D>((oss2.str()+"_unsmeared").c_str(), |
| 267 |
< |
(oss2.str()+"_unsmeared").c_str(), |
| 268 |
< |
60, max_x-rcone, max_x+rcone, |
| 269 |
< |
60, max_y-rcone, max_y+rcone); |
| 270 |
< |
for (int i = 0; i < particles.size(); i++) { |
| 271 |
< |
double N = particles[i]->energy(); |
| 272 |
< |
double x = particles[i]->eta(); |
| 273 |
< |
double y = particles[i]->phi(); |
| 274 |
< |
histo_unsmeared->Fill(x, y, N); |
| 275 |
< |
} |
| 187 |
> |
TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(), |
| 188 |
> |
30, |
| 189 |
> |
highest_e_jet->eta()-0.5, |
| 190 |
> |
highest_e_jet->eta()+0.5, |
| 191 |
> |
30, |
| 192 |
> |
highest_e_jet->phi()-0.5, |
| 193 |
> |
highest_e_jet->phi()+0.5); |
| 194 |
|
|
| 195 |
< |
// create a smeared reduced histo |
| 195 |
> |
// create a smeared histo |
| 196 |
|
// create a temporary 2D vector for smeared energies |
| 197 |
< |
XbinSize = (histo->GetXaxis()->GetXmax() |
| 197 |
> |
double XbinSize = (histo->GetXaxis()->GetXmax() |
| 198 |
|
- histo->GetXaxis()->GetXmin()) / |
| 199 |
|
static_cast<double>(histo->GetXaxis()->GetNbins()); |
| 200 |
< |
YbinSize = (histo->GetYaxis()->GetXmax() |
| 200 |
> |
double YbinSize = (histo->GetYaxis()->GetXmax() |
| 201 |
|
- histo->GetYaxis()->GetXmin()) / |
| 202 |
|
static_cast<double>(histo->GetYaxis()->GetNbins()); |
| 203 |
< |
vector< vector<double> > smeared(60, vector<double>(60, 0.0) ); |
| 203 |
> |
double Xlo = histo->GetXaxis()->GetXmin(); |
| 204 |
> |
double Xhi = histo->GetXaxis()->GetXmax(); |
| 205 |
> |
double Ylo = histo->GetYaxis()->GetXmin(); |
| 206 |
> |
double Yhi = histo->GetYaxis()->GetXmax(); |
| 207 |
> |
vector< vector<double> > smeared(30, vector<double>(30, 0.0) ); |
| 208 |
|
switch (smear_coord_) { |
| 209 |
|
case 1: |
| 210 |
|
for (int i = 0; i < particles.size(); i++) { |
| 211 |
|
double N = particles[i]->energy(); |
| 212 |
|
double x = particles[i]->eta(); |
| 213 |
|
double y = particles[i]->phi(); |
| 214 |
+ |
if (y < Ylo) { |
| 215 |
+ |
y += 2.0*PI; |
| 216 |
+ |
} |
| 217 |
+ |
if (y > Yhi) { |
| 218 |
+ |
y -= 2.0*PI; |
| 219 |
+ |
} |
| 220 |
|
// loop over bins and add Gaussian in proper angle to smeared |
| 221 |
< |
for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) { |
| 222 |
< |
for (vector< double >::size_type j2 = 0; j2 < 60; j2++) { |
| 221 |
> |
for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) { |
| 222 |
> |
for (vector< double >::size_type j2 = 0; j2 < 30; j2++) { |
| 223 |
|
double eta = static_cast<double>((signed int)i2) * XbinSize + |
| 224 |
< |
max_x - rcone - x; |
| 225 |
< |
double phi0 = static_cast<double>((signed int)j2) * YbinSize + |
| 226 |
< |
max_y - rcone - y; |
| 227 |
< |
double phi = acos(cos(phi0)); |
| 300 |
< |
phi = sin(phi0) > 0 ? phi : -phi; |
| 301 |
< |
|
| 224 |
> |
Xlo - x; |
| 225 |
> |
double phi = static_cast<double>((signed int)j2) * YbinSize + |
| 226 |
> |
Ylo - y; |
| 227 |
> |
|
| 228 |
|
// transform eta, phi to proper angle |
| 229 |
|
double theta = 2.0*atan(exp(-eta)); |
| 230 |
|
double iota = asin(sin(theta)*sin(phi)); |
| 231 |
|
|
| 232 |
|
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 233 |
|
* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_)); |
| 308 |
– |
|
| 309 |
– |
// correct for out-of-range phi |
| 310 |
– |
double transform = 0.0; |
| 311 |
– |
if (histo->GetYaxis()->GetXmin() < -PI) { |
| 312 |
– |
transform = -2.0*PI; |
| 313 |
– |
phi += transform; |
| 314 |
– |
iota = asin(sin(theta)*sin(phi)); |
| 315 |
– |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 316 |
– |
* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_)); |
| 317 |
– |
} |
| 318 |
– |
if (histo->GetYaxis()->GetXmax() > PI) { |
| 319 |
– |
transform = 2.0*PI; |
| 320 |
– |
phi += transform; |
| 321 |
– |
iota = asin(sin(theta)*sin(phi)); |
| 322 |
– |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 323 |
– |
* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_)); |
| 324 |
– |
} |
| 234 |
|
} |
| 235 |
|
} |
| 236 |
|
} |
| 241 |
|
double N = particles[i]->energy(); |
| 242 |
|
double x = particles[i]->eta(); |
| 243 |
|
double y = particles[i]->phi(); |
| 244 |
+ |
if (y < Ylo) { |
| 245 |
+ |
y += 2.0*PI; |
| 246 |
+ |
} |
| 247 |
+ |
if (y > Yhi) { |
| 248 |
+ |
y -= 2.0*PI; |
| 249 |
+ |
} |
| 250 |
|
// loop over bins and add Gaussian to smeared |
| 251 |
< |
for (vector< vector<double> >::size_type i2 = 0; i2 < 60; i2++) { |
| 252 |
< |
for (vector< double >::size_type j2 = 0; j2 < 60; j2++) { |
| 251 |
> |
for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) { |
| 252 |
> |
for (vector< double >::size_type j2 = 0; j2 < 30; j2++) { |
| 253 |
|
double eta = static_cast<double>((signed int)i2) * XbinSize |
| 254 |
< |
+ max_x - rcone - x; |
| 255 |
< |
double phi0 = static_cast<double>((signed int)j2) * YbinSize + |
| 256 |
< |
max_y - rcone - y; |
| 342 |
< |
double phi = acos(cos(phi0)); |
| 343 |
< |
phi = sin(phi0) > 0 ? phi : -phi; |
| 344 |
< |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 345 |
< |
* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_)); |
| 254 |
> |
+ Xlo - x; |
| 255 |
> |
double phi = static_cast<double>((signed int)j2) * YbinSize + |
| 256 |
> |
Ylo - y; |
| 257 |
|
|
| 258 |
< |
// correct for out-of-range phi |
| 348 |
< |
double transform = 0.0; |
| 349 |
< |
if (histo->GetYaxis()->GetXmin() < -PI) { |
| 350 |
< |
transform = -2.0*PI; |
| 351 |
< |
phi += transform; |
| 352 |
< |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 353 |
< |
* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_)); |
| 354 |
< |
} |
| 355 |
< |
if (histo->GetYaxis()->GetXmax() > PI) { |
| 356 |
< |
transform = 2.0*PI; |
| 357 |
< |
phi += transform; |
| 358 |
< |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 258 |
> |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 259 |
|
* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_)); |
| 360 |
– |
} |
| 260 |
|
} |
| 261 |
|
} |
| 262 |
|
} |
| 263 |
|
} |
| 264 |
|
// set histogram to match smear vector |
| 265 |
< |
for (int i = 1; i <= 60; i++) { |
| 266 |
< |
for (int j = 1; j <= 60; j++) { |
| 265 |
> |
for (int i = 1; i <= 30; i++) { |
| 266 |
> |
for (int j = 1; j <= 30; j++) { |
| 267 |
|
histo->SetBinContent(i, j, smeared[i-1][j-1]); |
| 268 |
|
} |
| 269 |
|
} |
| 271 |
|
return histo; |
| 272 |
|
} |
| 273 |
|
|
| 274 |
< |
void seed_with_CA(vector<reco::Candidate *> gParticles, TH2 *histo, |
| 275 |
< |
jetfit::model_def &_mdef) { |
| 276 |
< |
// create a PseudoJet vector |
| 277 |
< |
vector<PseudoJet> particles; |
| 278 |
< |
for (unsigned i = 0; i < gParticles.size(); i++) { |
| 380 |
< |
double x_max = (histo->GetXaxis()->GetXmax() |
| 381 |
< |
+ histo->GetXaxis()->GetXmin()) / 2.0; |
| 382 |
< |
double y_max = (histo->GetYaxis()->GetXmax() |
| 383 |
< |
+ histo->GetYaxis()->GetXmin()) / 2.0; |
| 384 |
< |
valarray<double> pmom(4); |
| 385 |
< |
pmom[0] = gParticles[i]->px(); |
| 386 |
< |
pmom[1] = gParticles[i]->py(); |
| 387 |
< |
pmom[2] = gParticles[i]->pz(); |
| 388 |
< |
pmom[3] = gParticles[i]->energy(); |
| 389 |
< |
double eta = gParticles[i]->eta(); |
| 390 |
< |
double phi = gParticles[i]->phi(); |
| 391 |
< |
// phi range corrections |
| 392 |
< |
if (phi - y_max > 2.0*PI) phi -= 2.0*PI; |
| 393 |
< |
if (phi - y_max < -2.0*PI) phi += 2.0*PI; |
| 394 |
< |
if ((eta - x_max)*(eta - x_max) + (phi - y_max)*(phi - y_max) < 0.25) { |
| 395 |
< |
PseudoJet j(pmom); |
| 396 |
< |
particles.push_back(j); |
| 397 |
< |
} |
| 398 |
< |
} |
| 399 |
< |
|
| 400 |
< |
// choose a jet definition |
| 401 |
< |
double R = 0.2; |
| 402 |
< |
JetDefinition jet_def(cambridge_algorithm, R); |
| 403 |
< |
|
| 404 |
< |
// run clustering and extract the jets |
| 405 |
< |
ClusterSequence cs(particles, jet_def); |
| 406 |
< |
vector<PseudoJet> jets = cs.inclusive_jets(); |
| 407 |
< |
|
| 408 |
< |
double XbinSize = (histo->GetXaxis()->GetXmax() |
| 409 |
< |
- histo->GetXaxis()->GetXmin()) / |
| 410 |
< |
static_cast<double>(histo->GetXaxis()->GetNbins()); |
| 411 |
< |
double YbinSize = (histo->GetYaxis()->GetXmax() |
| 412 |
< |
- histo->GetYaxis()->GetXmin()) / |
| 413 |
< |
static_cast<double>(histo->GetYaxis()->GetNbins()); |
| 414 |
< |
|
| 415 |
< |
// seed with C-A jets |
| 274 |
> |
template <class Jet> |
| 275 |
> |
void seed_with_jetcoll(vector<Jet> jets, |
| 276 |
> |
jetfit::model_def &_mdef, double phi_lo = -PI, |
| 277 |
> |
double phi_hi = PI) { |
| 278 |
> |
// seed with jet collection |
| 279 |
|
int ijset = 0; |
| 280 |
|
for (unsigned ij = 0; ij < jets.size(); ij++) { |
| 281 |
< |
double N = jets[ij].e(); |
| 282 |
< |
if (N > 50.0) { |
| 283 |
< |
cout << N << endl; |
| 281 |
> |
double N = jets[ij].energy(); |
| 282 |
> |
double eta = jets[ij].eta(); |
| 283 |
> |
double phi = jets[ij].phi(); |
| 284 |
> |
if (N > 10.0) { |
| 285 |
|
_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1, |
| 286 |
|
0.0, 1.0e6); |
| 287 |
< |
_mdef.set_special_par(ijset, 1, jets[ij].eta(), 0.01, |
| 287 |
> |
_mdef.set_special_par(ijset, 1, eta, 0.01, |
| 288 |
|
0.0, 0.0); |
| 289 |
< |
double mdef_phi = jets[ij].phi() > PI ? jets[ij].phi() - 2*PI |
| 290 |
< |
: jets[ij].phi(); |
| 289 |
> |
double mdef_phi = phi > PI ? phi - 2*PI |
| 290 |
> |
: phi; |
| 291 |
> |
if (mdef_phi < phi_lo) { |
| 292 |
> |
mdef_phi += 2.0*PI; |
| 293 |
> |
} |
| 294 |
> |
if (mdef_phi > phi_hi) { |
| 295 |
> |
mdef_phi -= 2.0*PI; |
| 296 |
> |
} |
| 297 |
|
_mdef.set_special_par(ijset, 2, mdef_phi, 0.01, |
| 298 |
|
0.0, 0.0); |
| 299 |
|
_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0); |
| 308 |
|
class jf_model_def : public jetfit::model_def { |
| 309 |
|
public: |
| 310 |
|
virtual double chisquare_error(double E) { |
| 311 |
< |
return E > 2.7 ? 0.55*E - 1.0 : 0.5; |
| 312 |
< |
// study from 09-03-09 |
| 311 |
> |
return exp(-(4.2 + 0.11*E)); |
| 312 |
> |
// study from 09-04-09 |
| 313 |
|
} |
| 314 |
|
}; |
| 315 |
|
|
| 325 |
|
_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0); |
| 326 |
|
_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0); |
| 327 |
|
|
| 328 |
< |
seed_with_CA(get_particles(evt), histo, *_mdef); |
| 328 |
> |
// get jetcoll from event file and select highest e jet |
| 329 |
> |
if (info_type_ == 0) { |
| 330 |
> |
edm::Handle< vector<reco::GenJet> > jet_collection; |
| 331 |
> |
evt.getByLabel(jet_algo_, jet_collection); |
| 332 |
> |
reco::GenJet highest_e_jet; |
| 333 |
> |
bool found_jet = false; |
| 334 |
> |
for (unsigned i = 0; i < jet_collection->size(); i++) { |
| 335 |
> |
if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy()) |
| 336 |
> |
{ |
| 337 |
> |
highest_e_jet = (*jet_collection)[i]; |
| 338 |
> |
} |
| 339 |
> |
} |
| 340 |
> |
vector<reco::GenJet> highest_e_jet_coll; |
| 341 |
> |
highest_e_jet_coll.push_back(highest_e_jet); |
| 342 |
> |
seed_with_jetcoll(highest_e_jet_coll, *_mdef, |
| 343 |
> |
histo->GetYaxis()->GetXmin(), |
| 344 |
> |
histo->GetYaxis()->GetXmax()); |
| 345 |
> |
} |
| 346 |
> |
if (info_type_ == 1) { |
| 347 |
> |
edm::Handle< vector<reco::PFJet> > jet_collection; |
| 348 |
> |
evt.getByLabel(jet_algo_, jet_collection); |
| 349 |
> |
reco::PFJet highest_e_jet; |
| 350 |
> |
bool found_jet = false; |
| 351 |
> |
for (unsigned i = 0; i < jet_collection->size(); i++) { |
| 352 |
> |
if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy()) |
| 353 |
> |
{ |
| 354 |
> |
highest_e_jet = (*jet_collection)[i]; |
| 355 |
> |
} |
| 356 |
> |
} |
| 357 |
> |
vector<reco::PFJet> highest_e_jet_coll; |
| 358 |
> |
highest_e_jet_coll.push_back(highest_e_jet); |
| 359 |
> |
seed_with_jetcoll(highest_e_jet_coll, *_mdef, |
| 360 |
> |
histo->GetYaxis()->GetXmin(), |
| 361 |
> |
histo->GetYaxis()->GetXmax()); |
| 362 |
> |
} |
| 363 |
|
|
| 364 |
|
jetfit::set_model_def(_mdef); |
| 365 |
|
|
| 406 |
|
} |
| 407 |
|
jetfit::set_ngauss(_mdef->get_n_special_par_sets()); |
| 408 |
|
init_fit_histo->SetBinContent(i+1, j+1, |
| 409 |
< |
jetfit::fit_fcn(x, y, pval)); |
| 409 |
> |
jetfit::fit_fcn(x, y, pval) |
| 410 |
> |
* XbinSize * YbinSize); |
| 411 |
|
} |
| 412 |
|
} |
| 413 |
|
|
| 460 |
|
else |
| 461 |
|
out << "Not calculated" << endl; |
| 462 |
|
} |
| 463 |
+ |
else { |
| 464 |
+ |
out << "Error matrix accurate" << endl; |
| 465 |
+ |
} |
| 466 |
|
|
| 467 |
|
return out; |
| 468 |
|
} |
| 471 |
|
std::vector<jetfit::trouble> t, |
| 472 |
|
TH2 *hist_orig) { |
| 473 |
|
ofs << "Histogram "<<hist_orig->GetName() << endl; |
| 474 |
< |
for (int i = 0; i < unique_jets[hist_orig].size(); i++) { |
| 474 |
> |
for (int i = unique_jets[hist_orig].size() - r.chisquare.size(); |
| 475 |
> |
i < unique_jets[hist_orig].size(); i++) { |
| 476 |
> |
int ir = i - unique_jets[hist_orig].size() + r.chisquare.size(); |
| 477 |
|
ofs << "For "<<i+1<<" gaussians: " << endl |
| 478 |
< |
<< t.at(i) << endl |
| 479 |
< |
<< unique_jets[hist_orig][i].size()<<" unique jets found" << endl; |
| 478 |
> |
<< t.at(i) << endl; |
| 479 |
> |
ofs << "chisquare="<<r.chisquare.at(ir) |
| 480 |
> |
<< endl; |
| 481 |
> |
ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl; |
| 482 |
|
for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) { |
| 483 |
|
jet _jet = unique_jets[hist_orig][i][j]; |
| 484 |
|
ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta |
