Merge branch 'master' of https://github.com/AliceO2Group/O2Physics into doublephianalysisv12

Author: skundu692

Common/CCDB/CMakeLists.txt

@@ -19,4 +19,5 @@ o2physics_target_root_dictionary(AnalysisCCDB
               HEADERS EventSelectionParams.h
               HEADERS TriggerAliases.h
               HEADERS ctpRateFetcher.h
+              HEADERS RCTSelectionFlags.h
               LINKDEF AnalysisCCDBLinkDef.h)

PWGDQ/Core/HistogramsLibrary.cxx

@@ -1654,6 +1654,8 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
         hm->AddHistogram(histClass, "MassD0region_Rapidity_AveragePt", "", true, 140, 1.5, 2.2, VarManager::kMass, 10, -0.8, 0.8, VarManager::kRap, 150, 0.0, 30.0, VarManager::kPt);
         hm->AddHistogram(histClass, "MassD0region_Pt_ITStrackOccupancy", "Pair mass vs pair Pt vs event ITS occupancy", false, 70, 1.5, 2.2, VarManager::kMass, 160, 0., 20., VarManager::kPt, 200, 0., 20000., VarManager::kTrackOccupancyInTimeRange);
         hm->AddHistogram(histClass, "MassD0region_TPCnSigKa_pIN", "Pair mass vs kaon cand. pIN vs kaon cand. TPC n-#sigma(K)", false, 140, 1.5, 2.2, VarManager::kMass, 100, 0.0, 10.0, VarManager::kPin_leg1, 20, -5.0, 5.0, VarManager::kTPCnSigmaKa_leg1);
