Femtoworld

PWGCF/CMakeLists.txt

@@ -13,6 +13,7 @@ add_subdirectory(Core)
 # add_subdirectory(DataModel)
 add_subdirectory(GenericFramework)
 add_subdirectory(FemtoDream)
+add_subdirectory(FemtoWorld)
 add_subdirectory(MultiparticleCorrelations)
 add_subdirectory(Tasks)
 add_subdirectory(TableProducer)

PWGCF/FemtoWorld/CMakeLists.txt

@@ -0,0 +1,16 @@
+# Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+# All rights not expressly granted are reserved.
+#
+# This software is distributed under the terms of the GNU General Public
+# License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+#
+# In applying this license CERN does not waive the privileges and immunities
+# granted to it by virtue of its status as an Intergovernmental Organization
+# or submit itself to any jurisdiction.
+
+add_subdirectory(Core)
+add_subdirectory(DataModel)
+add_subdirectory(TableProducer)
+add_subdirectory(Tasks)
+

PWGCF/FemtoWorld/Core/CMakeLists.txt

@@ -0,0 +1,10 @@
+# Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+# All rights not expressly granted are reserved.
+#
+# This software is distributed under the terms of the GNU General Public
+# License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+#
+# In applying this license CERN does not waive the privileges and immunities
+# granted to it by virtue of its status as an Intergovernmental Organization
+# or submit itself to any jurisdiction.

PWGCF/FemtoWorld/Core/FemtoWorldContainer.h

