Skip to content

Commit

Permalink
added switch for the invariant masses
Browse files Browse the repository at this point in the history
  • Loading branch information
@eeremenk
eeremenk committed Aug 11, 2023
1 parent 239d9b6 commit 94f2e1d
Show file tree
Hide file tree
Showing 2 changed files with 90 additions and 86 deletions.
1 change: 1 addition & 0 deletions Modules/ITS/include/ITS/ITSTrackTask.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,7 @@ class ITSTrackTask : public TaskInterface

float mPiInvMass = 0.14;
float mProtonInvMass = 0.938;
Int_t mInvMasses = 1; // switch for the V0 invariant mass computation, 1 (default) - on, 0 - off

float mVertexXYsize = 0.5;
float mVertexZsize = 15.;
Expand Down
175 changes: 89 additions & 86 deletions Modules/ITS/src/ITSTrackTask.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@ void ITSTrackTask::initialize(o2::framework::InitContext& /*ctx*/)
mDoTTree = o2::quality_control_modules::common::getFromConfig<int>(mCustomParameters, "doTTree", mDoTTree);
nBCbins = o2::quality_control_modules::common::getFromConfig<int>(mCustomParameters, "nBCbins", nBCbins);
mDoNorm = o2::quality_control_modules::common::getFromConfig<int>(mCustomParameters, "doNorm", mDoNorm);
mInvMasses = o2::quality_control_modules::common::getFromConfig<int>(mCustomParameters, "InvMasses", mInvMasses);

createAllHistos();
publishHistos();
Expand Down Expand Up @@ -182,7 +183,7 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx)
// DCAFitter2 class initialization for the v0 part
using Vec3D = ROOT::Math::SVector<double, 3>; // this is a type of the fitted vertex
o2::vertexing::DCAFitter2 ft;
ft.setBz(-5.0);
ft.setBz(5.0);
ft.setPropagateToPCA(true);
ft.setMaxR(30);
ft.setMaxDZIni(0.1);
Expand All @@ -193,15 +194,15 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx)
// prepare variables for v0
float vx = 0, vy = 0, vz = 0;
float dca[2]{ 0., 0. };
float bz = -5.0;
float bz = 5.0;

// loop on tracks per ROF
for (int iROF = 0; iROF < trackRofArr.size(); iROF++) {

vMap.clear();
vEta.clear();
vPhi.clear();

int nClusterCntTrack = 0;
int nTracks = trackRofArr[iROF].getNEntries();
int start = trackRofArr[iROF].getFirstEntry();
Expand Down Expand Up @@ -254,86 +255,88 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx)
}

// Find V0s
if (track.getSign() < 0) // choose only positive tracks
continue;
track.getImpactParams(vx, vy, vz, bz, dca);
if ((track.getNumberOfClusters() < 6) || (abs(track.getTgl()) > 1.5) || (abs(dca[0]) < 0.01)) // conditions for the track acceptance
continue;

