DAS::TPS Namespace Reference

Classes
struct	has_p4
struct	has_p4< T, std::void_t< decltype(std::declval< T >().p4.Eta())> >
struct	Hist
struct	MultiJets

Enumerations
enum	RemovalStrategy { DeltaR, AntikT, None }

Functions
void	CheckConsistency (DT::MetaInfo &metainfo1, TTree *tIn1, DT::MetaInfo &metainfo2, TTree *tIn2)
RemovalStrategy	parseStrategy (const std::string &strategyStr, TTree *tree)
template<class Jet >
vector< fjcore::PseudoJet >	cluster (const vector< Jet > &jets, const double R)
template<class Jet >
bool	overlappingJets (const vector< Jet > &jets)
void	applyEventMixing (const vector< fs::path > &inputs1, const vector< fs::path > &inputs2, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
void	applyNjetSkim (const vector< fs::path > &inputs, const fs::path &output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})
template<class Jet >
std::vector< Jet >	PhaseSpaceSelection (const std::vector< Jet > &jets, float minpt)
void	getTPSobservables (const vector< fs::path > inputs, const fs::path output, const pt::ptree &config, const int steering, const DT::Slice slice={1, 0})

Enumeration Type Documentation

RemovalStrategy

enum RemovalStrategy

strong

Enumerator
DeltaR
AntikT
None

{ DeltaR, AntikT, None };

Function Documentation

applyEventMixing()

