DAS  3.0
Das Analysis System
DAS::Unfolding::InclusiveJet1D Namespace Reference

Classes

struct  PtVariation
 

Functions

void getUnfHistPt (const vector< fs::path > &inputs, const fs::path &output, const int steering, const DT::Slice slice={1, 0})
 

Variables

static const std::vector< double > recBins {74,84,97,114,133,153,174,196,220,245,272,300,330,362,395,430,468,507,548,592,638,686,737,790,846,905,967,1032,1101,1172,1248,1327,1410,1497,1588,1684,1784,1890,2000,2116,2238,2366,2500,2640,2787,2941,3103,3273,3450,3637,3832}
 
static const std::vector< double > genBins {74 ,97 ,133 ,174 ,220 ,272 ,330 ,395 ,468 ,548 ,638 ,737 ,846 ,967 ,1101 ,1248 ,1410 ,1588 ,1784 ,2000 ,2238 ,2500 ,2787 ,3103 ,3450, 3832}
 
static const int nRecBins = recBins.size()-1
 
static const int nGenBins = genBins.size()-1
 
static const double minpt = recBins.front()
 
static const double maxpt = recBins.back()
 
static const double maxy = 2.0
 

Function Documentation

◆ getUnfHistPt()

void DAS::Unfolding::InclusiveJet1D::getUnfHistPt ( const vector< fs::path > &  inputs,
const fs::path &  output,
const int  steering,
const DT::Slice  slice = {1,0} 
)

Get histograms for unfolding only for pt (in d=1), without using TUnfoldBinning.

Then the output can be unfolded using unfold.

Todo:
try to use the underflow and overflow for matching out of the phase space...
Parameters
inputsinput ROOT files (n-tuple)
outputoutput ROOT file (n-tuple)
steeringparameters obtained from explicit options
slicenumber and index of slice
105  {1,0}
106  )
107 {
108  cout << __func__ << ' ' << slice << " start" << endl;
109 
110  DT::Flow flow(steering);
111  auto tIn = flow.GetInputTree(inputs, slice);
112  auto [fOut, tOut] = flow.GetOutput(output);
113 
114  DT::MetaInfo metainfo(tOut);
115  auto isMC = metainfo.Get<bool>("flags", "isMC");
116  int R = metainfo.Get<int>("flags", "R");
117 
118  auto gEv = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
119  auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
120  auto genJets = isMC ? flow.GetBranchReadOnly <vector<GenJet>>("genJets") : nullptr;
121  auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
122 
123  JMEmatching<>::maxDR = R/10./2.;
124  cout << "Radius for matching: " << JMEmatching<>::maxDR << endl;
125 
126  cout << "Setting variations" << endl;
127 
129  vector<PtVariation> variations { PtVariation("nominal") };
130  if (steering & DT::syst) {
131  TList * ev_wgts = metainfo.List("variations", RecEvent::WeightVar);
132  for (TObject * obj: *ev_wgts) {
133  auto name = dynamic_cast<TObjString*>(obj)->GetString();
134  static int i = 0; ++i; // JEC GJW RJW GEW REW
135  variations.push_back( PtVariation(name, 0, 0, 0, i) );
136  }
137 
138  TList * JECs = metainfo.List("variations", RecJet::ScaleVar);
139  for (TObject * obj: *JECs) {
140  auto name = dynamic_cast<TObjString*>(obj)->GetString();
141  static int i = 0; ++i; // JEC GJW RJW GEW REW
142  variations.push_back( PtVariation(name, i, 0, 0, 0) );
143  }
144  }
145 
146  for (DT::Looper looper(tIn); looper(); ++looper) {
147  [[ maybe_unused ]]
148  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
149 
150  for (const PtVariation& v: variations) {
151 
152  v.tmp->Reset();
153 
154  auto evWgt = rEv->weights.at(v.iRecEvtWgt);
155  if (isMC) evWgt *= gEv->weights.at(v.iGenEvtWgt);
156 
157  vector<int> binIDs; // we save the indices of the filled bins to avoid trivial operations
158  for (const RecJet& recjet: *recJets) {
159  auto y = recjet.AbsRap();
160  if (y >= maxy) continue;
161  auto pt = recjet.CorrPt(v.iJEC);
162  if (pt < minpt || pt >= maxpt) continue;
163 
164  auto irecbin = v.rec->FindBin(pt);
165  if (find(binIDs.begin(), binIDs.end(), irecbin) == binIDs.end()) binIDs.push_back(irecbin);
166 
167  auto jWgt = recjet.weights.at(v.iRecJetWgt);
168 
169  v.rec->Fill(pt, evWgt * jWgt);
170  v.tmp->Fill(pt, evWgt * jWgt); // for cov
171  }
172 
173  for (auto x: binIDs) for (auto y: binIDs) {
174  double cCov = v.cov->GetBinContent(x,y),
175  cTmp = v.tmp->GetBinContent(x)*v.tmp->GetBinContent(y);
176  v.cov->SetBinContent(x,y,cCov+cTmp);
177  }
178  }
179 
180  if (!isMC) continue;
181 
182  JMEmatching matching(*recJets, *genJets);
183  for (const auto& [rec_it,gen_it]: matching.match_its) {
184 
185  auto genpt = gen_it->p4.Pt(),
186  geny = gen_it->AbsRap();
187 
188  bool LowGenPt = genpt < minpt,
189  HighGenPt = genpt >= maxpt,
190  HighGenY = geny >= maxy;
191  bool goodGen = (!LowGenPt) && (!HighGenPt) && (!HighGenY);
192 
193  for (const PtVariation& v: variations) {
194  auto gEvW = gEv->weights.at(v.iGenEvtWgt),
195  gJetW = gen_it->weights.at(v.iGenJetWgt);
196 
197  if (goodGen) v.gen->Fill(genpt, gEvW * gJetW);
198 
199  auto rEvW = rEv->weights.at(v.iRecEvtWgt),
200  rJetW = rec_it->weights.at(v.iRecJetWgt);
201 
202  auto recpt = rec_it->CorrPt(v.iJEC),
203  recy = abs(rec_it->Rapidity());
204 
205  bool LowRecPt = recpt < minpt,
206  HighRecPt = recpt >= maxpt,
207  HighRecY = recy >= maxy;
208  bool goodRec = (!LowRecPt) && (!HighRecPt) && (!HighRecY);
209 
210  if ( goodRec && goodGen) { v.RM->Fill(genpt, recpt, gEvW * gJetW * rEvW * rJetW );
211  v.missNoMatch->Fill(genpt, gEvW * gJetW * (1 - rEvW * rJetW) ); }
212  else if (!goodRec && goodGen) v.missOut->Fill(genpt, gEvW * gJetW );
213  else if ( goodRec && !goodGen) v.fakeOut->Fill( recpt, gEvW * gJetW * rEvW * rJetW );
214  }
216  } // end of loop on matched jets
217 
218  for (const auto& gen_it: matching.miss_its) {
219  for (const PtVariation& v: variations) {
220 
221  double y = gen_it->AbsRap();
222  if (y >= maxy) continue;
223  double pt = gen_it->p4.Pt();
224  if (pt < minpt || pt > maxpt) continue;
225 
226  auto evWgt = gEv->weights.at(v.iGenEvtWgt);
227  auto jWgt = gen_it->weights.at(v.iGenJetWgt);
228 
229  v.gen ->Fill(pt, evWgt * jWgt);
230  v.missNoMatch->Fill(pt, evWgt * jWgt);
231  }
232  } // end of loop on missed jets
233 
234  for (const auto& rec_it: matching.fake_its) {
235  for (const auto& v: variations) {
236 
237  double y = rec_it->AbsRap();
238  if (y >= maxy) continue;
239  double pt = rec_it->CorrPt(v.iJEC);
240  if (pt < minpt || pt > maxpt) continue;
241 
242  auto evWgt = rEv->weights.at(v.iRecEvtWgt) * gEv->weights.at(v.iGenEvtWgt);
243  auto jWgt = rec_it->weights.at(v.iRecJetWgt);
244 
245  v.fakeNoMatch->Fill(pt, evWgt * jWgt);
246  }
247  } // end of loop on fake jets
248  }
249 
250  fOut->cd();
251  for (PtVariation& v: variations)
252  v.Write(fOut);
253 
254  cout << __func__ << ' ' << slice << " stop" << endl;
255 }

