Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions PWGDQ/Core/HistogramsLibrary.h
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,10 @@
//
// Contact: iarsene@cern.ch, i.c.arsene@fys.uio.no
//

#ifndef PWGDQ_CORE_HISTOGRAMSLIBRARY_H_
#define PWGDQ_CORE_HISTOGRAMSLIBRARY_H_

#include <TString.h>
#include "PWGDQ/Core/HistogramManager.h"
#include "PWGDQ/Core/VarManager.h"
Expand Down Expand Up @@ -159,6 +163,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
hm->AddHistogram(histClass, "TPCncls_Run", "Number of cluster in TPC", true, (VarManager::GetNRuns() > 0 ? VarManager::GetNRuns() : 1), 0.5, 0.5 + VarManager::GetNRuns(), VarManager::kRunId,
10, -0.5, 159.5, VarManager::kTPCncls, 10, 0., 1., VarManager::kNothing, VarManager::GetRunStr().Data());
hm->AddHistogram(histClass, "TPCnclsCR", "Number of crossed rows in TPC", false, 160, -0.5, 159.5, VarManager::kTPCnclsCR);
hm->AddHistogram(histClass, "TPCncls_TPCnclsCR", "Number of TPC cluster vs Number of crossed rows in TPC", false, 160, -0.5, 159.5, VarManager::kTPCncls, 160, -0.5, 159.5, VarManager::kTPCnclsCR);
hm->AddHistogram(histClass, "IsTPCrefit", "", false, 2, -0.5, 1.5, VarManager::kIsTPCrefit);
hm->AddHistogram(histClass, "IsGoldenChi2", "", false, 2, -0.5, 1.5, VarManager::kIsGoldenChi2);
hm->AddHistogram(histClass, "TPCchi2", "TPC chi2", false, 100, 0.0, 10.0, VarManager::kTPCchi2);
Expand Down Expand Up @@ -294,6 +299,8 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
hm->AddHistogram(histClass, "Eta_Pt", "", false, 125, -2.0, 2.0, VarManager::kEta, 100, 0.0, 20.0, VarManager::kPt);
hm->AddHistogram(histClass, "Mass_VtxZ", "", true, 30, -15.0, 15.0, VarManager::kVtxZ, 100, 0.0, 20.0, VarManager::kMass);
hm->AddHistogram(histClass, "cosThetaHE", "", false, 100, -1., 1., VarManager::kCosThetaHE);
hm->AddHistogram(histClass, "PhiV", "", false, 100, 0.0, 3.2, VarManager::kPairPhiv);
hm->AddHistogram(histClass, "Mass_Pt_PhiV", "", false, 125, 0.0, 5.0, VarManager::kMass, 100, 0.0, 20.0, VarManager::kPt, 100, 0.0, 3.2, VarManager::kPairPhiv);
if (subGroupStr.Contains("vertexing-barrel")) {
hm->AddHistogram(histClass, "Lxy", "", false, 100, 0.0, 10.0, VarManager::kVertexingLxy);
hm->AddHistogram(histClass, "Lxyz", "", false, 100, 0.0, 10.0, VarManager::kVertexingLxyz);
Expand Down Expand Up @@ -372,3 +379,5 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
hm->AddHistogram(histClass, "DeltaEta_DeltaPhiSym", "", false, 20, -2.0, 2.0, VarManager::kDeltaEta, 50, -8.0, 8.0, VarManager::kDeltaPhiSym);
}
}

#endif // PWGDQ_CORE_HISTOGRAMSLIBRARY_H_
4 changes: 3 additions & 1 deletion PWGDQ/Core/VarManager.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@ void VarManager::FillEventDerived(float* values)
// Fill event-wise derived quantities (these are all quantities which can be computed just based on the values already filled in the FillEvent() function)
//
if (fgUsedVars[kRunId]) {
values[kRunId] = (fgRunMap.size() > 0 ? fgRunMap[int(values[kRunNo])] : 0);
values[kRunId] = (fgRunMap.size() > 0 ? fgRunMap[static_cast<int>(values[kRunNo])] : 0);
}
}

Expand Down Expand Up @@ -439,6 +439,8 @@ void VarManager::SetDefaultVarNames()
fgVariableUnits[kPairEta] = "";
fgVariableNames[kPairPhi] = "#varphi";
fgVariableUnits[kPairPhi] = "rad.";
fgVariableNames[kPairPhiv] = "#varphi_{V}";
fgVariableUnits[kPairPhiv] = "rad.";
fgVariableNames[kDeltaEta] = "#Delta#eta";
fgVariableUnits[kDeltaEta] = "";
fgVariableNames[kDeltaPhi] = "#Delta#phi";
Expand Down
98 changes: 82 additions & 16 deletions PWGDQ/Core/VarManager.h
Original file line number Diff line number Diff line change
Expand Up @@ -14,20 +14,21 @@
// Class to handle analysis variables
//

#ifndef VarManager_H
#define VarManager_H

#include <TObject.h>
#include <TString.h>
#include <Math/Vector4D.h>
#include "Math/Vector3D.h"
#include "Math/GenVector/Boost.h"
#include <TRandom.h>
#ifndef PWGDQ_CORE_VARMANAGER_H_
#define PWGDQ_CORE_VARMANAGER_H_

#include <vector>
#include <map>
#include <cmath>
#include <iostream>
#include <utility>

