|
36 | 36 | #include <Framework/runDataProcessing.h> |
37 | 37 | #include <ReconstructionDataFormats/PID.h> |
38 | 38 |
|
| 39 | +#include "TDatabasePDG.h" |
39 | 40 | #include <TComplex.h> |
40 | 41 | #include <TF1.h> |
41 | 42 | #include <TH2.h> |
42 | 43 | #include <THn.h> |
43 | 44 | #include <TList.h> |
44 | 45 | #include <TMath.h> |
| 46 | +#include <TPDGCode.h> |
45 | 47 | #include <TProfile.h> |
46 | 48 | #include <TProfile2D.h> |
47 | 49 | #include <TRandom3.h> |
@@ -101,6 +103,7 @@ struct V0ptHadPiKaProt { |
101 | 103 | kFV0A |
102 | 104 | }; |
103 | 105 |
|
| 106 | + Configurable<bool> cfgIsMC{"cfgIsMC", false, "Run MC"}; |
104 | 107 | Configurable<float> cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; |
105 | 108 | Configurable<float> cfgCutTpcChi2NCl{"cfgCutTpcChi2NCl", 2.5f, "Maximum TPCchi2NCl"}; |
106 | 109 | Configurable<float> cfgCutItsChi2NCl{"cfgCutItsChi2NCl", 36.0f, "Maximum ITSchi2NCl"}; |
@@ -221,6 +224,10 @@ struct V0ptHadPiKaProt { |
221 | 224 | using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::CentFV0As, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFDDMs, aod::Mults>>; |
222 | 225 | using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::pidTPCFullPr, aod::pidTOFFullPr, aod::pidTPCFullKa, aod::pidTOFFullKa, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::pidTPCFullEl, aod::pidTOFFullEl>>; |
223 | 226 |
|
| 227 | + using MyMCRecCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::CentFV0As, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::Mults, aod::McCollisionLabels>>; |
| 228 | + using MyMCTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::McTrackLabels, aod::pidTPCFullPr, aod::pidTOFFullPr, aod::pidTPCFullKa, aod::pidTOFFullKa, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::pidTPCFullEl, aod::pidTOFFullEl>>; |
| 229 | + Preslice<aod::McParticles> perMcCollision = aod::mcparticle::mcCollisionId; |
| 230 | + |
224 | 231 | std::array<float, 6> tofNsigmaCut; |
225 | 232 | std::array<float, 6> itsNsigmaCut; |
226 | 233 | std::array<float, 6> tpcNsigmaCut; |
@@ -463,6 +470,30 @@ struct V0ptHadPiKaProt { |
463 | 470 | fPtDepDCAz = new TF1("ptDepDCAz", Form("%s", cfgDCAzFunc->c_str()), 0.001, 1000); |
464 | 471 | } |
465 | 472 |
|
| 473 | + if (cfgIsMC) { |
| 474 | + // MC event counts |
| 475 | + histos.add("MCGenerated/hMC", "MC Event statistics", kTH1F, {{10, 0.0f, 10.0f}}); |
| 476 | + histos.add("MCGenerated/hPtEtaPhiCharged", "MC charged particles' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 477 | + histos.add("MCGenerated/hPtEtaPhiPion", "MC charged pions' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 478 | + histos.add("MCGenerated/hPtEtaPhiKaon", "MC charged kaons' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 479 | + histos.add("MCGenerated/hPtEtaPhiProton", "MC charged protons' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 480 | + |
| 481 | + histos.add("MCReconstructed/hPtEtaPhiChargedParticle", "MC reconstructed charged particles' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 482 | + histos.add("MCReconstructed/hPtEtaPhiChargedTrack", "MC reconstructed charged tracks' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 483 | + |
| 484 | + histos.add("MCReconstructed/hPtEtaPhiPionParticle", "MC reconstructed pion particles' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 485 | + histos.add("MCReconstructed/hPtEtaPhiPionTrack", "MC reconstructed pion tracks' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 486 | + histos.add("MCReconstructed/hPtEtaPhiTruePionTrack", "MC reconstructed pdgcode matched pion tracks' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 487 | + |
| 488 | + histos.add("MCReconstructed/hPtEtaPhiKaonParticle", "MC reconstructed kaon particles' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 489 | + histos.add("MCReconstructed/hPtEtaPhiKaonTrack", "MC reconstructed kaon tracks' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 490 | + histos.add("MCReconstructed/hPtEtaPhiTrueKaonTrack", "MC reconstructed pdgcode matched kaon tracks' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 491 | + |
| 492 | + histos.add("MCReconstructed/hPtEtaPhiProtonParticle", "MC reconstructed proton particles' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 493 | + histos.add("MCReconstructed/hPtEtaPhiProtonTrack", "MC reconstructed proton tracks' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 494 | + histos.add("MCReconstructed/hPtEtaPhiTrueProtonTrack", "MC reconstructed pdgcode matched proton tracks' pt, eta, phi", kTH3D, {ptAxis, {100, 0., o2::constants::math::TwoPI}, {100, -2.01, 2.01}}); |
| 495 | + } |
| 496 | + |
466 | 497 | } // end init |
467 | 498 |
|
468 | 499 | //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
@@ -845,8 +876,204 @@ struct V0ptHadPiKaProt { |
845 | 876 | return ptweight; |
846 | 877 | } |
847 | 878 |
|
| 879 | + // process MC recosnstructed data |
| 880 | + void processMCRec(MyMCRecCollisions::iterator const& collision, MyMCTracks const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles) |
| 881 | + { |
| 882 | + histos.fill(HIST("MCGenerated/hMC"), 5.5); |
| 883 | + |
| 884 | + if (!collision.has_mcCollision()) { |
| 885 | + return; |
| 886 | + } |
| 887 | + histos.fill(HIST("MCGenerated/hMC"), 6.5); |
| 888 | + |
| 889 | + if (!eventSelectionDefaultCuts(collision)) { |
| 890 | + return; |
| 891 | + } |
| 892 | + histos.fill(HIST("MCGenerated/hMC"), 7.5); |
| 893 | + |
| 894 | + fillMultCorrPlotsBeforeSel(collision, tracks); |
| 895 | + |
| 896 | + const auto centralityFT0C = collision.centFT0C(); |
| 897 | + if (cfgUseSmallIonAdditionalEventCut && !eventSelectedSmallion(collision, tracks.size(), centralityFT0C)) |
| 898 | + return; |
| 899 | + |
| 900 | + if (cfgUseSmallIonAdditionalEventCut) { |
| 901 | + fillMultCorrPlotsAfterSel(collision, tracks); |
| 902 | + } |
| 903 | + |
| 904 | + histos.fill(HIST("MCGenerated/hMC"), 8.5); |
| 905 | + |
| 906 | + // Centrality |
| 907 | + double cent = 0.0; |
| 908 | + if (cfgCentralityChoice == kFT0C) |
| 909 | + cent = collision.centFT0C(); |
| 910 | + else if (cfgCentralityChoice == kFT0A) |
| 911 | + cent = collision.centFT0A(); |
| 912 | + else if (cfgCentralityChoice == kFT0M) |
| 913 | + cent = collision.centFT0M(); |
| 914 | + else if (cfgCentralityChoice == kFV0A) |
| 915 | + cent = collision.centFV0A(); |
| 916 | + |
| 917 | + histos.fill(HIST("hZvtx_after_sel"), collision.posZ()); |
| 918 | + histos.fill(HIST("hCentrality"), cent); |
| 919 | + histos.fill(HIST("Hist2D_globalTracks_PVTracks"), collision.multNTracksPV(), tracks.size()); |
| 920 | + histos.fill(HIST("Hist2D_cent_nch"), tracks.size(), centralityFT0C); |
| 921 | + |
| 922 | + // Calculating generated level pt distributions |
| 923 | + auto mcColl = collision.mcCollision(); |
| 924 | + // Slice particles belonging only to this MC collision |
| 925 | + auto particlesThisEvent = mcParticles.sliceBy(perMcCollision, mcColl.globalIndex()); |
| 926 | + |
| 927 | + for (const auto& mcParticle : particlesThisEvent) { |
| 928 | + if (!mcParticle.has_mcCollision()) |
| 929 | + continue; |
| 930 | + |
| 931 | + // charged check |
| 932 | + auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcParticle.pdgCode()); |
| 933 | + if (!pdgEntry) |
| 934 | + continue; |
| 935 | + if (pdgEntry->Charge() == 0) |
| 936 | + continue; |
| 937 | + |
| 938 | + if (mcParticle.isPhysicalPrimary()) { |
| 939 | + if ((mcParticle.pt() > cfgCutPtLower) && (mcParticle.pt() < cfgCutPtUpper) && (std::abs(mcParticle.eta()) < cfgCutEta)) { |
| 940 | + histos.fill(HIST("MCGenerated/hPtEtaPhiCharged"), mcParticle.pt(), mcParticle.eta(), mcParticle.phi()); |
| 941 | + |
| 942 | + auto pdgcode = std::abs(mcParticle.pdgCode()); |
| 943 | + |
| 944 | + if (pdgcode == PDG_t::kPiPlus) |
| 945 | + histos.fill(HIST("MCGenerated/hPtEtaPhiPion"), mcParticle.pt(), mcParticle.eta(), mcParticle.phi()); |
| 946 | + |
| 947 | + if (pdgcode == PDG_t::kKPlus) |
| 948 | + histos.fill(HIST("MCGenerated/hPtEtaPhiKaon"), mcParticle.pt(), mcParticle.eta(), mcParticle.phi()); |
| 949 | + |
| 950 | + if (pdgcode == PDG_t::kProton) |
| 951 | + histos.fill(HIST("MCGenerated/hPtEtaPhiProton"), mcParticle.pt(), mcParticle.eta(), mcParticle.phi()); |
| 952 | + } |
| 953 | + } |
| 954 | + } //! end mc particle loop |
| 955 | + |
| 956 | + for (const auto& track : tracks) { // Loop over tracks |
| 957 | + |
| 958 | + if (!track.has_collision()) { |
| 959 | + continue; |
| 960 | + } |
| 961 | + if (!track.has_mcParticle()) { //! check if track has corresponding MC particle |
| 962 | + continue; |
| 963 | + } |
| 964 | + |
| 965 | + auto particle = track.mcParticle(); |
| 966 | + if (!particle.has_mcCollision()) { |
| 967 | + continue; |
| 968 | + } |
| 969 | + |
| 970 | + if (!track.isPVContributor()) { |
| 971 | + continue; |
| 972 | + } |
| 973 | + |
| 974 | + if (!(track.itsNCls() > cfgITScluster) || !(track.tpcNClsFound() >= cfgTPCcluster) || !(track.tpcNClsCrossedRows() >= cfgTPCnCrossedRows)) { |
| 975 | + continue; |
| 976 | + } |
| 977 | + |
| 978 | + if (cfgUsePtDepDCAxy && !(std::abs(track.dcaXY()) < fPtDepDCAxy->Eval(track.pt()))) { |
| 979 | + continue; |
| 980 | + } |
| 981 | + if (cfgUsePtDepDCAz && !(std::abs(track.dcaZ()) < fPtDepDCAz->Eval(track.pt()))) { |
| 982 | + continue; |
| 983 | + } |
| 984 | + |
| 985 | + if (track.sign() == 0) |
| 986 | + continue; |
| 987 | + |
| 988 | + if (particle.isPhysicalPrimary()) { |
| 989 | + if ((particle.pt() > cfgCutPtLower) && (particle.pt() < cfgCutPtUpper) && (std::abs(particle.eta()) < cfgCutEta)) { |
| 990 | + |
| 991 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiChargedParticle"), particle.pt(), particle.eta(), particle.phi()); |
| 992 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiChargedTrack"), track.pt(), track.eta(), track.phi()); |
| 993 | + |
| 994 | + // PID QAs before selection |
| 995 | + double nSigmaTpcPi = track.tpcNSigmaPi(); |
| 996 | + double nSigmaTpcKa = track.tpcNSigmaKa(); |
| 997 | + double nSigmaTpcProt = track.tpcNSigmaPr(); |
| 998 | + double nSigmaTofPi = track.tofNSigmaPi(); |
| 999 | + double nSigmaTofKa = track.tofNSigmaKa(); |
| 1000 | + double nSigmaTofProt = track.tofNSigmaPr(); |
| 1001 | + histos.fill(HIST("h2DnsigmaPionTpcVsPtBeforeCut"), track.pt(), nSigmaTpcPi); |
| 1002 | + histos.fill(HIST("h2DnsigmaKaonTpcVsPtBeforeCut"), track.pt(), nSigmaTpcKa); |
| 1003 | + histos.fill(HIST("h2DnsigmaProtonTpcVsPtBeforeCut"), track.pt(), nSigmaTpcProt); |
| 1004 | + histos.fill(HIST("h2DnsigmaPionTofVsPtBeforeCut"), track.pt(), nSigmaTofPi); |
| 1005 | + histos.fill(HIST("h2DnsigmaKaonTofVsPtBeforeCut"), track.pt(), nSigmaTofKa); |
| 1006 | + histos.fill(HIST("h2DnsigmaProtonTofVsPtBeforeCut"), track.pt(), nSigmaTofProt); |
| 1007 | + histos.fill(HIST("h2DnsigmaPionTpcVsTofBeforeCut"), nSigmaTpcPi, nSigmaTofPi); |
| 1008 | + histos.fill(HIST("h2DnsigmaKaonTpcVsTofBeforeCut"), nSigmaTpcKa, nSigmaTofKa); |
| 1009 | + histos.fill(HIST("h2DnsigmaProtonTpcVsTofBeforeCut"), nSigmaTpcProt, nSigmaTofProt); |
| 1010 | + |
| 1011 | + // identified particles selection |
| 1012 | + bool isPion = false; |
| 1013 | + bool isKaon = false; |
| 1014 | + bool isProton = false; |
| 1015 | + |
| 1016 | + if (cfgUseRun3V2PID) { |
| 1017 | + int pidVal = getNsigmaPID(track); |
| 1018 | + if (pidVal == PIONS + 1) |
| 1019 | + isPion = true; |
| 1020 | + if (pidVal == KAONS + 1) |
| 1021 | + isKaon = true; |
| 1022 | + if (pidVal == PROTONS + 1) |
| 1023 | + isProton = true; |
| 1024 | + } else { |
| 1025 | + isPion = selectionPion(track); |
| 1026 | + isKaon = selectionKaon(track); |
| 1027 | + isProton = selectionProton(track); |
| 1028 | + } |
| 1029 | + |
| 1030 | + // PID QAs after selection |
| 1031 | + if (isPion) { |
| 1032 | + histos.fill(HIST("h2DnsigmaPionTpcVsPtAfterCut"), track.pt(), nSigmaTpcPi); |
| 1033 | + histos.fill(HIST("h2DnsigmaPionTofVsPtAfterCut"), track.pt(), nSigmaTofPi); |
| 1034 | + histos.fill(HIST("h2DnsigmaPionTpcVsTofAfterCut"), nSigmaTpcPi, nSigmaTofPi); |
| 1035 | + } |
| 1036 | + if (isKaon) { |
| 1037 | + histos.fill(HIST("h2DnsigmaKaonTpcVsPtAfterCut"), track.pt(), nSigmaTpcKa); |
| 1038 | + histos.fill(HIST("h2DnsigmaKaonTofVsPtAfterCut"), track.pt(), nSigmaTofKa); |
| 1039 | + histos.fill(HIST("h2DnsigmaKaonTpcVsTofAfterCut"), nSigmaTpcKa, nSigmaTofKa); |
| 1040 | + } |
| 1041 | + if (isProton) { |
| 1042 | + histos.fill(HIST("h2DnsigmaProtonTpcVsPtAfterCut"), track.pt(), nSigmaTpcProt); |
| 1043 | + histos.fill(HIST("h2DnsigmaProtonTofVsPtAfterCut"), track.pt(), nSigmaTofProt); |
| 1044 | + histos.fill(HIST("h2DnsigmaProtonTpcVsTofAfterCut"), nSigmaTpcProt, nSigmaTofProt); |
| 1045 | + } |
| 1046 | + |
| 1047 | + auto pdgcodeRec = std::abs(particle.pdgCode()); |
| 1048 | + |
| 1049 | + if (isPion) { |
| 1050 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiPionParticle"), particle.pt(), particle.eta(), particle.phi()); |
| 1051 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiPionTrack"), track.pt(), track.eta(), track.phi()); |
| 1052 | + if (pdgcodeRec == PDG_t::kPiPlus) |
| 1053 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiTruePionTrack"), track.pt(), track.eta(), track.phi()); |
| 1054 | + } |
| 1055 | + |
| 1056 | + if (isKaon) { |
| 1057 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiKaonParticle"), particle.pt(), particle.eta(), particle.phi()); |
| 1058 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiKaonTrack"), track.pt(), track.eta(), track.phi()); |
| 1059 | + if (pdgcodeRec == PDG_t::kKPlus) |
| 1060 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiTrueKaonTrack"), track.pt(), track.eta(), track.phi()); |
| 1061 | + } |
| 1062 | + |
| 1063 | + if (isProton) { |
| 1064 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiProtonParticle"), particle.pt(), particle.eta(), particle.phi()); |
| 1065 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiProtonTrack"), track.pt(), track.eta(), track.phi()); |
| 1066 | + if (pdgcodeRec == PDG_t::kProton) |
| 1067 | + histos.fill(HIST("MCReconstructed/hPtEtaPhiTrueProtonTrack"), track.pt(), track.eta(), track.phi()); |
| 1068 | + } |
| 1069 | + } |
| 1070 | + } |
| 1071 | + } // end track loop |
| 1072 | + } |
| 1073 | + PROCESS_SWITCH(V0ptHadPiKaProt, processMCRec, "Process Monte-carlo data", false); |
| 1074 | + |
848 | 1075 | // process Data |
849 | | - void process(AodCollisions::iterator const& coll, aod::BCsWithTimestamps const&, AodTracks const& inputTracks) |
| 1076 | + void processData(AodCollisions::iterator const& coll, aod::BCsWithTimestamps const&, AodTracks const& inputTracks) |
850 | 1077 | { |
851 | 1078 | if (!eventSelectionDefaultCuts(coll)) { |
852 | 1079 | return; |
@@ -1212,6 +1439,7 @@ struct V0ptHadPiKaProt { |
1212 | 1439 | fPtProfileProtInWinB->Delete(); |
1213 | 1440 |
|
1214 | 1441 | } // End process loop |
| 1442 | + PROCESS_SWITCH(V0ptHadPiKaProt, processData, "Process Real Data", true); |
1215 | 1443 | }; |
1216 | 1444 |
|
1217 | 1445 | WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) |
|
0 commit comments