[PWGCF] Fixes in femto framework

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

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

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<int> limitTrack1{"limitTrack1", 1, "At least this many tracks of type 1 need to be in the collision"};
-  o2::framework::Configurable<int> limitTrack2{"limitTrack2", 0, "At least this many tracks of type 2 need to be in the collision"};
+  o2::framework::Configurable<int> 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<int> 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,43 +768,51 @@ class TrackBuilderDerivedToDerived
   {
     mLimitTrack1 = config.limitTrack1.value;
     mLimitTrack2 = config.limitTrack2.value;
+
+    if (mLimitTrack1 == 0 && mLimitTrack2 == 0) {
+      LOG(fatal) << "Both track limits are 0. Breaking...";
+    }
   }
 
   template <typename T1, typename T2, typename T3, typename T4, typename T5>
   bool collisionHasTooFewTracks(T1& col, T2& /*trackTable*/, T3& partitionTrack1, T4& partitionTrack2, T5& cache)
   {
     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 <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7>
   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 <typename T1, typename T2, typename T3>
   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 <typename T1, typename T2, typename T3>
@@ -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;
     }
   }

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 <typename T1, typename T2, typename T3, typename T4>
-  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 <typename T1, typename T2, typename T3, typename T4>
-  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 <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
-  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<T3>();
-      auto negDaughter = lambda.template negDau_as<T3>();
+      // auto posDaughter = lambda.template posDau_as<T3>();
+      // auto negDaughter = lambda.template negDau_as<T3>();
+      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 <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
-  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<T3>();
-      auto negDaughter = k0short.template negDau_as<T3>();
+      // auto posDaughter = k0short.template posDau_as<T3>();
+      // auto negDaughter = k0short.template negDau_as<T3>();
+      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);

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)