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

inputs	input ROOT files (n-tuple)
output	output ROOT file (n-tuple)
steering	parameters obtained from explicit options
slice	number and index of slice

                               {1,0} 
        )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow(steering, inputs);
    auto tIn = flow.GetInputTree(slice);
    auto [fOut, tOut] = flow.GetOutput(output);
 
    DT::MetaInfo metainfo(tOut);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    int R = metainfo.Get<int>("flags", "R");
 
    auto gEv = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
    auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto genJets = isMC ? flow.GetBranchReadOnly <vector<GenJet>>("genJets") : nullptr;
    auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
 
    JMEmatching<>::maxDR = R/10./2.;
    cout << "Radius for matching: " << JMEmatching<>::maxDR << endl;
 
    cout << "Setting variations" << endl;
 
    PtVariation::isMC = isMC;
    vector<PtVariation> variations { PtVariation("nominal") };
    if (steering & DT::syst) {
 
        auto push = [&variations](auto... args) {
            variations.emplace_back(args...);
        };
 
        auto getvars = [&push,&metainfo](const char * category, auto indices) {
            TList * ev_wgts = metainfo.List("variations", category);
            for (TObject * obj: *ev_wgts) {
                auto param = dynamic_cast<TParameter<int>*>(obj);
                const int val = param->GetVal();
                static int i = 0;
                // in case of several members (e.g. NNPDF replicas)
                for (int j = 0; j < val; ++j, ++i) {
                    TString varname = param->GetName();
                    if (val > 1) varname += j;
                    apply(push, indices(varname, i));
                }
            }
        };
                                                                                  // JEC GJW RJW GEW REW
        getvars(GenEvent::WeightVar, [](const char * v, int i) { return make_tuple(v, 0 , 0 , 0 , i , 0); });
        getvars(RecEvent::WeightVar, [](const char * v, int i) { return make_tuple(v, 0 , 0 , 0 , 0 , i); });
        getvars(RecJet  :: ScaleVar, [](const char * v, int i) { return make_tuple(v, i , 0 , 0 , 0 , 0); });
    }
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        for (const PtVariation& v: variations) {
 
            v.tmp->Reset();
 
            auto evWgt = rEv->weights.at(v.iRecEvtWgt);
            if (isMC) evWgt *= gEv->weights.at(v.iGenEvtWgt);
 
            vector<int> binIDs; // we save the indices of the filled bins to avoid trivial operations
            for (const RecJet& recjet: *recJets) {
                auto y  = recjet.AbsRap();
                if (y >= maxy) continue;
                auto pt = recjet.CorrPt(v.iJEC);
                if (pt < minpt || pt >= maxpt) continue;
 
                auto irecbin = v.rec->FindBin(pt);
                if (find(binIDs.begin(), binIDs.end(), irecbin) == binIDs.end()) binIDs.push_back(irecbin);
 
                auto jWgt = recjet.weights.at(v.iRecJetWgt);
 
                v.rec->Fill(pt, evWgt * jWgt);
                v.tmp->Fill(pt, evWgt * jWgt); // for cov
            }
 
            for (auto x: binIDs) for (auto y: binIDs) {
                double cCov = v.cov->GetBinContent(x,y),
                       cTmp = v.tmp->GetBinContent(x)*v.tmp->GetBinContent(y);
                v.cov->SetBinContent(x,y,cCov+cTmp);
            }
        }
 
        if (!isMC) continue;
 
        JMEmatching matching(*recJets, *genJets);
        for (const auto& [rec_it,gen_it]: matching.match_its) {
 
            auto genpt = gen_it->p4.Pt(),
                 geny = gen_it->AbsRap();
 
            bool LowGenPt  = genpt < minpt,
                 HighGenPt = genpt >= maxpt,
                 HighGenY  =  geny >= maxy;
            bool goodGen = (!LowGenPt) && (!HighGenPt) && (!HighGenY);
 
            for (const PtVariation& v: variations) {
                auto gEvW = gEv->weights.at(v.iGenEvtWgt),
                     gJetW = gen_it->weights.at(v.iGenJetWgt);
 
                if (goodGen) v.gen->Fill(genpt, gEvW * gJetW);
 
                auto rEvW = rEv->weights.at(v.iRecEvtWgt),
                     rJetW = rec_it->weights.at(v.iRecJetWgt);
 
                auto recpt = rec_it->CorrPt(v.iJEC),
                     recy  = abs(rec_it->Rapidity());
 
                bool LowRecPt  = recpt < minpt,
                     HighRecPt = recpt >= maxpt,
                     HighRecY  =  recy >= maxy;
                bool goodRec = (!LowRecPt) && (!HighRecPt) && (!HighRecY);
 
                if ( goodRec &&  goodGen) {    v.RM->Fill(genpt, recpt, gEvW * gJetW *      rEvW * rJetW  );
                                      v.missNoMatch->Fill(genpt,        gEvW * gJetW * (1 - rEvW * rJetW) ); }
           else if (!goodRec &&  goodGen) v.missOut->Fill(genpt,        gEvW * gJetW                      );
           else if ( goodRec && !goodGen) v.fakeOut->Fill(       recpt, gEvW * gJetW *      rEvW * rJetW  );
            }
        } // end of loop on matched jets
 
        for (const auto& gen_it: matching.miss_its) {
            for (const PtVariation& v: variations) {
 
                double y = gen_it->AbsRap();
                if (y >= maxy) continue;
                double pt = gen_it->p4.Pt();
                if (pt < minpt || pt > maxpt) continue;
 
                auto evWgt = gEv->weights.at(v.iGenEvtWgt);
                auto jWgt = gen_it->weights.at(v.iGenJetWgt);
 
                v.gen        ->Fill(pt, evWgt * jWgt);
                v.missNoMatch->Fill(pt, evWgt * jWgt);
            }
        } // end of loop on missed jets
 
        for (const auto& rec_it: matching.fake_its) {
            for (const auto& v: variations) {
 
                double y = rec_it->AbsRap();
                if (y >= maxy) continue;
                double pt = rec_it->CorrPt(v.iJEC);
                if (pt < minpt || pt > maxpt) continue;
 
                auto evWgt = rEv->weights.at(v.iRecEvtWgt) * gEv->weights.at(v.iGenEvtWgt);
                auto jWgt = rec_it->weights.at(v.iRecJetWgt);
 
                v.fakeNoMatch->Fill(pt, evWgt * jWgt);
            }
        } // end of loop on fake jets
    }
 
    fOut->cd();
    for (PtVariation& v: variations)
        v.Write(fOut);
 
    cout << __func__ << ' ' << slice << " stop" << endl;
}

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