Add Neutrons within the ZDC acceptance in the strangeness data model (for UPC analysis)

Author: romainschotter

PWGLF/DataModel/LFStrangenessTables.h

@@ -1725,6 +1725,38 @@ DECLARE_SOA_TABLE(Tracked3BodyColls, "AOD", "TRA3BODYCOLL", //! Table joinable w
 using Tracked3BodyColl = Tracked3BodyColls::iterator;
 using AssignedTracked3Bodys = soa::Join<aod::Tracked3Bodys, aod::Tracked3BodyColls>;
 using AssignedTracked3Body = AssignedTracked3Bodys::iterator;
+
+namespace zdcneutrons
+{
+// FOR DERIVED
+DECLARE_SOA_INDEX_COLUMN(StraMCCollision, straMCCollision); //!
+// DYNAMIC COLUMNS
+DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, //! neutron transverse momentum (GeV/c)
+                           [](float px, float py) -> float { return RecoDecay::sqrtSumOfSquares(px, py); });
+DECLARE_SOA_DYNAMIC_COLUMN(P, p, //! neutron total momentum (GeV/c)
+                           [](float px, float py, float pz) -> float { return RecoDecay::sqrtSumOfSquares(px, py, pz); });
+DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, //! neutron phi in the range [0, 2pi)
+                           [](float px, float py) -> float { return RecoDecay::phi(px, py); });
+DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, //! neutron pseudorapidity
+                           [](float px, float py, float pz) -> float { return RecoDecay::eta(std::array{px, py, pz}); });
+DECLARE_SOA_DYNAMIC_COLUMN(Y, y, //! neutron rapidity
+                           [](float pz, float e) -> float { return std::atanh(pz/e); });
+} // namespace zdcneutrons
+
+DECLARE_SOA_TABLE(ZDCNeutrons, "AOD", "ZDCNEUTRON", //! MC properties of the neutrons within ZDC acceptance (for UPC analysis)
+                  mcparticle::PdgCode, mcparticle::StatusCode, mcparticle::Flags, 
+                  mcparticle::Vx, mcparticle::Vy, mcparticle::Vz, mcparticle::Vt,
+                  mcparticle::Px, mcparticle::Py, mcparticle::Pz, mcparticle::E,
+                  // Dynamic columns for manipulating information
+                  zdcneutrons::Pt<mcparticle::Px, mcparticle::Py>,
+                  zdcneutrons::P<mcparticle::Px, mcparticle::Py, mcparticle::Pz>,
+                  zdcneutrons::Phi<mcparticle::Px, mcparticle::Py>,
+                  zdcneutrons::Eta<mcparticle::Px, mcparticle::Py, mcparticle::Pz>,
+                  zdcneutrons::Y<mcparticle::Pz, mcparticle::E>);
+
+DECLARE_SOA_TABLE(ZDCNeutronMCCollRefs, "AOD", "ZDCNEUTRONMCCOLLREF", //! refers MC candidate back to proper MC Collision
+                  o2::soa::Index<>, zdcneutrons::StraMCCollisionId, o2::soa::Marker<4>);
+
 } // namespace o2::aod
 
 //______________________________________________________

PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx

@@ -111,6 +111,11 @@ struct strangederivedbuilder {
   Produces<aod::CascMCMothers> cascmothers;   // casc mother references
   Produces<aod::MotherMCParts> motherMCParts; // mc particles for mothers
 
+  //__________________________________________________
+  // UPC specific information
+  Produces<aod::ZDCNeutrons> zdcNeutrons;                    // Primary neutrons within ZDC acceptance
+  Produces<aod::ZDCNeutronMCCollRefs> zdcNeutronsMCCollRefs; // references collisions from ZDCNeutrons
+
   //__________________________________________________
   // Q-vectors
   Produces<aod::StraFT0AQVs> StraFT0AQVs;     // FT0A Q-vector
@@ -1053,6 +1058,30 @@ struct strangederivedbuilder {
     straOrigin(origin.dataframeID());
   }
 
+  void processSimulatedZDCNeutrons(soa::Join<aod::McCollisions, aod::McCollsExtra, aod::MultsExtraMC> const& mcCollisions, aod::McParticles const& mcParticlesEntireTable)
+  {
+    for (const auto& mccollision : mcCollisions) {
+      const uint64_t mcCollIndex = mccollision.globalIndex();
+      auto mcParticles = mcParticlesEntireTable.sliceBy(mcParticlePerMcCollision, mcCollIndex);
+
+      for (const auto& mcPart : mcParticles) {
+        if (!mcPart.isPhysicalPrimary()) {
+          continue;
+        }
+
+        if (std::abs(mcPart.pdgCode()) == kNeutron) { // check if it is a neutron or anti-neutron
+          if (std::abs(mcPart.eta()) > 8.7) { // check if it is within ZDC acceptance
+            zdcNeutrons(mcPart.pdgCode(), mcPart.statusCode(), mcPart.flags(), 
+                        mcPart.vx(), mcPart.vy(), mcPart.vz(), mcPart.vt(),
+                        mcPart.px(), mcPart.py(), mcPart.pz(), mcPart.e());
+
+            zdcNeutronsMCCollRefs(mcPart.mcCollisionId());
+          }
+        }
+      }
+    }
+  }
+
   // debug processing
   PROCESS_SWITCH(strangederivedbuilder, processDataframeIDs, "Produce data frame ID tags", false);
 
@@ -1076,6 +1105,7 @@ struct strangederivedbuilder {
   PROCESS_SWITCH(strangederivedbuilder, processPureSimulation, "Produce pure simulated information", true);
   PROCESS_SWITCH(strangederivedbuilder, processReconstructedSimulation, "Produce reco-ed simulated information", true);
   PROCESS_SWITCH(strangederivedbuilder, processBinnedGenerated, "Produce binned generated information", false);
+  PROCESS_SWITCH(strangederivedbuilder, processSimulatedZDCNeutrons, "Produce generated neutrons (within ZDC acceptance) table for UPC analysis", false);
 
   // event plane information
   PROCESS_SWITCH(strangederivedbuilder, processFT0AQVectors, "Produce FT0A Q-vectors table", false);