Skip to content

Commit

Permalink
Fixes for InteractionTag, TrackStudy, add different t,z statistics
Browse files Browse the repository at this point in the history
  • Loading branch information
@shahor02
shahor02 committed Nov 4, 2023
1 parent 84f7458 commit b1ffd4e
Show file tree
Hide file tree
Showing 2 changed files with 93 additions and 68 deletions.
2 changes: 1 addition & 1 deletion Detectors/FIT/FT0/reconstruction/include/FT0Reconstruction/InteractionTag.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ struct InteractionTag : public o2::conf::ConfigurableParamHelper<InteractionTag>
int minAmplitudeC = 1;
bool isSelected(const RecPoints& rp) const
{
return rp.getTrigger().getVertex() && rp.getTrigger().getAmplA() >= minAmplitudeA && rp.getTrigger().getAmplC() >= minAmplitudeC && (rp.getTrigger().getAmplA() + rp.getTrigger().getAmplC()) > minAmplitudeAC;
return rp.getTrigger().getVertex() && rp.getTrigger().getAmplA() >= minAmplitudeA && rp.getTrigger().getAmplC() >= minAmplitudeC && (rp.getTrigger().getAmplA() + rp.getTrigger().getAmplC()) >= minAmplitudeAC;
}

O2ParamDef(InteractionTag, "ft0tag");
Expand Down
159 changes: 92 additions & 67 deletions Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#include "DataFormatsGlobalTracking/RecoContainerCreateTracksVariadic.h"
#include "ReconstructionDataFormats/TrackTPCITS.h"
#include "ReconstructionDataFormats/GlobalTrackID.h"
#include "DataFormatsCalibration/MeanVertexObject.h"
#include "DetectorsBase/Propagator.h"
#include "DetectorsBase/GeometryManager.h"
#include "SimulationDataFormat/MCEventLabel.h"
Expand Down Expand Up @@ -77,6 +78,7 @@ class TrackingStudySpec : public Task
int mMaxNeighbours = 3;
float mMaxVTTimeDiff = 80.; // \mus
GTrackID::mask_t mTracksSrc{};
o2::dataformats::MeanVertexObject mMeanVtx{};
o2::steer::MCKinematicsReader mcReader; // reader of MC information
};

Expand Down Expand Up @@ -112,6 +114,7 @@ void TrackingStudySpec::updateTimeDependentParams(ProcessingContext& pc)
} else {
mITSROFrameLengthMUS = alpParams.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS * 1e-3; // ITS ROFrame duration in \mus
}
pc.inputs().get<o2::dataformats::MeanVertexObject*>("meanvtx");
}
}

Expand All @@ -123,58 +126,62 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
auto prop = o2::base::Propagator::Instance();
auto FITInfo = recoData.getFT0RecPoints();
static int TFCount = 0;
int nv = vtxRefs.size() - 1;
int nv = vtxRefs.size();
o2::dataformats::PrimaryVertexExt pveDummy;
std::vector<o2::dataformats::PrimaryVertexExt> pveVec;
float tBiasITS = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameBiasInBC * o2::constants::lhc::LHCBunchSpacingMUS;
o2::dataformats::PrimaryVertex vtxDummy(mMeanVtx.getPos(), {}, {}, 0);
std::vector<o2::dataformats::PrimaryVertexExt> pveVec(nv - 1);
const auto& alpParams = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance();
float tBiasITS = alpParams.roFrameBiasInBC * o2::constants::lhc::LHCBunchSpacingMUS;
std::vector<float> dtvec, dzvec;
const o2::ft0::InteractionTag& ft0Params = o2::ft0::InteractionTag::Instance();