Variable Documentation

◆ genBins

const std::vector<double> genBins {74 ,97 ,133 ,174 ,220 ,272 ,330 ,395 ,468 ,548 ,638 ,737 ,846 ,967 ,1101 ,1248 ,1410 ,1588 ,1784 ,2000 ,2238 ,2500 ,2787 ,3103 ,3450, 3832}
static

◆ maxpt

const double maxpt = recBins.back()
static

◆ maxy

const double maxy = 2.0
static

◆ minpt

const double minpt = recBins.front()
static

◆ nGenBins

const int nGenBins = genBins.size()-1
static

◆ nRecBins

const int nRecBins = recBins.size()-1
static

◆ recBins

const std::vector<double> recBins {74,84,97,114,133,153,174,196,220,245,272,300,330,362,395,430,468,507,548,592,638,686,737,790,846,905,967,1032,1101,1172,1248,1327,1410,1497,1588,1684,1784,1890,2000,2116,2238,2366,2500,2640,2787,2941,3103,3273,3450,3637,3832}
static
DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.name
name
Definition: DYToLL_M-50_13TeV_pythia8_cff_GEN_SIM_RECOBEFMIX_DIGI_L1_DIGI2RAW_L1Reco_RECO.py:48
jmarExample.pt
pt
Definition: jmarExample.py:19
Darwin::Tools::Flow
User-friendly handling of input and output n-tuples.
Definition: Flow.h:69
Step::verbose
static bool verbose
Definition: Step.h:40
DAS::Unfolding::InclusiveJet1D::maxpt
static const double maxpt
Definition: getUnfHistPt.cc:39
Darwin::Tools::Looper
Facility to loop over a n-tuple, including parallelisation and printing.
Definition: Looper.h:32
Darwin::Tools::syst
@ syst
activate -s to systematic uncertainties
Definition: Options.h:29
Darwin::Tools::MetaInfo
Generic meta-information for n-tuple (including speficities to Darwin).
Definition: MetaInfo.h:68
DAS::Unfolding::InclusiveJet1D::minpt
static const double minpt
Definition: getUnfHistPt.cc:39
DAS::PhysicsObject::AbsRap
float AbsRap(const Uncertainties::Variation &=Uncertainties::nominal) const final
absolute rapidity
Definition: PhysicsObject.h:58
recjet
DAS::RecJet recjet
Definition: classes.h:15
DAS::PhysicsObject::CorrPt
float CorrPt(size_t i=0) const
corrected transverse momentum
Definition: PhysicsObject.h:55
Ntupliser_cfg.isMC
string isMC
Definition: Ntupliser_cfg.py:59
DAS::PhysicsObject::weights
Weights weights
object weights
Definition: PhysicsObject.h:52
jercExample.inputs
def inputs
Definition: jercExample.py:118
Darwin::Tools::dev_null
static std::ostream dev_null(nullptr)
to redirect output stream to nowhere
DAS::Unfolding::InclusiveJet1D::maxy
static const double maxy
Definition: getUnfHistPt.cc:39