Skip to content

Commit 3e46f29

Browse files
dsekihatdsekihat
authored
[PWGEM/Dilepton] add a task to test bremsstrahlung in MC (#14704)
1 parent 7836101 commit 3e46f29

File tree

2 files changed

+232
-0
lines changed

2 files changed

+232
-0
lines changed

‎PWGEM/Dilepton/Tasks/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,11 @@ o2physics_add_dpl_workflow(create-resolution-map
4040
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::GlobalTracking
4141
COMPONENT_NAME Analysis)
4242

43+
o2physics_add_dpl_workflow(test-bremsstrahlung
44+
SOURCES testBremsstrahlung.cxx
45+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
46+
COMPONENT_NAME Analysis)
47+
4348
o2physics_add_dpl_workflow(table-reader-barrel
4449
SOURCES tableReaderBarrel.cxx
4550
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::PWGDQCore

‎PWGEM/Dilepton/Tasks/testBremsstrahlung.cxx‎

Lines changed: 227 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,227 @@
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+
//
13+
// Analysis task to produce resolution mapfor electrons/muons in dilepton analysis
14+
// Please write to: daiki.sekihata@cern.ch
15+
16+
#include "PWGEM/Dilepton/Utils/MCUtilities.h"
17+
18+
#include "Common/CCDB/RCTSelectionFlags.h"
19+
#include "Common/Core/RecoDecay.h"
20+
#include "Common/Core/fwdtrackUtilities.h"
21+
#include "Common/Core/trackUtilities.h"
22+
#include "Common/DataModel/Centrality.h"
23+
#include "Common/DataModel/CollisionAssociationTables.h"
24+
#include "Common/DataModel/EventSelection.h"
25+
#include "Common/DataModel/TrackSelectionTables.h"
26+
27+
#include "CCDB/BasicCCDBManager.h"
28+
#include "DataFormatsCalibration/MeanVertexObject.h"
29+
#include "DataFormatsParameters/GRPMagField.h"
30+
#include "DetectorsBase/Propagator.h"
31+
#include "Field/MagneticField.h"
32+
#include "Framework/ASoA.h"
33+
#include "Framework/ASoAHelpers.h"
34+
#include "Framework/AnalysisDataModel.h"
35+
#include "Framework/AnalysisTask.h"
36+
#include "Framework/DataTypes.h"
37+
#include "Framework/HistogramRegistry.h"
38+
#include "Framework/runDataProcessing.h"
39+
40+
// #include "TGeoGlobalMagField.h"
41+
42+
#include <array>
43+
#include <map>
44+
#include <set>
45+
#include <string>
46+
#include <tuple>
47+
#include <unordered_map>
48+
#include <utility>
49+
#include <vector>
50+
51+
using namespace o2;
52+
using namespace o2::framework;
53+
using namespace o2::framework::expressions;
54+
using namespace o2::aod;
55+
using namespace o2::soa;
56+
using namespace o2::aod::pwgem::dilepton::utils::mcutil;
57+
58+
struct testBremsstrahlung {
59+
Configurable<std::string> ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
60+
Configurable<std::string> grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"};
61+
Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
62+
Configurable<std::string> geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"};
63+
Configurable<std::string> lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"};
64+
Configurable<std::string> mVtxPath{"mVtxPath", "GLO/Calib/MeanVertex", "Path of the mean vertex file"};
65+
Configurable<float> d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"};
66+
Configurable<bool> skipGRPOquery{"skipGRPOquery", true, "skip grpo query"};
67+
68+
HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject};
69+
70+
o2::ccdb::CcdbApi ccdbApi;
71+
Service<o2::ccdb::BasicCCDBManager> ccdb;
72+
int mRunNumber = 0;
73+
float d_bz = 0;
74+
o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT;
75+
o2::dataformats::VertexBase mVtx;
76+
const o2::dataformats::MeanVertexObject* mMeanVtx = nullptr;
77+
o2::base::MatLayerCylSet* lut = nullptr;
78+
79+
~testBremsstrahlung() {}
80+
81+
void init(o2::framework::InitContext&)
82+
{
83+
ccdb->setURL(ccdburl);
84+
ccdb->setCaching(true);
85+
ccdb->setLocalObjectValidityChecking();
86+
ccdb->setFatalWhenNull(false);
87+
ccdbApi.init(ccdburl);
88+
89+
mRunNumber = 0;
90+
d_bz = 0;
91+
92+
if (doprocessGen) {
93+
registry.add("Event/hGenID", "generator ID;generator ID;Number of mc collisions", kTH1F, {{7, -1.5, 5.5}}, false);
94+
registry.add("Photon/hXY", "prod. vtx. of #gamma^{bremss};X (cm);Y (cm)", kTH2D, {{400, -100, +100}, {400, -100, +100}}, false);
95+
registry.add("Photon/hs", "kinetic var. at prod. vtx.;p_{T,#gamma}^{bremss} (GeV/c);#eta_{#gamma}^{bremss};#varphi_{#gamma}^{bremss} (rad.);", kTHnSparseD, {{100, 0, 1}, {40, -1, +1}, {72, 0, 2 * M_PI}}, false);
96+
registry.add("Photon/hsDelta", "hsDelta;(p_{T,e}^{after} - p_{T,e}^{before})/p_{T,e}^{before};p_{T,e}^{before} (GeV/c);r_{xy}^{bremss} (cm);#varphi^{bremss} (rad.);z^{bremss} (cm);", kTHnSparseD, {{100, -1, 0}, {100, 0, 1}, {200, 0, 100}, {72, 0, 2 * M_PI}, {200, -100, +100}}, false);
97+
}
98+
}
99+
100+
// void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
101+
// {
102+
// if (mRunNumber == bc.runNumber()) {
103+
// return;
104+
// }
105+
106+
// // load matLUT for this timestamp
107+
// if (!lut) {
108+
// LOG(info) << "Loading material look-up table for timestamp: " << bc.timestamp();
109+
// lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForTimeStamp<o2::base::MatLayerCylSet>(lutPath, bc.timestamp()));
110+
// } else {
111+
// LOG(info) << "Material look-up table already in place. Not reloading.";
112+
// }
113+
114+
// // In case override, don't proceed, please - no CCDB access required
115+
// if (d_bz_input > -990) {
116+
// d_bz = d_bz_input;
117+
// o2::parameters::GRPMagField grpmag;
118+
// if (std::fabs(d_bz) > 1e-5) {
119+
// grpmag.setL3Current(30000.f / (d_bz / 5.0f));
120+
// }
121+
// o2::base::Propagator::initFieldFromGRP(&grpmag);
122+
// o2::base::Propagator::Instance()->setMatLUT(lut);
123+
// mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
124+
// mRunNumber = bc.runNumber();
125+
126+
// if (!o2::base::GeometryManager::isGeometryLoaded()) {
127+
// ccdb->get<TGeoManager>(geoPath);
128+
// }
129+
// }
130+
131+
// auto run3grp_timestamp = bc.timestamp();
132+
// o2::parameters::GRPObject* grpo = 0x0;
133+
// o2::parameters::GRPMagField* grpmag = 0x0;
134+
// if (!skipGRPOquery) {
135+
// grpo = ccdb->getForTimeStamp<o2::parameters::GRPObject>(grpPath, run3grp_timestamp);
136+
// }
137+
// if (grpo) {
138+
// o2::base::Propagator::initFieldFromGRP(grpo);
139+
// o2::base::Propagator::Instance()->setMatLUT(lut);
140+
// mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
141+
// // Fetch magnetic field from ccdb for current collision
142+
// d_bz = grpo->getNominalL3Field();
143+
// LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
144+
// } else {
145+
// grpmag = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(grpmagPath, run3grp_timestamp);
146+
// if (!grpmag) {
147+
// LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp;
148+
// }
149+
// o2::base::Propagator::initFieldFromGRP(grpmag);
150+
// o2::base::Propagator::Instance()->setMatLUT(lut);
151+
// mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
152+
153+
// // Fetch magnetic field from ccdb for current collision
154+
// d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f);
155+
// LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
156+
// }
157+
158+
// mRunNumber = bc.runNumber();
159+
// }
160+
161+
SliceCache cache;
162+
Preslice<aod::TracksIU> perCollision_mid = o2::aod::track::collisionId;
163+
Preslice<aod::McParticles> perMcCollision = aod::mcparticle::mcCollisionId;
164+
165+
using MyCollisions = Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs>;
166+
using MyCollision = MyCollisions::iterator;
167+
168+
using MyTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::McTrackLabels>;
169+
using MyTrack = MyTracks::iterator;
170+
171+
void processGen(aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles)
172+
{
173+
for (const auto& mcCollision : mcCollisions) {
174+
registry.fill(HIST("Event/hGenID"), mcCollision.getSubGeneratorId());
175+
176+
auto mcParticles_per_collision = mcParticles.sliceBy(perMcCollision, mcCollision.globalIndex());
177+
// auto photons_per_collision = photons->sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache);
178+
179+
for (const auto& mcParticle : mcParticles_per_collision) {
180+
if (std::abs(mcParticle.pdgCode()) != 22) {
181+
continue;
182+
}
183+
if (!mcParticle.has_mothers()) { // select seconday mcParticles
184+
continue;
185+
}
186+
if (mcParticle.isPhysicalPrimary() || mcParticle.producedByGenerator()) { // select seconday mcParticles
187+
continue;
188+
}
189+
190+
auto mother = mcParticle.template mothers_as<aod::McParticles>()[0]; // to be primary electron
191+
if (mother.eta() < -0.9 || 0.9 < mother.eta()) {
192+
continue;
193+
}
194+
if (std::abs(mother.pdgCode()) != 11) { // select bremsstrahlung photon from primary electron
195+
continue;
196+
}
197+
if (!(mother.isPhysicalPrimary() || mother.producedByGenerator())) { // select bremsstrahlung photon from primary electron
198+
continue;
199+
}
200+
201+
registry.fill(HIST("Photon/hXY"), mcParticle.vx(), mcParticle.vy());
202+
registry.fill(HIST("Photon/hs"), mcParticle.pt(), mcParticle.eta(), mcParticle.phi() > 0 ? mcParticle.phi() : mcParticle.phi() + 2 * M_PI);
203+
// int ndau = mother.daughtersIds()[1] - mother.daughtersIds()[0] + 1;
204+
// LOGF(info, "ndau = %d", ndau);
205+
206+
for (const auto& daughter : mother.daughters_as<aod::McParticles>()) {
207+
// LOGF(info, "mother.globalIndex() = %d, daughter.globalIndex() = %d, daughter.pdgCode() = %d", mother.globalIndex(), daughter.globalIndex(), daughter.pdgCode());
208+
if (std::abs(daughter.pdgCode()) == 11 && (!daughter.isPhysicalPrimary() && !daughter.producedByGenerator())) {
209+
float reldpt = (daughter.pt() - mother.pt()) / mother.pt();
210+
float rxy = std::sqrt(std::pow(daughter.vx(), 2) + std::pow(daughter.vy(), 2));
211+
float phi_bremss = RecoDecay::phi(daughter.vy(), daughter.vx());
212+
o2::math_utils::bringTo02Pi(phi_bremss);
213+
registry.fill(HIST("Photon/hsDelta"), reldpt, mother.pt(), rxy, phi_bremss, daughter.vz());
214+
}
215+
}
216+
217+
} // end of mc particle per mc collision
218+
219+
} // end of mc collision loop
220+
}
221+
PROCESS_SWITCH(testBremsstrahlung, processGen, "process generated info", true);
222+
};
223+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
224+
{
225+
return WorkflowSpec{
226+
adaptAnalysisTask<testBremsstrahlung>(cfgc, TaskName{"test-bremsstrahlung"})};
227+
}

0 commit comments

Comments
 (0)