Skip to content

Commit e101263

Browse files
peressounkoperessounko
peressounko
authored and
peressounko
committed
event selection fixes
1 parent 09462de commit e101263

File tree

1 file changed

+108
-81
lines changed

1 file changed

+108
-81
lines changed

‎PWGEM/Tasks/phosPi0.cxx‎

Lines changed: 108 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -14,34 +14,32 @@
1414
/// \author Dmitri Peresunko <Dmitri.Peresunko@cern.ch>
1515
///
1616

17-
#include <algorithm>
18-
#include <climits>
19-
#include <cstdlib>
20-
#include <map>
21-
#include <memory>
22-
#include <string>
23-
#include <vector>
24-
2517
#include "Common/DataModel/CaloClusters.h"
2618
#include "Common/DataModel/Centrality.h"
2719
#include "Common/DataModel/EventSelection.h"
2820
#include "Common/DataModel/Multiplicity.h"
29-
30-
#include "Framework/ConfigParamSpec.h"
31-
#include "Framework/runDataProcessing.h"
32-
#include "Framework/AnalysisTask.h"
33-
#include "Framework/AnalysisDataModel.h"
34-
#include "Framework/ASoA.h"
35-
#include "Framework/ASoAHelpers.h"
36-
#include "Framework/HistogramRegistry.h"
37-
3821
#include "EventFiltering/Zorro.h"
3922
#include "EventFiltering/ZorroSummary.h"
4023

41-
#include "PHOSBase/Geometry.h"
42-
#include "CommonDataFormat/InteractionRecord.h"
4324
#include "CCDB/BasicCCDBManager.h"
25+
#include "CommonDataFormat/InteractionRecord.h"
4426
#include "DataFormatsParameters/GRPLHCIFData.h"
27+
#include "Framework/ASoA.h"
28+
#include "Framework/ASoAHelpers.h"
29+
#include "Framework/AnalysisDataModel.h"
30+
#include "Framework/AnalysisTask.h"
31+
#include "Framework/ConfigParamSpec.h"
32+
#include "Framework/HistogramRegistry.h"
33+
#include "Framework/runDataProcessing.h"
34+
#include "PHOSBase/Geometry.h"
35+
36+
#include <algorithm>
37+
#include <climits>
38+
#include <cstdlib>
39+
#include <map>
40+
#include <memory>
41+
#include <string>
42+
#include <vector>
4543