void DAS::TPS::applyEventMixing	(	const vector< fs::path > &	inputs1,
		const vector< fs::path > &	inputs2,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Mix the jets from second input to first input to simulate DPS or TPS events. Overlapping jets are properly treated. The input must be taken from applyNjetSkim and must have the same isMC flag and pTcut.

Todo:

fix output logic

Todo:

Merge metainfos

Parameters

inputs1	input ROOT files (n-tuples)
inputs2	input ROOT files (n-tuples)
output	output ROOT file (n-tuple)
config	config handled with `Darwin::Tools::options`
steering	parameters obtained from explicit options
slice	number and index of slice

                                  {1,0} 
            )
{
    cout << __func__ << ' ' << slice << " start" << endl;
 
    DT::Flow flow1(steering, inputs1);
    DT::Flow flow2(steering, inputs2);
    auto tIn1 = flow1.GetInputTree(slice);
    auto tIn2 = flow2.GetInputTree(slice);
    auto tOut2 = flow2.GetOutputTree(output);
    auto tOut1 = flow1.GetOutputTree(output); 
 
    DT::MetaInfo metainfo1(tOut1);
    DT::MetaInfo metainfo2(tOut2);
    metainfo1.Check(config);
    metainfo2.Check(config); 
    auto isMC = metainfo1.Get<bool>("flags", "isMC");
    auto minpt = metainfo1.Get<float>("skims","jets","minpt");
    auto R = metainfo1.Get<int>("flags","R");
 
    CheckConsistency(metainfo1, tIn1, metainfo2, tIn2);
 
    auto recEvent1 = flow1.GetBranchReadWrite<RecEvent>("recEvent");
    auto recEvent2 = flow2.GetBranchReadOnly <RecEvent>("recEvent");
    auto genEvent1 = isMC ? flow1.GetBranchReadWrite<GenEvent>("genEvent") : nullptr;
    auto genEvent2 = isMC ? flow2.GetBranchReadOnly <GenEvent>("genEvent") : nullptr;
    auto recJets1 = flow1.GetBranchReadWrite<vector<RecJet>>("recJets");
    auto recJets2 = flow2.GetBranchReadOnly <vector<RecJet>>("recJets");
    auto genJets1 = isMC ? flow1.GetBranchReadWrite<vector<GenJet>>("genJets") : nullptr;
    auto genJets2 = isMC ? flow2.GetBranchReadOnly <vector<GenJet>>("genJets") : nullptr;
 
    const auto strategyStr = config.get<string>("corrections.event_mixing.overlap");
    const auto strategy = parseStrategy(strategyStr, tIn1);
 
 
    for (DT::Looper looper1(tIn1), looper2(tIn2); looper1() && looper2(); ++looper1, ++looper2) {
        [[maybe_unused]] 
        static auto& cout = (steering & DT::verbose) ? ::cout : DT::dev_null;
 
        if (recJets1 && recJets2) {
            // Merge input2's jets into input1's jets
            recJets1->insert(recJets1->end(), recJets2->begin(), recJets2->end());
            recEvent1->weights.front() *= recEvent2->weights.front();
 
            sort(recJets1->begin(), recJets1->end(), pt_sort); // sort again the jets by pt
 
            switch (strategy) {
                case RemovalStrategy::DeltaR:
                    if (overlappingJets(*recJets1)) continue;
                    break;
                case RemovalStrategy::AntikT:{
                    // reclustering the jets with anti-kt algorithm
                    vector<fjcore::PseudoJet> inclusive_recJets = cluster(*recJets1, R/10.0);
                    // Jet multiplicity may vary due to the reclustering.
                    // We only leave the event with unchanged jet multiplicity.
                    int nJets1 = PhaseSpaceSelection(*recJets1, minpt).size();
                    int nJets2 = PhaseSpaceSelection(inclusive_recJets, minpt).size();
                    if (nJets1 != nJets2) {
                        cout << "Jet multiplicity changed: " << nJets1 << " -> " << nJets2 << endl;
                        cout << "recJets: ";
                        for (const auto& jet : *recJets1) cout << jet.CorrPt() << " ";
                        cout << "\nAnti-kT: ";
                        for (const auto& jet : inclusive_recJets) cout << jet.pt() << " ";
                        cout << endl;
                        continue;
                    }
                    break;
                }
                case RemovalStrategy::None:;
            }
        }
        if (isMC && genJets1 && genJets2) {
            // Merge input2's jets into input1's jets
            genJets1->insert(genJets1->end(), genJets2->begin(), genJets2->end());
            genEvent1->weights.front() *= genEvent2->weights.front();
 
            sort(genJets1->begin(), genJets1->end(), pt_sort); // sort again the jets by pt
 
            switch (strategy) {
                case RemovalStrategy::DeltaR:
                    if (overlappingJets(*genJets1)) continue;
                    break;
                case RemovalStrategy::AntikT:{
                    // reclustering the jets with anti-kt algorithm
                    vector<fjcore::PseudoJet> inclusive_genJets = cluster(*genJets1, R/10.0);
                    // Jet multiplicity may vary due to the reclustering.
                    // We only leave the event with unchanged jet multiplicity.
                    int nJets1 = PhaseSpaceSelection(*genJets1, minpt).size();
                    int nJets2 = PhaseSpaceSelection(inclusive_genJets, minpt).size();
                    if (nJets1 != nJets2) {
                        cout << "Jet multiplicity changed: " << nJets1 << " -> " << nJets2 << endl;
                        cout << "genJets: ";
                        for (const auto& jet : *genJets1) cout << jet.CorrPt() << " ";
                        cout << "\nAnti-kT: ";
                        for (const auto& jet : inclusive_genJets) cout << jet.pt() << " ";
                        cout << endl;
                        continue;
                    }
                    break;
                }
                case RemovalStrategy::None:;
            }
        }
 
        if (steering & DT::fill) tOut1->Fill();
    }
 
    metainfo1.Set<bool>("git", "complete", true);
    cout << __func__ << ' ' << slice << " stop" << endl;
}

applyNjetSkim()