@@ -0,0 +1,140 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file FemtoWorldContainer.h
+/// \brief Definition of the FemtoWorldContainer
+/// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de
+/// \author Valentina Mantovani Sarti, valentina.mantovani-sarti@tum.de
+
+#ifndef ANALYSIS_TASKS_PWGCF_FEMTOWORLD_FEMTOWORDLCONTAINER_H_
+#define ANALYSIS_TASKS_PWGCF_FEMTOWORLD_FEMTOWORDLCONTAINER_H_
+
+#include "Framework/HistogramRegistry.h"
+#include "FemtoWorldMath.h"
+
+#include "Math/Vector4D.h"
+#include "TMath.h"
+#include "TDatabasePDG.h"
+
+using namespace o2::framework;
+
+namespace o2::analysis::femtoWorld
+{
+
+namespace femtoWorldContainer
+{
+/// Femtoscopic observable to be computed
+enum Observable { kstar ///< kstar
+};
+
+/// Type of the event processind
+enum EventType { same, ///< Pair from same event
+                 mixed ///< Pair from mixed event
+};
+}; // namespace femtoWorldContainer
+
+/// \class FemtoWorldContainer
+/// \brief Container for all histogramming related to the correlation function. The two
+/// particles of the pair are passed here, and the correlation function and QA histograms
+/// are filled according to the specified observable
+/// \tparam eventType Type of the event (same/mixed)
+/// \tparam obs Observable to be computed (k*/Q_inv/...)
+template <femtoWorldContainer::EventType eventType, femtoWorldContainer::Observable obs>
+class FemtoWorldContainer
+{
+ public:
+  /// Destructor
+  virtual ~FemtoWorldContainer() = default;
+
+  /// Initializes histograms for the task
+  /// \tparam T Type of the configurable for the axis configuration
+  /// \param registry Histogram registry to be passed
+  /// \param kstarBins k* binning for the histograms
+  /// \param multBins multiplicity binning for the histograms
+  /// \param kTBins kT binning for the histograms
+  /// \param mTBins mT binning for the histograms
+  template <typename T>
+  void init(HistogramRegistry* registry, T& kstarBins, T& multBins, T& kTBins, T& mTBins)
+  {
+    mHistogramRegistry = registry;
+    std::string femtoObs;
+    if constexpr (mFemtoObs == femtoWorldContainer::Observable::kstar) {
+      femtoObs = "#it{k*} (GeV/#it{c})";
+    }
+    std::vector<double> tmpVecMult = multBins;
+    framework::AxisSpec multAxis = {tmpVecMult, "Multiplicity"};
+    framework::AxisSpec femtoObsAxis = {kstarBins, femtoObs.c_str()};
+    framework::AxisSpec kTAxis = {kTBins, "#it{k}_{T} (GeV/#it{c})"};
+    framework::AxisSpec mTAxis = {mTBins, "#it{m}_{T} (GeV/#it{c}^{2})"};
+
+    std::string folderName = static_cast<std::string>(mFolderSuffix[mEventType]);
+    mHistogramRegistry->add((folderName + "relPairDist").c_str(), ("; " + femtoObs + "; Entries").c_str(), kTH1F, {femtoObsAxis});
+    mHistogramRegistry->add((folderName + "relPairkT").c_str(), "; #it{k}_{T} (GeV/#it{c}); Entries", kTH1F, {kTAxis});
+    mHistogramRegistry->add((folderName + "relPairkstarkT").c_str(), ("; " + femtoObs + "; #it{k}_{T} (GeV/#it{c})").c_str(), kTH2F, {femtoObsAxis, kTAxis});
+    mHistogramRegistry->add((folderName + "relPairkstarmT").c_str(), ("; " + femtoObs + "; #it{m}_{T} (GeV/#it{c}^{2})").c_str(), kTH2F, {femtoObsAxis, mTAxis});
+    mHistogramRegistry->add((folderName + "relPairkstarMult").c_str(), ("; " + femtoObs + "; Multiplicity").c_str(), kTH2F, {femtoObsAxis, multAxis});
+    mHistogramRegistry->add((folderName + "kstarPtPart1").c_str(), ("; " + femtoObs + "; #it{p} _{T} Particle 1 (GeV/#it{c})").c_str(), kTH2F, {femtoObsAxis, {375, 0., 7.5}});
+    mHistogramRegistry->add((folderName + "kstarPtPart2").c_str(), ("; " + femtoObs + "; #it{p} _{T} Particle 2 (GeV/#it{c})").c_str(), kTH2F, {femtoObsAxis, {375, 0., 7.5}});
+    mHistogramRegistry->add((folderName + "MultPtPart1").c_str(), "; #it{p} _{T} Particle 1 (GeV/#it{c}); Multiplicity", kTH2F, {{375, 0., 7.5}, multAxis});
+    mHistogramRegistry->add((folderName + "MultPtPart2").c_str(), "; #it{p} _{T} Particle 2 (GeV/#it{c}); Multiplicity", kTH2F, {{375, 0., 7.5}, multAxis});
+    mHistogramRegistry->add((folderName + "PtPart1PtPart2").c_str(), "; #it{p} _{T} Particle 1 (GeV/#it{c}); #it{p} _{T} Particle 2 (GeV/#it{c})", kTH2F, {{375, 0., 7.5}, {375, 0., 7.5}});
+  }
+
+  /// Set the PDG codes of the two particles involved
+  /// \param pdg1 PDG code of particle one
+  /// \param pdg2 PDG code of particle two
+  void setPDGCodes(const int pdg1, const int pdg2)
+  {
+    mMassOne = TDatabasePDG::Instance()->GetParticle(pdg1)->Mass();
+    mMassTwo = TDatabasePDG::Instance()->GetParticle(pdg2)->Mass();
+  }
+
+  /// Pass a pair to the container and compute all the relevant observables
+  /// \tparam T type of the femtoworldparticle
+  /// \param part1 Particle one
+  /// \param part2 Particle two
+  /// \param mult Multiplicity of the event
+  template <typename T>
+  void setPair(T const& part1, T const& part2, const int mult)
+  {
+    float femtoObs;
+    if constexpr (mFemtoObs == femtoWorldContainer::Observable::kstar) {
+      femtoObs = FemtoWorldMath::getkstar(part1, mMassOne, part2, mMassTwo);
+    }
+    const float kT = FemtoWorldMath::getkT(part1, mMassOne, part2, mMassTwo);
+    const float mT = FemtoWorldMath::getmT(part1, mMassOne, part2, mMassTwo);
+
+    if (mHistogramRegistry) {
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairDist"), femtoObs);
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairkT"), kT);
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairkstarkT"), femtoObs, kT);
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairkstarmT"), femtoObs, mT);
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("relPairkstarMult"), femtoObs, mult);
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("kstarPtPart1"), femtoObs, part1.pt());
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("kstarPtPart2"), femtoObs, part2.pt());
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("MultPtPart1"), part1.pt(), mult);
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("MultPtPart2"), part2.pt(), mult);
+      mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST("PtPart1PtPart2"), part1.pt(), part2.pt());
+    }
+  }
+
+ protected:
+  HistogramRegistry* mHistogramRegistry = nullptr;                                    ///< For QA output
+  static constexpr std::string_view mFolderSuffix[2] = {"SameEvent/", "MixedEvent/"}; ///< Folder naming for the output according to mEventType
+  static constexpr femtoWorldContainer::Observable mFemtoObs = obs;                   ///< Femtoscopic observable to be computed (according to femtoWorldContainer::Observable)
+  static constexpr int mEventType = eventType;                                        ///< Type of the event (same/mixed, according to femtoWorldContainer::EventType)
+  float mMassOne = 0.f;                                                               ///< PDG mass of particle 1
+  float mMassTwo = 0.f;                                                               ///< PDG mass of particle 2
+};
+
+} // namespace o2::analysis::femtoWorld
+
+#endif /* ANALYSIS_TASKS_PWGCF_FEMTOWORLD_FEMTOWORDLCONTAINER_H_ */

