Skip to content

Commit 84f40fb

Browse files
nicola-wilsonwirthnialibuild
authored
[PWGJE] Add diffusion wake analysis (#16570)
Co-authored-by: wirthni <74405269+wirthni@users.noreply.github.com> Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent ba4e0a5 commit 84f40fb

2 files changed

Lines changed: 286 additions & 0 deletions

File tree

‎PWGJE/TableProducer/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -110,3 +110,8 @@ o2physics_add_dpl_workflow(slim-tables-producer
110110
SOURCES slimTablesProducer.cxx
111111
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
112112
COMPONENT_NAME Analysis)
113+
114+
o2physics_add_dpl_workflow(diff-wake-tree-producer
115+
SOURCES diffWakeTreeProducer.cxx
116+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
117+
COMPONENT_NAME Analysis)

‎PWGJE/TableProducer/diffWakeTreeProducer.cxx‎

Lines changed: 281 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,281 @@
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 diffWakeTreeProducer.cxx
13+
/// \brief This task writes a collision and track table which are further used in a diffusion wake analysis
14+
/// \authors Nicola Wilson <nicola.wilson@cern.ch> and Nicolas Wirth <nicolas.wirth@cern.ch>
15+
16+
#include "Common/CCDB/EventSelectionParams.h"
17+
#include "Common/Core/EventPlaneHelper.h"
18+
#include "Common/DataModel/Centrality.h"
19+
#include "Common/DataModel/EventSelection.h"
20+
#include "Common/DataModel/Multiplicity.h"
21+
#include "Common/DataModel/Qvectors.h"
22+
#include "Common/DataModel/TrackSelectionTables.h"
23+
24+
#include <Framework/ASoA.h>
25+
#include <Framework/AnalysisDataModel.h>
26+
#include <Framework/AnalysisHelpers.h>
27+
#include <Framework/AnalysisTask.h>
28+
#include <Framework/HistogramRegistry.h>
29+
#include <Framework/HistogramSpec.h>
30+
#include <Framework/InitContext.h>
31+
#include <Framework/OutputObjHeader.h>
32+
#include <Framework/runDataProcessing.h>
33+
34+
// Event selection: Only events that contain track above some threshold
35+
// -------------------------------------------------------------------------------------------
36+
// TRACK DATA
37+
// -------------------------------------------------------------------------------------------
38+
// BEFORE COMPRESSION AFTER COMPRESSION
39+
// Name Data Type Size(b) Name Data Type Size(b)
40+
// ColID int32_t 4 [same]
41+
// Charge short 2 [same]
42+
// Px, Py, Pz float 3x4 P unsigned long 8
43+
// DEdx float 4 DEdx unsigned short 2
44+
// DCAXY float 4 DCAXY short 2
45+
// DCAZ float 4 DCAZ short 2
46+
// Length float 4 Length unsigned short 2
47+
48+
// OVERALL COMPRESSION 34b->22b
49+
50+
// -------------------------------------------------------------------------------------------
51+
// EVENT DATA
52+
// -------------------------------------------------------------------------------------------
53+
// GI int64_t 8 [same]
54+
// RN int 4 [same]
55+
// Cent float 4 [same]
56+
// Mult int 4 [same]
57+
// VertexX float 4 [same]
58+
// VertexY float 4 [same]
59+
// VertexZ float 4 [same]
60+
// Psi2 float 4 Psi2 short 2
61+
// Psi3 float 4 Psi3 short 2
62+
63+
// OVERALL COMPRESSION 40b->36b
64+
//--------------------------------------------------------
65+
namespace o2::aod
66+
{
67+
namespace testcol
68+
{
69+
// Event properties
70+
// DECLARE_SOA_COLUMN(Gi, gi, int64_t);
71+
DECLARE_SOA_COLUMN(Rn, rn, int32_t); // run number
72+
DECLARE_SOA_COLUMN(Cent, cent, float); // FT0C centrality
73+
DECLARE_SOA_COLUMN(Mult, mult, int32_t); // TPC multiplicity
74+
DECLARE_SOA_COLUMN(Occu, occu, int32_t); // Occupancy ITS
75+
DECLARE_SOA_COLUMN(Occuft0, occuft0, float); // Occupancy FT0C amplitudes
76+
DECLARE_SOA_COLUMN(VertexX, vertexX, float);
77+
DECLARE_SOA_COLUMN(VertexY, vertexY, float);
78+
DECLARE_SOA_COLUMN(VertexZ, vertexZ, float);
79+
DECLARE_SOA_COLUMN(Psi2, psi2, int16_t);
80+
DECLARE_SOA_COLUMN(Psi3, psi3, int16_t);
81+
} // namespace testcol
82+
83+
DECLARE_SOA_TABLE(TableCols, "AOD", "TABLECOL",
84+
o2::soa::Index<>,
85+
testcol::Rn,
86+
testcol::Cent,
87+
testcol::Mult,
88+
testcol::Occu,
89+
testcol::Occuft0,
90+
testcol::VertexX,
91+
testcol::VertexY,
92+
testcol::VertexZ,
93+
testcol::Psi2,
94+
testcol::Psi3);
95+
using TableCol = TableCols::iterator;
96+
97+
namespace testtrack
98+
{
99+
100+
// Track properties
101+
DECLARE_SOA_INDEX_COLUMN(TableCol, tableCol);
102+
DECLARE_SOA_COLUMN(Charge, charge, int16_t);
103+
DECLARE_SOA_COLUMN(P, p, uint64_t);
104+
DECLARE_SOA_COLUMN(Dedx, dedx, uint16_t);
105+
DECLARE_SOA_COLUMN(Dcaxy, dcaxy, int16_t);
106+
DECLARE_SOA_COLUMN(Dcaz, dcaz, int16_t);
107+
} // namespace testtrack
108+
109+
DECLARE_SOA_TABLE(TableTrack, "AOD", "TABLETRACK",
110+
o2::soa::Index<>,
111+
testtrack::TableColId,
112+
testtrack::Charge,
113+
testtrack::P,
114+
testtrack::Dedx,
115+
testtrack::Dcaxy,
116+
testtrack::Dcaz);
117+
} // namespace o2::aod
118+
//--------------------------------------------------------
119+
using namespace o2;
120+
using namespace o2::framework;
121+
122+
struct DiffWakeTreeProducer {
123+
124+
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
125+
Configurable<int> nBinsPt{"nBinsPt", 100, "N bins in pT histo"};
126+
Configurable<double> ptThresh{"ptThresh", 20.0, "pT threshold"};
127+
Configurable<float> centMax{"centMax", 10, "centrality"};
128+
Configurable<float> zVertCut{"zVertCut", 10.0, "z_vertex cut"};
129+
Configurable<float> etaCut{"etaCut", 0.9, "eta cut"};
130+
131+
int64_t collisionCounter = 0;
132+
133+
Produces<o2::aod::TableCols> testcol;
134+
Produces<o2::aod::TableTrack> testtrack;
135+
136+
EventPlaneHelper helperEP;
137+
138+
void init(InitContext const&)
139+
{
140+
const AxisSpec axisEta{30, -1.5, +1.5, "#eta"};
141+
const AxisSpec axispT{nBinsPt, 0, 250, "p_{T}"};
142+
143+
histos.add("etaHistogram", "etaHistogram", kTH1F, {axisEta});
144+
histos.add("pTHistogram", "pTHistogram", kTH1F, {axispT});
145+
}
146+
147+
using Bcs = aod::BCs;
148+
void process(soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::TPCMults, aod::QvectorFT0Cs>::iterator const& col,
149+
soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection> const& tracks,
150+
Bcs const&)
151+
{
152+
const float maxMomentum = 173.0; // max for px, py, pz
153+
const float minMomentum = 0.1; // min for pT
154+
155+
// Event selection corresponds to sel8FullPbPb
156+
if (!col.sel8())
157+
return;
158+
if (col.centFT0C() > centMax)
159+
return; // Centrality 0 - 10 %
160+
if (std::abs(col.posZ()) > zVertCut)
161+
return; // z position < 10 cm
162+
if (!col.selection_bit(o2::aod::evsel::kNoCollInRofStandard))
163+
return;
164+
if (!col.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard))
165+
return;
166+
167+
//------ Get Run number ---------------------
168+
auto bc = col.bc_as<Bcs>();
169+
int run = bc.runNumber();
170+
//------------ EP ---------------------------
171+
double ep2 = 0.0;
172+
double ep3 = 0.0;
173+
ep2 = helperEP.GetEventPlane(col.qvecFT0CRe(), col.qvecFT0CIm(), 2);
174+
ep3 = helperEP.GetEventPlane(col.qvecFT0CRe(), col.qvecFT0CIm(), 3);
175+
176+
//------- Only events with track above some thresh ----------
177+
178+
bool eventHighpT = false;
179+
for (const auto& track : tracks) {
180+
181+
if (!track.isGlobalTrackWoPtEta())
182+
continue;
183+
if (track.pt() < minMomentum)
184+
continue;
185+
if (std::abs(track.eta()) > etaCut)
186+
continue;
187+
if (track.pt() > ptThresh) {
188+
eventHighpT = true;
189+
break;
190+
}
191+
}
192+
if (!eventHighpT)
193+
return;
194+
//------------------------------------------------------------
195+
// Translate values to less memory consuming values
196+
int16_t substituteEp2 = static_cast<int16_t>(ep2 * 1000);
197+
int16_t substituteEp3 = static_cast<int16_t>(ep3 * 1000);
198+
199+
testcol(run,
200+
col.centFT0C(),
201+
col.multTPC(),
202+
col.trackOccupancyInTimeRange(),
203+
col.ft0cOccupancyInTimeRange(),
204+
col.posX(),
205+
col.posY(),
206+
col.posZ(),
207+
substituteEp2,
208+
substituteEp3);
209+
210+
for (const auto& track : tracks) {
211+
212+
// Track cut
213+
if (!track.isGlobalTrackWoPtEta())
214+
continue; // General track cuts, but pT and eta cut are set manually
215+
if (track.pt() < minMomentum)
216+
continue;
217+
if (std::abs(track.eta()) > etaCut)
218+
continue;
219+
220+
if (std::abs(track.px()) > maxMomentum || std::abs(track.py()) > maxMomentum || std::abs(track.pz()) > maxMomentum)
221+
continue; // to avoid overflow in Substitute_p
222+
223+
histos.fill(HIST("etaHistogram"), track.eta());
224+
histos.fill(HIST("pTHistogram"), track.pt());
225+
226+
//------------ Translate values to less memory consuming values --------------------
227+
// Px, Py, Pz
228+
uint64_t substituteP = 0;
229+
uint8_t uppermostBit = 20;
230+
uint8_t lowermostBit = 0;
231+
uint64_t bitmask20Bits = 0b11111111111111111111;
232+
233+
int64_t particlePx = (track.px() * 6000);
234+
if (particlePx < 0)
235+
substituteP |= static_cast<uint64_t>(1) << uppermostBit;
236+
if (particlePx < 0)
237+
particlePx = (-1) * particlePx;
238+
substituteP |= (particlePx & bitmask20Bits) << lowermostBit;
239+
240+
uppermostBit = 41;
241+
lowermostBit = 21;
242+
int64_t particlePy = (track.py() * 6000);
243+
if (particlePy < 0)
244+
substituteP |= static_cast<uint64_t>(1) << uppermostBit;
245+
if (particlePy < 0)
246+
particlePy = (-1) * particlePy;
247+
substituteP |= (particlePy & bitmask20Bits) << lowermostBit;
248+
249+
uppermostBit = 62;
250+
lowermostBit = 42;
251+
int64_t particlePz = (track.pz() * 6000);
252+
if (particlePz < 0)
253+
substituteP |= static_cast<uint64_t>(1) << uppermostBit;
254+
if (particlePz < 0)
255+
particlePz = (-1) * particlePz;
256+
substituteP |= (particlePz & bitmask20Bits) << lowermostBit;
257+
258+
// dEdx
259+
uint16_t substituteDEDX = static_cast<uint16_t>(track.tpcSignal() * 10);
260+
261+
// DCA
262+
int16_t substituteDCAXY = static_cast<int16_t>(track.dcaXY() * 100);
263+
int16_t substituteDCAZ = static_cast<int16_t>(track.dcaZ() * 100);
264+
265+
//--------------- Fill track table ------------------
266+
testtrack(collisionCounter,
267+
track.sign(),
268+
substituteP,
269+
substituteDEDX,
270+
substituteDCAXY,
271+
substituteDCAZ);
272+
}
273+
collisionCounter++;
274+
}
275+
};
276+
277+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
278+
{
279+
return WorkflowSpec{
280+
adaptAnalysisTask<DiffWakeTreeProducer>(cfgc)};
281+
}

0 commit comments

Comments
 (0)