Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
169 changes: 89 additions & 80 deletions ALICE3/TableProducer/OTF/onTheFlyTracker.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@
#include <TGeoGlobalMagField.h>
#include <TH1.h>
#include <TH2.h>
#include <TLorentzVector.h>

Check failure on line 81 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
#include <TMCProcess.h>
#include <TMath.h>
#include <TPDGCode.h>
Expand All @@ -102,8 +102,6 @@
#include <utility>
#include <vector>

#include <math.h>

using namespace o2;
using namespace o2::framework;
using std::array;
Expand Down Expand Up @@ -291,9 +289,9 @@
nSiliconHits(nSiliconHitsInput),
nTPCHits(nTPCHitsInput),
trackType(trackTypeInput) {}
const TimeEst& getTimeMUS() const { return timeEst; }
[[nodiscard]] const TimeEst& getTimeMUS() const { return timeEst; }
int64_t mcLabel = -1; ///< MC label of the track
TimeEst timeEst{}; ///< time estimate in ns
TimeEst timeEst; ///< time estimate in ns
bool isDecayDau = false; ///< is a decay daughter
bool isPVContributor = false; ///< track is PV contributor
bool isWeakDecayDau = false; ///< is a weak decay daughter
Expand Down Expand Up @@ -450,7 +448,7 @@
// load LUTs for primaries
for (const auto& entry : globalConfiguration) {
int pdg = 0;
if (entry.first.find("lut") != 0) {
if (!entry.first.starts_with("lut")) {
continue;
}
if (entry.first.find("lutEl") != std::string::npos) {
Expand Down Expand Up @@ -487,8 +485,7 @@
if (filename.empty()) {
LOG(warning) << "No LUT file passed for pdg " << pdg << ", skipping.";
}
bool success = mSmearer[icfg]->loadTable(pdg, filename.c_str());
if (!success) {
if (!mSmearer[icfg]->loadTable(pdg, filename.c_str())) {
LOG(fatal) << "Having issue with loading the LUT " << pdg << " " << filename;
}
}
Expand All @@ -500,7 +497,7 @@
mSmearer[icfg]->skipUnreconstructed(!processUnreconstructedTracks.value);

insertHist(histPath + "hPtGenerated", "hPtGenerated;#it{p}_{T} (GeV/c);Counts", kTH1D, {{axes.axisMomentum}});
insertHist(histPath + "hPhiGenerated", "hPhiGenerated;#phi (rad);Counts", kTH1D, {{100, 0.0f, 2 * M_PI, "#phi (rad)"}});
insertHist(histPath + "hPhiGenerated", "hPhiGenerated;#phi (rad);Counts", kTH1D, {{100, 0.0f, o2::constants::math::TwoPI, "#phi (rad)"}});

insertHist(histPath + "hPtGeneratedEl", "hPtGeneratedEl;Gen #it{p}_{T} (GeV/c);Counts", kTH1D, {{axes.axisMomentum}});
insertHist(histPath + "hPtGeneratedPi", "hPtGeneratedPi;Gen #it{p}_{T} (GeV/c);Counts", kTH1D, {{axes.axisMomentum}});
Expand Down Expand Up @@ -740,7 +737,7 @@
// Cross-check
LOGF(info, "Check field at (0, 0, 0): %.1f kG, nominal: %.1f", static_cast<float>(fieldInstance->GetBz(0, 0, 0)), static_cast<float>(field));
LOGF(info, "Initializing empty material cylinder LUT - could be better in the future");
o2::base::MatLayerCylSet* lut = new o2::base::MatLayerCylSet();
auto* lut = new o2::base::MatLayerCylSet();
lut->addLayer(200, 200.1, 2, 1.0f, 100.0f);
LOGF(info, "MatLayerCylSet::optimizePhiSlices()");
lut->optimizePhiSlices();
Expand Down Expand Up @@ -814,7 +811,7 @@
/// \param xiDecayVertex the address of the xi decay vertex
/// \param laDecayVertex the address of the la decay vertex
template <typename McParticleType>
void decayCascade(const McParticleType& particle, o2::track::TrackParCov track, std::vector<TLorentzVector>& decayDaughters, std::vector<double>& xiDecayVertex, std::vector<double>& laDecayVertex)

Check failure on line 814 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
{
const double uXi = rand.Uniform(0, 1);
const double ctauXi = 4.91; // cm
Expand All @@ -837,12 +834,12 @@
xiDecayVertex.push_back(particle.vz() + rxyzXi * (particle.pz() / particle.p()));

std::vector<double> xiDaughters = {o2::constants::physics::MassLambda, o2::constants::physics::MassPionCharged};
TLorentzVector xi(newPx, newPy, particle.pz(), newE);

Check failure on line 837 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
TGenPhaseSpace xiDecay;
xiDecay.SetDecay(xi, 2, xiDaughters.data());
xiDecay.Generate();
decayDaughters.push_back(*xiDecay.GetDecay(1));
TLorentzVector la = *xiDecay.GetDecay(0);

Check failure on line 842 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.

const double uLa = rand.Uniform(0, 1);
const double ctauLa = 7.845; // cm
Expand All @@ -865,7 +862,7 @@
/// \param decayDaughters the address of resulting daughters
/// \param v0DecayVertex the address of the la decay vertex
template <typename McParticleType>
void decayV0Particle(const McParticleType& particle, std::vector<TLorentzVector>& decayDaughters, std::vector<double>& v0DecayVertex, int pdgCode)

Check failure on line 865 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
{
double u = rand.Uniform(0, 1);
double v0Mass = -1.;
Expand Down Expand Up @@ -899,7 +896,7 @@

const double v0BetaGamma = particle.p() / v0Mass;
const double v0rxyz = (-v0BetaGamma * ctau * std::log(1 - u));
TLorentzVector v0(particle.px(), particle.py(), particle.pz(), particle.e());

Check failure on line 899 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.

v0DecayVertex.push_back(particle.vx() + v0rxyz * (particle.px() / particle.p()));
v0DecayVertex.push_back(particle.vy() + v0rxyz * (particle.py() / particle.p()));
Expand Down Expand Up @@ -974,7 +971,7 @@
o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB);
const std::string histPath = "Configuration_" + std::to_string(icfg) + "/";

std::vector<TLorentzVector> cascadeDecayProducts;

Check failure on line 974 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
std::vector<double> xiDecayVertex, laDecayVertex;
static constexpr int kCascProngs = 3;
std::array<o2::track::TrackParCov, kCascProngs> xiDaughterTrackParCovsPerfect;
Expand Down Expand Up @@ -1004,11 +1001,11 @@
getHist<TH1>(histPath + "hXiBuilding")->Fill(0.0f);
}

o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[0], xiDecayVertex, xiDaughterTrackParCovsPerfect[0], pdgDB);

Check failure on line 1004 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
xiDaughterTrackParCovsPerfect[0].setPID(pdgCodeToPID(PDG_t::kPiMinus));
o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[1], laDecayVertex, xiDaughterTrackParCovsPerfect[1], pdgDB);

Check failure on line 1006 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
xiDaughterTrackParCovsPerfect[1].setPID(pdgCodeToPID(PDG_t::kPiMinus));
o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kProton, cascadeDecayProducts[2], laDecayVertex, xiDaughterTrackParCovsPerfect[2], pdgDB);