void DAS::TPS::applyNjetSkim	(	const vector< fs::path > &	inputs,
		const fs::path &	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Skims the input tree to keep only events with different jet multiplicity, with configurable pT thresholds for leading Njets.

Note

Drop the selection for genJets, as it will tighter the selection on rec level in MC samples than data. In case we only interested in the rec level in this stage, we can drop the gen level selection.

Parameters

inputs	input ROOT files (n-tuples)
output	output ROOT file (n-tuple)
config	config handled with `Darwin::Tools::options`
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);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
 
    const auto minpt = config.get<float>("skims.jets.minpt");
    metainfo.Set<float>("skims", "jets", "minpt", minpt);
 
    const auto Njets = config.get<int>("skims.jets.N");
    metainfo.Set<int>("skims", "jets", "N", Njets);
    cout << "Skimming events: Njets = " << Njets << "\t minpt > " << minpt << endl;
 
    auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
    auto genJets = isMC ? flow.GetBranchReadOnly<vector<GenJet>>("genJets") : nullptr;
    auto pileup  = flow.GetBranchReadOnly<PileUp>("pileup");
    
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        if (pileup->nVtx > 1) continue; // skip events with more than one vertex
        
        bool passRec = recJets && PhaseSpaceSelection(*recJets, minpt).size() == Njets;
        bool passGen = genJets && PhaseSpaceSelection(*genJets, minpt).size() == Njets;
 
        if ((steering & DT::fill) && passRec) tOut->Fill(); 
    }
 
    metainfo.Set<bool>("git", "complete", true);
    cout << __func__ << ' ' << slice << " stop" << endl;
}

CheckConsistency()

void DAS::TPS::CheckConsistency	(	DT::MetaInfo &	metainfo1,
		TTree *	tIn1,
		DT::MetaInfo &	metainfo2,
		TTree *	tIn2
	)

Check the consistency between two inputs make sure they have the same isMC flag and pTcut. The number of entries in the two input files doesn’t have to be the same, but the number of entries in the output file will be affected by the smaller one.

                                                                                              {
    struct Dataset {
        string name;
        bool isMC;
        int Njets;
        float minpt;
        long long entries;
    } datasets[] = {
        {"Input1", metainfo1.Get<bool>("flags", "isMC"), metainfo1.Get<int>("skims", "jets", "N"), metainfo1.Get<float>("skims", "jets", "minpt"), tIn1 ? tIn1->GetEntries() : 0},
        {"Input2", metainfo2.Get<bool>("flags", "isMC"), metainfo2.Get<int>("skims", "jets", "N"), metainfo2.Get<float>("skims", "jets", "minpt"), tIn2 ? tIn2->GetEntries() : 0}
    };
    int minIndex = (datasets[0].entries < datasets[1].entries) ? 0 : 1;
    cout << "+---------+------+-------+-------+-------------+\n"
         << "| Dataset | isMC | Njets | pTcut |   Entries   |\n"
         << "+---------+------+-------+-------+-------------+\n";
    for (int i = 0; i < 2; i++) {
        cout << "| " << setw(7) << left << datasets[i].name
             << " | " << setw(4) << left << (datasets[i].isMC ? "Yes" : "No")
             << " | " << setw(5) << left << datasets[i].Njets
             << " | " << setw(5) << left << datasets[i].minpt
             << " | " << setw(11 + (i == minIndex ? 9 : 0)) << right 
             << (i == minIndex ? orange + to_string(datasets[i].entries) + def : to_string(datasets[i].entries))
             << " |\n";
    }
    cout << "+---------+------+-------+-------+-------------+\n" << flush;
 
    if (datasets[0].isMC != datasets[1].isMC)
        BOOST_THROW_EXCEPTION( DE::BadInput("The isMC flag is inconsistent between inputs!", tIn1) );
    if (datasets[0].minpt != datasets[1].minpt)
        BOOST_THROW_EXCEPTION( DE::BadInput("The pTcut is inconsistent between inputs!", tIn1) );
}

cluster()

vector<fjcore::PseudoJet> DAS::TPS::cluster	(	const vector< Jet > &	jets,
		const double	R
	)

