Skip to content

Commit 0e18d19

Browse files
mherzer28mherzer28
authored
[PWGLF] mass and bethe bloch calculation (#11385)
1 parent ef4b9ab commit 0e18d19

File tree

1 file changed

+41
-20
lines changed

1 file changed

+41
-20
lines changed

‎PWGLF/TableProducer/Nuspex/trHeAnalysis.cxx‎

Lines changed: 41 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -110,13 +110,14 @@ static const std::vector<int> particlePdgCodes{
110110
static const std::vector<float> particleMasses{
111111
o2::constants::physics::MassTriton, o2::constants::physics::MassHelium3};
112112
static const std::vector<int> particleCharge{1, 2};
113+
static const std::vector<float> particleChargeFactor{2.3, 2.55};
113114
static const std::vector<std::string> betheBlochParNames{
114115
"p0", "p1", "p2", "p3", "p4", "resolution"};
115116
constexpr float betheBlochDefault[nParticles][nBetheParams]{
116-
{0.313129, 181.664226, 2779397163087.684082, 2.130773, 29.609643,
117-
0.09}, // triton
118-
{70.584685, 3.196364, 0.133878, 2.731736, 1.675617, 0.09}}; // Helion
119-
117+
{0.248753, 3.58634, 0.0167065, 2.29194, 0.774344,
118+
0.07}, // triton
119+
{0.0274556, 18.3054, 3.99987e-05, 3.17219, 11.1775,
120+
0.07}}; // Helion
120121
} // namespace
121122
using namespace o2;
122123
using namespace o2::framework;
@@ -135,16 +136,18 @@ class Particle
135136
float mass;
136137
int charge;
137138
float resolution;
139+
float chargeFactor;
138140
std::vector<float> betheParams;
139141
static constexpr int NNumBetheParams = 5;
140142