Check failure on line 1008 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector]

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
xiDaughterTrackParCovsPerfect[2].setPID(pdgCodeToPID(PDG_t::kProton));
o2::track::TrackParCov perfectCascadeTrack;
o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectCascadeTrack, pdgDB);
Expand Down Expand Up @@ -1070,10 +1067,10 @@
if (TMath::IsNaN(xiDaughterTrackParCovsTracked[i].getZ())) {
isReco[i] = false;
continue;
} else {
getHist<TH1>(histPath + "hXiBuilding")->Fill(4.0f);
histos.fill(HIST("hNaNBookkeeping"), i + 1, 1.0f);
}
getHist<TH1>(histPath + "hXiBuilding")->Fill(4.0f);
histos.fill(HIST("hNaNBookkeeping"), i + 1, 1.0f);

trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus;
// TODO: add flag for whether it's a ghost track or not, currently assuming all are reconstructed tracks if they pass the fast tracker requirements
TrackType trackType = isReco[i] ? TrackType::kRecoCascDaug : TrackType::kGenCascDaug;
Expand Down Expand Up @@ -1298,7 +1295,7 @@
o2::track::TrackParCov trackedCascade;
const o2::track::TrackParCov& trackedBach = xiDaughterTrackParCovsTracked[0];
const int nCascHits = fastTracker[icfg]->FastTrack(perfectCascadeTrack, trackedCascade, dNdEta, xiDecayRadius2D);
reconstructedCascade = (fastTrackerSettings.minSiliconHitsForKinkReco < nCascHits) ? true : false;
reconstructedCascade = fastTrackerSettings.minSiliconHitsForKinkReco < nCascHits;
if (reconstructedCascade) {
std::array<float, 3> pCasc{};
std::array<float, 3> pBach{};
Expand Down Expand Up @@ -1358,10 +1355,9 @@
std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda});
newCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas
float trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus;
if (reconstructedCascade) {
tracksCascadeProngs[kCascProngs + 1] = TrackAlice3{newCascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, false, false, false, 1, thisCascade.foundClusters, TrackType::kRecoCascDaug};
}
tracksCascadeProngs[kCascProngs + 1] = TrackAlice3{newCascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, false, false, false, 1, thisCascade.foundClusters, 0, TrackType::kRecoCascDaug};