Anti-kt clustering algorithm to recluster the jets. The jets are sorted by pT and the reclustered jets are returned.

                                                                          {
    vector<fjcore::PseudoJet> input_jets;
    input_jets.reserve(jets.size()+1);
    for (const auto& jet: jets) {
        const auto& p4 = jet.CorrP4();
        input_jets.push_back(fjcore::PseudoJet(p4.px(),p4.py(),p4.pz(),p4.E()));
    }
    fjcore::JetDefinition jet_def(fjcore::antikt_algorithm, R);
    fjcore::ClusterSequence clust_seq(input_jets, jet_def);
    vector<fjcore::PseudoJet> inclusive_jets = sorted_by_pt(clust_seq.inclusive_jets(0.0));
    return inclusive_jets;
}

getTPSobservables()

void DAS::TPS::getTPSobservables	(	const vector< fs::path >	inputs,
		const fs::path	output,
		const pt::ptree &	config,
		const int	steering,
		const DT::Slice	slice = {1,0}
	)

Calculate the observables that can be used to caracterise 6 jets events.

Parameters

inputs	input ROOT files (n-tuples)
output	name of output root file containing the histograms
config	config file from `DTOptions`
steering	steering parameters from `DTOptions`
slice	slices for running

                                   {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);
    metainfo.Check(config);
    auto isMC = metainfo.Get<bool>("flags", "isMC");
    
    const auto minpt = config.get<float>("skims.jets.minpt");
    metainfo.Set<float>("skims", "jets", "minpt", minpt);
 
    auto rEv = flow.GetBranchReadOnly<RecEvent>("recEvent");
    auto recJets = flow.GetBranchReadOnly<vector<RecJet>>("recJets");
    auto gEv = isMC ? flow.GetBranchReadOnly<GenEvent>("genEvent") : nullptr;
    auto genJets = isMC ? flow.GetBranchReadOnly<vector<GenJet>>("genJets") : nullptr;
 
    Hist genHist("gen"),
         recHist("rec");
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        MultiJets recMultijets(*recJets, minpt);
        if (recMultijets) {
            double evweight = rEv->weights.front();
            if (isMC) evweight *= gEv->weights.front();
            recHist.Fill(recMultijets, evweight);
        }
 
        if (!isMC) continue;
        MultiJets genMultijets(*genJets, minpt);
        if (genMultijets) {
            double evweight = gEv->weights.front();
            genHist.Fill(genMultijets, evweight);
        }
    }
 
    recHist.Write(fOut);
    if (isMC)
        genHist.Write(fOut);
    metainfo.Set<bool>("git", "complete", true);
 
    cout << __func__ << ' ' << slice << " end" << endl;
}

overlappingJets()

bool DAS::TPS::overlappingJets

(

const vector< Jet > &

jets

)

Loop all possible jets combinations checking azimuthal difference to avioding overlapping jets.

                                               {
    for (size_t i = 0; i < jets.size(); ++i) {
        for (size_t j = i + 1; j < jets.size(); ++j) {
            if (DeltaR(jets[i].p4, jets[j].p4) < 0.4)
                return true;
        }
    }
    return false;
}

parseStrategy()

RemovalStrategy DAS::TPS::parseStrategy	(	const std::string &	strategyStr,
		TTree *	tree
	)

Define the strategy to remove overlapping jets.

                                                                         {
    cout << "Removal strategy = " << strategyStr << endl;
    if (strategyStr == "DeltaR")      return RemovalStrategy::DeltaR;
    else if (strategyStr == "antikt") return RemovalStrategy::AntikT;
    else if (strategyStr == "none")   return RemovalStrategy::None;
 
    BOOST_THROW_EXCEPTION(DE::BadInput(("Unrecognized strategy: " + strategyStr).c_str(), tree));
}

PhaseSpaceSelection()

std::vector<Jet> DAS::TPS::PhaseSpaceSelection	(	const std::vector< Jet > &	jets,
		float	minpt
	)

                                                                            {
    std::vector<Jet> tempjets;
    for (const auto& jet : jets) {
        double eta;
        if constexpr(has_p4<Jet>::value) {
            eta = jet.p4.Eta();
        } else {
            eta = jet.eta();
        }
        if (jet.CorrPt() > minpt && eta < 5.0)
            tempjets.push_back(jet);
    }
    return tempjets;
}