4644
using namespace o2;
4745
using namespace o2::aod::evsel;
@@ -52,7 +50,7 @@ struct PhosPi0 {
5250
Configurable<bool> skimmedProcessing{"skimmedProcessing", false, "Skimmed dataset processing"};
5351
Configurable<std::string> trigName{"trigName", "fPHOSPhoton", "name of offline trigger"};
5452
Configurable<std::string> zorroCCDBpath{"zorroCCDBpath", "/Users/m/mpuccio/EventFiltering/OTS/", "path to the zorro ccdb objects"};
55-
Configurable<int> evSelTrig{"evSelTrig", aod::evsel::kIsTriggerTVX, "Select events with this trigger"};
53+
Configurable<int> evSelTrig{"evSelTrig", kTVXinPHOS, "Select events with this trigger"};
5654
Configurable<bool> isMC{"isMC", false, "to fill MC histograms"};
5755
Configurable<float> minCluE{"minCluE", 0.3, "Minimum cluster energy for analysis"};
5856
Configurable<float> minCluTime{"minCluTime", -25.e-9, "Min. cluster time"};
@@ -107,14 +105,19 @@ struct PhosPi0 {
107105
int mRunNumber = 0; // Current run number
108106
int mixedEventBin = 0; // Which list of Mixed use for mixing
109107
std::vector<Photon> mCurEvent;
110-
static constexpr int kMaxMixBins = 20; // maximal number of kinds of events for mixing
108+
static constexpr int kMixBinsZ = 20;
109+
static constexpr int kMixBinsPhi = 6;
110+
static constexpr int kMaxMixBins = kMixBinsZ * kMixBinsPhi + 1; // maximal number of bins of events for mixing: vtx*event plane + hard
111+
static constexpr double kEmixCut = 10.; // minimal clu energy for special hard mixing bin
112+
static constexpr int kHardMixBin = kMaxMixBins - 1; // special hard mixing bin
113+
111114
std::array<std::deque<std::vector<Photon>>, kMaxMixBins> mMixedEvents;
112115
std::array<std::deque<std::vector<Photon>>, kMaxMixBins> mAmbMixedEvents;
113116

114117
int mPrevMCColId = -1; // mark MC collissions already scanned
115118
// fast access to histos
116119
TH1* hColl;
117-
TH3 *hReMod, *hMiMod;
120+
TH3 *hReMod, *hMiMod, *hReAsym, *hMiAsym;
118121
TH2 *hReAll, *hReDisp, *hReCPV, *hReBoth, *hSignalAll, *hPi0SignalAll, *hPi0SignalCPV, *hPi0SignalDisp,
119122
*hPi0SignalBoth, *hMiAll, *hMiDisp, *hMiCPV, *hMiBoth;
120123
TH2 *hReOneAll, *hReOneDisp, *hReOneCPV, *hReOneBoth, *hMiOneAll, *hMiOneDisp, *hMiOneCPV, *hMiOneBoth;
@@ -154,8 +157,8 @@ struct PhosPi0 {
154157
hColl->GetXaxis()->SetBinLabel(4, "T0a&&T0c");
155158
hColl->GetXaxis()->SetBinLabel(5, "kTVXinPHOS");
156159
hColl->GetXaxis()->SetBinLabel(6, "kIsTriggerTVX");
157-
hColl->GetXaxis()->SetBinLabel(7, "PHOSClu");
158-
hColl->GetXaxis()->SetBinLabel(8, "PHOSClu&&kTVXinPHOS");
160+
hColl->GetXaxis()->SetBinLabel(7, "kPHOS");
161+
hColl->GetXaxis()->SetBinLabel(8, "kPHOS&&kTVXinPHOS");
159162
hColl->GetXaxis()->SetBinLabel(9, "Accepted");
160163

161164
auto h2{std::get<std::shared_ptr<TH1>>(mHistManager.add("eventsBC", "Number of events per trigger", HistType::kTH1F, {{8, 0., 8.}}))};
@@ -186,6 +189,7 @@ struct PhosPi0 {
186189
mHistManager.add("cluCPVOcc", "Cluster with CPV occupancy", HistType::kTH3F, {phiAxis, zAxis, modAxis});
187190
mHistManager.add("cluDispOcc", "Cluster with Disp occupancy", HistType::kTH3F, {phiAxis, zAxis, modAxis});
188191
mHistManager.add("cluBothOcc", "Cluster with Both occupancy", HistType::kTH3F, {phiAxis, zAxis, modAxis});
192+
mHistManager.add("qvec", "qvector", HistType::kTH2F, {{10, 0, o2::constants::math::PI, "#phi", "#phi (rad)"}, {10, 0., 1., "|q|", "|q|"}});
189193
}
190194

191195
hReMod = std::get<std::shared_ptr<TH3>>(mHistManager.add("mggReModComb", "inv mass per module",
@@ -203,6 +207,9 @@ struct PhosPi0 {
203207
hReBoth = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReBoth", "inv mass for centrality",
204208
HistType::kTH2F, {mggAxis, ptAxis}))
205209
.get();
210+
hReAsym = std::get<std::shared_ptr<TH3>>(mHistManager.add("mggReAsym", "inv mass vs pt vs asym",
211+
HistType::kTH3F, {mggAxis, ptAxis, {10, 0., 1., "asym", "a"}}))
212+
.get();
206213
hReOneAll = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReOneAll", "inv mass for centrality",
207214
HistType::kTH2F, {mggAxis, ptAxis}))
208215
.get();
@@ -262,6 +269,9 @@ struct PhosPi0 {
262269
hMiBoth = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggMiBoth", "inv mass for centrality",
263270
HistType::kTH2F, {mggAxis, ptAxis}))
264271
.get();
272+
hMiAsym = std::get<std::shared_ptr<TH3>>(mHistManager.add("mggMiAsym", "inv mass vs pt vs asym",
273+
HistType::kTH3F, {mggAxis, ptAxis, {10, 0., 1., "asym", "a"}}))
274+
.get();
265275
hMiOneAll = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggMiOneAll", "inv mass for centrality",
266276
HistType::kTH2F, {mggAxis, ptAxis}))
267277
.get();
@@ -306,25 +316,28 @@ struct PhosPi0 {
306316
/// \brief Process PHOS data
307317
void processData(SelCollisions::iterator const& col,
308318
aod::CaloClusters const& clusters,
319+
aod::FullTracks const& tracks,
309320
aod::BCsWithTimestamps const&)
310321
{
311322
aod::McParticles const* mcPart = nullptr;
312-
scanAll<false>(col, clusters, mcPart);
323+
scanAll<false>(col, clusters, tracks, mcPart);
313324
}
314325
PROCESS_SWITCH(PhosPi0, processData, "processData", true);
315326
void processMC(SelCollisionsMC::iterator const& col,
316327
McClusters const& clusters,
328+
aod::FullTracks const& tracks,
317329
aod::McParticles const& mcPart,
318330
aod::McCollisions const& /*mcCol*/,
319331
aod::BCsWithTimestamps const&)
320332
{
321-
scanAll<true>(col, clusters, &mcPart);
333+
scanAll<true>(col, clusters, tracks, &mcPart);
322334
}
323335
PROCESS_SWITCH(PhosPi0, processMC, "processMC", false);
324336

325337
template <bool isMC, typename TCollision, typename TClusters>
326338
void scanAll(TCollision& col,
327339
TClusters& clusters,
340+
aod::FullTracks const& tracks,
328341
aod::McParticles const* mcPart)
329342
{
330343
mixedEventBin = 0;
@@ -341,7 +354,7 @@ struct PhosPi0 {
341354
return; ///
342355
}
343356
} else {
344-
if (!col.selection_bit(evSelTrig)) {
357+
if (!col.alias_bit(evSelTrig)) {
345358
return;
346359
}
347360
}
@@ -350,13 +363,14 @@ struct PhosPi0 {
350363
const double vtxCut = 10.;
351364
double vtxZ = col.posZ();
352365
mHistManager.fill(HIST("vertex"), vtxZ);
353-
bool isColSelected = false;
354-
if constexpr (isMC) {
355-
isColSelected = (col.selection_bit(kIsTriggerTVX) && (clusters.size() > 0));
356-
} else {
357-
isColSelected = col.selection_bit(evSelTrig) && std::abs(vtxZ) < vtxCut; // col.alias_bit(evSelTrig)
358-
// collision.selection_bit(aod::evsel::kNoTimeFrameBorder);
359-
}
366+
bool isColSelected = col.alias_bit(evSelTrig) && std::abs(vtxZ) < vtxCut;
367+
;
368+
// if constexpr (isMC) {
369+
// isColSelected = (col.selection_bit(kIsTriggerTVX) && (clusters.size() > 0));
370+
// } else {
371+
// isColSelected = col.alias_bit(evSelTrig) && std::abs(vtxZ) < vtxCut; //
372+
// // collision.selection_bit(aod::evsel::kNoTimeFrameBorder);
373+
// }
360374

361375
if (col.selection_bit(kIsBBT0A) || col.selection_bit(kIsBBT0C)) {
362376
hColl->Fill(2.5);
@@ -367,55 +381,28 @@ struct PhosPi0 {
367381
if (col.alias_bit(kTVXinPHOS)) {
368382
hColl->Fill(4.5);
369383
}
370-
if (col.selection_bit(kIsTriggerTVX)) {
384+
if (col.alias_bit(kIsTriggerTVX)) {
371385
hColl->Fill(5.5);
372386
}
373-
if (clusters.size() > 0) {
387+
if (col.alias_bit(kPHOS)) {
374388
hColl->Fill(6.5);
375389
if (col.alias_bit(kTVXinPHOS)) {
376390
hColl->Fill(7.5);
377391
}
378392
}
379-
// //Event Plane| jet orientation
380-
// if (flag & (kProton | kDeuteron | kTriton | kHe3 | kHe4) || doprocessMC) { /// ignore PID pre-selections for the MC
381-
// if constexpr (std::is_same<Tcoll, CollWithEP>::value) {
382-
// nuclei::candidates_flow.emplace_back(NucleusCandidateFlow{
383-
// collision.centFV0A(),
384-
// collision.centFT0M(),
385-
// collision.centFT0A(),
386-
// collision.centFT0C(),
387-
// collision.psiFT0A(),
388-
// collision.multFT0A(),
389-
// collision.psiFT0C(),
390-
// collision.multFT0C(),
391-
// collision.psiTPC(),
392-
// collision.psiTPCL(),
393-
// collision.psiTPCR(),
394-
// collision.multTPC()});
395-
// } else if constexpr (std::is_same<Tcoll, CollWithQvec>::value) {
396-
// nuclei::candidates_flow.emplace_back(NucleusCandidateFlow{
397-
// collision.centFV0A(),
398-
// collision.centFT0M(),
399-
// collision.centFT0A(),
400-
// collision.centFT0C(),
401-
// 0.5 * std::atan2(collision.qvecFT0AIm(), collision.qvecFT0ARe()),
402-
// collision.multFT0A(),
403-
// 0.5 * std::atan2(collision.qvecFT0CIm(), collision.qvecFT0CRe()),
404-
// collision.multFT0C(),
405-
// -999.,
406-
// 0.5 * std::atan2(collision.qvecBNegIm(), collision.qvecBNegRe()),
407-
// 0.5 * std::atan2(collision.qvecBPosIm(), collision.qvecBPosRe()),
408-
// collision.multTPC()});
409-
// }
410-
411-
int mult = 1.; // multiplicity TODO!!!
412-
mixedEventBin = findMixedEventBin(vtxZ, mult);
413393

414394
if (!isColSelected) {
415395
return;
416396
}
417397
hColl->Fill(8.5);
418398

399+
int mult = 1.; // multiplicity TODO!!!
400+
std::pair<double, double> q = evalQvec(tracks);
401+
402+
// find proper event for mixing
403+
// note, that if one find hard cluster in PHOS later, it will change bin to special one
404+
mixedEventBin = findMixedEventBin(vtxZ, mult, q);
405+
419406
// Fill MC distributions
420407
// pion rapidity, pt, phi
421408
// secondary pi0s
@@ -482,6 +469,9 @@ struct PhosPi0 {
482469
clu.m02() < minM02) {
483470
continue;
484471
}
472+
if (clu.e() > kEmixCut) {
473+
mixedEventBin = kHardMixBin;
474+
}
485475
if (fillQC) {
486476
mHistManager.fill(HIST("cluSp"), clu.e(), clu.mod());
487477
if (clu.e() > minOccE) {
@@ -529,6 +519,7 @@ struct PhosPi0 {
529519
}
530520
hReMod->Fill(m, pt, modComb, w);
531521
hReAll->Fill(m, pt, w);
522+
hReAsym->Fill(m, pt, std::abs((ph1.e - ph2.e) / (ph1.e + ph2.e)), w);
532523
hReOneAll->Fill(m, ph1.pt(), w);
533524
hReOneAll->Fill(m, ph2.pt(), w);
534525
if (ph1.isCPVOK()) {
@@ -633,6 +624,7 @@ struct PhosPi0 {
633624
}
634625
hMiMod->Fill(m, pt, modComb, w);
635626
hMiAll->Fill(m, pt, w);
627+
hMiAsym->Fill(m, pt, std::abs((ph1.e - ph2.e) / (ph1.e + ph2.e)), w);
636628
hMiOneAll->Fill(m, ph1.pt(), w);
637629
hMiOneAll->Fill(m, ph2.pt(), w);
638630
if (ph1.isCPVOK()) {
@@ -681,9 +673,10 @@ struct PhosPi0 {
681673
mixedEventBin = 0;
682674

683675
mHistManager.fill(HIST("eventsBC"), 0.);
684-
double vtxZ = 0; // no vtx info
685-
int mult = 1.; // multiplicity TODO!!!
686-
mixedEventBin = findMixedEventBin(vtxZ, mult);
676+
double vtxZ = 0; // no vtx info
677+
int mult = 1.; // multiplicity TODO!!!
678+
std::pair<double, double> q{0, 0}; // fake q-vector as no tracks
679+
mixedEventBin = findMixedEventBin(vtxZ, mult, q);
687680

688681
if (!isSelected) {
689682
return;
@@ -831,17 +824,27 @@ struct PhosPi0 {
831824
4.;
832825
}
833826
//_____________________________________________________________________________
834-
int findMixedEventBin(double vtxZ, double /*mult */)
827+
int findMixedEventBin(double vtxZ, double /*mult*/, std::pair<double, double>& q)
835828
{
836829
// calculate index for event mixing
837830
const double zwidth = 1.; // Width of zvtx bin
838-
int res = static_cast<int>((vtxZ + 10.) / zwidth);
839-
840-
if (res < 0)
841-
return 0;
842-
if (res >= kMaxMixBins)
843-
return kMaxMixBins - 1;
844-
return res;
831+
int iz = static_cast<int>((vtxZ + 10.) / zwidth);
832+
833+
if (iz < 0)
834+
iz = 0;
835+
if (iz >= kMixBinsZ)
836+
iz = kMixBinsZ - 1;
837+
838+
// event plane orientation
839+
double phi = 0.5 * std::atan2(q.second, q.first); // 1/2 due to second order flow harmonic
840+
while (phi < 0)
841+
phi += o2::constants::math::PI;
842+
while (phi > o2::constants::math::PI)
843+
phi -= o2::constants::math::PI;
844+
int iphi = static_cast<int>(kMixBinsPhi * phi / o2::constants::math::PI);
845+
mHistManager.fill(HIST("qvec"), phi, std::sqrt(q.first * q.first + q.second * q.second));
846+
847+
return iz * iphi;
845848
}
846849
//----------------------------------------
847850
int commonParentPDG(int lab1, int lab2, aod::McParticles const* mcParticles)
@@ -907,6 +910,30 @@ struct PhosPi0 {
907910
}
908911
return 1.;
909912
}
913+
//_____________________________________________________________________________
914+
std::pair<double, double> evalQvec(aod::FullTracks const& tracks)
915+
{
916+
// calculate approximate q-vector for event
917+
const int ord = 2; // flow order
918+
std::pair<double, double> q{0, 0};
919+
int ntr = 0;
920+
for (const auto& track : tracks) {
921+
if (!track.has_collision()) { // ignore orphan tracks without collision
922+
continue;
923+
}
924+
// if (!track.isGlobalTrack()) { // only global tracks
925+
// continue;
926+
// }
927+
q.first += std::cos(ord * track.phi());
928+
q.second += std::sin(ord * track.phi());
929+
ntr++;
930+
}
931+
if (ntr > 0) {
932+
q.first /= ntr;
933+
q.second /= ntr;
934+
}
935+
return q;
936+
}
910937
};
911938

912939
o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc)

0 commit comments

Comments
 (0)