tracksCascadeProngs[kCascProngs + 1] = TrackAlice3{newCascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, false, false, false, 1, thisCascade.foundClusters, TrackType::kRecoCascDaug};

fillCascadeTable = true;
} // end fitter OK
} // end cascade found
Expand Down Expand Up @@ -1840,57 +1836,58 @@
/// \param trackParCov track of the particle to compute bremsstrahlung for
void computeBremsstrahlungLoss(const int icfg, const auto& mcParticle, o2::track::TrackParCov& trackParCov)
{
if (brSettings.radiateBR) {
const o2::fastsim::GeometryEntry geoEntry = mGeoContainer.getEntry(icfg);

for (auto const& layerName : geoEntry.getLayerNames()) {
if (layerName.find("global") != std::string::npos) { // Layers with global tag are skipped
continue;
}
if (!brSettings.radiateBR) {
return;
}
const o2::fastsim::GeometryEntry geoEntry = mGeoContainer.getEntry(icfg);

float mass = o2::constants::physics::MassElectron;
for (auto const& layerName : geoEntry.getLayerNames()) {
if (layerName.find("global") != std::string::npos) { // Layers with global tag are skipped
continue;
}

switch (std::abs(mcParticle.pdgCode())) {
case kElectron:
mass = o2::constants::physics::MassElectron;
break;
case kMuonMinus:
mass = o2::constants::physics::MassMuon;
break;
case kPiPlus:
mass = o2::constants::physics::MassPionCharged;
break;
case kKPlus:
mass = o2::constants::physics::MassKaonCharged;
break;
case kProton:
mass = o2::constants::physics::MassProton;
break;
default:
break;
}
float mass = o2::constants::physics::MassElectron;

switch (std::abs(mcParticle.pdgCode())) {
case kElectron:
mass = o2::constants::physics::MassElectron;
break;
case kMuonMinus:
mass = o2::constants::physics::MassMuon;
break;
case kPiPlus:
mass = o2::constants::physics::MassPionCharged;
break;
case kKPlus:
mass = o2::constants::physics::MassKaonCharged;
break;
case kProton:
mass = o2::constants::physics::MassProton;
break;
default:
break;
}

float lambda = brSettings.radiationStrength * mcParticle.e() * geoEntry.getFloatValue(layerName, "x0") / (mass * mass);
ULong64_t nPhotons = gRandom->Poisson(lambda);
float lambda = brSettings.radiationStrength * mcParticle.e() * geoEntry.getFloatValue(layerName, "x0") / (mass * mass);
ULong64_t nPhotons = gRandom->Poisson(lambda);

double initialMomentum = trackParCov.getP();
double initialMomentum = trackParCov.getP();

for (ULong64_t photon = 0; photon < nPhotons; ++photon) {
float radiativeLoss = 1.0f - brSettings.minBREnergyFraction * std::pow(brSettings.maxBREnergyFraction / brSettings.minBREnergyFraction, gRandom->Rndm());
trackParCov.setQ2Pt(trackParCov.getQ2Pt() / radiativeLoss);
}
for (ULong64_t photon = 0; photon < nPhotons; ++photon) {
float radiativeLoss = 1.0f - brSettings.minBREnergyFraction * std::pow(brSettings.maxBREnergyFraction / brSettings.minBREnergyFraction, gRandom->Rndm());
trackParCov.setQ2Pt(trackParCov.getQ2Pt() / radiativeLoss);
}

double afterRadiationMomentum = trackParCov.getP();
double afterRadiationMomentum = trackParCov.getP();

if (brSettings.doBRQA) {
const std::string histPath = "Configuration_" + std::to_string(icfg) + "/";
if (brSettings.doBRQA) {
const std::string histPath = "Configuration_" + std::to_string(icfg) + "/";

getHist<TH1>(histPath + "h1dNBRPhotons")->Fill(static_cast<double>(nPhotons));
getHist<TH1>(histPath + "h1dBREnergyLoss")->Fill((initialMomentum - afterRadiationMomentum) / afterRadiationMomentum);
getHist<TH1>(histPath + "h1dNBRPhotons")->Fill(static_cast<double>(nPhotons));
getHist<TH1>(histPath + "h1dBREnergyLoss")->Fill((initialMomentum - afterRadiationMomentum) / afterRadiationMomentum);

getHist<TH2>(histPath + "h2dBRPtRes")->Fill(trackParCov.getPt(), (trackParCov.getPt() - mcParticle.pt()) / trackParCov.getPt());
getHist<TH2>(histPath + "h2dBRPtResAbs")->Fill(trackParCov.getPt(), trackParCov.getPt() - mcParticle.pt());
}
getHist<TH2>(histPath + "h2dBRPtRes")->Fill(trackParCov.getPt(), (trackParCov.getPt() - mcParticle.pt()) / trackParCov.getPt());
getHist<TH2>(histPath + "h2dBRPtResAbs")->Fill(trackParCov.getPt(), trackParCov.getPt() - mcParticle.pt());
}
}
}
Expand Down Expand Up @@ -1980,28 +1977,44 @@
}
getHist<TH1>(histPath + "hPtGenerated")->Fill(mcParticle.pt());
getHist<TH1>(histPath + "hPhiGenerated")->Fill(mcParticle.phi());
if (std::abs(mcParticle.pdgCode()) == kElectron)
getHist<TH1>(histPath + "hPtGeneratedEl")->Fill(mcParticle.pt());
if (std::abs(mcParticle.pdgCode()) == kPiPlus)
getHist<TH1>(histPath + "hPtGeneratedPi")->Fill(mcParticle.pt());
if (std::abs(mcParticle.pdgCode()) == kKPlus)
getHist<TH1>(histPath + "hPtGeneratedKa")->Fill(mcParticle.pt());
if (std::abs(mcParticle.pdgCode()) == kProton)
getHist<TH1>(histPath + "hPtGeneratedPr")->Fill(mcParticle.pt());
switch (std::abs(mcParticle.pdgCode())) {
case kElectron:
getHist<TH1>(histPath + "hPtGeneratedEl")->Fill(mcParticle.pt());
break;
case kPiPlus:
getHist<TH1>(histPath + "hPtGeneratedPi")->Fill(mcParticle.pt());
break;
case kKPlus:
getHist<TH1>(histPath + "hPtGeneratedKa")->Fill(mcParticle.pt());
break;
case kProton:
getHist<TH1>(histPath + "hPtGeneratedPr")->Fill(mcParticle.pt());
break;
default:
break;
}

if (!reconstructed && !processUnreconstructedTracks) {
continue;
}

getHist<TH1>(histPath + "hPtReconstructed")->Fill(trackParCov.getPt());
if (std::abs(mcParticle.pdgCode()) == kElectron)
getHist<TH1>(histPath + "hPtReconstructedEl")->Fill(trackParCov.getPt());
if (std::abs(mcParticle.pdgCode()) == kPiPlus)
getHist<TH1>(histPath + "hPtReconstructedPi")->Fill(trackParCov.getPt());
if (std::abs(mcParticle.pdgCode()) == kKPlus)
getHist<TH1>(histPath + "hPtReconstructedKa")->Fill(trackParCov.getPt());
if (std::abs(mcParticle.pdgCode()) == kProton)
getHist<TH1>(histPath + "hPtReconstructedPr")->Fill(trackParCov.getPt());
switch (std::abs(mcParticle.pdgCode())) {
case kElectron:
getHist<TH1>(histPath + "hPtReconstructedEl")->Fill(mcParticle.pt());
break;
case kPiPlus:
getHist<TH1>(histPath + "hPtReconstructedPi")->Fill(mcParticle.pt());
break;
case kKPlus:
getHist<TH1>(histPath + "hPtReconstructedKa")->Fill(mcParticle.pt());
break;
case kProton:
getHist<TH1>(histPath + "hPtReconstructedPr")->Fill(mcParticle.pt());
break;
default:
break;
}
}
if (doExtraQA) {
getHist<TH2>(histPath + "h2dPtRes")->Fill(trackParCov.getPt(), (trackParCov.getPt() - mcParticle.pt()) / trackParCov.getPt());
Expand Down Expand Up @@ -2189,11 +2202,7 @@
computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov);
perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta);
if (nTrkHits < fastTrackerSettings.minSiliconHits) {
reconstructed = false;
} else {
reconstructed = true;
}
reconstructed = nTrkHits >= fastTrackerSettings.minSiliconHits;
}

if (!reconstructed && !processUnreconstructedTracks) {
Expand Down
Loading