diff --git a/ALICE3/Core/ALICE3CoreLinkDef.h b/ALICE3/Core/ALICE3CoreLinkDef.h index 05d2a78fc00..076c21f8a79 100644 --- a/ALICE3/Core/ALICE3CoreLinkDef.h +++ b/ALICE3/Core/ALICE3CoreLinkDef.h @@ -14,3 +14,5 @@ #pragma link off all functions; #pragma link C++ class o2::pid::tof::TOFResoALICE3 + ; +#pragma link C++ class std::vector < std::vector < unsigned int>> + ; +#pragma link C++ class std::vector < std::vector < std::uint32_t>> + ; diff --git a/ALICE3/DataModel/tracksAlice3.h b/ALICE3/DataModel/tracksAlice3.h index 280c5ccb110..6e03e5e68ba 100644 --- a/ALICE3/DataModel/tracksAlice3.h +++ b/ALICE3/DataModel/tracksAlice3.h @@ -21,6 +21,7 @@ // O2 includes #include "Framework/AnalysisDataModel.h" +#include namespace o2::aod { @@ -29,16 +30,30 @@ namespace track_alice3 DECLARE_SOA_COLUMN(IsReconstructed, isReconstructed, bool); //! is reconstructed or not DECLARE_SOA_COLUMN(NSiliconHits, nSiliconHits, int); //! number of silicon hits DECLARE_SOA_COLUMN(NTPCHits, nTPCHits, int); //! number of tpc hits +DECLARE_SOA_COLUMN(PdgCode, pdgCode, int); //! PDG code of the linked truth MC particle } // namespace track_alice3 DECLARE_SOA_TABLE(TracksAlice3, "AOD", "TRACKSALICE3", track_alice3::IsReconstructed); using TrackAlice3 = TracksAlice3::iterator; +DECLARE_SOA_TABLE(TracksAlice3Pdg, "AOD", "TRACKSALICE3PDG", + track_alice3::PdgCode); +using TrackAlice3Pdg = TracksAlice3Pdg::iterator; + DECLARE_SOA_TABLE(TracksExtraA3, "AOD", "TracksExtraA3", track_alice3::NSiliconHits, track_alice3::NTPCHits); using TrackExtraA3 = TracksExtraA3::iterator; +namespace mcparticle_alice3 +{ +DECLARE_SOA_COLUMN(NHits, nHits, int); //! number of silicon hits +DECLARE_SOA_COLUMN(Charge, charge, float); //! particle charge +} // namespace mcparticle_alice3 +DECLARE_SOA_TABLE(MCParticlesExtraA3, "AOD", "MCParticlesExtraA3", + mcparticle_alice3::NHits, + mcparticle_alice3::Charge); +using MCParticleExtraA3 = MCParticlesExtraA3::iterator; } // namespace o2::aod #endif // ALICE3_DATAMODEL_TRACKSALICE3_H_ diff --git a/ALICE3/TableProducer/CMakeLists.txt b/ALICE3/TableProducer/CMakeLists.txt index 280a15f59ce..233088b6d5e 100644 --- a/ALICE3/TableProducer/CMakeLists.txt +++ b/ALICE3/TableProducer/CMakeLists.txt @@ -58,7 +58,7 @@ o2physics_add_dpl_workflow(alice3-hf-tree-creator-3prong o2physics_add_dpl_workflow(alice3-tracking-translator SOURCES alice3TrackingTranslator.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::ALICE3Core COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(alice3-dq-table-maker @@ -68,5 +68,5 @@ o2physics_add_dpl_workflow(alice3-dq-table-maker o2physics_add_dpl_workflow(alice3strangenessfinder SOURCES alice3strangenessFinder.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::ALICE3Core COMPONENT_NAME Analysis) diff --git a/ALICE3/TableProducer/alice3TrackingTranslator.cxx b/ALICE3/TableProducer/alice3TrackingTranslator.cxx index 5fc4728de98..706d5cba462 100644 --- a/ALICE3/TableProducer/alice3TrackingTranslator.cxx +++ b/ALICE3/TableProducer/alice3TrackingTranslator.cxx @@ -42,11 +42,6 @@ #include #include -#ifdef __CLING__ -#pragma link C++ class std::vector < std::vector < unsigned int>> + ; -#pragma link C++ class std::vector < std::vector < std::uint32_t>> + ; -#endif - TString inputPath; struct Alice3TrackingTranslator { @@ -61,7 +56,9 @@ struct Alice3TrackingTranslator { o2::framework::Produces tableTracksDCACov; o2::framework::Produces tableCollisionsAlice3; o2::framework::Produces tableTracksAlice3; + o2::framework::Produces tableTracksAlice3Pdg; o2::framework::Produces tableTracksExtraA3; + o2::framework::Produces tableMCParticlesExtraA3; o2::framework::Produces tableStoredTracksExtra; o2::framework::Produces tableTrackSelection; @@ -71,13 +68,12 @@ struct Alice3TrackingTranslator { o2::framework::Produces tableOTFLUTConfigId; o2::framework::Configurable maxCollisions{"maxCollisions", -1000, "Maximum number of collisions translated"}; + o2::framework::Configurable addDaughterInfo{"addDaughterInfo", false, "Add daughter particle information to the MC truth output tables"}; void init(o2::framework::InitContext&) { // Initialization if needed LOG(info) << "Alice3TrackingTranslator init called"; - // Load dictionary for nested vector - gInterpreter->GenerateDictionary("vector>", "vector"); } #define SETADDRESS(branchname, branchvar) \ @@ -123,6 +119,12 @@ struct Alice3TrackingTranslator { SETADDRESS("pz", m_pz); SETADDRESS("m", m_m); SETADDRESS("p", m_p); + SETADDRESS("q", m_q); + SETADDRESS("number_of_hits", m_number_of_hits); + SETADDRESS("particle_id", m_particleId); + if (mTree->GetBranchStatus("mother_particle_id")) { + SETADDRESS("mother_particle_id", m_motherId); + } } std::vector* m_particle_type = nullptr; std::vector* m_vx = nullptr; @@ -134,6 +136,10 @@ struct Alice3TrackingTranslator { std::vector* m_pz = nullptr; std::vector* m_m = nullptr; std::vector* m_p = nullptr; + std::vector* m_number_of_hits = nullptr; + std::vector* m_particleId = nullptr; + std::vector* m_q = nullptr; + std::vector* m_motherId = nullptr; }; struct TrackStruct : public FileStruct { @@ -155,8 +161,8 @@ struct Alice3TrackingTranslator { SETADDRESS("eQOP_fit", m_eQOP_fit); SETADDRESS("eT_fit", m_eT_fit); SETADDRESS("nMajorityHits", m_nMajorityHits); - // SETADDRESS("majorityParticleId", m_majorityParticleId); - mTree->SetBranchAddress("majorityParticleId", &m_majorityParticleId); + SETADDRESS("majorityParticleId", m_majorityParticleId); + // mTree->SetBranchAddress("majorityParticleId", &m_majorityParticleId); SETADDRESS("t_charge", m_t_charge); SETADDRESS("t_vx", m_t_vx); SETADDRESS("t_vy", m_t_vy); @@ -187,20 +193,20 @@ struct Alice3TrackingTranslator { std::vector* m_eT_fit = nullptr; // time // The majority truth particle info - std::vector* m_nMajorityHits = nullptr; /// The number of hits from majority particle - std::vector>* m_majorityParticleId = nullptr; /// The particle Id of the majority particle - std::vector* m_t_charge = nullptr; /// Charge of majority particle - std::vector* m_t_time = nullptr; /// Time of majority particle - std::vector* m_t_vx = nullptr; /// Vertex x positions of majority particle - std::vector* m_t_vy = nullptr; /// Vertex y positions of majority particle - std::vector* m_t_vz = nullptr; /// Vertex z positions of majority particle - std::vector* m_t_px = nullptr; /// Initial momenta m_px of majority particle - std::vector* m_t_py = nullptr; /// Initial momenta m_py of majority particle - std::vector* m_t_pz = nullptr; /// Initial momenta m_pz of majority particle - std::vector* m_t_theta = nullptr; /// Initial momenta theta of majority particle - std::vector* m_t_phi = nullptr; /// Initial momenta phi of majority particle - std::vector* m_t_pT = nullptr; /// Initial momenta pT of majority particle - std::vector* m_t_eta = nullptr; /// Initial momenta eta of majority particle + std::vector* m_nMajorityHits = nullptr; /// The number of hits from majority particle + std::vector* m_majorityParticleId = nullptr; /// The particle Id of the majority particle + std::vector* m_t_charge = nullptr; /// Charge of majority particle + std::vector* m_t_time = nullptr; /// Time of majority particle + std::vector* m_t_vx = nullptr; /// Vertex x positions of majority particle + std::vector* m_t_vy = nullptr; /// Vertex y positions of majority particle + std::vector* m_t_vz = nullptr; /// Vertex z positions of majority particle + std::vector* m_t_px = nullptr; /// Initial momenta m_px of majority particle + std::vector* m_t_py = nullptr; /// Initial momenta m_py of majority particle + std::vector* m_t_pz = nullptr; /// Initial momenta m_pz of majority particle + std::vector* m_t_theta = nullptr; /// Initial momenta theta of majority particle + std::vector* m_t_phi = nullptr; /// Initial momenta phi of majority particle + std::vector* m_t_pT = nullptr; /// Initial momenta pT of majority particle + std::vector* m_t_eta = nullptr; /// Initial momenta eta of majority particle std::vector* m_majorityParticlePDG = nullptr; // IA }; @@ -213,6 +219,28 @@ struct Alice3TrackingTranslator { } std::vector* barcode = nullptr; }; + void addMCParticle(int collIndex, ParticleStruct& fileParticles, int iParticle, uint8_t flags, int firstMother, int firstDaughter, int numberOfHits) + { + int mothers[2] = {firstMother, -1}; + int daughters[2] = {firstDaughter, -1}; + tableStoredMcParticles(collIndex, // mcCollisionId + fileParticles.m_particle_type->at(iParticle), // pdgCode + 0, // statusCode + flags, // flags + mothers, // mothersIds + daughters, // daughtersIdSlice + 1.0f, // weight + fileParticles.m_px->at(iParticle), // m_px + fileParticles.m_py->at(iParticle), // m_py + fileParticles.m_pz->at(iParticle), // m_pz + std::hypot(fileParticles.m_p->at(iParticle), fileParticles.m_m->at(iParticle)), // e + fileParticles.m_vx->at(iParticle), // m_vx + fileParticles.m_vy->at(iParticle), // m_vy + fileParticles.m_vz->at(iParticle), // m_vz + fileParticles.m_vt->at(iParticle)); // m_vt + tableMCParticlesExtraA3(numberOfHits, // number_of_hits + fileParticles.m_q->at(iParticle)); // charge + }; void process(o2::aod::BCs const&) { @@ -248,23 +276,34 @@ struct Alice3TrackingTranslator { LOG(info) << "Processing file: " << justFilename.Data(); files[justFilename.Data()] = filename; } - + LOG(info) << "All files loaded successfully"; // Now open the files to translate and read the trees - ParticleStruct fileParticles(files["particles_simulation.root"], "particles"); + ParticleStruct fileParticles(files["particles.root"], "particles"); + LOG(info) << "Particles loaded successfully"; + ParticleStruct fileParticlesSim(files["particles_simulation.root"], "particles"); + LOG(info) << "Particles Sim loaded successfully"; + std::string daughterFileName = addDaughterInfo ? "particles_decay.root" : "particles_simulation.root"; + ParticleStruct fileDaughterParticles(files[daughterFileName], "particles"); + LOG(info) << "Daughter particles loaded successfully from file " << daughterFileName; // FileStruct fileVertices(files["performance_vertexing.root"], "vertexing"); TrackStruct fileTracksummary(files["tracksummary_ambi.root"], "tracksummary"); // HitsStruct fileHits(files["hits.root"], "hits"); + LOG(info) << "Tracks loaded successfully"; const Long64_t kEvents = fileParticles.getEntries(); + int indexOfLastParticleAfterEvent = -1; for (Long64_t iEvent = 0; iEvent < kEvents; ++iEvent) { if (iEvent > 0 && maxCollisions.value > 0 && (iEvent % maxCollisions) == 0) { - LOG(info) << "Processing event " << iEvent << "/" << kEvents; + LOG(info) << "Stopping at event " << iEvent << "/" << kEvents; break; } fileParticles.setEventEntry(iEvent); // fileVertices.setEventEntry(iEvent); fileTracksummary.setEventEntry(iEvent); // fileHits.setEventEntry(iEvent); + if (addDaughterInfo) + fileDaughterParticles.setEventEntry(iEvent); + fileParticlesSim.setEventEntry(iEvent); LOG(info) << "Processing event " << iEvent << "/" << kEvents; @@ -275,6 +314,12 @@ struct Alice3TrackingTranslator { float collisionZ = 0.0f; tableOTFLUTConfigId(0); // dummy for the moment + + // Determine the collision ID for the new entry. + // If the table is empty, lastIndex() returns -1, so we start at 0. + // If it has entries, lastIndex() returns the index of the last element, so we use lastIndex() + 1. + int collisionId = tableCollisions.lastIndex() + 1; + tableCollisions(0, // bcId collisionX, // posX collisionY, // posY @@ -296,10 +341,6 @@ struct Alice3TrackingTranslator { tableCollisionsAlice3(0.f); // multDensity - // Fill MC particles - int mothers[2] = {-1, -1}; - int daughters[2] = {-1, -1}; - struct addedParticle { float px; float py; @@ -308,51 +349,92 @@ struct Alice3TrackingTranslator { float vy; float vz; }; - std::map> addedParticles; // Convert tracks from ACTS to ALICE format const size_t nParticlesGen = fileParticles.m_vx->size(); - const size_t nParticles = fileTracksummary.m_t_vx->size(); + const size_t nParticlesSim = fileParticlesSim.m_vx->size(); + const size_t nDaughterParticles = fileDaughterParticles.m_vx->size(); const size_t nTracks = fileTracksummary.m_eLOC0_fit->size(); + std::vector idMCparticles; + for (size_t iTrack = 0; iTrack < nTracks; ++iTrack) { - LOG(info) << "Processing track " << iTrack << "/" << nTracks << " (nParticles=" << nParticles << ") nParticlesGen=" << nParticlesGen; + LOG(info) << "Processing track " << iTrack << "/" << nTracks << " (nParticlesSim=" << nParticlesSim << ") nParticlesGen=" << nParticlesGen; const size_t iParticle = iTrack; + if (iParticle == 0) { - tableMcCollisions(0, // mccollision::BCId, - 0, // mccollision::GeneratorsID, - fileTracksummary.m_t_vx->at(iParticle), // mccollision::PosX, - fileTracksummary.m_t_vy->at(iParticle), // mccollision::PosY, - fileTracksummary.m_t_vz->at(iParticle), // mccollision::PosZ - fileTracksummary.m_t_time->at(iParticle), // mccollision::T - 1.0f, // mccollision::Weight - 0.0f, // mccollision::ImpactParameter, - 0.f); // mccollision::EventPlaneAngle, + tableMcCollisions(0, // mccollision::BCId, + 0, // mccollision::GeneratorsID, + fileParticles.m_vx->at(iParticle), // mccollision::PosX, + fileParticles.m_vy->at(iParticle), // mccollision::PosY, + fileParticles.m_vz->at(iParticle), // mccollision::PosZ + fileParticles.m_vt->at(iParticle), // mccollision::T + 1.0f, // mccollision::Weight + 0.0f, // mccollision::ImpactParameter, + 0.f); // mccollision::EventPlaneAngle, } uint8_t flags = 0; - flags |= o2::aod::mcparticle::enums::PhysicalPrimary; - - addedParticles[fileTracksummary.m_majorityParticlePDG->at(iParticle)].push_back({fileTracksummary.m_t_px->at(iParticle), - fileTracksummary.m_t_py->at(iParticle), - fileTracksummary.m_t_pz->at(iParticle), - fileTracksummary.m_t_vx->at(iParticle), - fileTracksummary.m_t_vy->at(iParticle), - fileTracksummary.m_t_vz->at(iParticle)}); - - tableStoredMcParticles(tableMcCollisions.lastIndex(), // mcCollisionId - fileTracksummary.m_majorityParticlePDG->at(iParticle), // pdgCode - 0, // statusCode - flags, // flags - mothers, // mothersIds - daughters, // daughtersIdSlice - 1.0f, // weight - fileTracksummary.m_t_px->at(iParticle), // m_px - fileTracksummary.m_t_py->at(iParticle), // m_py - fileTracksummary.m_t_pz->at(iParticle), // m_pz - 0, // e - fileTracksummary.m_t_vx->at(iParticle), // m_vx - fileTracksummary.m_t_vy->at(iParticle), // m_vy - fileTracksummary.m_t_vz->at(iParticle), // m_vz - fileTracksummary.m_t_time->at(iParticle)); // m_vt + ULong64_t idMCTrueParticle = fileTracksummary.m_majorityParticleId->at(iParticle); + int32_t mcParticleId = -1; + int pdgCode = -1; + + for (size_t iMC = 0; iMC < nParticlesGen; ++iMC) { + if (fileParticles.m_particleId->at(iMC) == idMCTrueParticle) { + if (count(idMCparticles.begin(), idMCparticles.end(), fileParticles.m_particleId->at(iMC)) > 0) { + continue; + } + idMCparticles.push_back(fileParticles.m_particleId->at(iMC)); + flags |= o2::aod::mcparticle::enums::PhysicalPrimary; + int nHits = 0; + for (size_t iPartSim = 0; iPartSim < nParticlesSim; ++iPartSim) { + if (fileParticlesSim.m_particleId->at(iPartSim) == fileParticles.m_particleId->at(iMC)) { + nHits = fileParticlesSim.m_number_of_hits->at(iPartSim); + break; + } + } + addMCParticle(tableMcCollisions.lastIndex(), fileParticles, iMC, flags, -1, -1, nHits); + mcParticleId = tableStoredMcParticles.lastIndex(); + pdgCode = fileParticles.m_particle_type->at(iMC); + break; + } + } + if (addDaughterInfo) { + for (size_t iMC = 0; iMC < nDaughterParticles; ++iMC) { + if (fileDaughterParticles.m_particleId->at(iMC) == idMCTrueParticle) { + if (count(idMCparticles.begin(), idMCparticles.end(), fileDaughterParticles.m_particleId->at(iMC)) > 0) { + break; + } + + int nHits = 0; + for (size_t iPartSim = 0; iPartSim < nParticlesSim; ++iPartSim) { + if (fileParticlesSim.m_particleId->at(iPartSim) == fileDaughterParticles.m_particleId->at(iMC)) { + nHits = fileParticlesSim.m_number_of_hits->at(iPartSim); + break; + } + } + for (size_t iMother = 0; iMother < nParticlesGen; ++iMother) { + if (fileDaughterParticles.m_motherId->at(iMC) == fileParticles.m_particleId->at(iMother)) { + if (count(idMCparticles.begin(), idMCparticles.end(), fileParticles.m_particleId->at(iMother)) > 0) { + break; + } + idMCparticles.push_back(fileParticles.m_particleId->at(iMother)); + uint8_t flagsMother = o2::aod::mcparticle::enums::PhysicalPrimary; + addMCParticle(tableMcCollisions.lastIndex(), fileParticles, iMother, flagsMother, -1, tableStoredMcParticles.lastIndex() + 2, 0); + break; + } + } + int motherId = -1; + if (count(idMCparticles.begin(), idMCparticles.end(), fileDaughterParticles.m_motherId->at(iMC)) > 0) { + auto it = find(idMCparticles.begin(), idMCparticles.end(), fileDaughterParticles.m_motherId->at(iMC)); + motherId = it - idMCparticles.begin() + indexOfLastParticleAfterEvent + 1; + } + idMCparticles.push_back(fileDaughterParticles.m_particleId->at(iMC)); + addMCParticle(tableMcCollisions.lastIndex(), fileDaughterParticles, iMC, flags, motherId, -1, nHits); + mcParticleId = tableStoredMcParticles.lastIndex(); + pdgCode = fileDaughterParticles.m_particle_type->at(iMC); + break; + } + } + } // Extract ACTS track parameters const float phi = fileTracksummary.m_ePHI_fit->at(iTrack); const float theta = fileTracksummary.m_eTHETA_fit->at(iTrack); @@ -409,7 +491,7 @@ struct Alice3TrackingTranslator { o2::track::TrackParCov trackParCov(x, alpha, trackParams, trackCov, charge); // Fill StoredTracks table (basic track parameters) - tableStoredTracks(tableCollisions.lastIndex(), // collisionId + tableStoredTracks(collisionId, // collisionId o2::aod::track::TrackTypeEnum::Track, // trackType trackParCov.getX(), // x trackParCov.getAlpha(), // alpha @@ -460,8 +542,6 @@ struct Alice3TrackingTranslator { // Fill MC track labels // Get particle linkage from hits using the majority hit index - int32_t mcParticleId = -1; // Default to invalid particle ID - mcParticleId = tableStoredMcParticles.lastIndex(); // Temporary: link all tracks to the last added MC particle // if (fileTracksummary.nMajorityHits && iTrack < fileTracksummary.nMajorityHits->size()) { // unsigned int hitIndex = fileTracksummary.nMajorityHits->at(iTrack); // if (fileHits.barcode && hitIndex < fileHits.barcode->size()) { @@ -489,7 +569,8 @@ struct Alice3TrackingTranslator { 0.0f); // sigmaDcaZ2 // Fill ALICE3 specific tables - tableTracksAlice3(true); // isReconstructed + tableTracksAlice3(true); // isReconstructed + tableTracksAlice3Pdg(pdgCode); // PdgCode to the linked MC truth particle tableTracksExtraA3(m_nMeasurements, // nSiliconHits (using m_nMeasurements as proxy) 0); // nTPCHits @@ -546,43 +627,64 @@ struct Alice3TrackingTranslator { } for (size_t iParticle = 0; iParticle < nParticlesGen; ++iParticle) { - + if (iParticle == 0 && nTracks == 0) { + tableMcCollisions(0, // mccollision::BCId, + 0, // mccollision::GeneratorsID, + fileParticles.m_vx->at(iParticle), // mccollision::PosX, + fileParticles.m_vy->at(iParticle), // mccollision::PosY, + fileParticles.m_vz->at(iParticle), // mccollision::PosZ + fileParticles.m_vt->at(iParticle), // mccollision::T + 1.0f, // mccollision::Weight + 0.0f, // mccollision::ImpactParameter, + 0.f); // mccollision::EventPlaneAngle, + } + if (idMCparticles.end() != std::find(idMCparticles.begin(), idMCparticles.end(), fileParticles.m_particleId->at(iParticle))) { + // Already added via track + continue; + } uint8_t flags = 0; flags |= o2::aod::mcparticle::enums::PhysicalPrimary; - bool alreadyAdded = false; - for (const auto& ap : addedParticles[fileParticles.m_particle_type->at(iParticle)]) { - if (std::abs(ap.px - fileParticles.m_px->at(iParticle)) <= 1.e-5 && - std::abs(ap.py - fileParticles.m_py->at(iParticle)) <= 1.e-5 && - std::abs(ap.pz - fileParticles.m_pz->at(iParticle)) <= 1.e-5 && - std::abs(ap.vx - fileParticles.m_vx->at(iParticle)) <= 1.e-5 && - std::abs(ap.vy - fileParticles.m_vy->at(iParticle)) <= 1.e-5 && - std::abs(ap.vz - fileParticles.m_vz->at(iParticle)) <= 1.e-5) { - alreadyAdded = true; + + int nHits = 0; + for (size_t iPartSim = 0; iPartSim < nParticlesSim; ++iPartSim) { + if (fileParticlesSim.m_particleId->at(iPartSim) == fileParticles.m_particleId->at(iParticle)) { + nHits = fileParticlesSim.m_number_of_hits->at(iPartSim); break; } } - if (alreadyAdded) { - continue; + addMCParticle(tableMcCollisions.lastIndex(), fileParticles, iParticle, flags, -1, -1, nHits); + idMCparticles.push_back(fileParticles.m_particleId->at(iParticle)); + } + if (addDaughterInfo) { + for (size_t iParticle = 0; iParticle < nDaughterParticles; ++iParticle) { + if (idMCparticles.end() != std::find(idMCparticles.begin(), idMCparticles.end(), fileDaughterParticles.m_particleId->at(iParticle))) { + // Already added via track + continue; + } + uint8_t flags = 0; + int nHits = 0; + for (size_t iPartSim = 0; iPartSim < nParticlesSim; ++iPartSim) { + if (fileParticlesSim.m_particleId->at(iPartSim) == fileDaughterParticles.m_particleId->at(iParticle)) { + nHits = fileParticlesSim.m_number_of_hits->at(iPartSim); + break; + } + } + int motherId = -1; + for (size_t iMother = 0; iMother < nParticlesGen; ++iMother) { + if (fileDaughterParticles.m_motherId->at(iParticle) == fileParticles.m_particleId->at(iMother)) { + if (count(idMCparticles.begin(), idMCparticles.end(), fileDaughterParticles.m_motherId->at(iParticle)) > 0) { + auto it = find(idMCparticles.begin(), idMCparticles.end(), fileDaughterParticles.m_motherId->at(iParticle)); + motherId = it - idMCparticles.begin() + indexOfLastParticleAfterEvent + 1; + } + } + } + addMCParticle(tableMcCollisions.lastIndex(), fileDaughterParticles, iParticle, flags, motherId, -1, nHits); } - - tableStoredMcParticles(tableMcCollisions.lastIndex(), // mcCollisionId - fileParticles.m_particle_type->at(iParticle), // pdgCode - 0, // statusCode - flags, // flags - mothers, // mothersIds - daughters, // daughtersIdSlice - 1.0f, // weight - fileParticles.m_px->at(iParticle), // m_px - fileParticles.m_py->at(iParticle), // m_py - fileParticles.m_pz->at(iParticle), // m_pz - std::hypot(fileParticles.m_p->at(iParticle), fileParticles.m_m->at(iParticle)), // e - fileParticles.m_vx->at(iParticle), // m_vx - fileParticles.m_vy->at(iParticle), // m_vy - fileParticles.m_vz->at(iParticle), // m_vz - fileParticles.m_vt->at(iParticle)); // m_vt } - LOG(info) << "Event " << iEvent << ": has " << nTracks << " tracks and " << nParticles << " particles."; + LOG(info) << "Event " << iEvent << ": has " << nTracks << " tracks, " << nParticlesGen << " particles " << nDaughterParticles << " daughter particles, " << nParticlesSim << " propagated particles."; + LOG(info) << "Total numbers of stored MC particles: " << tableStoredMcParticles.lastIndex() + 1; + indexOfLastParticleAfterEvent = tableStoredMcParticles.lastIndex(); } } }; diff --git a/ALICE3/TableProducer/alice3strangenessFinder.cxx b/ALICE3/TableProducer/alice3strangenessFinder.cxx index 5419509d951..74cc26d6b49 100644 --- a/ALICE3/TableProducer/alice3strangenessFinder.cxx +++ b/ALICE3/TableProducer/alice3strangenessFinder.cxx @@ -21,11 +21,13 @@ #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "ALICE3/Core/TrackUtilities.h" #include "ALICE3/DataModel/OTFPIDTrk.h" #include "ALICE3/DataModel/OTFRICH.h" #include "ALICE3/DataModel/OTFStrangeness.h" #include "ALICE3/DataModel/OTFTOF.h" #include "ALICE3/DataModel/tracksAlice3.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -36,6 +38,10 @@ #include "ReconstructionDataFormats/Track.h" #include #include +#include +#include + +#include #include @@ -45,8 +51,8 @@ using namespace o2::framework; using namespace o2::constants::physics; using Alice3TracksWPid = soa::Join; -using Alice3Tracks = soa::Join; - +using Alice3Tracks = soa::Join; +using Alice3MCParticles = soa::Join; struct Alice3strangenessFinder { SliceCache cache; @@ -57,9 +63,12 @@ struct Alice3strangenessFinder { Produces tableCascadeCores; Produces tableCascadeIndices; + Configurable buildCascade{"buildCascade", false, "build cascade candidates"}; + Configurable nSigmaTOF{"nSigmaTOF", 5.0f, "Nsigma for TOF PID (if enabled)"}; Configurable dcaXYconstant{"dcaXYconstant", -1.0f, "[0] in |DCAxy| > [0]+[1]/pT"}; Configurable dcaXYpTdep{"dcaXYpTdep", 0.0, "[1] in |DCAxy| > [0]+[1]/pT"}; + ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.025f, 0.05f, 0.075f, 0.1f, 0.125f, 0.15f, 0.175f, 0.2f, 0.225f, 0.25f, 0.275f, 0.3f, 0.325f, 0.35f, 0.375f, 0.4f, 0.425f, 0.45f, 0.475f, 0.5f, 0.525f, 0.55f, 0.575f, 0.6f, 0.625f, 0.65f, 0.675f, 0.7f, 0.725f, 0.75f, 0.775f, 0.8f, 0.82f, 0.85f, 0.875f, 0.9f, 0.925f, 0.95f, 0.975f, 1.0f, 1.05f, 1.1f}, "pt axis for QA histograms"}; Configurable bachMinConstDCAxy{"bachMinConstDCAxy", -1.0f, "[0] in |DCAxy| > [0]+[1]/pT"}; Configurable bachMinPtDepDCAxy{"bachMinPtDepDCAxy", 0.0, "[1] in |DCAxy| > [0]+[1]/pT"}; @@ -73,26 +82,31 @@ struct Alice3strangenessFinder { Configurable maxR{"maxR", 150., "reject PCA's above this radius"}; Configurable maxDZIni{"maxDZIni", 5, "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; Configurable maxDXYIni{"maxDXYIni", 4, "reject (if>0) PCA candidate if tracks DXY exceeds threshold"}; - Configurable maxVtxChi2{"maxVtxChi2", 2, "reject (if>0) vtx. chi2 above this value"}; + Configurable maxVtxChi2{"maxVtxChi2", 10, "reject (if>0) vtx. chi2 above this value"}; Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; Configurable acceptedLambdaMassWindow{"acceptedLambdaMassWindow", 0.2f, "accepted Lambda mass window around PDG mass"}; - // Operation and minimisation criteria + // Operation Configurable magneticField{"magneticField", 20.0f, "Magnetic field (in kilogauss)"}; - Configurable doDCAplotsD{"doDCAplotsD", true, "do daughter prong DCA plots for D mesons"}; - Configurable doDCAplots3Prong{"doDCAplots3Prong", true, "do daughter prong DCA plots for Lc baryons"}; - Configurable doTopoPlotsForSAndB{"doTopoPlotsForSAndB", true, "do topological variable distributions for S and B separately"}; - Configurable dcaDaughtersSelection{"dcaDaughtersSelection", 1000.0f, "DCA between daughters (cm)"}; Configurable mcSameMotherCheck{"mcSameMotherCheck", true, "check if tracks come from the same MC mother"}; // propagation options Configurable usePropagator{"usePropagator", false, "use external propagator"}; Configurable refitWithMatCorr{"refitWithMatCorr", false, "refit V0 applying material corrections"}; Configurable useCollinearV0{"useCollinearV0", true, "use collinear approximation for V0 fitting"}; + Configurable maxIter{"maxIter", 30, "maximum number of iterations for vertex fitter"}; + + // for the ACTS study + Configurable isK0Gun{"isK0Gun", false, "is K0s Monte Carlo gun used"}; + Configurable isLambdaGun{"isLambdaGun", true, "is Lambda Monte Carlo gun used"}; + Configurable skipFitter{"skipFitter", false, "calculate V0 properties without calling the DCA fitter, using only the track parameters "}; + Configurable useOriginalTrackParams{"useOriginalTrackParams", false, "use original track parameters instead of the ones propagated to PCA (effective only if skipFitter is false) and for MC truth info"}; o2::vertexing::DCAFitterN<2> fitter; o2::vertexing::DCAFitterN<3> fitter3; + Service pdgDB; + // partitions for D mesons Partition positiveSecondaryTracks = aod::track::signed1Pt > 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); @@ -100,7 +114,8 @@ struct Alice3strangenessFinder { aod::track::signed1Pt < 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); Partition bachelorTracks = nabs(aod::track::dcaXY) > bachMinConstDCAxy + bachMinPtDepDCAxy* nabs(aod::track::signed1Pt) && nabs(aod::track::dcaZ) > bachMinConstDCAz + bachMinPtDepDCAz* nabs(aod::track::signed1Pt); - + Partition positiveMCParticles = aod::mcparticle_alice3::charge > 0.0f; + Partition negativeMCParticles = aod::mcparticle_alice3::charge < 0.0f; // Partition negativeSecondaryPions = nabs(aod::upgrade_tof::nSigmaPionInnerTOF) < nSigmaTOF && nabs(aod::upgrade_tof::nSigmaPionOuterTOF) < nSigmaTOF && aod::track::signed1Pt < 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); // Partition positiveSecondaryPions = nabs(aod::upgrade_tof::nSigmaPionInnerTOF) < nSigmaTOF && nabs(aod::upgrade_tof::nSigmaPionOuterTOF) < nSigmaTOF && aod::track::signed1Pt > 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); // Partition secondaryProtons = nabs(aod::upgrade_tof::nSigmaProtonInnerTOF) < nSigmaTOF && nabs(aod::upgrade_tof::nSigmaProtonOuterTOF) < nSigmaTOF && aod::track::signed1Pt > 0.0f && nabs(aod::track::dcaXY) > dcaXYconstant + dcaXYpTdep* nabs(aod::track::signed1Pt); @@ -134,8 +149,15 @@ struct Alice3strangenessFinder { fitter.setUsePropagator(usePropagator); fitter.setRefitWithMatCorr(refitWithMatCorr); fitter.setCollinear(useCollinearV0); + fitter.setMaxIter(maxIter); fitter.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrNONE); + histos.add("hFitterQA", "", kTH1D, {{10, 0, 10}}); // For QA reasons, counting found candidates at different stages + histos.add("hPtPosDau", "", kTH1D, {axisPt}); + histos.add("hPtNegDau", "", kTH1D, {axisPt}); + histos.add("hPtPosDauAfterV0Finding", "", kTH2D, {axisPt, axisPt}); + histos.add("hPtNegDauAfterV0Finding", "", kTH2D, {axisPt, axisPt}); + histos.add("hEventCounter", "", kTH1D, {{1, 0, 2}}); // counting processed events auto hV0Counter = histos.add("hV0Counter", "hV0Counter", kTH1D, {{4, 0, 4}}); hV0Counter->GetXaxis()->SetBinLabel(1, "K0S"); hV0Counter->GetXaxis()->SetBinLabel(2, "Lambda"); @@ -148,6 +170,9 @@ struct Alice3strangenessFinder { hCascadeCounter->GetXaxis()->SetBinLabel(3, "Omega"); hCascadeCounter->GetXaxis()->SetBinLabel(4, "AntiOmega"); hCascadeCounter->GetXaxis()->SetBinLabel(5, "Misidentified"); + histos.add("hRadiusVsHitsNeg", "", kTH2D, {{400, 0, 400}, {12, 0.5, 12.5}}); // radius vs hist for MC studies + histos.add("hRadiusVsHitsPos", "", kTH2D, {{400, 0, 400}, {12, 0.5, 12.5}}); // radius vs hist for MC studies + histos.print(); } float calculateDCAStraightToPV(float X, float Y, float Z, float Px, float Py, float Pz, float pvX, float pvY, float pvZ) @@ -180,65 +205,104 @@ struct Alice3strangenessFinder { template bool buildDecayCandidateTwoBody(TTrackType const& t0, TTrackType const& t1, std::array vtx, Candidate& thisCandidate) { - //}-{}-{}-{}-{}-{}-{}-{}-{}-{} - // Move close to minima - int nCand = 0; - try { - nCand = fitter.process(t0, t1); - } catch (...) { - return false; - } - if (nCand == 0) { - return false; - } - //}-{}-{}-{}-{}-{}-{}-{}-{}-{} - if (!fitter.isPropagateTracksToVertexDone() && !fitter.propagateTracksToVertex()) { - LOG(debug) << "RejProp failed"; - return false; - } - o2::track::TrackParCov t0New = fitter.getTrack(0); - o2::track::TrackParCov t1New = fitter.getTrack(1); - t0New.getPxPyPzGlo(thisCandidate.pDau0); - t1New.getPxPyPzGlo(thisCandidate.pDau1); - - thisCandidate.dcaDau = std::sqrt(fitter.getChi2AtPCACandidate()); - thisCandidate.p[0] = thisCandidate.pDau0[0] + thisCandidate.pDau1[0]; - thisCandidate.p[1] = thisCandidate.pDau0[1] + thisCandidate.pDau1[1]; - thisCandidate.p[2] = thisCandidate.pDau0[2] + thisCandidate.pDau1[2]; - - const auto posSV = fitter.getPCACandidatePos(); - thisCandidate.posSV[0] = posSV[0]; - thisCandidate.posSV[1] = posSV[1]; - thisCandidate.posSV[2] = posSV[2]; - - std::array covA = {0}; - std::array covB = {0}; - fitter.getTrack(0).getCovXYZPxPyPzGlo(covA); - fitter.getTrack(1).getCovXYZPxPyPzGlo(covB); - - static constexpr std::array MomentumIndices = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component - for (size_t i = 0; i < MomentumIndices.size(); i++) { - int j = MomentumIndices[i]; - thisCandidate.parentTrackCovMatrix[j] = covA[j] + covB[j]; - } + histos.fill(HIST("hPtNegDau"), t1.getPt()); + histos.fill(HIST("hPtPosDau"), t0.getPt()); + + if (!skipFitter) { - auto covVtx = fitter.calcPCACovMatrix(); - thisCandidate.parentTrackCovMatrix[0] = covVtx(0, 0); - thisCandidate.parentTrackCovMatrix[1] = covVtx(1, 0); - thisCandidate.parentTrackCovMatrix[2] = covVtx(1, 1); - thisCandidate.parentTrackCovMatrix[3] = covVtx(2, 0); - thisCandidate.parentTrackCovMatrix[4] = covVtx(2, 1); - thisCandidate.parentTrackCovMatrix[5] = covVtx(2, 2); - - thisCandidate.eta = RecoDecay::eta(std::array{thisCandidate.p[0], thisCandidate.p[1], thisCandidate.p[2]}); - thisCandidate.cosPA = RecoDecay::cpa(vtx, std::array{thisCandidate.posSV[0], thisCandidate.posSV[1], thisCandidate.posSV[2]}, - std::array{thisCandidate.p[0], thisCandidate.p[1], thisCandidate.p[2]}); - thisCandidate.dcaToPV = calculateDCAStraightToPV(thisCandidate.posSV[0], thisCandidate.posSV[1], thisCandidate.posSV[2], - thisCandidate.p[0], thisCandidate.p[1], thisCandidate.p[2], - vtx[0], vtx[1], vtx[2]); - - return true; + histos.fill(HIST("hFitterQA"), 0.5); + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + // Move close to minima + int nCand = 0; + try { + nCand = fitter.process(t0, t1); + } catch (...) { + return false; + } + histos.fill(HIST("hFitterQA"), 1.5); + if (nCand == 0) { + LOG(info) << "0 candidates found by fitter"; + return false; + } + histos.fill(HIST("hFitterQA"), 2.5); + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + if (!fitter.isPropagateTracksToVertexDone() && !fitter.propagateTracksToVertex()) { + LOG(info) << "RejProp failed"; + return false; + } + histos.fill(HIST("hFitterQA"), 3.5); + o2::track::TrackParCov t0New = fitter.getTrack(0); + o2::track::TrackParCov t1New = fitter.getTrack(1); + t0New.getPxPyPzGlo(thisCandidate.pDau0); + t1New.getPxPyPzGlo(thisCandidate.pDau1); + if (useOriginalTrackParams) { + t0.getPxPyPzGlo(thisCandidate.pDau0); + t1.getPxPyPzGlo(thisCandidate.pDau1); + } + histos.fill(HIST("hPtNegDauAfterV0Finding"), std::sqrt(thisCandidate.pDau1[0] * thisCandidate.pDau1[0] + thisCandidate.pDau1[1] + thisCandidate.pDau1[1]), t1.getPt()); + histos.fill(HIST("hPtPosDauAfterV0Finding"), std::sqrt(thisCandidate.pDau0[0] * thisCandidate.pDau0[0] + thisCandidate.pDau0[1] + thisCandidate.pDau0[1]), t0.getPt()); + + thisCandidate.dcaDau = std::sqrt(fitter.getChi2AtPCACandidate()); + thisCandidate.p[0] = thisCandidate.pDau0[0] + thisCandidate.pDau1[0]; + thisCandidate.p[1] = thisCandidate.pDau0[1] + thisCandidate.pDau1[1]; + thisCandidate.p[2] = thisCandidate.pDau0[2] + thisCandidate.pDau1[2]; + const auto posSV = fitter.getPCACandidatePos(); + thisCandidate.posSV[0] = posSV[0]; + thisCandidate.posSV[1] = posSV[1]; + thisCandidate.posSV[2] = posSV[2]; + + std::array covA = {0}; + std::array covB = {0}; + fitter.getTrack(0).getCovXYZPxPyPzGlo(covA); + fitter.getTrack(1).getCovXYZPxPyPzGlo(covB); + + static constexpr std::array MomentumIndices = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (size_t i = 0; i < MomentumIndices.size(); i++) { + int j = MomentumIndices[i]; + thisCandidate.parentTrackCovMatrix[j] = covA[j] + covB[j]; + } + + auto covVtx = fitter.calcPCACovMatrix(); + thisCandidate.parentTrackCovMatrix[0] = covVtx(0, 0); + thisCandidate.parentTrackCovMatrix[1] = covVtx(1, 0); + thisCandidate.parentTrackCovMatrix[2] = covVtx(1, 1); + thisCandidate.parentTrackCovMatrix[3] = covVtx(2, 0); + thisCandidate.parentTrackCovMatrix[4] = covVtx(2, 1); + thisCandidate.parentTrackCovMatrix[5] = covVtx(2, 2); + thisCandidate.parentTrackCovMatrix[0] = 0; + thisCandidate.parentTrackCovMatrix[1] = 0; + thisCandidate.parentTrackCovMatrix[2] = 0; + thisCandidate.parentTrackCovMatrix[3] = 0; + thisCandidate.parentTrackCovMatrix[4] = 0; + thisCandidate.parentTrackCovMatrix[5] = 0; + + thisCandidate.eta = RecoDecay::eta(std::array{thisCandidate.p[0], thisCandidate.p[1], thisCandidate.p[2]}); + thisCandidate.cosPA = RecoDecay::cpa(vtx, std::array{thisCandidate.posSV[0], thisCandidate.posSV[1], thisCandidate.posSV[2]}, + std::array{thisCandidate.p[0], thisCandidate.p[1], thisCandidate.p[2]}); + thisCandidate.dcaToPV = calculateDCAStraightToPV(thisCandidate.posSV[0], thisCandidate.posSV[1], thisCandidate.posSV[2], + thisCandidate.p[0], thisCandidate.p[1], thisCandidate.p[2], + vtx[0], vtx[1], vtx[2]); + + return true; + } else { + t0.getPxPyPzGlo(thisCandidate.pDau0); + t1.getPxPyPzGlo(thisCandidate.pDau1); + thisCandidate.dcaDau = 0; + thisCandidate.p[0] = thisCandidate.pDau0[0] + thisCandidate.pDau1[0]; + thisCandidate.p[1] = thisCandidate.pDau0[1] + thisCandidate.pDau1[1]; + thisCandidate.p[2] = thisCandidate.pDau0[2] + thisCandidate.pDau1[2]; + thisCandidate.posSV[0] = 0; + thisCandidate.posSV[1] = 0; + thisCandidate.posSV[2] = 0; + thisCandidate.eta = RecoDecay::eta(std::array{thisCandidate.p[0], thisCandidate.p[1], thisCandidate.p[2]}); + thisCandidate.cosPA = RecoDecay::cpa(vtx, std::array{thisCandidate.posSV[0], thisCandidate.posSV[1], thisCandidate.posSV[2]}, + std::array{thisCandidate.p[0], thisCandidate.p[1], thisCandidate.p[2]}); + thisCandidate.dcaToPV = calculateDCAStraightToPV(thisCandidate.posSV[0], thisCandidate.posSV[1], thisCandidate.posSV[2], + thisCandidate.p[0], thisCandidate.p[1], thisCandidate.p[2], + vtx[0], vtx[1], vtx[2]); + return true; + } } void processFindV0CandidateNoPid(aod::Collision const& collision, Alice3Tracks const&, aod::McParticles const&) @@ -248,12 +312,15 @@ struct Alice3strangenessFinder { auto bachelorTracksGrouped = bachelorTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); const std::array vtx = {collision.posX(), collision.posY(), collision.posZ()}; + histos.fill(HIST("hEventCounter"), 1.0); + for (auto const& posTrack : positiveSecondaryTracksGrouped) { if (!posTrack.isReconstructed()) { continue; // no ghost tracks } o2::track::TrackParCov pos = getTrackParCov(posTrack); + for (auto const& negTrack : negativeSecondaryTracksGrouped) { if (!negTrack.isReconstructed()) { continue; // no ghost tracks @@ -262,28 +329,32 @@ struct Alice3strangenessFinder { if (mcSameMotherCheck && !checkSameMother(posTrack, negTrack)) { continue; // keep only if same mother } - + if ((posTrack.pdgCode() != kPiPlus && negTrack.pdgCode() != kPiMinus) && isK0Gun) + continue; + if ((posTrack.pdgCode() != kProton && negTrack.pdgCode() != kPiMinus) && isLambdaGun) + continue; o2::track::TrackParCov neg = getTrackParCov(negTrack); Candidate v0cand; if (!buildDecayCandidateTwoBody(pos, neg, vtx, v0cand)) { continue; // failed at building candidate } - auto mcParticle1 = posTrack.template mcParticle_as(); - if (mcParticle1.pdgCode() == kK0Short) { - histos.fill(HIST("hV0Counter"), 0.5); - } else if (mcParticle1.pdgCode() == kLambda0) { - histos.fill(HIST("hV0Counter"), 1.5); - } else if (mcParticle1.pdgCode() == kLambda0Bar) { - histos.fill(HIST("hV0Counter"), 2.5); - } else { - histos.fill(HIST("hV0Counter"), 3.5); - } + // TODO: not all ACTS tracks have MC association, so this check is not possible for all candidates, fix is needed + // auto mcParticle1 = posTrack.template mcParticle_as(); + // if (mcParticle1.pdgCode() == kK0Short) { + // histos.fill(HIST("hV0Counter"), 0.5); + // } else if (mcParticle1.pdgCode() == kLambda0) { + // histos.fill(HIST("hV0Counter"), 1.5); + // } else if (mcParticle1.pdgCode() == kLambda0Bar) { + // histos.fill(HIST("hV0Counter"), 2.5); + // } else { + // histos.fill(HIST("hV0Counter"), 3.5); + // } v0CandidateIndices(collision.globalIndex(), posTrack.globalIndex(), negTrack.globalIndex(), - mcParticle1.globalIndex()); + -1); v0CandidateCores(v0cand.posSV[0], v0cand.posSV[1], v0cand.posSV[2], v0cand.pDau0[0], v0cand.pDau0[1], v0cand.pDau0[2], @@ -300,12 +371,14 @@ struct Alice3strangenessFinder { std::array{v0cand.pDau1[0], v0cand.pDau1[1], v0cand.pDau1[2]}}, std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton}); - const bool inLambdaMassWindow = std::abs(lambdaMassHypothesis - o2::constants::physics::MassLambda0) > acceptedLambdaMassWindow; - const bool inAntiLambdaMassWindow = std::abs(antiLambdaMassHypothesis - o2::constants::physics::MassLambda0) > acceptedLambdaMassWindow; - if (inLambdaMassWindow || inAntiLambdaMassWindow) { + const bool inLambdaMassWindow = std::abs(lambdaMassHypothesis - o2::constants::physics::MassLambda0) < acceptedLambdaMassWindow; + const bool inAntiLambdaMassWindow = std::abs(antiLambdaMassHypothesis - o2::constants::physics::MassLambda0) < acceptedLambdaMassWindow; + if (!buildCascade) { + continue; // not building cascades, so skip the rest + } + if (!inLambdaMassWindow && !inAntiLambdaMassWindow) { continue; // Likely not a lambda, should not be considered for cascade building } - for (const auto& bachTrack : bachelorTracksGrouped) { if (bachTrack.globalIndex() == posTrack.globalIndex() || bachTrack.globalIndex() == negTrack.globalIndex()) { continue; // avoid using any track that was already used @@ -380,6 +453,83 @@ struct Alice3strangenessFinder { } // end negTrack } // end posTrack } + void processMCTrueFromACTS(aod::McCollision const& collision, Alice3MCParticles const&) + { + + auto negativeMCParticlesGrouped = negativeMCParticles->sliceByCached(aod::mcparticle::mcCollisionId, collision.globalIndex(), cache); + auto positiveMCParticlesGrouped = positiveMCParticles->sliceByCached(aod::mcparticle::mcCollisionId, collision.globalIndex(), cache); + const std::array vtx = {collision.posX(), collision.posY(), collision.posZ()}; + + float radiusPos = 0.0f; + float radiusNeg = 0.0f; + bool isK0s = false; + bool isLambda = false; + bool isAntiLambda = false; + int v0PdgCode = 0; + int iPosPart = 0; + for (auto const& posParticle : positiveMCParticlesGrouped) { + radiusPos = std::hypot(posParticle.vx(), posParticle.vy()); + histos.fill(HIST("hRadiusVsHitsPos"), radiusPos, posParticle.nHits()); + for (auto const& negParticle : negativeMCParticlesGrouped) { + if (negParticle.pdgCode() != kPiMinus && negParticle.pdgCode() != kProtonBar) { + continue; + } + radiusNeg = std::hypot(negParticle.vx(), negParticle.vy()); + if (iPosPart == 0) { + histos.fill(HIST("hRadiusVsHitsNeg"), radiusNeg, negParticle.nHits()); + } + if (radiusPos == radiusNeg) { + isK0s = (posParticle.pdgCode() == kPiPlus && negParticle.pdgCode() == kPiMinus); + if (isK0s) + v0PdgCode = kK0Short; + isLambda = (posParticle.pdgCode() == kProton && negParticle.pdgCode() == kPiMinus); + if (isLambda) + v0PdgCode = kLambda0; + isAntiLambda = (posParticle.pdgCode() == kPiPlus && negParticle.pdgCode() == kProtonBar); + if (isAntiLambda) + v0PdgCode = kLambda0Bar; + if (isK0s || isLambda || isAntiLambda) { + if (!isK0s && isK0Gun) + continue; + if (!isLambda && isLambdaGun) + continue; + Candidate v0cand; + std::vector v0DecayVertex; + v0DecayVertex.push_back(negParticle.vx()); + v0DecayVertex.push_back(negParticle.vy()); + v0DecayVertex.push_back(negParticle.vz()); + TLorentzVector posLorVector = {posParticle.px(), posParticle.py(), posParticle.pz(), posParticle.e()}; + TLorentzVector negLorVector = {negParticle.px(), negParticle.py(), negParticle.pz(), negParticle.e()}; + o2::track::TrackParCov posParCov; + o2::track::TrackParCov negParCov; + o2::upgrade::convertTLorentzVectorToO2Track(1, posLorVector, v0DecayVertex, posParCov); + o2::upgrade::convertTLorentzVectorToO2Track(-1, negLorVector, v0DecayVertex, negParCov); + if (!buildDecayCandidateTwoBody(posParCov, negParCov, vtx, v0cand)) + continue; + v0CandidateIndices(collision.globalIndex(), + posParticle.globalIndex(), + negParticle.globalIndex(), + 0); + v0CandidateCores(v0cand.posSV[0], v0cand.posSV[1], v0cand.posSV[2], + v0cand.pDau0[0], v0cand.pDau0[1], v0cand.pDau0[2], + v0cand.pDau1[0], v0cand.pDau1[1], v0cand.pDau1[2], + v0cand.dcaDau, 0, 0, + v0cand.cosPA, v0cand.dcaToPV); + if (isK0s) { + histos.fill(HIST("hV0Counter"), 0.5); + } else if (isLambda) { + histos.fill(HIST("hV0Counter"), 1.5); + } else if (isAntiLambda) { + histos.fill(HIST("hV0Counter"), 2.5); + } else { + histos.fill(HIST("hV0Counter"), 3.5); + } + } + } + } + iPosPart++; + } + } // void processFindV0CandidateWithPid(aod::Collision const& collision, aod::McParticles const& mcParticles, Alice3TracksWPid const&) // { // auto negativeSecondaryPionsGrouped = negativeSecondaryPions->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); @@ -388,6 +538,7 @@ struct Alice3strangenessFinder { // auto secondaryAntiProtonsGrouped = secondaryAntiProtons->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); // } PROCESS_SWITCH(Alice3strangenessFinder, processFindV0CandidateNoPid, "find V0 without PID", true); + PROCESS_SWITCH(Alice3strangenessFinder, processMCTrueFromACTS, "process MC truth from ACTS", false); // PROCESS_SWITCH(alice3strangenessFinder, processFindV0CandidateWithPid, "find V0 with PID", false); }; diff --git a/ALICE3/Tasks/alice3Strangeness.cxx b/ALICE3/Tasks/alice3Strangeness.cxx index 8d005b59ebd..3b539616fab 100644 --- a/ALICE3/Tasks/alice3Strangeness.cxx +++ b/ALICE3/Tasks/alice3Strangeness.cxx @@ -246,10 +246,12 @@ struct Alice3Strangeness { { // if(collision.lutConfigId()!=idGeometry) // return; + float collisionZ = collision.posZ(); + histos.fill(HIST("hPVz"), collisionZ); for (auto const& v0 : v0Candidates) { - bool isK0 = (v0.mK0Short() - o2::constants::physics::MassK0Short) < selectionValues.acceptedK0MassWindow; - bool isLambda = (v0.mLambda() - o2::constants::physics::MassLambda0) < selectionValues.acceptedLambdaMassWindow; - bool isAntiLambda = (v0.mAntiLambda() - o2::constants::physics::MassLambda0) < selectionValues.acceptedLambdaMassWindow; + bool isK0 = std::abs(v0.mK0Short() - o2::constants::physics::MassK0Short) < selectionValues.acceptedK0MassWindow; + bool isLambda = std::abs(v0.mLambda() - o2::constants::physics::MassLambda0) < selectionValues.acceptedLambdaMassWindow; + bool isAntiLambda = std::abs(v0.mAntiLambda() - o2::constants::physics::MassLambda0) < selectionValues.acceptedLambdaMassWindow; histos.fill(HIST("reconstructedCandidates/hArmeterosBeforeAllSelections"), v0.alpha(), v0.qtArm()); histos.fill(HIST("hV0CandidateCounter"), 0.5);