4646using namespace o2 ;
4747using namespace o2 ::framework;
4848using namespace o2 ::framework::expressions;
49- using selectedClusters = o2::soa::Filtered<o2::aod::JClusters>;
49+ using selectedClusters = o2::soa::Filtered<o2::soa::Join<o2:: aod::JClusters,o2::aod::JClusterTracks> >;
5050
5151#include " Framework/runDataProcessing.h"
5252
@@ -72,12 +72,11 @@ struct GammaJetTreeProducer {
7272 trackSelections{" trackSelections" , " globalTracks" , " set track selections" };
7373 Configurable<float > trackMinPt{" trackMinPt" , 0.15 , " minimum track pT cut" };
7474 Configurable<float > jetPtMin{" jetPtMin" , 5.0 , " minimum jet pT cut" };
75- Configurable<float > jetR{" jetR" , 0.4 , " jet resolution parameter" };
7675 Configurable<float > isoR{" isoR" , 0.4 , " isolation cone radius" };
77-
76+ Configurable<float > perpConeJetR{" perpConeJetR" , 0.4 , " perpendicular cone radius used to calculate perp cone rho for jet" };
77+ Configurable<float > trackMatchingEoverP{" trackMatchingEoverP" , 2.0 , " closest track is required to have E/p < value" };
7878 // cluster cuts
7979 Configurable<int > mClusterDefinition {" clusterDefinition" , 10 , " cluster definition to be selected, e.g. 10=kV3Default" };
80- // Preslice<o2::aod::JClusterTracks> perClusterMatchedTracks = o2::aod::jcluster::clusterId;
8180
8281 int mRunNumber = 0 ;
8382 int eventSelection = -1 ;
@@ -98,14 +97,20 @@ struct GammaJetTreeProducer {
9897 // create histograms
9998 LOG (info) << " Creating histograms" ;
10099
101- const o2Axis ptAxis{100 , 0 , 100 , " p_{T} (GeV/c)" };
100+ const o2Axis ptAxis{100 , 0 , 200 , " p_{T} (GeV/c)" };
102101 const o2Axis energyAxis{100 , 0 , 100 , " E (GeV)" };
103102 const o2Axis m02Axis{100 , 0 , 3 , " m02" };
103+ const o2Axis etaAxis{100 , -1 , 1 , " #eta" };
104+ const o2Axis phiAxis{100 , 0 , 2 *TMath::Pi (), " #phi" };
104105
105106 mHistograms .add (" clusterE" , " Energy of cluster" , o2HistType::kTH1F , {energyAxis});
106107 mHistograms .add (" trackPt" , " pT of track" , o2HistType::kTH1F , {ptAxis});
107108 mHistograms .add (" chjetPt" , " pT of charged jet" , o2HistType::kTH1F , {ptAxis});
109+ mHistograms .add (" chjetPtEtaPhi" , " pT of charged jet" , o2HistType::kTHnSparseF , {ptAxis, etaAxis, phiAxis});
108110 mHistograms .add (" chjetpt_vs_constpt" , " pT of charged jet vs pT of constituents" , o2HistType::kTH2F , {ptAxis, ptAxis});
111+
112+ // track QA THnSparse
113+ mHistograms .add (" trackPtEtaPhi" , " Track QA" , o2HistType::kTHnSparseF , {ptAxis, etaAxis, phiAxis});
109114 }
110115
111116 // ---------------------
@@ -161,7 +166,7 @@ struct GammaJetTreeProducer {
161166 double ptSumLeft = 0 ;
162167 double ptSumRight = 0 ;
163168
164- double cPhi = TVector2::Phi_0_2pi (cluster .phi ());
169+ double cPhi = TVector2::Phi_0_2pi (object .phi ());
165170
166171 // rotate cone left by 90 degrees
167172 float cPhiLeft = cPhi - TMath::Pi () / 2 ;
@@ -173,8 +178,8 @@ struct GammaJetTreeProducer {
173178 if (!isTrackSelected (track)) {
174179 continue ;
175180 }
176- dRLeft = jetutilities::deltaR (cluster .eta (), cPhiLeft, track.eta (), track.phi ());
177- dRRight = jetutilities::deltaR (cluster .eta (), cPhiRight, track.eta (), track.phi ());
181+ dRLeft = jetutilities::deltaR (object .eta (), cPhiLeft, track.eta (), track.phi ());
182+ dRRight = jetutilities::deltaR (object .eta (), cPhiRight, track.eta (), track.phi ());
178183
179184 if (dRLeft < radius) {
180185 ptSumLeft += track.pt ();
@@ -208,9 +213,12 @@ struct GammaJetTreeProducer {
208213 return ;
209214 }
210215
211- eventsTable (collision.multiplicity (), collision.centrality (), collision.rho (), collision.eventSel (), collision.alias_raw ());
216+ eventsTable (collision.multiplicity (), collision.centrality (), collision.rho (), collision.eventSel (), collision.trackOccupancyInTimeRange (), collision. alias_raw ());
212217 collisionMapping[collision.globalIndex ()] = eventsTable.lastIndex ();
213218
219+ // loop over tracks one time for QA
220+ runTrackQA (tracks);
221+
214222 // loop over clusters
215223 for (auto cluster : clusters) {
216224
@@ -226,26 +234,24 @@ struct GammaJetTreeProducer {
226234 // double dRMin = 100;
227235 double p = -1 ;
228236
229- // auto tracksofcluster = matchedtracks.sliceBy(perClusterMatchedTracks, cluster.globalIndex());
230- // for (const auto& match : tracksofcluster) {
231- // // ask the jtracks table for track with ID trackID
232- // double dR = deltaR(cluster.eta(), cluster.phi(), match.tracks_as<o2::aod::JTracks>().Eta(), match.tracks_as<o2::aod::JTracks>().Phi());
233- // if (dR < dRMin) {
234- // dRMin = dR;
235- // dEta = cluster.eta() - match.tracks_as<o2::aod::JTracks>().eta();
236- // dPhi = TVector2::Phi_0_2pi(cluster.phi()) - TVector2::Phi_0_2pi(match.tracks_as<o2::aod::JTracks>().phi());
237- // if (abs(dPhi) > M_PI) {
238- // dPhi = 2 * M_PI - abs(dPhi);
239- // }
240- // p = match.tracks_as<o2::aod::JTracks>().p();
241- // }
242- // }
243-
244- // // for compression reasons make dPhi and dEta 0 if no match is found
245- // if (p == -1) {
246- // dPhi = 0;
247- // dEta = 0;
248- // }
237+ // do track matching
238+ auto tracksofcluster = cluster.matchedTracks_as <JetTracks>();
239+
240+ for (auto track : tracksofcluster) {
241+ if (!isTrackSelected (track)) {
242+ continue ;
243+ }
244+ // find closest track that still has E/p < trackMatchingEoverP
245+ if (cluster.energy ()/track.p () > trackMatchingEoverP) {
246+ continue ;
247+ } else {
248+ // TODO make it eta on emcal surface and phi on emcal surface
249+ dEta = cluster.eta () - track.eta ();
250+ dPhi = RecoDecay::constrainAngle (RecoDecay::constrainAngle (track.phi (), -M_PI) - RecoDecay::constrainAngle (cluster.phi (), -M_PI), -M_PI);
251+ p = track.p ();
252+ break ;
253+ }
254+ }
249255
250256 gammasTable (eventsTable.lastIndex (), cluster.energy (), cluster.eta (), cluster.phi (), cluster.m02 (), cluster.m20 (), cluster.nCells (), cluster.time (), cluster.isExotic (), cluster.distanceToBadChannel (), cluster.nlm (), isoraw, perpconerho, dPhi, dEta, p);
251257 }
@@ -257,30 +263,38 @@ struct GammaJetTreeProducer {
257263 }
258264 PROCESS_SWITCH (GammaJetTreeProducer, processClusters, " Process EMCal clusters" , true );
259265
260- Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100 . 0f ) ;
266+ Filter jetCuts = aod::jet::pt > jetPtMin;
261267 // Process charged jets
262268 void processChargedJets (soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::JCollisionBCs>::iterator const & collision, soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>> const & chargedJets, aod::JetTracks const &)
263269 {
264270 // event selection
265271 if (!isEventAccepted (collision)) {
266272 return ;
267273 }
268-
274+ float leadingTrackPt = 0 ;
275+ ushort nconst = 0 ;
269276 // loop over charged jets
270277 for (auto jet : chargedJets) {
271278 if (jet.pt () < jetPtMin)
272279 continue ;
273- ushort nconst = 0 ;
280+ nconst = 0 ;
281+ leadingTrackPt = 0 ;
274282 // loop over constituents
275283 for (auto & constituent : jet.template tracks_as <aod::JetTracks>()) {
276284 mHistograms .fill (HIST (" chjetpt_vs_constpt" ), jet.pt (), constituent.pt ());
277285 nconst++;
286+ if (constituent.pt () > leadingTrackPt) {
287+ leadingTrackPt = constituent.pt ();
288+ }
278289 }
279290 int32_t storedColIndex = -1 ;
280291 if (auto foundCol = collisionMapping.find (collision.globalIndex ()); foundCol != collisionMapping.end ()) {
281292 storedColIndex = foundCol->second ;
282293 }
283- chargedJetsTable (storedColIndex, jet.pt (), jet.eta (), jet.phi (), jet.energy (), jet.mass (), jet.area (), nconst);
294+ // calculate perp cone rho
295+ double perpconerho = ch_perp_cone_rho (jet, tracks, perpConeJetR);
296+ mHistograms .fill (HIST (" chjetPtEtaPhi" ), jet.pt (), jet.eta (), jet.phi ());
297+ chargedJetsTable (storedColIndex, jet.pt (), jet.eta (), jet.phi (), jet.r (), jet.energy (), jet.mass (), jet.area (), leadingTrackPt, perpconerho,nconst);
284298 // fill histograms
285299 mHistograms .fill (HIST (" chjetPt" ), jet.pt ());
286300 }
0 commit comments