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
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})
float	ShiftJetOrigin (const FourVector &p4_orig, const float &vertexZorig)
template<typename ObjT >
void	ShiftP4EtaAndFillHist (std::vector< ObjT > &objects, float vertexZ, TH2D &hist)
void	applyPVshift (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 R = metainfo1.Get<int>("flags","R");
 
    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 recMuons1 = flow1.GetBranchReadWrite<vector<RecMuon>>("recMuons", DT::facultative);
    auto recMuons2 = flow2.GetBranchReadOnly <vector<RecMuon>>("recMuons", DT::facultative);
    auto genJets1 = isMC ? flow1.GetBranchReadWrite<vector<GenJet>>("genJets") : nullptr;
    auto genJets2 = isMC ? flow2.GetBranchReadOnly <vector<GenJet>>("genJets") : nullptr;
    auto genMuons1 = isMC ? flow1.GetBranchReadWrite<vector<GenMuon>>("genMuons", DT::facultative) : nullptr;
    auto genMuons2 = isMC ? flow2.GetBranchReadOnly <vector<GenMuon>>("genMuons", DT::facultative) : nullptr;
    auto recPVs1 = flow1.GetBranchReadOnly<vector<PrimaryVertex>>("recPVs");
    auto recPVs2 = flow2.GetBranchReadOnly<vector<PrimaryVertex>>("recPVs");
    auto genPVs1 = isMC ? flow1.GetBranchReadOnly<vector<PrimaryVertex>>("genPVs") : nullptr;
    auto genPVs2 = isMC ? flow2.GetBranchReadOnly<vector<PrimaryVertex>>("genPVs") : nullptr;
 
    const auto minpt = config.get<float>("skims.jets.minpt");
    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;
 
        const auto& matchingEntries = runToEntries[run1];
        size_t currentIndex = lastUsedIndex[run1];
        if (currentIndex >= matchingEntries.size()) continue; // no more entries to match
        Long64_t entry2 = matchingEntries[currentIndex];
        tIn2->GetEntry(entry2);
        lastUsedIndex[run1] = currentIndex + 1;
 
        // sanity check
        if (recEvent1->runNo != recEvent2->runNo)
            cerr << orange << "Run number mismatch: " << recEvent1->runNo << " != " << recEvent2->runNo << def << endl;
        if (recPVs1->front().z != recPVs2->front().z)
            cerr << orange << "Unmatched rec-level primary vertex: " << recPVs1->front().z << " != " << recPVs2->front().z << def << endl;
 
        if (recJets1 && recJets2) {
            // Merge input2's jets into input1's jets
            recJets1->insert(recJets1->end(), recJets2->begin(), recJets2->end());
            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 (recMuons1 && recMuons2) {
            // Merge input2's muons into input1's muons
            recMuons1->insert(recMuons1->end(), recMuons2->begin(), recMuons2->end());
            sort(recMuons1->begin(), recMuons1->end(), pt_sort);
        }
 
        if ((recJets1 && recJets2) || (recMuons1 && recMuons2))
             recEvent1->weights.front() *= recEvent2->weights.front();
 
        if (isMC) {
            if (genJets1 && genJets2) {
                // sanity check
                if (genPVs1->front().z != genPVs2->front().z)
                    BOOST_THROW_EXCEPTION(DE::BadInput(("Unmatched gen-level primary vertex: " + to_string(genPVs1->front().z) + " != " + to_string(genPVs2->front().z)).c_str(), tIn1));
                // Merge input2's jets into input1's jets
                genJets1->insert(genJets1->end(), genJets2->begin(), genJets2->end());
 
                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 (genMuons1 && genMuons2) {
                // Merge input2's muons into input1's muons
                genMuons1->insert(genMuons1->end(), genMuons2->begin(), genMuons2->end());
                sort(genMuons1->begin(), genMuons1->end(), pt_sort);
            }
 
            if ((genJets1 && genJets2) || (genMuons1 && genMuons2))
                genEvent1->weights.front() *= genEvent2->weights.front();
        }
        ++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 / min(tIn1->GetEntries(), tIn2->GetEntries()) << "%)." << endl;
 
    metainfo1.Set<bool>("git", "complete", true);
    cout << __func__ << ' ' << slice << " stop" << endl;
}

applyPVshift()

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

Translation of the primary vertex to the center of the detector. The jet four-momentum is adjusted to account for the vertex shift.

Todo:

Shift photons too

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");
 
    auto recJets = flow.GetBranchReadWrite<vector<RecJet>>("recJets");
    auto genJets = isMC ? flow.GetBranchReadWrite<vector<GenJet>>("genJets") : nullptr;
    auto recMuons = flow.GetBranchReadWrite<vector<RecMuon>>("recMuons", DT::facultative);
    auto genMuons = isMC ? flow.GetBranchReadWrite<vector<GenMuon>>("genMuons", DT::facultative) : nullptr;
 
    auto recPVs = flow.GetBranchReadWrite<vector<PrimaryVertex>>("recPVs");
    auto genPVs = isMC ? flow.GetBranchReadWrite<vector<PrimaryVertex>>("genPVs") : nullptr;
 
    TH2D rec_Deta_PVshift_jet("recDeta_PVshift_jet", "#Delta#eta vs V_{z};V_{z} (cm);#Delta#eta", 100, -20, 20, 100, -0.2, 0.2);
    TH2D gen_Deta_PVshift_jet("genDeta_PVshift_jet", "#Delta#eta vs V_{z};V_{z} (cm);#Delta#eta", 100, -20, 20, 100, -0.2, 0.2);
    TH2D rec_Deta_PVshift_muon("recDeta_PVshift_muon", "#Delta#eta vs V_{z} (muons);V_{z} (cm);#Delta#eta", 100, -20, 20, 100, -0.2, 0.2);
    TH2D gen_Deta_PVshift_muon("genDeta_PVshift_muon", "#Delta#eta vs V_{z} (gen muons);V_{z} (cm);#Delta#eta", 100, -20, 20, 100, -0.2, 0.2);
 
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[maybe_unused]]
        static auto& cout = (steering & DT::verbose) ? ::cout : DT::dev_null;
        
        auto& vertexZ = recPVs->front().z;
        ShiftP4EtaAndFillHist(*recJets, vertexZ, rec_Deta_PVshift_jet);
        if (recMuons != nullptr)
            ShiftP4EtaAndFillHist(*recMuons, vertexZ, rec_Deta_PVshift_muon);
        vertexZ = 0;
        if (isMC) {
            auto& vertexZ = genPVs ->front().z;
            ShiftP4EtaAndFillHist(*genJets, vertexZ, gen_Deta_PVshift_jet);
            if (genMuons != nullptr)
                ShiftP4EtaAndFillHist(*genMuons, vertexZ, gen_Deta_PVshift_muon);
            vertexZ = 0;
        }
 
        if (steering & DT::fill) tOut->Fill();
    }
    fOut->cd();
    rec_Deta_PVshift_jet.Write();
    gen_Deta_PVshift_jet.Write();
    rec_Deta_PVshift_muon.Write();
    gen_Deta_PVshift_muon.Write();
    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;
}

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);
        } else {
            break;
        }
    }
    return tempjets;
}

