| 70 |
|
return sum1 * sum2; |
| 71 |
|
} |
| 72 |
|
|
| 73 |
< |
static void jetfinder(TMinuit *gMinuit, TH2 *hist, int ngauss) { |
| 74 |
< |
double dist_sq = numeric_limits<double>::infinity(); |
| 75 |
< |
unique_jets[hist].resize(ngauss); |
| 76 |
< |
int nbinsX = hist->GetXaxis()->GetNbins(); |
| 77 |
< |
int nbinsY = hist->GetYaxis()->GetNbins(); |
| 78 |
< |
double XbinSize = (hist->GetXaxis()->GetXmax() |
| 79 |
< |
- hist->GetXaxis()->GetXmin()) |
| 80 |
< |
/ static_cast<double>(nbinsX); |
| 81 |
< |
double YbinSize = (hist->GetYaxis()->GetXmax() |
| 82 |
< |
- hist->GetYaxis()->GetXmin()) |
| 83 |
< |
/ static_cast<double>(nbinsY); |
| 84 |
< |
for (int i = 0; i < ngauss; i++) { |
| 85 |
< |
double N, mu_x, mu_y, sig, err, lo, hi; |
| 86 |
< |
int iuint; |
| 87 |
< |
TString name; |
| 88 |
< |
gMinuit->mnpout(4*i, name, N, err, lo, hi, iuint); |
| 89 |
< |
gMinuit->mnpout(4*i + 1, name, mu_x, err, lo, hi, iuint); |
| 90 |
< |
gMinuit->mnpout(4*i + 2, name, mu_y, err, lo, hi, iuint); |
| 91 |
< |
gMinuit->mnpout(4*i + 3, name, sig, err, lo, hi, iuint); |
| 92 |
< |
for (int j = 0; j < i; j++) { |
| 93 |
< |
double N2, mu_x2, mu_y2, sig2; |
| 94 |
< |
gMinuit->mnpout(4*j, name, N2, err, lo, hi, iuint); |
| 95 |
< |
gMinuit->mnpout(4*j + 1, name, mu_x2, err, lo, hi, iuint); |
| 96 |
< |
gMinuit->mnpout(4*j + 2, name, mu_y2, err, lo, hi, iuint); |
| 97 |
< |
gMinuit->mnpout(4*j + 3, name, sig2, err, lo, hi, iuint); |
| 98 |
< |
double _dist_sq = (mu_x2 - mu_x)*(mu_x2 - mu_x) |
| 99 |
< |
+ (mu_y2 - mu_y)*(mu_y2 - mu_y); |
| 100 |
< |
if (_dist_sq < dist_sq) |
| 101 |
< |
dist_sq = _dist_sq; |
| 102 |
< |
} |
| 103 |
< |
|
| 104 |
< |
jet j; |
| 105 |
< |
j.energy = N; |
| 106 |
< |
j.eta = mu_x; j.phi = mu_y; |
| 107 |
< |
unique_jets[hist][ngauss-1].push_back(j); |
| 108 |
< |
} |
| 109 |
< |
} |
| 110 |
< |
|
| 111 |
< |
virtual void beginJob(const edm::EventSetup&); |
| 112 |
< |
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 *); |
| 78 |
|
virtual void analyze_results(jetfit::results r, |
| 79 |
|
std::vector<jetfit::trouble> t, TH2 *); |
| 80 |
|
vector<reco::Candidate *> get_particles(const edm::Event&); |
| 119 |
– |
void fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>&, |
| 120 |
– |
const edm::InputTag&, const edm::Event&) const; |
| 121 |
– |
void fetchCandidateCollection(edm::Handle<reco::GenParticleCollection>&, |
| 122 |
– |
const edm::InputTag&, const edm::Event&) const; |
| 81 |
|
|
| 82 |
|
fstream ofs; |
| 83 |
|
edm::InputTag inputTagPFCandidates_; |
| 86 |
|
double smear_; |
| 87 |
|
int smear_coord_; |
| 88 |
|
string jet_algo_; |
| 89 |
+ |
string reclustered_jet_algo_; |
| 90 |
|
}; |
| 91 |
|
|
| 92 |
|
map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > > |
| 109 |
|
// 0 = eta-phi smear |
| 110 |
|
// 1 = proper angle smear |
| 111 |
|
jet_algo_ = pSet.getParameter<string>("jet_algo"); |
| 112 |
< |
set_user_minuit(jetfinder); |
| 154 |
< |
} |
| 155 |
< |
|
| 156 |
< |
void |
| 157 |
< |
JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::PFCandidateCollection>& c, |
| 158 |
< |
const edm::InputTag& tag, |
| 159 |
< |
const edm::Event& iEvent) const { |
| 160 |
< |
|
| 161 |
< |
bool found = iEvent.getByLabel(tag, c); |
| 162 |
< |
|
| 163 |
< |
if(!found ) { |
| 164 |
< |
ostringstream err; |
| 165 |
< |
err<<" cannot get PFCandidates: " |
| 166 |
< |
<<tag<<endl; |
| 167 |
< |
edm::LogError("PFCandidates")<<err.str(); |
| 168 |
< |
throw cms::Exception( "MissingProduct", err.str()); |
| 169 |
< |
} |
| 170 |
< |
|
| 171 |
< |
} |
| 172 |
< |
|
| 173 |
< |
void |
| 174 |
< |
JetFinderAnalyzer::fetchCandidateCollection(edm::Handle<reco::GenParticleCollection>& c, |
| 175 |
< |
const edm::InputTag& tag, |
| 176 |
< |
const edm::Event& iEvent) const { |
| 177 |
< |
|
| 178 |
< |
bool found = iEvent.getByLabel(tag, c); |
| 179 |
< |
|
| 180 |
< |
if(!found ) { |
| 181 |
< |
ostringstream err; |
| 182 |
< |
err<<" cannot get GenParticles: " |
| 183 |
< |
<<tag<<endl; |
| 184 |
< |
edm::LogError("GenParticles")<<err.str(); |
| 185 |
< |
throw cms::Exception( "MissingProduct", err.str()); |
| 186 |
< |
} |
| 187 |
< |
|
| 188 |
< |
} |
| 189 |
< |
|
| 190 |
< |
vector<reco::Candidate *> JetFinderAnalyzer::get_particles(const edm::Event &evt) { |
| 191 |
< |
// fill unreduced histo |
| 192 |
< |
edm::Handle<reco::GenParticleCollection> hRaw; |
| 193 |
< |
edm::Handle<reco::PFCandidateCollection> hPFlow; |
| 194 |
< |
if (info_type_ == 0) { |
| 195 |
< |
fetchCandidateCollection(hRaw, |
| 196 |
< |
inputTagGenParticles_, |
| 197 |
< |
evt); |
| 198 |
< |
} |
| 199 |
< |
if (info_type_ == 1) { |
| 200 |
< |
fetchCandidateCollection(hPFlow, |
| 201 |
< |
inputTagPFCandidates_, |
| 202 |
< |
evt); |
| 203 |
< |
} |
| 204 |
< |
|
| 205 |
< |
vector<reco::Candidate *> particles; |
| 206 |
< |
|
| 207 |
< |
switch (info_type_) { |
| 208 |
< |
case 0: |
| 209 |
< |
for (unsigned i = 0; i < hRaw->size(); i++) { |
| 210 |
< |
if ((*hRaw)[i].status() == 1) { |
| 211 |
< |
particles.push_back((reco::Candidate *)&((*hRaw)[i])); |
| 212 |
< |
} |
| 213 |
< |
} |
| 214 |
< |
break; |
| 215 |
< |
case 1: |
| 216 |
< |
for (unsigned i = 0; i < hPFlow->size(); i++) { |
| 217 |
< |
particles.push_back((reco::Candidate *)&((*hPFlow)[i])); |
| 218 |
< |
} |
| 219 |
< |
break; |
| 220 |
< |
default: |
| 221 |
< |
cerr << "Unknown event type" << endl; // TODO use MessageLogger |
| 222 |
< |
} |
| 223 |
< |
|
| 224 |
< |
return particles; |
| 112 |
> |
set_user_minuit(0); // we don't need this feature--it's deprecated |
| 113 |
|
} |
| 114 |
|
|
| 115 |
|
TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) { |
| 116 |
< |
const reco::Jet * highest_e_jet; |
| 116 |
> |
const reco::Jet * highest_e_jet = 0; |
| 117 |
|
|
| 118 |
|
if (info_type_ == 0) { |
| 119 |
|
edm::Handle< vector<reco::GenJet> > evtJets; |
| 120 |
|
evt.getByLabel(jet_algo_, evtJets); |
| 121 |
|
for (unsigned i = 0; i < evtJets->size(); i++) { |
| 122 |
< |
if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) { |
| 122 |
> |
if (highest_e_jet == 0 |
| 123 |
> |
|| (*evtJets)[i].energy() > highest_e_jet->energy()) { |
| 124 |
|
highest_e_jet = &((*evtJets)[i]); |
| 125 |
|
} |
| 126 |
|
} |
| 129 |
|
edm::Handle< vector<reco::PFJet> > evtJets; |
| 130 |
|
evt.getByLabel(jet_algo_, evtJets); |
| 131 |
|
for (unsigned i = 0; i < evtJets->size(); i++) { |
| 132 |
< |
if (i == 0 || (*evtJets)[i].energy() > highest_e_jet->energy()) { |
| 132 |
> |
if (highest_e_jet == 0 |
| 133 |
> |
|| (*evtJets)[i].energy() > highest_e_jet->energy()) { |
| 134 |
|
highest_e_jet = &((*evtJets)[i]); |
| 135 |
|
} |
| 136 |
|
} |
| 137 |
|
} |
| 248 |
– |
cout << "found highest e jet" << endl; |
| 138 |
|
|
| 139 |
|
if (highest_e_jet == 0) { |
| 140 |
< |
cerr << "No fatjets found!" << endl; exit(1); |
| 140 |
> |
cerr << "No fat jets found!" << endl; |
| 141 |
> |
return 0; |
| 142 |
|
} |
| 143 |
|
|
| 144 |
|
vector<const reco::Candidate *> particles = |
| 145 |
|
highest_e_jet->getJetConstituentsQuick(); |
| 146 |
< |
cout << "found highest e jet" << endl; |
| 257 |
< |
|
| 146 |
> |
|
| 147 |
|
ostringstream oss; |
| 148 |
|
oss << "eta_phi_energy"<<evt.id().event() << flush; |
| 149 |
< |
TH2D *histo = new TH2D(oss.str().c_str(), oss.str().c_str(), |
| 149 |
> |
edm::Service<TFileService> fs; |
| 150 |
> |
TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(), |
| 151 |
|
30, |
| 152 |
|
highest_e_jet->eta()-0.5, |
| 153 |
|
highest_e_jet->eta()+0.5, |
| 155 |
|
highest_e_jet->phi()-0.5, |
| 156 |
|
highest_e_jet->phi()+0.5); |
| 157 |
|
|
| 158 |
< |
for (int i = 0; i < particles.size(); i++) { |
| 269 |
< |
histo->Fill(particles[i]->eta(), |
| 270 |
< |
particles[i]->phi(), |
| 271 |
< |
particles[i]->energy()); |
| 272 |
< |
} |
| 273 |
< |
|
| 158 |
> |
if (smear_ > 0.0) { |
| 159 |
|
// create a smeared histo |
| 160 |
|
// create a temporary 2D vector for smeared energies |
| 161 |
|
double XbinSize = (histo->GetXaxis()->GetXmax() |
| 165 |
|
- histo->GetYaxis()->GetXmin()) / |
| 166 |
|
static_cast<double>(histo->GetYaxis()->GetNbins()); |
| 167 |
|
double Xlo = histo->GetXaxis()->GetXmin(); |
| 168 |
+ |
double Xhi = histo->GetXaxis()->GetXmax(); |
| 169 |
|
double Ylo = histo->GetYaxis()->GetXmin(); |
| 170 |
+ |
double Yhi = histo->GetYaxis()->GetXmax(); |
| 171 |
|
vector< vector<double> > smeared(30, vector<double>(30, 0.0) ); |
| 172 |
|
switch (smear_coord_) { |
| 173 |
|
case 1: |
| 175 |
|
double N = particles[i]->energy(); |
| 176 |
|
double x = particles[i]->eta(); |
| 177 |
|
double y = particles[i]->phi(); |
| 178 |
+ |
if (y < Ylo) { |
| 179 |
+ |
y += 2.0*PI; |
| 180 |
+ |
} |
| 181 |
+ |
if (y > Yhi) { |
| 182 |
+ |
y -= 2.0*PI; |
| 183 |
+ |
} |
| 184 |
|
// loop over bins and add Gaussian in proper angle to smeared |
| 185 |
|
for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) { |
| 186 |
|
for (vector< double >::size_type j2 = 0; j2 < 30; j2++) { |
| 187 |
|
double eta = static_cast<double>((signed int)i2) * XbinSize + |
| 188 |
< |
Xlo; |
| 189 |
< |
double phi0 = static_cast<double>((signed int)j2) * YbinSize + |
| 190 |
< |
Ylo; |
| 191 |
< |
double phi = acos(cos(phi0)); |
| 299 |
< |
phi = sin(phi0) > 0 ? phi : -phi; |
| 300 |
< |
|
| 188 |
> |
Xlo - x; |
| 189 |
> |
double phi = static_cast<double>((signed int)j2) * YbinSize + |
| 190 |
> |
Ylo - y; |
| 191 |
> |
|
| 192 |
|
// transform eta, phi to proper angle |
| 193 |
|
double theta = 2.0*atan(exp(-eta)); |
| 194 |
|
double iota = asin(sin(theta)*sin(phi)); |
| 195 |
|
|
| 196 |
|
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 197 |
|
* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_)); |
| 307 |
– |
|
| 308 |
– |
// correct for out-of-range phi |
| 309 |
– |
double transform = 0.0; |
| 310 |
– |
if (histo->GetYaxis()->GetXmin() < -PI) { |
| 311 |
– |
transform = -2.0*PI; |
| 312 |
– |
phi += transform; |
| 313 |
– |
iota = asin(sin(theta)*sin(phi)); |
| 314 |
– |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 315 |
– |
* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_)); |
| 316 |
– |
} |
| 317 |
– |
if (histo->GetYaxis()->GetXmax() > PI) { |
| 318 |
– |
transform = 2.0*PI; |
| 319 |
– |
phi += transform; |
| 320 |
– |
iota = asin(sin(theta)*sin(phi)); |
| 321 |
– |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 322 |
– |
* exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_)); |
| 323 |
– |
} |
| 198 |
|
} |
| 199 |
|
} |
| 200 |
|
} |
| 205 |
|
double N = particles[i]->energy(); |
| 206 |
|
double x = particles[i]->eta(); |
| 207 |
|
double y = particles[i]->phi(); |
| 208 |
+ |
if (y < Ylo) { |
| 209 |
+ |
y += 2.0*PI; |
| 210 |
+ |
} |
| 211 |
+ |
if (y > Yhi) { |
| 212 |
+ |
y -= 2.0*PI; |
| 213 |
+ |
} |
| 214 |
|
// loop over bins and add Gaussian to smeared |
| 215 |
|
for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) { |
| 216 |
|
for (vector< double >::size_type j2 = 0; j2 < 30; j2++) { |
| 217 |
|
double eta = static_cast<double>((signed int)i2) * XbinSize |
| 218 |
< |
+ Xlo; |
| 219 |
< |
double phi0 = static_cast<double>((signed int)j2) * YbinSize + |
| 220 |
< |
Ylo; |
| 341 |
< |
double phi = acos(cos(phi0)); |
| 342 |
< |
phi = sin(phi0) > 0 ? phi : -phi; |
| 343 |
< |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 344 |
< |
* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_)); |
| 218 |
> |
+ Xlo - x; |
| 219 |
> |
double phi = static_cast<double>((signed int)j2) * YbinSize + |
| 220 |
> |
Ylo - y; |
| 221 |
|
|
| 222 |
< |
// correct for out-of-range phi |
| 347 |
< |
double transform = 0.0; |
| 348 |
< |
if (histo->GetYaxis()->GetXmin() < -PI) { |
| 349 |
< |
transform = -2.0*PI; |
| 350 |
< |
phi += transform; |
| 351 |
< |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 352 |
< |
* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_)); |
| 353 |
< |
} |
| 354 |
< |
if (histo->GetYaxis()->GetXmax() > PI) { |
| 355 |
< |
transform = 2.0*PI; |
| 356 |
< |
phi += transform; |
| 357 |
< |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 222 |
> |
smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_)) |
| 223 |
|
* exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_)); |
| 359 |
– |
} |
| 224 |
|
} |
| 225 |
|
} |
| 226 |
|
} |
| 231 |
|
histo->SetBinContent(i, j, smeared[i-1][j-1]); |
| 232 |
|
} |
| 233 |
|
} |
| 234 |
+ |
} |
| 235 |
+ |
else { |
| 236 |
+ |
// don't smear--just fill with particles |
| 237 |
+ |
for (int i = 0; i < particles.size(); i++) { |
| 238 |
+ |
histo->Fill(particles[i]->eta(), particles[i]->phi(), |
| 239 |
+ |
particles[i]->energy()); |
| 240 |
+ |
} |
| 241 |
+ |
} |
| 242 |
|
|
| 243 |
|
return histo; |
| 244 |
|
} |
| 245 |
|
|
| 246 |
|
template <class Jet> |
| 247 |
|
void seed_with_jetcoll(vector<Jet> jets, |
| 248 |
< |
jetfit::model_def &_mdef) { |
| 248 |
> |
jetfit::model_def &_mdef, double phi_lo = -PI, |
| 249 |
> |
double phi_hi = PI) { |
| 250 |
|
// seed with jet collection |
| 251 |
|
int ijset = 0; |
| 252 |
|
for (unsigned ij = 0; ij < jets.size(); ij++) { |
| 258 |
|
0.0, 1.0e6); |
| 259 |
|
_mdef.set_special_par(ijset, 1, eta, 0.01, |
| 260 |
|
0.0, 0.0); |
| 261 |
< |
double mdef_phi = phi > PI ? phi - 2*PI |
| 262 |
< |
: phi; |
| 263 |
< |
_mdef.set_special_par(ijset, 2, mdef_phi, 0.01, |
| 261 |
> |
if (phi < phi_lo) { |
| 262 |
> |
phi += 2.0*PI; |
| 263 |
> |
} |
| 264 |
> |
if (phi > phi_hi) { |
| 265 |
> |
phi -= 2.0*PI; |
| 266 |
> |
} |
| 267 |
> |
_mdef.set_special_par(ijset, 2, phi, 0.01, |
| 268 |
|
0.0, 0.0); |
| 269 |
|
_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0); |
| 270 |
|
ijset++; |
| 305 |
|
if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy()) |
| 306 |
|
{ |
| 307 |
|
highest_e_jet = (*jet_collection)[i]; |
| 308 |
+ |
found_jet = true; |
| 309 |
|
} |
| 310 |
|
} |
| 311 |
|
vector<reco::GenJet> highest_e_jet_coll; |
| 312 |
|
highest_e_jet_coll.push_back(highest_e_jet); |
| 313 |
< |
seed_with_jetcoll(highest_e_jet_coll, *_mdef); |
| 313 |
> |
seed_with_jetcoll(highest_e_jet_coll, *_mdef, |
| 314 |
> |
histo->GetYaxis()->GetXmin(), |
| 315 |
> |
histo->GetYaxis()->GetXmax()); |
| 316 |
|
} |
| 317 |
|
if (info_type_ == 1) { |
| 318 |
|
edm::Handle< vector<reco::PFJet> > jet_collection; |
| 323 |
|
if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy()) |
| 324 |
|
{ |
| 325 |
|
highest_e_jet = (*jet_collection)[i]; |
| 326 |
+ |
found_jet = true; |
| 327 |
|
} |
| 328 |
|
} |
| 329 |
|
vector<reco::PFJet> highest_e_jet_coll; |
| 330 |
|
highest_e_jet_coll.push_back(highest_e_jet); |
| 331 |
< |
seed_with_jetcoll(highest_e_jet_coll, *_mdef); |
| 331 |
> |
seed_with_jetcoll(highest_e_jet_coll, *_mdef, |
| 332 |
> |
histo->GetYaxis()->GetXmin(), |
| 333 |
> |
histo->GetYaxis()->GetXmax()); |
| 334 |
|
} |
| 335 |
|
|
| 336 |
|
jetfit::set_model_def(_mdef); |
| 465 |
|
for (vector< vector<double> >::size_type i = 0; |
| 466 |
|
i < r.pval.size(); i++) { |
| 467 |
|
jetfit::set_ngauss(r.pval[i].size() / 4); |
| 585 |
– |
TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2, |
| 586 |
– |
hist_orig->GetXaxis()->GetXmin(), |
| 587 |
– |
hist_orig->GetXaxis()->GetXmax(), |
| 588 |
– |
hist_orig->GetYaxis()->GetXmin(), |
| 589 |
– |
hist_orig->GetYaxis()->GetXmax(), |
| 590 |
– |
r.pval[i].size()); |
| 591 |
– |
for (vector<double>::size_type j = 0; j < r.pval[i].size(); j++) { |
| 592 |
– |
tf2->SetParameter(j, r.pval[i][j]); |
| 593 |
– |
} |
| 468 |
|
ostringstream fit_histo_oss; |
| 469 |
|
fit_histo_oss << hist_orig->GetName()<<"_fit_"<<i << flush; |
| 596 |
– |
tf2->SetNpx(hist_orig->GetXaxis()->GetNbins()); |
| 597 |
– |
tf2->SetNpy(hist_orig->GetYaxis()->GetNbins()); |
| 470 |
|
TH2D *fit_histo = fs->make<TH2D>(fit_histo_oss.str().c_str(), |
| 471 |
|
fit_histo_oss.str().c_str(), |
| 472 |
|
hist_orig->GetXaxis()->GetNbins(), |
| 475 |
|
hist_orig->GetYaxis()->GetNbins(), |
| 476 |
|
hist_orig->GetYaxis()->GetXmin(), |
| 477 |
|
hist_orig->GetYaxis()->GetXmax()); |
| 478 |
< |
TH1 *tf2_histo = tf2->CreateHistogram(); |
| 479 |
< |
double XbinSize = (fit_histo->GetXaxis()->GetXmax() |
| 480 |
< |
- fit_histo->GetXaxis()->GetXmin()) |
| 481 |
< |
/ static_cast<double>(fit_histo->GetXaxis()->GetNbins()); |
| 482 |
< |
double YbinSize = (fit_histo->GetYaxis()->GetXmax() |
| 483 |
< |
- fit_histo->GetYaxis()->GetXmin()) |
| 484 |
< |
/ static_cast<double>(fit_histo->GetYaxis()->GetNbins()); |
| 485 |
< |
for (int ih = 0; ih < tf2->GetNpx(); ih++) { |
| 486 |
< |
for (int jh = 0; jh < tf2->GetNpy(); jh++) { |
| 487 |
< |
fit_histo->SetBinContent(ih+1, jh+1, |
| 488 |
< |
tf2_histo->GetBinContent(ih+1, jh+1) |
| 489 |
< |
* XbinSize * YbinSize); |
| 478 |
> |
double XbinSize = (hist_orig->GetXaxis()->GetXmax() |
| 479 |
> |
- hist_orig->GetXaxis()->GetXmin()) |
| 480 |
> |
/ static_cast<double>(hist_orig->GetXaxis()->GetNbins()); |
| 481 |
> |
double YbinSize = (hist_orig->GetYaxis()->GetXmax() |
| 482 |
> |
- hist_orig->GetYaxis()->GetXmin()) |
| 483 |
> |
/ static_cast<double>(hist_orig->GetYaxis()->GetNbins()); |
| 484 |
> |
|
| 485 |
> |
double Xlo = hist_orig->GetXaxis()->GetXmin(); |
| 486 |
> |
double Xhi = hist_orig->GetXaxis()->GetXmax(); |
| 487 |
> |
double Ylo = hist_orig->GetYaxis()->GetXmin(); |
| 488 |
> |
double Yhi = hist_orig->GetYaxis()->GetXmax(); |
| 489 |
> |
for (int bi = 0; bi < 30; bi++) { |
| 490 |
> |
for (int bj = 0; bj < 30; bj++) { |
| 491 |
> |
double x = (static_cast<double>(bi) + 0.5)*XbinSize + Xlo; |
| 492 |
> |
double y = (static_cast<double>(bj) + 0.5)*YbinSize + Ylo; |
| 493 |
> |
double pval[256]; |
| 494 |
> |
for (int ii = 0; ii < r.pval[i].size(); ii++) { |
| 495 |
> |
pval[ii] = r.pval[i][ii]; |
| 496 |
> |
} |
| 497 |
> |
fit_histo->SetBinContent(bi+1, bj+1, |
| 498 |
> |
jetfit::fit_fcn(x, y, pval) |
| 499 |
> |
* XbinSize * YbinSize); |
| 500 |
|
} |
| 501 |
|
} |
| 502 |
|
} |
