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})
float	ShiftJetOrigin (const FourVector &p4_orig, const float &vertexZorig)
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 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;
    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 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());
            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) {
            // 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());
            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 / min(tIn1->GetEntries(), tIn2->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;
    
    for (DT::Looper looper(tIn); looper(); ++looper) {
        [[ maybe_unused ]]
        static auto& cout = steering & DT::verbose ? ::cout : DT::dev_null;
 
        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;
}

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.

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 recPVs = flow.GetBranchReadWrite<vector<PrimaryVertex>>("recPVs");
    auto genPVs = isMC ? flow.GetBranchReadWrite<vector<PrimaryVertex>>("genPVs") : nullptr;
 
    TH2D rec_Deta_PVshift("recDeta_PVshift", "#Delta#eta vs V_{z};V_{z} (cm);#Delta#eta", 100, -20, 20, 100, -0.2, 0.2);
    TH2D gen_Deta_PVshift("genDeta_PVshift", "#Delta#eta vs V_{z};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;
        
        for (auto& jet : *recJets) {
            auto jet_orig_eta = jet.p4.Eta();
            jet.p4.SetEta(ShiftJetOrigin(jet.p4, recPVs->front().z));
            rec_Deta_PVshift.Fill(-recPVs->front().z, (jet.p4.Eta() - jet_orig_eta));
            recPVs->front().z = 0;
        }
        if (isMC) {
            for (auto& jet : *genJets) {
                auto jet_orig_eta = jet.p4.Eta();
                jet.p4.SetEta(ShiftJetOrigin(jet.p4, genPVs->front().z));
                gen_Deta_PVshift.Fill(-genPVs->front().z, (jet.p4.Eta() - jet_orig_eta));
                genPVs->front().z = 0;
            }
        }
 
        if (steering & DT::fill) tOut->Fill();
    }
    fOut->cd();
    rec_Deta_PVshift.Write();
    gen_Deta_PVshift.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;
}

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

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