DAS::TPS Namespace Reference

Classes
struct	Hist
struct	MultiJets

Functions
void	CheckConsistency (DT::MetaInfo &metainfo1, TTree *tIn1, DT::MetaInfo &metainfo2, TTree *tIn2)
template<class Jet >
bool	overlappingJets (const vector< Jet > &jets1, const vector< Jet > &jets2)
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})

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");
 
    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;
 
    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) {
            if (overlappingJets(*recJets1, *recJets2)) continue;
 
            // 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
        }
        if (isMC && genJets1 && genJets2) {
            if (overlappingJets(*genJets1, *genJets2)) continue;
 
            // 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
        }
 
        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) );
}

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 > &	jets1,
		const vector< Jet > &	jets2
	)

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

Todo:

Implement the anti-kt clustering algorithm.

                                                                          {
    for (const auto& jet1 : jets1) {
        for (const auto& jet2 : jets2) {
            if (DeltaR(jet1.p4, jet2.p4) < 0.4) return true;
        }
    }
    return false;
}

PhaseSpaceSelection()

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

return Jets pass the phase space selection

Todo:

Check pseudorapidty coverage

                                                                            {
    std::vector<Jet> tempjets;
    for (const auto& jet : jets) 
        if (jet.CorrPt() > minpt && jet.p4.Eta() < 5.0) 
            tempjets.push_back(jet); 
    return tempjets;
}