for (int iv = 0; iv < nv; iv++) {
LOGP(debug, "processing PV {} of {}", iv, nv);
const auto& vtref = vtxRefs[iv];
auto pv = pvvec[iv];
auto& pve = pveVec.emplace_back();
static_cast<o2::dataformats::PrimaryVertex&>(pve) = pv;
dtvec.clear();
dzvec.clear();
dtvec.reserve(pv.getNContributors());
dzvec.reserve(pv.getNContributors());

float bestTimeDiff = 1000, bestTime = -999;
int bestFTID = -1;
if (mTracksSrc[GTrackID::FT0]) {
for (int ift0 = vtref.getFirstEntryOfSource(GTrackID::FT0); ift0 < vtref.getFirstEntryOfSource(GTrackID::FT0) + vtref.getEntriesOfSource(GTrackID::FT0); ift0++) {
const auto& ft0 = FITInfo[trackIndex[ift0]];
if (ft0.isValidTime(o2::ft0::RecPoints::TimeMean) && ft0.getTrigger().getAmplA() + ft0.getTrigger().getAmplC() > 20) {
auto fitTime = ft0.getInteractionRecord().differenceInBCMUS(recoData.startIR);
if (std::abs(fitTime - pv.getTimeStamp().getTimeStamp()) < bestTimeDiff) {
bestTimeDiff = fitTime - pv.getTimeStamp().getTimeStamp();
bestFTID = trackIndex[ift0];
if (iv != nv - 1) {
auto& pve = pveVec[iv];
static_cast<o2::dataformats::PrimaryVertex&>(pve) = pvvec[iv];
dtvec.clear();
dzvec.clear();
dtvec.reserve(pve.getNContributors());
dzvec.reserve(pve.getNContributors());
float bestTimeDiff = 1000, bestTime = -999;
int bestFTID = -1;
if (mTracksSrc[GTrackID::FT0]) {
for (int ift0 = vtref.getFirstEntryOfSource(GTrackID::FT0); ift0 < vtref.getFirstEntryOfSource(GTrackID::FT0) + vtref.getEntriesOfSource(GTrackID::FT0); ift0++) {
const auto& ft0 = FITInfo[trackIndex[ift0]];
if (ft0Params.isSelected(ft0)) {
auto fitTime = ft0.getInteractionRecord().differenceInBCMUS(recoData.startIR);
if (std::abs(fitTime - pve.getTimeStamp().getTimeStamp()) < bestTimeDiff) {
bestTimeDiff = fitTime - pve.getTimeStamp().getTimeStamp();
bestFTID = trackIndex[ift0];
}
}
}
} else {
LOGP(warn, "FT0 is not requested, cannot set complete vertex info");
}
} else {
LOGP(warn, "FT0 is not requested, cannot set complete vertex info");
}
if (bestFTID >= 0) {
pve.FT0A = FITInfo[bestFTID].getTrigger().getAmplA();
pve.FT0C = FITInfo[bestFTID].getTrigger().getAmplC();
pve.FT0Time = FITInfo[bestFTID].getInteractionRecord().differenceInBCMUS(recoData.startIR);
if (bestFTID >= 0) {
pve.FT0A = FITInfo[bestFTID].getTrigger().getAmplA();
pve.FT0C = FITInfo[bestFTID].getTrigger().getAmplC();
pve.FT0Time = FITInfo[bestFTID].getInteractionRecord().differenceInBCMUS(recoData.startIR);
}
pve.VtxID = iv;
}
pve.VtxID = iv;
float meanT = 0, meanZ = 0, rmsT = 0, rmsZ = 0;
float meanTW = 0, meanZW = 0, rmsTW = 0, rmsZW = 0, WT = 0, WZ = 0;
float meanT0 = 0, rmsT0 = 0;
float meanTW0 = 0, rmsTW0 = 0, WT0 = 0;
int nContAdd = 0, nContAdd0 = 0, ntITS = 0;
int nAdjusted = 0;
for (int is = 0; is < GTrackID::NSources; is++) {
bool skipTracks = (!GTrackID::getSourceDetectorsMask(is)[GTrackID::ITS] || !mTracksSrc[is] || !recoData.isTrackSourceLoaded(is));
int idMin = vtref.getFirstEntryOfSource(is), idMax = idMin + vtref.getEntriesOfSource(is);
for (int i = idMin; i < idMax; i++) {
auto vid = trackIndex[i];
bool pvCont = vid.isPVContributor();
if (pvCont) {
pve.nSrc[is]++;
pveVec[iv].nSrc[is]++;
}
if (skipTracks) {
continue;
Expand All @@ -183,7 +190,7 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
auto trc = recoData.getTrackParam(vid);
float xmin = trc.getX();
o2::dataformats::DCA dca;
if (!prop->propagateToDCA(pv, trc, prop->getNominalBz(), 2., o2::base::PropagatorF::MatCorrType::USEMatCorrLUT, &dca)) {
if (!prop->propagateToDCA(iv == nv - 1 ? vtxDummy : pvvec[iv], trc, prop->getNominalBz(), 2., o2::base::PropagatorF::MatCorrType::USEMatCorrLUT, &dca)) {
continue;
}
float ttime = 0, ttimeE = 0;
Expand All @@ -196,11 +203,20 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
if (++ntITS > 0) { // do not allow ITS in the MAD
acceptForPV0 = false;
}
} else if (vid.getSource() == GTrackID::ITSTPC) {
float bcf = ttime / o2::constants::lhc::LHCBunchSpacingMUS + o2::constants::lhc::LHCMaxBunches;
int bcWrtROF = int(bcf - alpParams.roFrameBiasInBC) % alpParams.roFrameLengthInBC;
if (bcWrtROF == 0) {
float dbc = bcf - (int(bcf / alpParams.roFrameBiasInBC)) * alpParams.roFrameBiasInBC;
if (std::abs(dbc) < 1e-6 && (++nAdjusted) > 1) {
acceptForPV0 = false; // do not allow more than 1 adjusted track MAD
}
}
} else if (vid.getSource() == GTrackID::TPC) {
ttimeE *= o2::constants::lhc::LHCBunchSpacingMUS * 8;
}
if (pvCont) {
float dt = ttime - pve.getTimeStamp().getTimeStamp();
float dt = ttime - pveVec[iv].getTimeStamp().getTimeStamp();
float tW = 1. / (ttimeE * ttimeE), zW = 1. / trc.getSigmaZ2();
dzvec.push_back(dca.getZ());
WT += tW;
Expand Down Expand Up @@ -228,45 +244,48 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
}
(*mDBGOut) << "dca"
<< "tfID=" << TFCount << "ttime=" << ttime << "ttimeE=" << ttimeE
<< "gid=" << vid << "pv=" << pv << "trc=" << trc << "pvCont=" << pvCont << "ambig=" << ambig << "dca=" << dca << "xmin=" << xmin << "\n";
<< "gid=" << vid << "pv=" << (iv == nv - 1 ? vtxDummy : pvvec[iv]) << "trc=" << trc << "pvCont=" << pvCont << "ambig=" << ambig << "dca=" << dca << "xmin=" << xmin << "\n";
}
}

if (nContAdd) {
rmsT /= nContAdd;
rmsZ /= nContAdd;
meanT /= nContAdd;
meanZ /= nContAdd;
pve.rmsT = (rmsT - meanT * meanT);
pve.rmsT = pve.rmsT > 0 ? std::sqrt(pve.rmsT) : 0;
pve.rmsZ = rmsZ - meanZ * meanZ;
pve.rmsZ = pve.rmsZ > 0 ? std::sqrt(pve.rmsZ) : 0;
}
if (nContAdd0) {
rmsT0 /= nContAdd0;
meanT0 /= nContAdd0;
pve.rmsT0 = (rmsT0 - meanT0 * meanT0);
pve.rmsT0 = pve.rmsT0 > 0 ? std::sqrt(pve.rmsT0) : 0;
}
if (WT0 > 0) {
rmsTW0 /= WT0;
meanTW0 /= WT0;
pve.rmsTW0 = (rmsTW0 - meanTW0 * meanTW0);
pve.rmsTW0 = pve.rmsTW0 > 0 ? std::sqrt(pve.rmsTW0) : 0;
}
//
if (WT > 0 && WZ > 0) {
rmsTW /= WT;
meanTW /= WT;
pve.rmsTW = (rmsTW - meanTW * meanTW);
pve.rmsTW = pve.rmsTW > 0 ? std::sqrt(pve.rmsTW) : 0;
rmsZW /= WZ;
meanZW /= WZ;
pve.rmsZW = rmsZW - meanZW * meanZW;
pve.rmsZW = pve.rmsZ > 0 ? std::sqrt(pve.rmsZ) : 0;
if (iv != nv - 1) {
auto& pve = pveVec[iv];
if (nContAdd) {
rmsT /= nContAdd;
rmsZ /= nContAdd;
meanT /= nContAdd;
meanZ /= nContAdd;
pve.rmsT = (rmsT - meanT * meanT);
pve.rmsT = pve.rmsT > 0 ? std::sqrt(pve.rmsT) : 0;
pve.rmsZ = rmsZ - meanZ * meanZ;
pve.rmsZ = pve.rmsZ > 0 ? std::sqrt(pve.rmsZ) : 0;
}
if (nContAdd0) {
rmsT0 /= nContAdd0;
meanT0 /= nContAdd0;
pve.rmsT0 = (rmsT0 - meanT0 * meanT0);
pve.rmsT0 = pve.rmsT0 > 0 ? std::sqrt(pve.rmsT0) : 0;
}
if (WT0 > 0) {
rmsTW0 /= WT0;
meanTW0 /= WT0;
pve.rmsTW0 = (rmsTW0 - meanTW0 * meanTW0);
pve.rmsTW0 = pve.rmsTW0 > 0 ? std::sqrt(pve.rmsTW0) : 0;
}
//
if (WT > 0 && WZ > 0) {
rmsTW /= WT;
meanTW /= WT;
pve.rmsTW = (rmsTW - meanTW * meanTW);
pve.rmsTW = pve.rmsTW > 0 ? std::sqrt(pve.rmsTW) : 0;
rmsZW /= WZ;
meanZW /= WZ;
pve.rmsZW = rmsZW - meanZW * meanZW;
pve.rmsZW = pve.rmsZ > 0 ? std::sqrt(pve.rmsZ) : 0;
}
pve.tMAD = o2::math_utils::MAD2Sigma(dtvec.size(), dtvec.data());
pve.zMAD = o2::math_utils::MAD2Sigma(dzvec.size(), dzvec.data());
}
pve.tMAD = o2::math_utils::MAD2Sigma(dtvec.size(), dtvec.data());
pve.zMAD = o2::math_utils::MAD2Sigma(dzvec.size(), dzvec.data());
}
int nvtot = mMaxNeighbours < 0 ? -1 : (int)pveVec.size();

Expand Down Expand Up @@ -368,6 +387,11 @@ void TrackingStudySpec::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) {
return;
}
if (matcher == ConcreteDataMatcher("GLO", "MEANVERTEX", 0)) {
LOG(info) << "Imposing new MeanVertex: " << ((const o2::dataformats::MeanVertexObject*)obj)->asString();
mMeanVtx = *(const o2::dataformats::MeanVertexObject*)obj;
return;
}
}

DataProcessorSpec getTrackingStudySpec(GTrackID::mask_t srcTracks, GTrackID::mask_t srcClusters, bool useMC)
Expand All @@ -378,6 +402,7 @@ DataProcessorSpec getTrackingStudySpec(GTrackID::mask_t srcTracks, GTrackID::mas
dataRequest->requestTracks(srcTracks, useMC);
dataRequest->requestClusters(srcClusters, useMC);
dataRequest->requestPrimaryVertertices(useMC);
dataRequest->inputs.emplace_back("meanvtx", "GLO", "MEANVERTEX", 0, Lifetime::Condition, ccdbParamSpec("GLO/Calib/MeanVertex", {}, 1));
auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(false, // orbitResetTime
true, // GRPECS=true
true, // GRPLHCIF
Expand Down

0 comments on commit b1ffd4e

Please sign in to comment.