Project
Loading...
Searching...
No Matches
AODProducerWorkflowSpec.cxx
Go to the documentation of this file.
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
13
41#include "MathUtils/Utils.h"
51#include "Framework/DataTypes.h"
55#include "FT0Base/Geometry.h"
80#include "ZDCBase/Constants.h"
83#include "TOFBase/Utils.h"
84#include "O2Version.h"
85#include "TMath.h"
86#include "MathUtils/Utils.h"
87#include "Math/SMatrix.h"
88#include "TString.h"
89#include <limits>
90#include <map>
91#include <numeric>
92#include <type_traits>
93#include <unordered_map>
94#include <set>
95#include <string>
96#include <vector>
97#include <thread>
98#include "TLorentzVector.h"
99#include "TVector3.h"
100#include "MathUtils/Tsallis.h"
101#include <random>
102#ifdef WITH_OPENMP
103#include <omp.h>
104#endif
105#include <filesystem>
106#include <nlohmann/json.hpp>
107
108using namespace o2::framework;
109using namespace o2::math_utils::detail;
116
117namespace o2::aodproducer
118{
119
120void AODProducerWorkflowDPL::createCTPReadout(const o2::globaltracking::RecoContainer& recoData, std::vector<o2::ctp::CTPDigit>& ctpDigits, ProcessingContext& pc)
121{
122 // Extraxt CTP Config from CCDB
123 const auto ctpcfg = pc.inputs().get<o2::ctp::CTPConfiguration*>("ctpconfig");
124 ctpcfg->printStream(std::cout);
125 // o2::ctp::CTPConfiguration ctpcfg = o2::ctp::CTPRunManager::getConfigFromCCDB(-1, std::to_string(runNumber)); // how to get run
126 // Extract inputs from recoData
127 uint64_t classMaskEMCAL = 0, classMaskTRD = 0, classMaskPHOSCPV = 0;
128 for (const auto& trgclass : ctpcfg->getCTPClasses()) {
129 if (trgclass.cluster->getClusterDetNames().find("EMC") != std::string::npos) {
130 classMaskEMCAL = trgclass.classMask;
131 }
132 if (trgclass.cluster->getClusterDetNames().find("PHS") != std::string::npos) {
133 classMaskPHOSCPV = trgclass.classMask;
134 }
135 if (trgclass.cluster->getClusterDetNames().find("TRD") != std::string::npos) {
136 classMaskTRD = trgclass.classMask;
137 }
138 }
139 LOG(info) << "createCTPReadout: Class Mask EMCAL -> " << classMaskEMCAL;
140 LOG(info) << "createCTPReadout: Class Mask PHOS/CPV -> " << classMaskPHOSCPV;
141 LOG(info) << "createCTPReadout: Class Mask TRD -> " << classMaskTRD;
142
143 // const auto& fddRecPoints = recoData.getFDDRecPoints();
144 // const auto& fv0RecPoints = recoData.getFV0RecPoints();
145 const auto& triggerrecordEMCAL = recoData.getEMCALTriggers();
146 const auto& triggerrecordPHOSCPV = recoData.getPHOSTriggers();
147 const auto& triggerrecordTRD = recoData.getTRDTriggerRecords();
148 // For EMCAL filter remove calibration triggers
149 std::vector<o2::emcal::TriggerRecord> triggerRecordEMCALPhys;
150 for (const auto& trg : triggerrecordEMCAL) {
151 if (trg.getTriggerBits() & o2::trigger::Cal) {
152 continue;
153 }
154 triggerRecordEMCALPhys.push_back(trg);
155 }
156 // const auto& triggerrecordTRD =recoData.getITSTPCTRDTriggers()
157 //
158
159 // Find TVX triggers, only TRD/EMCAL/PHOS/CPV triggers in coincidence will be accepted
160 std::set<uint64_t> bcsMapT0triggers;
161 const auto& ft0RecPoints = recoData.getFT0RecPoints();
162 for (auto& ft0RecPoint : ft0RecPoints) {
163 auto t0triggers = ft0RecPoint.getTrigger();
164 if (t0triggers.getVertex()) {
165 uint64_t globalBC = ft0RecPoint.getInteractionRecord().toLong();
166 bcsMapT0triggers.insert(globalBC);
167 }
168 }
169
170 auto genericCTPDigitizer = [&bcsMapT0triggers, &ctpDigits](auto triggerrecords, uint64_t classmask) -> int {
171 // Strategy:
172 // find detector trigger based on trigger record from readout and add CTPDigit if trigger there
173 int cntwarnings = 0;
174 uint32_t orbitPrev = 0;
175 uint16_t bcPrev = 0;
176 for (auto& trigger : triggerrecords) {
177 auto orbitPrevT = orbitPrev;
178 auto bcPrevT = bcPrev;
179 bcPrev = trigger.getBCData().bc;
180 orbitPrev = trigger.getBCData().orbit;
181 // dedicated for TRD: remove bogus triggers
182 if (orbitPrev < orbitPrevT || bcPrev >= o2::constants::lhc::LHCMaxBunches || (orbitPrev == orbitPrevT && bcPrev < bcPrevT)) {
183 cntwarnings++;
184 // LOGP(warning, "Bogus TRD trigger at bc:{}/orbit:{} (previous was {}/{}), with {} tracklets and {} digits",bcPrev, orbitPrev, bcPrevT, orbitPrevT, trig.getNumberOfTracklets(), trig.getNumberOfDigits());
185 } else {
186 uint64_t globalBC = trigger.getBCData().toLong();
187 auto t0entry = bcsMapT0triggers.find(globalBC);
188 if (t0entry != bcsMapT0triggers.end()) {
189 auto ctpdig = std::find_if(ctpDigits.begin(), ctpDigits.end(), [globalBC](const o2::ctp::CTPDigit& dig) { return static_cast<uint64_t>(dig.intRecord.toLong()) == globalBC; });
190 if (ctpdig != ctpDigits.end()) {
191 // CTP digit existing from other trigger, merge detector class mask
192 ctpdig->CTPClassMask |= std::bitset<64>(classmask);
193 LOG(debug) << "createCTPReadout: Merging " << classmask << " CTP digits with existing digit, CTP mask " << ctpdig->CTPClassMask;
194 } else {
195 // New CTP digit needed
196 LOG(debug) << "createCTPReadout: New CTP digit needed for class " << classmask << std::endl;
197 auto& ctpdigNew = ctpDigits.emplace_back();
198 ctpdigNew.intRecord.setFromLong(globalBC);
199 ctpdigNew.CTPClassMask = classmask;
200 }
201 } else {
202 LOG(warning) << "createCTPReadout: Found " << classmask << " and no MTVX:" << globalBC;
203 }
204 }
205 }
206 return cntwarnings;
207 };
208
209 auto warningsTRD = genericCTPDigitizer(triggerrecordTRD, classMaskTRD);
210 auto warningsEMCAL = genericCTPDigitizer(triggerRecordEMCALPhys, classMaskEMCAL);
211 auto warningsPHOSCPV = genericCTPDigitizer(triggerrecordPHOSCPV, classMaskPHOSCPV);
212
213 LOG(info) << "createCTPReadout:# of TRD bogus triggers:" << warningsTRD;
214 LOG(info) << "createCTPReadout:# of EMCAL bogus triggers:" << warningsEMCAL;
215 LOG(info) << "createCTPReadout:# of PHOS/CPV bogus triggers:" << warningsPHOSCPV;
216}
217
218void AODProducerWorkflowDPL::collectBCs(const o2::globaltracking::RecoContainer& data,
219 const std::vector<o2::InteractionTimeRecord>& mcRecords,
220 std::map<uint64_t, int>& bcsMap)
221{
222 const auto& primVertices = data.getPrimaryVertices();
223 const auto& fddRecPoints = data.getFDDRecPoints();
224 const auto& ft0RecPoints = data.getFT0RecPoints();
225 const auto& fv0RecPoints = data.getFV0RecPoints();
226 const auto& caloEMCCellsTRGR = data.getEMCALTriggers();
227 const auto& caloPHOSCellsTRGR = data.getPHOSTriggers();
228 const auto& cpvTRGR = data.getCPVTriggers();
229 const auto& ctpDigits = data.getCTPDigits();
230 const auto& zdcBCRecData = data.getZDCBCRecData();
231
232 bcsMap[mStartIR.toLong()] = 1; // store the start of TF
233
234 // collecting non-empty BCs and enumerating them
235 for (auto& rec : mcRecords) {
236 uint64_t globalBC = rec.toLong();
237 bcsMap[globalBC] = 1;
238 }
239
240 for (auto& fddRecPoint : fddRecPoints) {
241 uint64_t globalBC = fddRecPoint.getInteractionRecord().toLong();
242 bcsMap[globalBC] = 1;
243 }
244
245 for (auto& ft0RecPoint : ft0RecPoints) {
246 uint64_t globalBC = ft0RecPoint.getInteractionRecord().toLong();
247 bcsMap[globalBC] = 1;
248 }
249
250 for (auto& fv0RecPoint : fv0RecPoints) {
251 uint64_t globalBC = fv0RecPoint.getInteractionRecord().toLong();
252 bcsMap[globalBC] = 1;
253 }
254
255 for (auto& zdcRecData : zdcBCRecData) {
256 uint64_t globalBC = zdcRecData.ir.toLong();
257 bcsMap[globalBC] = 1;
258 }
259
260 for (auto& vertex : primVertices) {
261 auto& timeStamp = vertex.getTimeStamp();
262 double tsTimeStamp = timeStamp.getTimeStamp() * 1E3; // mus to ns
263 uint64_t globalBC = relativeTime_to_GlobalBC(tsTimeStamp);
264 bcsMap[globalBC] = 1;
265 }
266
267 for (auto& emcaltrg : caloEMCCellsTRGR) {
268 uint64_t globalBC = emcaltrg.getBCData().toLong();
269 bcsMap[globalBC] = 1;
270 }
271
272 for (auto& phostrg : caloPHOSCellsTRGR) {
273 uint64_t globalBC = phostrg.getBCData().toLong();
274 bcsMap[globalBC] = 1;
275 }
276
277 for (auto& cpvtrg : cpvTRGR) {
278 uint64_t globalBC = cpvtrg.getBCData().toLong();
279 bcsMap[globalBC] = 1;
280 }
281
282 for (auto& ctpDigit : ctpDigits) {
283 uint64_t globalBC = ctpDigit.intRecord.toLong();
284 bcsMap[globalBC] = 1;
285 }
286
287 int bcID = 0;
288 for (auto& item : bcsMap) {
289 item.second = bcID;
290 bcID++;
291 }
292}
293
294template <typename TracksCursorType, typename TracksCovCursorType>
295void AODProducerWorkflowDPL::addToTracksTable(TracksCursorType& tracksCursor, TracksCovCursorType& tracksCovCursor,
296 const o2::track::TrackParCov& track, int collisionID, aod::track::TrackTypeEnum type)
297{
298 tracksCursor(collisionID,
299 type,
300 truncateFloatFraction(track.getX(), mTrackX),
301 truncateFloatFraction(track.getAlpha(), mTrackAlpha),
302 track.getY(),
303 track.getZ(),
304 truncateFloatFraction(track.getSnp(), mTrackSnp),
305 truncateFloatFraction(track.getTgl(), mTrackTgl),
306 truncateFloatFraction(track.getQ2Pt(), mTrack1Pt));
307 // trackscov
308 float sY = TMath::Sqrt(track.getSigmaY2()), sZ = TMath::Sqrt(track.getSigmaZ2()), sSnp = TMath::Sqrt(track.getSigmaSnp2()),
309 sTgl = TMath::Sqrt(track.getSigmaTgl2()), sQ2Pt = TMath::Sqrt(track.getSigma1Pt2());
310 tracksCovCursor(truncateFloatFraction(sY, mTrackCovDiag),
311 truncateFloatFraction(sZ, mTrackCovDiag),
312 truncateFloatFraction(sSnp, mTrackCovDiag),
313 truncateFloatFraction(sTgl, mTrackCovDiag),
314 truncateFloatFraction(sQ2Pt, mTrackCovDiag),
315 (Char_t)(128. * track.getSigmaZY() / (sZ * sY)),
316 (Char_t)(128. * track.getSigmaSnpY() / (sSnp * sY)),
317 (Char_t)(128. * track.getSigmaSnpZ() / (sSnp * sZ)),
318 (Char_t)(128. * track.getSigmaTglY() / (sTgl * sY)),
319 (Char_t)(128. * track.getSigmaTglZ() / (sTgl * sZ)),
320 (Char_t)(128. * track.getSigmaTglSnp() / (sTgl * sSnp)),
321 (Char_t)(128. * track.getSigma1PtY() / (sQ2Pt * sY)),
322 (Char_t)(128. * track.getSigma1PtZ() / (sQ2Pt * sZ)),
323 (Char_t)(128. * track.getSigma1PtSnp() / (sQ2Pt * sSnp)),
324 (Char_t)(128. * track.getSigma1PtTgl() / (sQ2Pt * sTgl)));
325}
326
327template <typename TracksExtraCursorType>
328void AODProducerWorkflowDPL::addToTracksExtraTable(TracksExtraCursorType& tracksExtraCursor, TrackExtraInfo& extraInfoHolder)
329{
330 // In case of TPC-only tracks, do not truncate the time error since we encapsulate there a special encoding of
331 // the deltaFwd/Bwd times
332 auto trackTimeRes = extraInfoHolder.trackTimeRes;
333 if (!extraInfoHolder.isTPConly) {
334 trackTimeRes = truncateFloatFraction(trackTimeRes, mTrackTimeError);
335 }
336
337 // extra
338 tracksExtraCursor(truncateFloatFraction(extraInfoHolder.tpcInnerParam, mTrack1Pt),
339 extraInfoHolder.flags,
340 extraInfoHolder.itsClusterSizes,
341 extraInfoHolder.tpcNClsFindable,
342 extraInfoHolder.tpcNClsFindableMinusFound,
343 extraInfoHolder.tpcNClsFindableMinusPID,
344 extraInfoHolder.tpcNClsFindableMinusCrossedRows,
345 extraInfoHolder.tpcNClsShared,
346 extraInfoHolder.trdPattern,
347 truncateFloatFraction(extraInfoHolder.itsChi2NCl, mTrackChi2),
348 truncateFloatFraction(extraInfoHolder.tpcChi2NCl, mTrackChi2),
349 truncateFloatFraction(extraInfoHolder.trdChi2, mTrackChi2),
350 truncateFloatFraction(extraInfoHolder.tofChi2, mTrackChi2),
351 truncateFloatFraction(extraInfoHolder.tpcSignal, mTrackSignal),
352 truncateFloatFraction(extraInfoHolder.trdSignal, mTrackSignal),
353 truncateFloatFraction(extraInfoHolder.length, mTrackSignal),
354 truncateFloatFraction(extraInfoHolder.tofExpMom, mTrack1Pt),
355 truncateFloatFraction(extraInfoHolder.trackEtaEMCAL, mTrackPosEMCAL),
356 truncateFloatFraction(extraInfoHolder.trackPhiEMCAL, mTrackPosEMCAL),
357 truncateFloatFraction(extraInfoHolder.trackTime, mTrackTime),
358 trackTimeRes);
359}
360
361template <typename TracksQACursorType>
362void AODProducerWorkflowDPL::addToTracksQATable(TracksQACursorType& tracksQACursor, TrackQA& trackQAInfoHolder)
363{
364 tracksQACursor(
365 trackQAInfoHolder.trackID,
366 truncateFloatFraction(trackQAInfoHolder.tpcTime0, mTPCTime0),
367 truncateFloatFraction(trackQAInfoHolder.tpcdEdxNorm, mTrackSignal),
368 trackQAInfoHolder.tpcdcaR,
369 trackQAInfoHolder.tpcdcaZ,
370 trackQAInfoHolder.tpcClusterByteMask,
371 trackQAInfoHolder.tpcdEdxMax0R,
372 trackQAInfoHolder.tpcdEdxMax1R,
373 trackQAInfoHolder.tpcdEdxMax2R,
374 trackQAInfoHolder.tpcdEdxMax3R,
375 trackQAInfoHolder.tpcdEdxTot0R,
376 trackQAInfoHolder.tpcdEdxTot1R,
377 trackQAInfoHolder.tpcdEdxTot2R,
378 trackQAInfoHolder.tpcdEdxTot3R,
379 trackQAInfoHolder.dRefContY,
380 trackQAInfoHolder.dRefContZ,
381 trackQAInfoHolder.dRefContSnp,
382 trackQAInfoHolder.dRefContTgl,
383 trackQAInfoHolder.dRefContQ2Pt,
384 trackQAInfoHolder.dRefGloY,
385 trackQAInfoHolder.dRefGloZ,
386 trackQAInfoHolder.dRefGloSnp,
387 trackQAInfoHolder.dRefGloTgl,
388 trackQAInfoHolder.dRefGloQ2Pt,
389 trackQAInfoHolder.dTofdX,
390 trackQAInfoHolder.dTofdZ);
391}
392
393template <typename mftTracksCursorType, typename AmbigMFTTracksCursorType>
394void AODProducerWorkflowDPL::addToMFTTracksTable(mftTracksCursorType& mftTracksCursor, AmbigMFTTracksCursorType& ambigMFTTracksCursor,
395 GIndex trackID, const o2::globaltracking::RecoContainer& data, int collisionID,
396 std::uint64_t collisionBC, const std::map<uint64_t, int>& bcsMap)
397{
398 // mft tracks
399 int bcSlice[2] = {-1, -1};
400 const auto& track = data.getMFTTrack(trackID);
401 const auto& rof = data.getMFTTracksROFRecords()[mMFTROFs[trackID.getIndex()]];
402 float trackTime = rof.getBCData().differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingNS + mMFTROFrameHalfLengthNS + mMFTROFBiasNS;
403 float trackTimeRes = mMFTROFrameHalfLengthNS;
404 bool needBCSlice = collisionID < 0;
405 std::uint64_t bcOfTimeRef;
406 if (needBCSlice) {
407 double error = mTimeMarginTrackTime + trackTimeRes;
408 bcOfTimeRef = fillBCSlice(bcSlice, trackTime - error, trackTime + error, bcsMap);
409 } else {
410 bcOfTimeRef = collisionBC - mStartIR.toLong(); // by default (unambiguous) track time is wrt collision BC
411 }
412 trackTime -= bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS;
413
414 // the Cellular Automaton track-finding algorithm flag is stored in first of the 4 bits not used for the cluster size
415 uint64_t mftClusterSizesAndTrackFlags = track.getClusterSizes();
416 mftClusterSizesAndTrackFlags |= (track.isCA()) ? (1ULL << (60)) : 0;
417
418 mftTracksCursor(collisionID,
419 track.getX(),
420 track.getY(),
421 truncateFloatFraction(track.getZ(), mTrackX), // for the forward tracks Z has the same role as X in barrel
422 truncateFloatFraction(track.getPhi(), mTrackAlpha),
423 truncateFloatFraction(track.getTanl(), mTrackTgl),
424 truncateFloatFraction(track.getInvQPt(), mTrack1Pt),
425 mftClusterSizesAndTrackFlags,
426 truncateFloatFraction(track.getTrackChi2(), mTrackChi2),
427 truncateFloatFraction(trackTime, mTrackTime),
428 truncateFloatFraction(trackTimeRes, mTrackTimeError));
429 if (needBCSlice) {
430 ambigMFTTracksCursor(mTableTrMFTID, bcSlice);
431 }
432}
433
434template <typename TracksCursorType, typename TracksCovCursorType, typename TracksExtraCursorType, typename TracksQACursorType, typename AmbigTracksCursorType,
435 typename MFTTracksCursorType, typename MFTTracksCovCursorType, typename AmbigMFTTracksCursorType,
436 typename FwdTracksCursorType, typename FwdTracksCovCursorType, typename AmbigFwdTracksCursorType, typename FwdTrkClsCursorType>
437void AODProducerWorkflowDPL::fillTrackTablesPerCollision(int collisionID,
438 std::uint64_t collisionBC,
439 const o2::dataformats::VtxTrackRef& trackRef,
440 const gsl::span<const GIndex>& GIndices,
442 TracksCursorType& tracksCursor,
443 TracksCovCursorType& tracksCovCursor,
444 TracksExtraCursorType& tracksExtraCursor,
445 TracksQACursorType& tracksQACursor,
446 AmbigTracksCursorType& ambigTracksCursor,
447 MFTTracksCursorType& mftTracksCursor,
448 MFTTracksCovCursorType& mftTracksCovCursor,
449 AmbigMFTTracksCursorType& ambigMFTTracksCursor,
450 FwdTracksCursorType& fwdTracksCursor,
451 FwdTracksCovCursorType& fwdTracksCovCursor,
452 AmbigFwdTracksCursorType& ambigFwdTracksCursor,
453 FwdTrkClsCursorType& fwdTrkClsCursor,
454 const std::map<uint64_t, int>& bcsMap)
455{
456 for (int src = GIndex::NSources; src--;) {
457 if (!GIndex::isTrackSource(src)) {
458 continue;
459 }
460 int start = trackRef.getFirstEntryOfSource(src);
461 int end = start + trackRef.getEntriesOfSource(src);
462 int nToReserve = end - start; // + last index for a given table
463 if (src == GIndex::Source::MFT) {
464 mftTracksCursor.reserve(nToReserve + mftTracksCursor.lastIndex());
466 fwdTracksCursor.reserve(nToReserve + fwdTracksCursor.lastIndex());
467 fwdTracksCovCursor.reserve(nToReserve + fwdTracksCovCursor.lastIndex());
469 mftTracksCovCursor.reserve(nToReserve + mftTracksCovCursor.lastIndex());
470 }
471 } else {
472 tracksCursor.reserve(nToReserve + tracksCursor.lastIndex());
473 tracksCovCursor.reserve(nToReserve + tracksCovCursor.lastIndex());
474 tracksExtraCursor.reserve(nToReserve + tracksExtraCursor.lastIndex());
475 }
476 for (int ti = start; ti < end; ti++) {
477 const auto& trackIndex = GIndices[ti];
478 if (GIndex::includesSource(src, mInputSources)) {
479 if (src == GIndex::Source::MFT) { // MFT tracks are treated separately since they are stored in a different table
480 if (trackIndex.isAmbiguous() && mGIDToTableMFTID.find(trackIndex) != mGIDToTableMFTID.end()) { // was it already stored ?
481 continue;
482 }
483 addToMFTTracksTable(mftTracksCursor, ambigMFTTracksCursor, trackIndex, data, collisionID, collisionBC, bcsMap);
484 mGIDToTableMFTID.emplace(trackIndex, mTableTrMFTID);
485 mTableTrMFTID++;
486 } else if (src == GIndex::Source::MCH || src == GIndex::Source::MFTMCH || src == GIndex::Source::MCHMID) { // FwdTracks tracks are treated separately since they are stored in a different table
487 if (trackIndex.isAmbiguous() && mGIDToTableFwdID.find(trackIndex) != mGIDToTableFwdID.end()) { // was it already stored ?
488 continue;
489 }
490 addToFwdTracksTable(fwdTracksCursor, fwdTracksCovCursor, ambigFwdTracksCursor, mftTracksCovCursor, trackIndex, data, collisionID, collisionBC, bcsMap);
491 mGIDToTableFwdID.emplace(trackIndex, mTableTrFwdID);
492 addClustersToFwdTrkClsTable(data, fwdTrkClsCursor, trackIndex, mTableTrFwdID);
493 mTableTrFwdID++;
494 } else {
495 // barrel track: normal tracks table
496 if (trackIndex.isAmbiguous() && mGIDToTableID.find(trackIndex) != mGIDToTableID.end()) { // was it already stored ?
497 continue;
498 }
499
500 float weight = 0;
501 static std::uniform_real_distribution<> distr(0., 1.);
502 bool writeQAData = o2::math_utils::Tsallis::downsampleTsallisCharged(data.getTrackParam(trackIndex).getPt(), mTrackQCFraction, mSqrtS, weight, distr(mGenerator));
503 auto extraInfoHolder = processBarrelTrack(collisionID, collisionBC, trackIndex, data, bcsMap);
504
505 if (writeQAData) {
506 auto trackQAInfoHolder = processBarrelTrackQA(collisionID, collisionBC, trackIndex, data, bcsMap);
507 if (std::bitset<8>(trackQAInfoHolder.tpcClusterByteMask).count() >= mTrackQCNTrCut) {
508 trackQAInfoHolder.trackID = mTableTrID;
509 // LOGP(info, "orig time0 in bc: {} diffBCRef: {}, ttime: {} -> {}", trackQAInfoHolder.tpcTime0*8, extraInfoHolder.diffBCRef, extraInfoHolder.trackTime, (trackQAInfoHolder.tpcTime0 * 8 - extraInfoHolder.diffBCRef) * o2::constants::lhc::LHCBunchSpacingNS - extraInfoHolder.trackTime);
510 trackQAInfoHolder.tpcTime0 = (trackQAInfoHolder.tpcTime0 * 8 - extraInfoHolder.diffBCRef) * o2::constants::lhc::LHCBunchSpacingNS - extraInfoHolder.trackTime;
511 // difference between TPC track time0 and stored track nominal time in ns instead of TF start
512 addToTracksQATable(tracksQACursor, trackQAInfoHolder);
513 } else {
514 writeQAData = false;
515 }
516 }
517
518 // include specific selection of tpc standalone tracks if thinning is active
519 if (mThinTracks && extraInfoHolder.isTPConly && !writeQAData) { // if trackQA is written then no check has to be done
520 auto trk = data.getTPCTrack(trackIndex);
521 if (trk.getNClusters() >= mTrackQCNCls && trk.getPt() >= mTrackQCPt) {
522 o2::dataformats::DCA dcaInfo{999.f, 999.f, 999.f, 999.f, 999.f};
523 o2::dataformats::VertexBase v = mVtx.getMeanVertex(collisionID < 0 ? 0.f : data.getPrimaryVertex(collisionID).getZ());
524 if (o2::base::Propagator::Instance()->propagateToDCABxByBz(v, trk, 2., mMatCorr, &dcaInfo) && std::abs(dcaInfo.getY()) < mTrackQCDCAxy) {
525 writeQAData = true; // just setting this to not thin the track
526 }
527 }
528 }
529
530 // Skip thinning if not enabled or track is not tpc standalone or assoc. to a V0 or qa'ed
531 if (mThinTracks && src == GIndex::Source::TPC && mGIDUsedBySVtx.find(trackIndex) == mGIDUsedBySVtx.end() && mGIDUsedByStr.find(trackIndex) == mGIDUsedByStr.end() && !writeQAData) {
532 mGIDToTableID.emplace(trackIndex, -1); // pretend skipped tracks are stored; this is safe since they are are not written to disk and -1 indicates to all users to not use this track
533 continue;
534 }
535
536 if (!extraInfoHolder.isTPConly && extraInfoHolder.trackTimeRes < 0.f) { // failed or rejected?
537 LOG(warning) << "Barrel track " << trackIndex << " has no time set, rejection is not expected : time=" << extraInfoHolder.trackTime
538 << " timeErr=" << extraInfoHolder.trackTimeRes << " BCSlice: " << extraInfoHolder.bcSlice[0] << ":" << extraInfoHolder.bcSlice[1];
539 continue;
540 }
541 const auto& trOrig = data.getTrackParam(trackIndex);
542 bool isProp = false;
543 if (mPropTracks && trOrig.getX() < mMinPropR &&
544 mGIDUsedBySVtx.find(trackIndex) == mGIDUsedBySVtx.end() &&
545 mGIDUsedByStr.find(trackIndex) == mGIDUsedByStr.end()) { // Do not propagate track assoc. to V0s and str. tracking
546 auto trackPar(trOrig);
547 isProp = propagateTrackToPV(trackPar, data, collisionID);
548 if (isProp) {
549 addToTracksTable(tracksCursor, tracksCovCursor, trackPar, collisionID, aod::track::Track);
550 }
551 }
552 if (!isProp) {
553 addToTracksTable(tracksCursor, tracksCovCursor, trOrig, collisionID, aod::track::TrackIU);
554 }
555 addToTracksExtraTable(tracksExtraCursor, extraInfoHolder);
556
557 // collecting table indices of barrel tracks for V0s table
558 if (extraInfoHolder.bcSlice[0] >= 0 && collisionID < 0) {
559 ambigTracksCursor(mTableTrID, extraInfoHolder.bcSlice);
560 }
561 mGIDToTableID.emplace(trackIndex, mTableTrID);
562 mTableTrID++;
563 }
564 }
565 }
566 }
567 if (collisionID < 0) {
568 return;
569 }
571 auto sTracks = data.getStrangeTracks();
572 tracksCursor.reserve(mVertexStrLUT[collisionID + 1] + tracksCursor.lastIndex());
573 tracksCovCursor.reserve(mVertexStrLUT[collisionID + 1] + tracksCovCursor.lastIndex());
574 tracksExtraCursor.reserve(mVertexStrLUT[collisionID + 1] + tracksExtraCursor.lastIndex());
575 for (int iS{mVertexStrLUT[collisionID]}; iS < mVertexStrLUT[collisionID + 1]; ++iS) {
576 auto& collStrTrk = mCollisionStrTrk[iS];
577 auto& sTrk = sTracks[collStrTrk.second];
578 TrackExtraInfo extraInfo;
579 extraInfo.itsChi2NCl = sTrk.mTopoChi2; // TODO: this is the total chi2 of adding the ITS clusters, the topology chi2 meaning might change in the future
580 extraInfo.itsClusterSizes = sTrk.getClusterSizes();
581 addToTracksTable(tracksCursor, tracksCovCursor, sTrk.mMother, collisionID, aod::track::StrangeTrack);
582 addToTracksExtraTable(tracksExtraCursor, extraInfo);
583 mStrTrkIndices[collStrTrk.second] = mTableTrID;
584 mTableTrID++;
585 }
586}
587
588void AODProducerWorkflowDPL::fillIndexTablesPerCollision(const o2::dataformats::VtxTrackRef& trackRef, const gsl::span<const GIndex>& GIndices, const o2::globaltracking::RecoContainer& data)
589{
590 const auto& mchmidMatches = data.getMCHMIDMatches();
591
593 int start = trackRef.getFirstEntryOfSource(src);
594 int end = start + trackRef.getEntriesOfSource(src);
595 for (int ti = start; ti < end; ti++) {
596 auto& trackIndex = GIndices[ti];
597 if (GIndex::includesSource(src, mInputSources)) {
598 if (src == GIndex::Source::MFT) {
599 if (trackIndex.isAmbiguous() && mGIDToTableMFTID.find(trackIndex) != mGIDToTableMFTID.end()) {
600 continue;
601 }
602 mGIDToTableMFTID.emplace(trackIndex, mIndexMFTID);
603 mIndexTableMFT[trackIndex.getIndex()] = mIndexMFTID;
604 mIndexMFTID++;
606 if (trackIndex.isAmbiguous() && mGIDToTableFwdID.find(trackIndex) != mGIDToTableFwdID.end()) {
607 continue;
608 }
609 mGIDToTableFwdID.emplace(trackIndex, mIndexFwdID);
610 if (src == GIndex::Source::MCH) {
611 mIndexTableFwd[trackIndex.getIndex()] = mIndexFwdID;
612 } else if (src == GIndex::Source::MCHMID) {
613 const auto& mchmidMatch = mchmidMatches[trackIndex.getIndex()];
614 const auto mchTrackID = mchmidMatch.getMCHRef().getIndex();
615 mIndexTableFwd[mchTrackID] = mIndexFwdID;
616 }
617 mIndexFwdID++;
618 }
619 }
620 }
621 }
622}
623
624template <typename FwdTracksCursorType, typename FwdTracksCovCursorType, typename AmbigFwdTracksCursorType, typename mftTracksCovCursorType>
625void AODProducerWorkflowDPL::addToFwdTracksTable(FwdTracksCursorType& fwdTracksCursor, FwdTracksCovCursorType& fwdTracksCovCursor,
626 AmbigFwdTracksCursorType& ambigFwdTracksCursor, mftTracksCovCursorType& mftTracksCovCursor, GIndex trackID,
627 const o2::globaltracking::RecoContainer& data, int collisionID, std::uint64_t collisionBC,
628 const std::map<uint64_t, int>& bcsMap)
629{
630 const auto& mchTracks = data.getMCHTracks();
631 const auto& midTracks = data.getMIDTracks();
632 const auto& mchmidMatches = data.getMCHMIDMatches();
633 const auto& mchClusters = data.getMCHTrackClusters();
634
635 FwdTrackInfo fwdInfo;
636 FwdTrackCovInfo fwdCovInfo;
637 int bcSlice[2] = {-1, -1};
638
639 // helper lambda for mch bitmap -- common for global and standalone tracks
640 auto getMCHBitMap = [&](int mchTrackID) {
641 if (mchTrackID != -1) { // check matching just in case
642 const auto& mchTrack = mchTracks[mchTrackID];
643 int first = mchTrack.getFirstClusterIdx();
644 int last = mchTrack.getLastClusterIdx();
645 for (int i = first; i <= last; i++) { // check chamberIds of all clusters
646 const auto& cluster = mchClusters[i];
647 int chamberId = cluster.getChamberId();
648 fwdInfo.mchBitMap |= 1 << chamberId;
649 }
650 }
651 };
652
653 auto getMIDBitMapBoards = [&](int midTrackID) {
654 if (midTrackID != -1) { // check matching just in case
655 const auto& midTrack = midTracks[midTrackID];
656 fwdInfo.midBitMap = midTrack.getHitMap();
657 fwdInfo.midBoards = midTrack.getEfficiencyWord();
658 }
659 };
660
661 auto extrapMCHTrack = [&](int mchTrackID) {
662 const auto& track = mchTracks[mchTrackID];
663
664 // mch standalone tracks extrapolated to vertex
665 // compute 3 sets of tracks parameters :
666 // - at vertex
667 // - at DCA
668 // - at the end of the absorber
669 // extrapolate to vertex
670 float vx = 0, vy = 0, vz = 0;
671 if (collisionID >= 0) {
672 const auto& v = data.getPrimaryVertex(collisionID);
673 vx = v.getX();
674 vy = v.getY();
675 vz = v.getZ();
676 }
677
678 o2::mch::TrackParam trackParamAtVertex(track.getZ(), track.getParameters(), track.getCovariances());
679 if (mPropMuons) {
680 double errVtx{0.0}; // FIXME: get errors associated with vertex if available
681 double errVty{0.0};
682 if (!o2::mch::TrackExtrap::extrapToVertex(trackParamAtVertex, vx, vy, vz, errVtx, errVty)) {
683 return false;
684 }
685 }
686
687 // extrapolate to DCA
688 o2::mch::TrackParam trackParamAtDCA(track.getZ(), track.getParameters());
689 if (!o2::mch::TrackExtrap::extrapToVertexWithoutBranson(trackParamAtDCA, vz)) {
690 return false;
691 }
692
693 // extrapolate to the end of the absorber
694 o2::mch::TrackParam trackParamAtRAbs(track.getZ(), track.getParameters());
695 if (!o2::mch::TrackExtrap::extrapToZ(trackParamAtRAbs, -505.)) { // FIXME: replace hardcoded 505
696 return false;
697 }
698
699 double dcaX = trackParamAtDCA.getNonBendingCoor() - vx;
700 double dcaY = trackParamAtDCA.getBendingCoor() - vy;
701 double dca = std::sqrt(dcaX * dcaX + dcaY * dcaY);
702
703 double xAbs = trackParamAtRAbs.getNonBendingCoor();
704 double yAbs = trackParamAtRAbs.getBendingCoor();
705
706 double dpdca = track.getP() * dca;
707 double dchi2 = track.getChi2OverNDF();
708
709 auto fwdmuon = mMatching.MCHtoFwd(trackParamAtVertex);
710
711 fwdInfo.x = fwdmuon.getX();
712 fwdInfo.y = fwdmuon.getY();
713 fwdInfo.z = fwdmuon.getZ();
714 fwdInfo.phi = fwdmuon.getPhi();
715 fwdInfo.tanl = fwdmuon.getTgl();
716 fwdInfo.invqpt = fwdmuon.getInvQPt();
717 fwdInfo.rabs = std::sqrt(xAbs * xAbs + yAbs * yAbs);
718 fwdInfo.chi2 = dchi2;
719 fwdInfo.pdca = dpdca;
720 fwdInfo.nClusters = track.getNClusters();
721
722 fwdCovInfo.sigX = TMath::Sqrt(fwdmuon.getCovariances()(0, 0));
723 fwdCovInfo.sigY = TMath::Sqrt(fwdmuon.getCovariances()(1, 1));
724 fwdCovInfo.sigPhi = TMath::Sqrt(fwdmuon.getCovariances()(2, 2));
725 fwdCovInfo.sigTgl = TMath::Sqrt(fwdmuon.getCovariances()(3, 3));
726 fwdCovInfo.sig1Pt = TMath::Sqrt(fwdmuon.getCovariances()(4, 4));
727 fwdCovInfo.rhoXY = (Char_t)(128. * fwdmuon.getCovariances()(0, 1) / (fwdCovInfo.sigX * fwdCovInfo.sigY));
728 fwdCovInfo.rhoPhiX = (Char_t)(128. * fwdmuon.getCovariances()(0, 2) / (fwdCovInfo.sigPhi * fwdCovInfo.sigX));
729 fwdCovInfo.rhoPhiY = (Char_t)(128. * fwdmuon.getCovariances()(1, 2) / (fwdCovInfo.sigPhi * fwdCovInfo.sigY));
730 fwdCovInfo.rhoTglX = (Char_t)(128. * fwdmuon.getCovariances()(0, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigX));
731 fwdCovInfo.rhoTglY = (Char_t)(128. * fwdmuon.getCovariances()(1, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigY));
732 fwdCovInfo.rhoTglPhi = (Char_t)(128. * fwdmuon.getCovariances()(2, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigPhi));
733 fwdCovInfo.rho1PtX = (Char_t)(128. * fwdmuon.getCovariances()(0, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigX));
734 fwdCovInfo.rho1PtY = (Char_t)(128. * fwdmuon.getCovariances()(1, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigY));
735 fwdCovInfo.rho1PtPhi = (Char_t)(128. * fwdmuon.getCovariances()(2, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigPhi));
736 fwdCovInfo.rho1PtTgl = (Char_t)(128. * fwdmuon.getCovariances()(3, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigTgl));
737
738 return true;
739 };
740
741 if (trackID.getSource() == GIndex::MCH) { // This is an MCH track
742 int mchTrackID = trackID.getIndex();
743 getMCHBitMap(mchTrackID);
744 if (!extrapMCHTrack(mchTrackID)) {
745 LOGF(warn, "Unable to extrapolate MCH track with ID %d! Dummy parameters will be used", mchTrackID);
746 }
747 fwdInfo.trackTypeId = o2::aod::fwdtrack::MCHStandaloneTrack;
748 const auto& rof = data.getMCHTracksROFRecords()[mMCHROFs[mchTrackID]];
749 auto time = rof.getTimeMUS(mStartIR).first;
750 fwdInfo.trackTime = time.getTimeStamp() * 1.e3;
751 fwdInfo.trackTimeRes = time.getTimeStampError() * 1.e3;
752 } else if (trackID.getSource() == GIndex::MCHMID) { // This is an MCH-MID track
753 fwdInfo.trackTypeId = o2::aod::fwdtrack::MuonStandaloneTrack;
754 auto mchmidMatch = mchmidMatches[trackID.getIndex()];
755 auto mchTrackID = mchmidMatch.getMCHRef().getIndex();
756 if (!extrapMCHTrack(mchTrackID)) {
757 LOGF(warn, "Unable to extrapolate MCH track with ID %d! Dummy parameters will be used", mchTrackID);
758 }
759 auto midTrackID = mchmidMatch.getMIDRef().getIndex();
760 fwdInfo.chi2matchmchmid = mchmidMatch.getMatchChi2OverNDF();
761 getMCHBitMap(mchTrackID);
762 getMIDBitMapBoards(midTrackID);
763 auto time = mchmidMatch.getTimeMUS(mStartIR).first;
764 fwdInfo.trackTime = time.getTimeStamp() * 1.e3;
765 fwdInfo.trackTimeRes = time.getTimeStampError() * 1.e3;
766 } else { // This is a GlobalMuonTrack or a GlobalForwardTrack
767 const auto& track = data.getGlobalFwdTrack(trackID);
768 const auto& mftTracks = data.getMFTTracks();
769 const auto& mfttrack = mftTracks[track.getMFTTrackID()];
770 if (!extrapMCHTrack(track.getMCHTrackID())) {
771 LOGF(warn, "Unable to extrapolate MCH track with ID %d! Dummy parameters will be used", track.getMCHTrackID());
772 }
773 fwdInfo.x = track.getX();
774 fwdInfo.y = track.getY();
775 fwdInfo.z = track.getZ();
776 fwdInfo.phi = track.getPhi();
777 fwdInfo.tanl = track.getTanl();
778 fwdInfo.invqpt = track.getInvQPt();
779 fwdInfo.chi2 = track.getTrackChi2();
780 // fwdInfo.nClusters = track.getNumberOfPoints();
781 fwdInfo.chi2matchmchmid = track.getMIDMatchingChi2();
782 fwdInfo.chi2matchmchmft = track.getMFTMCHMatchingChi2();
783 fwdInfo.matchscoremchmft = track.getMFTMCHMatchingScore();
784 fwdInfo.matchmfttrackid = mIndexTableMFT[track.getMFTTrackID()];
785 fwdInfo.matchmchtrackid = mIndexTableFwd[track.getMCHTrackID()];
786 fwdInfo.trackTime = track.getTimeMUS().getTimeStamp() * 1.e3;
787 fwdInfo.trackTimeRes = track.getTimeMUS().getTimeStampError() * 1.e3;
788
789 getMCHBitMap(track.getMCHTrackID());
790 getMIDBitMapBoards(track.getMIDTrackID());
791
792 fwdCovInfo.sigX = TMath::Sqrt(track.getCovariances()(0, 0));
793 fwdCovInfo.sigY = TMath::Sqrt(track.getCovariances()(1, 1));
794 fwdCovInfo.sigPhi = TMath::Sqrt(track.getCovariances()(2, 2));
795 fwdCovInfo.sigTgl = TMath::Sqrt(track.getCovariances()(3, 3));
796 fwdCovInfo.sig1Pt = TMath::Sqrt(track.getCovariances()(4, 4));
797 fwdCovInfo.rhoXY = (Char_t)(128. * track.getCovariances()(0, 1) / (fwdCovInfo.sigX * fwdCovInfo.sigY));
798 fwdCovInfo.rhoPhiX = (Char_t)(128. * track.getCovariances()(0, 2) / (fwdCovInfo.sigPhi * fwdCovInfo.sigX));
799 fwdCovInfo.rhoPhiY = (Char_t)(128. * track.getCovariances()(1, 2) / (fwdCovInfo.sigPhi * fwdCovInfo.sigY));
800 fwdCovInfo.rhoTglX = (Char_t)(128. * track.getCovariances()(0, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigX));
801 fwdCovInfo.rhoTglY = (Char_t)(128. * track.getCovariances()(1, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigY));
802 fwdCovInfo.rhoTglPhi = (Char_t)(128. * track.getCovariances()(2, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigPhi));
803 fwdCovInfo.rho1PtX = (Char_t)(128. * track.getCovariances()(0, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigX));
804 fwdCovInfo.rho1PtY = (Char_t)(128. * track.getCovariances()(1, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigY));
805 fwdCovInfo.rho1PtPhi = (Char_t)(128. * track.getCovariances()(2, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigPhi));
806 fwdCovInfo.rho1PtTgl = (Char_t)(128. * track.getCovariances()(3, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigTgl));
807
808 fwdInfo.trackTypeId = (fwdInfo.chi2matchmchmid >= 0) ? o2::aod::fwdtrack::GlobalMuonTrack : o2::aod::fwdtrack::GlobalForwardTrack;
809
810 float sX = TMath::Sqrt(mfttrack.getSigma2X()), sY = TMath::Sqrt(mfttrack.getSigma2Y()), sPhi = TMath::Sqrt(mfttrack.getSigma2Phi()),
811 sTgl = TMath::Sqrt(mfttrack.getSigma2Tanl()), sQ2Pt = TMath::Sqrt(mfttrack.getSigma2InvQPt());
812
813 mftTracksCovCursor(fwdInfo.matchmfttrackid,
814 truncateFloatFraction(sX, mTrackCovDiag),
815 truncateFloatFraction(sY, mTrackCovDiag),
816 truncateFloatFraction(sPhi, mTrackCovDiag),
817 truncateFloatFraction(sTgl, mTrackCovDiag),
818 truncateFloatFraction(sQ2Pt, mTrackCovDiag),
819 (Char_t)(128. * mfttrack.getCovariances()(0, 1) / (sX * sY)),
820 (Char_t)(128. * mfttrack.getCovariances()(0, 2) / (sPhi * sX)),
821 (Char_t)(128. * mfttrack.getCovariances()(1, 2) / (sPhi * sY)),
822 (Char_t)(128. * mfttrack.getCovariances()(0, 3) / (sTgl * sX)),
823 (Char_t)(128. * mfttrack.getCovariances()(1, 3) / (sTgl * sY)),
824 (Char_t)(128. * mfttrack.getCovariances()(2, 3) / (sTgl * sPhi)),
825 (Char_t)(128. * mfttrack.getCovariances()(0, 4) / (sQ2Pt * sX)),
826 (Char_t)(128. * mfttrack.getCovariances()(1, 4) / (sQ2Pt * sY)),
827 (Char_t)(128. * mfttrack.getCovariances()(2, 4) / (sQ2Pt * sPhi)),
828 (Char_t)(128. * mfttrack.getCovariances()(3, 4) / (sQ2Pt * sTgl)));
829 }
830
831 std::uint64_t bcOfTimeRef;
832 bool needBCSlice = collisionID < 0;
833 if (needBCSlice) { // need to store BC slice
834 float err = mTimeMarginTrackTime + fwdInfo.trackTimeRes;
835 bcOfTimeRef = fillBCSlice(bcSlice, fwdInfo.trackTime - err, fwdInfo.trackTime + err, bcsMap);
836 } else {
837 bcOfTimeRef = collisionBC - mStartIR.toLong(); // by default track time is wrt collision BC (unless no collision assigned)
838 }
839 fwdInfo.trackTime -= bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS;
840
841 fwdTracksCursor(collisionID,
842 fwdInfo.trackTypeId,
843 fwdInfo.x,
844 fwdInfo.y,
845 truncateFloatFraction(fwdInfo.z, mTrackX), // for the forward tracks Z has the same role as X in the barrel
846 truncateFloatFraction(fwdInfo.phi, mTrackAlpha),
847 truncateFloatFraction(fwdInfo.tanl, mTrackTgl),
848 truncateFloatFraction(fwdInfo.invqpt, mTrack1Pt),
849 fwdInfo.nClusters,
850 truncateFloatFraction(fwdInfo.pdca, mTrackX),
851 truncateFloatFraction(fwdInfo.rabs, mTrackX),
852 truncateFloatFraction(fwdInfo.chi2, mTrackChi2),
853 truncateFloatFraction(fwdInfo.chi2matchmchmid, mTrackChi2),
854 truncateFloatFraction(fwdInfo.chi2matchmchmft, mTrackChi2),
855 truncateFloatFraction(fwdInfo.matchscoremchmft, mTrackChi2),
856 fwdInfo.matchmfttrackid,
857 fwdInfo.matchmchtrackid,
858 fwdInfo.mchBitMap,
859 fwdInfo.midBitMap,
860 fwdInfo.midBoards,
861 truncateFloatFraction(fwdInfo.trackTime, mTrackTime),
862 truncateFloatFraction(fwdInfo.trackTimeRes, mTrackTimeError));
863
864 fwdTracksCovCursor(truncateFloatFraction(fwdCovInfo.sigX, mTrackCovDiag),
865 truncateFloatFraction(fwdCovInfo.sigY, mTrackCovDiag),
866 truncateFloatFraction(fwdCovInfo.sigPhi, mTrackCovDiag),
867 truncateFloatFraction(fwdCovInfo.sigTgl, mTrackCovDiag),
868 truncateFloatFraction(fwdCovInfo.sig1Pt, mTrackCovDiag),
869 fwdCovInfo.rhoXY,
870 fwdCovInfo.rhoPhiX,
871 fwdCovInfo.rhoPhiY,
872 fwdCovInfo.rhoTglX,
873 fwdCovInfo.rhoTglY,
874 fwdCovInfo.rhoTglPhi,
875 fwdCovInfo.rho1PtX,
876 fwdCovInfo.rho1PtY,
877 fwdCovInfo.rho1PtPhi,
878 fwdCovInfo.rho1PtTgl);
879
880 if (needBCSlice) {
881 ambigFwdTracksCursor(mTableTrFwdID, bcSlice);
882 }
883}
884
885//------------------------------------------------------------------
886void AODProducerWorkflowDPL::updateMCHeader(MCCollisionCursor& collisionCursor,
887 XSectionCursor& xSectionCursor,
888 PdfInfoCursor& pdfInfoCursor,
889 HeavyIonCursor& heavyIonCursor,
890 const MCEventHeader& header,
891 int collisionID,
892 int bcID,
893 float time,
894 short generatorID,
895 int sourceID)
896{
901
902 auto genID = updateMCCollisions(collisionCursor,
903 bcID,
904 time,
905 header,
906 generatorID,
907 sourceID,
908 mCollisionPosition);
909 mXSectionUpdate = (updateHepMCXSection(xSectionCursor, //
910 collisionID, //
911 genID, //
912 header, //
913 mXSectionUpdate)
914 ? HepMCUpdate::always
915 : HepMCUpdate::never);
916 mPdfInfoUpdate = (updateHepMCPdfInfo(pdfInfoCursor, //
917 collisionID, //
918 genID, //
919 header, //
920 mPdfInfoUpdate)
921 ? HepMCUpdate::always
922 : HepMCUpdate::never);
923 mHeavyIonUpdate = (updateHepMCHeavyIon(heavyIonCursor, //
924 collisionID, //
925 genID, //
926 header,
927 mHeavyIonUpdate)
928 ? HepMCUpdate::always
929 : HepMCUpdate::never);
930}
931
932void dimensionMCKeepStore(std::vector<std::vector<std::unordered_map<int, int>>>& store, int Nsources, int NEvents)
933{
934 store.resize(Nsources);
935 for (int s = 0; s < Nsources; ++s) {
936 store[s].resize(NEvents);
937 }
938}
939
940void clearMCKeepStore(std::vector<std::vector<std::unordered_map<int, int>>>& store)
941{
942 for (auto s = 0U; s < store.size(); ++s) {
943 for (auto e = 0U; e < store[s].size(); ++e) {
944 store[s][e].clear();
945 }
946 }
947}
948
949// helper function to add a particle/track to the MC keep store
950void keepMCParticle(std::vector<std::vector<std::unordered_map<int, int>>>& store, int source, int event, int track, int value = 1)
951{
952 if (track < 0) {
953 LOG(warn) << "trackID is smaller than 0. Neglecting";
954 return;
955 }
956 store[source][event][track] = value;
957}
958
959void AODProducerWorkflowDPL::fillMCParticlesTable(o2::steer::MCKinematicsReader& mcReader,
960 MCParticlesCursor& mcParticlesCursor,
961 const gsl::span<const o2::dataformats::VtxTrackRef>& primVer2TRefs,
962 const gsl::span<const GIndex>& GIndices,
964 const std::vector<std::vector<int>>& mcColToEvSrc)
965{
966 int NSources = 0;
967 int NEvents = 0;
968 for (auto& p : mcColToEvSrc) {
969 NSources = std::max(p[1], NSources);
970 NEvents = std::max(p[2], NEvents);
971 }
972 NSources++; // 0 - indexed
973 NEvents++;
974 LOG(info) << " number of events " << NEvents;
975 LOG(info) << " number of sources " << NSources;
976 dimensionMCKeepStore(mToStore, NSources, NEvents);
977
978 std::vector<int> particleIDsToKeep;
979
980 auto markMCTrackForSrc = [&](std::array<GID, GID::NSources>& contributorsGID, uint8_t src) {
981 auto mcLabel = data.getTrackMCLabel(contributorsGID[src]);
982 if (!mcLabel.isValid()) {
983 return;
984 }
985 keepMCParticle(mToStore, mcLabel.getSourceID(), mcLabel.getEventID(), mcLabel.getTrackID());
986 };
987
988 // mark reconstructed MC particles to store them into the table
989 for (auto& trackRef : primVer2TRefs) {
990 for (int src = GIndex::NSources; src--;) {
991 int start = trackRef.getFirstEntryOfSource(src);
992 int end = start + trackRef.getEntriesOfSource(src);
993 for (int ti = start; ti < end; ti++) {
994 auto& trackIndex = GIndices[ti];
995 if (GIndex::includesSource(src, mInputSources)) {
996 auto mcTruth = data.getTrackMCLabel(trackIndex);
997 if (!mcTruth.isValid()) {
998 continue;
999 }
1000 keepMCParticle(mToStore, mcTruth.getSourceID(), mcTruth.getEventID(), mcTruth.getTrackID());
1001 // treating contributors of global tracks
1002 auto contributorsGID = data.getSingleDetectorRefs(trackIndex);
1003 if (contributorsGID[GIndex::Source::TPC].isIndexSet()) {
1004 markMCTrackForSrc(contributorsGID, GIndex::Source::TPC);
1005 }
1006 if (contributorsGID[GIndex::Source::ITS].isIndexSet()) {
1007 markMCTrackForSrc(contributorsGID, GIndex::Source::ITS);
1008 }
1009 if (contributorsGID[GIndex::Source::TOF].isIndexSet()) {
1010 const auto& labelsTOF = data.getTOFClustersMCLabels()->getLabels(contributorsGID[GIndex::Source::TOF]);
1011 for (auto& mcLabel : labelsTOF) {
1012 if (!mcLabel.isValid()) {
1013 continue;
1014 }
1015 keepMCParticle(mToStore, mcLabel.getSourceID(), mcLabel.getEventID(), mcLabel.getTrackID());
1016 }
1017 }
1018 }
1019 }
1020 }
1021 }
1022 // mark calorimeter signals as reconstructed particles
1023 if (mInputSources[GIndex::EMC]) {
1024 auto& mcCaloEMCCellLabels = data.getEMCALCellsMCLabels()->getTruthArray();
1025 for (auto& mcTruth : mcCaloEMCCellLabels) {
1026 if (!mcTruth.isValid()) {
1027 continue;
1028 }
1029 keepMCParticle(mToStore, mcTruth.getSourceID(), mcTruth.getEventID(), mcTruth.getTrackID());
1030 }
1031 }
1032 if (mInputSources[GIndex::PHS]) {
1033 auto& mcCaloPHOSCellLabels = data.getPHOSCellsMCLabels()->getTruthArray();
1034 for (auto& mcTruth : mcCaloPHOSCellLabels) {
1035 if (!mcTruth.isValid()) {
1036 continue;
1037 }
1038 keepMCParticle(mToStore, mcTruth.getSourceID(), mcTruth.getEventID(), mcTruth.getTrackID());
1039 }
1040 }
1041 using namespace aodmchelpers;
1042 using MCTrackNavigator = o2::mcutils::MCTrackNavigator;
1043
1044 size_t offset = 0;
1045 for (auto& colInfo : mcColToEvSrc) { // loop over "<eventID, sourceID> <-> combined MC col. ID" key pairs
1046 int event = colInfo[2];
1047 int source = colInfo[1];
1048 int mcColId = colInfo[0];
1049 std::vector<MCTrack> const& mcParticles = mcReader.getTracks(source, event);
1050 LOG(debug) << "Event=" << event << " source=" << source << " collision=" << mcColId;
1051 auto& preselect = mToStore[source][event];
1052
1053 offset = updateParticles(mcParticlesCursor,
1054 mcColId,
1055 mcParticles,
1056 preselect,
1057 offset,
1058 mRecoOnly,
1059 source == 0, // background
1060 mMcParticleW,
1061 mMcParticleMom,
1062 mMcParticlePos);
1063
1065 }
1066}
1067
1068template <typename MCTrackLabelCursorType, typename MCMFTTrackLabelCursorType, typename MCFwdTrackLabelCursorType>
1069void AODProducerWorkflowDPL::fillMCTrackLabelsTable(MCTrackLabelCursorType& mcTrackLabelCursor,
1070 MCMFTTrackLabelCursorType& mcMFTTrackLabelCursor,
1071 MCFwdTrackLabelCursorType& mcFwdTrackLabelCursor,
1072 const o2::dataformats::VtxTrackRef& trackRef,
1073 const gsl::span<const GIndex>& primVerGIs,
1075 int vertexId)
1076{
1077 // labelMask (temporary) usage:
1078 // bit 13 -- ITS/TPC with ITS label (track of AB tracklet) different from TPC
1079 // bit 14 -- isNoise() == true
1080 // bit 15 -- isFake() == true (defined by the fakeness of the top level global track, i.e. if TOF is present, fake means that the track of the TPC label does not contribute to TOF cluster)
1081 // labelID = -1 -- label is not set
1082
1083 for (int src = GIndex::NSources; src--;) {
1084 int start = trackRef.getFirstEntryOfSource(src);
1085 int end = start + trackRef.getEntriesOfSource(src);
1086 mcMFTTrackLabelCursor.reserve(end - start + mcMFTTrackLabelCursor.lastIndex());
1087 mcFwdTrackLabelCursor.reserve(end - start + mcFwdTrackLabelCursor.lastIndex());
1088 mcTrackLabelCursor.reserve(end - start + mcTrackLabelCursor.lastIndex());
1089 for (int ti = start; ti < end; ti++) {
1090 const auto trackIndex = primVerGIs[ti];
1091
1092 // check if the label was already stored (or the track was rejected for some reason in the fillTrackTablesPerCollision)
1093 auto needToStore = [trackIndex](std::unordered_map<GIndex, int>& mp) {
1094 auto entry = mp.find(trackIndex);
1095 if (entry == mp.end() || entry->second == -1) {
1096 return false;
1097 }
1098 entry->second = -1;
1099 return true;
1100 };
1101
1102 if (GIndex::includesSource(src, mInputSources)) {
1103 auto mcTruth = data.getTrackMCLabel(trackIndex);
1104 MCLabels labelHolder{};
1105 if ((src == GIndex::Source::MFT) || (src == GIndex::Source::MFTMCH) || (src == GIndex::Source::MCH) || (src == GIndex::Source::MCHMID)) { // treating mft and fwd labels separately
1106 if (!needToStore(src == GIndex::Source::MFT ? mGIDToTableMFTID : mGIDToTableFwdID)) {
1107 continue;
1108 }
1109 if (mcTruth.isValid()) { // if not set, -1 will be stored
1110 labelHolder.labelID = (mToStore[mcTruth.getSourceID()][mcTruth.getEventID()])[mcTruth.getTrackID()];
1111 }
1112 if (mcTruth.isFake()) {
1113 labelHolder.fwdLabelMask |= (0x1 << 7);
1114 }
1115 if (mcTruth.isNoise()) {
1116 labelHolder.fwdLabelMask |= (0x1 << 6);
1117 }
1118 if (src == GIndex::Source::MFT) {
1119 mcMFTTrackLabelCursor(labelHolder.labelID,
1120 labelHolder.fwdLabelMask);
1121 } else {
1122 mcFwdTrackLabelCursor(labelHolder.labelID,
1123 labelHolder.fwdLabelMask);
1124 }
1125 } else {
1126 if (!needToStore(mGIDToTableID)) {
1127 continue;
1128 }
1129 if (mcTruth.isValid()) { // if not set, -1 will be stored
1130 labelHolder.labelID = (mToStore[mcTruth.getSourceID()][mcTruth.getEventID()])[mcTruth.getTrackID()]; // defined by TPC if it contributes, otherwise: by ITS
1131 if (mcTruth.isFake()) {
1132 labelHolder.labelMask |= (0x1 << 15);
1133 }
1134 if (trackIndex.includesDet(DetID::TPC) && trackIndex.getSource() != GIndex::Source::TPC) { // this is global track
1135 auto contributorsGID = data.getSingleDetectorRefs(trackIndex);
1136 if (contributorsGID[GIndex::Source::ITSTPC].isIndexSet()) { // there is a match to ITS tracks or ITSAB tracklet!
1137 if (data.getTrackMCLabel(contributorsGID[GIndex::Source::ITSTPC]).isFake()) {
1138 labelHolder.labelMask |= (0x1 << 13);
1139 }
1140 }
1141 }
1142 } else if (mcTruth.isNoise()) {
1143 labelHolder.labelMask |= (0x1 << 14);
1144 }
1145 mcTrackLabelCursor(labelHolder.labelID, labelHolder.labelMask);
1146 }
1147 }
1148 }
1149 }
1150
1151 // filling the tables with the strangeness tracking labels
1152 auto sTrackLabels = data.getStrangeTracksMCLabels();
1153 // check if vertexId and vertexId + 1 maps into mVertexStrLUT
1154 if (!(vertexId < 0 || vertexId >= mVertexStrLUT.size() - 1)) {
1155 mcTrackLabelCursor.reserve(mVertexStrLUT[vertexId + 1] + mcTrackLabelCursor.lastIndex());
1156 for (int iS{mVertexStrLUT[vertexId]}; iS < mVertexStrLUT[vertexId + 1]; ++iS) {
1157 auto& collStrTrk = mCollisionStrTrk[iS];
1158 auto& label = sTrackLabels[collStrTrk.second];
1159 MCLabels labelHolder;
1160 labelHolder.labelID = label.isValid() ? (mToStore[label.getSourceID()][label.getEventID()])[label.getTrackID()] : -1;
1161 labelHolder.labelMask = (label.isFake() << 15) | (label.isNoise() << 14);
1162 mcTrackLabelCursor(labelHolder.labelID, labelHolder.labelMask);
1163 }
1164 }
1165}
1166
1167template <typename V0CursorType, typename CascadeCursorType, typename Decay3BodyCursorType>
1168void AODProducerWorkflowDPL::fillSecondaryVertices(const o2::globaltracking::RecoContainer& recoData, V0CursorType& v0Cursor, CascadeCursorType& cascadeCursor, Decay3BodyCursorType& decay3BodyCursor)
1169{
1170
1171 auto v0s = recoData.getV0sIdx();
1172 auto cascades = recoData.getCascadesIdx();
1173 auto decays3Body = recoData.getDecays3BodyIdx();
1174
1175 v0Cursor.reserve(v0s.size());
1176 // filling v0s table
1177 for (size_t iv0 = 0; iv0 < v0s.size(); iv0++) {
1178 const auto& v0 = v0s[iv0];
1179 auto trPosID = v0.getProngID(0);
1180 auto trNegID = v0.getProngID(1);
1181 uint8_t v0flags = v0.getBits();
1182 int posTableIdx = -1, negTableIdx = -1, collID = -1;
1183 auto item = mGIDToTableID.find(trPosID);
1184 if (item != mGIDToTableID.end()) {
1185 posTableIdx = item->second;
1186 } else {
1187 LOG(warn) << "Could not find a positive track index for prong ID " << trPosID;
1188 }
1189 item = mGIDToTableID.find(trNegID);
1190 if (item != mGIDToTableID.end()) {
1191 negTableIdx = item->second;
1192 } else {
1193 LOG(warn) << "Could not find a negative track index for prong ID " << trNegID;
1194 }
1195 auto itemV = mVtxToTableCollID.find(v0.getVertexID());
1196 if (itemV == mVtxToTableCollID.end()) {
1197 LOG(warn) << "Could not find V0 collisionID for the vertex ID " << v0.getVertexID();
1198 } else {
1199 collID = itemV->second;
1200 }
1201 if (posTableIdx != -1 and negTableIdx != -1 and collID != -1) {
1202 v0Cursor(collID, posTableIdx, negTableIdx, v0flags);
1203 mV0ToTableID[int(iv0)] = mTableV0ID++;
1204 }
1205 }
1206
1207 // filling cascades table
1208 cascadeCursor.reserve(cascades.size());
1209 for (auto& cascade : cascades) {
1210 auto itemV0 = mV0ToTableID.find(cascade.getV0ID());
1211 if (itemV0 == mV0ToTableID.end()) {
1212 continue;
1213 }
1214 int v0tableID = itemV0->second, bachTableIdx = -1, collID = -1;
1215 auto bachelorID = cascade.getBachelorID();
1216 auto item = mGIDToTableID.find(bachelorID);
1217 if (item != mGIDToTableID.end()) {
1218 bachTableIdx = item->second;
1219 } else {
1220 LOG(warn) << "Could not find a bachelor track index";
1221 continue;
1222 }
1223 auto itemV = mVtxToTableCollID.find(cascade.getVertexID());
1224 if (itemV != mVtxToTableCollID.end()) {
1225 collID = itemV->second;
1226 } else {
1227 LOG(warn) << "Could not find cascade collisionID for the vertex ID " << cascade.getVertexID();
1228 continue;
1229 }
1230 cascadeCursor(collID, v0tableID, bachTableIdx);
1231 }
1232
1233 // filling 3 body decays table
1234 decay3BodyCursor.reserve(decays3Body.size());
1235 for (size_t i3Body = 0; i3Body < decays3Body.size(); i3Body++) {
1236 const auto& decay3Body = decays3Body[i3Body];
1237 GIndex trIDs[3]{
1238 decay3Body.getProngID(0),
1239 decay3Body.getProngID(1),
1240 decay3Body.getProngID(2)};
1241 int tableIdx[3]{-1, -1, -1}, collID = -1;
1242 bool missing{false};
1243 for (int i{0}; i < 3; ++i) {
1244 auto item = mGIDToTableID.find(trIDs[i]);
1245 if (item != mGIDToTableID.end()) {
1246 tableIdx[i] = item->second;
1247 } else {
1248 LOG(warn) << fmt::format("Could not find a track index for prong ID {}", (int)trIDs[i]);
1249 missing = true;
1250 }
1251 }
1252 auto itemV = mVtxToTableCollID.find(decay3Body.getVertexID());
1253 if (itemV == mVtxToTableCollID.end()) {
1254 LOG(warn) << "Could not find 3 body collisionID for the vertex ID " << decay3Body.getVertexID();
1255 missing = true;
1256 } else {
1257 collID = itemV->second;
1258 }
1259 if (missing) {
1260 continue;
1261 }
1262 decay3BodyCursor(collID, tableIdx[0], tableIdx[1], tableIdx[2]);
1263 }
1264}
1265
1266template <typename FwdTrkClsCursorType>
1267void AODProducerWorkflowDPL::addClustersToFwdTrkClsTable(const o2::globaltracking::RecoContainer& recoData, FwdTrkClsCursorType& fwdTrkClsCursor, GIndex trackID, int fwdTrackId)
1268{
1269 const auto& mchTracks = recoData.getMCHTracks();
1270 const auto& mchmidMatches = recoData.getMCHMIDMatches();
1271 const auto& mchClusters = recoData.getMCHTrackClusters();
1272
1273 int mchTrackID = -1;
1274 if (trackID.getSource() == GIndex::MCH) { // This is an MCH track
1275 mchTrackID = trackID.getIndex();
1276 } else if (trackID.getSource() == GIndex::MCHMID) { // This is an MCH-MID track
1277 auto mchmidMatch = mchmidMatches[trackID.getIndex()];
1278 mchTrackID = mchmidMatch.getMCHRef().getIndex();
1279 } // Others are Global Forward Tracks, their clusters will be or were added with the corresponding MCH track
1280
1281 if (mchTrackID > -1 && mchTrackID < mchTracks.size()) {
1282 const auto& mchTrack = mchTracks[mchTrackID];
1283 fwdTrkClsCursor.reserve(mchTrack.getNClusters() + fwdTrkClsCursor.lastIndex());
1284 int first = mchTrack.getFirstClusterIdx();
1285 int last = mchTrack.getLastClusterIdx();
1286 for (int i = first; i <= last; i++) {
1287 const auto& cluster = mchClusters[i];
1288 fwdTrkClsCursor(fwdTrackId,
1289 truncateFloatFraction(cluster.x, mMuonCl),
1290 truncateFloatFraction(cluster.y, mMuonCl),
1291 truncateFloatFraction(cluster.z, mMuonCl),
1292 (((cluster.ey < 5.) & 0x1) << 12) | (((cluster.ex < 5.) & 0x1) << 11) | cluster.getDEId());
1293 }
1294 }
1295}
1296
1297template <typename HMPCursorType>
1298void AODProducerWorkflowDPL::fillHMPID(const o2::globaltracking::RecoContainer& recoData, HMPCursorType& hmpCursor)
1299{
1300 auto hmpMatches = recoData.getHMPMatches();
1301
1302 hmpCursor.reserve(hmpMatches.size());
1303
1304 // filling HMPs table
1305 for (size_t iHmp = 0; iHmp < hmpMatches.size(); iHmp++) {
1306
1307 const auto& match = hmpMatches[iHmp];
1308
1309 float xTrk, yTrk, theta, phi;
1310 float xMip, yMip;
1311 int charge, nph;
1312
1313 match.getHMPIDtrk(xTrk, yTrk, theta, phi);
1314 match.getHMPIDmip(xMip, yMip, charge, nph);
1315
1316 auto photChargeVec = match.getPhotCharge();
1317
1318 float photChargeVec2[10]; // = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
1319
1320 for (Int_t i = 0; i < 10; i++) {
1321 photChargeVec2[i] = photChargeVec[i];
1322 }
1323 auto tref = mGIDToTableID.find(match.getTrackRef());
1324 if (tref != mGIDToTableID.end()) {
1325 hmpCursor(tref->second, match.getHMPsignal(), xTrk, yTrk, xMip, yMip, nph, charge, match.getIdxHMPClus(), match.getHmpMom(), photChargeVec2);
1326 } else {
1327 LOG(error) << "Could not find AOD track table entry for HMP-matched track " << match.getTrackRef().asString();
1328 }
1329 }
1330}
1331
1332void AODProducerWorkflowDPL::prepareStrangenessTracking(const o2::globaltracking::RecoContainer& recoData)
1333{
1334 auto v0s = recoData.getV0sIdx();
1335 auto cascades = recoData.getCascadesIdx();
1336 auto decays3Body = recoData.getDecays3BodyIdx();
1337
1338 int sTrkID = 0;
1339 mCollisionStrTrk.clear();
1340 mCollisionStrTrk.reserve(recoData.getStrangeTracks().size());
1341 mVertexStrLUT.clear();
1342 mVertexStrLUT.resize(recoData.getPrimaryVertices().size() + 1, 0);
1343 for (auto& sTrk : recoData.getStrangeTracks()) {
1344 auto ITSIndex = GIndex{sTrk.mITSRef, GIndex::ITS};
1345 int vtxId{0};
1346 if (sTrk.mPartType == dataformats::kStrkV0) {
1347 vtxId = v0s[sTrk.mDecayRef].getVertexID();
1348 } else if (sTrk.mPartType == dataformats::kStrkCascade) {
1349 vtxId = cascades[sTrk.mDecayRef].getVertexID();
1350 } else {
1351 vtxId = decays3Body[sTrk.mDecayRef].getVertexID();
1352 }
1353 mCollisionStrTrk.emplace_back(vtxId, sTrkID++);
1354 mVertexStrLUT[vtxId]++;
1355 }
1356 std::exclusive_scan(mVertexStrLUT.begin(), mVertexStrLUT.end(), mVertexStrLUT.begin(), 0);
1357
1358 // sort by collision ID
1359 std::sort(mCollisionStrTrk.begin(), mCollisionStrTrk.end(), [](const auto& a, const auto& b) { return a.first < b.first; });
1360 mStrTrkIndices.clear();
1361 mStrTrkIndices.resize(mCollisionStrTrk.size(), -1);
1362}
1363
1364template <typename V0C, typename CC, typename D3BC>
1365void AODProducerWorkflowDPL::fillStrangenessTrackingTables(const o2::globaltracking::RecoContainer& recoData, V0C& v0Curs, CC& cascCurs, D3BC& d3BodyCurs)
1366{
1367 int itsTableIdx = -1;
1368 int sTrkID = 0;
1369 int nV0 = 0;
1370 int nCasc = 0;
1371 int nD3Body = 0;
1372
1373 for (const auto& sTrk : recoData.getStrangeTracks()) {
1374 if (sTrk.mPartType == dataformats::kStrkV0) {
1375 nV0++;
1376 } else if (sTrk.mPartType == dataformats::kStrkCascade) {
1377 nCasc++;
1378 } else {
1379 nD3Body++;
1380 }
1381 }
1382
1383 v0Curs.reserve(nV0);
1384 cascCurs.reserve(nCasc);
1385 d3BodyCurs.reserve(nD3Body);
1386
1387 for (const auto& sTrk : recoData.getStrangeTracks()) {
1388 auto ITSIndex = GIndex{sTrk.mITSRef, GIndex::ITS};
1389 auto item = mGIDToTableID.find(ITSIndex);
1390 if (item != mGIDToTableID.end()) {
1391 itsTableIdx = item->second;
1392 } else {
1393 LOG(warn) << "Could not find a ITS strange track index " << ITSIndex;
1394 continue;
1395 }
1396 if (sTrk.mPartType == dataformats::kStrkV0) {
1397 v0Curs(mStrTrkIndices[sTrkID++],
1398 itsTableIdx,
1399 sTrk.mDecayRef,
1400 sTrk.mDecayVtx[0],
1401 sTrk.mDecayVtx[1],
1402 sTrk.mDecayVtx[2],
1403 sTrk.mMasses[0],
1404 sTrk.mMasses[1],
1405 sTrk.mMatchChi2,
1406 sTrk.mTopoChi2,
1407 sTrk.getAverageClusterSize());
1408 } else if (sTrk.mPartType == dataformats::kStrkCascade) {
1409 cascCurs(mStrTrkIndices[sTrkID++],
1410 itsTableIdx,
1411 sTrk.mDecayRef,
1412 sTrk.mDecayVtx[0],
1413 sTrk.mDecayVtx[1],
1414 sTrk.mDecayVtx[2],
1415 sTrk.mMasses[0],
1416 sTrk.mMasses[1],
1417 sTrk.mMatchChi2,
1418 sTrk.mTopoChi2,
1419 sTrk.getAverageClusterSize());
1420 } else {
1421 d3BodyCurs(mStrTrkIndices[sTrkID++],
1422 itsTableIdx,
1423 sTrk.mDecayRef,
1424 sTrk.mDecayVtx[0],
1425 sTrk.mDecayVtx[1],
1426 sTrk.mDecayVtx[2],
1427 sTrk.mMasses[0],
1428 sTrk.mMasses[1],
1429 sTrk.mMatchChi2,
1430 sTrk.mTopoChi2,
1431 sTrk.getAverageClusterSize());
1432 }
1433 }
1434}
1435
1436void AODProducerWorkflowDPL::countTPCClusters(const o2::globaltracking::RecoContainer& data)
1437{
1438 const auto& tpcTracks = data.getTPCTracks();
1439 const auto& tpcClusRefs = data.getTPCTracksClusterRefs();
1440 const auto& tpcClusShMap = data.clusterShMapTPC;
1441 const auto& tpcClusAcc = data.getTPCClusters();
1442 constexpr int maxRows = 152;
1443 constexpr int neighbour = 2;
1444 int ntr = tpcTracks.size();
1445 mTPCCounters.clear();
1446 mTPCCounters.resize(ntr);
1447#ifdef WITH_OPENMP
1448 int ngroup = std::min(50, std::max(1, ntr / mNThreads));
1449#pragma omp parallel for schedule(dynamic, ngroup) num_threads(mNThreads)
1450#endif
1451 for (int itr = 0; itr < ntr; itr++) {
1452 std::array<bool, maxRows> clMap{}, shMap{};
1453 uint8_t sectorIndex, rowIndex;
1454 uint32_t clusterIndex;
1455 auto& counters = mTPCCounters[itr];
1456 const auto& track = tpcTracks[itr];
1457 for (int i = 0; i < track.getNClusterReferences(); i++) {
1458 o2::tpc::TrackTPC::getClusterReference(tpcClusRefs, i, sectorIndex, rowIndex, clusterIndex, track.getClusterRef());
1459 unsigned int absoluteIndex = tpcClusAcc.clusterOffset[sectorIndex][rowIndex] + clusterIndex;
1460 clMap[rowIndex] = true;
1461 if (tpcClusShMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) {
1462 if (!shMap[rowIndex]) {
1463 counters.shared++;
1464 }
1465 shMap[rowIndex] = true;
1466 }
1467 }
1468 int last = -1;
1469 for (int i = 0; i < maxRows; i++) {
1470 if (clMap[i]) {
1471 counters.crossed++;
1472 counters.found++;
1473 last = i;
1474 } else if ((i - last) <= neighbour) {
1475 counters.crossed++;
1476 } else {
1477 int lim = std::min(i + 1 + neighbour, maxRows);
1478 for (int j = i + 1; j < lim; j++) {
1479 if (clMap[j]) {
1480 counters.crossed++;
1481 break;
1482 }
1483 }
1484 }
1485 }
1486 }
1487}
1488
1489uint8_t AODProducerWorkflowDPL::getTRDPattern(const o2::trd::TrackTRD& track)
1490{
1491 uint8_t pattern = 0;
1492 for (int il = o2::trd::TrackTRD::EGPUTRDTrack::kNLayers - 1; il >= 0; il--) {
1493 if (track.getTrackletIndex(il) != -1) {
1494 pattern |= 0x1 << il;
1495 }
1496 }
1497 if (track.getHasNeighbor()) {
1498 pattern |= 0x1 << 6;
1499 }
1500 if (track.getHasPadrowCrossing()) {
1501 pattern |= 0x1 << 7;
1502 }
1503 return pattern;
1504}
1505
1506template <typename TCaloHandler, typename TCaloCursor, typename TCaloTRGCursor, typename TMCCaloLabelCursor>
1507void AODProducerWorkflowDPL::addToCaloTable(TCaloHandler& caloHandler, TCaloCursor& caloCellCursor, TCaloTRGCursor& caloTRGCursor,
1508 TMCCaloLabelCursor& mcCaloCellLabelCursor, int eventID, int bcID, int8_t caloType)
1509{
1510 auto inputEvent = caloHandler.buildEvent(eventID);
1511 auto cellsInEvent = inputEvent.mCells; // get cells belonging to current event
1512 auto cellMClabels = inputEvent.mMCCellLabels; // get MC labels belonging to current event (only implemented for EMCal currently!)
1513 caloCellCursor.reserve(cellsInEvent.size() + caloCellCursor.lastIndex());
1514 caloTRGCursor.reserve(cellsInEvent.size() + caloTRGCursor.lastIndex());
1515 if (mUseMC) {
1516 mcCaloCellLabelCursor.reserve(cellsInEvent.size() + mcCaloCellLabelCursor.lastIndex());
1517 }
1518 for (auto iCell = 0U; iCell < cellsInEvent.size(); iCell++) {
1519 caloCellCursor(bcID,
1520 CellHelper::getCellNumber(cellsInEvent[iCell]),
1521 truncateFloatFraction(CellHelper::getAmplitude(cellsInEvent[iCell]), mCaloAmp),
1522 truncateFloatFraction(CellHelper::getTimeStamp(cellsInEvent[iCell]), mCaloTime),
1523 cellsInEvent[iCell].getType(),
1524 caloType); // 1 = emcal, -1 = undefined, 0 = phos
1525
1526 // todo: fix dummy values in CellHelper when it is clear what is filled for trigger information
1527 if (CellHelper::isTRU(cellsInEvent[iCell])) { // Write only trigger cells into this table
1528 caloTRGCursor(bcID,
1529 CellHelper::getFastOrAbsID(cellsInEvent[iCell]),
1530 CellHelper::getLnAmplitude(cellsInEvent[iCell]),
1531 CellHelper::getTriggerBits(cellsInEvent[iCell]),
1532 caloType);
1533 }
1534 if (mUseMC) {
1535 // Common for PHOS and EMCAL
1536 // loop over all MC Labels for the current cell
1537 std::vector<int32_t> particleIds;
1538 std::vector<float> amplitudeFraction;
1539 if (!mEMCselectLeading) {
1540 particleIds.reserve(cellMClabels.size());
1541 amplitudeFraction.reserve(cellMClabels.size());
1542 }
1543 float tmpMaxAmplitude = 0;
1544 int32_t tmpindex = 0;
1545 for (auto& mclabel : cellMClabels[iCell]) {
1546 // do not fill noise lables!
1547 if (mclabel.isValid()) {
1548 if (mEMCselectLeading) {
1549 if (mclabel.getAmplitudeFraction() > tmpMaxAmplitude) {
1550 // Check if this MCparticle added to be kept?
1551 if (mToStore.at(mclabel.getSourceID()).at(mclabel.getEventID()).find(mclabel.getTrackID()) !=
1552 mToStore.at(mclabel.getSourceID()).at(mclabel.getEventID()).end()) {
1553 tmpMaxAmplitude = mclabel.getAmplitudeFraction();
1554 tmpindex = (mToStore.at(mclabel.getSourceID()).at(mclabel.getEventID())).at(mclabel.getTrackID());
1555 }
1556 }
1557 } else {
1558 auto trackStore = mToStore.at(mclabel.getSourceID()).at(mclabel.getEventID());
1559 auto iter = trackStore.find(mclabel.getTrackID());
1560 if (iter != trackStore.end()) {
1561 amplitudeFraction.emplace_back(mclabel.getAmplitudeFraction());
1562 particleIds.emplace_back(iter->second);
1563 } else {
1564 particleIds.emplace_back(-1); // should the mc particle not be in mToStore make sure something (e.g. -1) is saved in particleIds so the length of particleIds is the same es amplitudeFraction!
1565 amplitudeFraction.emplace_back(0.f);
1566 LOG(warn) << "CaloTable: Could not find track for mclabel (" << mclabel.getSourceID() << "," << mclabel.getEventID() << "," << mclabel.getTrackID() << ") in the AOD MC store";
1567 if (mMCKineReader) {
1568 auto mctrack = mMCKineReader->getTrack(mclabel);
1569 TVector3 vec;
1570 mctrack->GetStartVertex(vec);
1571 LOG(warn) << " ... this track is of PDG " << mctrack->GetPdgCode() << " produced by " << mctrack->getProdProcessAsString() << " at (" << vec.X() << "," << vec.Y() << "," << vec.Z() << ")";
1572 }
1573 }
1574 }
1575 }
1576 } // end of loop over all MC Labels for the current cell
1577 if (mEMCselectLeading) {
1578 amplitudeFraction.emplace_back(tmpMaxAmplitude);
1579 particleIds.emplace_back(tmpindex);
1580 }
1581 if (particleIds.size() == 0) {
1582 particleIds.emplace_back(-1);
1583 amplitudeFraction.emplace_back(0.f);
1584 }
1585 mcCaloCellLabelCursor(particleIds,
1586 amplitudeFraction);
1587 }
1588 } // end of loop over cells in current event
1589}
1590
1591// fill calo related tables (cells and calotrigger table)
1592template <typename TCaloCursor, typename TCaloTRGCursor, typename TMCCaloLabelCursor>
1593void AODProducerWorkflowDPL::fillCaloTable(TCaloCursor& caloCellCursor, TCaloTRGCursor& caloTRGCursor,
1594 TMCCaloLabelCursor& mcCaloCellLabelCursor, const std::map<uint64_t, int>& bcsMap,
1596{
1597 // get calo information
1598 auto caloEMCCells = data.getEMCALCells();
1599 auto caloEMCCellsTRGR = data.getEMCALTriggers();
1600 auto mcCaloEMCCellLabels = data.getEMCALCellsMCLabels();
1601
1602 auto caloPHOSCells = data.getPHOSCells();
1603 auto caloPHOSCellsTRGR = data.getPHOSTriggers();
1604 auto mcCaloPHOSCellLabels = data.getPHOSCellsMCLabels();
1605
1606 if (!mInputSources[GIndex::PHS]) {
1607 caloPHOSCells = {};
1608 caloPHOSCellsTRGR = {};
1609 mcCaloPHOSCellLabels = {};
1610 }
1611
1612 if (!mInputSources[GIndex::EMC]) {
1613 caloEMCCells = {};
1614 caloEMCCellsTRGR = {};
1615 mcCaloEMCCellLabels = {};
1616 }
1617
1620
1621 // get cell belonging to an eveffillnt instead of timeframe
1622 emcEventHandler.reset();
1623 emcEventHandler.setCellData(caloEMCCells, caloEMCCellsTRGR);
1624 emcEventHandler.setCellMCTruthContainer(mcCaloEMCCellLabels);
1625
1626 phsEventHandler.reset();
1627 phsEventHandler.setCellData(caloPHOSCells, caloPHOSCellsTRGR);
1628 phsEventHandler.setCellMCTruthContainer(mcCaloPHOSCellLabels);
1629
1630 int emcNEvents = emcEventHandler.getNumberOfEvents();
1631 int phsNEvents = phsEventHandler.getNumberOfEvents();
1632
1633 std::vector<std::tuple<uint64_t, int8_t, int>> caloEvents; // <bc, caloType, eventID>
1634
1635 caloEvents.reserve(emcNEvents + phsNEvents);
1636
1637 for (int iev = 0; iev < emcNEvents; ++iev) {
1638 uint64_t bc = emcEventHandler.getInteractionRecordForEvent(iev).toLong();
1639 caloEvents.emplace_back(std::make_tuple(bc, 1, iev));
1640 }
1641
1642 for (int iev = 0; iev < phsNEvents; ++iev) {
1643 uint64_t bc = phsEventHandler.getInteractionRecordForEvent(iev).toLong();
1644 caloEvents.emplace_back(std::make_tuple(bc, 0, iev));
1645 }
1646
1647 std::sort(caloEvents.begin(), caloEvents.end(),
1648 [](const auto& left, const auto& right) { return std::get<0>(left) < std::get<0>(right); });
1649
1650 // loop over events
1651 for (int i = 0; i < emcNEvents + phsNEvents; ++i) {
1652 uint64_t globalBC = std::get<0>(caloEvents[i]);
1653 int8_t caloType = std::get<1>(caloEvents[i]);
1654 int eventID = std::get<2>(caloEvents[i]);
1655 auto item = bcsMap.find(globalBC);
1656 int bcID = -1;
1657 if (item != bcsMap.end()) {
1658 bcID = item->second;
1659 } else {
1660 LOG(warn) << "Error: could not find a corresponding BC ID for a calo point; globalBC = " << globalBC << ", caloType = " << (int)caloType;
1661 }
1662 if (caloType == 0) { // phos
1663 addToCaloTable(phsEventHandler, caloCellCursor, caloTRGCursor, mcCaloCellLabelCursor, eventID, bcID, caloType);
1664 }
1665 if (caloType == 1) { // emc
1666 addToCaloTable(emcEventHandler, caloCellCursor, caloTRGCursor, mcCaloCellLabelCursor, eventID, bcID, caloType);
1667 }
1668 }
1669
1670 caloEvents.clear();
1671}
1672
1674{
1675 mTimer.Stop();
1677 mLPMProdTag = ic.options().get<string>("lpmp-prod-tag");
1678 mAnchorPass = ic.options().get<string>("anchor-pass");
1679 mAnchorProd = ic.options().get<string>("anchor-prod");
1680 mUser = ic.options().get<string>("created-by");
1681 mRecoPass = ic.options().get<string>("reco-pass");
1682 mTFNumber = ic.options().get<int64_t>("aod-timeframe-id");
1683 mRecoOnly = ic.options().get<int>("reco-mctracks-only");
1684 mTruncate = ic.options().get<int>("enable-truncation");
1685 mRunNumber = ic.options().get<int>("run-number");
1686 mCTPReadout = ic.options().get<int>("ctpreadout-create");
1687 mNThreads = std::max(1, ic.options().get<int>("nthreads"));
1688 mEMCselectLeading = ic.options().get<bool>("emc-select-leading");
1689 mThinTracks = ic.options().get<bool>("thin-tracks");
1690 mPropTracks = ic.options().get<bool>("propagate-tracks");
1691 mPropMuons = ic.options().get<bool>("propagate-muons");
1692 if (auto s = ic.options().get<std::string>("with-streamers"); !s.empty()) {
1693 mStreamerFlags.set(s);
1694 if (mStreamerFlags) {
1695 LOGP(info, "Writing streamer data with mask:");
1696 LOG(info) << mStreamerFlags;
1697 } else {
1698 LOGP(warn, "Specified non-default empty streamer mask!");
1699 }
1700 }
1701 mTrackQCFraction = ic.options().get<float>("trackqc-fraction");
1702 mTrackQCNTrCut = ic.options().get<int64_t>("trackqc-NTrCut");
1703 mTrackQCDCAxy = ic.options().get<float>("trackqc-tpc-dca");
1704 mTrackQCPt = ic.options().get<float>("trackqc-tpc-pt");
1705 mTrackQCNCls = ic.options().get<int>("trackqc-tpc-cls");
1706 if (auto seed = ic.options().get<int>("seed"); seed == 0) {
1707 LOGP(info, "Using random device for seeding");
1708 std::random_device rd;
1709 std::array<int, std::mt19937::state_size> seed_data{};
1710 std::generate(std::begin(seed_data), std::end(seed_data), std::ref(rd));
1711 std::seed_seq seq(std::begin(seed_data), std::end(seed_data));
1712 mGenerator = std::mt19937(seq);
1713 } else {
1714 LOGP(info, "Using seed {} for sampling", seed);
1715 mGenerator.seed(seed);
1716 }
1717#ifdef WITH_OPENMP
1718 LOGP(info, "Multi-threaded parts will run with {} OpenMP threads", mNThreads);
1719#else
1720 mNThreads = 1;
1721 LOG(info) << "OpenMP is disabled";
1722#endif
1723 if (mTFNumber == -1L) {
1724 LOG(info) << "TFNumber will be obtained from CCDB";
1725 }
1726 if (mRunNumber == -1L) {
1727 LOG(info) << "The Run number will be obtained from DPL headers";
1728 }
1729
1730 // set no truncation if selected by user
1731 if (mTruncate != 1) {
1732 LOG(info) << "Truncation is not used!";
1733 mCollisionPosition = 0xFFFFFFFF;
1734 mCollisionPositionCov = 0xFFFFFFFF;
1735 mTrackX = 0xFFFFFFFF;
1736 mTrackAlpha = 0xFFFFFFFF;
1737 mTrackSnp = 0xFFFFFFFF;
1738 mTrackTgl = 0xFFFFFFFF;
1739 mTrack1Pt = 0xFFFFFFFF;
1740 mTrackChi2 = 0xFFFFFFFF;
1741 mTrackCovDiag = 0xFFFFFFFF;
1742 mTrackCovOffDiag = 0xFFFFFFFF;
1743 mTrackSignal = 0xFFFFFFFF;
1744 mTrackTime = 0xFFFFFFFF;
1745 mTPCTime0 = 0xFFFFFFFF;
1746 mTrackTimeError = 0xFFFFFFFF;
1747 mTrackPosEMCAL = 0xFFFFFFFF;
1748 mTracklets = 0xFFFFFFFF;
1749 mMcParticleW = 0xFFFFFFFF;
1750 mMcParticlePos = 0xFFFFFFFF;
1751 mMcParticleMom = 0xFFFFFFFF;
1752 mCaloAmp = 0xFFFFFFFF;
1753 mCaloTime = 0xFFFFFFFF;
1754 mCPVPos = 0xFFFFFFFF;
1755 mCPVAmpl = 0xFFFFFFFF;
1756 mMuonTr1P = 0xFFFFFFFF;
1757 mMuonTrThetaX = 0xFFFFFFFF;
1758 mMuonTrThetaY = 0xFFFFFFFF;
1759 mMuonTrZmu = 0xFFFFFFFF;
1760 mMuonTrBend = 0xFFFFFFFF;
1761 mMuonTrNonBend = 0xFFFFFFFF;
1762 mMuonTrCov = 0xFFFFFFFF;
1763 mMuonCl = 0xFFFFFFFF;
1764 mMuonClErr = 0xFFFFFFFF;
1765 mV0Time = 0xFFFFFFFF;
1766 mV0ChannelTime = 0xFFFFFFFF;
1767 mFDDTime = 0xFFFFFFFF;
1768 mFDDChannelTime = 0xFFFFFFFF;
1769 mT0Time = 0xFFFFFFFF;
1770 mT0ChannelTime = 0xFFFFFFFF;
1771 mV0Amplitude = 0xFFFFFFFF;
1772 mFDDAmplitude = 0xFFFFFFFF;
1773 mT0Amplitude = 0xFFFFFFFF;
1774 }
1775 // Initialize ZDC helper maps
1776 for (int ic = 0; ic < o2::zdc::NChannels; ic++) {
1777 mZDCEnergyMap[ic] = -std::numeric_limits<float>::infinity();
1778 }
1779 for (int ic = 0; ic < o2::zdc::NTDCChannels; ic++) {
1780 mZDCTDCMap[ic] = -std::numeric_limits<float>::infinity();
1781 }
1782
1783 std::string hepmcUpdate = ic.options().get<std::string>("hepmc-update");
1784 HepMCUpdate when = (hepmcUpdate == "never" ? HepMCUpdate::never : hepmcUpdate == "always" ? HepMCUpdate::always
1785 : hepmcUpdate == "all" ? HepMCUpdate::allKeys
1786 : HepMCUpdate::anyKey);
1787 mXSectionUpdate = when;
1788 mPdfInfoUpdate = when;
1789 mHeavyIonUpdate = when;
1790
1791 mTimer.Reset();
1792
1793 if (mStreamerFlags) {
1794 mStreamer = std::make_unique<o2::utils::TreeStreamRedirector>("AO2DStreamer.root", "RECREATE");
1795 }
1796}
1797
1798namespace
1799{
1800void add_additional_meta_info(std::vector<TString>& keys, std::vector<TString>& values)
1801{
1802 // see if we should put additional meta info (e.g. from MC)
1803 auto aod_external_meta_info_file = getenv("AOD_ADDITIONAL_METADATA_FILE");
1804 if (aod_external_meta_info_file != nullptr) {
1805 LOG(info) << "Trying to inject additional AOD meta-data from " << aod_external_meta_info_file;
1806 if (std::filesystem::exists(aod_external_meta_info_file)) {
1807 std::ifstream input_file(aod_external_meta_info_file);
1808 if (input_file) {
1809 nlohmann::json json_data;
1810 try {
1811 input_file >> json_data;
1812 } catch (nlohmann::json::parse_error& e) {
1813 std::cerr << "JSON Parse Error: " << e.what() << "\n";
1814 std::cerr << "Exception ID: " << e.id << "\n";
1815 std::cerr << "Byte position: " << e.byte << "\n";
1816 return;
1817 }
1818 // If parsing succeeds, iterate over key-value pairs
1819 for (const auto& [key, value] : json_data.items()) {
1820 LOG(info) << "Adding AOD MetaData" << key << " : " << value;
1821 keys.push_back(key.c_str());
1822 values.push_back(value.get<std::string>());
1823 }
1824 }
1825 }
1826 }
1827}
1828} // namespace
1829
1831{
1832 mTimer.Start(false);
1834 recoData.collectData(pc, *mDataRequest);
1835 updateTimeDependentParams(pc); // Make sure that this is called after the RecoContainer collect data, since some condition objects are fetched there
1836
1837 mStartIR = recoData.startIR;
1838
1839 auto primVertices = recoData.getPrimaryVertices();
1840 auto primVer2TRefs = recoData.getPrimaryVertexMatchedTrackRefs();
1841 auto primVerGIs = recoData.getPrimaryVertexMatchedTracks();
1842 auto primVerLabels = recoData.getPrimaryVertexMCLabels();
1843
1844 auto fddChData = recoData.getFDDChannelsData();
1845 auto fddRecPoints = recoData.getFDDRecPoints();
1846 auto ft0ChData = recoData.getFT0ChannelsData();
1847 auto ft0RecPoints = recoData.getFT0RecPoints();
1848 auto fv0ChData = recoData.getFV0ChannelsData();
1849 auto fv0RecPoints = recoData.getFV0RecPoints();
1850
1851 auto zdcEnergies = recoData.getZDCEnergy();
1852 auto zdcBCRecData = recoData.getZDCBCRecData();
1853 auto zdcTDCData = recoData.getZDCTDCData();
1854
1855 auto cpvClusters = recoData.getCPVClusters();
1856 auto cpvTrigRecs = recoData.getCPVTriggers();
1857
1858 auto ctpDigits = recoData.getCTPDigits();
1859 const auto& tinfo = pc.services().get<o2::framework::TimingInfo>();
1860 std::vector<o2::ctp::CTPDigit> ctpDigitsCreated;
1861 if (mCTPReadout == 1) {
1862 LOG(info) << "CTP : creating ctpreadout in AOD producer";
1863 createCTPReadout(recoData, ctpDigitsCreated, pc);
1864 LOG(info) << "CTP : ctpreadout created from AOD";
1865 ctpDigits = gsl::span<o2::ctp::CTPDigit>(ctpDigitsCreated);
1866 }
1867 LOG(debug) << "FOUND " << primVertices.size() << " primary vertices";
1868 LOG(debug) << "FOUND " << ft0RecPoints.size() << " FT0 rec. points";
1869 LOG(debug) << "FOUND " << fv0RecPoints.size() << " FV0 rec. points";
1870 LOG(debug) << "FOUND " << fddRecPoints.size() << " FDD rec. points";
1871 LOG(debug) << "FOUND " << cpvClusters.size() << " CPV clusters";
1872 LOG(debug) << "FOUND " << cpvTrigRecs.size() << " CPV trigger records";
1873
1874 LOG(info) << "FOUND " << primVertices.size() << " primary vertices";
1875
1876 using namespace o2::aodhelpers;
1877
1878 auto bcCursor = createTableCursor<o2::aod::BCs>(pc);
1879 auto bcFlagsCursor = createTableCursor<o2::aod::BCFlags>(pc);
1880 auto cascadesCursor = createTableCursor<o2::aod::Cascades>(pc);
1881 auto collisionsCursor = createTableCursor<o2::aod::Collisions>(pc);
1882 auto decay3BodyCursor = createTableCursor<o2::aod::Decay3Bodys>(pc);
1883 auto trackedCascadeCursor = createTableCursor<o2::aod::TrackedCascades>(pc);
1884 auto trackedV0Cursor = createTableCursor<o2::aod::TrackedV0s>(pc);
1885 auto tracked3BodyCurs = createTableCursor<o2::aod::Tracked3Bodys>(pc);
1886 auto fddCursor = createTableCursor<o2::aod::FDDs>(pc);
1887 auto fddExtraCursor = createTableCursor<o2::aod::FDDsExtra>(pc);
1888 auto ft0Cursor = createTableCursor<o2::aod::FT0s>(pc);
1889 auto ft0ExtraCursor = createTableCursor<o2::aod::FT0sExtra>(pc);
1890 auto fv0aCursor = createTableCursor<o2::aod::FV0As>(pc);
1891 auto fv0aExtraCursor = createTableCursor<o2::aod::FV0AsExtra>(pc);
1892 auto fwdTracksCursor = createTableCursor<o2::aod::StoredFwdTracks>(pc);
1893 auto fwdTracksCovCursor = createTableCursor<o2::aod::StoredFwdTracksCov>(pc);
1894 auto fwdTrkClsCursor = createTableCursor<o2::aod::FwdTrkCls>(pc);
1895 auto mftTracksCursor = createTableCursor<o2::aod::StoredMFTTracks>(pc);
1896 auto mftTracksCovCursor = createTableCursor<o2::aod::StoredMFTTracksCov>(pc);
1897 auto tracksCursor = createTableCursor<o2::aod::StoredTracksIU>(pc);
1898 auto tracksCovCursor = createTableCursor<o2::aod::StoredTracksCovIU>(pc);
1899 auto tracksExtraCursor = createTableCursor<o2::aod::StoredTracksExtra>(pc);
1900 auto tracksQACursor = createTableCursor<o2::aod::TracksQAVersion>(pc);
1901 auto ambigTracksCursor = createTableCursor<o2::aod::AmbiguousTracks>(pc);
1902 auto ambigMFTTracksCursor = createTableCursor<o2::aod::AmbiguousMFTTracks>(pc);
1903 auto ambigFwdTracksCursor = createTableCursor<o2::aod::AmbiguousFwdTracks>(pc);
1904 auto v0sCursor = createTableCursor<o2::aod::V0s>(pc);
1905 auto zdcCursor = createTableCursor<o2::aod::Zdcs>(pc);
1906 auto hmpCursor = createTableCursor<o2::aod::HMPIDs>(pc);
1907 auto caloCellsCursor = createTableCursor<o2::aod::Calos>(pc);
1908 auto caloCellsTRGTableCursor = createTableCursor<o2::aod::CaloTriggers>(pc);
1909 auto cpvClustersCursor = createTableCursor<o2::aod::CPVClusters>(pc);
1910 auto originCursor = createTableCursor<o2::aod::Origins>(pc);
1911
1912 // Declare MC cursors type without adding the output for a table
1923 if (mUseMC) { // This creates the actual writercursor
1924 mcColLabelsCursor = createTableCursor<o2::aod::McCollisionLabels>(pc);
1925 mcCollisionsCursor = createTableCursor<o2::aod::McCollisions>(pc);
1926 hepmcXSectionsCursor = createTableCursor<o2::aod::HepMCXSections>(pc);
1927 hepmcPdfInfosCursor = createTableCursor<o2::aod::HepMCPdfInfos>(pc);
1928 hepmcHeavyIonsCursor = createTableCursor<o2::aod::HepMCHeavyIons>(pc);
1929 mcMFTTrackLabelCursor = createTableCursor<o2::aod::McMFTTrackLabels>(pc);
1930 mcFwdTrackLabelCursor = createTableCursor<o2::aod::McFwdTrackLabels>(pc);
1931 mcParticlesCursor = createTableCursor<o2::aod::StoredMcParticles_001>(pc);
1932 mcTrackLabelCursor = createTableCursor<o2::aod::McTrackLabels>(pc);
1933 mcCaloLabelsCursor = createTableCursor<o2::aod::McCaloLabels_001>(pc);
1934 }
1935
1936 std::unique_ptr<o2::steer::MCKinematicsReader> mcReader;
1937 if (mUseMC) {
1938 mcReader = std::make_unique<o2::steer::MCKinematicsReader>("collisioncontext.root");
1939 }
1940 mMCKineReader = mcReader.get(); // for use in different functions
1941 std::map<uint64_t, int> bcsMap;
1942 collectBCs(recoData, mUseMC ? mcReader->getDigitizationContext()->getEventRecords() : std::vector<o2::InteractionTimeRecord>{}, bcsMap);
1943 if (!primVer2TRefs.empty()) { // if the vertexing was done, the last slot refers to orphan tracks
1944 addRefGlobalBCsForTOF(primVer2TRefs.back(), primVerGIs, recoData, bcsMap);
1945 }
1946 // initialize the bunch crossing container for further use below
1947 mBCLookup.init(bcsMap);
1948
1949 uint64_t tfNumber;
1950 const int runNumber = (mRunNumber == -1) ? int(tinfo.runNumber) : mRunNumber;
1951 if (mTFNumber == -1L) {
1952 // TODO has to use absolute time of TF
1953 tfNumber = uint64_t(tinfo.firstTForbit) + (uint64_t(tinfo.runNumber) << 32); // getTFNumber(mStartIR, runNumber);
1954 } else {
1955 tfNumber = mTFNumber;
1956 }
1957
1958 std::vector<float> aAmplitudes, aTimes;
1959 std::vector<uint8_t> aChannels;
1960 fv0aCursor.reserve(fv0RecPoints.size());
1961 for (auto& fv0RecPoint : fv0RecPoints) {
1962 aAmplitudes.clear();
1963 aChannels.clear();
1964 aTimes.clear();
1965 const auto channelData = fv0RecPoint.getBunchChannelData(fv0ChData);
1966 for (auto& channel : channelData) {
1967 if (channel.charge > 0) {
1968 aAmplitudes.push_back(truncateFloatFraction(channel.charge, mV0Amplitude));
1969 aTimes.push_back(truncateFloatFraction(channel.time * 1.E-3, mV0ChannelTime));
1970 aChannels.push_back(channel.channel);
1971 }
1972 }
1973 uint64_t bc = fv0RecPoint.getInteractionRecord().toLong();
1974 auto item = bcsMap.find(bc);
1975 int bcID = -1;
1976 if (item != bcsMap.end()) {
1977 bcID = item->second;
1978 } else {
1979 LOG(fatal) << "Error: could not find a corresponding BC ID for a FV0 rec. point; BC = " << bc;
1980 }
1981 fv0aCursor(bcID,
1982 aAmplitudes,
1983 aChannels,
1984 truncateFloatFraction(fv0RecPoint.getCollisionGlobalMeanTime() * 1E-3, mV0Time), // ps to ns
1985 fv0RecPoint.getTrigger().getTriggersignals());
1986
1987 if (mEnableFITextra) {
1988 fv0aExtraCursor(bcID,
1989 aTimes);
1990 }
1991 }
1992
1993 std::vector<float> zdcEnergy, zdcAmplitudes, zdcTime;
1994 std::vector<uint8_t> zdcChannelsE, zdcChannelsT;
1995 zdcCursor.reserve(zdcBCRecData.size());
1996 for (auto zdcRecData : zdcBCRecData) {
1997 uint64_t bc = zdcRecData.ir.toLong();
1998 auto item = bcsMap.find(bc);
1999 int bcID = -1;
2000 if (item != bcsMap.end()) {
2001 bcID = item->second;
2002 } else {
2003 LOG(fatal) << "Error: could not find a corresponding BC ID for a ZDC rec. point; BC = " << bc;
2004 }
2005 int fe, ne, ft, nt, fi, ni;
2006 zdcRecData.getRef(fe, ne, ft, nt, fi, ni);
2007 zdcEnergy.clear();
2008 zdcChannelsE.clear();
2009 zdcAmplitudes.clear();
2010 zdcTime.clear();
2011 zdcChannelsT.clear();
2012 for (int ie = 0; ie < ne; ie++) {
2013 auto& zdcEnergyData = zdcEnergies[fe + ie];
2014 zdcEnergy.emplace_back(zdcEnergyData.energy());
2015 zdcChannelsE.emplace_back(zdcEnergyData.ch());
2016 }
2017 for (int it = 0; it < nt; it++) {
2018 auto& tdc = zdcTDCData[ft + it];
2019 zdcAmplitudes.emplace_back(tdc.amplitude());
2020 zdcTime.emplace_back(tdc.value());
2021 zdcChannelsT.emplace_back(o2::zdc::TDCSignal[tdc.ch()]);
2022 }
2023 zdcCursor(bcID,
2024 zdcEnergy,
2025 zdcChannelsE,
2026 zdcAmplitudes,
2027 zdcTime,
2028 zdcChannelsT);
2029 }
2030
2031 // keep track event/source id for each mc-collision
2032 // using map and not unordered_map to ensure
2033 // correct ordering when iterating over container elements
2034 std::vector<std::vector<int>> mcColToEvSrc;
2035
2036 if (mUseMC) {
2037 using namespace o2::aodmchelpers;
2038
2039 // filling mcCollision table
2040 int nMCCollisions = mcReader->getDigitizationContext()->getNCollisions();
2041 const auto& mcRecords = mcReader->getDigitizationContext()->getEventRecords();
2042 const auto& mcParts = mcReader->getDigitizationContext()->getEventParts();
2043
2044 // count all parts
2045 int totalNParts = 0;
2046 for (int iCol = 0; iCol < nMCCollisions; iCol++) {
2047 totalNParts += mcParts[iCol].size();
2048 }
2049 mcCollisionsCursor.reserve(totalNParts);
2050
2051 for (int iCol = 0; iCol < nMCCollisions; iCol++) {
2052 const auto time = mcRecords[iCol].getTimeOffsetWrtBC();
2053 auto globalBC = mcRecords[iCol].toLong();
2054 auto item = bcsMap.find(globalBC);
2055 int bcID = -1;
2056 if (item != bcsMap.end()) {
2057 bcID = item->second;
2058 } else {
2059 LOG(fatal) << "Error: could not find a corresponding BC ID "
2060 << "for MC collision; BC = " << globalBC
2061 << ", mc collision = " << iCol;
2062 }
2063 auto& colParts = mcParts[iCol];
2064 auto nParts = colParts.size();
2065 for (auto colPart : colParts) {
2066 auto eventID = colPart.entryID;
2067 auto sourceID = colPart.sourceID;
2068 // enable embedding: if several colParts exist, then they are
2069 // saved as one collision
2070 if (nParts == 1 || sourceID == 0) {
2071 // FIXME:
2072 // use generators' names for generatorIDs (?)
2073 auto& header = mcReader->getMCEventHeader(sourceID, eventID);
2074 updateMCHeader(mcCollisionsCursor.cursor,
2075 hepmcXSectionsCursor.cursor,
2076 hepmcPdfInfosCursor.cursor,
2077 hepmcHeavyIonsCursor.cursor,
2078 header,
2079 iCol,
2080 bcID,
2081 time,
2082 0,
2083 sourceID);
2084 }
2085 mcColToEvSrc.emplace_back(std::vector<int>{iCol, sourceID, eventID}); // point background and injected signal events to one collision
2086 }
2087 }
2088 }
2089
2090 std::sort(mcColToEvSrc.begin(), mcColToEvSrc.end(),
2091 [](const std::vector<int>& left, const std::vector<int>& right) { return (left[0] < right[0]); });
2092
2093 // vector of FDD amplitudes
2094 int16_t aFDDAmplitudesA[8] = {0u}, aFDDAmplitudesC[8] = {0u};
2095 float aFDDTimesA[8] = {0.f}, aFDDTimesC[8] = {0.f};
2096 // filling FDD table
2097 fddCursor.reserve(fddRecPoints.size());
2098 for (const auto& fddRecPoint : fddRecPoints) {
2099 for (int i = 0; i < 8; i++) {
2100 aFDDAmplitudesA[i] = 0;
2101 aFDDAmplitudesC[i] = 0;
2102 aFDDTimesA[i] = 0.f;
2103 aFDDTimesC[i] = 0.f;
2104 }
2105 uint64_t globalBC = fddRecPoint.getInteractionRecord().toLong();
2106 uint64_t bc = globalBC;
2107 auto item = bcsMap.find(bc);
2108 int bcID = -1;
2109 if (item != bcsMap.end()) {
2110 bcID = item->second;
2111 } else {
2112 LOG(fatal) << "Error: could not find a corresponding BC ID for a FDD rec. point; BC = " << bc;
2113 }
2114 const auto channelData = fddRecPoint.getBunchChannelData(fddChData);
2115 for (const auto& channel : channelData) {
2116 if (channel.mPMNumber < 8) {
2117 aFDDAmplitudesC[channel.mPMNumber] = channel.mChargeADC; // amplitude
2118 aFDDTimesC[channel.mPMNumber] = truncateFloatFraction(channel.mTime * 1E-3, mFDDChannelTime); // time
2119 } else {
2120 aFDDAmplitudesA[channel.mPMNumber - 8] = channel.mChargeADC; // amplitude
2121 aFDDTimesA[channel.mPMNumber - 8] = truncateFloatFraction(channel.mTime * 1E-3, mFDDChannelTime); // time
2122 }
2123 }
2124
2125 fddCursor(bcID,
2126 aFDDAmplitudesA,
2127 aFDDAmplitudesC,
2128 truncateFloatFraction(fddRecPoint.getCollisionTimeA() * 1E-3, mFDDTime), // ps to ns
2129 truncateFloatFraction(fddRecPoint.getCollisionTimeC() * 1E-3, mFDDTime), // ps to ns
2130 fddRecPoint.getTrigger().getTriggersignals());
2131 if (mEnableFITextra) {
2132 fddExtraCursor(bcID,
2133 aFDDTimesA,
2134 aFDDTimesC);
2135 }
2136 }
2137
2138 // filling FT0 table
2139 std::vector<float> aAmplitudesA, aAmplitudesC, aTimesA, aTimesC;
2140 std::vector<uint8_t> aChannelsA, aChannelsC;
2141 ft0Cursor.reserve(ft0RecPoints.size());
2142 for (auto& ft0RecPoint : ft0RecPoints) {
2143 aAmplitudesA.clear();
2144 aAmplitudesC.clear();
2145 aTimesA.clear();
2146 aTimesC.clear();
2147 aChannelsA.clear();
2148 aChannelsC.clear();
2149 const auto channelData = ft0RecPoint.getBunchChannelData(ft0ChData);
2150 for (auto& channel : channelData) {
2151 // TODO: switch to calibrated amplitude
2152 if (channel.QTCAmpl > 0) {
2153 constexpr int nFT0ChannelsAside = o2::ft0::Geometry::NCellsA * 4;
2154 if (channel.ChId < nFT0ChannelsAside) {
2155 aChannelsA.push_back(channel.ChId);
2156 aAmplitudesA.push_back(truncateFloatFraction(channel.QTCAmpl, mT0Amplitude));
2157 aTimesA.push_back(truncateFloatFraction(channel.CFDTime * 1E-3, mT0ChannelTime));
2158 } else {
2159 aChannelsC.push_back(channel.ChId - nFT0ChannelsAside);
2160 aAmplitudesC.push_back(truncateFloatFraction(channel.QTCAmpl, mT0Amplitude));
2161 aTimesC.push_back(truncateFloatFraction(channel.CFDTime * 1E-3, mT0ChannelTime));
2162 }
2163 }
2164 }
2165 uint64_t globalBC = ft0RecPoint.getInteractionRecord().toLong();
2166 uint64_t bc = globalBC;
2167 auto item = bcsMap.find(bc);
2168 int bcID = -1;
2169 if (item != bcsMap.end()) {
2170 bcID = item->second;
2171 } else {
2172 LOG(fatal) << "Error: could not find a corresponding BC ID for a FT0 rec. point; BC = " << bc;
2173 }
2174 ft0Cursor(bcID,
2175 aAmplitudesA,
2176 aChannelsA,
2177 aAmplitudesC,
2178 aChannelsC,
2179 truncateFloatFraction(ft0RecPoint.getCollisionTimeA() * 1E-3, mT0Time), // ps to ns
2180 truncateFloatFraction(ft0RecPoint.getCollisionTimeC() * 1E-3, mT0Time), // ps to ns
2181 ft0RecPoint.getTrigger().getTriggersignals());
2182 if (mEnableFITextra) {
2183 ft0ExtraCursor(bcID,
2184 aTimesA,
2185 aTimesC);
2186 }
2187 }
2188
2189 if (mUseMC) {
2190 // filling MC collision labels
2191 mcColLabelsCursor.reserve(primVerLabels.size());
2192 for (auto& label : primVerLabels) {
2193 auto it = std::find_if(mcColToEvSrc.begin(), mcColToEvSrc.end(),
2194 [&label](const auto& mcColInfo) { return mcColInfo[1] == label.getSourceID() && mcColInfo[2] == label.getEventID(); });
2195 int32_t mcCollisionID = -1;
2196 if (it != mcColToEvSrc.end()) {
2197 mcCollisionID = it->at(0);
2198 }
2199 uint16_t mcMask = 0; // todo: set mask using normalized weights?
2200 mcColLabelsCursor(mcCollisionID, mcMask);
2201 }
2202 }
2203
2204 cacheTriggers(recoData);
2205 countTPCClusters(recoData);
2206
2207 int collisionID = 0;
2208 mIndexTableMFT.resize(recoData.getMFTTracks().size());
2209 mIndexTableFwd.resize(recoData.getMCHTracks().size());
2210
2211 auto& trackReffwd = primVer2TRefs.back();
2212 fillIndexTablesPerCollision(trackReffwd, primVerGIs, recoData);
2213 collisionID = 0;
2214 for (auto& vertex : primVertices) {
2215 auto& trackReffwd = primVer2TRefs[collisionID];
2216 fillIndexTablesPerCollision(trackReffwd, primVerGIs, recoData); // this function must follow the same track order as 'fillTrackTablesPerCollision' to fill the map of track indices
2217 collisionID++;
2218 }
2219
2221 prepareStrangenessTracking(recoData);
2222
2223 mGIDToTableFwdID.clear(); // reset the tables to be used by 'fillTrackTablesPerCollision'
2224 mGIDToTableMFTID.clear();
2225
2226 if (mPropTracks || mThinTracks) {
2227 auto v0s = recoData.getV0sIdx();
2228 auto cascades = recoData.getCascadesIdx();
2229 auto decays3Body = recoData.getDecays3BodyIdx();
2230 mGIDUsedBySVtx.reserve(v0s.size() * 2 + cascades.size() + decays3Body.size() * 3);
2231 for (const auto& v0 : v0s) {
2232 mGIDUsedBySVtx.insert(v0.getProngID(0));
2233 mGIDUsedBySVtx.insert(v0.getProngID(1));
2234 }
2235 for (const auto& cascade : cascades) {
2236 mGIDUsedBySVtx.insert(cascade.getBachelorID());
2237 }
2238 for (const auto& id3Body : decays3Body) {
2239 mGIDUsedBySVtx.insert(id3Body.getProngID(0));
2240 mGIDUsedBySVtx.insert(id3Body.getProngID(1));
2241 mGIDUsedBySVtx.insert(id3Body.getProngID(2));
2242 }
2243
2244 mGIDUsedByStr.reserve(recoData.getStrangeTracks().size());
2245 for (const auto& sTrk : recoData.getStrangeTracks()) {
2246 mGIDUsedByStr.emplace(sTrk.mITSRef, GIndex::ITS);
2247 }
2248 }
2249
2250 // filling unassigned tracks first
2251 // so that all unassigned tracks are stored in the beginning of the table together
2252 auto& trackRef = primVer2TRefs.back(); // references to unassigned tracks are at the end
2253 // fixme: interaction time is undefined for unassigned tracks (?)
2254 fillTrackTablesPerCollision(-1, std::uint64_t(-1), trackRef, primVerGIs, recoData, tracksCursor, tracksCovCursor, tracksExtraCursor, tracksQACursor,
2255 ambigTracksCursor, mftTracksCursor, mftTracksCovCursor, ambigMFTTracksCursor,
2256 fwdTracksCursor, fwdTracksCovCursor, ambigFwdTracksCursor, fwdTrkClsCursor, bcsMap);
2257
2258 // filling collisions and tracks into tables
2259 collisionID = 0;
2260 collisionsCursor.reserve(primVertices.size());
2261 for (auto& vertex : primVertices) {
2262 auto& cov = vertex.getCov();
2263 auto& timeStamp = vertex.getTimeStamp(); // this is a relative time
2264 const double interactionTime = timeStamp.getTimeStamp() * 1E3; // mus to ns
2265 uint64_t globalBC = relativeTime_to_GlobalBC(interactionTime);
2266 uint64_t localBC = relativeTime_to_LocalBC(interactionTime);
2267 LOG(debug) << "global BC " << globalBC << " local BC " << localBC << " relative interaction time " << interactionTime;
2268 // collision timestamp in ns wrt the beginning of collision BC
2269 const float relInteractionTime = static_cast<float>(localBC * o2::constants::lhc::LHCBunchSpacingNS - interactionTime);
2270 auto item = bcsMap.find(globalBC);
2271 int bcID = -1;
2272 if (item != bcsMap.end()) {
2273 bcID = item->second;
2274 } else {
2275 LOG(fatal) << "Error: could not find a corresponding BC ID for a collision; BC = " << globalBC << ", collisionID = " << collisionID;
2276 }
2277 collisionsCursor(bcID,
2278 truncateFloatFraction(vertex.getX(), mCollisionPosition),
2279 truncateFloatFraction(vertex.getY(), mCollisionPosition),
2280 truncateFloatFraction(vertex.getZ(), mCollisionPosition),
2281 truncateFloatFraction(cov[0], mCollisionPositionCov),
2282 truncateFloatFraction(cov[1], mCollisionPositionCov),
2283 truncateFloatFraction(cov[2], mCollisionPositionCov),
2284 truncateFloatFraction(cov[3], mCollisionPositionCov),
2285 truncateFloatFraction(cov[4], mCollisionPositionCov),
2286 truncateFloatFraction(cov[5], mCollisionPositionCov),
2287 vertex.getFlags(),
2288 truncateFloatFraction(vertex.getChi2(), mCollisionPositionCov),
2289 vertex.getNContributors(),
2290 truncateFloatFraction(relInteractionTime, mCollisionPosition),
2291 truncateFloatFraction(timeStamp.getTimeStampError() * 1E3, mCollisionPositionCov));
2292 mVtxToTableCollID[collisionID] = mTableCollID++;
2293
2294 auto& trackRef = primVer2TRefs[collisionID];
2295 // passing interaction time in [ps]
2296 fillTrackTablesPerCollision(collisionID, globalBC, trackRef, primVerGIs, recoData, tracksCursor, tracksCovCursor, tracksExtraCursor, tracksQACursor, ambigTracksCursor,
2297 mftTracksCursor, mftTracksCovCursor, ambigMFTTracksCursor,
2298 fwdTracksCursor, fwdTracksCovCursor, ambigFwdTracksCursor, fwdTrkClsCursor, bcsMap);
2299 collisionID++;
2300 }
2301
2302 fillSecondaryVertices(recoData, v0sCursor, cascadesCursor, decay3BodyCursor);
2303 fillHMPID(recoData, hmpCursor);
2304 fillStrangenessTrackingTables(recoData, trackedV0Cursor, trackedCascadeCursor, tracked3BodyCurs);
2305
2306 // helper map for fast search of a corresponding class mask for a bc
2307 auto emcalIncomplete = filterEMCALIncomplete(recoData.getEMCALTriggers());
2308 std::unordered_map<uint64_t, std::pair<uint64_t, uint64_t>> bcToClassMask;
2309 if (mInputSources[GID::CTP]) {
2310 LOG(debug) << "CTP input available";
2311 for (auto& ctpDigit : ctpDigits) {
2312 uint64_t bc = ctpDigit.intRecord.toLong();
2313 uint64_t classMask = ctpDigit.CTPClassMask.to_ulong();
2314 uint64_t inputMask = ctpDigit.CTPInputMask.to_ulong();
2315 if (emcalIncomplete.find(bc) != emcalIncomplete.end()) {
2316 // reject EMCAL triggers as BC was rejected as incomplete at readout level
2317 auto classMaskOrig = classMask;
2318 classMask = classMask & ~mEMCALTrgClassMask;
2319 LOG(debug) << "Found EMCAL incomplete event, mask before " << std::bitset<64>(classMaskOrig) << ", after " << std::bitset<64>(classMask);
2320 }
2321 bcToClassMask[bc] = {classMask, inputMask};
2322 // LOG(debug) << Form("classmask:0x%llx", classMask);
2323 }
2324 }
2325
2326 // filling BC table
2327 bcCursor.reserve(bcsMap.size());
2328 for (auto& item : bcsMap) {
2329 uint64_t bc = item.first;
2330 std::pair<uint64_t, uint64_t> masks{0, 0};
2331 if (mInputSources[GID::CTP]) {
2332 auto bcClassPair = bcToClassMask.find(bc);
2333 if (bcClassPair != bcToClassMask.end()) {
2334 masks = bcClassPair->second;
2335 }
2336 }
2337 bcCursor(runNumber,
2338 bc,
2339 masks.first,
2340 masks.second);
2341 }
2342
2343 bcToClassMask.clear();
2344
2345 // filling BC flags table:
2346 auto bcFlags = fillBCFlags(recoData, bcsMap);
2347 bcFlagsCursor.reserve(bcFlags.size());
2348 for (auto f : bcFlags) {
2349 bcFlagsCursor(f);
2350 }
2351
2352 // fill cpvcluster table
2353 if (mInputSources[GIndex::CPV]) {
2354 float posX, posZ;
2355 cpvClustersCursor.reserve(cpvTrigRecs.size());
2356 for (auto& cpvEvent : cpvTrigRecs) {
2357 uint64_t bc = cpvEvent.getBCData().toLong();
2358 auto item = bcsMap.find(bc);
2359 int bcID = -1;
2360 if (item != bcsMap.end()) {
2361 bcID = item->second;
2362 } else {
2363 LOG(fatal) << "Error: could not find a corresponding BC ID for a CPV Trigger Record; BC = " << bc;
2364 }
2365 for (int iClu = cpvEvent.getFirstEntry(); iClu < cpvEvent.getFirstEntry() + cpvEvent.getNumberOfObjects(); iClu++) {
2366 auto& clu = cpvClusters[iClu];
2367 clu.getLocalPosition(posX, posZ);
2368 cpvClustersCursor(bcID,
2369 truncateFloatFraction(posX, mCPVPos),
2370 truncateFloatFraction(posZ, mCPVPos),
2371 truncateFloatFraction(clu.getEnergy(), mCPVAmpl),
2372 clu.getPackedClusterStatus());
2373 }
2374 }
2375 }
2376
2377 if (mUseMC) {
2378 TStopwatch timer;
2379 timer.Start();
2380 // filling mc particles table
2381 fillMCParticlesTable(*mcReader,
2382 mcParticlesCursor.cursor,
2383 primVer2TRefs,
2384 primVerGIs,
2385 recoData,
2386 mcColToEvSrc);
2387 timer.Stop();
2388 LOG(info) << "FILL MC took " << timer.RealTime() << " s";
2389 mcColToEvSrc.clear();
2390
2391 // ------------------------------------------------------
2392 // filling track labels
2393
2394 // need to go through labels in the same order as for tracks
2395 fillMCTrackLabelsTable(mcTrackLabelCursor, mcMFTTrackLabelCursor, mcFwdTrackLabelCursor, primVer2TRefs.back(), primVerGIs, recoData);
2396 for (auto iref = 0U; iref < primVer2TRefs.size() - 1; iref++) {
2397 auto& trackRef = primVer2TRefs[iref];
2398 fillMCTrackLabelsTable(mcTrackLabelCursor, mcMFTTrackLabelCursor, mcFwdTrackLabelCursor, trackRef, primVerGIs, recoData, iref);
2399 }
2400 }
2401
2402 // Fill calo tables and if MC also the MCCaloTable, therefore, has to be after fillMCParticlesTable call!
2403 if (mInputSources[GIndex::PHS] || mInputSources[GIndex::EMC]) {
2404 fillCaloTable(caloCellsCursor, caloCellsTRGTableCursor, mcCaloLabelsCursor, bcsMap, recoData);
2405 }
2406
2407 bcsMap.clear();
2408 clearMCKeepStore(mToStore);
2409 mGIDToTableID.clear();
2410 mTableTrID = 0;
2411 mGIDToTableFwdID.clear();
2412 mTableTrFwdID = 0;
2413 mGIDToTableMFTID.clear();
2414 mTableTrMFTID = 0;
2415 mVtxToTableCollID.clear();
2416 mTableCollID = 0;
2417 mV0ToTableID.clear();
2418 mTableV0ID = 0;
2419
2420 mIndexTableFwd.clear();
2421 mIndexFwdID = 0;
2422 mIndexTableMFT.clear();
2423 mIndexMFTID = 0;
2424
2425 mBCLookup.clear();
2426
2427 mGIDUsedBySVtx.clear();
2428 mGIDUsedByStr.clear();
2429
2430 originCursor(tfNumber);
2431
2432 // sending metadata to writer
2433 TString dataType = mUseMC ? "MC" : "RAW";
2434 TString O2Version = o2::fullVersion();
2435 TString ROOTVersion = ROOT_RELEASE;
2436 mMetaDataKeys = {"DataType", "Run", "O2Version", "ROOTVersion", "RecoPassName", "AnchorProduction", "AnchorPassName", "LPMProductionTag", "CreatedBy"};
2437 mMetaDataVals = {dataType, "3", O2Version, ROOTVersion, mRecoPass, mAnchorProd, mAnchorPass, mLPMProdTag, mUser};
2438 add_additional_meta_info(mMetaDataKeys, mMetaDataVals);
2439
2440 pc.outputs().snapshot(Output{"AMD", "AODMetadataKeys", 0}, mMetaDataKeys);
2441 pc.outputs().snapshot(Output{"AMD", "AODMetadataVals", 0}, mMetaDataVals);
2442
2443 pc.outputs().snapshot(Output{"TFN", "TFNumber", 0}, tfNumber);
2444 pc.outputs().snapshot(Output{"TFF", "TFFilename", 0}, "");
2445
2446 mTimer.Stop();
2447}
2448
2449void AODProducerWorkflowDPL::cacheTriggers(const o2::globaltracking::RecoContainer& recoData)
2450{
2451 // ITS tracks->ROF
2452 {
2453 mITSROFs.clear();
2454 const auto& rofs = recoData.getITSTracksROFRecords();
2455 uint16_t count = 0;
2456 for (const auto& rof : rofs) {
2457 int first = rof.getFirstEntry(), last = first + rof.getNEntries();
2458 for (int i = first; i < last; i++) {
2459 mITSROFs.push_back(count);
2460 }
2461 count++;
2462 }
2463 }
2464 // MFT tracks->ROF
2465 {
2466 mMFTROFs.clear();
2467 const auto& rofs = recoData.getMFTTracksROFRecords();
2468 uint16_t count = 0;
2469 for (const auto& rof : rofs) {
2470 int first = rof.getFirstEntry(), last = first + rof.getNEntries();
2471 for (int i = first; i < last; i++) {
2472 mMFTROFs.push_back(count);
2473 }
2474 count++;
2475 }
2476 }
2477 // ITSTPCTRD tracks -> TRD trigger
2478 {
2479 mITSTPCTRDTriggers.clear();
2480 const auto& itstpctrigs = recoData.getITSTPCTRDTriggers();
2481 int count = 0;
2482 for (const auto& trig : itstpctrigs) {
2483 int first = trig.getFirstTrack(), last = first + trig.getNumberOfTracks();
2484 for (int i = first; i < last; i++) {
2485 mITSTPCTRDTriggers.push_back(count);
2486 }
2487 count++;
2488 }
2489 }
2490 // TPCTRD tracks -> TRD trigger
2491 {
2492 mTPCTRDTriggers.clear();
2493 const auto& tpctrigs = recoData.getTPCTRDTriggers();
2494 int count = 0;
2495 for (const auto& trig : tpctrigs) {
2496 int first = trig.getFirstTrack(), last = first + trig.getNumberOfTracks();
2497 for (int i = first; i < last; i++) {
2498 mTPCTRDTriggers.push_back(count);
2499 }
2500 count++;
2501 }
2502 }
2503 // MCH tracks->ROF
2504 {
2505 mMCHROFs.clear();
2506 const auto& rofs = recoData.getMCHTracksROFRecords();
2507 uint16_t count = 0;
2508 for (const auto& rof : rofs) {
2509 int first = rof.getFirstIdx(), last = first + rof.getNEntries();
2510 for (int i = first; i < last; i++) {
2511 mMCHROFs.push_back(count);
2512 }
2513 count++;
2514 }
2515 }
2516}
2517
2518AODProducerWorkflowDPL::TrackExtraInfo AODProducerWorkflowDPL::processBarrelTrack(int collisionID, std::uint64_t collisionBC, GIndex trackIndex,
2519 const o2::globaltracking::RecoContainer& data, const std::map<uint64_t, int>& bcsMap)
2520{
2521 TrackExtraInfo extraInfoHolder;
2522 if (collisionID < 0) {
2523 extraInfoHolder.flags |= o2::aod::track::OrphanTrack;
2524 }
2525 bool needBCSlice = collisionID < 0; // track is associated to multiple vertices
2526 uint64_t bcOfTimeRef = collisionBC - mStartIR.toLong(); // by default track time is wrt collision BC (unless no collision assigned)
2527
2528 auto setTrackTime = [&](double t, double terr, bool gaussian) {
2529 // set track time and error, for ambiguous tracks define the bcSlice as it was used in vertex-track association
2530 // provided track time (wrt TF start) and its error should be in ns, gaussian flag tells if the error is assumed to be gaussin or half-interval
2531 if (!gaussian) {
2532 extraInfoHolder.flags |= o2::aod::track::TrackTimeResIsRange;
2533 }
2534 extraInfoHolder.trackTimeRes = terr;
2535 if (needBCSlice) { // need to define BC slice
2536 double error = this->mTimeMarginTrackTime + (gaussian ? extraInfoHolder.trackTimeRes * this->mNSigmaTimeTrack : extraInfoHolder.trackTimeRes);
2537 bcOfTimeRef = fillBCSlice(extraInfoHolder.bcSlice, t - error, t + error, bcsMap);
2538 }
2539 extraInfoHolder.trackTime = float(t - bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS);
2540 extraInfoHolder.diffBCRef = int(bcOfTimeRef);
2541 LOGP(debug, "time : {}/{} -> {}/{} -> trunc: {}/{} CollID: {} Amb: {}", t, terr, t - bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS, terr,
2542 truncateFloatFraction(extraInfoHolder.trackTime, mTrackTime), truncateFloatFraction(extraInfoHolder.trackTimeRes, mTrackTimeError),
2543 collisionID, trackIndex.isAmbiguous());
2544 };
2545 auto contributorsGID = data.getSingleDetectorRefs(trackIndex);
2546 const auto& trackPar = data.getTrackParam(trackIndex);
2547 extraInfoHolder.flags |= trackPar.getPID() << 28;
2548 auto src = trackIndex.getSource();
2549 if (contributorsGID[GIndex::Source::TOF].isIndexSet()) { // ITS-TPC-TRD-TOF, ITS-TPC-TOF, TPC-TRD-TOF, TPC-TOF
2550 const auto& tofMatch = data.getTOFMatch(trackIndex);
2551 extraInfoHolder.tofChi2 = tofMatch.getChi2();
2552 const auto& tofInt = tofMatch.getLTIntegralOut();
2553 float intLen = tofInt.getL();
2554 extraInfoHolder.length = intLen;
2555 const float mass = o2::constants::physics::MassPionCharged; // default pid = pion
2556 if (tofInt.getTOF(o2::track::PID::Pion) > 0.f) {
2557 float expBeta = (intLen / (tofInt.getTOF(o2::track::PID::Pion) * cSpeed));
2558 if (expBeta > o2::constants::math::Almost1) {
2560 }
2561 extraInfoHolder.tofExpMom = mass * expBeta / std::sqrt(1.f - expBeta * expBeta);
2562 }
2563 // correct the time of the track
2564 const double massZ = o2::track::PID::getMass2Z(trackPar.getPID());
2565 const double energy = sqrt((massZ * massZ) + (extraInfoHolder.tofExpMom * extraInfoHolder.tofExpMom));
2566 const double exp = extraInfoHolder.length * energy / (cSpeed * extraInfoHolder.tofExpMom);
2567 auto tofSignal = (tofMatch.getSignal() - exp) * 1e-3; // time in ns wrt TF start
2568 setTrackTime(tofSignal, 0.2, true); // FIXME: calculate actual resolution (if possible?)
2569 }
2570 if (contributorsGID[GIndex::Source::TRD].isIndexSet()) { // ITS-TPC-TRD-TOF, TPC-TRD-TOF, TPC-TRD, ITS-TPC-TRD
2571 const auto& trdOrig = data.getTrack<o2::trd::TrackTRD>(contributorsGID[GIndex::Source::TRD]); // refitted TRD trac
2572 extraInfoHolder.trdChi2 = trdOrig.getChi2();
2573 extraInfoHolder.trdSignal = trdOrig.getSignal();
2574 extraInfoHolder.trdPattern = getTRDPattern(trdOrig);
2575 if (extraInfoHolder.trackTimeRes < 0.) { // time is not set yet, this is possible only for TPC-TRD and ITS-TPC-TRD tracks, since those with TOF are set upstream
2576 // TRD is triggered: time uncertainty is within a BC
2577 const auto& trdTrig = (src == GIndex::Source::TPCTRD) ? data.getTPCTRDTriggers()[mTPCTRDTriggers[trackIndex.getIndex()]] : data.getITSTPCTRDTriggers()[mITSTPCTRDTriggers[trackIndex.getIndex()]];
2578 double ttrig = trdTrig.getBCData().differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingNS; // 1st get time wrt TF start
2579 setTrackTime(ttrig, 1., true); // FIXME: calculate actual resolution (if possible?)
2580 }
2581 }
2582 if (contributorsGID[GIndex::Source::ITS].isIndexSet()) {
2583 const auto& itsTrack = data.getITSTrack(contributorsGID[GIndex::ITS]);
2584 int nClusters = itsTrack.getNClusters();
2585 float chi2 = itsTrack.getChi2();
2586 extraInfoHolder.itsChi2NCl = nClusters != 0 ? chi2 / (float)nClusters : 0;
2587 extraInfoHolder.itsClusterSizes = itsTrack.getClusterSizes();
2588 if (src == GIndex::ITS) { // standalone ITS track should set its time from the ROF
2589 const auto& rof = data.getITSTracksROFRecords()[mITSROFs[trackIndex.getIndex()]];
2590 double t = rof.getBCData().differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingNS + mITSROFrameHalfLengthNS + mITSROFBiasNS;
2591 setTrackTime(t, mITSROFrameHalfLengthNS, false);
2592 }
2593 } else if (contributorsGID[GIndex::Source::ITSAB].isIndexSet()) { // this is an ITS-TPC afterburner contributor
2594 extraInfoHolder.itsClusterSizes = data.getITSABRefs()[contributorsGID[GIndex::Source::ITSAB].getIndex()].getClusterSizes();
2595 }
2596 if (contributorsGID[GIndex::Source::TPC].isIndexSet()) {
2597 const auto& tpcOrig = data.getTPCTrack(contributorsGID[GIndex::TPC]);
2598 const auto& tpcClData = mTPCCounters[contributorsGID[GIndex::TPC]];
2599 const auto& dEdx = tpcOrig.getdEdx().dEdxTotTPC > 0 ? tpcOrig.getdEdx() : tpcOrig.getdEdxAlt();
2600 if (tpcOrig.getdEdx().dEdxTotTPC == 0) {
2601 extraInfoHolder.flags |= o2::aod::track::TPCdEdxAlt;
2602 }
2603 extraInfoHolder.tpcInnerParam = tpcOrig.getP() / tpcOrig.getAbsCharge();
2604 extraInfoHolder.tpcChi2NCl = tpcOrig.getNClusters() ? tpcOrig.getChi2() / tpcOrig.getNClusters() : 0;
2605 extraInfoHolder.tpcSignal = dEdx.dEdxTotTPC;
2606 extraInfoHolder.tpcNClsFindable = tpcOrig.getNClusters();
2607 extraInfoHolder.tpcNClsFindableMinusFound = tpcOrig.getNClusters() - tpcClData.found;
2608 extraInfoHolder.tpcNClsFindableMinusCrossedRows = tpcOrig.getNClusters() - tpcClData.crossed;
2609 extraInfoHolder.tpcNClsShared = tpcClData.shared;
2610 uint32_t clsUsedForPID = dEdx.NHitsIROC + dEdx.NHitsOROC1 + dEdx.NHitsOROC2 + dEdx.NHitsOROC3;
2611 extraInfoHolder.tpcNClsFindableMinusPID = tpcOrig.getNClusters() - clsUsedForPID;
2612 if (src == GIndex::TPC) { // standalone TPC track should set its time from their timebins range
2613 if (needBCSlice) {
2614 double t = (tpcOrig.getTime0() + 0.5 * (tpcOrig.getDeltaTFwd() - tpcOrig.getDeltaTBwd())) * mTPCBinNS; // central value
2615 double terr = 0.5 * (tpcOrig.getDeltaTFwd() + tpcOrig.getDeltaTBwd()) * mTPCBinNS;
2616 double err = mTimeMarginTrackTime + terr;
2617 bcOfTimeRef = fillBCSlice(extraInfoHolder.bcSlice, t - err, t + err, bcsMap);
2618 }
2620 p.setDeltaTFwd(tpcOrig.getDeltaTFwd());
2621 p.setDeltaTBwd(tpcOrig.getDeltaTBwd());
2622 extraInfoHolder.trackTimeRes = p.getTimeErr();
2623 extraInfoHolder.trackTime = float(tpcOrig.getTime0() * mTPCBinNS - bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS);
2624 extraInfoHolder.diffBCRef = int(bcOfTimeRef);
2625 extraInfoHolder.isTPConly = true; // no truncation
2626 extraInfoHolder.flags |= o2::aod::track::TrackTimeAsym;
2627 } else if (src == GIndex::ITSTPC) { // its-tpc matched tracks have gaussian time error and the time was not set above
2628 const auto& trITSTPC = data.getTPCITSTrack(trackIndex);
2629 auto ts = trITSTPC.getTimeMUS();
2630 setTrackTime(ts.getTimeStamp() * 1.e3, ts.getTimeStampError() * 1.e3, true);
2631 }
2632 }
2633
2634 extrapolateToCalorimeters(extraInfoHolder, data.getTrackParamOut(trackIndex));
2635 // set bit encoding for PVContributor property as part of the flag field
2636 if (trackIndex.isPVContributor()) {
2637 extraInfoHolder.flags |= o2::aod::track::PVContributor;
2638 }
2639 return extraInfoHolder;
2640}
2641
2642AODProducerWorkflowDPL::TrackQA AODProducerWorkflowDPL::processBarrelTrackQA(int collisionID, std::uint64_t collisionBC, GIndex trackIndex,
2643 const o2::globaltracking::RecoContainer& data, const std::map<uint64_t, int>& bcsMap)
2644{
2645 TrackQA trackQAHolder;
2646 auto contributorsGID = data.getTPCContributorGID(trackIndex);
2647 const auto& trackPar = data.getTrackParam(trackIndex);
2648 if (contributorsGID.isIndexSet()) {
2649 auto prop = o2::base::Propagator::Instance();
2650 const auto& tpcOrig = data.getTPCTrack(contributorsGID);
2653 const o2::base::Propagator::MatCorrType mMatType = o2::base::Propagator::MatCorrType::USEMatCorrLUT;
2654 const o2::dataformats::VertexBase v = mVtx.getMeanVertex(collisionID < 0 ? 0.f : data.getPrimaryVertex(collisionID).getZ());
2655 std::array<float, 2> dcaInfo{-999., -999.};
2656 if (prop->propagateToDCABxByBz({v.getX(), v.getY(), v.getZ()}, tpcTMP, 2.f, mMatType, &dcaInfo)) {
2657 trackQAHolder.tpcdcaR = 100. * dcaInfo[0] / sqrt(1. + trackPar.getQ2Pt() * trackPar.getQ2Pt());
2658 trackQAHolder.tpcdcaZ = 100. * dcaInfo[1] / sqrt(1. + trackPar.getQ2Pt() * trackPar.getQ2Pt());
2659 }
2660 // This allows to safely clamp any float to one byte, using the
2661 // minmal/maximum values as under-/overflow borders and rounding to the nearest integer
2662 auto safeInt8Clamp = [](auto value) -> int8_t {
2663 using ValType = decltype(value);
2664 return static_cast<int8_t>(TMath::Nint(std::clamp(value, static_cast<ValType>(std::numeric_limits<int8_t>::min()), static_cast<ValType>(std::numeric_limits<int8_t>::max()))));
2665 };
2666 auto safeUInt8Clamp = [](auto value) -> uint8_t {
2667 using ValType = decltype(value);
2668 return static_cast<uint8_t>(TMath::Nint(std::clamp(value, static_cast<ValType>(std::numeric_limits<uint8_t>::min()), static_cast<ValType>(std::numeric_limits<uint8_t>::max()))));
2669 };
2670
2672 uint8_t clusterCounters[8] = {0};
2673 {
2674 uint8_t sectorIndex, rowIndex;
2675 uint32_t clusterIndex;
2676 const auto& tpcClusRefs = data.getTPCTracksClusterRefs();
2677 for (int i = 0; i < tpcOrig.getNClusterReferences(); i++) {
2678 o2::tpc::TrackTPC::getClusterReference(tpcClusRefs, i, sectorIndex, rowIndex, clusterIndex, tpcOrig.getClusterRef());
2679 char indexTracklet = (rowIndex % 152) / 19;
2680 clusterCounters[indexTracklet]++;
2681 }
2682 }
2683 uint8_t byteMask = 0;
2684 for (int i = 0; i < 8; i++) {
2685 if (clusterCounters[i] > 5) {
2686 byteMask |= (1 << i);
2687 }
2688 }
2689 trackQAHolder.tpcTime0 = tpcOrig.getTime0();
2690 trackQAHolder.tpcClusterByteMask = byteMask;
2691 const auto& dEdxInfoAlt = tpcOrig.getdEdxAlt(); // tpcOrig.getdEdx()
2692 const float dEdxNorm = (dEdxInfoAlt.dEdxTotTPC > 0) ? 100. / dEdxInfoAlt.dEdxTotTPC : 0;
2693 trackQAHolder.tpcdEdxNorm = dEdxInfoAlt.dEdxTotTPC;
2694 trackQAHolder.tpcdEdxMax0R = safeUInt8Clamp(dEdxInfoAlt.dEdxMaxIROC * dEdxNorm);
2695 trackQAHolder.tpcdEdxMax1R = safeUInt8Clamp(dEdxInfoAlt.dEdxMaxOROC1 * dEdxNorm);
2696 trackQAHolder.tpcdEdxMax2R = safeUInt8Clamp(dEdxInfoAlt.dEdxMaxOROC2 * dEdxNorm);
2697 trackQAHolder.tpcdEdxMax3R = safeUInt8Clamp(dEdxInfoAlt.dEdxMaxOROC3 * dEdxNorm);
2698 //
2699 trackQAHolder.tpcdEdxTot0R = safeUInt8Clamp(dEdxInfoAlt.dEdxTotIROC * dEdxNorm);
2700 trackQAHolder.tpcdEdxTot1R = safeUInt8Clamp(dEdxInfoAlt.dEdxTotOROC1 * dEdxNorm);
2701 trackQAHolder.tpcdEdxTot2R = safeUInt8Clamp(dEdxInfoAlt.dEdxTotOROC2 * dEdxNorm);
2702 trackQAHolder.tpcdEdxTot3R = safeUInt8Clamp(dEdxInfoAlt.dEdxTotOROC3 * dEdxNorm);
2704 float scaleTOF{0};
2705 auto contributorsGIDA = data.getSingleDetectorRefs(trackIndex);
2706 if (contributorsGIDA[GIndex::Source::TOF].isIndexSet()) { // ITS-TPC-TRD-TOF, ITS-TPC-TOF, TPC-TRD-TOF, TPC-TOF
2707 const auto& tofMatch = data.getTOFMatch(trackIndex);
2708 const float qpt = trackPar.getQ2Pt();
2710 trackQAHolder.dTofdX = safeInt8Clamp(tofMatch.getDXatTOF() / scaleTOF);
2711 trackQAHolder.dTofdZ = safeInt8Clamp(tofMatch.getDZatTOF() / scaleTOF);
2712 }
2713
2714 // Add matching information at a reference point (defined by
2715 // o2::aod::track::trackQARefRadius) in the same frame as the global track
2716 // without material corrections and error propagation
2717 if (auto itsContGID = data.getITSContributorGID(trackIndex); itsContGID.isIndexSet() && itsContGID.getSource() != GIndex::ITSAB) {
2718 const auto& itsOrig = data.getITSTrack(itsContGID);
2719 o2::track::TrackPar gloCopy = trackPar;
2720 o2::track::TrackPar itsCopy = itsOrig;
2721 o2::track::TrackPar tpcCopy = tpcOrig;
2722 if (prop->propagateToX(gloCopy, o2::aod::track::trackQARefRadius, prop->getNominalBz(), o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, mMatCorr) &&
2723 prop->propagateToAlphaX(tpcCopy, gloCopy.getAlpha(), o2::aod::track::trackQARefRadius, false, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, 1, mMatCorr) &&
2724 prop->propagateToAlphaX(itsCopy, gloCopy.getAlpha(), o2::aod::track::trackQARefRadius, false, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, 1, mMatCorr)) {
2725 // All tracks are now at the same radius and in the same frame and we can calculate the deltas wrt. to the global track
2726 // The scale is defined by the global track scaling depending on beta0
2727 const float beta0 = std::sqrt(std::min(50.f / tpcOrig.getdEdx().dEdxMaxTPC, 1.f));
2728 const float qpt = gloCopy.getQ2Pt();
2729 const float x = qpt / beta0;
2730 // scaling is defined as sigmaBins/sqrt(p0^2 + (p1 * q/pt / beta)^2)
2731 auto scaleCont = [&x](int i) -> float {
2733 };
2734 auto scaleGlo = [&x](int i) -> float {
2736 };
2737
2738 // Calculate deltas for contributors
2739 trackQAHolder.dRefContY = safeInt8Clamp((itsCopy.getY() - tpcCopy.getY()) * scaleCont(0));
2740 trackQAHolder.dRefContZ = safeInt8Clamp((itsCopy.getZ() - tpcCopy.getZ()) * scaleCont(1));
2741 trackQAHolder.dRefContSnp = safeInt8Clamp((itsCopy.getSnp() - tpcCopy.getSnp()) * scaleCont(2));
2742 trackQAHolder.dRefContTgl = safeInt8Clamp((itsCopy.getTgl() - tpcCopy.getTgl()) * scaleCont(3));
2743 trackQAHolder.dRefContQ2Pt = safeInt8Clamp((itsCopy.getQ2Pt() - tpcCopy.getQ2Pt()) * scaleCont(4));
2744 // Calculate deltas for global track against averaged contributors
2745 trackQAHolder.dRefGloY = safeInt8Clamp(((itsCopy.getY() + tpcCopy.getY()) * 0.5f - gloCopy.getY()) * scaleGlo(0));
2746 trackQAHolder.dRefGloZ = safeInt8Clamp(((itsCopy.getZ() + tpcCopy.getZ()) * 0.5f - gloCopy.getZ()) * scaleGlo(1));
2747 trackQAHolder.dRefGloSnp = safeInt8Clamp(((itsCopy.getSnp() + tpcCopy.getSnp()) * 0.5f - gloCopy.getSnp()) * scaleGlo(2));
2748 trackQAHolder.dRefGloTgl = safeInt8Clamp(((itsCopy.getTgl() + tpcCopy.getTgl()) * 0.5f - gloCopy.getTgl()) * scaleGlo(3));
2749 trackQAHolder.dRefGloQ2Pt = safeInt8Clamp(((itsCopy.getQ2Pt() + tpcCopy.getQ2Pt()) * 0.5f - gloCopy.getQ2Pt()) * scaleGlo(4));
2750 //
2751
2752 if (mStreamerFlags[AODProducerStreamerFlags::TrackQA]) {
2753 (*mStreamer) << "trackQA"
2754 << "trackITSOrig=" << itsOrig
2755 << "trackTPCOrig=" << tpcOrig
2756 << "trackITSTPCOrig=" << trackPar
2757 << "trackITSProp=" << itsCopy
2758 << "trackTPCProp=" << tpcCopy
2759 << "trackITSTPCProp=" << gloCopy
2760 << "refRadius=" << o2::aod::track::trackQARefRadius
2761 << "scaleBins=" << o2::aod::track::trackQAScaleBins
2762 << "scaleCont0=" << scaleCont(0)
2763 << "scaleCont1=" << scaleCont(1)
2764 << "scaleCont2=" << scaleCont(2)
2765 << "scaleCont3=" << scaleCont(3)
2766 << "scaleCont4=" << scaleCont(4)
2767 << "scaleGlo0=" << scaleGlo(0)
2768 << "scaleGlo1=" << scaleGlo(1)
2769 << "scaleGlo2=" << scaleGlo(2)
2770 << "scaleGlo3=" << scaleGlo(3)
2771 << "scaleGlo4=" << scaleGlo(4)
2772 << "trackQAHolder.tpcTime0=" << trackQAHolder.tpcTime0
2773 << "trackQAHolder.tpcdEdxNorm=" << trackQAHolder.tpcdEdxNorm
2774 << "trackQAHolder.tpcdcaR=" << trackQAHolder.tpcdcaR
2775 << "trackQAHolder.tpcdcaZ=" << trackQAHolder.tpcdcaZ
2776 << "trackQAHolder.tpcdcaClusterByteMask=" << trackQAHolder.tpcClusterByteMask
2777 << "trackQAHolder.tpcdEdxMax0R=" << trackQAHolder.tpcdEdxMax0R
2778 << "trackQAHolder.tpcdEdxMax1R=" << trackQAHolder.tpcdEdxMax1R
2779 << "trackQAHolder.tpcdEdxMax2R=" << trackQAHolder.tpcdEdxMax2R
2780 << "trackQAHolder.tpcdEdxMax3R=" << trackQAHolder.tpcdEdxMax3R
2781 << "trackQAHolder.tpcdEdxTot0R=" << trackQAHolder.tpcdEdxTot0R
2782 << "trackQAHolder.tpcdEdxTot1R=" << trackQAHolder.tpcdEdxTot1R
2783 << "trackQAHolder.tpcdEdxTot2R=" << trackQAHolder.tpcdEdxTot2R
2784 << "trackQAHolder.tpcdEdxTot3R=" << trackQAHolder.tpcdEdxTot3R
2785 << "trackQAHolder.dRefContY=" << trackQAHolder.dRefContY
2786 << "trackQAHolder.dRefContZ=" << trackQAHolder.dRefContZ
2787 << "trackQAHolder.dRefContSnp=" << trackQAHolder.dRefContSnp
2788 << "trackQAHolder.dRefContTgl=" << trackQAHolder.dRefContTgl
2789 << "trackQAHolder.dRefContQ2Pt=" << trackQAHolder.dRefContQ2Pt
2790 << "trackQAHolder.dRefGloY=" << trackQAHolder.dRefGloY
2791 << "trackQAHolder.dRefGloZ=" << trackQAHolder.dRefGloZ
2792 << "trackQAHolder.dRefGloSnp=" << trackQAHolder.dRefGloSnp
2793 << "trackQAHolder.dRefGloTgl=" << trackQAHolder.dRefGloTgl
2794 << "trackQAHolder.dRefGloQ2Pt=" << trackQAHolder.dRefGloQ2Pt
2795 << "trackQAHolder.dTofdX=" << trackQAHolder.dTofdX
2796 << "trackQAHolder.dTofdZ=" << trackQAHolder.dTofdZ
2797 << "scaleTOF=" << scaleTOF
2798 << "\n";
2799 }
2800 }
2801 }
2802 }
2803
2804 return trackQAHolder;
2805}
2806
2807bool AODProducerWorkflowDPL::propagateTrackToPV(o2::track::TrackParametrizationWithError<float>& trackPar,
2809 int colID)
2810{
2811 o2::dataformats::DCA dcaInfo;
2812 dcaInfo.set(999.f, 999.f, 999.f, 999.f, 999.f);
2813 o2::dataformats::VertexBase v = mVtx.getMeanVertex(colID < 0 ? 0.f : data.getPrimaryVertex(colID).getZ());
2814 return o2::base::Propagator::Instance()->propagateToDCABxByBz(v, trackPar, 2.f, mMatCorr, &dcaInfo);
2815}
2816
2817void AODProducerWorkflowDPL::extrapolateToCalorimeters(TrackExtraInfo& extraInfoHolder, const o2::track::TrackPar& track)
2818{
2819 constexpr float XEMCAL = 440.f, XPHOS = 460.f, XEMCAL2 = XEMCAL * XEMCAL;
2820 constexpr float ETAEMCAL = 0.75; // eta of EMCAL/DCAL with margin
2821 constexpr float ZEMCALFastCheck = 460.; // Max Z (with margin to check with straightline extrapolarion)
2822 constexpr float ETADCALINNER = 0.22; // eta of the DCAL PHOS Hole (at XEMCAL)
2823 constexpr float ETAPHOS = 0.13653194; // nominal eta of the PHOS acceptance (at XPHOS): -log(tan((TMath::Pi()/2 - atan2(63, 460))/2))
2824 constexpr float ETAPHOSMARGIN = 0.17946979; // etat of the PHOS acceptance with 20 cm margin (at XPHOS): -log(tan((TMath::Pi()/2 + atan2(63+20., 460))/2)), not used, for the ref only
2825 constexpr float ETADCALPHOSSWITCH = (ETADCALINNER + ETAPHOS) / 2; // switch to DCAL to PHOS check if eta < this value
2826 constexpr short SNONE = 0, SEMCAL = 0x1, SPHOS = 0x2;
2827 constexpr short SECTORTYPE[18] = {
2828 SNONE, SNONE, SNONE, SNONE, // 0:3
2829 SEMCAL, SEMCAL, SEMCAL, SEMCAL, SEMCAL, SEMCAL, // 3:9 EMCAL only
2830 SNONE, SNONE, // 10:11
2831 SPHOS, // 12 PHOS only
2832 SPHOS | SEMCAL, SPHOS | SEMCAL, SPHOS | SEMCAL, // 13:15 PHOS & DCAL
2833 SEMCAL, // 16 DCAL only
2834 SNONE // 17
2835 };
2836
2837 o2::track::TrackPar outTr{track};
2838 auto prop = o2::base::Propagator::Instance();
2839 // 1st propagate to EMCAL nominal radius
2840 float xtrg = 0;
2841 // quick check with straight line propagtion
2842 if (!outTr.getXatLabR(XEMCAL, xtrg, prop->getNominalBz(), o2::track::DirType::DirOutward) ||
2843 (std::abs(outTr.getZAt(xtrg, 0)) > ZEMCALFastCheck) ||
2844 !prop->PropagateToXBxByBz(outTr, xtrg, 0.95, 10, o2::base::Propagator::MatCorrType::USEMatCorrLUT)) {
2845 LOGP(debug, "preliminary step: does not reach R={} {}", XEMCAL, outTr.asString());
2846 return;
2847 }
2848 // we do not necessarilly reach wanted radius in a single propagation
2849 if ((outTr.getX() * outTr.getX() + outTr.getY() * outTr.getY() < XEMCAL2) &&
2850 (!outTr.rotateParam(outTr.getPhi()) ||
2851 !outTr.getXatLabR(XEMCAL, xtrg, prop->getNominalBz(), o2::track::DirType::DirOutward) ||
2852 !prop->PropagateToXBxByBz(outTr, xtrg, 0.95, 10, o2::base::Propagator::MatCorrType::USEMatCorrLUT))) {
2853 LOGP(debug, "does not reach R={} {}", XEMCAL, outTr.asString());
2854 return;
2855 }
2856 // rotate to proper sector
2857 int sector = o2::math_utils::angle2Sector(outTr.getPhiPos());
2858
2859 auto propExactSector = [&outTr, &sector, prop](float xprop) -> bool { // propagate exactly to xprop in the proper sector frame
2860 int ntri = 0;
2861 while (ntri < 2) {
2862 auto outTrTmp = outTr;
2863 float alpha = o2::math_utils::sector2Angle(sector);
2864 if ((std::abs(outTr.getZ()) > ZEMCALFastCheck) || !outTrTmp.rotateParam(alpha) ||
2865 !prop->PropagateToXBxByBz(outTrTmp, xprop, 0.95, 10, o2::base::Propagator::MatCorrType::USEMatCorrLUT)) {
2866 LOGP(debug, "failed on rotation to {} (sector {}) or propagation to X={} {}", alpha, sector, xprop, outTrTmp.asString());
2867 return false;
2868 }
2869 // make sure we are still in the target sector
2870 int sectorTmp = o2::math_utils::angle2Sector(outTrTmp.getPhiPos());
2871 if (sectorTmp == sector) {
2872 outTr = outTrTmp;
2873 break;
2874 }
2875 sector = sectorTmp;
2876 ntri++;
2877 }
2878 if (ntri == 2) {
2879 LOGP(debug, "failed to rotate to sector, {}", outTr.asString());
2880 return false;
2881 }
2882 return true;
2883 };
2884
2885 // we are at the EMCAL X, check if we are in the good sector
2886 if (!propExactSector(XEMCAL) || SECTORTYPE[sector] == SNONE) { // propagation failed or neither EMCAL not DCAL/PHOS
2887 return;
2888 }
2889
2890 // check if we are in a good eta range
2891 float r = std::sqrt(outTr.getX() * outTr.getX() + outTr.getY() * outTr.getY()), tg = std::atan2(r, outTr.getZ());
2892 float eta = -std::log(std::tan(0.5f * tg)), etaAbs = std::abs(eta);
2893 if (etaAbs > ETAEMCAL) {
2894 LOGP(debug, "eta = {} is off at EMCAL radius", eta, outTr.asString());
2895 return;
2896 }
2897 // are we in the PHOS hole (with margin)?
2898 if ((SECTORTYPE[sector] & SPHOS) && etaAbs < ETADCALPHOSSWITCH) { // propagate to PHOS radius
2899 if (!propExactSector(XPHOS)) {
2900 return;
2901 }
2902 r = std::sqrt(outTr.getX() * outTr.getX() + outTr.getY() * outTr.getY());
2903 tg = std::atan2(r, outTr.getZ());
2904 eta = -std::log(std::tan(0.5f * tg));
2905 } else if (!(SECTORTYPE[sector] & SEMCAL)) { // are in the sector with PHOS only
2906 return;
2907 }
2908 extraInfoHolder.trackPhiEMCAL = outTr.getPhiPos();
2909 extraInfoHolder.trackEtaEMCAL = eta;
2910 LOGP(debug, "eta = {} phi = {} sector {} for {}", extraInfoHolder.trackEtaEMCAL, extraInfoHolder.trackPhiEMCAL, sector, outTr.asString());
2911 //
2912}
2913
2914std::set<uint64_t> AODProducerWorkflowDPL::filterEMCALIncomplete(const gsl::span<const o2::emcal::TriggerRecord> triggers)
2915{
2916 std::set<uint64_t> emcalIncompletes;
2917 for (const auto& trg : triggers) {
2918 if (trg.getTriggerBits() & o2::emcal::triggerbits::Inc) {
2919 // trigger record masked at incomplete at readout level
2920 emcalIncompletes.insert(trg.getBCData().toLong());
2921 }
2922 }
2923 return emcalIncompletes;
2924}
2925
2926void AODProducerWorkflowDPL::updateTimeDependentParams(ProcessingContext& pc)
2927{
2929 static bool initOnceDone = false;
2930 if (!initOnceDone) { // this params need to be queried only once
2931 initOnceDone = true;
2932 // Note: DPLAlpideParam for ITS and MFT will be loaded by the RecoContainer
2933 mSqrtS = o2::base::GRPGeomHelper::instance().getGRPLHCIF()->getSqrtS();
2934 // apply settings
2937 std::bitset<3564> bs = bcf.getBCPattern();
2938 for (auto i = 0U; i < bs.size(); i++) {
2939 if (bs.test(i)) {
2941 }
2942 }
2943
2945 mITSROFrameHalfLengthNS = 0.5 * (grpECS->isDetContinuousReadOut(o2::detectors::DetID::ITS) ? alpParamsITS.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS : alpParamsITS.roFrameLengthTrig);
2948 mMFTROFrameHalfLengthNS = 0.5 * (grpECS->isDetContinuousReadOut(o2::detectors::DetID::MFT) ? alpParamsMFT.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS : alpParamsMFT.roFrameLengthTrig);
2950
2951 // RS FIXME: this is not yet fetched from the CCDB
2953 mTPCBinNS = elParam.ZbinWidth * 1.e3;
2954
2955 const auto& pvParams = o2::vertexing::PVertexerParams::Instance();
2956 mNSigmaTimeTrack = pvParams.nSigmaTimeTrack;
2957 mTimeMarginTrackTime = pvParams.timeMarginTrackTime * 1.e3;
2958 mFieldON = std::abs(o2::base::Propagator::Instance()->getNominalBz()) > 0.01;
2959
2960 pc.inputs().get<o2::ctp::CTPConfiguration*>("ctpconfig");
2961 }
2962 if (mPropTracks) {
2964 }
2965}
2966
2967//_______________________________________
2969{
2970 // Note: strictly speaking, for Configurable params we don't need finaliseCCDB check, the singletons are updated at the CCDB fetcher level
2971 if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) {
2972 if (matcher == ConcreteDataMatcher("GLO", "GRPMAGFIELD", 0)) {
2974 }
2975 return;
2976 }
2977 if (matcher == ConcreteDataMatcher("ITS", "ALPIDEPARAM", 0)) {
2978 LOG(info) << "ITS Alpide param updated";
2980 par.printKeyValues();
2981 return;
2982 }
2983 if (matcher == ConcreteDataMatcher("MFT", "ALPIDEPARAM", 0)) {
2984 LOG(info) << "MFT Alpide param updated";
2986 par.printKeyValues();
2987 return;
2988 }
2989 if (matcher == ConcreteDataMatcher("GLO", "MEANVERTEX", 0)) {
2990 LOG(info) << "Imposing new MeanVertex: " << ((const o2::dataformats::MeanVertexObject*)obj)->asString();
2991 mVtx = *(const o2::dataformats::MeanVertexObject*)obj;
2992 return;
2993 }
2994 if (matcher == ConcreteDataMatcher("CTP", "CTPCONFIG", 0)) {
2995 // construct mask with EMCAL trigger classes for rejection of incomplete triggers
2996 auto ctpconfig = *(const o2::ctp::CTPConfiguration*)obj;
2997 mEMCALTrgClassMask = 0;
2998 for (const auto& trgclass : ctpconfig.getCTPClasses()) {
2999 if (trgclass.cluster->maskCluster[o2::detectors::DetID::EMC]) {
3000 mEMCALTrgClassMask |= trgclass.classMask;
3001 }
3002 }
3003 LOG(info) << "Loaded EMCAL trigger class mask: " << std::bitset<64>(mEMCALTrgClassMask);
3004 }
3005}
3006
3007void AODProducerWorkflowDPL::addRefGlobalBCsForTOF(const o2::dataformats::VtxTrackRef& trackRef, const gsl::span<const GIndex>& GIndices,
3008 const o2::globaltracking::RecoContainer& data, std::map<uint64_t, int>& bcsMap)
3009{
3010 // Orphan tracks need to refer to some globalBC and for tracks with TOF this BC should be whithin an orbit
3011 // from the track abs time (to guarantee time precision). Therefore, we may need to insert some dummy globalBCs
3012 // to guarantee proper reference.
3013 // complete globalBCs by dummy entries necessary to provide BC references for TOF tracks with requested precision
3014 // to provide a reference for the time of the orphan tracks we should make sure that there are no gaps longer
3015 // than needed to store the time with sufficient precision
3016
3017 // estimate max distance in BCs between TOF time and eventual reference BC
3018 int nbitsFrac = 24 - (32 - o2::math_utils::popcount(mTrackTime)); // number of bits used to encode the fractional part of float truncated by the mask
3019 int nbitsLoss = std::max(0, int(std::log2(TOFTimePrecPS))); // allowed bit loss guaranteeing needed precision in PS
3020 assert(nbitsFrac > 1);
3021 std::uint64_t maxRangePS = std::uint64_t(0x1) << (nbitsFrac + nbitsLoss);
3022 int maxGapBC = maxRangePS / (o2::constants::lhc::LHCBunchSpacingNS * 1e3); // max gap in BCs allowing to store time with required precision
3023 LOG(info) << "Max gap of " << maxGapBC << " BCs to closest globalBC reference is needed for TOF tracks to provide precision of "
3024 << TOFTimePrecPS << " ps";
3025
3026 // check if there are tracks orphan tracks at all
3027 if (!trackRef.getEntries()) {
3028 return;
3029 }
3030 // the bscMap has at least TF start BC
3031 std::uint64_t maxBC = mStartIR.toLong();
3032 const auto& tofClus = data.getTOFClusters();
3033 for (int src = GIndex::NSources; src--;) {
3034 if (!GIndex::getSourceDetectorsMask(src)[o2::detectors::DetID::TOF]) { // check only tracks with TOF contribution
3035 continue;
3036 }
3037 int start = trackRef.getFirstEntryOfSource(src);
3038 int end = start + trackRef.getEntriesOfSource(src);
3039 for (int ti = start; ti < end; ti++) {
3040 auto& trackIndex = GIndices[ti];
3041 const auto& tofMatch = data.getTOFMatch(trackIndex);
3042 const auto& tofInt = tofMatch.getLTIntegralOut();
3043 float intLen = tofInt.getL();
3044 float tofExpMom = 0.;
3045 if (tofInt.getTOF(o2::track::PID::Pion) > 0.f) {
3046 float expBeta = (intLen / (tofInt.getTOF(o2::track::PID::Pion) * cSpeed));
3047 if (expBeta > o2::constants::math::Almost1) {
3049 }
3050 tofExpMom = o2::constants::physics::MassPionCharged * expBeta / std::sqrt(1.f - expBeta * expBeta);
3051 } else {
3052 continue;
3053 }
3054 double massZ = o2::track::PID::getMass2Z(data.getTrackParam(trackIndex).getPID());
3055 double energy = sqrt((massZ * massZ) + (tofExpMom * tofExpMom));
3056 double exp = intLen * energy / (cSpeed * tofExpMom);
3057 auto tofSignal = (tofMatch.getSignal() - exp) * 1e-3; // time in ns wrt TF start
3058 auto bc = relativeTime_to_GlobalBC(tofSignal);
3059
3060 auto it = bcsMap.lower_bound(bc);
3061 if (it == bcsMap.end() || it->first > bc + maxGapBC) {
3062 bcsMap.emplace_hint(it, bc, 1);
3063 LOG(debug) << "adding dummy BC " << bc;
3064 }
3065 if (bc > maxBC) {
3066 maxBC = bc;
3067 }
3068 }
3069 }
3070 // make sure there is a globalBC exceeding the max encountered bc
3071 if ((--bcsMap.end())->first <= maxBC) {
3072 bcsMap.emplace_hint(bcsMap.end(), maxBC + 1, 1);
3073 }
3074 // renumber BCs
3075 int bcID = 0;
3076 for (auto& item : bcsMap) {
3077 item.second = bcID;
3078 bcID++;
3079 }
3080}
3081
3082std::uint64_t AODProducerWorkflowDPL::fillBCSlice(int (&slice)[2], double tmin, double tmax, const std::map<uint64_t, int>& bcsMap) const
3083{
3084 // for ambiguous tracks (no or multiple vertices) we store the BC slice corresponding to track time window used for track-vertex matching,
3085 // see VertexTrackMatcher::extractTracks creator method, i.e. central time estimated +- uncertainty defined as:
3086 // 1) for tracks having a gaussian time error: PVertexerParams.nSigmaTimeTrack * trackSigma + PVertexerParams.timeMarginTrackTime
3087 // 2) for tracks having time uncertainty in a fixed time interval (TPC,ITS,MFT..): half of the interval + PVertexerParams.timeMarginTrackTime
3088 // The track time in the TrackExtraInfo is stored in ns wrt the collision BC for unambigous tracks and wrt bcSlice[0] for ambiguous ones,
3089 // with convention for errors: trackSigma in case (1) and half of the time interval for case (2) above.
3090
3091 // find indices of widest slice of global BCs in the map compatible with provided BC range. bcsMap is guaranteed to be non-empty.
3092 // We also assume that tmax >= tmin.
3093
3094 uint64_t bcMin = relativeTime_to_GlobalBC(tmin), bcMax = relativeTime_to_GlobalBC(tmax);
3095
3096 /*
3097 // brute force way of searching bcs via direct binary search in the map
3098 auto lower = bcsMap.lower_bound(bcMin), upper = bcsMap.upper_bound(bcMax);
3099
3100 if (lower == bcsMap.end()) {
3101 --lower;
3102 }
3103 if (upper != lower) {
3104 --upper;
3105 }
3106 slice[0] = std::distance(bcsMap.begin(), lower);
3107 slice[1] = std::distance(bcsMap.begin(), upper);
3108 */
3109
3110 // faster way to search in bunch crossing via the accelerated bunch crossing lookup structure
3111 auto p = mBCLookup.lower_bound(bcMin);
3112 // assuming that bcMax will be >= bcMin and close to bcMin; we can find
3113 // the upper bound quickly by lineary iterating from p.first to the point where
3114 // the time becomes larger than bcMax.
3115 // (if this is not the case we could determine it with a similar call to mBCLookup)
3116 auto& bcvector = mBCLookup.getBCTimeVector();
3117 auto upperindex = p.first;
3118 while (upperindex < bcvector.size() && bcvector[upperindex] <= bcMax) {
3119 upperindex++;
3120 }
3121 if (upperindex != p.first) {
3122 upperindex--;
3123 }
3124 slice[0] = p.first;
3125 slice[1] = upperindex;
3126
3127 auto bcOfTimeRef = p.second - this->mStartIR.toLong();
3128 LOG(debug) << "BC slice t:" << tmin << " " << slice[0]
3129 << " t: " << tmax << " " << slice[1]
3130 << " bcref: " << bcOfTimeRef;
3131 return bcOfTimeRef;
3132}
3133
3134std::vector<uint8_t> AODProducerWorkflowDPL::fillBCFlags(const o2::globaltracking::RecoContainer& data, std::map<uint64_t, int>& bcsMap) const
3135{
3136 std::vector<uint8_t> flags(bcsMap.size());
3137
3138 // flag BCs belonging to UPC mode ITS ROFs
3139 auto bcIt = bcsMap.begin();
3140 auto itsrofs = data.getITSTracksROFRecords();
3143 for (auto& rof : itsrofs) {
3144 if (!rof.getFlag(o2::itsmft::ROFRecord::VtxUPCMode)) {
3145 continue;
3146 }
3147 uint64_t globalBC0 = rof.getBCData().toLong() + bROF, globalBC1 = globalBC0 + lROF - 1;
3148 // BCs are sorted, iterate until the start of ROF
3149 while (bcIt != bcsMap.end()) {
3150 if (bcIt->first < globalBC0) {
3151 ++bcIt;
3152 continue;
3153 }
3154 if (bcIt->first > globalBC1) {
3155 break;
3156 }
3157 flags[bcIt->second] |= o2::aod::bc::ITSUPCMode;
3158 ++bcIt;
3159 }
3160 }
3161 return flags;
3162}
3163
3165{
3166 LOGF(info, "aod producer dpl total timing: Cpu: %.3e Real: %.3e s in %d slots",
3167 mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1);
3168
3169 mStreamer.reset();
3170}
3171
3172DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, bool enableStrangenessTracking, bool useMC, bool CTPConfigPerRun, bool enableFITextra)
3173{
3174 auto dataRequest = std::make_shared<DataRequest>();
3175 dataRequest->inputs.emplace_back("ctpconfig", "CTP", "CTPCONFIG", 0, Lifetime::Condition, ccdbParamSpec("CTP/Config/Config", CTPConfigPerRun));
3176
3177 dataRequest->requestTracks(src, useMC);
3178 dataRequest->requestPrimaryVertices(useMC);
3179 if (src[GID::CTP]) {
3180 dataRequest->requestCTPDigits(useMC);
3181 }
3182 if (enableSV) {
3183 dataRequest->requestSecondaryVertices(useMC);
3184 }
3185 if (enableStrangenessTracking) {
3186 dataRequest->requestStrangeTracks(useMC);
3187 LOGF(info, "requestStrangeTracks Finish");
3188 }
3189 if (src[GID::ITS]) {
3190 dataRequest->requestClusters(GIndex::getSourcesMask("ITS"), false);
3191 }
3192 if (src[GID::TPC]) {
3193 dataRequest->requestClusters(GIndex::getSourcesMask("TPC"), false); // no need to ask for TOF clusters as they are requested with TOF tracks
3194 }
3195 if (src[GID::TOF]) {
3196 dataRequest->requestTOFClusters(useMC);
3197 }
3198 if (src[GID::PHS]) {
3199 dataRequest->requestPHOSCells(useMC);
3200 }
3201 if (src[GID::TRD]) {
3202 dataRequest->requestTRDTracklets(false);
3203 }
3204 if (src[GID::EMC]) {
3205 dataRequest->requestEMCALCells(useMC);
3206 }
3207 if (src[GID::CPV]) {
3208 dataRequest->requestCPVClusters(useMC);
3209 }
3210
3211 auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(true, // orbitResetTime
3212 true, // GRPECS=true
3213 true, // GRPLHCIF
3214 true, // GRPMagField
3215 true, // askMatLUT
3217 dataRequest->inputs,
3218 true); // query only once all objects except mag.field
3219
3220 dataRequest->inputs.emplace_back("meanvtx", "GLO", "MEANVERTEX", 0, Lifetime::Condition, ccdbParamSpec("GLO/Calib/MeanVertex", {}, 1));
3221
3222 using namespace o2::aod;
3223 using namespace o2::aodproducer;
3224
3225 std::vector<OutputSpec> outputs{
3259 OutputSpec{"TFN", "TFNumber"},
3260 OutputSpec{"TFF", "TFFilename"},
3261 OutputSpec{"AMD", "AODMetadataKeys"},
3262 OutputSpec{"AMD", "AODMetadataVals"}};
3263
3264 if (useMC) {
3265 outputs.insert(outputs.end(),
3266 {OutputForTable<McCollisions>::spec(),
3267 OutputForTable<HepMCXSections>::spec(),
3268 OutputForTable<HepMCPdfInfos>::spec(),
3269 OutputForTable<HepMCHeavyIons>::spec(),
3270 OutputForTable<McMFTTrackLabels>::spec(),
3271 OutputForTable<McFwdTrackLabels>::spec(),
3272 OutputForTable<StoredMcParticles_001>::spec(),
3273 OutputForTable<McTrackLabels>::spec(),
3274 OutputForTable<McCaloLabels_001>::spec(),
3275 // todo: use addTableToOuput helper?
3276 // currently the description is MCCOLLISLABEL, so
3277 // the name in AO2D would be O2mccollislabel
3278 // addTableToOutput<McCollisionLabels>(outputs);
3279 {OutputLabel{"McCollisionLabels"}, "AOD", "MCCOLLISIONLABEL", 0, Lifetime::Timeframe}});
3280 }
3281
3282 return DataProcessorSpec{
3283 "aod-producer-workflow",
3284 dataRequest->inputs,
3285 outputs,
3286 AlgorithmSpec{adaptFromTask<AODProducerWorkflowDPL>(src, dataRequest, ggRequest, enableSV, useMC, enableFITextra)},
3287 Options{
3288 ConfigParamSpec{"run-number", VariantType::Int64, -1L, {"The run-number. If left default we try to get it from DPL header."}},
3289 ConfigParamSpec{"aod-timeframe-id", VariantType::Int64, -1L, {"Set timeframe number"}},
3290 ConfigParamSpec{"fill-calo-cells", VariantType::Int, 1, {"Fill calo cells into cell table"}},
3291 ConfigParamSpec{"enable-truncation", VariantType::Int, 1, {"Truncation parameter: 1 -- on, != 1 -- off"}},
3292 ConfigParamSpec{"lpmp-prod-tag", VariantType::String, "", {"LPMProductionTag"}},
3293 ConfigParamSpec{"anchor-pass", VariantType::String, "", {"AnchorPassName"}},
3294 ConfigParamSpec{"anchor-prod", VariantType::String, "", {"AnchorProduction"}},
3295 ConfigParamSpec{"reco-pass", VariantType::String, "", {"RecoPassName"}},
3296 ConfigParamSpec{"created-by", VariantType::String, "", {"Who created this AO2D"}},
3297 ConfigParamSpec{"nthreads", VariantType::Int, std::max(1, int(std::thread::hardware_concurrency() / 2)), {"Number of threads"}},
3298 ConfigParamSpec{"reco-mctracks-only", VariantType::Int, 0, {"Store only reconstructed MC tracks and their mothers/daughters. 0 -- off, != 0 -- on"}},
3299 ConfigParamSpec{"ctpreadout-create", VariantType::Int, 0, {"Create CTP digits from detector readout and CTP inputs. !=1 -- off, 1 -- on"}},
3300 ConfigParamSpec{"emc-select-leading", VariantType::Bool, false, {"Flag to select if only the leading contributing particle for an EMCal cell should be stored"}},
3301 ConfigParamSpec{"propagate-tracks", VariantType::Bool, false, {"Propagate tracks (not used for secondary vertices) to IP"}},
3302 ConfigParamSpec{"hepmc-update", VariantType::String, "always", {"When to update HepMC Aux tables: always - force update, never - never update, all - if all keys are present, any - when any key is present (not valid yet)"}},
3303 ConfigParamSpec{"propagate-muons", VariantType::Bool, false, {"Propagate muons to IP"}},
3304 ConfigParamSpec{"thin-tracks", VariantType::Bool, false, {"Produce thinned track tables"}},
3305 ConfigParamSpec{"trackqc-fraction", VariantType::Float, float(0.1), {"Fraction of tracks to QC"}},
3306 ConfigParamSpec{"trackqc-NTrCut", VariantType::Int64, 4L, {"Minimal length of the track - in amount of tracklets"}},
3307 ConfigParamSpec{"trackqc-tpc-dca", VariantType::Float, 3.f, {"Keep TPC standalone track with this DCAxy to the PV"}},
3308 ConfigParamSpec{"trackqc-tpc-cls", VariantType::Int, 80, {"Keep TPC standalone track with this #clusters"}},
3309 ConfigParamSpec{"trackqc-tpc-pt", VariantType::Float, 0.2f, {"Keep TPC standalone track with this pt"}},
3310 ConfigParamSpec{"with-streamers", VariantType::String, "", {"Bit-mask to steer writing of intermediate streamer files"}},
3311 ConfigParamSpec{"seed", VariantType::Int, 0, {"Set seed for random generator used for sampling (0 (default) means using a random_device)"}},
3312 }};
3313}
3314
3315} // namespace o2::aodproducer
Class to refer to the reconstructed information.
Definition of the 32 Central Trigger System (CTS) Trigger Types defined in https://twiki....
General auxilliary methods.
definition of CTPConfiguration and related CTP structures
Wrapper container for different reconstructed object types.
Definition of the MCH cluster minimal structure.
Reconstructed MID track.
Base track model for the Barrel, params only, w/o covariance.
uint64_t exp(uint64_t base, uint8_t exp) noexcept
definition of CTPDigit, CTPInputDigit
int16_t charge
Definition RawEventData.h:5
uint64_t vertex
Definition RawEventData.h:9
uint64_t bc
Definition RawEventData.h:5
int16_t time
Definition RawEventData.h:4
Definition of the FIT RecPoints class.
Definition of the FV0 RecPoints class.
int32_t i
std::string getType()
Definition Utils.h:138
Header of the AggregatedRunInfo struct.
Global Forward Muon tracks.
Global index for barrel track: provides provenance (detectors combination), index in respective array...
Definition of the MCTrack class.
Definition of a container to keep Monte Carlo truth external to simulation objects.
Utility functions for MC particles.
Definition of the MCH ROFrame record.
Class to perform MFT MCH (and MID) matching.
Class to store the output of the matching to HMPID.
Class to perform TOF matching to global tracks.
Definition of the parameter class for the detector electronics.
Header to collect physics constants.
uint32_t j
Definition RawData.h:0
Definition of the FDD RecPoint class.
Definition of tools for track extrapolation.
Definition of the ITS track.
Definition of the MUON track.
Definition of the MCH track.
Definition of the MCH track parameters for internal use.
Result of refitting TPC-ITS matched track.
std::ostringstream debug
Extention of GlobalTrackID by flags relevant for verter-track association.
Referenc on track indices contributing to the vertex, with possibility chose tracks from specific sou...
Container class to store energy released in the ZDC.
Container class to store a TDC hit in a ZDC channel.
int nClusters
StringRef key
const auto & getBCPattern() const
void GetStartVertex(TVector3 &vertex) const
Definition MCTrack.h:326
void endOfStream(framework::EndOfStreamContext &ec) final
This is invoked whenever we have an EndOfStream event.
void finaliseCCDB(ConcreteDataMatcher &matcher, void *obj) final
std::pair< size_t, uint64_t > lower_bound(uint64_t timestamp) const
void init(std::map< uint64_t, int > const &bcs)
initialize this container (to be ready for lookup/search queries)
void clear()
clear/reset this container
std::vector< uint64_t > const & getBCTimeVector() const
return the sorted vector of increaing BC times
static float getAmplitude(const o2::emcal::Cell &cell)
static int16_t getLnAmplitude(const o2::emcal::Cell &)
static int8_t getTriggerBits(const o2::emcal::Cell &)
static int16_t getCellNumber(const o2::emcal::Cell &cell)
static int16_t getFastOrAbsID(const o2::emcal::Cell &)
static bool isTRU(const o2::emcal::Cell &cell)
static float getTimeStamp(const o2::emcal::Cell &cell)
void checkUpdates(o2::framework::ProcessingContext &pc)
static GRPGeomHelper & instance()
void setRequest(std::shared_ptr< GRPGeomRequest > req)
static constexpr float MAX_SIN_PHI
Definition Propagator.h:72
static constexpr float MAX_STEP
Definition Propagator.h:73
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:143
void printStream(std::ostream &stream) const
static mask_t getSourcesMask(const std::string_view srcList)
VertexBase getMeanVertex(float z) const
int getEntriesOfSource(int s) const
Definition VtxTrackRef.h:62
int getFirstEntryOfSource(int s) const
Definition VtxTrackRef.h:55
Static class with identifiers, bitmasks and names for ALICE detectors.
Definition DetID.h:58
static constexpr ID ITS
Definition DetID.h:63
static constexpr ID MFT
Definition DetID.h:71
static constexpr ID TPC
Definition DetID.h:64
static constexpr ID EMC
Definition DetID.h:69
static constexpr ID TOF
Definition DetID.h:66
Handler for EMCAL event data.
void reset()
Reset containers with empty ranges.
void setCellData(CellRange cells, TriggerRange triggers)
Setting the data at cell level.
void setCellMCTruthContainer(const o2::dataformats::MCTruthContainer< o2::emcal::MCLabel > *mclabels)
Setting the pointer for the MCTruthContainer for cells.
InteractionRecord getInteractionRecordForEvent(int eventID) const
Get the interaction record for the given event.
int getNumberOfEvents() const
Get the number of events handled by the event handler.
void snapshot(const Output &spec, T const &object)
ConfigParamRegistry const & options()
Definition InitContext.h:33
decltype(auto) get(R binding, int part=0) const
DataAllocator & outputs()
The data allocator is used to allocate memory for the output data.
InputRecord & inputs()
The inputs associated with this processing context.
ServiceRegistryRef services()
The services registry associated with this processing context.
static constexpr int NCellsA
Definition Geometry.h:52
o2::dataformats::GlobalFwdTrack MCHtoFwd(const o2::mch::TrackParam &mchTrack)
Converts mchTrack parameters to Forward coordinate system.
static bool extrapToZ(TrackParam &trackParam, double zEnd)
static void setField()
static bool extrapToVertex(TrackParam &trackParam, double xVtx, double yVtx, double zVtx, double errXVtx, double errYVtx)
Definition TrackExtrap.h:62
static bool extrapToVertexWithoutBranson(TrackParam &trackParam, double zVtx, double xUpstream=0., double yUpstream=0., std::optional< double > zUpstream=std::nullopt)
Definition TrackExtrap.h:74
track parameters for internal use
Definition TrackParam.h:34
Double_t getNonBendingCoor() const
return non bending coordinate (cm)
Definition TrackParam.h:51
Double_t getBendingCoor() const
return bending coordinate (cm)
Definition TrackParam.h:59
void setCellMCTruthContainer(const o2::dataformats::MCTruthContainer< o2::phos::MCLabel > *mclabels)
Setting the pointer for the MCTruthContainer for cells.
void setCellData(CellRange cells, TriggerRange triggers)
Setting the data at cell level.
void reset()
Reset containers with empty ranges.
InteractionRecord getInteractionRecordForEvent(int eventID) const
int getNumberOfEvents() const
MCTrack const * getTrack(o2::MCCompLabel const &) const
void releaseTracksForSourceAndEvent(int source, int event)
API to ask releasing tracks (freeing memory) for source + event.
std::vector< MCTrack > const & getTracks(int source, int event) const
variant returning all tracks for source and event at once
static void addInteractionBC(int bc, bool fromCollisonCotext=false)
Definition Utils.cxx:52
static constexpr ID Pion
Definition PID.h:96
Double_t getX() const
Definition TrackFwd.h:49
void set(const std::string &s="", int base=2)
Definition EnumFlags.h:401
bool match(const std::vector< std::string > &queries, const char *pattern)
Definition dcs-ccdb.cxx:229
struct _cl_event * event
Definition glcorearb.h:2982
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
GLint GLenum GLint x
Definition glcorearb.h:403
GLenum src
Definition glcorearb.h:1767
GLint GLsizei count
Definition glcorearb.h:399
GLuint entry
Definition glcorearb.h:5735
GLint GLint GLsizei GLuint * counters
Definition glcorearb.h:3985
GLuint GLuint end
Definition glcorearb.h:469
const GLdouble * v
Definition glcorearb.h:832
GLdouble GLdouble right
Definition glcorearb.h:4077
GLdouble f
Definition glcorearb.h:310
GLuint GLuint GLfloat weight
Definition glcorearb.h:5477
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLsizei GLsizei GLchar * source
Definition glcorearb.h:798
GLsizei const GLfloat * value
Definition glcorearb.h:819
GLint GLint GLsizei GLint GLenum GLenum type
Definition glcorearb.h:275
GLenum GLsizei GLsizei GLint * values
Definition glcorearb.h:1576
GLboolean * data
Definition glcorearb.h:298
GLintptr offset
Definition glcorearb.h:660
GLuint GLsizei const GLchar * label
Definition glcorearb.h:2519
GLfloat v0
Definition glcorearb.h:811
GLbitfield flags
Definition glcorearb.h:1570
GLboolean r
Definition glcorearb.h:1233
GLuint start
Definition glcorearb.h:469
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
constexpr std::array< float, 2 > trackQAScaledTOF
Definition DataTypes.h:135
constexpr std::array< float, 5 > trackQAScaleContP1
Definition DataTypes.h:132
uint8_t itsSharedClusterMap uint8_t
constexpr std::array< float, 5 > trackQAScaleContP0
Definition DataTypes.h:131
constexpr std::array< float, 5 > trackQAScaleGloP0
Definition DataTypes.h:133
constexpr std::array< float, 5 > trackQAScaleGloP1
Definition DataTypes.h:134
constexpr float trackQAScaleBins
Definition DataTypes.h:129
constexpr float trackQARefRadius
Definition DataTypes.h:128
bool updateHepMCHeavyIon(const HeavyIonCursor &cursor, int collisionID, short generatorID, o2::dataformats::MCEventHeader const &header, HepMCUpdate when=HepMCUpdate::anyKey)
short updateMCCollisions(const CollisionCursor &cursor, int bcId, float time, o2::dataformats::MCEventHeader const &header, short generatorId=0, int sourceId=0, unsigned int mask=0xFFFFFFF0)
bool updateHepMCPdfInfo(const PdfInfoCursor &cursor, int collisionID, short generatorID, o2::dataformats::MCEventHeader const &header, HepMCUpdate when=HepMCUpdate::anyKey)
uint32_t updateParticles(const ParticleCursor &cursor, int collisionID, std::vector< MCTrack > const &tracks, TrackToIndex &preselect, uint32_t offset=0, bool filter=false, bool background=false, uint32_t weightMask=0xFFFFFFF0, uint32_t momentumMask=0xFFFFFFF0, uint32_t positionMask=0xFFFFFFF0)
bool updateHepMCXSection(const XSectionCursor &cursor, int collisionID, short generatorID, o2::dataformats::MCEventHeader const &header, HepMCUpdate when=HepMCUpdate::anyKey)
void keepMCParticle(std::vector< std::vector< std::unordered_map< int, int > > > &store, int source, int event, int track, int value=1)
framework::DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, bool enableST, bool useMC, bool CTPConfigPerRun, bool enableFITextra)
create a processor spec
void dimensionMCKeepStore(std::vector< std::vector< std::unordered_map< int, int > > > &store, int Nsources, int NEvents)
void clearMCKeepStore(std::vector< std::vector< std::unordered_map< int, int > > > &store)
constexpr int LHCMaxBunches
constexpr double LHCBunchSpacingNS
constexpr float Almost1
constexpr double MassPionCharged
Defining PrimaryVertex explicitly as messageable.
Definition TFIDInfo.h:20
std::vector< ConfigParamSpec > ccdbParamSpec(std::string const &path, int runDependent, std::vector< CCDBMetadata > metadata={}, int qrate=0)
int angle2Sector(float phi)
Definition Utils.h:183
float sector2Angle(int sect)
Definition Utils.h:193
TrackParCovF TrackParCov
Definition Track.h:33
constexpr uint32_t Cal
Definition Triggers.h:32
const int TDCSignal[NTDCChannels]
Definition Constants.h:181
constexpr int NTDCChannels
Definition Constants.h:90
constexpr int NChannels
Definition Constants.h:65
std::string fullVersion()
get full version information (official O2 release and git commit)
GPUReconstruction * rec
static OutputSpec const spec()
decltype(FFL(std::declval< cursor_t >())) cursor
gsl::span< const o2::trd::TriggerRecord > getTRDTriggerRecords() const
void collectData(o2::framework::ProcessingContext &pc, const DataRequest &request)
static bool downsampleTsallisCharged(float pt, float factorPt, float sqrts, float &weight, float rnd, float mass=0.13957)
Definition Tsallis.cxx:31
std::vector< o2::ctf::BufferType > vec
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
Cluster clu
std::uniform_int_distribution< unsigned long long > distr
std::random_device rd
std::array< uint16_t, 5 > pattern