for (int intrack = start; intrack < end; intrack++) { // goes through the tracks one more time
auto& ntrack = trackArr[intrack];
if (ntrack.getSign() > 0) // choose only negative tracks
continue;
ntrack.getImpactParams(vx, vy, vz, bz, dca);
if ((ntrack.getNumberOfClusters() < 6) || (abs(ntrack.getTgl()) > 1.5) || (abs(dca[0]) < 0.01)) // conditions for the track acceptance
continue;

int nc = 0;
try {
nc = ft.process(track, ntrack);
} catch (...) {
continue;
}
if (nc == 0)
continue;

int ibest = 0;
float bestChi2 = 1e7;
for (int i = 0; i < nc; i++) {
auto chi2 = ft.getChi2AtPCACandidate(i);
if (chi2 > bestChi2)
continue;
bestChi2 = chi2;
ibest = i;
}
// conditions for v0
auto vtx = ft.getPCACandidate(ibest);
auto x = vtx[0] + 0.02985;
auto y = vtx[1] + 0.01949;
auto r = sqrt(x * x + y * y);
if ((r < 0.5) || (r > 3.5))
continue;

const auto& t0 = ft.getTrack(0, ibest); // Positive daughter track
const auto& t1 = ft.getTrack(1, ibest); // Negative daughter track

auto r0 = t0.getXYZGlo();
auto r1 = t1.getXYZGlo();
auto dx = r0.X() - r1.X();
auto dy = r0.Y() - r1.Y();
auto dz = r0.Z() - r1.Z();
auto d = sqrt(dx * dx + dy * dy + dz * dz);
if (d > 0.02)
continue;
std::array<float, 3> p0; // Positive daughter momentum
t0.getPxPyPzGlo(p0);
std::array<float, 3> p1; // Negative daughter momentum
t1.getPxPyPzGlo(p1);
std::array<float, 3> v0p; // V0 particle momentum
v0p = { p0[0] + p1[0], p0[1] + p1[1], p0[2] + p1[2] };

// Strangness inv mass calculation
auto pV0 = sqrt(v0p[0] * v0p[0] + v0p[1] * v0p[1] + v0p[2] * v0p[2]); // Particle momentum
auto p2DaughterPos = p0[0] * p0[0] + p0[1] * p0[1] + p0[2] * p0[2]; // Positive daughter momentum
auto p2DaughterNeg = p1[0] * p1[0] + p1[1] * p1[1] + p1[2] * p1[2]; // Negative daughter momentum
// K0s
auto enDaughterPos = sqrt(mPiInvMass * mPiInvMass + p2DaughterPos); // Positive daughter energy
auto enDaughterNeg = sqrt(mPiInvMass * mPiInvMass + p2DaughterNeg); // Negative daughter energy
auto enV0 = enDaughterPos + enDaughterNeg;
auto K0sInvMass = sqrt(enV0 * enV0 - pV0 * pV0);
hInvMassK0s->Fill(K0sInvMass);
// Lambda
enDaughterPos = sqrt(mProtonInvMass * mProtonInvMass + p2DaughterPos); // Positive daughter energy
enDaughterNeg = sqrt(mPiInvMass * mPiInvMass + p2DaughterNeg); // Negative daughter energy
enV0 = enDaughterPos + enDaughterNeg;
auto LambdaInvMass = sqrt(enV0 * enV0 - pV0 * pV0);
hInvMassLambda->Fill(LambdaInvMass);
// LambdaBar
enDaughterPos = sqrt(mPiInvMass * mPiInvMass + p2DaughterPos); // Positive daughter energy
enDaughterNeg = sqrt(mProtonInvMass * mProtonInvMass + p2DaughterNeg); // Negative daughter energy
enV0 = enDaughterPos + enDaughterNeg;
auto LambdaBarInvMass = sqrt(enV0 * enV0 - pV0 * pV0);
hInvMassLambdaBar->Fill(LambdaBarInvMass);
if (mInvMasses == 1){
if (track.getSign() < 0) // choose only positive tracks
continue;
track.getImpactParams(vx, vy, vz, bz, dca);
if ((track.getNumberOfClusters() < 6) || (abs(track.getTgl()) > 1.5) || (abs(dca[0]) < 0.06)) // conditions for the track acceptance
continue;

for (int intrack = start; intrack < end; intrack++) { // goes through the tracks one more time
auto& ntrack = trackArr[intrack];
if (ntrack.getSign() > 0) // choose only negative tracks
continue;
ntrack.getImpactParams(vx, vy, vz, bz, dca);
if ((ntrack.getNumberOfClusters() < 6) || (abs(ntrack.getTgl()) > 1.5) || (abs(dca[0]) < 0.06)) // conditions for the track acceptance
continue;

int nc = 0;
try {
nc = ft.process(track, ntrack);
} catch (...) {
continue;
}
if (nc == 0)
continue;

int ibest = 0;
float bestChi2 = 1e7;
for (int i = 0; i < nc; i++) {
auto chi2 = ft.getChi2AtPCACandidate(i);
if (chi2 > bestChi2)
continue;
bestChi2 = chi2;
ibest = i;
}
// conditions for v0
auto vtx = ft.getPCACandidate(ibest);
auto x = vtx[0] + 0.02985;
auto y = vtx[1] + 0.01949;
auto r = sqrt(x * x + y * y);
if ((r < 0.5) || (r > 3.5))
continue;

const auto& t0 = ft.getTrack(0, ibest); // Positive daughter track
const auto& t1 = ft.getTrack(1, ibest); // Negative daughter track

auto r0 = t0.getXYZGlo();
auto r1 = t1.getXYZGlo();
auto dx = r0.X() - r1.X();
auto dy = r0.Y() - r1.Y();
auto dz = r0.Z() - r1.Z();
auto d = sqrt(dx * dx + dy * dy + dz * dz);
if (d > 0.02)
continue;
std::array<float, 3> p0; // Positive daughter momentum
t0.getPxPyPzGlo(p0);
std::array<float, 3> p1; // Negative daughter momentum
t1.getPxPyPzGlo(p1);
std::array<float, 3> v0p; // V0 particle momentum
v0p = { p0[0] + p1[0], p0[1] + p1[1], p0[2] + p1[2] };

// Strangness inv mass calculation
auto pV0 = sqrt(v0p[0] * v0p[0] + v0p[1] * v0p[1] + v0p[2] * v0p[2]); // Particle momentum
auto p2DaughterPos = p0[0] * p0[0] + p0[1] * p0[1] + p0[2] * p0[2]; // Positive daughter momentum
auto p2DaughterNeg = p1[0] * p1[0] + p1[1] * p1[1] + p1[2] * p1[2]; // Negative daughter momentum
// K0s
auto enDaughterPos = sqrt(mPiInvMass * mPiInvMass + p2DaughterPos); // Positive daughter energy
auto enDaughterNeg = sqrt(mPiInvMass * mPiInvMass + p2DaughterNeg); // Negative daughter energy
auto enV0 = enDaughterPos + enDaughterNeg;
auto K0sInvMass = sqrt(enV0 * enV0 - pV0 * pV0);
hInvMassK0s->Fill(K0sInvMass);
// Lambda
enDaughterPos = sqrt(mProtonInvMass * mProtonInvMass + p2DaughterPos); // Positive daughter energy
enDaughterNeg = sqrt(mPiInvMass * mPiInvMass + p2DaughterNeg); // Negative daughter energy
enV0 = enDaughterPos + enDaughterNeg;
auto LambdaInvMass = sqrt(enV0 * enV0 - pV0 * pV0);
hInvMassLambda->Fill(LambdaInvMass);
// LambdaBar
enDaughterPos = sqrt(mPiInvMass * mPiInvMass + p2DaughterPos); // Positive daughter energy
enDaughterNeg = sqrt(mProtonInvMass * mProtonInvMass + p2DaughterNeg); // Negative daughter energy
enV0 = enDaughterPos + enDaughterNeg;
auto LambdaBarInvMass = sqrt(enV0 * enV0 - pV0 * pV0);
hInvMassLambdaBar->Fill(LambdaBarInvMass);
}
}
}

