diff --git a/PWGCF/Femto/Core/collisionBuilder.h b/PWGCF/Femto/Core/collisionBuilder.h index 311f5b14e0e..2d0aa6b532d 100644 --- a/PWGCF/Femto/Core/collisionBuilder.h +++ b/PWGCF/Femto/Core/collisionBuilder.h @@ -395,20 +395,10 @@ class CollisionBuilder void init(o2::framework::HistogramRegistry* registry, T1& confFilter, T2& confBits, T3& confRct, T4& confCcdb, T5& confTable, T6& initContext) { LOG(info) << "Initialize femto collision builder..."; - mProducedCollisions = utils::enableTable("FCols_001", confTable.produceCollisions.value, initContext); - mProducedCollisionMasks = utils::enableTable("FColMasks_001", confTable.produceCollisionMasks.value, initContext); - mProducedPositions = utils::enableTable("FColPos_001", confTable.producePositions.value, initContext); - mProducedSphericities = utils::enableTable("FColSphericities_001", confTable.produceSphericities.value, initContext); - mProducedMultiplicities = utils::enableTable("FColMults_001", confTable.produceMults.value, initContext); - mProducedCentralities = utils::enableTable("FColCents_001", confTable.produceCents.value, initContext); - mProduceQns = utils::enableTable("FColQnBins_001", confTable.produceQns.value, initContext); - if (mProducedCollisions || mProducedCollisionMasks || mProducedPositions || mProducedSphericities || mProducedMultiplicities || mProducedCentralities) { - mFillAnyTable = true; - } else { - LOG(info) << "No tables configured, Selection object will not be configured..."; - LOG(info) << "Initialization done..."; - return; - } + + mMagFieldForced = confCcdb.magFieldForced.value; + mGrpPath = confCcdb.grpPath.value; + mSubGeneratorId = confFilter.subGeneratorId.value; if (!confBits.triggers.value.empty()) { mUseTrigger = true; @@ -425,9 +415,21 @@ class CollisionBuilder LOG(info) << "Init RCT flag checker with label: " << confRct.label.value << "; use ZDC: " << confRct.useZdc.value << "; Limimted acceptance is bad: " << confRct.treatLimitedAcceptanceAsBad.value; mRctFlagsChecker.init(confRct.label.value, confRct.useZdc.value, confRct.treatLimitedAcceptanceAsBad.value); } - mMagFieldForced = confCcdb.magFieldForced.value; - mGrpPath = confCcdb.grpPath.value; - mSubGeneratorId = confFilter.subGeneratorId.value; + + mProducedCollisions = utils::enableTable("FCols_001", confTable.produceCollisions.value, initContext); + mProducedCollisionMasks = utils::enableTable("FColMasks_001", confTable.produceCollisionMasks.value, initContext); + mProducedPositions = utils::enableTable("FColPos_001", confTable.producePositions.value, initContext); + mProducedSphericities = utils::enableTable("FColSphericities_001", confTable.produceSphericities.value, initContext); + mProducedMultiplicities = utils::enableTable("FColMults_001", confTable.produceMults.value, initContext); + mProducedCentralities = utils::enableTable("FColCents_001", confTable.produceCents.value, initContext); + mProduceQns = utils::enableTable("FColQnBins_001", confTable.produceQns.value, initContext); + if (mProducedCollisions || mProducedCollisionMasks || mProducedPositions || mProducedSphericities || mProducedMultiplicities || mProducedCentralities) { + mFillAnyTable = true; + } else { + LOG(info) << "No tables configured, Selection object will not be configured..."; + LOG(info) << "Initialization done..."; + return; + } mCollisionSelection.configure(registry, confFilter, confBits); mCollisionSelection.printSelections(colSelsName); diff --git a/PWGCF/Femto/Core/femtoUtils.h b/PWGCF/Femto/Core/femtoUtils.h index d7d0a61f9aa..26e7649776d 100644 --- a/PWGCF/Femto/Core/femtoUtils.h +++ b/PWGCF/Femto/Core/femtoUtils.h @@ -18,6 +18,7 @@ #include "Common/Core/TableHelper.h" +#include "CommonConstants/MathConstants.h" #include "CommonConstants/PhysicsConstants.h" #include "Framework/InitContext.h" @@ -136,14 +137,15 @@ float qn(T const& col) /// Recalculate pT for Kinks (Sigmas) using kinematic constraints inline float calcPtnew(float pxMother, float pyMother, float pzMother, float pxDaughter, float pyDaughter, float pzDaughter) { + float almost0 = 1e-6f; // Particle masses in GeV/c^2 - const float massPion = 0.13957f; - const float massNeutron = 0.93957f; - const float massSigmaMinus = 1.19745f; + auto massPion = o2::constants::physics::MassPionCharged; + auto massNeutron = o2::constants::physics::MassNeutron; + auto massSigmaMinus = o2::constants::physics::MassSigmaMinus; // Calculate mother momentum and direction versor float pMother = std::sqrt(pxMother * pxMother + pyMother * pyMother + pzMother * pzMother); - if (pMother < 1e-6f) + if (pMother < almost0) return -999.f; float versorX = pxMother / pMother; @@ -154,39 +156,39 @@ inline float calcPtnew(float pxMother, float pyMother, float pzMother, float pxD float ePi = std::sqrt(massPion * massPion + pxDaughter * pxDaughter + pyDaughter * pyDaughter + pzDaughter * pzDaughter); // Scalar product of versor with daughter momentum - float a = versorX * pxDaughter + versorY * pyDaughter + versorZ * pzDaughter; + float scalarProduct = versorX * pxDaughter + versorY * pyDaughter + versorZ * pzDaughter; // Solve quadratic equation for momentum magnitude - float K = massSigmaMinus * massSigmaMinus + massPion * massPion - massNeutron * massNeutron; - float A = 4.f * (ePi * ePi - a * a); - float B = -4.f * a * K; - float C = 4.f * ePi * ePi * massSigmaMinus * massSigmaMinus - K * K; + float k = massSigmaMinus * massSigmaMinus + massPion * massPion - massNeutron * massNeutron; + float a = 4.f * (ePi * ePi - scalarProduct * scalarProduct); + float b = -4.f * scalarProduct * k; + float c = 4.f * ePi * ePi * massSigmaMinus * massSigmaMinus - k * k; - if (std::abs(A) < 1e-6f) + if (std::abs(a) < almost0) return -999.f; - float D = B * B - 4.f * A * C; - if (D < 0.f) + float d = b * b - 4.f * a * c; + if (d < 0.f) return -999.f; - float sqrtD = std::sqrt(D); - float P1 = (-B + sqrtD) / (2.f * A); - float P2 = (-B - sqrtD) / (2.f * A); + float sqrtD = std::sqrt(d); + float p1 = (-b + sqrtD) / (2.f * a); + float p2 = (-b - sqrtD) / (2.f * a); // Pick physical solution: prefer P2 if positive, otherwise P1 - if (P2 < 0.f && P1 < 0.f) + if (p2 < 0.f && p1 < 0.f) return -999.f; - if (P2 < 0.f) - return P1; + if (p2 < 0.f) + return p1; // Choose solution closest to original momentum - float p1Diff = std::abs(P1 - pMother); - float p2Diff = std::abs(P2 - pMother); - float P = (p1Diff < p2Diff) ? P1 : P2; + float p1Diff = std::abs(p1 - pMother); + float p2Diff = std::abs(p2 - pMother); + float p = (p1Diff < p2Diff) ? p1 : p2; // Calculate pT from recalibrated momentum - float pxS = versorX * P; - float pyS = versorY * P; + float pxS = versorX * p; + float pyS = versorY * p; return std::sqrt(pxS * pxS + pyS * pyS); } diff --git a/PWGCF/Femto/Core/trackBuilder.h b/PWGCF/Femto/Core/trackBuilder.h index b3356b8f1d1..c33b5850468 100644 --- a/PWGCF/Femto/Core/trackBuilder.h +++ b/PWGCF/Femto/Core/trackBuilder.h @@ -753,8 +753,8 @@ struct TrackBuilderDerivedToDerivedProducts : o2::framework::ProducesGroup { struct ConfTrackTablesDerivedToDerived : o2::framework::ConfigurableGroup { std::string prefix = std::string("TrackTables"); - o2::framework::Configurable limitTrack1{"limitTrack1", 1, "At least this many tracks of type 1 need to be in the collision"}; - o2::framework::Configurable limitTrack2{"limitTrack2", 0, "At least this many tracks of type 2 need to be in the collision"}; + o2::framework::Configurable limitTrack1{"limitTrack1", 1, "At least this many tracks of type 1 need to be in the collision. Ignored if set to 0."}; + o2::framework::Configurable limitTrack2{"limitTrack2", 0, "At least this many tracks of type 2 need to be in the collision. Ignored if set to 0."}; }; class TrackBuilderDerivedToDerived @@ -768,6 +768,10 @@ class TrackBuilderDerivedToDerived { mLimitTrack1 = config.limitTrack1.value; mLimitTrack2 = config.limitTrack2.value; + + if (mLimitTrack1 == 0 && mLimitTrack2 == 0) { + LOG(fatal) << "Both track limits are 0. Breaking..."; + } } template @@ -775,36 +779,40 @@ class TrackBuilderDerivedToDerived { auto trackSlice1 = partitionTrack1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); auto trackSlice2 = partitionTrack2->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); - if (trackSlice1.size() >= mLimitTrack1 && trackSlice2.size() >= mLimitTrack2) { - return false; - } - return true; + return trackSlice1.size() < mLimitTrack1 || trackSlice2.size() < mLimitTrack2; } template void processTracks(T1& col, T2& /*trackTable*/, T3& partitionTrack1, T4& partitionTrack2, T5& cache, T6& newTrackTable, T7& newCollisionTable) { indexMap.clear(); - auto trackSlice1 = partitionTrack1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); - auto trackSlice2 = partitionTrack2->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); - for (auto const& track : trackSlice1) { - this->fillTrack(track, newTrackTable, newCollisionTable); + if (mLimitTrack1 > 0) { + auto trackSlice1 = partitionTrack1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + for (auto const& track : trackSlice1) { + this->fillTrack(track, newTrackTable, newCollisionTable); + } } - for (auto const& track : trackSlice2) { - this->fillTrack(track, newTrackTable, newCollisionTable); + + if (mLimitTrack2 > 0) { + auto trackSlice2 = partitionTrack2->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + for (auto const& track : trackSlice2) { + this->fillTrack(track, newTrackTable, newCollisionTable); + } } } template void fillTrack(T1 const& track, T2& trackProducts, T3& collisionProducts) { - trackProducts.producedTracks(collisionProducts.producedCollision.lastIndex(), - track.signedPt(), - track.eta(), - track.phi()); - trackProducts.producedTrackMasks(track.mask()); - indexMap.emplace(track.globalIndex(), trackProducts.producedTracks.lastIndex()); + if (indexMap.find(track.globalIndex()) == indexMap.end()) { // protect against double filling + trackProducts.producedTracks(collisionProducts.producedCollision.lastIndex(), + track.signedPt(), + track.eta(), + track.phi()); + trackProducts.producedTrackMasks(track.mask()); + indexMap.emplace(track.globalIndex(), trackProducts.producedTracks.lastIndex()); + } } template @@ -816,7 +824,6 @@ class TrackBuilderDerivedToDerived } else { this->fillTrack(daughter, trackProducts, collisionProducts); int64_t idx = trackProducts.producedTracks.lastIndex(); - indexMap.emplace(daughter.globalIndex(), idx); return idx; } } diff --git a/PWGCF/Femto/Core/v0Builder.h b/PWGCF/Femto/Core/v0Builder.h index 868e17e4629..16de0ea04cc 100644 --- a/PWGCF/Femto/Core/v0Builder.h +++ b/PWGCF/Femto/Core/v0Builder.h @@ -555,10 +555,14 @@ class V0BuilderDerivedToDerived { mLimitLambda = config.limitLambda.value; mLimitK0short = config.limitK0short.value; + + if (mLimitLambda == 0 && mLimitK0short == 0) { + LOG(fatal) << "Both lambda limit and k0short limit are 0. Breaking..."; + } } template - bool collisionHasTooFewLambdas(T1& col, T2& /*lambdaTable*/, T3& partitionLambda, T4& cache) + bool collisionHasTooFewLambdas(T1 const& col, T2 const& /*lambdaTable*/, T3& partitionLambda, T4& cache) { auto lambdaSlice = partitionLambda->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); if (lambdaSlice.size() >= mLimitLambda) { @@ -568,7 +572,7 @@ class V0BuilderDerivedToDerived } template - bool collisionHasTooFewK0shorts(T1& col, T2& /*k0shortTable*/, T3& partitionK0short, T4& cache) + bool collisionHasTooFewK0shorts(T1 const& col, T2 const& /*k0shortTable*/, T3& partitionK0short, T4& cache) { auto k0shortSlice = partitionK0short->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); if (k0shortSlice.size() >= mLimitK0short) { @@ -578,14 +582,16 @@ class V0BuilderDerivedToDerived } template - void processLambdas(T1& col, T2& /*lambdaTable*/, T3& /*oldTrackTable*/, T4& partitionLambda, T5& trackBuilder, T6& cache, T7& newLambdaTable, T8& newTrackTable, T9& newCollisionTable) + void processLambdas(T1 const& col, T2 const& /*lambdaTable*/, T3 const& oldTrackTable, T4& partitionLambda, T5& trackBuilder, T6& cache, T7& newLambdaTable, T8& newTrackTable, T9& newCollisionTable) { auto lambdaSlice = partitionLambda->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); for (auto const& lambda : lambdaSlice) { - auto posDaughter = lambda.template posDau_as(); - auto negDaughter = lambda.template negDau_as(); + // auto posDaughter = lambda.template posDau_as(); + // auto negDaughter = lambda.template negDau_as(); + auto posDaughter = oldTrackTable.rawIteratorAt(lambda.posDauId() - oldTrackTable.offset()); + auto negDaughter = oldTrackTable.rawIteratorAt(lambda.negDauId() - oldTrackTable.offset()); int posDaughterIndex = trackBuilder.getDaughterIndex(posDaughter, newTrackTable, newCollisionTable); int negDaughterIndex = trackBuilder.getDaughterIndex(negDaughter, newTrackTable, newCollisionTable); @@ -602,14 +608,16 @@ class V0BuilderDerivedToDerived } template - void processK0shorts(T1& col, T2& /*k0shortTable*/, T3& /*oldTrackTable*/, T4& partitionK0short, T5& trackBuilder, T6& cache, T7& newK0shortTable, T8& newTrackTable, T9& newCollisionTable) + void processK0shorts(T1 const& col, T2 const& /*k0shortTable*/, T3 const& oldTrackTable, T4& partitionK0short, T5& trackBuilder, T6& cache, T7& newK0shortTable, T8& newTrackTable, T9& newCollisionTable) { auto k0shortSlice = partitionK0short->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); for (auto const& k0short : k0shortSlice) { - auto posDaughter = k0short.template posDau_as(); - auto negDaughter = k0short.template negDau_as(); + // auto posDaughter = k0short.template posDau_as(); + // auto negDaughter = k0short.template negDau_as(); + auto posDaughter = oldTrackTable.rawIteratorAt(k0short.posDauId() - oldTrackTable.offset()); + auto negDaughter = oldTrackTable.rawIteratorAt(k0short.negDauId() - oldTrackTable.offset()); int posDaughterIndex = trackBuilder.getDaughterIndex(posDaughter, newTrackTable, newCollisionTable); int negDaughterIndex = trackBuilder.getDaughterIndex(negDaughter, newTrackTable, newCollisionTable); diff --git a/PWGCF/Femto/TableProducer/femtoProducerDerivedToDerived.cxx b/PWGCF/Femto/TableProducer/femtoProducerDerivedToDerived.cxx index ec684af5f62..85534fecfe3 100644 --- a/PWGCF/Femto/TableProducer/femtoProducerDerivedToDerived.cxx +++ b/PWGCF/Femto/TableProducer/femtoProducerDerivedToDerived.cxx @@ -79,7 +79,7 @@ struct FemtoProducerDerivedToDerived { trackBuilder.init(confTrackBuilder); v0Builder.init(confV0Builder); - if ((doprocessTracks + doprocessLambdas + doprocessK0shorts) > 1) { + if ((doprocessTracks + doprocessTracksLambdas + doprocessTracksK0shorts) > 1) { LOG(fatal) << "Only one proccess function can be activated"; } } @@ -95,31 +95,27 @@ struct FemtoProducerDerivedToDerived { } PROCESS_SWITCH(FemtoProducerDerivedToDerived, processTracks, "Process tracks", true); - void processLambdas(FilteredCollision const& col, Tracks const& tracks, Lambdas const& lambdas) + void processTracksLambdas(FilteredCollision const& col, Tracks const& tracks, Lambdas const& lambdas) { - if (trackBuilder.collisionHasTooFewTracks(col, tracks, trackPartition1, trackPartition2, cache) && v0Builder.collisionHasTooFewLambdas(col, lambdas, lambdaPartition, cache)) { + if (trackBuilder.collisionHasTooFewTracks(col, tracks, trackPartition1, trackPartition2, cache) || v0Builder.collisionHasTooFewLambdas(col, lambdas, lambdaPartition, cache)) { return; } - if (trackBuilder.collisionHasTooFewTracks(col, tracks, trackPartition1, trackPartition2, cache)) { - collisionBuilder.processCollision(col, collisionBuilderProducts); - trackBuilder.processTracks(col, tracks, trackPartition1, trackPartition2, cache, trackBuilderProducts, collisionBuilderProducts); - v0Builder.processLambdas(col, lambdas, tracks, lambdaPartition, trackBuilder, cache, v0BuilderProducts, trackBuilderProducts, collisionBuilderProducts); - } + collisionBuilder.processCollision(col, collisionBuilderProducts); + trackBuilder.processTracks(col, tracks, trackPartition1, trackPartition2, cache, trackBuilderProducts, collisionBuilderProducts); + v0Builder.processLambdas(col, lambdas, tracks, lambdaPartition, trackBuilder, cache, v0BuilderProducts, trackBuilderProducts, collisionBuilderProducts); } - PROCESS_SWITCH(FemtoProducerDerivedToDerived, processLambdas, "Process lambdas and tracks", false); + PROCESS_SWITCH(FemtoProducerDerivedToDerived, processTracksLambdas, "Process lambdas and tracks", false); - void processK0shorts(FilteredCollision const& col, Tracks const& tracks, K0shorts const& k0shorts) + void processTracksK0shorts(FilteredCollision const& col, Tracks const& tracks, K0shorts const& k0shorts) { - if (trackBuilder.collisionHasTooFewTracks(col, tracks, trackPartition1, trackPartition2, cache) && v0Builder.collisionHasTooFewK0shorts(col, k0shorts, k0shortPartition, cache)) { + if (trackBuilder.collisionHasTooFewTracks(col, tracks, trackPartition1, trackPartition2, cache) || v0Builder.collisionHasTooFewK0shorts(col, k0shorts, k0shortPartition, cache)) { return; } - if (trackBuilder.collisionHasTooFewTracks(col, tracks, trackPartition1, trackPartition2, cache)) { - collisionBuilder.processCollision(col, collisionBuilderProducts); - trackBuilder.processTracks(col, tracks, trackPartition1, trackPartition2, cache, trackBuilderProducts, collisionBuilderProducts); - v0Builder.processK0shorts(col, k0shorts, tracks, k0shortPartition, trackBuilder, cache, v0BuilderProducts, trackBuilderProducts, collisionBuilderProducts); - } + collisionBuilder.processCollision(col, collisionBuilderProducts); + trackBuilder.processTracks(col, tracks, trackPartition1, trackPartition2, cache, trackBuilderProducts, collisionBuilderProducts); + v0Builder.processK0shorts(col, k0shorts, tracks, k0shortPartition, trackBuilder, cache, v0BuilderProducts, trackBuilderProducts, collisionBuilderProducts); } - PROCESS_SWITCH(FemtoProducerDerivedToDerived, processK0shorts, "Process k0short and tracks", false); + PROCESS_SWITCH(FemtoProducerDerivedToDerived, processTracksK0shorts, "Process k0short and tracks", false); }; o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc)