diff --git a/Analysis/DataModel/include/AnalysisDataModel/StrangenessTables.h b/Analysis/DataModel/include/AnalysisDataModel/StrangenessTables.h index d71b1c3e33256..40cc1d5d411c9 100644 --- a/Analysis/DataModel/include/AnalysisDataModel/StrangenessTables.h +++ b/Analysis/DataModel/include/AnalysisDataModel/StrangenessTables.h @@ -72,6 +72,18 @@ DECLARE_SOA_DYNAMIC_COLUMN(QtArm, qtarm, //! return std::sqrt(RecoDecay::P2(pxneg, pyneg, pzneg) - dp * dp / momTot); //qtarm }); +// Psi pair angle: angle between the plane defined by the electron and positron momenta and the xy plane +DECLARE_SOA_DYNAMIC_COLUMN(PsiPair, psipair, //! + [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) { + auto clipToPM1 = [](float x) { return x < -1.f ? -1.f : (x > 1.f ? 1.f : x); }; + float ptot2 = RecoDecay::P2(pxpos, pypos, pzpos) * RecoDecay::P2(pxneg, pyneg, pzneg); + float argcos = RecoDecay::dotProd(array{pxpos, pypos, pzpos}, array{pxneg, pyneg, pzneg}) / std::sqrt(ptot2); + float thetaPos = std::atan2(RecoDecay::sqrtSumOfSquares(pxpos, pypos), pzpos); + float thetaNeg = std::atan2(RecoDecay::sqrtSumOfSquares(pxneg, pyneg), pzneg); + float argsin = (thetaNeg - thetaPos) / std::acos(clipToPM1(argcos)); + return std::asin(clipToPM1(argsin)); + }); + //Calculated on the fly with mass assumption + dynamic tables DECLARE_SOA_DYNAMIC_COLUMN(MLambda, mLambda, //! [](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) -> float { return RecoDecay::M(array{array{pxpos, pypos, pzpos}, array{pxneg, pyneg, pzneg}}, array{RecoDecay::getMassPDG(kProton), RecoDecay::getMassPDG(kPiPlus)}); }); @@ -114,6 +126,7 @@ DECLARE_SOA_TABLE_FULL(StoredV0Datas, "V0Datas", "AOD", "V0DATA", //! v0data::DCAV0ToPV, v0data::Alpha, v0data::QtArm, + v0data::PsiPair, //Invariant masses v0data::MLambda, diff --git a/Analysis/Tasks/CMakeLists.txt b/Analysis/Tasks/CMakeLists.txt index 5d4dab653051c..b24d18ad592ef 100644 --- a/Analysis/Tasks/CMakeLists.txt +++ b/Analysis/Tasks/CMakeLists.txt @@ -11,6 +11,7 @@ add_subdirectory(PWGCF) add_subdirectory(PWGDQ) add_subdirectory(PWGHF) +add_subdirectory(PWGGA) add_subdirectory(PWGJE) add_subdirectory(PWGLF) add_subdirectory(PWGUD) diff --git a/Analysis/Tasks/PWGGA/CMakeLists.txt b/Analysis/Tasks/PWGGA/CMakeLists.txt new file mode 100644 index 0000000000000..9f42842030240 --- /dev/null +++ b/Analysis/Tasks/PWGGA/CMakeLists.txt @@ -0,0 +1,4 @@ +o2_add_dpl_workflow(gammaconversionsmc + SOURCES gammaConversionsMC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) diff --git a/Analysis/Tasks/PWGGA/gammaConversionsMC.cxx b/Analysis/Tasks/PWGGA/gammaConversionsMC.cxx new file mode 100644 index 0000000000000..42ecc1316599a --- /dev/null +++ b/Analysis/Tasks/PWGGA/gammaConversionsMC.cxx @@ -0,0 +1,335 @@ + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" + +#include "AnalysisDataModel/StrangenessTables.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" // for BigTracks + +#include "AnalysisDataModel/PID/PIDResponse.h" +#include "AnalysisDataModel/PID/PIDTPC.h" + +#include +#include + +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +using tracksAndTPCInfo = soa::Join; + +struct GammaConversionsMc { + + Configurable fSinglePtCut{"fSinglePtCut", 0.04, "minimum daughter track pt"}; + + Configurable fEtaCut{"fEtaCut", 0.8, "accepted eta range"}; + + Configurable fMinR{"fMinR", 5., "minimum conversion radius of the V0s"}; + Configurable fMaxR{"fMaxR", 180., "maximum conversion radius of the V0s"}; + + Configurable fPIDnSigmaBelowElectronLine{"fPIDnSigmaBelowElectronLine", -3., "minimum sigma electron PID for V0 daughter tracks"}; + + Configurable fPIDnSigmaAboveElectronLine{"fPIDnSigmaAboveElectronLine", 3., "maximum sigma electron PID for V0 daughter tracks"}; + + Configurable fPIDnSigmaAbovePionLine{"fPIDnSigmaAbovePionLine", 3., "minimum sigma to be over the pion line for low momentum tracks"}; //case 4: 3.0sigma, 1.0 sigma at high momentum + + Configurable fPIDnSigmaAbovePionLineHighP{"fPIDnSigmaAbovePionLineHighP", 1., "minimum sigma to be over the pion line for high momentum tracks"}; + + Configurable fPIDMinPnSigmaAbovePionLine{"fPIDMinPnSigmaAbovePionLine", 0.4, "minimum track momentum to apply any pion rejection"}; //case 7: // 0.4 GeV + + Configurable fPIDMaxPnSigmaAbovePionLine{"fPIDMaxPnSigmaAbovePionLine", 8., "border between low and high momentum pion rejection"}; //case 7: // 8. GeV + + Configurable fMinTPCFoundOverFindableCls{"fMinTPCNClsFoundOverFindable", 0.6, "minimum ratio found tpc clusters over findable"}; //case 9: // 0.6 + + Configurable fMinTPCCrossedRowsOverFindableCls{"fMinTPCCrossedRowsOverFindableCls", 0.0, "minimum ratio TPC crossed rows over findable clusters"}; + + Configurable fQtPtMax{"fQtPtMax", 0.11, "up to fQtMax, multiply the pt of the V0s by this value to get the maximum qt "}; + + Configurable fQtMax{"fQtMax", 0.040, "maximum qt"}; + Configurable fMaxPhotonAsymmetry{"fMaxPhotonAsymmetry", 0.95, "maximum photon asymetry"}; + Configurable fPsiPairCut{"fPsiPairCut", 0.1, "maximum psi angle of the track pair"}; + + Configurable fCosPAngleCut{"fCosPAngleCut", 0.85, "mimimum cosinus of the pointing angle"}; // case 4 + + HistogramRegistry registry{ + "registry", + { + {"IsPhotonSelected", "IsPhotonSelected", {HistType::kTH1F, {{13, -0.5f, 11.5f}}}}, + + {"beforeCuts/hPtRec", "hPtRec_before", {HistType::kTH1F, {{100, 0.0f, 25.0f}}}}, + {"beforeCuts/hEtaRec", "hEtaRec_before", {HistType::kTH1F, {{1000, -2.f, 2.f}}}}, + {"beforeCuts/hPhiRec", "hEtaRec_before", {HistType::kTH1F, {{1000, 0.f, 2.f * M_PI}}}}, + {"beforeCuts/hConvPointR", "hConvPointR_before", {HistType::kTH1F, {{800, 0.f, 200.f}}}}, + + {"hPtRec", "hPtRec", {HistType::kTH1F, {{100, 0.0f, 25.0f}}}}, + {"hEtaRec", "hEtaRec", {HistType::kTH1F, {{1000, -2.f, 2.f}}}}, + {"hPhiRec", "hEtaRec", {HistType::kTH1F, {{1000, 0.f, 2.f * M_PI}}}}, + {"hConvPointR", "hConvPointR", {HistType::kTH1F, {{800, 0.f, 200.f}}}}, + + {"hTPCdEdxSigEl", "hTPCdEdxSigEl", {HistType::kTH2F, {{150, 0.03f, 20.f}, {400, -10.f, 10.f}}}}, + {"hTPCdEdxSigPi", "hTPCdEdxSigPi", {HistType::kTH2F, {{150, 0.03f, 20.f}, {400, -10.f, 10.f}}}}, + {"hTPCdEdx", "hTPCdEdx", {HistType::kTH2F, {{150, 0.03f, 20.f}, {800, 0.f, 200.f}}}}, + + {"hTPCFoundOverFindableCls", "hTPCFoundOverFindableCls", {HistType::kTH1F, {{100, 0.f, 1.f}}}}, + {"hTPCCrossedRowsOverFindableCls", "hTPCCrossedRowsOverFindableCls", {HistType::kTH1F, {{100, 0.f, 1.5f}}}}, + + {"hArmenteros", "hArmenteros", {HistType::kTH2F, {{200, -1.f, 1.f}, {250, 0.f, 25.f}}}}, + {"hPsiPtRec", "hPsiPtRec", {HistType::kTH2F, {{500, -2.f, 2.f}, {100, 0.f, 25.f}}}}, + + {"hCosPAngle", "hCosPAngle", {HistType::kTH1F, {{1000, -1.f, 1.f}}}}, + + // resolution histos + {"resolutions/hPtRes", "hPtRes_Rec-MC", {HistType::kTH1F, {{100, -0.5f, 0.5f}}}}, + {"resolutions/hEtaRes", "hEtaRes_Rec-MC", {HistType::kTH1F, {{100, -2.f, 2.f}}}}, + {"resolutions/hPhiRes", "hPhiRes_Rec-MC", {HistType::kTH1F, {{100, -M_PI, M_PI}}}}, + + {"resolutions/hConvPointRRes", "hConvPointRRes_Rec-MC", {HistType::kTH1F, {{100, -200.f, 200.f}}}}, + {"resolutions/hConvPointAbsoluteDistanceRes", "hConvPointAbsoluteDistanceRes", {HistType::kTH1F, {{100, -0.0f, 200.f}}}}, + }, + }; + + enum photonCuts { + kPhotonIn = 0, + kTrackEta, + kTrackPt, + kElectronPID, + kPionRejLowMom, + kPionRejHighMom, + kTPCFoundOverFindableCls, + kTPCCrossedRowsOverFindableCls, + kV0Radius, + kArmenteros, + kPsiPair, + kCosinePA, + kPhotonOut + }; + + std::vector fPhotCutsLabels{ + "kPhotonIn", + "kTrackEta", + "kTrackPt", + "kElectronPID", + "kPionRejLowMom", + "kPionRejHighMom", + "kTPCFoundOverFindableCls", + "kTPCCrossedRowsOverFindableCls", + "kV0Radius", + "kArmenteros", + "kPsiPair", + "kCosinePA", + "kPhotonOut"}; + + void init(InitContext const&) + { + TAxis* lXaxis = registry.get(HIST("IsPhotonSelected"))->GetXaxis(); + for (size_t i = 0; i < fPhotCutsLabels.size(); ++i) { + lXaxis->SetBinLabel(i + 1, fPhotCutsLabels[i]); + } + } + + void process(aod::Collision const& theCollision, + aod::V0Datas const& theV0s, + tracksAndTPCInfo const& theTracks, + aod::McParticles const& theMcParticles) + { + for (auto& lV0 : theV0s) { + + fillHistogramsBeforeCuts(lV0); + + auto lTrackPos = lV0.template posTrack_as(); //positive daughter + auto lTrackNeg = lV0.template negTrack_as(); //negative daughter + + // apply track cuts + if (!(trackPassesCuts(lTrackPos) && trackPassesCuts(lTrackNeg))) { + continue; + } + + float lV0CosinePA = RecoDecay::CPA(array{theCollision.posX(), theCollision.posY(), theCollision.posZ()}, array{lV0.x(), lV0.y(), lV0.z()}, array{lV0.px(), lV0.py(), lV0.pz()}); + + // apply photon cuts + if (!passesPhotonCuts(lV0, lV0CosinePA)) { + continue; + } + + fillHistogramsAfterCuts(lV0, lTrackPos, lTrackNeg, lV0CosinePA); + + processTruePhotons(lV0, lTrackPos, lTrackNeg, theMcParticles); + } + } + + template + void fillHistogramsBeforeCuts(const T& theV0) + { + // fill some QA histograms before any cuts + registry.fill(HIST("beforeCuts/hPtRec"), theV0.pt()); + registry.fill(HIST("beforeCuts/hEtaRec"), theV0.eta()); + registry.fill(HIST("beforeCuts/hPhiRec"), theV0.phi()); + registry.fill(HIST("beforeCuts/hConvPointR"), theV0.v0radius()); + registry.fill(HIST("IsPhotonSelected"), kPhotonIn); + } + + template + bool trackPassesCuts(const T& theTrack) + { + + // single track eta cut + if (TMath::Abs(theTrack.eta()) > fEtaCut) { + registry.fill(HIST("IsPhotonSelected"), kTrackEta); + return kFALSE; + } + + // single track pt cut + if (theTrack.pt() < fSinglePtCut) { + registry.fill(HIST("IsPhotonSelected"), kTrackPt); + return kFALSE; + } + + if (!(selectionPIDTPC_track(theTrack))) { + return kFALSE; + } + + if (theTrack.tpcFoundOverFindableCls() < fMinTPCFoundOverFindableCls) { + registry.fill(HIST("IsPhotonSelected"), kTPCFoundOverFindableCls); + return kFALSE; + } + + if (theTrack.tpcCrossedRowsOverFindableCls() < fMinTPCCrossedRowsOverFindableCls) { + registry.fill(HIST("IsPhotonSelected"), kTPCCrossedRowsOverFindableCls); + return kFALSE; + } + + return kTRUE; + } + + template + bool passesPhotonCuts(const T& theV0, float theV0CosinePA) + { + if (theV0.v0radius() < fMinR || theV0.v0radius() > fMaxR) { + registry.fill(HIST("IsPhotonSelected"), kV0Radius); + return kFALSE; + } + + if (!ArmenterosQtCut(theV0.alpha(), theV0.qtarm(), theV0.pt())) { + registry.fill(HIST("IsPhotonSelected"), kArmenteros); + return kFALSE; + } + + if (TMath::Abs(theV0.psipair()) > fPsiPairCut) { + registry.fill(HIST("IsPhotonSelected"), kPsiPair); + return kFALSE; + } + + if (theV0CosinePA < fCosPAngleCut) { + registry.fill(HIST("IsPhotonSelected"), kCosinePA); + return kFALSE; + } + + return kTRUE; + } + + template + void fillHistogramsAfterCuts(const TV0& theV0, const TTRACK& theTrackPos, const TTRACK& theTrackNeg, float theV0CosinePA) + { + registry.fill(HIST("IsPhotonSelected"), kPhotonOut); + + registry.fill(HIST("hPtRec"), theV0.pt()); + registry.fill(HIST("hEtaRec"), theV0.eta()); + registry.fill(HIST("hPhiRec"), theV0.phi()); + registry.fill(HIST("hConvPointR"), theV0.v0radius()); + + registry.fill(HIST("hTPCdEdxSigEl"), theTrackNeg.p(), theTrackNeg.tpcNSigmaEl()); + registry.fill(HIST("hTPCdEdxSigEl"), theTrackPos.p(), theTrackPos.tpcNSigmaEl()); + registry.fill(HIST("hTPCdEdxSigPi"), theTrackNeg.p(), theTrackNeg.tpcNSigmaPi()); + registry.fill(HIST("hTPCdEdxSigPi"), theTrackPos.p(), theTrackPos.tpcNSigmaPi()); + + registry.fill(HIST("hTPCdEdx"), theTrackNeg.p(), theTrackNeg.tpcSignal()); + registry.fill(HIST("hTPCdEdx"), theTrackPos.p(), theTrackPos.tpcSignal()); + + registry.fill(HIST("hTPCFoundOverFindableCls"), theTrackNeg.tpcFoundOverFindableCls()); + registry.fill(HIST("hTPCFoundOverFindableCls"), theTrackPos.tpcFoundOverFindableCls()); + + registry.fill(HIST("hTPCCrossedRowsOverFindableCls"), theTrackNeg.tpcCrossedRowsOverFindableCls()); + registry.fill(HIST("hTPCCrossedRowsOverFindableCls"), theTrackPos.tpcCrossedRowsOverFindableCls()); + + registry.fill(HIST("hArmenteros"), theV0.alpha(), theV0.qtarm()); + registry.fill(HIST("hPsiPtRec"), theV0.psipair(), theV0.pt()); + + registry.fill(HIST("hCosPAngle"), theV0CosinePA); + } + + template + void processTruePhotons(const TV0& theV0, const TTRACK& theTrackPos, const TTRACK& theTrackNeg, const TMC& theMcParticles) + { + // todo: verify it is enough to check only mother0 being equal + if (theTrackPos.mcParticle().mother0() > -1 && + theTrackPos.mcParticle().mother0() == theTrackNeg.mcParticle().mother0()) { + auto lMother = theMcParticles.iteratorAt(theTrackPos.mcParticle().mother0()); + + if (lMother.pdgCode() == 22) { + + registry.fill(HIST("resolutions/hPtRes"), theV0.pt() - lMother.pt()); + registry.fill(HIST("resolutions/hEtaRes"), theV0.eta() - lMother.eta()); + registry.fill(HIST("resolutions/hPhiRes"), theV0.phi() - lMother.phi()); + + TVector3 lConvPointRec(theV0.x(), theV0.y(), theV0.z()); + TVector3 lPosTrackVtxMC(theTrackPos.mcParticle().vx(), theTrackPos.mcParticle().vy(), theTrackPos.mcParticle().vz()); + // take the origin of the positive mc track as conversion point (should be identical with negative, verified this on a few photons) + TVector3 lConvPointMC(lPosTrackVtxMC); + + registry.fill(HIST("resolutions/hConvPointRRes"), lConvPointRec.Perp() - lConvPointMC.Perp()); + registry.fill(HIST("resolutions/hConvPointAbsoluteDistanceRes"), TVector3(lConvPointRec - lConvPointMC).Mag()); + } + } + } + + Bool_t ArmenterosQtCut(Double_t theAlpha, Double_t theQt, Double_t thePt) + { + // in AliPhysics this is the cut for if fDo2DQt && fDoQtGammaSelection == 2 + Float_t lQtMaxPtDep = fQtPtMax * thePt; + if (lQtMaxPtDep > fQtMax) { + lQtMaxPtDep = fQtMax; + } + if (!(TMath::Power(theAlpha / fMaxPhotonAsymmetry, 2) + TMath::Power(theQt / lQtMaxPtDep, 2) < 1)) { + return kFALSE; + } + return kTRUE; + } + + template + bool selectionPIDTPC_track(const T& theTrack) + { + // TPC Electron Line + if (theTrack.tpcNSigmaEl() < fPIDnSigmaBelowElectronLine || theTrack.tpcNSigmaEl() > fPIDnSigmaAboveElectronLine) { + registry.fill(HIST("IsPhotonSelected"), kElectronPID); + return kFALSE; + } + + // TPC Pion Line + if (theTrack.p() > fPIDMinPnSigmaAbovePionLine) { + // low pt Pion rej + if (theTrack.p() < fPIDMaxPnSigmaAbovePionLine) { + if (theTrack.tpcNSigmaEl() > fPIDnSigmaBelowElectronLine && theTrack.tpcNSigmaEl() < fPIDnSigmaAboveElectronLine && theTrack.tpcNSigmaPi() < fPIDnSigmaAbovePionLine) { + registry.fill(HIST("IsPhotonSelected"), kPionRejLowMom); + return kFALSE; + } + } + // High Pt Pion rej + else { + if (theTrack.tpcNSigmaEl() > fPIDnSigmaBelowElectronLine && theTrack.tpcNSigmaEl() < fPIDnSigmaAboveElectronLine && theTrack.tpcNSigmaPi() < fPIDnSigmaAbovePionLineHighP) { + registry.fill(HIST("IsPhotonSelected"), kPionRejHighMom); + return kFALSE; + } + } + } + return kTRUE; + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/Framework/Core/include/Framework/AnalysisDataModel.h b/Framework/Core/include/Framework/AnalysisDataModel.h index 04e4c12d0b9ed..dd37035c9b9af 100644 --- a/Framework/Core/include/Framework/AnalysisDataModel.h +++ b/Framework/Core/include/Framework/AnalysisDataModel.h @@ -239,7 +239,13 @@ DECLARE_SOA_DYNAMIC_COLUMN(ITSNClsInnerBarrel, itsNClsInnerBarrel, //! Number of return itsNclsInnerBarrel; }); -DECLARE_SOA_DYNAMIC_COLUMN(TPCCrossedRowsOverFindableCls, tpcCrossedRowsOverFindableCls, //! Ratio of crossed rows over findable clusters +DECLARE_SOA_DYNAMIC_COLUMN(TPCFoundOverFindableCls, tpcFoundOverFindableCls, //! Ratio of found over findable clusters + [](uint8_t tpcNClsFindable, int8_t tpcNClsFindableMinusFound) -> float { + int16_t tpcNClsFound = (int16_t)tpcNClsFindable - tpcNClsFindableMinusFound; + return (float)tpcNClsFound / (float)tpcNClsFindable; + }); + +DECLARE_SOA_DYNAMIC_COLUMN(TPCCrossedRowsOverFindableCls, tpcCrossedRowsOverFindableCls, //! Ratio crossed rows over findable clusters [](uint8_t tpcNClsFindable, int8_t tpcNClsFindableMinusCrossedRows) -> float { int16_t tpcNClsCrossedRows = (int16_t)tpcNClsFindable - tpcNClsFindableMinusCrossedRows; return (float)tpcNClsCrossedRows / (float)tpcNClsFindable; @@ -303,6 +309,7 @@ DECLARE_SOA_TABLE(TracksExtra, "AOD", "TRACKEXTRA", //! Additional track informa track::TPCNClsCrossedRows, track::ITSNCls, track::ITSNClsInnerBarrel, track::TPCCrossedRowsOverFindableCls, + track::TPCFoundOverFindableCls, track::TPCFractionSharedCls, track::TrackEtaEMCAL, track::TrackPhiEMCAL);