Skip to content

Commit f1bd16f

Browse files
bhagyasahoalibuild
bhagyasaho
and
alibuild
authored
[PWGDQ] Adding muon DCA Task (#7062)
* Add files via upload * Update ReducedInfoTables.h * [PWGDQ] Update muonDCA.cxx * Update muonDCA.cxx * Please consider the following formatting changes --------- Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent d0cde60 commit f1bd16f

2 files changed

Lines changed: 201 additions & 0 deletions

File tree

‎PWGDQ/Tasks/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -69,6 +69,11 @@ o2physics_add_dpl_workflow(task-j-psi-hf
6969
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore
7070
COMPONENT_NAME Analysis)
7171

72+
o2physics_add_dpl_workflow(task-muon-dca
73+
SOURCES muonDCA.cxx
74+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore
75+
COMPONENT_NAME Analysis)
76+
7277
o2physics_add_dpl_workflow(correlation
7378
SOURCES dqCorrelation.cxx
7479
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::GFWCore

‎PWGDQ/Tasks/muonDCA.cxx‎

Lines changed: 196 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,196 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
//
12+
/// \file muonDCA.cxx
13+
/// \brief Task to compute and evaluate DCA quantities
14+
/// \author Nicolas Bizé <nicolas.bize@cern.ch>, SUBATECH
15+
//
16+
#include "Framework/runDataProcessing.h"
17+
#include "Framework/AnalysisTask.h"
18+
#include "Framework/ASoAHelpers.h"
19+
#include "GlobalTracking/MatchGlobalFwd.h"
20+
#include "CCDB/BasicCCDBManager.h"
21+
#include "DataFormatsParameters/GRPMagField.h"
22+
#include "PWGDQ/DataModel/ReducedInfoTables.h"
23+
#include "PWGDQ/Core/VarManager.h"
24+
25+
using namespace o2;
26+
using namespace o2::framework;
27+
using namespace o2::aod;
28+
29+
using MyMuons = soa::Join<aod::ReducedMuons, aod::ReducedMuonsExtra, aod::ReducedMuonsCov>;
30+
using MyEvents = soa::Join<aod::ReducedEvents, aod::ReducedEventsExtended>;
31+
using MyEventsVtxCov = soa::Join<aod::ReducedEvents, aod::ReducedEventsExtended, aod::ReducedEventsVtxCov>;
32+
33+
// constexpr static uint32_t gkMuonDCAFillMapWithCov = VarManager::ObjTypes::ReducedMuon | VarManager::ObjTypes::ReducedMuonExtra | VarManager::ObjTypes::ReducedMuonCov | VarManager::ObjTypes::MuonDCA;
34+
35+
constexpr static int toVertex = VarManager::kToVertex;
36+
constexpr static int toDCA = VarManager::kToDCA;
37+
constexpr static int toRabs = VarManager::kToRabs;
38+
39+
static o2::globaltracking::MatchGlobalFwd mExtrap;
40+
template <typename T>
41+
bool isSelected(const T& muon);
42+
43+
struct muonExtrap {
44+
Produces<ReducedMuonsDca> dcaTable;
45+
Configurable<std::string> geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"};
46+
Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
47+
Configurable<string> fConfigCcdbUrl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
48+
49+
Service<o2::ccdb::BasicCCDBManager> fCCDB;
50+
o2::parameters::GRPMagField* grpmag = nullptr; // for run 3, we access GRPMagField from GLO/Config/GRPMagField
51+
int fCurrentRun; // needed to detect if the run changed and trigger update of magnetic field
52+
53+
HistogramRegistry registry{
54+
"registry",
55+
{}};
56+
57+
void init(o2::framework::InitContext&)
58+
{
59+
// Load geometry
60+
fCCDB->setURL(fConfigCcdbUrl);
61+
fCCDB->setCaching(true);
62+
fCCDB->setLocalObjectValidityChecking();
63+
64+
if (!o2::base::GeometryManager::isGeometryLoaded()) {
65+
LOGF(info, "Load geometry from CCDB");
66+
fCCDB->get<TGeoManager>(geoPath);
67+
}
68+
69+
AxisSpec pdcaAxis = {5000, 0.0, 5000.0, "p #times DCA"};
70+
AxisSpec dcaAxis = {200, 0.0, 200.0, "DCA"};
71+
AxisSpec dcaxAxis = {200, -100.0, 100.0, "DCA_x"};
72+
AxisSpec dcayAxis = {200, -100.0, 100.0, "DCA_y"};
73+
AxisSpec rabsAxis = {100, 0., 100.0, "R_{abs}"};
74+
AxisSpec xAxis = {200, -100., 100.0, "x (cm)"};
75+
AxisSpec yAxis = {200, -100., 100.0, "y (cm)"};
76+
AxisSpec zAxis = {200, -100., 100.0, "z (cm)"};
77+
78+
HistogramConfigSpec pdcaSpec({HistType::kTH1F, {pdcaAxis}});
79+
HistogramConfigSpec dcaSpec({HistType::kTH1F, {dcaAxis}});
80+
HistogramConfigSpec dcaxSpec({HistType::kTH1F, {dcaxAxis}});
81+
HistogramConfigSpec dcaySpec({HistType::kTH1F, {dcayAxis}});
82+
HistogramConfigSpec rabsSpec({HistType::kTH1F, {rabsAxis}});
83+
HistogramConfigSpec xSpec({HistType::kTH1F, {xAxis}});
84+
HistogramConfigSpec ySpec({HistType::kTH1F, {yAxis}});
85+
HistogramConfigSpec zSpec({HistType::kTH1F, {zAxis}});
86+
87+
registry.add("pdca", "pDCA", pdcaSpec);
88+
registry.add("dca", "DCA", dcaSpec);
89+
registry.add("dcax", "DCA_x", dcaxSpec);
90+
registry.add("dcay", "DCA_y", dcaySpec);
91+
registry.add("rabs", "R_{abs}", rabsSpec);
92+
registry.add("xAtVtx", "x at vertex", xSpec);
93+
registry.add("xAtDCA", "x at DCA", xSpec);
94+
registry.add("xAtRabs", "x at end abs", xSpec);
95+
registry.add("yAtVtx", "y at vertex", ySpec);
96+
registry.add("yAtDCA", "y at DCA", ySpec);
97+
registry.add("yAtRabs", "y at end abs", ySpec);
98+
registry.add("zAtVtx", "z at vertex", zSpec);
99+
registry.add("zAtDCA", "z at DCA", zSpec);
100+
registry.add("zAtRabs", "z at end abs", zSpec);
101+
}
102+
103+
void processExtrapolation(MyEventsVtxCov::iterator const& collision, MyMuons const& muons)
104+
{
105+
if (fCurrentRun != collision.runNumber()) {
106+
grpmag = fCCDB->getForTimeStamp<o2::parameters::GRPMagField>(grpmagPath, collision.timestamp());
107+
if (grpmag != nullptr) {
108+
LOGF(info, "Init field from GRP");
109+
o2::base::Propagator::initFieldFromGRP(grpmag);
110+
}
111+
LOGF(info, "Set field for muons");
112+
VarManager::SetupMuonMagField();
113+
fCurrentRun = collision.runNumber();
114+
}
115+
116+
for (auto& muon : muons) {
117+
if (static_cast<int>(muon.trackType()) < 2) {
118+
continue; // Make sure to remove global muon tracks
119+
}
120+
// propagate muon track to vertex
121+
o2::dataformats::GlobalFwdTrack muonTrackAtVertex = VarManager::PropagateMuon(muon, collision, toVertex);
122+
123+
// propagate muon track to DCA
124+
o2::dataformats::GlobalFwdTrack muonTrackAtDCA = VarManager::PropagateMuon(muon, collision, toDCA);
125+
126+
// propagate to Rabs
127+
o2::dataformats::GlobalFwdTrack muonTrackAtRabs = VarManager::PropagateMuon(muon, collision, toRabs);
128+
129+
// Calculate DCA quantities (preferable to do it with VarManager)
130+
double dcax = muonTrackAtDCA.getX() - collision.posX();
131+
double dcay = muonTrackAtDCA.getY() - collision.posY();
132+
double dca = std::sqrt(dcax * dcax + dcay * dcay);
133+
double pdca = muonTrackAtVertex.getP() * dca;
134+
double xAtVtx = muonTrackAtVertex.getX();
135+
double yAtVtx = muonTrackAtVertex.getY();
136+
double zAtVtx = muonTrackAtVertex.getZ();
137+
double xAtDCA = muonTrackAtDCA.getX();
138+
double yAtDCA = muonTrackAtDCA.getY();
139+
double zAtDCA = muonTrackAtDCA.getZ();
140+
double xAbs = muonTrackAtRabs.getX();
141+
double yAbs = muonTrackAtRabs.getY();
142+
double zAbs = muonTrackAtRabs.getZ();
143+
144+
double rabs = std::sqrt(xAbs * xAbs + yAbs * yAbs);
145+
146+
// QA histograms
147+
registry.get<TH1>(HIST("pdca"))->Fill(pdca);
148+
registry.get<TH1>(HIST("dca"))->Fill(dca);
149+
registry.get<TH1>(HIST("dcax"))->Fill(dcax);
150+
registry.get<TH1>(HIST("dcay"))->Fill(dcay);
151+
registry.get<TH1>(HIST("rabs"))->Fill(rabs);
152+
153+
registry.get<TH1>(HIST("xAtDCA"))->Fill(xAtDCA);
154+
registry.get<TH1>(HIST("xAtRabs"))->Fill(xAbs);
155+
registry.get<TH1>(HIST("xAtVtx"))->Fill(xAtVtx);
156+
157+
registry.get<TH1>(HIST("yAtDCA"))->Fill(yAtDCA);
158+
registry.get<TH1>(HIST("yAtRabs"))->Fill(yAbs);
159+
registry.get<TH1>(HIST("yAtVtx"))->Fill(yAtVtx);
160+
161+
registry.get<TH1>(HIST("zAtDCA"))->Fill(zAtDCA);
162+
registry.get<TH1>(HIST("zAtRabs"))->Fill(zAbs);
163+
registry.get<TH1>(HIST("zAtVtx"))->Fill(zAtVtx);
164+
165+
// Fill DCA table
166+
dcaTable(pdca,
167+
dca,
168+
dcax,
169+
dcay,
170+
rabs,
171+
muonTrackAtVertex.getPt(),
172+
muonTrackAtVertex.getEta(),
173+
muonTrackAtVertex.getPhi(),
174+
muon.sign(),
175+
muon.isAmbiguous(),
176+
muonTrackAtVertex.getPx(),
177+
muonTrackAtVertex.getPy(),
178+
muonTrackAtVertex.getPz());
179+
}
180+
}
181+
182+
PROCESS_SWITCH(muonExtrap, processExtrapolation, "process extrapolation", false);
183+
184+
void processDummy(MyEventsVtxCov&)
185+
{
186+
// do nothing
187+
}
188+
189+
PROCESS_SWITCH(muonExtrap, processDummy, "do nothing", false);
190+
};
191+
192+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
193+
{
194+
return WorkflowSpec{
195+
adaptAnalysisTask<muonExtrap>(cfgc)};
196+
};

0 commit comments

Comments
 (0)