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 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)
map< Long64_t, vector< Long64_t > >	buildRunToEntriesMap (TTree *tree, RecEvent *recEvent)
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. Only events with matching run numbers are mixed. 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);
 
    // Build the map from run numbers to entry indices and track the last used index for each run in the 2nd tree
    map<Long64_t, vector<Long64_t>> runToEntries = buildRunToEntriesMap(tIn2, recEvent2);
    map<Long64_t, size_t> lastUsedIndex;
    
    Long64_t matchedEvents = 0;
    
    for (DT::Looper looper(tIn1); looper(); ++looper) {
        [[maybe_unused]] 
        static auto& cout = (steering & DT::verbose) ? ::cout : DT::dev_null;
 
        int run1 = recEvent1->runNo;
        // Skip if no matching run number
        if (runToEntries.find(run1) == runToEntries.end())
            continue;
        
        // Get all entries with matching run number
        const auto& matchingEntries = runToEntries[run1];
        // Get the next event to use in the 2nd tree
        size_t currentIndex = lastUsedIndex[run1];
        Long64_t entry2 = matchingEntries[currentIndex];
        tIn2->GetEntry(entry2);
 
        // sanity check
        if (recEvent1->runNo != recEvent2->runNo) 
            BOOST_THROW_EXCEPTION(DE::BadInput(("Run number mismatch: " + to_string(recEvent1->runNo) + " != " + to_string(recEvent2->runNo)).c_str(), tIn1));
 
        // Update index, wrap around if we reach the end
        lastUsedIndex[run1] = (currentIndex + 1) % matchingEntries.size();
 
        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:;
            }
        }
        ++matchedEvents;
        if (steering & DT::fill) tOut1->Fill();
    }
    
    // Calculate final match rate
    cout << "Processed " << tIn1->GetEntries() << " events, matched " << matchedEvents << " events (" 
         << fixed << setprecision(2) << matchedEvents * 100.0 / tIn1->GetEntries() << "%)." << endl;
 
    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;
}

buildRunToEntriesMap()

map<Long64_t, vector<Long64_t> > DAS::TPS::buildRunToEntriesMap	(	TTree *	tree,
		RecEvent *	recEvent
	)

Build a map that groups event entries by run number.

                                                                                      {
    map<Long64_t, vector<Long64_t>> runToEntries;
    cout << "Creating run-to-entries map..." << endl;
    for (DT::Looper looper(tree); looper(); ++looper) {
        runToEntries[recEvent->runNo].push_back(*looper);
    }
    cout << "Found " << runToEntries.size() << " unique run numbers in tree" << endl;
    return runToEntries;
}

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 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;
}