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  unique_ptr<TChain> tIn = DT::GetChain(inputs);
111  unique_ptr<TFile> fOut(DT_GetOutput(output));
112  auto tOut = unique_ptr<TTree>(tIn->CloneTree(0));
113 
114  DT::MetaInfo metainfo(tOut);
115  auto isMC = metainfo.Get<bool>("flags", "isMC");
116  int R = metainfo.Get<int>("flags", "R");
117 
118  // event information
119  RecEvent * rEv = nullptr;
120  GenEvent * gEv = nullptr;
121  tIn->SetBranchAddress("recEvent", &rEv);
122  if (isMC)
123  tIn->SetBranchAddress("genEvent", &gEv);
124 
125  // jet information
126  vector<RecJet> * recJets = nullptr;
127  tIn->SetBranchAddress("recJets", &recJets);
128  vector<GenJet> * genJets = nullptr;
129  if (isMC)
130  tIn->SetBranchAddress("genJets", &genJets);
131 
132  JMEmatching<>::maxDR = R/10./2.;
133  cout << "Radius for matching: " << JMEmatching<>::maxDR << endl;
134 
135  cout << "Setting variations" << endl;
136 
138  vector<PtVariation> variations { PtVariation("nominal") };
139  if (steering & DT::syst) {
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; // JEC GJW RJW GEW REW
144  variations.push_back( PtVariation(name, 0, 0, 0, i) );
145  }
146 
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; // JEC GJW RJW GEW REW
151  variations.push_back( PtVariation(name, i, 0, 0, 0) );
152  }
153  }
154 
155  for (DT::Looper looper(tIn, slice); looper(); ++looper) {
156  [[ maybe_unused ]]
157  static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
158 
159  for (const PtVariation& v: variations) {
160 
161  v.tmp->Reset();
162 
163  auto evWgt = rEv->weights.at(v.iRecEvtWgt);
164  if (isMC) evWgt *= gEv->weights.at(v.iGenEvtWgt);
165 
166  vector<int> binIDs; // we save the indices of the filled bins to avoid trivial operations
167  for (const RecJet& recjet: *recJets) {
168  auto y = recjet.AbsRap();
169  if (y >= maxy) continue;
170  auto pt = recjet.CorrPt(v.iJEC);
171  if (pt < minpt || pt >= maxpt) continue;
172 
173  auto irecbin = v.rec->FindBin(pt);
174  if (find(binIDs.begin(), binIDs.end(), irecbin) == binIDs.end()) binIDs.push_back(irecbin);
175 
176  auto jWgt = recjet.weights.at(v.iRecJetWgt);
177 
178  v.rec->Fill(pt, evWgt * jWgt);
179  v.tmp->Fill(pt, evWgt * jWgt); // for cov
180  }
181 
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);
186  }
187  }
188 
189  if (!isMC) continue;
190 
191  JMEmatching matching(*recJets, *genJets);
192  for (const auto& [rec_it,gen_it]: matching.match_its) {
193 
194  auto genpt = gen_it->p4.Pt(),
195  geny = gen_it->AbsRap();
196 
197  bool LowGenPt = genpt < minpt,
198  HighGenPt = genpt >= maxpt,
199  HighGenY = geny >= maxy;
200  bool goodGen = (!LowGenPt) && (!HighGenPt) && (!HighGenY);
201 
202  for (const PtVariation& v: variations) {
203  auto gEvW = gEv->weights.at(v.iGenEvtWgt),
204  gJetW = gen_it->weights.at(v.iGenJetWgt);
205 
206  if (goodGen) v.gen->Fill(genpt, gEvW * gJetW);
207 
208  auto rEvW = rEv->weights.at(v.iRecEvtWgt),
209  rJetW = rec_it->weights.at(v.iRecJetWgt);
210 
211  auto recpt = rec_it->CorrPt(v.iJEC),
212  recy = abs(rec_it->Rapidity());
213 
214  bool LowRecPt = recpt < minpt,
215  HighRecPt = recpt >= maxpt,
216  HighRecY = recy >= maxy;
217  bool goodRec = (!LowRecPt) && (!HighRecPt) && (!HighRecY);
218 
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 );
223  }
225  } // end of loop on matched jets
226 
227  for (const auto& gen_it: matching.miss_its) {
228  for (const PtVariation& v: variations) {
229 
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;
234 
235  auto evWgt = gEv->weights.at(v.iGenEvtWgt);
236  auto jWgt = gen_it->weights.at(v.iGenJetWgt);
237 
238  v.gen ->Fill(pt, evWgt * jWgt);
239  v.missNoMatch->Fill(pt, evWgt * jWgt);
240  }
241  } // end of loop on missed jets
242 
243  for (const auto& rec_it: matching.fake_its) {
244  for (const auto& v: variations) {
245 
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;
250 
251  auto evWgt = rEv->weights.at(v.iRecEvtWgt) * gEv->weights.at(v.iGenEvtWgt);
252  auto jWgt = rec_it->weights.at(v.iRecJetWgt);
253 
254  v.fakeNoMatch->Fill(pt, evWgt * jWgt);
255  }
256  } // end of loop on fake jets
257  }
258 
259  fOut->cd();
260  for (PtVariation& v: variations)
261  v.Write(fOut.get());
262 
263  fOut->cd();
264  tOut->Write();
265 
266  cout << __func__ << ' ' << slice << " stop" << endl;
267 }

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::syst
@ syst
activate -s to systematic uncertainties
Definition: Options.h:30
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:33
Darwin::Tools::GetChain
std::unique_ptr< TChain > GetChain(std::vector< std::filesystem::path > inputs, const char *name="events")
Load chain from a list of files.
Definition: FileUtils.cc:67
Darwin::Tools::MetaInfo
Generic meta-information for n-tuple (including speficities to Darwin).
Definition: MetaInfo.h:65
DAS::Unfolding::InclusiveJet1D::minpt
static const double minpt
Definition: getUnfHistPt.cc:39
DAS::PhysicsObject::AbsRap
float AbsRap(const Uncertainties::Variation &=Uncertainties::nominal) const
absolute rapidity
Definition: PhysicsObject.h:54
recjet
DAS::RecJet recjet
Definition: classes.h:15
DAS::PhysicsObject::CorrPt
float CorrPt(size_t i=0) const
corrected transverse momentum
Definition: PhysicsObject.h:52
Ntupliser_cfg.isMC
string isMC
Definition: Ntupliser_cfg.py:59
DAS::PhysicsObject::weights
Weights weights
object weights
Definition: PhysicsObject.h:49
jercExample.inputs
def inputs
Definition: jercExample.py:118
Darwin::Tools::dev_null
static std::ostream dev_null(nullptr)
to redirect output stream to nowhere
DT_GetOutput
#define DT_GetOutput(output)
Definition: FileUtils.h:96
DAS::Unfolding::InclusiveJet1D::maxy
static const double maxy
Definition: getUnfHistPt.cc:39