Skip to content

Commit 2bf0277

Browse files
sashingosashingo
committed
added UE study for Z->ee and multiplicity analysis in pp
1 parent f106b3b commit 2bf0277

File tree

1 file changed

+105
-8
lines changed

1 file changed

+105
-8
lines changed

‎PWGHF/HFL/Tasks/taskElectronWeakBoson.cxx‎

Lines changed: 105 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525
#include "Common/Core/ZorroSummary.h"
2626
#include "Common/DataModel/Centrality.h"
2727
#include "Common/DataModel/EventSelection.h"
28+
#include "Common/DataModel/Multiplicity.h"
2829
#include "Common/DataModel/PIDResponseTPC.h"
2930
#include "Common/DataModel/TrackSelectionTables.h"
3031
#include "Tools/KFparticle/KFUtilities.h"
@@ -46,6 +47,7 @@
4647
#include <Framework/runDataProcessing.h>
4748

4849
#include <TPDGCode.h>
50+
#include <TRandom3.h>
4951

5052
#include <KFPTrack.h>
5153
#include <KFParticle.h>
@@ -130,6 +132,8 @@ struct HfTaskElectronWeakBoson {
130132
// Centrality estimator configuration
131133
Configurable<int> centralityEstimator{"centralityEstimator", CentralityEstimator::FT0M, "Centrality estimator. See CentralityEstimator for valid values."};
132134
Configurable<bool> enableCentralityAnalysis{"enableCentralityAnalysis", true, "Enable centrality-dependent analysis"};
135+
Configurable<bool> enableMultiplicityPVAnalysis{"enableMultiplicityPVAnalysis", true, "Enable centrality-dependent analysis"};
136+
Configurable<bool> enableMultiplicityFT0MAnalysis{"enableMultiplicityFT0MAnalysis", true, "Enable centrality-dependent analysis"};
133137
Configurable<float> centralityMin{"centralityMin", -1, "minimum cut on centrality selection"};
134138
Configurable<float> centralityMax{"centralityMax", 101, "maximum cut on centrality selection"};
135139
Configurable<std::vector<double>> centralityBins{"centralityBins", {0, 20, 60, 100}, "centrality bins"};
@@ -139,6 +143,13 @@ struct HfTaskElectronWeakBoson {
139143
Configurable<float> massZMinQA{"massZMinQA", 0.1, "minimum mass cut for Zee Reco QA"};
140144
// CCDB service object
141145
Service<o2::ccdb::BasicCCDBManager> ccdb{};
146+
// UE
147+
Configurable<int> nRandomCones{"nRandomCones", 100, "number of random cones"};
148+
Configurable<float> rcHardE{"rcHardE", 5.0, "hard cluster veto energy"};
149+
Configurable<float> rcVetoR{"rcVetoR", 0.4, "veto radius"};
150+
Configurable<bool> useUEsub{"useUEsub", true, "apply UE subtraction in isolation"};
151+
152+
TRandom3* rnd = new TRandom3(0);
142153

143154
struct HfElectronCandidate {
144155
float pt, eta, phi, dcaxyTrk, dcazTrk, eop, energyIso, momIso;
@@ -157,7 +168,7 @@ struct HfTaskElectronWeakBoson {
157168
: pt(ptr), eta(e), phi(ph), mass(m), ptchild0(ptzee0), ptchild1(ptzee1), charge(ch) {}
158169
};
159170
std::vector<HfZeeCandidate> reconstructedZ;
160-
using CollisionsWithCent = soa::Join<aod::Collisions, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms, aod::CentFV0As>;
171+
using CollisionsWithCent = soa::Join<aod::Collisions, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms, aod::CentFV0As, aod::Mults>;
161172
using SelectedClusters = o2::aod::EMCALClusters;
162173
// PbPb
163174
// using TrackEle = o2::soa::Join<o2::aod::Tracks, o2::aod::FullTracks, o2::aod::TracksExtra, o2::aod::TracksDCA, o2::aod::TrackSelection, o2::aod::pidTPCFullEl>;
@@ -200,6 +211,8 @@ struct HfTaskElectronWeakBoson {
200211
ConfigurableAxis confaxisInvMassZgamma{"confaxisInvMassZgamma", {150, 0, 150}, "M_{ee} (GeV/c^{2})"};
201212
ConfigurableAxis confaxisInvMassZ{"confaxisInvMassZ", {130, 20, 150}, "M_{ee} (GeV/c^{2})"};
202213
ConfigurableAxis confaxisZfrag{"confaxisZfrag", {200, 0, 2.0}, "p_{T,h}/p_{T,Z}"};
214+
ConfigurableAxis confaxisMultPV{"confaxisMultPV", {200, 0, 200.0}, "multiplicity"};
215+
ConfigurableAxis confaxisMultFT0{"confaxisMultFT0", {1000, 0, 1000.0}, "multiplicity"};
203216

204217
// Histogram registry: an object to hold your registrygrams
205218
HistogramRegistry registry{"registry"};
@@ -281,12 +294,18 @@ struct HfTaskElectronWeakBoson {
281294
const AxisSpec axisInvMassZgamma{confaxisInvMassZgamma, "M_{ee} (GeV/c^{2})"};
282295
const AxisSpec axisInvMassZ{confaxisInvMassZ, "M_{ee} (GeV/c^{2})"};
283296
const AxisSpec axisZfrag{confaxisZfrag, "p_{T,h}/p_{T,Z}"};
297+
const AxisSpec axisMultPV{confaxisMultPV, "multiplicity"};
298+
const AxisSpec axisMultFT0{confaxisMultFT0, "multiplicity"};
284299

285300
// create registrygrams
286301
registry.add("hZvtx", "Z vertex", kTH1D, {axisZvtx});
287302
registry.add("hEventCounterInit", "hEventCounterInit", kTH1D, {axisCounter});
288303
registry.add("hEventCounter", "hEventCounter", kTH1D, {axisCounter});
289304
registry.add("hCentrality", "Centrality distribution", kTH1D, {axisCentrality});
305+
registry.add("hCentMultCorr", "Centrality distribution", kTH2D, {{axisCentrality}, {axisMultFT0}});
306+
registry.add("hMultPV", "multiplicity distribution for PV", kTH1D, {axisMultPV});
307+
registry.add("hMultFT0", "multiplicity distribution for FT0", kTH1D, {axisMultFT0});
308+
registry.add("hMultFT0PV", "multiplicity distribution", kTH2D, {{axisMultFT0}, {axisMultPV}});
290309
registry.add("hITSchi2", "ITS #chi^{2}", kTH1F, {axisChi2});
291310
registry.add("hTPCchi2", "TPC #chi^{2}", kTH1F, {axisChi2});
292311
registry.add("hTPCnCls", "TPC NCls", kTH1F, {axisCluster});
@@ -320,12 +339,19 @@ struct HfTaskElectronWeakBoson {
320339

321340
// hisotgram for EMCal trigger
322341
registry.add("hEMCalTrigger", "EMCal trigger", kTH1D, {axisTrigger});
342+
343+
// histogram for UE
344+
registry.add("hRho", "rho UE density", kTH1F, {axisE});
345+
registry.add("hSumERC", "RC sumE", kTH1F, {axisE});
346+
registry.add("hUE", "UE vs. centrality", kTH2F, {{axisCentrality}, {axisE}});
323347
}
324348

325349
double getIsolatedCluster(const o2::aod::EMCALCluster& cluster,
326-
const SelectedClusters& clusters)
350+
const SelectedClusters& clusters,
351+
float UE)
327352
{
328353
double energySum = 0.0;
354+
double energySum_excl = 0.0;
329355
double isoEnergy = 10.0;
330356
double const etaAssCluster = cluster.eta();
331357
double const phiAssCluster = cluster.phi();
@@ -346,11 +372,12 @@ struct HfTaskElectronWeakBoson {
346372
energySum += associateCluster.energy();
347373
}
348374
}
349-
375+
energySum_excl = energySum - cluster.energy();
350376
if (energySum > 0) {
351-
isoEnergy = energySum / cluster.energy() - 1.0;
377+
isoEnergy = (energySum_excl - UE) / cluster.energy();
352378
}
353379

380+
// LOG(info) <<"clustE = " << cluster.energy() << " ; energySum = " << energySum << " ; nclust in Cone = " << nclustSum - 1 << " ; UE = " << UE << " ; isoEnergy = " << isoEnergy;
354381
registry.fill(HIST("hIsolationEnergy"), cluster.energy(), isoEnergy);
355382

356383
return (isoEnergy);
@@ -387,6 +414,50 @@ struct HfTaskElectronWeakBoson {
387414
// LOG(info) << "isop = " << isoMomentum;
388415
return std::make_pair(trackCount - 1, isoMomentum);
389416
}
417+
float estimateRhoRC(const SelectedClusters& clusters)
418+
{
419+
const float R = rIsolation;
420+
const float A = o2::constants::math::PI * R * R;
421+
422+
std::vector<float> sumErc;
423+
sumErc.reserve(nRandomCones);
424+
425+
for (int i = 0; i < nRandomCones; i++) {
426+
427+
float etarc = rnd->Uniform(-etaEmcMax, etaEmcMax); // in EMCal acceptance
428+
float phirc = rnd->Uniform(phiEmcMin, phiEmcMax); // in EMCal acceptance
429+
430+
float s = 0;
431+
432+
for (auto& c : clusters) {
433+
if (c.energy() > rcHardE)
434+
continue;
435+
double dEtarc = etarc - c.eta();
436+
double dPhirc = phirc - c.phi();
437+
dPhirc = RecoDecay::constrainAngle(dPhirc, -o2::constants::math::PI);
438+
double const deltaRrc = std::sqrt(dEtarc * dEtarc + dPhirc * dPhirc);
439+
if (deltaRrc < R) {
440+
s += c.energy();
441+
}
442+
}
443+
444+
registry.fill(HIST("hSumERC"), s);
445+
sumErc.push_back(s);
446+
}
447+
448+
if (sumErc.empty())
449+
return 0;
450+
451+
std::nth_element(sumErc.begin(),
452+
sumErc.begin() + sumErc.size() / 2,
453+
sumErc.end());
454+
455+
float median = sumErc[sumErc.size() / 2];
456+
// LOG(info) << "median = " << median;
457+
registry.fill(HIST("hRho"), median / A);
458+
459+
return median / A;
460+
}
390461

391462
void recoMassZee(const KFParticle& kfpIsoEle,
392463
int charge,
@@ -527,13 +598,39 @@ struct HfTaskElectronWeakBoson {
527598
float centrality = 1.0;
528599
if (enableCentralityAnalysis) {
529600
centrality = o2::hf_centrality::getCentralityColl(collision, centralityEstimator);
530-
// LOG(info) << centrality;
601+
// LOG(info) << "centrality = " << o2::hf_centrality::getCentralityColl(collision, centralityEstimator) << " ; FTC = " << collision.multFT0C();
531602
if (centrality < centralityMin || centrality > centralityMax) {
532603
return;
533604
}
534-
registry.fill(HIST("hCentrality"), centrality);
605+
registry.fill(HIST("hCentMultCorr"), centrality, collision.multFT0M());
606+
}
607+
608+
if (enableMultiplicityFT0MAnalysis || enableMultiplicityPVAnalysis) {
609+
if (enableMultiplicityFT0MAnalysis)
610+
centrality = collision.multFT0M();
611+
if (enableMultiplicityPVAnalysis)
612+
centrality = collision.multNTracksPV();
613+
// LOG(info) << "raw mult PV = " << collision.multNTracksPV();
614+
// LOG(info) << "raw mult FT0M = " << collision.multFT0M();
615+
registry.fill(HIST("hMultPV"), collision.multNTracksPV());
616+
registry.fill(HIST("hMultFT0"), collision.multFT0M());
617+
registry.fill(HIST("hMultFT0PV"), collision.multFT0M(), collision.multNTracksPV());
618+
}
619+
620+
registry.fill(HIST("hCentrality"), centrality);
621+
622+
// UE estimate
623+
float rho = 0.f;
624+
float UE = 0.f;
625+
626+
if (useUEsub) {
627+
rho = estimateRhoRC(emcClusters);
628+
UE = rho * (float)(o2::constants::math::PI * rIsolation * rIsolation);
629+
registry.fill(HIST("hUE"), centrality, UE);
630+
// LOG(info) << "UE = " << UE;
535631
}
536632

633+
// track loop
537634
for (const auto& track : tracks) {
538635

539636
if (std::abs(track.eta()) > etaTrMax) {
@@ -568,7 +665,7 @@ struct HfTaskElectronWeakBoson {
568665
registry.fill(HIST("hTPCNsigma"), track.p(), track.tpcNSigmaEl());
569666

570667
float eop = -0.01;
571-
float isoEnergy = 1.0;
668+
float isoEnergy = 99.0;
572669
// track isolation
573670
auto [trackCount, isoMomentum] = getIsolatedTrack(track.eta(), track.phi(), track.p(), tracks);
574671
// LOG(info) << "isoMomentum = " << isoMomentum;
@@ -657,7 +754,7 @@ struct HfTaskElectronWeakBoson {
657754
eop = energyEmc / match.track_as<TrackEle>().p();
658755
// LOG(info) << "eop = " << eop;
659756

660-
isoEnergy = getIsolatedCluster(cluster, emcClusters);
757+
isoEnergy = getIsolatedCluster(cluster, emcClusters, UE);
661758

662759
if (match.track_as<TrackEle>().pt() > ptTHnThresh && isTHnElectron) {
663760
registry.fill(HIST("hTHnElectrons"), match.track_as<TrackEle>().pt(), match.track_as<TrackEle>().tpcNSigmaEl(), m02Emc, eop, isoEnergy, isoMomentum, trackCount, track.eta(), track.tpcSignal());

0 commit comments

Comments
 (0)