Skip to content

Commit 2c33ec0

Browse files
smaff92alibuild
smaff92
and
alibuild
authored
[PWGJE] Adding MC weights for JEJE injections, and adding mask for s… (#11643)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent d2d3acd commit 2c33ec0

File tree

1 file changed

+67
-77
lines changed

1 file changed

+67
-77
lines changed

‎PWGJE/Tasks/statPromptPhoton.cxx‎

Lines changed: 67 additions & 77 deletions
Original file line numberDiff line numberDiff line change
@@ -13,30 +13,43 @@
1313
/// \brief Reconstruction of Phi yield through track-track Minv correlations for resonance hadrochemistry analysis.
1414
///
1515
///
16-
/// \author Adrian Fereydon Nassirpour <adrian.fereydon.nassirpour@cern.ch>
16+
/// \author Adrian Fereydon Nassirpour <adrian.fereydon.nassirpour@cern.ch>
1717

18+
#include "PWGJE/Core/FastJetUtilities.h"
1819
#include "PWGJE/Core/JetDerivedDataUtilities.h"
1920
#include "PWGJE/DataModel/EMCALClusters.h"
20-
#include "PWGJE/DataModel/JetReducedData.h"
21+
#include "PWGJE/DataModel/Jet.h"
2122

23+
#include "Common/Core/RecoDecay.h"
24+
#include "Common/Core/TrackSelection.h"
25+
#include "Common/Core/TrackSelectionDefaults.h"
26+
#include "Common/Core/trackUtilities.h"
2227
#include "Common/DataModel/EventSelection.h"
28+
#include "Common/DataModel/Multiplicity.h"
29+
#include "Common/DataModel/PIDResponse.h"
2330
#include "Common/DataModel/TrackSelectionTables.h"
2431

32+
#include "CommonConstants/PhysicsConstants.h"
33+
#include "CommonDataFormat/InteractionRecord.h"
34+
#include "DataFormatsEMCAL/AnalysisCluster.h"
35+
#include "DataFormatsEMCAL/Cell.h"
36+
#include "DataFormatsEMCAL/Constants.h"
37+
#include "DataFormatsParameters/GRPMagField.h"
38+
#include "DataFormatsParameters/GRPObject.h"
39+
#include "DetectorsBase/Propagator.h"
40+
#include "EMCALBase/Geometry.h"
41+
#include "EMCALCalib/BadChannelMap.h"
2542
#include "Framework/ASoA.h"
2643
#include "Framework/AnalysisDataModel.h"
2744
#include "Framework/AnalysisTask.h"
2845
#include "Framework/HistogramRegistry.h"
29-
#include <Framework/Configurable.h>
30-
#include <Framework/HistogramSpec.h>
31-
#include <Framework/InitContext.h>
32-
#include <Framework/OutputObjHeader.h>
33-
#include <Framework/runDataProcessing.h>
46+
#include "Framework/runDataProcessing.h"
47+
#include "ReconstructionDataFormats/Track.h"
48+
#include <CCDB/BasicCCDBManager.h>
3449

3550
#include <TLorentzVector.h>
36-
#include <TMath.h>
3751
#include <TVector2.h>
3852

39-
#include <cmath>
4053
#include <iostream>
4154
#include <optional>
4255
#include <string>
@@ -47,6 +60,7 @@ using namespace o2::framework;
4760
using namespace o2::framework::expressions;
4861

4962
struct statPromptPhoton {
63+
5064
SliceCache cache;
5165
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
5266
Configurable<double> cfgMaxDCArToPVcut{"cfgMaxDCArToPVcut", 0.5, "Track DCAr cut to PV Maximum"};
@@ -85,15 +99,20 @@ struct statPromptPhoton {
8599
Configurable<bool> cfgGenHistograms{"cfgGenHistograms", false, "Enables Generated histograms"};
86100
Configurable<bool> cfgRecHistograms{"cfgRecHistograms", false, "Enables Reconstructed histograms"};
87101
Configurable<bool> cfgDataHistograms{"cfgDataHistograms", false, "Enables Data histograms"};
102+
Configurable<bool> cfgSkimmedTrigger{"cfgSkimmedTrigger", false, "Enables trigger for skimmied datasets (2023 onwards)"};
103+
Configurable<std::string> cfgTriggerMasks{"cfgTriggerMasks", "", "possible JE Trigger masks: fJetChLowPt,fJetChHighPt,fTrackLowPt,fTrackHighPt,fJetD0ChLowPt,fJetD0ChHighPt,fJetLcChLowPt,fJetLcChHighPt,fEMCALReadout,fJetFullHighPt,fJetFullLowPt,fJetNeutralHighPt,fJetNeutralLowPt,fGammaVeryHighPtEMCAL,fGammaVeryHighPtDCAL,fGammaHighPtEMCAL,fGammaHighPtDCAL,fGammaLowPtEMCAL,fGammaLowPtDCAL,fGammaVeryLowPtEMCAL,fGammaVeryLowPtDCAL"};
88104
Configurable<bool> cfgDebug{"cfgDebug", false, "Enables debug information for local running"};
105+
89106
int trackFilter = -1;
107+
std::vector<int> triggerMaskBits;
90108

91109
// INIT
92110
void init(InitContext const&)
93111
{
94112
std::vector<double> ptBinning = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 12.0, 16.0, 20.0, 25.0, 30.0, 40.0, 50.0, 75.0, 100.0, 150.0, 200.0, 300.0, 500.0};
95113
AxisSpec pthadAxis = {ptBinning, "#it{p}_{T}^{had sum} [GeV/c]"};
96114

115+
triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(cfgTriggerMasks);
97116
if (cfgJETracks) {
98117
trackFilter = jetderiveddatautilities::initialiseTrackSelection(static_cast<std::string>(cfgTrackFilter));
99118
}
@@ -235,8 +254,11 @@ struct statPromptPhoton {
235254

236255
using jMCClusters = o2::soa::Join<o2::aod::JMcClusterLbs, o2::aod::JClusters, o2::aod::JClusterTracks>;
237256
using jClusters = o2::soa::Join<o2::aod::JClusters, o2::aod::JClusterTracks>;
238-
using jselectedCollisions = soa::Join<aod::JCollisions, aod::JCollisionBCs, aod::JCollisionPIs, aod::EvSels, aod::JEMCCollisionLbs>;
257+
using jselectedCollisions = soa::Join<aod::JCollisions, aod::JCollisionBCs, aod::JCollisionPIs, aod::EvSels, aod::JEMCCollisionLbs, aod::JMcCollisionLbs>;
258+
using jselectedDataCollisions = soa::Join<aod::JCollisions, aod::JCollisionBCs, aod::JCollisionPIs, aod::EvSels, aod::JEMCCollisionLbs>;
259+
// using jselectedDataCollisions = soa::Join<aod::JCollisions, aod::JCollisionBCs, aod::JCollisionPIs, aod::JCollisionMcInfos, aod::EvSels, aod::JEMCCollisionLbs>;
239260
using jfilteredCollisions = soa::Filtered<jselectedCollisions>;
261+
using jfilteredDataCollisions = soa::Filtered<jselectedDataCollisions>;
240262
using jfilteredMCClusters = soa::Filtered<jMCClusters>;
241263
using jfilteredClusters = soa::Filtered<jClusters>;
242264

@@ -588,7 +610,7 @@ struct statPromptPhoton {
588610

589611
PresliceUnsorted<jEMCtracks> EMCTrackPerTrack = aod::jemctrack::trackId;
590612
int nEventsRecMC_JE = 0;
591-
void processMCRec_JE(jfilteredCollisions::iterator const& collision, jfilteredMCClusters const& mcclusters, jTrackCandidates const& tracks, soa::Join<aod::JTracks, aod::JTrackExtras, aod::JTrackPIs> const&, TrackCandidates const&, aod::JMcParticles const&, BcCandidates const&, jEMCtracks const& emctracks)
613+
void processMCRec_JE(jfilteredCollisions::iterator const& collision, jfilteredMCClusters const& mcclusters, jTrackCandidates const& tracks, soa::Join<aod::JTracks, aod::JTrackExtras, aod::JTrackPIs> const&, TrackCandidates const&, aod::JMcParticles const&, BcCandidates const&, jEMCtracks const& emctracks, aod::JetMcCollisions const&)
592614
{
593615

594616
nEventsRecMC_JE++;
@@ -613,6 +635,19 @@ struct statPromptPhoton {
613635
}
614636
histos.fill(HIST("REC_nEvents"), 2.5);
615637

638+
if (cfgSkimmedTrigger) {
639+
if (!jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) {
640+
return;
641+
}
642+
} // JE Software Triggers
643+
644+
histos.fill(HIST("REC_nEvents"), 3.5);
645+
646+
double weight = 1;
647+
if (collision.has_mcCollision()) {
648+
weight = collision.mcCollision().weight();
649+
}
650+
616651
bool noTrk = true;
617652
for (auto& track : tracks) {
618653
if (cfgJETracks) {
@@ -828,9 +863,9 @@ struct statPromptPhoton {
828863
histos.fill(HIST("REC_Cluster_QA"), 4.5);
829864
clustertrigger = true;
830865
double pthadsum = GetPtHadSum(tracks, mccluster, cfgMinR, cfgMaxR, false, false, true);
831-
histos.fill(HIST("REC_Trigger_V_PtHadSum_Photon"), mccluster.energy(), pthadsum);
832-
histos.fill(HIST("REC_PtHadSum_Photon"), pthadsum);
833-
histos.fill(HIST("REC_Trigger_Energy"), mccluster.energy());
866+
histos.fill(HIST("REC_Trigger_V_PtHadSum_Photon"), mccluster.energy(), pthadsum, weight);
867+
histos.fill(HIST("REC_PtHadSum_Photon"), pthadsum, weight);
868+
histos.fill(HIST("REC_Trigger_Energy"), mccluster.energy(), weight);
834869
}
835870

836871
auto ClusterParticles = mccluster.mcParticles_as<aod::JMcParticles>();
@@ -839,11 +874,6 @@ struct statPromptPhoton {
839874
bool goodgentrigger = true;
840875
double chPe = 0;
841876
for (auto& clusterparticle : ClusterParticles) {
842-
// double etaP = clusterparticle.eta();
843-
// double etaC = mccluster.eta();
844-
// double phiP = clusterparticle.phi();
845-
// double phiC = mccluster.phi();
846-
// double ptP = clusterparticle.pt();
847877
int cindex = clusterparticle.globalIndex();
848878
double pdgcode = fabs(clusterparticle.pdgCode());
849879
if (!clusterparticle.isPhysicalPrimary()) {
@@ -884,8 +914,6 @@ struct statPromptPhoton {
884914
histos.fill(HIST("REC_Cluster_ParticleWITHtrack_Phi"), clusterparticle.phi());
885915
histos.fill(HIST("REC_Cluster_ParticleWITHtrack_Eta"), clusterparticle.eta());
886916
histos.fill(HIST("REC_Cluster_ParticleWITHtrack_Pt_Phi"), clusterparticle.pt(), clusterparticle.phi());
887-
// if (phiPrimeP > (0.12/ptP + TMath::Pi()/18. + 0.035) ||
888-
// phiPrimeP < (0.1/ptP/ptP + TMath::Pi()/18. - 0.025) ) {
889917
histos.fill(HIST("REC_Cluster_ParticleWITHtrack_Pt_PhiPrime"), ptP, phiPrimeP);
890918
if (photontrigger) {
891919
histos.fill(HIST("REC_Impurity_ParticleWITHtrack_Pt_PhiPrime"), ptP, phiPrimeP);
@@ -904,8 +932,6 @@ struct statPromptPhoton {
904932
histos.fill(HIST("REC_Cluster_ParticleWITHOUTtrack_Phi"), clusterparticle.phi());
905933
histos.fill(HIST("REC_Cluster_ParticleWITHOUTtrack_Eta"), clusterparticle.eta());
906934
histos.fill(HIST("REC_Cluster_ParticleWITHOUTtrack_Pt_Phi"), clusterparticle.pt(), clusterparticle.phi());
907-
// if (phiPrimeP > (0.12/ptP + TMath::Pi()/18. + 0.035) ||
908-
// phiPrimeP < (0.1/ptP/ptP + TMath::Pi()/18. - 0.025) ) {
909935
histos.fill(HIST("REC_Cluster_ParticleWITHOUTtrack_Pt_PhiPrime"), ptP, phiPrimeP);
910936
if (photontrigger) {
911937
histos.fill(HIST("REC_Impurity_ParticleWITHOUTtrack_Pt_PhiPrime"), ptP, phiPrimeP);
@@ -964,12 +990,12 @@ struct statPromptPhoton {
964990
std::cout << "Photon mom 2: " << mom2 << std::endl;
965991
}
966992
if (std::abs(clusterparticle.getGenStatusCode()) > 19 && std::abs(clusterparticle.getGenStatusCode()) < 90) {
967-
histos.fill(HIST("REC_True_Prompt_Trigger_Energy"), clusterparticle.e());
993+
histos.fill(HIST("REC_True_Prompt_Trigger_Energy"), clusterparticle.e(), weight);
968994
TLorentzVector lRealPhoton;
969995
lRealPhoton.SetPxPyPzE(clusterparticle.px(), clusterparticle.py(), clusterparticle.pz(), clusterparticle.e());
970996
double truepthadsum = GetPtHadSum(tracks, lRealPhoton, cfgMinR, cfgMaxR, false, false, false);
971997
truephotonPt = clusterparticle.e();
972-
histos.fill(HIST("REC_TrueTrigger_V_PtHadSum_Photon"), truephotonPt, truepthadsum);
998+
histos.fill(HIST("REC_TrueTrigger_V_PtHadSum_Photon"), truephotonPt, truepthadsum, weight);
973999
}
9741000
} // photon check
9751001
} // clusterparticle loop
@@ -1000,39 +1026,6 @@ struct statPromptPhoton {
10001026
}
10011027
} // cluster loop
10021028

1003-
// auto bc = collision.bc_as<BcCandidates>();
1004-
// int rnr = bc.runNumber();
1005-
1006-
// std::string rnrstring = std::to_string(rnr);
1007-
// if (runs.find(rnrstring) == std::string::npos) {
1008-
// std::cout<<"++++++++++++++++++++++++++++++++"<<std::endl;
1009-
// std::cout<<"FETCHING NEW RUN NUMBER FOR RUN: "<<rnr<<std::endl;
1010-
// std::cout<<"++++++++++++++++++++++++++++++++"<<std::endl;
1011-
1012-
// std::string ccdbpath="GLO/Config/GRPMagField";
1013-
// static o2::parameters::GRPMagField* grpmag = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(ccdbpath, bc.timestamp());
1014-
// if(grpmag) {
1015-
// bfield = std::lround(5.f * grpmag->getL3Current() / 30000.f);
1016-
// std::cout<<"++++++++++++++++++++++++++++++++"<<std::endl;
1017-
// std::cout<<"MAG FIELD IS: "<<bfield<<std::endl;
1018-
// std::cout<<"++++++++++++++++++++++++++++++++"<<std::endl;
1019-
// }
1020-
// else {
1021-
// ccdbpath="GLO/GRP/GRP";
1022-
// static o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp<o2::parameters::GRPObject>(ccdbpath, bc.timestamp());
1023-
// if(!grpo) {
1024-
// std::cout<<"WE CAN NEITHER FETCH GRPMAG OR GRPO!!! SHIT IS SCREWED"<<std::endl;
1025-
// }
1026-
// bfield = grpo->getNominalL3Field();
1027-
// }
1028-
// bfield = 5;
1029-
// runs += rnrstring;
1030-
// std::cout << "++++++++++++++++++++++++++++++++" << std::endl;
1031-
// std::cout << "Run is now appended to string: " << runs << std::endl;
1032-
// std::cout << "++++++++++++++++++++++++++++++++" << std::endl;
1033-
1034-
// } // check mag field for current run number: done!
1035-
10361029
// clusters done, now we do the sternheimer tracks
10371030
for (auto& track : tracks) {
10381031
bool sterntrigger = false;
@@ -1057,19 +1050,9 @@ struct statPromptPhoton {
10571050
phiPrime = 2 * TMath::Pi() - phiPrime;
10581051
}
10591052

1060-
// if (bfield < 0) {
1061-
// phiPrime = 2 * TMath::Pi() - phiPrime;
1062-
// }
1063-
10641053
phiPrime = phiPrime + TMath::Pi() / 18.;
10651054
phiPrime = fmod(phiPrime, 2 * TMath::Pi() / 18.);
1066-
// double pt = track.pt();
1067-
// if (phiPrime > (0.12/pt + TMath::Pi()/18. + 0.035) ||
1068-
// phiPrime < (0.1/pt/pt + TMath::Pi()/18. - 0.025) ) {
10691055
histos.fill(HIST("REC_Track_PhiPrime_Pt"), phiPrime, track.pt());
1070-
// }//geo cut
1071-
// Done with geometric cuts
1072-
10731056
histos.fill(HIST("REC_Track_Pt"), track.pt());
10741057
histos.fill(HIST("REC_Track_Phi"), track.phi());
10751058
if (clustertrigger) {
@@ -1083,12 +1066,12 @@ struct statPromptPhoton {
10831066
}
10841067
}
10851068
double pthadsum = GetPtHadSum(tracks, track, cfgMinR, cfgMaxR, true, false, true);
1086-
histos.fill(HIST("REC_Trigger_V_PtHadSum_Nch"), sternPt, pthadsum);
1069+
histos.fill(HIST("REC_Trigger_V_PtHadSum_Nch"), sternPt, pthadsum, weight);
10871070
if (sterntrigger) {
10881071
bool doStern = true;
10891072
double sterncount = 1.0;
10901073
while (doStern) {
1091-
histos.fill(HIST("REC_Trigger_V_PtHadSum_Stern"), sterncount, pthadsum, 2.0 / sternPt);
1074+
histos.fill(HIST("REC_Trigger_V_PtHadSum_Stern"), sterncount, pthadsum, (2.0 / sternPt) * weight);
10921075
if (sterncount < sternPt) {
10931076
sterncount++;
10941077
} else {
@@ -1103,18 +1086,21 @@ struct statPromptPhoton {
11031086
PROCESS_SWITCH(statPromptPhoton, processMCRec_JE, "processJE MC data", false);
11041087

11051088
int nEventsData = 0;
1106-
void processData(jfilteredCollisions::iterator const& collision, jfilteredClusters const& clusters, jDataTrackCandidates const& tracks, soa::Join<aod::JTracks, aod::JTrackExtras, aod::JTrackPIs> const&, TrackCandidates const&, BcCandidates const&, jEMCtracks const& emctracks)
1089+
void processData(jfilteredDataCollisions::iterator const& collision, jfilteredClusters const& clusters, jDataTrackCandidates const& tracks, soa::Join<aod::JTracks, aod::JTrackExtras, aod::JTrackPIs> const&, TrackCandidates const&, BcCandidates const&, jEMCtracks const& emctracks)
11071090
{
1108-
11091091
nEventsData++;
11101092
if (cfgDebug) {
11111093
if (nEventsData == 1) {
11121094
std::cout << "Starting Data Processing: " << nEventsData << std::endl;
11131095
}
11141096
if ((nEventsData + 1) % 10000 == 0) {
11151097
std::cout << "Processed Data Events: " << nEventsData << std::endl;
1098+
std::cout << "Events Trigger Bit: " << collision.triggerSel() << std::endl;
1099+
std::cout << "Trigger Mask Bit: " << triggerMaskBits[0] << std::endl;
1100+
std::cout << "Trigger Mask Cfg Line: " << cfgTriggerMasks << std::endl;
11161101
}
11171102
}
1103+
11181104
histos.fill(HIST("DATA_nEvents"), 0.5);
11191105

11201106
// required cuts
@@ -1124,13 +1110,21 @@ struct statPromptPhoton {
11241110
return;
11251111

11261112
histos.fill(HIST("DATA_nEvents"), 1.5);
1127-
11281113
if (cfgEmcTrigger) {
11291114
if (!collision.isEmcalReadout())
11301115
return;
11311116
}
1117+
11321118
histos.fill(HIST("DATA_nEvents"), 2.5);
11331119

1120+
if (cfgSkimmedTrigger) {
1121+
if (!jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) {
1122+
return;
1123+
}
1124+
} // JE Software Triggers
1125+
1126+
histos.fill(HIST("DATA_nEvents"), 3.5);
1127+
11341128
bool noTrk = true;
11351129
for (auto& track : tracks) {
11361130

@@ -1331,10 +1325,6 @@ struct statPromptPhoton {
13311325
phiPrime = 2 * TMath::Pi() - phiPrime;
13321326
}
13331327

1334-
// if (bfield < 0) {
1335-
// phiPrime = 2 * TMath::Pi() - phiPrime;
1336-
// }
1337-
13381328
phiPrime = phiPrime + TMath::Pi() / 18.;
13391329
phiPrime = fmod(phiPrime, 2 * TMath::Pi() / 18.);
13401330
double pt = track.pt();

0 commit comments

Comments
 (0)