PWGCF/FemtoWorld/Core/FemtoWorldDetaDphiStar.h

@@ -0,0 +1,203 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file FemtoWorldDetaDphiStar.h
+/// \brief FemtoWorldDetaDphiStar - Checks particles for the close pair rejection.
+/// \author Laura Serksnyte, TU München, laura.serksnyte@tum.de
+
+#ifndef ANALYSIS_TASKS_PWGCF_O2FEMTO_O2FEMTODETADPHISTAR_H_
+#define ANALYSIS_TASKS_PWGCF_O2FEMTO_O2FEMTODETADPHISTAR_H_
+
+#include "PWGCF/DataModel/FemtoDerived.h"
+
+#include "Framework/HistogramRegistry.h"
+#include <string>
+
+namespace o2::analysis
+{
+namespace femtoWorld
+{
+
+/// \class FemtoWorldDetaDphiStar
+/// \brief Class to check particles for the close pair rejection.
+/// \tparam partOne Type of particle 1 (Track/V0/Cascade/...)
+/// \tparam partTwo Type of particle 2 (Track/V0/Cascade/...)
+template <o2::aod::femtodreamparticle::ParticleType partOne, o2::aod::femtodreamparticle::ParticleType partTwo>
+class FemtoWorldDetaDphiStar
+{
+ public:
+  /// Destructor
+  virtual ~FemtoWorldDetaDphiStar() = default;
+  /// Initalization of the histograms and setting required values
+  void init(HistogramRegistry* registry, HistogramRegistry* registryQA, float ldeltaPhiMax, float ldeltaEtaMax, bool lplotForEveryRadii)
+  {
+    deltaPhiMax = ldeltaPhiMax;
+    deltaEtaMax = ldeltaEtaMax;
+    plotForEveryRadii = lplotForEveryRadii;
+    mHistogramRegistry = registry;
+    mHistogramRegistryQA = registryQA;
+
+    if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kTrack) {
+      std::string dirName = static_cast<std::string>(dirNames[0]);
+      histdetadpi[0][0] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[0][0])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
+      histdetadpi[0][1] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[1][0])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
+      if (plotForEveryRadii) {
+        for (int i = 0; i < 9; i++) {
+          histdetadpiRadii[0][i] = mHistogramRegistryQA->add<TH2>((dirName + static_cast<std::string>(histNamesRadii[0][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
+        }
+      }
+    }
+    if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kV0) {
+      for (int i = 0; i < 2; i++) {
+        std::string dirName = static_cast<std::string>(dirNames[1]);
+        histdetadpi[i][0] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[0][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
+        histdetadpi[i][1] = mHistogramRegistry->add<TH2>((dirName + static_cast<std::string>(histNames[1][i])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
+        if (plotForEveryRadii) {
+          for (int j = 0; j < 9; j++) {
+            histdetadpiRadii[i][j] = mHistogramRegistryQA->add<TH2>((dirName + static_cast<std::string>(histNamesRadii[i][j])).c_str(), "; #Delta #eta; #Delta #phi", kTH2F, {{100, -0.15, 0.15}, {100, -0.15, 0.15}});
+          }
+        }
+      }
+    }
+  }
+  ///  Check if pair is close or not
+  template <typename Part, typename Parts>
+  bool isClosePair(Part const& part1, Part const& part2, Parts const& particles, float lmagfield)
+  {
+    magfield = lmagfield;
+
+    if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kTrack) {
+      /// Track-Track combination
+      // check if provided particles are in agreement with the class instantiation
+      if (part1.partType() != o2::aod::femtodreamparticle::ParticleType::kTrack || part2.partType() != o2::aod::femtodreamparticle::ParticleType::kTrack) {
+        LOG(fatal) << "FemtoWorldDetaDphiStar: passed arguments don't agree with FemtoWorldDetaDphiStar instantiation! Please provide kTrack,kTrack candidates.";
+        return false;
+      }
+      auto deta = part1.eta() - part2.eta();
+      auto dphiAvg = AveragePhiStar(part1, part2, 0);
+      histdetadpi[0][0]->Fill(deta, dphiAvg);
+      if (pow(dphiAvg, 2) / pow(deltaPhiMax, 2) + pow(deta, 2) / pow(deltaEtaMax, 2) < 1.) {
+        return true;
+      } else {
+        histdetadpi[0][1]->Fill(deta, dphiAvg);
+        return false;
+      }
+
+    } else if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kV0) {
+      /// Track-V0 combination
+      // check if provided particles are in agreement with the class instantiation
+      if (part1.partType() != o2::aod::femtodreamparticle::ParticleType::kTrack || part2.partType() != o2::aod::femtodreamparticle::ParticleType::kV0) {
+        LOG(fatal) << "FemtoWorldDetaDphiStar: passed arguments don't agree with FemtoWorldDetaDphiStar instantiation! Please provide kTrack,kV0 candidates.";
+        return false;
+      }
+
+      bool pass = false;
+      for (int i = 0; i < 2; i++) {
+        auto indexOfDaughter = part2.index() - 2 + i;
+        auto daughter = particles.begin() + indexOfDaughter;
+        auto deta = part1.eta() - daughter.eta();
+        auto dphiAvg = AveragePhiStar(part1, *daughter, i);
+        histdetadpi[i][0]->Fill(deta, dphiAvg);
+        if (pow(dphiAvg, 2) / pow(deltaPhiMax, 2) + pow(deta, 2) / pow(deltaEtaMax, 2) < 1.) {
+          pass = true;
+        } else {
+          histdetadpi[i][1]->Fill(deta, dphiAvg);
+        }
+      }
+      return pass;
+    } else {
+      LOG(fatal) << "FemtoWorldPairCleaner: Combination of objects not defined - quitting!";
+      return false;
+    }
+  }
+
+ private:
+  HistogramRegistry* mHistogramRegistry = nullptr;   ///< For main output
+  HistogramRegistry* mHistogramRegistryQA = nullptr; ///< For QA output
+  static constexpr std::string_view dirNames[2] = {"kTrack_kTrack/", "kTrack_kV0/"};
+
+  static constexpr std::string_view histNames[2][2] = {{"detadphidetadphi0Before_0", "detadphidetadphi0Before_1"},
+                                                       {"detadphidetadphi0After_0", "detadphidetadphi0After_1"}};
+  static constexpr std::string_view histNamesRadii[2][9] = {{"detadphidetadphi0Before_0_0", "detadphidetadphi0Before_0_1", "detadphidetadphi0Before_0_2",
+                                                             "detadphidetadphi0Before_0_3", "detadphidetadphi0Before_0_4", "detadphidetadphi0Before_0_5",
+                                                             "detadphidetadphi0Before_0_6", "detadphidetadphi0Before_0_7", "detadphidetadphi0Before_0_8"},
+                                                            {"detadphidetadphi0Before_1_0", "detadphidetadphi0Before_1_1", "detadphidetadphi0Before_1_2",
+                                                             "detadphidetadphi0Before_1_3", "detadphidetadphi0Before_1_4", "detadphidetadphi0Before_1_5",
+                                                             "detadphidetadphi0Before_1_6", "detadphidetadphi0Before_1_7", "detadphidetadphi0Before_1_8"}};
+
+  static constexpr o2::aod::femtodreamparticle::ParticleType mPartOneType = partOne; ///< Type of particle 1
+  static constexpr o2::aod::femtodreamparticle::ParticleType mPartTwoType = partTwo; ///< Type of particle 2
+
+  static constexpr float tmpRadiiTPC[9] = {85., 105., 125., 145., 165., 185., 205., 225., 245.};
+
+  static constexpr uint32_t kSignMinusMask = 1;
+  static constexpr uint32_t kSignPlusMask = 1 << 1;
+  static constexpr uint32_t kValue0 = 0;
+
+  float deltaPhiMax;
+  float deltaEtaMax;
+  float magfield;
+  bool plotForEveryRadii = false;
+
+  std::array<std::array<std::shared_ptr<TH2>, 2>, 2> histdetadpi{};
+  std::array<std::array<std::shared_ptr<TH2>, 9>, 2> histdetadpiRadii{};
+
+  ///  Calculate phi at all required radii stored in tmpRadiiTPC
+  /// Magnetic field to be provided in Tesla
+  template <typename T>
+  void PhiAtRadiiTPC(const T& part, std::vector<float>& tmpVec)
+  {
+
+    float phi0 = part.phi();
+    // Start: Get the charge from cutcontainer using masks
+    float charge = 0.;
+    if ((part.cut() & kSignMinusMask) == kValue0 && (part.cut() & kSignPlusMask) == kValue0) {
+      charge = 0;
+    } else if ((part.cut() & kSignPlusMask) == kSignPlusMask) {
+      charge = 1;
+    } else if ((part.cut() & kSignMinusMask) == kSignMinusMask) {
+      charge = -1;
+    } else {
+      LOG(fatal) << "FemtoWorldDetaDphiStar: Charge bits are set wrong!";
+    }
+    // End: Get the charge from cutcontainer using masks
+    float pt = part.pt();
+    for (size_t i = 0; i < 9; i++) {
+      tmpVec.push_back(phi0 - std::asin(0.3 * charge * 0.1 * magfield * tmpRadiiTPC[i] * 0.01 / (2. * pt)));
+    }
+  }
+
+  ///  Calculate average phi
+  template <typename T1, typename T2>
+  float AveragePhiStar(const T1& part1, const T2& part2, int iHist)
+  {
+    std::vector<float> tmpVec1;
+    std::vector<float> tmpVec2;
+    PhiAtRadiiTPC(part1, tmpVec1);
+    PhiAtRadiiTPC(part2, tmpVec2);
+    const int num = tmpVec1.size();
+    float dPhiAvg = 0;
+    for (int i = 0; i < num; i++) {
+      float dphi = tmpVec1.at(i) - tmpVec2.at(i);
+      dphi = TVector2::Phi_mpi_pi(dphi);
+      dPhiAvg += dphi;
+      if (plotForEveryRadii) {
+        histdetadpiRadii[iHist][i]->Fill(part1.eta() - part2.eta(), dphi);
+      }
+    }
+    return (dPhiAvg / (float)num);
+  }
+};
+
+} /* namespace femtoWorld */
+} /* namespace o2::analysis */
+
+#endif /* ANALYSIS_TASKS_PWGCF_O2FEMTO_O2FEMTODETADPHISTAR_H_ */