Expand Down Expand Up @@ -639,19 +642,19 @@ void ITSTrackTask::createAllHistos()
}

// Invariant mass K0s, Lambda, LambdaBar
hInvMassK0s = new TH1D("hInvMassK0s", "K0s invariant mass", 40, 0.0, 1.0);
hInvMassK0s = new TH1D("hInvMassK0s", "K0s invariant mass", 80, 0.0, 1.0);
hInvMassK0s->SetTitle(Form("Invariant mass of K0s"));
addObject(hInvMassK0s);
formatAxes(hInvMassK0s, "m_{inv} (Gev/c)", "Counts", 1, 1.10);
hInvMassK0s->SetStats(0);

hInvMassLambda = new TH1D("hInvMassLambda", "Lambda invariant mass", 40, 1.0, 2.0);
hInvMassLambda = new TH1D("hInvMassLambda", "Lambda invariant mass", 80, 1.0, 2.0);
hInvMassLambda->SetTitle(Form("Invariant mass of Lambda"));
addObject(hInvMassLambda);
formatAxes(hInvMassLambda, "m_{inv} (Gev/c)", "Counts", 1, 1.10);
hInvMassLambda->SetStats(0);

hInvMassLambdaBar = new TH1D("hInvMassLambdaBar", "LambdaBar invariant mass", 40, 1.0, 2.0);
hInvMassLambdaBar = new TH1D("hInvMassLambdaBar", "LambdaBar invariant mass", 80, 1.0, 2.0);
hInvMassLambdaBar->SetTitle(Form("Invariant mass of LambdaBar"));
addObject(hInvMassLambdaBar);
formatAxes(hInvMassLambdaBar, "m_{inv} (Gev/c)", "Counts", 1, 1.10);
Expand Down

0 comments on commit 94f2e1d

Please sign in to comment.