ShiftJetOrigin()

float DAS::TPS::ShiftJetOrigin	(	const FourVector &	p4_orig,
		const float &	vertexZorig
	)

Adjust jet four-momentum to account for a vertex shift

1. Jet direction: $\vec{d} = (\rho=1, \eta_0, \phi_0)$ (invariant)
2. Hit point from original vertex $(0,0,z_v)$:
   • Barrel: $\vec{h}_1 = (R_{cal}\cos\phi, R_{cal}\sin\phi, s \cdot d_z)$ where $s = R_{cal}/\rho_d$
   • Endcap: $\vec{h}_1 = (s\rho_d\cos\phi, s\rho_d\sin\phi, Z_{cal})$ where $s = (Z_{cal})/d_z$
3. Shifted hit point: $\vec{h}_2 = \vec{h}_1 - (0,0,z_v)$ (detector shift)
4. Return: $(p_T, \eta_{new}, \phi_0, m)$ where $\eta_{new} = \vec{h}_2.\text{Eta()}$

                                                                          {
    // CMS calorimeter geometry (units: cm)
    // Only use as reference, not for actual calculations
    const float R_CAL_BARREL = 129.0f;   // ECAL Barrel inner radius
    const float Z_CAL_ENDCAP = 314.0f;   // ECAL Endcap z-position
    const float ETA_TRANSITION = 1.479f; // Barrel/endcap transition
 
    ROOT::Math::RhoEtaPhiVector dir(1.0, p4_orig.Eta(), p4_orig.Phi());
    ROOT::Math::XYZPoint hitPoint;
    
    if (abs(p4_orig.Eta()) < ETA_TRANSITION) {
        // Barrel: hit point on cylindrical surface at R_CAL_BARREL
        float s = R_CAL_BARREL / dir.Rho();
        float hit_z = s * dir.Z();
        hitPoint.SetXYZ(R_CAL_BARREL * cos(p4_orig.Phi()), R_CAL_BARREL * sin(p4_orig.Phi()), hit_z - vertexZorig);
    } else {
        // Endcap: hit point on planar surface at target_z
        float target_z = copysign(Z_CAL_ENDCAP, p4_orig.Eta());
        float s = target_z / dir.Z();
        float hit_rho = s * dir.Rho();
        hitPoint.SetXYZ(hit_rho * cos(p4_orig.Phi()), hit_rho * sin(p4_orig.Phi()), target_z - vertexZorig);
    }
    return hitPoint.Eta();
}

ShiftP4EtaAndFillHist()

void DAS::TPS::ShiftP4EtaAndFillHist	(	std::vector< ObjT > &	objects,
		float	vertexZ,
		TH2D &	hist
	)

Adjusts objects' eta and fills histogram with vertex shift effect.

This function:

1. Loops over objects and gets original eta from four-momentum p4.
2. Computes shifted eta using ShiftJetOrigin with vertexZ.
3. Updates object's p4 with new eta.
4. Fills histogram with (-vertexZ) and eta difference (new_eta - orig_eta).

                                                                                {
    for (auto& obj: objects) {
        float orig_eta = obj.p4.Eta();
        float new_eta = ShiftJetOrigin(obj.p4, vertexZ);
        obj.p4.SetEta(new_eta);
        hist.Fill(-vertexZ, new_eta - orig_eta);
    }
}