+        hm->AddHistogram(histClass, "Mass_Pt_Ft0cOccupancy", "", false, 150, 0.0, 5.0, VarManager::kMass, 10, 0., 10., VarManager::kPt, 20, 0., 20000., VarManager::kFT0COccupancyInTimeRange);
+        hm->AddHistogram(histClass, "Mass_Pt_ITStrackOccupancy", "", false, 150, 0.0, 5.0, VarManager::kMass, 10, 0., 10., VarManager::kPt, 20, 0., 20000., VarManager::kTrackOccupancyInTimeRange);
       }
       if (subGroupStr.Contains("lambdac")) {
         hm->AddHistogram(histClass, "MassLambdacRegion", "", false, 50, 2.15, 2.4, VarManager::kMass);

PWGDQ/Tasks/dqEfficiency_withAssoc.cxx

@@ -1252,7 +1252,7 @@ struct AnalysisSameEventPairing {
     Configurable<std::string> recSignals{"cfgBarrelMCRecSignals", "", "Comma separated list of MC signals (reconstructed)"};
     Configurable<std::string> recSignalsJSON{"cfgMCRecSignalsJSON", "", "Comma separated list of MC signals (reconstructed) via JSON"};
     Configurable<bool> skimSignalOnly{"cfgSkimSignalOnly", false, "Configurable to select only matched candidates"};
-    Configurable<bool> runMCGenPair{"cfgRunMCGenPair", false, "Do pairing of true MC particles"};
+    // Configurable<bool> runMCGenPair{"cfgRunMCGenPair", false, "Do pairing of true MC particles"};
   } fConfigMC;
 
   struct : ConfigurableGroup {
@@ -1267,6 +1267,7 @@ struct AnalysisSameEventPairing {
   Service<o2::ccdb::BasicCCDBManager> fCCDB;
 
   // Filter filterEventSelected = aod::dqanalysisflags::isEventSelected & uint32_t(1);
+  Filter eventFilter = aod::dqanalysisflags::isEventSelected > static_cast<uint32_t>(0);
 
   HistogramManager* fHistMan;
 
@@ -1298,6 +1299,7 @@ struct AnalysisSameEventPairing {
     if (context.mOptions.get<bool>("processDummy")) {
       return;
     }
+    bool isMCGen = context.mOptions.get<bool>("processMCGen") || context.mOptions.get<bool>("processMCGenWithEventSelection");
     VarManager::SetDefaultVarNames();
 
     fEnableBarrelHistos = context.mOptions.get<bool>("processAllSkimmed") || context.mOptions.get<bool>("processBarrelOnlySkimmed") || context.mOptions.get<bool>("processBarrelOnlyWithCollSkimmed");
@@ -1555,12 +1557,16 @@ struct AnalysisSameEventPairing {
       }
     }
 
-    for (auto& sig : fGenMCSignals) {
-      if (sig->GetNProngs() == 1) {
-        histNames += Form("MCTruthGen_%s;", sig->GetName()); // TODO: Add these names to a std::vector to avoid using Form in the process function
-      } else if (sig->GetNProngs() == 2) {
-        histNames += Form("MCTruthGenPair_%s;", sig->GetName());
-        fHasTwoProngGenMCsignals = true;
+    if (isMCGen) {
+      for (auto& sig : fGenMCSignals) {
+        if (sig->GetNProngs() == 1) {
+          histNames += Form("MCTruthGen_%s;", sig->GetName()); // TODO: Add these names to a std::vector to avoid using Form in the process function
+          histNames += Form("MCTruthGenSel_%s;", sig->GetName());
+        } else if (sig->GetNProngs() == 2) {
+          histNames += Form("MCTruthGenPair_%s;", sig->GetName());
+          histNames += Form("MCTruthGenPairSel_%s;", sig->GetName());
+          fHasTwoProngGenMCsignals = true;
+        }
       }
     }
 
@@ -2044,9 +2050,10 @@ struct AnalysisSameEventPairing {
   {
     runSameEventPairing<true, VarManager::kDecayToEE, gkEventFillMapWithCov, gkTrackFillMapWithCov>(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks);
     runSameEventPairing<true, VarManager::kDecayToMuMu, gkEventFillMapWithCov, gkMuonFillMapWithCov>(events, muonAssocsPerCollision, muonAssocs, muons, mcEvents, mcTracks);
-    if (fConfigMC.runMCGenPair) {
+    // Feature replaced by processMCGen and processMCGenWithEventSelection
+    /*if (fConfigMC.runMCGenPair) {
       runMCGen(mcEvents, mcTracks);
-    }
+    }*/
     // runSameEventPairing<true, VarManager::kElectronMuon, gkEventFillMap, gkTrackFillMap>(event, tracks, muons);
   }
 
@@ -2055,30 +2062,81 @@ struct AnalysisSameEventPairing {
                                 MyBarrelTracksWithCovWithAmbiguities const& barrelTracks, ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks)
   {
     runSameEventPairing<true, VarManager::kDecayToEE, gkEventFillMapWithCov, gkTrackFillMapWithCov>(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks);
-    if (fConfigMC.runMCGenPair) {
+    // Feature replaced by processMCGen and processMCGenWithEventSelection
+    /*if (fConfigMC.runMCGenPair) {
       runMCGen(mcEvents, mcTracks);
-    }
+    }*/
   }
 
   void processBarrelOnlyWithCollSkimmed(MyEventsVtxCovSelected const& events,
                                         soa::Join<aod::ReducedTracksAssoc, aod::BarrelTrackCuts, aod::Prefilter> const& barrelAssocs,
                                         MyBarrelTracksWithCovWithAmbiguitiesWithColl const& barrelTracks, ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks)
   {
     runSameEventPairing<true, VarManager::kDecayToEE, gkEventFillMapWithCov, gkTrackFillMapWithCovWithColl>(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks);
-    if (fConfigMC.runMCGenPair) {
+    // Feature replaced by processMCGen and processMCGenWithEventSelection
+    /*if (fConfigMC.runMCGenPair) {
       runMCGen(mcEvents, mcTracks);
-    }
+    }*/
   }
 
   void processMuonOnlySkimmed(MyEventsVtxCovSelected const& events,
                               soa::Join<aod::ReducedMuonsAssoc, aod::MuonTrackCuts> const& muonAssocs, MyMuonTracksWithCovWithAmbiguities const& muons, ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks)
   {
     runSameEventPairing<true, VarManager::kDecayToMuMu, gkEventFillMapWithCov, gkMuonFillMapWithCov>(events, muonAssocsPerCollision, muonAssocs, muons, mcEvents, mcTracks);
-    if (fConfigMC.runMCGenPair) {
+    // Feature replaced by processMCGen and processMCGenWithEventSelection
+    /*if (fConfigMC.runMCGenPair) {
       runMCGen(mcEvents, mcTracks);
+    }*/
+  }
+
+  void processMCGen(ReducedMCTracks const& mcTracks)
+  {
+    // loop over mc stack and fill histograms for pure MC truth signals
+    // group all the MC tracks which belong to the MC event corresponding to the current reconstructed event
+    // auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex());
+    for (auto& mctrack : mcTracks) {
+
+      VarManager::FillTrackMC(mcTracks, mctrack);
+      // NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete.
+      // NOTE: However, the working model is that the decisions on MC signals are precomputed during skimming and are stored in the mcReducedFlags member.
+      // TODO:  Use the mcReducedFlags to select signals
+      for (auto& sig : fGenMCSignals) {
+        if (sig->CheckSignal(true, mctrack)) {
+          fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), VarManager::fgValues);
+        }
+      }
     }
   }
 
+  PresliceUnsorted<ReducedMCTracks> perReducedMcGenEvent = aod::reducedtrackMC::reducedMCeventId;
+
+  void processMCGenWithEventSelection(soa::Filtered<MyEventsVtxCovSelected> const& events,
+                                      ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& mcTracks)
+  {
+    for (auto& event : events) {
+      if (!event.isEventSelected_bit(0)) {
+        continue;
+      }
+      if (!event.has_reducedMCevent()) {
+        continue;
+      }
+
+      auto groupedMCTracks = mcTracks.sliceBy(perReducedMcGenEvent, event.reducedMCeventId());
+      groupedMCTracks.bindInternalIndicesTo(&mcTracks);
+      for (auto& track : groupedMCTracks) {
+
+        VarManager::FillTrackMC(mcTracks, track);
+
+        auto track_raw = groupedMCTracks.rawIteratorAt(track.globalIndex());
+        for (auto& sig : fGenMCSignals) {
+          if (sig->CheckSignal(true, track_raw)) {
+            fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), VarManager::fgValues);
+          }
+        }
+      }
+    } // end loop over reconstructed events
+  }
+
   void processDummy(MyEvents&)
   {
     // do nothing
@@ -2088,6 +2146,8 @@ struct AnalysisSameEventPairing {
   PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlySkimmed, "Run barrel only pairing, with skimmed tracks", false);
   PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlyWithCollSkimmed, "Run barrel only pairing, with skimmed tracks and with collision information", false);
   PROCESS_SWITCH(AnalysisSameEventPairing, processMuonOnlySkimmed, "Run muon only pairing, with skimmed tracks", false);
+  PROCESS_SWITCH(AnalysisSameEventPairing, processMCGen, "Loop over MC particle stack and fill generator level histograms", false);
+  PROCESS_SWITCH(AnalysisSameEventPairing, processMCGenWithEventSelection, "Loop over MC particle stack and fill generator level histograms with  event selection", false);
   PROCESS_SWITCH(AnalysisSameEventPairing, processDummy, "Dummy function, enabled only if none of the others are enabled", false);
 };
 

PWGEM/Dilepton/Tasks/lmeeHFCocktail.cxx

@@ -160,11 +160,49 @@ struct MyConfigs : ConfigurableGroup {
 
 struct lmeehfcocktailprefilter {
 
+  HistogramRegistry registry{"registry", {}};
+  std::vector<std::shared_ptr<TH1>> hRapQuark;
   Produces<o2::aod::HfTable> hfTable;
+  ConfigurableAxis fConfigRapBins{"cfgRapBins", {200, -10.f, 10.f}, "Quark rapidity binning"};
+
+  void init(o2::framework::InitContext&)
+  {
+    const int Nchannels = 2;
+    const char* typeNamesSingle[Nchannels] = {"b", "c"};
+    const char* typeTitlesSingle[Nchannels] = {"b", "c"};
+
+    AxisSpec rap_axis = {fConfigRapBins, "y_{b}"};
+
+    // quark histograms
+    for (int i = 0; i < Nchannels; i++) {
+      hRapQuark.push_back(registry.add<TH1>(Form("Quark_Rap_%s", typeNamesSingle[i]), Form("Rap Quark %s", typeTitlesSingle[i]), HistType::kTH1F, {rap_axis}, true));
+    }
+  }
+
   void process(aod::McParticles const& mcParticles)
   {
     for (auto const& p : mcParticles) {
+      // Look at quarks which fragment
+      if (abs(p.pdgCode()) == 5 || abs(p.pdgCode()) == 4) {
+        bool foundhadrons = kFALSE;
+        if (p.has_daughters()) {
+          const auto& daughtersSlice = p.daughters_as<aod::McParticles>();
+          for (auto& d : daughtersSlice) {
+            int pdgfragment = d.pdgCode();
+            if (static_cast<int>(abs(pdgfragment) / 100.) == abs(p.pdgCode()) || static_cast<int>(abs(pdgfragment) / 1000.) == abs(p.pdgCode())) {
+              foundhadrons = kTRUE;
+            }
+          }
+        }
+        if (foundhadrons) {
+          if (abs(p.pdgCode()) == 4)
+            hRapQuark[1]->Fill(p.y());
+          else if (abs(p.pdgCode()) == 5)
+            hRapQuark[0]->Fill(p.y());
+        }
+      }
 
+      // Look at electrons
       if (abs(p.pdgCode()) != 11 || o2::mcgenstatus::getHepMCStatusCode(p.statusCode()) != 1 || !p.has_mothers()) {
         hfTable(EFromHFType::kNoE, -1, -1, -1, -1, -1, -1, -999., -999.);
         continue;

PWGEM/Dilepton/Utils/MCUtilities.h

@@ -481,14 +481,16 @@ int searchMothers(T& p, U& mcParticles, int pdg, bool equal)
     for (int i : allmothersids) {
       auto mother = mcParticles.iteratorAt(i);
       int mpdg = mother.pdgCode();
-      if (abs(mpdg) == pdg && mpdg * p.pdgCode() > 0) { // check for quark
-        if (quark_id > -1 || next_mother_id > -1) {     // we already found a possible candidate in the list of mothers, so now we have (at least) two
+      // if (abs(mpdg) == pdg && mpdg * p.pdgCode() > 0) { // check for quark
+      if (abs(mpdg) == pdg) {                       // check for quark to allow for beauty and charm + oscillation
+        if (quark_id > -1 || next_mother_id > -1) { // we already found a possible candidate in the list of mothers, so now we have (at least) two
           // LOG(warning) << "Flavour tracking is ambiguous. Stopping here.";
           return -1;
         }
         quark_id = i;
-      } else if ((static_cast<int>(abs(mpdg) / 100) == pdg || static_cast<int>(abs(mpdg) / 1000) == pdg) && mpdg * p.pdgCode() > 0) { // check for other mothers with flavour content
-        if (quark_id > -1 || next_mother_id > -1) {                                                                                   // we already found a possible candidate in the list of mothers, so now we have (at least) two
+        //} else if ((static_cast<int>(abs(mpdg) / 100) == pdg || static_cast<int>(abs(mpdg) / 1000) == pdg) && mpdg * p.pdgCode() > 0) { // check for other mothers with flavour content
+      } else if ((static_cast<int>(abs(mpdg) / 100) == pdg || static_cast<int>(abs(mpdg) / 1000) == pdg)) { // check for other mothers with flavour content to allow for beauty and charm
+        if (quark_id > -1 || next_mother_id > -1) {                                                         // we already found a possible candidate in the list of mothers, so now we have (at least) two
           // LOG(warning) << "Flavour tracking is ambiguous. Stopping here.";
           return -1;
         }

PWGEM/PhotonMeson/DataModel/bcWiseTables.h

@@ -79,29 +79,29 @@ DECLARE_SOA_COLUMN(StoredCentrality, storedCentrality, uint8_t); //! FT0M centra
 DECLARE_SOA_COLUMN(StoredZVtx, storedZVtx, int16_t);             //! Z-vertex position (x1000)
 
 DECLARE_SOA_DYNAMIC_COLUMN(Centrality, centrality, [](uint8_t storedcentrality) -> float { return storedcentrality / emdownscaling::downscalingFactors[emdownscaling::kFT0MCent]; }); //! Centrality (0-100)
-DECLARE_SOA_DYNAMIC_COLUMN(ZVtx, zVtx, [](uint8_t storedzvtx) -> float { return storedzvtx / emdownscaling::downscalingFactors[emdownscaling::kZVtx]; });                             //! Centrality (0-100)
+DECLARE_SOA_DYNAMIC_COLUMN(ZVtx, zVtx, [](int16_t storedzvtx) -> float { return storedzvtx / emdownscaling::downscalingFactors[emdownscaling::kZVtx]; });                             //! Centrality (0-100)
 } // namespace bcwisecollision
 DECLARE_SOA_TABLE(BCWiseCollisions, "AOD", "BCWISECOLL", //! table of skimmed EMCal clusters
                   o2::soa::Index<>, BCWiseBCId, bcwisecollision::StoredCentrality, bcwisecollision::StoredZVtx,
                   bcwisecollision::Centrality<bcwisecollision::StoredCentrality>, bcwisecollision::ZVtx<bcwisecollision::StoredZVtx>);
 
 namespace bcwisecluster
 {
-DECLARE_SOA_COLUMN(StoredDefinition, storedDefinition, uint8_t); //! cluster definition, see EMCALClusterDefinition.h
-DECLARE_SOA_COLUMN(StoredE, storedE, uint16_t);                  //! cluster energy (1 MeV -> Maximum cluster energy of ~65 GeV)
-DECLARE_SOA_COLUMN(StoredEta, storedEta, int16_t);               //! cluster pseudorapidity (x10,000)
-DECLARE_SOA_COLUMN(StoredPhi, storedPhi, uint16_t);              //! cluster azimuthal angle (x10 000) from 0 to 2pi
-DECLARE_SOA_COLUMN(StoredNCells, storedNCells, uint8_t);         //! number of cells in cluster
-DECLARE_SOA_COLUMN(StoredM02, storedM02, uint16_t);              //! shower shape long axis (x10 000)
-DECLARE_SOA_COLUMN(StoredTime, storedTime, int16_t);             //! cluster time (10 ps resolution)
-DECLARE_SOA_COLUMN(StoredIsExotic, storedIsExotic, bool);        //! flag to mark cluster as exotic
-
-DECLARE_SOA_DYNAMIC_COLUMN(Definition, definition, [](uint8_t storedDefinition) -> uint8_t { return storedDefinition; });                                 //! cluster definition, see EMCALClusterDefinition.h
-DECLARE_SOA_DYNAMIC_COLUMN(E, e, [](uint16_t storedE) -> float { return storedE / emdownscaling::downscalingFactors[emdownscaling::kEnergy]; });          //! cluster energy (GeV)
+DECLARE_SOA_COLUMN(StoredDefinition, storedDefinition, int8_t); //! cluster definition, see EMCALClusterDefinition.h
+DECLARE_SOA_COLUMN(StoredE, storedE, int16_t);                  //! cluster energy (1 MeV -> Maximum cluster energy of ~32 GeV)
+DECLARE_SOA_COLUMN(StoredEta, storedEta, int16_t);              //! cluster pseudorapidity (x10,000)
+DECLARE_SOA_COLUMN(StoredPhi, storedPhi, uint16_t);             //! cluster azimuthal angle (x10 000) from 0 to 2pi
+DECLARE_SOA_COLUMN(StoredNCells, storedNCells, int8_t);         //! number of cells in cluster
+DECLARE_SOA_COLUMN(StoredM02, storedM02, int16_t);              //! shower shape long axis (x10 000)
+DECLARE_SOA_COLUMN(StoredTime, storedTime, int16_t);            //! cluster time (10 ps resolution)
+DECLARE_SOA_COLUMN(StoredIsExotic, storedIsExotic, bool);       //! flag to mark cluster as exotic
+
+DECLARE_SOA_DYNAMIC_COLUMN(Definition, definition, [](int8_t storedDefinition) -> int8_t { return storedDefinition; });                                   //! cluster definition, see EMCALClusterDefinition.h
+DECLARE_SOA_DYNAMIC_COLUMN(E, e, [](int16_t storedE) -> float { return storedE / emdownscaling::downscalingFactors[emdownscaling::kEnergy]; });           //! cluster energy (GeV)
 DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, [](int16_t storedEta) -> float { return storedEta / emdownscaling::downscalingFactors[emdownscaling::kEta]; });      //! cluster pseudorapidity
 DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, [](uint16_t storedPhi) -> float { return storedPhi / emdownscaling::downscalingFactors[emdownscaling::kPhi]; });     //! cluster azimuthal angle (0 to 2pi)
-DECLARE_SOA_DYNAMIC_COLUMN(NCells, nCells, [](uint16_t storedNCells) -> uint16_t { return storedNCells; });                                               //! number of cells in cluster
-DECLARE_SOA_DYNAMIC_COLUMN(M02, m02, [](uint16_t storedM02) -> float { return storedM02 / emdownscaling::downscalingFactors[emdownscaling::kM02]; });     //! shower shape long axis
+DECLARE_SOA_DYNAMIC_COLUMN(NCells, nCells, [](int16_t storedNCells) -> int16_t { return storedNCells; });                                                 //! number of cells in cluster
+DECLARE_SOA_DYNAMIC_COLUMN(M02, m02, [](int16_t storedM02) -> float { return storedM02 / emdownscaling::downscalingFactors[emdownscaling::kM02]; });      //! shower shape long axis
 DECLARE_SOA_DYNAMIC_COLUMN(Time, time, [](int16_t storedTime) -> float { return storedTime / emdownscaling::downscalingFactors[emdownscaling::kTime]; }); //! cluster time (ns)
 DECLARE_SOA_DYNAMIC_COLUMN(IsExotic, isExotic, [](bool storedIsExotic) -> bool { return storedIsExotic; });                                               //! flag to mark cluster as exotic