Skip to content

Commit bafc417

Browse files
ddobrigkchiarazampolli
ddobrigk
authored and
chiarazampolli
committed
Add Flow Mapper for synthetic v2 addition (#13288)
* Add Flow Mapper for synthetic v2 addition This PR, related also to #13287, adds a Flow Mapping utility that works via rotating a set of phi angles with respect to an event plane in a controlled manner such as to provide a realistic azimuthal anisotropy (to order n = 2) in generators such as PYTHIA. It works by afterburning any particle sample after hadronization, preserving the number of particles but laying them out differently in azimuth. The implementation itself is based around the calculation of a phiOld -> phiNew look-up table that has as dimensions phiOld, impact parameter and transverse momentum. The LUT is calculated using the inversion of the cumulative function that is desired at each (impact parameter, pT) bin and intermediate values of any of the three relevant variables in the LUT are always calculated using linear interpolation. Special care is taken to ensure that v2 -> 0 as b -> 0 and pT -> 0. The LUT calculation takes as input a v2 vs pT curve (typically from measured data) and eccentricity vs b curve (typically from glauber MC). The v2 vs pT is provided in a base centrality (typically 40-50%) and the eccentricity values are used as scaling variable for other centralities. More details here: https://www.dropbox.com/scl/fi/ss57me45ykp2f6960a2gt/DDChinellato-SynthV2Gen-01.pdf?rlkey=pph4d1286giy3sd65rz136010&dl=0 (cherry picked from commit 1511172)
1 parent d10870c commit bafc417

File tree

4 files changed

+234
-0
lines changed

4 files changed

+234
-0
lines changed

‎Generators/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ o2_add_library(Generators
4040
src/GeneratorTParticle.cxx
4141
src/GeneratorTParticleParam.cxx
4242
src/GeneratorService.cxx
43+
src/FlowMapper.cxx
4344
$<$<BOOL:${pythia_FOUND}>:src/GeneratorPythia8.cxx>
4445
$<$<BOOL:${pythia_FOUND}>:src/DecayerPythia8.cxx>
4546
$<$<BOOL:${pythia_FOUND}>:src/GeneratorPythia8Param.cxx>
@@ -82,6 +83,7 @@ set(headers
8283
include/Generators/GeneratorTParticleParam.h
8384
include/Generators/GeneratorService.h
8485
include/Generators/BoxGenerator.h
86+
include/Generators/FlowMapper.h
8587
)
8688

8789
if(pythia_FOUND)

‎Generators/include/Generators/FlowMapper.h‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
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+
#ifndef ALICEO2_EVENTGEN_FLOWMAPPER_H_
13+
#define ALICEO2_EVENTGEN_FLOWMAPPER_H_
14+
15+
#include "TH1D.h"
16+
#include "TH3D.h"
17+
#include "TF1.h"
18+
19+
namespace o2
20+
{
21+
namespace eventgen
22+
{
23+
24+
/*****************************************************************/
25+
/*****************************************************************/
26+
27+
// this class implements a mapper that introduces a synthetic v2
28+
// into an otherwise uniform initial distribution of phis. It can
29+
// be used, for instance, to create artificial flow of realistic
30+
// intensity in PYTHIA simulations.
31+
32+
// The input histograms are to be read from the CCDB and can
33+
// be customized if necessary. Multiple copies of this mapper
34+
// could be used in case different species should have different
35+
// additional flow.
36+
37+
// N.B.: the main advantages of this mapper is that:
38+
// 1) it preserves total number of particles
39+
// 2) it retains a (distorted) event structure from
40+
// an original event generator (e.g. PYTHIA)
41+
42+
class FlowMapper
43+
{
44+
public:
45+
// Constructor
46+
FlowMapper();
47+
48+
void Setv2VsPt(TH1D hv2VsPtProvided);
49+
void SetEccVsB(TH1D hEccVsBProvided);
50+
51+
void SetNBinsPhi(long mBinsPhiProvided) { mBinsPhi = mBinsPhiProvided; };
52+
void SetPrecision(long mPrecisionProvided) { mPrecision = mPrecisionProvided; };
53+
void SetDerivative(long mDerivativeProvided) { mDerivative = mDerivativeProvided; };
54+
55+
void CreateLUT(); // to be called if all is set
56+
57+
Double_t MapPhi(Double_t lPhiInput, Double_t b, Double_t pt);
58+
59+
long mBinsPhi; // number of phi bins to use
60+
double mPrecision = 1e-6; // precision threshold for numerical inversion success
61+
double mDerivative = 1e-4; // delta-X for derivative calculation
62+
63+
std::unique_ptr<TH1D> mhv2vsPt; // input v2 vs pT from measurement
64+
std::unique_ptr<TH1D> mhEccVsB; // ecc vs B (from Glauber MC or elsewhere)
65+
66+
// Cumulative function to be inverted
67+
std::unique_ptr<TF1> mCumulative;
68+
69+
// the look-up table
70+
std::unique_ptr<TH3D> mhLUT;
71+
72+
ClassDef(FlowMapper, 1);
73+
};
74+
75+
/*****************************************************************/
76+
/*****************************************************************/
77+
78+
} // namespace eventgen
79+
} // namespace o2
80+
81+
#endif /* ALICEO2_EVENTGEN_FLOWMAPPER_H_ */

‎Generators/src/FlowMapper.cxx‎

Lines changed: 150 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,150 @@
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+
#include "Generators/FlowMapper.h"
13+
#include "TH1D.h"
14+
#include "TH3D.h"
15+
#include "TF1.h"
16+
#include <fairlogger/Logger.h>
17+
18+
namespace o2
19+
{
20+
namespace eventgen
21+
{
22+
23+
/*****************************************************************/
24+
/*****************************************************************/
25+
26+
FlowMapper::FlowMapper()
27+
{
28+
// base constructor. Creates cumulative function only so that's already in place but not much else.
29+
mCumulative = std::make_unique<TF1>("mCumulative", "x+[0]*TMath::Sin(2*x)", 0, 2 * TMath::Pi());
30+
mBinsPhi = 4000; // a first guess
31+
}
32+
33+
void FlowMapper::Setv2VsPt(TH1D hv2VsPtProvided)
34+
{
35+
// Sets the v2 vs pT to be used.
36+
mhv2vsPt = std::make_unique<TH1D>(hv2VsPtProvided);
37+
}
38+
void FlowMapper::SetEccVsB(TH1D hEccVsBProvided)
39+
{
40+
// Sets the v2 vs pT to be used.
41+
mhEccVsB = std::make_unique<TH1D>(hEccVsBProvided);
42+
}
43+
void FlowMapper::CreateLUT()
44+
{
45+
if (!mhv2vsPt) {
46+
LOG(fatal) << "You did not specify a v2 vs pT histogram!";
47+
return;
48+
}
49+
if (!mhEccVsB) {
50+
LOG(fatal) << "You did not specify an ecc vs B histogram!";
51+
return;
52+
}
53+
LOG(info) << "Proceeding to creating a look-up table...";
54+
const Long_t nbinsB = mhEccVsB->GetNbinsX();
55+
const Long_t nbinsPt = mhv2vsPt->GetNbinsX();
56+
const Long_t nbinsPhi = mBinsPhi; // constant in this context necessary
57+
std::vector<double> binsB(nbinsB + 1, 0);
58+
std::vector<double> binsPt(nbinsPt + 1, 0);
59+
std::vector<double> binsPhi(nbinsPhi + 1, 0);
60+
61+
for (int ii = 0; ii < nbinsB + 1; ii++) {
62+
binsB[ii] = mhEccVsB->GetBinLowEdge(ii + 1);
63+
}
64+
for (int ii = 0; ii < nbinsPt + 1; ii++) {
65+
binsPt[ii] = mhv2vsPt->GetBinLowEdge(ii + 1);
66+
}
67+
for (int ii = 0; ii < nbinsPhi + 1; ii++) {
68+
binsPhi[ii] = static_cast<Double_t>(ii) * 2 * TMath::Pi() / static_cast<Double_t>(nbinsPhi);
69+
}
70+
71+
// std::make_unique<TH1F>("hSign", "Sign of electric charge;charge sign", 3, -1.5, 1.5);
72+
73+
mhLUT = std::make_unique<TH3D>("mhLUT", "", nbinsB, binsB.data(), nbinsPt, binsPt.data(), nbinsPhi, binsPhi.data());
74+
75+
// loop over each centrality (b) bin
76+
for (int ic = 0; ic < nbinsB; ic++) {
77+
// loop over each pt bin
78+
for (int ip = 0; ip < nbinsPt; ip++) {
79+
// find target v2 value and set cumulative for inversion
80+
double v2target = mhv2vsPt->GetBinContent(ip + 1) * mhEccVsB->GetBinContent(ic + 1);
81+
LOG(info) << "At b ~ " << mhEccVsB->GetBinCenter(ic + 1) << ", pt ~ " << mhv2vsPt->GetBinCenter(ip + 1) << ", ref v2 is " << mhv2vsPt->GetBinContent(ip + 1) << ", scale is " << mhEccVsB->GetBinContent(ic + 1) << ", target v2 is " << v2target << ", inverting...";
82+
mCumulative->SetParameter(0, v2target); // set up
83+
for (Int_t ia = 0; ia < nbinsPhi; ia++) {
84+
// Look systematically for the X value that gives this Y
85+
// There are probably better ways of doing this, but OK
86+
Double_t lY = mhLUT->GetZaxis()->GetBinCenter(ia + 1);
87+
Double_t lX = lY; // a first reasonable guess
88+
Bool_t lConverged = kFALSE;
89+
while (!lConverged) {
90+
Double_t lDistance = mCumulative->Eval(lX) - lY;
91+
if (TMath::Abs(lDistance) < mPrecision) {
92+
lConverged = kTRUE;
93+
break;
94+
}
95+
Double_t lDerivativeValue = mDerivative / (mCumulative->Eval(lX + mDerivative) - mCumulative->Eval(lX));
96+
lX = lX - lDistance * lDerivativeValue * 0.25; // 0.5: speed factor, don't overshoot but control reasonable
97+
}
98+
mhLUT->SetBinContent(ic + 1, ip + 1, ia + 1, lX);
99+
}
100+
}
101+
}
102+
}
103+
104+
Double_t FlowMapper::MapPhi(Double_t lPhiInput, Double_t b, Double_t pt)
105+
{
106+
Int_t lLowestPeriod = TMath::Floor(lPhiInput / (2 * TMath::Pi()));
107+
Double_t lPhiOld = lPhiInput - 2 * lLowestPeriod * TMath::Pi();
108+
Double_t lPhiNew = lPhiOld;
109+
110+
// Avoid interpolation problems in dimension: pT
111+
Double_t lMaxPt = mhLUT->GetYaxis()->GetBinCenter(mhLUT->GetYaxis()->GetNbins());
112+
Double_t lMinPt = mhLUT->GetYaxis()->GetBinCenter(1);
113+
if (pt > lMaxPt) {
114+
pt = lMaxPt; // avoid interpolation problems at edge
115+
}
116+
117+
Double_t phiWidth = mhLUT->GetZaxis()->GetBinWidth(1); // any bin, assume constant
118+
119+
// Valid if not at edges. If at edges, that's ok, do not map
120+
bool validPhi = lPhiNew > phiWidth / 2.0f && lPhiNew < 2.0 * TMath::Pi() - phiWidth / 2.0f;
121+
122+
// If at very high b, do not map
123+
bool validB = b < mhLUT->GetXaxis()->GetBinCenter(mhLUT->GetXaxis()->GetNbins());
124+
Double_t minB = mhLUT->GetXaxis()->GetBinCenter(1);
125+
126+
if (validPhi && validB) {
127+
128+
Double_t scaleFactor = 1.0; // no need if not special conditions
129+
if (pt < lMinPt) {
130+
scaleFactor *= pt / lMinPt; // downscale the difference, zero at zero pT
131+
pt = lMinPt;
132+
}
133+
if (b < minB) {
134+
scaleFactor *= b / minB; // downscale the difference, zero at zero b
135+
b = minB;
136+
}
137+
lPhiNew = mhLUT->Interpolate(b, pt, lPhiOld);
138+
139+
lPhiNew = scaleFactor * lPhiNew + (1.0 - scaleFactor) * lPhiOld;
140+
}
141+
return lPhiNew + 2.0 * lLowestPeriod * TMath::Pi();
142+
}
143+
144+
/*****************************************************************/
145+
/*****************************************************************/
146+
147+
} /* namespace eventgen */
148+
} /* namespace o2 */
149+
150+
ClassImp(o2::eventgen::FlowMapper);

‎Generators/src/GeneratorsLinkDef.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -74,5 +74,6 @@
7474
#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::GeneratorFileOrCmdParam> + ;
7575

7676
#pragma link C++ class o2::eventgen::BoxGenerator + ;
77+
#pragma link C++ class o2::eventgen::FlowMapper + ;
7778

7879
#endif

0 commit comments

Comments
 (0)