141143
Particle(const std::string name_, int pdgCode_, float mass_, int charge_,
142-
LabeledArray<float> bethe)
144+
LabeledArray<float> bethe, float chargeFactor_)
143145
{
144146
name = TString(name_);
145147
pdgCode = pdgCode_;
146148
mass = mass_;
147149
charge = charge_;
150+
chargeFactor = chargeFactor_;
148151

149152
resolution =
150153
bethe.get(name, "resolution"); // Access the "resolution" parameter
@@ -213,7 +216,7 @@ struct TrHeAnalysis {
213216
Configurable<float> cfgCutMaxTofMassH3{"cfgCutMaxTofMassH3", 3.32f, "Maximum TOF mass H3"};
214217
// Set the kinematic and PID cuts for tracks
215218
struct : ConfigurableGroup {
216-
Configurable<float> pCut{"pCut", 0.3f, "Value of the p selection for spectra (default 0.3)"};
219+
Configurable<float> pCut{"pCut", 0.6f, "Value of the p selection for spectra (default 0.3)"};
217220
Configurable<float> etaCut{"etaCut", 0.8f, "Value of the eta selection for spectra (default 0.8)"};
218221
Configurable<float> yLowCut{"yLowCut", -1.0f, "Value of the low rapidity selection for spectra (default -1.0)"};
219222
Configurable<float> yHighCut{"yHighCut", 1.0f, "Value of the high rapidity selection for spectra (default 1.0)"};
@@ -294,7 +297,7 @@ struct TrHeAnalysis {
294297
for (int i = 0; i < nParticles; i++) {
295298
particles.push_back(Particle(particleNames.at(i), particlePdgCodes.at(i),
296299
particleMasses.at(i), particleCharge.at(i),
297-
cfgBetheBlochParams));
300+
cfgBetheBlochParams, particleChargeFactor.at(i)));
298301
}
299302
}
300303
void process(soa::Join<aod::Collisions, aod::EvSels>::iterator const& event,
@@ -387,7 +390,7 @@ struct TrHeAnalysis {
387390
histos.fill(HIST("histogram/cuts"), 12);
388391
continue;
389392
}
390-
if (track.mass() < cfgCutMinTofMassH3 || track.mass() > cfgCutMaxTofMassH3) {
393+
if (getMass(track) < cfgCutMinTofMassH3 || getMass(track) > cfgCutMaxTofMassH3) {
391394
histos.fill(HIST("histogram/cuts"), 13);
392395
continue;
393396
}
@@ -402,8 +405,8 @@ struct TrHeAnalysis {
402405
float tPhi = track.phi();
403406
int8_t tCharge = track.sign();
404407
float tH3DeDx = track.tpcSignal();
405-
float tnSigmaTpc = track.tpcNSigmaTr();
406-
float tTofSignalH3 = track.mass();
408+
float tnSigmaTpc = getTPCnSigma(track, particles.at(0));
409+
float tTofSignalH3 = getMass(track);
407410
float tDcaXY = track.dcaXY();
408411
float tDcaZ = track.dcaZ();
409412
float tSigmaYX = track.sigmaY();
@@ -449,8 +452,8 @@ struct TrHeAnalysis {
449452
float tPhi = track.phi();
450453
int8_t tCharge = 2.f * track.sign();
451454
float tHeDeDx = track.tpcSignal();
452-
float tnSigmaTpc = track.tpcNSigmaHe();
453-
float tTofSignalHe = track.mass();
455+
float tnSigmaTpc = getTPCnSigma(track, particles.at(1));
456+
float tTofSignalHe = getMass(track);
454457
float tDcaXY = track.dcaXY();
455458
float tDcaZ = track.dcaZ();
456459
float tSigmaYX = track.sigmaY();
@@ -542,7 +545,7 @@ struct TrHeAnalysis {
542545
histos.fill(HIST("histogram/cuts"), 12);
543546
continue;
544547
}
545-
if (track.mass() < cfgCutMinTofMassH3 || track.mass() > cfgCutMaxTofMassH3) {
548+
if (getMass(track) < cfgCutMinTofMassH3 || getMass(track) > cfgCutMaxTofMassH3) {
546549
histos.fill(HIST("histogram/cuts"), 13);
547550
continue;
548551
}
@@ -558,7 +561,7 @@ struct TrHeAnalysis {
558561
int8_t tCharge = track.sign();
559562
float tH3DeDx = track.tpcSignal();
560563
float tnSigmaTpc = track.tpcNSigmaTr();
561-
float tTofSignalH3 = track.mass();
564+
float tTofSignalH3 = getMass(track);
562565
float tDcaXY = track.dcaXY();
563566
float tDcaZ = track.dcaZ();
564567
float tSigmaYX = track.sigmaY();
@@ -604,7 +607,7 @@ struct TrHeAnalysis {
604607
int8_t tCharge = 2.f * track.sign();
605608
float tHeDeDx = track.tpcSignal();
606609
float tnSigmaTpc = track.tpcNSigmaHe();
607-
float tTofSignalHe = track.mass();
610+
float tTofSignalHe = getMass(track);
608611
float tDcaXY = track.dcaXY();
609612
float tDcaZ = track.dcaZ();
610613
float tSigmaYX = track.sigmaY();
@@ -634,17 +637,26 @@ struct TrHeAnalysis {
634637
if (!track.hasTPC())
635638
return -999;
636639

637-
float expBethe{tpc::BetheBlochAleph(
638-
static_cast<float>(particle.charge * rigidity / particle.mass),
639-
particle.betheParams[0], particle.betheParams[1],
640-
particle.betheParams[2], particle.betheParams[3],
641-
particle.betheParams[4])};
640+
float expBethe{betheBlochAleph(particle, rigidity)};
642641
float expSigma{expBethe * particle.resolution};
643642
float sigmaTPC =
644643
static_cast<float>((track.tpcSignal() - expBethe) / expSigma);
645644
return sigmaTPC;
646645
}
647646

647+
template <class T>
648+
float betheBlochAleph(Particle const& particle, T const& rigidity)
649+
{
650+
double bg = particle.charge * rigidity / particle.mass;
651+
double beta = bg / std::sqrt(1. + bg * bg);
652+
double aa = std::pow(beta, particle.betheParams[3]);
653+
double bb = std::pow(1. / bg, particle.betheParams[4]);
654+
if ((particle.betheParams[2] + bb) <= 0)
655+
return 0;
656+
bb = std::log(particle.betheParams[2] + bb);
657+
return std::pow(particle.charge, particle.chargeFactor) * 50 * (particle.betheParams[1] - aa - bb) * particle.betheParams[0] / aa;
658+
}
659+
648660
template <class T>
649661
float getMeanItsClsSize(T const& track)
650662
{
@@ -668,6 +680,15 @@ struct TrHeAnalysis {
668680
bool hePID = track.pidForTracking() == o2::track::PID::Helium3 || track.pidForTracking() == o2::track::PID::Alpha;
669681
return hePID ? track.tpcInnerParam() / 2 : track.tpcInnerParam();
670682
}
683+
template <class T>
684+
float getMass(T const& track)
685+
{
686+
const float beta = track.beta();
687+
const float rigidity = getRigidity(track);
688+
float gamma = 1 / std::sqrt(1 - beta * beta);
689+
float mass = (rigidity / std::sqrt(gamma * gamma - 1));
690+
return mass;
691+
}
671692
};
672693

673694
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)