#include <TObject.h>
#include <TString.h>
#include "TRandom.h"
#include "Math/Vector4D.h"
#include "Math/Vector3D.h"
#include "Math/GenVector/Boost.h"

#include "Framework/DataTypes.h"
#include "ReconstructionDataFormats/Track.h"
Expand Down Expand Up @@ -315,6 +316,7 @@ class VarManager : public TObject
kPairPtDau,
kPairEta,
kPairPhi,
kPairPhiv,
kDeltaEta,
kDeltaPhi,
kDeltaPhiSym,
Expand Down Expand Up @@ -671,13 +673,13 @@ void VarManager::FillTrack(T const& track, float* values)
values[kIsGlobalTrack] = track.filteringFlags() & (uint64_t(1) << 0);
values[kIsGlobalTrackSDD] = track.filteringFlags() & (uint64_t(1) << 1);

values[kIsLegFromGamma] = bool(track.filteringFlags() & (uint64_t(1) << 2));
values[kIsLegFromK0S] = bool(track.filteringFlags() & (uint64_t(1) << 3));
values[kIsLegFromLambda] = bool(track.filteringFlags() & (uint64_t(1) << 4));
values[kIsLegFromAntiLambda] = bool(track.filteringFlags() & (uint64_t(1) << 5));
values[kIsLegFromOmega] = bool(track.filteringFlags() & (uint64_t(1) << 6));
values[kIsLegFromGamma] = static_cast<bool>(track.filteringFlags() & (uint64_t(1) << 2));
values[kIsLegFromK0S] = static_cast<bool>(track.filteringFlags() & (uint64_t(1) << 3));
values[kIsLegFromLambda] = static_cast<bool>(track.filteringFlags() & (uint64_t(1) << 4));
values[kIsLegFromAntiLambda] = static_cast<bool>(track.filteringFlags() & (uint64_t(1) << 5));
values[kIsLegFromOmega] = static_cast<bool>(track.filteringFlags() & (uint64_t(1) << 6));

values[kIsProtonFromLambdaAndAntiLambda] = bool((values[kIsLegFromLambda] * track.sign() > 0) || (values[kIsLegFromAntiLambda] * (-track.sign()) > 0));
values[kIsProtonFromLambdaAndAntiLambda] = static_cast<bool>((values[kIsLegFromLambda] * track.sign() > 0) || (values[kIsLegFromAntiLambda] * (-track.sign()) > 0));
}
}

Expand Down Expand Up @@ -922,6 +924,70 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values)
values[kQuadDCAsigXY] = std::sqrt((dca1sig * dca1sig + dca2sig * dca2sig) / 2);
}
}
if (fgUsedVars[kPairPhiv]) {
// cos(phiv) = w*a /|w||a|
// with w = u x v
// and a = u x z / |u x z| , unit vector perpendicular to v12 and z-direction (magnetic field)
// u = v12 / |v12| , the unit vector of v12
// v = v1 x v2 / |v1 x v2| , unit vector perpendicular to v1 and v2

float bz = fgFitterTwoProngBarrel.getBz();

// ordering of tracks, so v1 has larger momentum
if (v1.P() < v2.P()) {
ROOT::Math::PtEtaPhiMVector v3 = v1;
v1 = v2;
v2 = v3;
}

// momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by definition.
// vector product of pep X pem
float vpx = 0, vpy = 0, vpz = 0;
if (t1.sign() * t2.sign() > 0) { // Like Sign
if (bz * t1.sign() < 0) {
vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py();
vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz();
vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px();
} else {
vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py();
vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz();
vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px();
}
} else { // Unlike Sign
if (bz * t1.sign() > 0) {
vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py();
vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz();
vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px();
} else {
vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py();
vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz();
vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px();
}
}

// unit vector of pep X pem
float vx = vpx / TMath::Sqrt(vpx * vpx + vpy * vpy + vpz * vpz);
float vy = vpy / TMath::Sqrt(vpx * vpx + vpy * vpy + vpz * vpz);
float vz = vpz / TMath::Sqrt(vpx * vpx + vpy * vpy + vpz * vpz);

float px = v12.Px();
float py = v12.Py();
float pz = v12.Pz();

// unit vector of (pep+pem)
float ux = px / TMath::Sqrt(px * px + py * py + pz * pz);
float uy = py / TMath::Sqrt(px * px + py * py + pz * pz);
float uz = pz / TMath::Sqrt(px * px + py * py + pz * pz);
float ax = uy / TMath::Sqrt(ux * ux + uy * uy);
float ay = -ux / TMath::Sqrt(ux * ux + uy * uy);

// The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
float wx = uy * vz - uz * vy;
float wy = uz * vx - ux * vz;
// by construction, (wx,wy,wz) must be a unit vector. Measure angle between (wx,wy,wz) and (ax,ay,0).
// The angle between them should be small if the pair is conversion. This function then returns values close to pi!
values[kPairPhiv] = TMath::ACos(wx * ax + wy * ay); // phiv in [0,pi] //cosPhiV = wx * ax + wy * ay;
}
}

template <int pairType, typename T1, typename T2>
Expand Down Expand Up @@ -1399,4 +1465,4 @@ void VarManager::FillDileptonHadron(T1 const& dilepton, T2 const& hadron, float*
values[kDeltaEta] = dilepton.eta() - hadron.eta();
}
}
#endif
#endif // PWGDQ_CORE_VARMANAGER_H_