108 cout << __func__ <<
' ' << slice <<
" start" << endl;
112 auto tOut = unique_ptr<TTree>(tIn->CloneTree(0));
115 auto isMC = metainfo.Get<
bool>(
"flags",
"isMC");
116 int R = metainfo.Get<
int>(
"flags",
"R");
119 RecEvent * rEv =
nullptr;
120 GenEvent * gEv =
nullptr;
121 tIn->SetBranchAddress(
"recEvent", &rEv);
123 tIn->SetBranchAddress(
"genEvent", &gEv);
126 vector<RecJet> * recJets =
nullptr;
127 tIn->SetBranchAddress(
"recJets", &recJets);
128 vector<GenJet> * genJets =
nullptr;
130 tIn->SetBranchAddress(
"genJets", &genJets);
132 JMEmatching<>::maxDR = R/10./2.;
133 cout <<
"Radius for matching: " << JMEmatching<>::maxDR << endl;
135 cout <<
"Setting variations" << endl;
138 vector<PtVariation> variations { PtVariation(
"nominal") };
140 TList * ev_wgts = metainfo.List(
"variations", RecEvent::WeightVar);
141 for (TObject * obj: *ev_wgts) {
142 auto name =
dynamic_cast<TObjString*
>(obj)->GetString();
143 static int i = 0; ++i;
144 variations.push_back( PtVariation(
name, 0, 0, 0, i) );
147 TList * JECs = metainfo.List(
"variations", RecJet::ScaleVar);
148 for (TObject * obj: *JECs) {
149 auto name =
dynamic_cast<TObjString*
>(obj)->GetString();
150 static int i = 0; ++i;
151 variations.push_back( PtVariation(
name, i, 0, 0, 0) );
155 for (
DT::Looper looper(tIn, slice); looper(); ++looper) {
159 for (
const PtVariation& v: variations) {
163 auto evWgt = rEv->weights.at(v.iRecEvtWgt);
164 if (
isMC) evWgt *= gEv->weights.at(v.iGenEvtWgt);
167 for (
const RecJet&
recjet: *recJets) {
169 if (y >=
maxy)
continue;
171 if (pt < minpt || pt >=
maxpt)
continue;
173 auto irecbin = v.rec->FindBin(
pt);
174 if (find(binIDs.begin(), binIDs.end(), irecbin) == binIDs.end()) binIDs.push_back(irecbin);
178 v.rec->Fill(
pt, evWgt * jWgt);
179 v.tmp->Fill(
pt, evWgt * jWgt);
182 for (
auto x: binIDs)
for (
auto y: binIDs) {
183 double cCov = v.cov->GetBinContent(x,y),
184 cTmp = v.tmp->GetBinContent(x)*v.tmp->GetBinContent(y);
185 v.cov->SetBinContent(x,y,cCov+cTmp);
191 JMEmatching matching(*recJets, *genJets);
192 for (
const auto& [rec_it,gen_it]: matching.match_its) {
194 auto genpt = gen_it->p4.Pt(),
195 geny = gen_it->AbsRap();
197 bool LowGenPt = genpt <
minpt,
198 HighGenPt = genpt >=
maxpt,
199 HighGenY = geny >=
maxy;
200 bool goodGen = (!LowGenPt) && (!HighGenPt) && (!HighGenY);
202 for (
const PtVariation& v: variations) {
203 auto gEvW = gEv->weights.at(v.iGenEvtWgt),
204 gJetW = gen_it->weights.at(v.iGenJetWgt);
206 if (goodGen) v.gen->Fill(genpt, gEvW * gJetW);
208 auto rEvW = rEv->weights.at(v.iRecEvtWgt),
209 rJetW = rec_it->weights.at(v.iRecJetWgt);
211 auto recpt = rec_it->CorrPt(v.iJEC),
212 recy = abs(rec_it->Rapidity());
214 bool LowRecPt = recpt <
minpt,
215 HighRecPt = recpt >=
maxpt,
216 HighRecY = recy >=
maxy;
217 bool goodRec = (!LowRecPt) && (!HighRecPt) && (!HighRecY);
219 if ( goodRec && goodGen) { v.RM->Fill(genpt, recpt, gEvW * gJetW * rEvW * rJetW );
220 v.missNoMatch->Fill(genpt, gEvW * gJetW * (1 - rEvW * rJetW) ); }
221 else if (!goodRec && goodGen) v.missOut->Fill(genpt, gEvW * gJetW );
222 else if ( goodRec && !goodGen) v.fakeOut->Fill( recpt, gEvW * gJetW * rEvW * rJetW );
227 for (
const auto& gen_it: matching.miss_its) {
228 for (
const PtVariation& v: variations) {
230 double y = gen_it->AbsRap();
231 if (y >=
maxy)
continue;
232 double pt = gen_it->p4.Pt();
233 if (pt < minpt || pt >
maxpt)
continue;
235 auto evWgt = gEv->weights.at(v.iGenEvtWgt);
236 auto jWgt = gen_it->weights.at(v.iGenJetWgt);
238 v.gen ->Fill(
pt, evWgt * jWgt);
239 v.missNoMatch->Fill(
pt, evWgt * jWgt);
243 for (
const auto& rec_it: matching.fake_its) {
244 for (
const auto& v: variations) {
246 double y = rec_it->AbsRap();
247 if (y >=
maxy)
continue;
248 double pt = rec_it->CorrPt(v.iJEC);
249 if (pt < minpt || pt >
maxpt)
continue;
251 auto evWgt = rEv->weights.at(v.iRecEvtWgt) * gEv->weights.at(v.iGenEvtWgt);
252 auto jWgt = rec_it->weights.at(v.iRecJetWgt);
254 v.fakeNoMatch->Fill(
pt, evWgt * jWgt);
260 for (PtVariation& v: variations)
266 cout << __func__ <<
' ' << slice <<
" stop" << endl;